package lehre;

import java.util.*;
import java.io.*;



public class RS_Blatt3 {

	static Random r = new Random();
	
	// Erzeugt eine mit Parameter lambda Poisson-verteilte ZV
	// als Summe von Exp-verteilten ZV 
	public static int poisson(double lambda)
	{
		int count = 0;
		double e = -Math.log(r.nextDouble());
		while (e <= lambda) 
		{
		  	count++;
			e -= Math.log(r.nextDouble());
		}
		return count;
	}
	
	// Erzeugt eine Realisierung eines Poisson-Prozesses mit Intensitaetsfunktion f
	// im Fenster (x_axis[0],x_axis[1]) X (y_axis[0], y_axis[1])
	public static Vector PoissonProcess(double[] x_axis, double[] y_axis, intensityFunction f) 
	{
		Vector points = new Vector();
		double x_length = x_axis[1] - x_axis[0];
		double y_length = y_axis[1] - y_axis[0];
		
		// Bestimme Anzahl der Punkte 
		int number_of_points = poisson(f.max()*x_length*y_length);
		
		for(int i=0; i<number_of_points; i++)
		{
			// Bedingte Gleichverteilung
			double x = x_length * r.nextDouble() + x_axis[0];
			double y = y_length * r.nextDouble() + y_axis[0];
			
			// Unabhaengige Verduennung (Thinning)
			if(r.nextDouble()*f.max()<f.value(x,y))
			{
			   points.add(new double[]{x,y});
			}
		}
	    return points;
	}
	
	public static interface intensityFunction
	{
		public double value(double x, double y);
		public double max();
	}
	
	// Konstante Intensitaetsfunktion f(x,y)=lambda
	public static class HomogeneousIntensityFunction implements intensityFunction
	{
		double lambda = 0;
		
		public HomogeneousIntensityFunction(double l)
		{
			lambda = l;
		}
		
		public double value(double x, double y)
		{
			return lambda;
		}
		
		public double max()
		{
			return lambda;
		}
		
	}
	
	// Inhomgenene Intensitaetsfunktion der speziellen Form f(x,y)=lambda*x*y
	public static class InhomogeneousIntensityFunction implements intensityFunction
	{
		double lambda = 0;
		double max = 0;
		
		public InhomogeneousIntensityFunction(double l, double max_x, double max_y)
		{
			lambda = l;
			max = l * max_x * max_y;
		}
		
		public double value(double x, double y)
		{
			return lambda * x * y;
		}
		
		public double max()
		{
			return max;
		}
		
	}
	
	// Teststatistik T vom Uebungsblatt
	public static double calculateTestStatistic(Vector process, double window_volume, double hypothetical_value)
	{
	    double help = Math.sqrt(window_volume * window_volume / process.size()) * ((process.size() / window_volume) - hypothetical_value);
		return help;
	}
	
	// Randkorrigierter Schaetzer aus dem Skript 
	public static double edgeCorrectedEstimator(Vector P, double h, double x, double y, double x_max, double y_max)
	{
		int count = 0;
		
		
		// Zaehle alle Punkte, die in der Maximumsnorm naeher als h/2 an (x,y) und auch in W liegen. 
		// Da alle Punkte in P immer in W liegen, kann die Ueberpruefung,
		// ob sie in W liegen, entfallen.
		for(int i=0; i<P.size(); i++)
		{
			double x_new = ((double[]) P.get(i))[0];
			double y_new = ((double[]) P.get(i))[1];
			
			double distance = Math.max(Math.abs(x-x_new),Math.abs(y-y_new));
			
			if(distance<0.5*h) count++;
		}
		
		// Bestimme Grenzen und Volumen des Schnittes aus Beobachtungsfenster
		// und Schaetzbereich
		double x_upper_border = Math.min(x+0.5*h, x_max);
		double x_lower_border = Math.max(x-0.5*h, 0);
		
		double y_upper_border = Math.min(y+0.5*h, y_max);
		double y_lower_border = Math.max(y-0.5*h, 0);
		
		double volume = (x_upper_border - x_lower_border) * (y_upper_border -y_lower_border);
		
		
		return count / volume;
	}
	
	
	// Finde die aufsummierte LogLikelihood fuer die Punkte des Poisson-Prozesses 
	// bei einem randkorrigierten Schaetzer, der mit der Bandbreite h bestimmt wurde.
	public static double likelihoodCrossValidation(Vector P, double h, double x_max, double y_max)
	{
		double value = 0;
		for(int i=0; i<P.size(); i++)
		{
			// Lasse einen Punkt aus ...
			double[] leave_out = (double[]) P.remove(i);
			
			//... und bestimme seine LogLikelihood aus dem randkorrigierten Schaetzer mit Bandbreite h
			value += Math.log(edgeCorrectedEstimator(P,h,leave_out[0],leave_out[1],x_max,y_max));
			
			// Fuege den Punkt fuer die naechsten Durchlaeufe wieder ein.
			P.add(i,leave_out);
		}
		
		return value;
	}
	
    
	
public static void main(String[] args) {
		
		long sysBeg = System.currentTimeMillis();
	    try{
		
	    /*	
		{ //begin Aufgabe 1
			System.out.println("Blatt 3, Aufgabe 1a)");
			double lambda = 0.001;
			int n = 100;
			double[] m = {100,300,500};
			
			// Fuer jede der drei Fenstergroesen ...
			for(int m_index=0; m_index<m.length; m_index++)
			{
			
				FileOutputStream file = new FileOutputStream("C:\\Dokumente und Einstellungen\\Sebastian\\Eigene Dateien\\Lehre\\räumliche_statistik\\java_programs\\s_plus_commands\\commands"+m[m_index]+".txt"); 
			    BufferedWriter w = new BufferedWriter(new OutputStreamWriter(file));
				
			    //..werden 100...
				for(int i=0; i<n; i++)
			    {
					
					// ... Realisierungen des Poisson-Prozesses im Beobachtungsfenster erzeugt ...
					 Vector P=new Vector();
//					wenn der Poissonprozess keine Punkte im Beobachtungsfenster hat, ist die Testgröße nicht definiert
					 while(P.size()==0){ 
					  P = PoissonProcess(new double[]{0,m[m_index]}, new double[]{0,m[m_index]}, new RS_Blatt3.HomogeneousIntensityFunction(lambda));
					 }
				    
					//... und die zugehoerige Statistik T in eine Datei geschrieben ...
					w.write(""+calculateTestStatistic(P, m[m_index]*m[m_index], lambda));
				   
					
						w.write("\n");
				    	
				   }
				
				// ... gemeinsam mit den passenden Befehlen fuer den R - Aufruf.
				
				
				w.flush();
				w.close();
				file.close();
				
			}
			
			
				
		}  // End Aufgabe 1a)
		*/
		/*
     	{  System.out.println("Blatt 3, Aufgabe 1b) bzw. 1c)");
		
     	    // Fuer 10 verschiedene lambdas ...
     	    double[] lambda = {0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01};
			double[] rejections = new double[lambda.length];
     		double number_of_mc_tests = 100.0;
     		
     		// .. und eine Fenstergroesse von 100 X 100 (Aufgabe b) bzw. 1000 X 1000 (Aufgabe c) ...
     		double X_MAX = 100;
			
     		// ... fuehren wir 100 Mal den MC-Test durch ..
     		for(int mc_test=0; mc_test<number_of_mc_tests; mc_test++)
     		{
     			// ... indem wir fuer jedes lambda ...
				for(int k=0; k<lambda.length; k++)
				{ 
				   int m = 99;
				   double[] test_statistics = new double[m+1];
				
				   // ... 99 Realisierungen ...
				   for(int i=0; i<m; i++)
			       {
					   // ... eines entsprechenden Poisson-Prozesses erzeugen, ...
					   Vector P=new Vector();
					   while(P.size()==0){
				        P = PoissonProcess(new double[]{0,X_MAX}, new double[]{0,X_MAX}, new RS_Blatt3.HomogeneousIntensityFunction(0.001));
					   }
				        // ... die zugehoerige Teststatistik fuer den Vergleich mit lambda=0.001 berechnen ...
				       test_statistics[i] = Math.abs(calculateTestStatistic(P, X_MAX*X_MAX, 0.001));
				   }
				
				   // ... und dann noch eine Realisierung eines Poisson-Prozesses mit lambda=0.001 erzeugen,
				   // die wir als "die echten Daten" betrachten, ...
				   Vector P=new Vector();
				   while(P.size()==0){
				    P = PoissonProcess(new double[]{0,X_MAX}, new double[]{0,X_MAX}, new RS_Blatt3.HomogeneousIntensityFunction(lambda[k]));
				   }
				    // .. und ebenfalls die passende Statistik T berechnen.
				   double real_data_test_statistic = Math.abs(calculateTestStatistic(P, X_MAX*X_MAX, 0.001));
				   test_statistics[m] = real_data_test_statistic;
				
				   // Dann werden die 100 Statistiken T der Groesse nach geordnet ...
				   Arrays.sort(test_statistics);
				   
				   int position=-1;
				   
				   // ... und die Position der Statistik der "echten Daten" bestimmt.
				   for(int i=0; i<test_statistics.length; i++)
				   {
				      if(test_statistics[i]==real_data_test_statistic)
				      {
				    	  position=i;
				    	  break;
				      }
				   }  // End for i
				   
				   // Ist diese Position groesser als 95, lehnen wir die Hypothese ab.
				   if(position+1>95) rejections[k]++;
				   
				}  // End for k	
			} // End for mc_test;	
     		
     		for(int k=0; k<rejections.length; k++)
     		{
     			// Die empirische Guetefunktion an einer Stelle lambda ist dann einfach der Anteil der Tests,
     			// bei denen die Hypothese, dass lambda = lambda_0 (= 0.001) abgelehnt wurde.
     			System.out.println("Guetefunktion des MC-Tests fuer lambda="+lambda[k]+": "+(rejections[k]/number_of_mc_tests));   			
     		}
				
				
		}  // End Aufgabe 1b) bzw. 1c)
     	
     	*/
	    
     	{    System.out.println("Blatt 3, Aufgabe 2");
     		
     	     // Wir betrachten viele verschiedene Bandbreiten h.
     	     double[] h = {11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40};
  		    
     	     // Wir erzeugen 10 Mal ...
     	     for(int k=0; k<10; k++)
     		 {
     	       // ... eine Realisierung des inhomogenen Poisson-Prozesess... 
     		   RS_Blatt3.InhomogeneousIntensityFunction I = new RS_Blatt3.InhomogeneousIntensityFunction(0.00001,100,100);
     		   Vector P = PoissonProcess(new double[]{0,100}, new double[]{0,100}, I);
     		    
     		   // Aufgabe 2a)
     		   double maximum = Double.NEGATIVE_INFINITY;
     		   double h_max = -1;
     		   
     		   // ... und bestimmen dann fuer jedes der h ...
     		   for(int l=0; l<h.length; l++)
     		   {
     			   // ... die Summe der LogLikelihoods wie im Skript.
     			   double L = likelihoodCrossValidation(P, h[l], 100, 100);  
     			   //System.out.println("h="+h[l]+", L="+L);
     			   
     			   // Das h, bei dem die LogLikelihood maximal wird, ...
     			  if(L>maximum)
     			   {
     				   maximum = L;
     				   h_max = h[l];
     			   }
     		   }
     		   
     		   // .. ist dann die "optimale" Bandbreite.
     		   System.out.println("Nr. "+(k+1)+", Optimale Bandbreite: "+h_max);
     		   
     		    
     		   // Aufgabe 2b)
     		   
     		   // Genauso bestimmen wir fuer jedes h ...
     		   for(int l=0; l<h.length; l++)
     		   {
     			   double mean_square_dev = 0;
     			   for(int i=0; i<=99; i++)
	     		   {
	     			   for(int j=0; j<=99; j++)
	     			   {
	     				   // ... auch diese Annaeherung an die mittlere quadratische Abweichung...
	     				   mean_square_dev += Math.pow(edgeCorrectedEstimator(P, h[l], i+0.5,j+0.5,100,100) - I.value(i+0.5,j+0.5),2.0);
	     			   } // end for j
	     			   
	     			   //System.out.println(""+i);
	     		   } // End for i
     			   
     			   mean_square_dev = mean_square_dev / 10000.0;
     			   
     			   // ... die beim RANDKORRIGIERTEN Schaetzer fuer lambda natuerlich immer kleiner werden sollte,
     			   // wenn h groesser wird (vgl. hierzu Skript).
     			   System.out.println("Fuer Poisson-Process Nr."+(k+1)+" mit Bandbreite "+h[l]+" ist die mittlere quadr. Abweichung "+mean_square_dev);
     		   } // End for l
     		   
     			   
     		 }  // End for k
     		 
     		 
     		
     		
     	} // End Aufgabe 2
		
     	
     	
	    }catch(Exception e) { System.err.println("RSBlatt3.main(): "+e); }
	
     	long sysEnd = System.currentTimeMillis();
        int min = (int)((sysEnd - sysBeg) / 60000);
        int sek = (int)(((sysEnd - sysBeg)-60000*min) / 1000);
        System.out.println( "Zeit: " + min +" Min und " + sek + " Sekunden.");
        System.out.println( "Fertig!");
        System.exit(0);
	}

}