package ole.raeumlicheStatistik;

import java.util.Arrays;
import java.util.Random;
import java.util.Vector;


public class RSBlatt6 {

	static final Geometry2D.Point zero = new Geometry2D.Point(0,0);
	static final Random r = new Random();
	
	// Fuer die radiale Simulation des Poisson-
	// Prozesses in Aufgabe 1
	double radial_distance=0;
	double intensity=-1;
	
	// Fuer die simulation der typischen 
	// Zelle in Aufgabe 1
	Vector points = new Vector(3,1);
	Geometry2D.ConvexPolygon initialCell;
	Geometry2D.ConvexPolygon typicalCell;
	
	public RSBlatt6(double lambda)
	{
		radial_distance=0;
		intensity = lambda;
	}
		
	
	/* *
	 *  Aufgabe 1                              *
	 ***/	
		
	// Bestimme die Polar-Koordinaten des naechsten Punktes 
	// in der radialen Simulation eines Poisson-Prozesses 
	// mit der im Konstruktor festgelegten Intensitaet
	// (Algorithmus siehe Skript)
	public double[] nextRadialPoissonPoint()
	{
		
		double e = -1.0 * Math.log(r.nextDouble()) / intensity;
		double theta = 2 * Math.PI * r.nextDouble();
		radial_distance += e ;
		
		return new double[] {Math.sqrt(radial_distance/Math.PI), theta};
	}

	// Bestimme die Punkte eines Poisson-Prozesses, die noetig sind 
	// um die initiale Zelle zu bestimmen, die bei der Simulation der 
	// typischen Zelle noetig sind.
	public Geometry2D.Point[] simulateInitialCellPoints() throws Exception
	{
		
	try{
		
		double[] angles = new double[2];
		
		// Bestimme zunaechst die 3 Punkte, die 
		// mindestens noetig sind.
		double[] point = nextRadialPoissonPoint();
		points.add(point);	
		angles[0] = point[1];
		
		point = nextRadialPoissonPoint();
		points.add(point);	
		angles[1] = point[1];
		
		point = nextRadialPoissonPoint();
		points.add(point);	
			
		while(true)
		{  /* So geht's auf jeden Fall:
			boolean[] quads = {false, false, false, false};
			
			for(int i=0; i<points.size(); i++)
			{
				double[] point = (double[]) points.get(i);
				if((point[1]>0) && (point[1]<=Math.PI / 2.0)) { quads[0]=true; continue;}
				if((point[1]>Math.PI / 2.0) && (point[1]<=Math.PI)) { quads[1]=true; continue;}
				if((point[1]>Math.PI) && (point[1]<=1.5 * Math.PI )) { quads[2]=true; continue;}
				if((point[1]>1.5 * Math.PI) && (point[1]<=2 * Math.PI)) { quads[3]=true; continue;}
			}
			
			if(quads[0] && quads[1] && quads[2] && quads[3]) 
			{ 
				return PolarToG2D(points);	
			}
			else
			{
				double[] point = nextRadialPoissonPoint();
			    points.add(points.size(),point); 
			}
			*/
			
			Arrays.sort(angles);
			
			// Ob ein weiterer Punkt benoetigt wird, haengt davon ab, wo 
			// der letzte Punkt liegt, d.h. insbesondere davon, wie sein
			// Winkel zur x-Achse ist.
			double[] last_added = (double[]) points.lastElement();
			double angle = last_added[1];
			
			// Dafuer gibt es 3 Moeglichkeiten:
			// a) Der Punkt liegt zwischen den bisherigen Punkten. Damit liegt er 
			//    nicht im fraglichen Sektor, veraendert den fraglichen Sektor aber 
			//    auch nicht. Es muss also ein weiterer Punkt simuliert werden. 
			if  
			(  ((angles[1]-angles[0]>Math.PI) && ((angle>angles[1]) || (angle<angles[0])))
			   ||((angles[1]-angles[0]<=Math.PI) && ((angle<angles[1]) && (angle>angles[0]))))
			{
				points.add(nextRadialPoissonPoint());
				continue;
			}
			
			// b) Der Punkt liegt im fraglichen Sektor. Damit ist das Abbruchkriterium
			//    erfuellt und es sind genug Punkte zur Bestimmung der initialen Zelle 
			//    simuliert. 
			double[] anti_angles = {angles[0], angles[1]};
			anti_angles[0] += (angles[0]>Math.PI)?(-Math.PI):(Math.PI);
			anti_angles[1] += (angles[1]>Math.PI)?(-Math.PI):(Math.PI);
			Arrays.sort(anti_angles);
			
			if
			(  ((anti_angles[1]-anti_angles[0]>Math.PI) && ((angle>anti_angles[1]) || (angle<anti_angles[0])))
			   ||((anti_angles[1]-anti_angles[0]<=Math.PI) && ((angle<anti_angles[1]) && (angle>anti_angles[0]))))
			{
				return PolarToG2D(points);	
			}
			
			// c) Ansonsten liegt der Punkt nicht im fraglichen Sektor, aber auch nicht 
			//    zwischen den bisherigen Punkten. Also muss der fragliche Sektor durch
			//    die zugehoerige Gerade erweitert werden und dann ein weiterer Punkt
			//    simuliert werden.
			double[] diff_angles = {Math.abs(angle-angles[0]), Math.abs(angle-angles[1])};
			diff_angles[0] += (diff_angles[0]>Math.PI)?(-Math.PI):(Math.PI);
			diff_angles[1] += (diff_angles[1]>Math.PI)?(-Math.PI):(Math.PI);
			
			if(diff_angles[0]<diff_angles[1])	
			{
				angles[0] = angle;
			}
			else
			{
				angles[1] = angle;
			}
			points.add(nextRadialPoissonPoint());
			continue;	
		}
	}
	catch(Exception e){ System.err.println("RSBlatt7.simulateInitialCellPoints: "+e); throw e;}
	}

	// Bestimme aus den uebergebenen Punkten die initiale Zelle 
	// fuer die Simulation der typischen Zelle 
	public void createInitialCell(Geometry2D.Point[] points)
	{
	try{
		// Bestimme dazu alle Bisektoren der Verbindungsstrecken
		// zwischen dem Null-Punkt und allen simulierten Punkten.
		Vector lines = new Vector(points.length*points.length,1);
		for(int i=0; i<points.length; i++)
		{
			Geometry2D.Line b = Geometry2D.bisector(zero,points[i]);
			lines.add(b);
		}
		
	    // Dann bestimme alle deren Schnittpunkte als 
		// theoretisch moegliche Ecken der initialen Zelle.
		double x_max = Double.NEGATIVE_INFINITY;
		double x_min = Double.POSITIVE_INFINITY;
		double y_max = Double.NEGATIVE_INFINITY;
		double y_min = Double.POSITIVE_INFINITY;
		
		for(int i=0; i<points.length; i++)
		{
			Geometry2D.Line line_i = (Geometry2D.Line) lines.get(i);
			
			for(int j=i+1; j<points.length; j++)
			{
				Geometry2D.Line line_j = (Geometry2D.Line) lines.get(j);
				Geometry2D.Point intersect = (Geometry2D.Point) Geometry2D.intersect(line_i, line_j);
				
				if(intersect.x>x_max) x_max = intersect.x;
				if(intersect.x<x_min) x_min = intersect.x;
				if(intersect.y>y_max) y_max = intersect.y;
				if(intersect.y<y_min) y_min = intersect.y;
			}
		}
		
		// Jeweils die minimalen bzw. maximalen x- bzw. y-Werte der 
		// Schnittpunkte der Bisektoren koennen als Rand eines 
		// Ausgangsfensters w genommen werden, in dem die initiale Zelle 
		// garantiert enthalten ist. 
		Geometry2D.Point a = new Geometry2D.Point(x_min-1, y_min-1);
		Geometry2D.Point b = new Geometry2D.Point(x_max+1, y_min-1);
		Geometry2D.Point c = new Geometry2D.Point(x_max+1, y_max+1);
		Geometry2D.Point d = new Geometry2D.Point(x_min-1, y_max+1);
		Geometry2D.ConvexPolygon w = new Geometry2D.ConvexPolygon(new Geometry2D.Point[] {a,b,c,d});
		
		// Schneide dieses Ausgangsfenster mit den Halbebenen, die jeweils durch 
		// die Bisektoren der Verbindungsstrecken zwischen Null-Punkt und 
		// den Poisson-Punkten definiert werden.
		for(int i=0; i<lines.size(); i++)
		{
		   w = Geometry2D.intersectWithHalfplane(w, ((Geometry2D.Line) lines.get(i)) );
		}
		
		// So erhaelt man die initiale Zelle.
		initialCell = w;
		
	}
	catch(Exception e){ System.err.println("RSBlatt7.createInitialCell: "+e); }

	}

	// Bestimme aus der bereits vorliegenden initialen Zelle 
	// die typische Zelle einer PVT
	public void createTypicalCell()
	{
		
	try{
		// Dazu finde die am weitesten vom Null-Punkt entfernte
		// Ecke der initialen Zelle 
		typicalCell = new Geometry2D.ConvexPolygon(initialCell);
		double r_max = findMaxCornerDistance(typicalCell);
		
		// Nur, wenn diese Entfernung mehr als halb so gross ist 
		// wie die Entfernung des naechsten simulierten Punktes 
		// vom Null-Punkt, kann dieser Punkt die Zelle evtl. 
		// beeinflussen und die Zelle muss mit der durch den 
		// zugehoerigen Bisektor definierten Halbebene 
		// geschnitten werden.
		double[] p = nextRadialPoissonPoint();
		
		while( p[0]  <= 2*r_max)
		{
			Geometry2D.Point help = new Geometry2D.Point(p[0] * Math.cos(p[1]), p[0] * Math.sin(p[1]));
			Geometry2D.Line b = Geometry2D.bisector(zero,help);
			typicalCell = Geometry2D.intersectWithHalfplane(typicalCell, b );
			points.add(p);
			r_max = findMaxCornerDistance(typicalCell);
			p = nextRadialPoissonPoint();
		}
		
	}
	catch(Exception e){ System.err.println("RSBlatt7.createTypicalCell: "+e); }
	  	
	}
		
	public double[] getTypicalCellCharacteristics()
	{
		double[] return_values = {0,0,0};
		
		return_values[0] = typicalCell.p.length;
		return_values[1] = typicalCell.getPerimeter();
		return_values[2] = typicalCell.getArea();
		
		return return_values;
	}

	//Rechne einen Vector von Punkten, gegeben als double[]
	//mit Polar-Koordinaten um in ein Geometry2D.Point- Array
	public static Geometry2D.Point[] PolarToG2D(Vector points)
	{
		Geometry2D.Point[] return_value = new Geometry2D.Point[points.size()];
		for(int i=0; i<points.size(); i++)
		{
			double[] p = (double[]) points.get(i);
			return_value[i] = new Geometry2D.Point(p[0] * Math.cos(p[1]), p[0] * Math.sin(p[1]));
		}

	return return_value;
	}

	//Bestimme den maximalen Abstand einer Ecke 
	//der uebergebenen Zelle vom Null-Punkt
	public static double findMaxCornerDistance(Geometry2D.ConvexPolygon poly)
	{
		double max = Double.NEGATIVE_INFINITY;

		for(int i=0; i<poly.p.length; i++)
		{
		    double d = poly.p[i].distanceTo(zero);
			if(d>max) max = d;
		}

	  return max;
	}
		
//	 Fuehre einen asymptotischen Test der Hypothese durch, dass eine i.i.d. Zufallsstichprobe 
//	 vom Umfang n mit Stichprobenmittel mean und Standardabweichung sd als Realisierungen
//	 einer Zufallsgroesse mit Erwartungswert theoretical_parameter aufgefasst werden kann
	public static void AsymptoticTestForExpectation(double theoretical_parameter, double mean, double sd, int n)
	{
		double T = Math.abs(Math.sqrt(n) * (mean - theoretical_parameter) / sd );
		
		System.out.println("Test-Statistic T="+T);
		if(T > 1.65) // 0.95 - Quantil von N(0.1)
		{
			System.out.println("Hypothesis rejected at the 10%-Level.");
			if(T > 1.96)  // 0.975 - Quantil von N(0.1)
			{
				System.out.println("Hypothesis rejected at the  5%-Level.");
				if(T > 2.58) // 0.995 - Quantil von N(0,1)
				{
					System.out.println("Hypothesis rejected at the  1%-Level.");
				}
				else
				{
					System.out.println("Hypothesis not rejected at the  1%-Level.");
				}
			}
			else
			{
				System.out.println("Hypothesis not rejected at the  5%-Level.");
				System.out.println("Hypothesis not rejected at the  1%-Level.");
			}
		}
		else
		{
			System.out.println("Hypothesis not rejected at the 10%-Level.");
			System.out.println("Hypothesis not rejected at the  5%-Level.");
			System.out.println("Hypothesis not rejected at the  1%-Level.");
		}
	}
	
	public static double[] Aufgabe1cUnd1d()
	{
		double[] mean_values = {0,0,0};
		double[] sum_of_squares = {0,0,0};
		double[] sd = {0,0,0};

	try{
		   // Fuehre die vorgegebene Anzahl von Durchlaeufen 
		   // mit der vorgegebenen Intensitaet durch.
		   final int number_of_runs = 10000;
		   final double intensity = 0.0005;
				   
		   for(int i=0; i<number_of_runs; i++)
		   {
			      // Simuliere die Punkte fuer die initiale Zelle...
			      // ... bestimme daraus die initiale Zelle und ...
			      // und daraus wiederum die typische Zelle.
				  RSBlatt6 r = new RSBlatt6(intensity); 
				  r.createInitialCell(r.simulateInitialCellPoints());
				  r.createTypicalCell();
				   
				  // Bestimme die Charakteristiken dieser Zelle...
				  double[] result = r.getTypicalCellCharacteristics();
				  
				  // ... und daraus Mittelwert und Standardabweichung. 
				  for(int j=0; j<result.length; j++)
				  {
					  mean_values[j] += result[j];
					  sum_of_squares[j] += result[j] * result[j];
				  }
				  
				  System.gc();
		 }
			
		 for(int j=0; j<mean_values.length; j++)
		 {
		    mean_values[j] /= number_of_runs;
			sd[j] = Math.sqrt((sum_of_squares[j] - number_of_runs * mean_values[j] * mean_values[j])/(number_of_runs-1));
		 }
			 
		 // Fuehre die Tests auf die theoretischen Werte fuer
		 // die typische Zelle einer PVT durch.
		 System.out.println("Mean Number Of Vertices: "+mean_values[0]);
		 System.out.println("Test for Expectation=6:");
		 AsymptoticTestForExpectation(6, mean_values[0], sd[0], number_of_runs);
		 System.out.println("--------------------------------------------------");
		  
	     System.out.println("Mean Perimeter: "+mean_values[1]);
		 System.out.println("Test for Expectation="+(4 / Math.sqrt(intensity))+":");
		 AsymptoticTestForExpectation(4 / Math.sqrt(intensity), mean_values[1], sd[1], number_of_runs);
		 System.out.println("--------------------------------------------------");
		 
		 System.out.println("Mean Area: "+mean_values[2]);
		 System.out.println("Test for Expectation="+(1 / intensity)+":");
		 AsymptoticTestForExpectation((1 / intensity), mean_values[2], sd[2], number_of_runs);
		 System.out.println("--------------------------------------------------");
	}
	catch(Exception e){ System.err.println("RSBlatt7.Aufgabe2d: "+e); } 
		   
	return mean_values;
	}

	public static void main(String[] args)
	{
	   long sysBeg = System.currentTimeMillis();

	   Aufgabe1cUnd1d();
	    
	   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);
	}

	/* Resultate:
	 
	Mean Number Of Vertices: 5.9923
	Test for Expectation=6:
	Test-Statistic T=0.5793540724210869
	Hypothesis not rejected at the 10%-Level.
	Hypothesis not rejected at the  5%-Level.
	Hypothesis not rejected at the  1%-Level.
	--------------------------------------------------
	Mean Perimeter: 178.72775870151793
	Test for Expectation=178.88543819998318:
	Test-Statistic T=0.3628134439724059
	Hypothesis not rejected at the 10%-Level.
	Hypothesis not rejected at the  5%-Level.
	Hypothesis not rejected at the  1%-Level.
	--------------------------------------------------
	Mean Area: 1998.547350009426
	Test for Expectation=2000.0:
	Test-Statistic T=0.13761627786083533
	Hypothesis not rejected at the 10%-Level.
	Hypothesis not rejected at the  5%-Level.
	Hypothesis not rejected at the  1%-Level.
	--------------------------------------------------
	Zeit: 38 Min und 54 Sekunden.
	Fertig!
	*/
}