package RS06;

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




public class Blatt_5 {

	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;
		}
		
	}
	/*
	 * Diese Methode generiert eine randeffektfreie Realisierung eines Matern-Cluster-Prozesses auf einem 
	 * rechteckigen Beobachtungsfenster in R^2. Die Punkte
	 * werden in einem Vector zurückgegeben, der Arrays der Form double[2] enthält. Das rechteckige 
	 * Beobachtungsfenster wird über die double[2] x_axis und y_axis spezifiziert. lambda_0 und lambda_1 sind 
	 * die Intensitäten des Primär- bzw. Sekundär-Poisson-Prozesses. R ist der Radius der Kreise, auf denen die
	 * Sekundärprozesse realisiert werden.
	 *  
	 */
	public static Vector Matern_Cluster_Process(double[] x_axis, double[] y_axis, double lambda_0, double lambda_1, double R){
		
		Vector points=new Vector(); //enthält nachher die Punkte des Prozesses als double[2] 
		
		/*
		 * Die Realisierung des Matern-Cluster-Prozesses auf dem Fenster 
		 * [x_axis[0], x_axis[1]] x [y_axis[0], y_axis[1]] wird genau von all jenen Clustern beeinflusst, 
		 * deren Elternpunkte einen Abstand vonweniger als R vom Fenster haben. Deshalb muss man diese Elternpunkte
		 * und ihre Cluster zusätzlich zu denen im Fenster realisieren.
		 */
		double[] edge_corrected_x_axis= new double[]{x_axis[0]-R, x_axis[1]+R};
		double[] edge_corrected_y_axis= new double[]{y_axis[0]-R, y_axis[1]+R};
		
		//realisiere den Poisson Prozess auf dem vergrößerten Fenster
		Vector parent_points=PoissonProcess(edge_corrected_x_axis, edge_corrected_y_axis, new HomogeneousIntensityFunction(lambda_0) );
		
		
		double[] child_x_axis= new double[2];
		double[] child_y_axis= new double[2];
		
		for(int i=0; i<parent_points.size(); i++){
			double[] parent=(double[]) parent_points.get(i);
			/*
			 * definiere für jeden Elternpunkt ein quadratisches Fenster, das den 
			 * Kreis mit Radius R um den Elternpunkt enthält
			 */
			child_x_axis[0]=parent[0]-R;
			child_x_axis[1]=parent[0]+R;
			child_y_axis[0]=parent[1]-R;
			child_y_axis[1]=parent[1]+R;
			 
			//realisiere den Tochterprozess auf dem quadratischen Fenster (Akzeptanz- bzw. Verwerfung für den Kreis ist unten implementiert)
			Vector child_points=PoissonProcess(child_x_axis, child_y_axis, new HomogeneousIntensityFunction(lambda_1) );
			
			double[] child;
			
			/*
			 * Tochterpunkte sollen genau dann in den Ergebnisvektor geschrieben werden, wenn sie sowohl im Kreis
			 * um den Elternpunkt mit Radius R als auch im Beobachtungsfenster liegen.
			 */
			for(int j=0; j<child_points.size(); j++){
				child=(double[]) child_points.get(j);
				if(x_axis[0]<=child[0] && child[0]<=x_axis[1] && 
						y_axis[0]<=child[1] && child[1]<=y_axis[1] &&
						Math.sqrt(Math.pow(child[0]-parent[0], 2)+Math.pow(child[1]-parent[1], 2))<R){
					
					points.add(child);
				}
			}
			
		}
		
		return points;
	}
	
/*
 * Schätzt die nächste-Nachbar-Abstandsverteilungsfunktion eines Punktprozesses, der auf dem durch double[] x_axis und double[] y_axis 
 * spezifizierten
 *  Fenster realisiert wurde. Der Vektor process enthält die Punkte als double[2]. double[] r_values enthält die Argumente, für die
 *  die NNDDF geschätzt werden soll. Die Ergebnisse werden in einem double[] result zurückgegeben, so dass result[i] den Schätzer der NNDDF
 *  an der Stelle r_values[i] enthält. 
 */	

	public static double[] NNDDF(double[] r_values, Vector process, double[] x_axis, double[] y_axis){
	
		double[] result=new double[r_values.length];
		/* Zunächst wird für jeden Punkt des Prozesses sein Abstand zum Fensterrand ermittelt. Der Randabstand des i-ten Punktes
		 * wird in boundary_distances[i] gespeichert.
		 */
		double[] point;
		double[] boundary_distances=new double[process.size()];
		double x_dist;
		double y_dist;
		                               
		for(int i=0; i<boundary_distances.length; i++){
			
			point=(double[]) process.get(i);
			x_dist=Math.min(point[0]-x_axis[0], x_axis[1]-point[0]);
			y_dist=Math.min(point[1]-y_axis[0], y_axis[1]-point[1]);
			
			boundary_distances[i]=Math.min(x_dist, y_dist);
			
		}
		
		
		/*Bestimme für jeden Punkt seinen Abstand zum nächsten Nachbarn und speichere 
		 * das Resultat für den i-ten Punkt in next_neighbor_distances[i]
		 * Das ist der rechenaufwändigste Teil der Prozedur. Da die Abstände zu den nächsten 
		 * Nachbarn dann nur einmal berechnet werden müssen und dann zu Berechnung der Werte D(r) 
		 * für alle r verwendet werden können, ist es sinnvoll, einen ganzen Array 
		 * von Werten für r zu übergeben, an denen die NNDDF geschätzt wird.
		 */		
		double[] next_neighbor_distances=new double[process.size()];
		
		for(int i=0; i<next_neighbor_distances.length; i++){
			next_neighbor_distances[i]=Double.POSITIVE_INFINITY;
		}
		
		double dist=0.0;
		double[] point_i;
		double[] point_j;
		for(int i=0; i<next_neighbor_distances.length; i++){
			point_i=(double[]) process.get(i);
			for(int j=0; j<next_neighbor_distances.length; j++){
				if(i!=j){
				point_j=(double[]) process.get(j);
				
				dist=Math.sqrt(Math.pow(point_i[0]-point_j[0],2)+Math.pow(point_i[1]-point_j[1],2));
				
				next_neighbor_distances[i]=Math.min(next_neighbor_distances[i], dist );
				
				}
			}
		}
		
		
		double D=0.0, lambda_H=0.0;
		double x_length=x_axis[1]-x_axis[0];
		double y_length=y_axis[1]-y_axis[0];
		double next_neighbor_dist;
		double boundary_dist;
		
		/*
		 * Berechne lambda_H und den Schätzer D(r) für alle Werte r in r_values
		 */
		
		for(int j=0; j<r_values.length; j++){
		
		for(int i=0; i<process.size(); i++ ){
			next_neighbor_dist=next_neighbor_distances[i];
			boundary_dist=boundary_distances[i];
			
			if(next_neighbor_dist <boundary_dist && next_neighbor_dist<=r_values[j]){
				
				D+=1.0/(x_length-2.0*next_neighbor_dist)/(y_length-2.0*next_neighbor_dist);
			}
			if(next_neighbor_dist <boundary_dist ){
				lambda_H+=1.0/(x_length-2.0*next_neighbor_dist)/(y_length-2.0*next_neighbor_dist);
			}
		}
		
		if(lambda_H>0){
			
			result[j]=D/lambda_H;
			D=0;
			lambda_H=0;
		}
		else{
		result[j]= -1;
		D=0;
		lambda_H=0;
		}
		}
		return result;
	}
	
	
	/*
	 * Prozedur liest eine Liste von doubles, die in einer Textdatei gespeichert sind, in einen Array ein und 
	 * gibt diesen zurück 
	 */
	public static double[] get_from_file(String path){ 
		double[] dummy= new double[1];
		Vector objects=new Vector();
		try{
		FileReader r=new FileReader(path);
		BufferedReader in = new BufferedReader(r);
		
		String s;
		s=in.readLine();
		while(s!=null){
			objects.add(s);
			s=in.readLine();
		}
		double[] result= new double[objects.size()];
		for(int i=0; i<objects.size(); i++){
			result[i]=Double.parseDouble((String) objects.get(i));
		}
		
		r.close();
		return result;
		}
		catch (IOException e) {
		    System.out.println("File not found");
		  }
		return dummy;
	}
    
	
public static void main(String[] args) {
		
		long sysBeg = System.currentTimeMillis();
	    try{
		/*
	    	{//begin Aufgabe 1 (a)
	    		
	    		FileOutputStream file_x = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\x"+".txt"); 
			    BufferedWriter wx = new BufferedWriter(new OutputStreamWriter(file_x));
			    
			    FileOutputStream file_y = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\y"+".txt"); 
			    BufferedWriter wy = new BufferedWriter(new OutputStreamWriter(file_y));
	    		
			    double x_max=1000;
			    double lambda=0.001;
			    double lambda_0=0.001;
			    double R=50;
			    double lambda_1=lambda/(R*R*Math.PI*lambda_0);
			    
			    
			    double[] x_axis=new double[]{0, x_max};
			    double[] y_axis= new double[]{0, x_max};
			    
			    Vector points=Matern_Cluster_Process(x_axis, y_axis, lambda_0, lambda_1, R);
			    
			    double[] point;
			    for(int i=0; i<points.size(); i++) {
			    	point=(double[]) points.get(i);
			    	wx.write(point[0]+"\n");
			    	wy.write(point[1]+"\n");
			    }
			    wx.flush();
			    wy.flush();
			    wx.close();
			    wy.close();
			    
	    	}	
	    	
	    	*/
	    	/*
	    	{ //begin Aufgabe 1(c)
	    		
	    		FileOutputStream file_x = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\rvalues"+".txt"); 
			    BufferedWriter wx = new BufferedWriter(new OutputStreamWriter(file_x));
			    
			    FileOutputStream file_Dy = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\dvalues"+".txt"); 
			    BufferedWriter wDy = new BufferedWriter(new OutputStreamWriter(file_Dy));
			    
			    FileOutputStream file_theo_Dy = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\dtheoreticalpoisson"+".txt"); 
			    BufferedWriter w_theo_Dy = new BufferedWriter(new OutputStreamWriter(file_theo_Dy));
			    
			    double x_max=1000;
			    double lambda=0.001;
			    double lambda_0=0.001;
			    double R=50;
			    double lambda_1=lambda/(R*R*Math.PI*lambda_0);
			    
			    
			    double[] x_axis=new double[]{0, x_max};
			    double[] y_axis= new double[]{0, x_max};
			    
			    Vector points=Matern_Cluster_Process(x_axis, y_axis, lambda_0, lambda_1, R);
			     
			     double[] D;
			     double[] r_val=new double[50];
			     
			     for(int j=0; j<r_val.length; j++){
			    	 r_val[j]=(j+1);
			     }
			     
			     D=NNDDF(r_val, points, x_axis, y_axis);
			     
			    for(int i=1; i<=D.length; i++){
			    	
			    	
			    	wx.write(i+"\n");
			    	wDy.write(D[i-1]+"\n");
			    	w_theo_Dy.write((1-Math.exp(-lambda*Math.PI*i*i))+"\n");
			    	
			    	
			    }
			    wx.flush();
			    wDy.flush();
			    wx.close();
			    wDy.close();
			    w_theo_Dy.flush();
			    w_theo_Dy.close();
	    	
	    	}
	    	*/
	    	
	    	/*
	    	{ //Ein kleines zusätzliches Programm: Simulation D bei Poisson empirisch vs. theoretisch
	    		
	    		FileOutputStream file_x = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\r_values"+".txt"); 
			    BufferedWriter wx = new BufferedWriter(new OutputStreamWriter(file_x));
			    
			    FileOutputStream file_Dy = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\D_values"+".txt"); 
			    BufferedWriter wDy = new BufferedWriter(new OutputStreamWriter(file_Dy));
			    
			    FileOutputStream file_theo_Dy = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\theoretical_D_values"+".txt"); 
			    BufferedWriter w_theo_Dy = new BufferedWriter(new OutputStreamWriter(file_theo_Dy));
			    
			    double x_max=1000;
	    		double[] x_axis=new double[]{0, x_max};
				double[] y_axis= new double[]{0, x_max};
			    double lambda=0.001;
				
			    Vector points=PoissonProcess(x_axis, y_axis, new HomogeneousIntensityFunction(lambda) );
			     double r=1;
			     double[] D=new double[1];
			     
			     double[] r_val=new double[1];
			     
			     for(int i=1; i<=50; i++){
			    	r_val[0]=i*r;
			    	 D=NNDDF(r_val, points, x_axis, y_axis);
			    	
			    	wx.write(i*r+"\n");
			    	wDy.write(D[0]+"\n");
			    	w_theo_Dy.write((1-Math.exp(-lambda*Math.PI*i*r*i*r))+"\n");
			    	
			    	System.out.println("difference: "+ (D[0]-(1-Math.exp(-lambda*Math.PI*i*r*i*r))));
			    }
			    wx.flush();
			    wDy.flush();
			    wx.close();
			    wDy.close();
			    w_theo_Dy.flush();
			    w_theo_Dy.close();
	    	
	    	}
	    	
	    	*/
	    	
	    	/*
	    	{	// Teil (d)
	    		
	    		 //Berechnung eines punktweisen empirischen 94% Konfidenzintervalls der NNDDAF für einen Poisson-Prozess
	    		 //mit Intensität lambda=0.001 für r=5,10,15,...,50
	    		 
	    		
	    		FileOutputStream file_x = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\r_values_(d)"+".txt"); 
			    BufferedWriter wr = new BufferedWriter(new OutputStreamWriter(file_x));
			    
			    FileOutputStream file_lower = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\lower"+".txt"); 
			    BufferedWriter low= new BufferedWriter(new OutputStreamWriter(file_lower));
			    
			    FileOutputStream file_upper = new FileOutputStream("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\upper"+".txt"); 
			    BufferedWriter up= new BufferedWriter(new OutputStreamWriter(file_upper));
			    
	    		
			    double x_max=1000;
	    		double[] x_axis=new double[]{0, x_max};
				double[] y_axis= new double[]{0, x_max};
			    double lambda=0.001;
	    		
	    		
	    		
	    		double[] r_values=new double[]{5,10,15,20,25,30,35,40,45,50};
	    		//Die Ergebnisse der NNDDF sollen an allen Stellen für r für jede Realisierung gespeichert werden
	    		double[][] results= new double[r_values.length][200];
	    		double[] specific_result;
	    		
	    		//NNDDF soll für 200 Realisierungen des Poisson Prozesses generiert werden
	    		
	    		for(int j=0; j<200; j++){
    			System.out.println(j);
   				 Vector points=PoissonProcess(x_axis, y_axis, new HomogeneousIntensityFunction(lambda) );
   				 specific_result =NNDDF(r_values, points, x_axis, y_axis);
   				 
   				 for(int i=0; i<r_values.length; i++){
   					 
   					 results[i][j]=specific_result[i];
   					 
   				 	}
   			
	    		}
	    		
	    		double r=5.0;
	    		// Jetzt müssen an allen Stellen r die Ergebnisse sortiert werden und die 6 kleinsten und die 6 größten Werte
	    		//werden gestrichen, d.h. der siebt-kleinste bzw. siebt-größte Wert der NNDDAF bilden das Konfidenzintervall
	    		 for(int i=0; i<r_values.length; i++){
	    			Arrays.sort(results[i]);
	    			
	    			System.out.println();
	    			wr.write((i+1)*r+"\n");
		    		
	    			low.write(results[i][6]+"\n");
		    		
		    		up.write(results[i][193]+"\n");
		    		
		    		System.out.println((r*(i+1))+", "+results[i][6]+", "+results[i][193]+", "+(1-Math.exp(-lambda*Math.PI*(i+1)*r*(i+1)*r)));
	    			 
	    			 
	    		 }
	    		
	    	
	    		wr.flush();
	    		wr.close();
	    		low.flush();
	    		up.flush();
	    		low.close();
	    		up.close();
	    	}
	    	
	   	
	    */
	    	{// Teil (e)
	    		
	    		double[] upper=get_from_file("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\upper.txt");
	    		double[] lower=get_from_file("C:\\Dokumente und Einstellungen\\Gast\\Eigene Dateien\\Lehre\\Sebastian\\data_blatt_5\\lower.txt");
	    	
	    	double x_max=1000;
	    	
	    	
	    	double intensity=0.001;
	    	double lambda_0=0.001;
	    	double R=50;
	    	double lambda_1=intensity/(R*R*Math.PI*lambda_0);
	    	 
	    	double[] x_axis=new double[]{0, x_max};
			double[] y_axis= new double[]{0, x_max};
	    	
			boolean contained;
			
			double[] D;
			
			int number_of_simulations=100;
			int count=0;
	    	
			double[] r_val=new double[10];
			
			for(int j=0; j<=9; j++){
				r_val[j]=5.0*(j+1);
			}
			
			for(int i=0; i<number_of_simulations; i++){
				contained=true;
				Vector points=Matern_Cluster_Process(x_axis, y_axis, lambda_0, lambda_1, R);
	    		
				
				
				 D=NNDDF(r_val, points, x_axis, y_axis);
				
				 for(int j=1; j<=10; j++){  
				   
				    if(lower[j-1]>D[j-1] || upper[j-1]<D[j-1]){
				    	contained=false;
				    
				    System.out.println("Realisierung "+i+" nicht im Schlauch bei "+(5*j));
				    break;
				    }
	    		}
				if(contained){
					count++;
					
				}
	    	}
			
			System.out.println("lambda_0="+lambda_0+": Wkt für NNAVF im Schlauch: " +(((double) count)/((double) number_of_simulations)));
	    }
	    	
	    
	    
	    
     	
     	
	    }catch(Exception e) { System.err.println("Blatt_5.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);
	}

}