package uk.ac.leeds.ccg.cluster;
import java.awt.*;
import java.awt.event.*;
import java.util.*;
//import com.sun.java.util.collections.*; //change for 1.2!
import uk.ac.leeds.ccg.geotools.*;

/** A class that implements Besag and Newell Cluster Hunt
 * The main class BandN creates the main frame and proivides 3 buttons
 * modified to use the new results object - 31/01/01
 */
public class Spattimescan extends SimpleMethod implements Method, Runnable
{
    // Members
		SpattimescanParameters parameters;
    SpattimescanParametersViewer parametersViewer;
    Database database;
    Results results;
    Viewer viewer;                        // Geotools viewer
    GeoRectangle geoRectangle = null;     // Geotools geoRectangle
    Shader shader;                        // Geotools shader
    
    // Methods
    public static void main(String args[]){
        Spattimescan spattimescan = new Spattimescan();
        spattimescan.run();
    }
    
    public Spattimescan(){
	    parameters = new SpattimescanParameters();
	    parametersViewer = new SpattimescanParametersViewer(parameters);
    }

    public String toString() {return "Spattimescan";}
    public Parameters getParameters() {return parameters;}
    public ParametersViewer getParametersViewer() {return parametersViewer;}
    public void setShader(Shader s) {shader=s;}
    public void setViewer(Viewer v) {viewer=v;}
    public void setData(Database d) {database=d;}

    public SignificanceTest getSignificanceTest(){return significanceTest;}
    
    public void setSearch(GeoRectangle r) {geoRectangle=r;}
    public Results getResults(){return results;}
    
    public void run(){
	results = new Results();
	double radMin = parameters.getRadMin();
	double radMax = parameters.getRadMax();
	double radInc = parameters.getRadInc();
	double timeMin = parameters.getTimeMin();
	double timeMax = parameters.getTimeMax();
	double timeInc = parameters.getTimeInc();
	int numRuns = parameters.getNumRuns();
	double thresh = parameters.getSignificanceThreshold();
        MyPoint p, p2;
	double n_g = (double) 0.0;
	double mu_g = (double) 0.0;
	double lambda = (double) 0.0;
	double maxLam = (double) 0.0;
	double maxLam2 = (double) 0.0;
	int nge = 0;
        double stat = (double) 0.0;
	    System.out.println("Spattimescan "+significanceTest);
        if (significanceTest == null) {
	    System.out.println("Spattimescan started with a null sig test using Poisson");
        significanceTest = new PoissonTest();
        significanceTest.setParameters(parameters);
        }
        System.out.println(""+this+" "+significanceTest);
        if(geoRectangle==null) geoRectangle=database.getBounds();
        Vector tcases = database.getCases();
        Vector cases= new Vector();
        for(int i=0;i<tcases.size();i++){
            p=(MyPoint)tcases.elementAt(i);
            if(geoRectangle.contains(p)) cases.addElement(p);
	    n_g += p.cases;
// System.out.println(p +" "+ n_g);
        }
Vector tcase = database.getPoints();
// System.out.println(tcase.elementAt(0));
for(int i=0;i<tcase.size();i++){
  p=(MyPoint)tcase.elementAt(i);
  mu_g += p.pop;
}
	// mu_g = n_g;  // is this right????
	p2 = (MyPoint) tcases.elementAt(0);  //need p2 initialised
System.out.println("# cases " + tcase.size() +" "+ cases.size());
System.out.println("# cases " + n_g +" "+ mu_g);
        // Go through the shapes vector that the canvas holds & extract the points
        // Calculate the overall rate of the disease

	GeoCircle gc = new GeoCircle();
	int rTimes = (int) ((radMax - radMin)/radInc +1.0);
	int tTimes = (int) ((timeMax - timeMin)/timeInc +1.0);
	double radius = radMin - radInc;
	double tidius = timeMin - timeInc;
	for (int loop=0; loop<rTimes; loop++){
	  radius += radInc;
	  for (int toop=0; toop<tTimes; toop++){
	    tidius += timeInc;
            for(int i=0;i<cases.size();i++){
                notifyProgressEventListener((double)i*100.0/(double)cases.size());
                // Set point p to be the element of the case vector at i
                p=(MyPoint)cases.elementAt(i);
// System.out.println(""+ p +" "+ p.x +" "+ p.y);
                gc.setCentre(p.x, p.y);
                gc.setRadius(radius);
                Vector points = database.getPoints(gc);
		double obsP = (double) 0.0;
		double obsC = (double) 0.0;
		int l2 = points.size();
		for (int i2=0; i2<l2; i2++){
	if ( ((MyPoint)points.elementAt(i2)).getTime() < timeMin  || 
	    ((MyPoint)points.elementAt(i2)).getTime() > timeMax ) continue;
	            obsP += ((MyPoint)points.elementAt(i2)).getPop();  //mu_z
	    	    obsC += ((MyPoint)points.elementAt(i2)).getCases();  //n_z
		}
// System.out.println("loop "+loop+" L2 "+l2+" obsP "+obsP+" obsC "+obsC);
if (obsC == 0) {
    // System.out.println("No cases in this interval");
    continue;
}
		double num1 = obsC / obsP;
		double num2 = (n_g - obsC) / (mu_g - obsP);
		double denom = n_g / mu_g;
		double fac1 = num1 /num2;
		double fac2 = num2 / denom;
// extra bits for binomial rather than Poisson
// denominator equals 1 as tot popn rescaled to equal tot no of cases!
		if (num1 < num2) continue;
		// double lambda = Math.pow(num1,obsC) *Math.pow(num2, (n_g-obsC));
		lambda = Math.pow(fac1,obsC) *Math.pow(fac2, n_g);
		// lambda /= Math.pow(denom, n_g);
// System.out.println("fac1 "+fac1+" fac2 "+ fac2 +" "+  Math.pow(fac1,obsC) +" "+  Math.pow(fac2, n_g) );
		if (lambda > maxLam){
    		    maxLam = lambda;
    		    p2 = p;
System.out.println("num1 "+num1+" num2 "+num2+" lambda "+lambda);
		}

// what's all this lot??
// this method just finds one cluster (at present) so I have commented out 
// the stuff that marks it(them) on a map - might have to check what sort
// of value is returned as 'stat' if/when I allow more than one cluster!
		// try{
			// if(significanceTest.isSignificant(points))
			// {
				// stat = significanceTest.getStatistic();
				// // results.addResult(cc,stat);
				// results.addResult(gc,stat);
			// }
		// }catch(SignificanceException e){
			// System.out.println(""+e);
		// }
            }
	  }
	}
if (lambda != 0) System.out.println(p2);
else System.out.println("No cases! - check time/circle radius parameters?");

// Here starts the Monte Carlo significance part
	Vector orig = database.getPoints();
// System.out.println(orig.size());
	int no_points = orig.size();
	// numRuns = 2; //just for the time being!
	for (int k=0; k<numRuns; k++){
if (lambda == 0) break;
	    Vector pts = new Vector();
	    Vector ncase = new Vector();
	    Vector newcases = new Vector();
	    for (int i=0; i<no_points; i++){
		ncase.add(new Integer(((MyPoint)orig.elementAt(i)).getCases() ));
            }
            Collections.shuffle(ncase);
            for (int i=0; i<no_points; i++){
		int c = ((Integer) ncase.elementAt(i)).intValue();
		MyPoint np = (MyPoint) ((MyPoint) orig.elementAt(i)).clone();
		np.setCases(c);
		pts.add(np);
if (np.getCases() != 0) newcases.addElement(np);
	    }

// that was random re-assigning of cases (without replacement!!), now calc lambda
	    radius = radMin - radInc;
	    tidius = timeMin - timeInc;
	    maxLam2 = 0.0;
            for (int loop=0; loop<rTimes; loop++){
              radius += radInc;
	      for (int toop=0; toop<tTimes; toop++){
		tidius += timeInc;
                for(int i=0;i<newcases.size();i++){
                    // Set point p to be the element of the case vector at i
                    p=(MyPoint)newcases.elementAt(i);
    // System.out.println(""+ p +" "+ p.x +" "+ p.y);
                    gc.setCentre(p.x, p.y);
                    gc.setRadius(radius);
                    // Vector points = database.getPoints(gc);
                    Vector newpoints = new Vector();
                    for(int j=0;j<no_points;j++){
                        p = (MyPoint) pts.elementAt(j);
                        if (gc.contains(p)) newpoints.addElement(p);
                    }
                    double obsP = (double) 0.0;
                    double obsC = (double) 0.0;
                    int l2 = newpoints.size();
                    for (int i2=0; i2<l2; i2++){
        if ( ((MyPoint)newpoints.elementAt(i2)).getTime() < timeMin  ||
            ((MyPoint)newpoints.elementAt(i2)).getTime() > timeMax ) continue;
                        obsP += ((MyPoint)newpoints.elementAt(i2)).getPop();  //mu_z
                        obsC += ((MyPoint)newpoints.elementAt(i2)).getCases();  //n_z
                    }
if (obsC == 0) {
    // System.out.println("No cases in this interval");
    continue;
}
    // System.out.println("loop "+loop+" l2 "+l2+" obsP "+obsP+" obsC "+obsC);
                    double num1 = obsC / obsP;
                    double num2 = (n_g - obsC) / (mu_g - obsP);
		    double denom = n_g / mu_g;
		    double fac1 = num1 / num2;
		    double fac2 = num2 / denom;
                    if (num1 < num2) continue;
                    double lamb2 = Math.pow(fac1,obsC) *Math.pow(fac2,n_g);
                    if (lamb2 > maxLam2){
                        maxLam2 = lamb2;
    // System.out.println("num1 "+num1+" num2 "+num2+" lamb2 "+lamb2);
                    }
                }
	      }
            }
    	    if (maxLam2 >= maxLam) nge++;
	    if ( (double) (nge+1)/(numRuns+1) > thresh){
    	        System.out.println("Will not be significant, nge = " + nge +","+ k +" out of NRuns");
    	        break;
	    }
	}

	if (lambda != 0){
	    if ( (double) (nge+1)/(numRuns+1) > thresh) System.out.println("Prob. of such a result is greater than " + (double) (nge+1) /(numRuns+1));
 	    else System.out.println("Prob. of such a result is less than, or equal to " + (double) (nge+1)/(numRuns+1));
	}
        notifyProgressEventListener(0.0);
    }
}