package WT;

import java.util.Random;
import java.awt.event.*;
import java.awt.*;
import javax.swing.JFrame;
import javax.swing.JPanel;
import java.io.*;

import WT.tools.XYPlot;
import WT.tools.Utilities;

public class Wiener1 {


   /** 
    * Aufruf und Einstellung von XYPlot
    *
    * @param	xyPairsToPlot	arrays von xy-Koordinaten, die als 
    *                           Stuetzstellen/-werte des Funktionsplots dienen
    * @param    title	        Titel des Ausgabefensters
    */
    
    private static void plot(double[][] xyPairsToPlot, int[] numOfPairs, String title){
        
	XYPlot plotter;
        plotter = new XYPlot(xyPairsToPlot, numOfPairs);

	double minX = 0.0;
        double maxX = 1.0;
	double minY = Utilities.getMin(xyPairsToPlot[0],2);
	double maxY = Utilities.getMax(xyPairsToPlot[0],2);
	for(int i = 1; i < xyPairsToPlot.length; i++){
	    minY = Math.min(minY,Utilities.getMin(xyPairsToPlot[i],2));
	    maxY = Math.max(maxY,Utilities.getMax(xyPairsToPlot[i],2));
	}
	plotter.setMinX(minX);	plotter.setMaxX(maxX);
	plotter.setMinY(minY);	plotter.setMaxY(maxY);

        /* die x- und y-Koordinaten werden gestreckt (bzw. gestaucht), so dass
	   die Darstellung im 600 x 600 Pixel Bild ausreichend groß ist. */
	
        double scaleX = 700.0/(maxX-minX);
	double scaleY = 400.0/(maxY-minY);
	plotter.setScaleX(scaleX); plotter.setScaleY(scaleY);

	JFrame f = new JFrame(title);
	f.addWindowListener(new WindowAdapter() {
	  public void windowClosing(WindowEvent e){ System.exit(0); }
	});		
	f.getContentPane().add(plotter);
	f.pack();
	f.setSize(800,500);
	f.setVisible(true);
    }
    
/**
 *  Bestimmt den Wert der Schauder-Funktionen S_n(t)
 *  @param t  betrachtete Stelle
 *  @param n  Ordnung der Schauder-Funktion
 *
 *  @return Funtkionswert der Schauderfunktion n-ter Ordnung an Stelle t
 */
    private static double schauderFunction(double t, int n){
	if (n==1) { return t;}

	/* Bestimme die Darstellung n=2^m+k mit k=1,...,2^m */
	int m = 0, k;
	double c = Math.pow(2,m);
	while (2*c < n){
	    m++;
	    c = Math.pow(2,m);
	}
	k = (int)(n - c);

        /*  Funktionswert wie in Vorlesung bzw. Musterlösung beschrieben */
	if( t < (k-1)/c || t >= k/c ){ 
	    return 0; }
	
	if( t < ((2*k)-1)/(2.0*c) ){
	    return ( Math.sqrt(c)*(t - (k-1)/c) ); }
	else{
	    return (-Math.sqrt(c)*(t - k/c) );     }
    }

/**
 *  Verwendung des ersten Algorithmus zur Approximation eines Wiener Prozesses X
 *  @param m  2^m = Ordnung der Approximation
 *  @param seed  Startwert für Zufallszahlengenerator
 *  @return Realisierungen von X an den Stützstellen
 */
    private static double[] algorithm1(int m, long seed){
	double n_d = Math.pow(2,m);
	int n = (int)n_d, k=0;
	double[] result = new double[2*(n+1)];
	/* Start in (0,0) */
	result[0] = 0;
	result[1] = 0;

	/* erzeuge n=2^m N(0,1)-Zufallszahlen */
	Random rand = new Random(seed);
	double[] y = new double[n+1];
	for(k = 0; k < n; k++){
	    y[k] = rand.nextGaussian();
	}

	/* Berechnung der Funktionswerte an den Stützstellen */
	double sum;
	for(int i = 1; i <= n; i++){
	    sum = 0;
	    result[2*i] = i/n_d;
	    for(k = 0; k < n; k++){
		sum = sum + y[k]*schauderFunction(i/n_d, k+1);
	    }
	    result[2*i+1] = sum;
	}
	return result;
    }
 

/**
 *  Verwendung des zweiten Algorithmus zur Approximation eines Wiener Prozesses X
 *  @param m  2^m = Ordnung der Approximation
 *  @param seed  Startwert für Zufallszahlengenerator
 *  @return Realisierungen von X an den Stützstellen
 */
    private static double[] algorithm2(int m, long seed){
	double n_d = Math.pow(2,m);
	int n = (int)n_d;
	double[] result = new double[2*(n+1)];

	/* Start in (0,0) */
	result[0] = 0;
	result[1] = 0;

	/* Berechnung der Funktionswerte an den Stützstellen 
	   unter Verwendung von U({-1,1})-Zufallszahlen */
	double sum = 0;
	double z = 0;
	Random rand = new Random(seed);
	for(int i = 1; i <= n; i++){
	    result[2*i] = i/n_d;
	    z = rand.nextDouble();
	    if( z < 0.5){ sum = sum - 1.0/Math.sqrt(n_d);  }
	    else{ 	  sum = sum + 1.0/Math.sqrt(n_d);  }
	    result[2*i+1] = sum;
	}
	return result;
    }
    
    /*
    public static void export(double[][] xyPairs, int[] numOfPairs, String filename)
    {
        for(int i=0;i<xyPairs.length;i++)
        {
            try{
                PrintWriter writer = new PrintWriter(new FileWriter(new File(filename+i+".txt")));
                for(int j=0;j<xyPairs[i].length;j++)
                {
                    System.out.println(j);
                    writer.println((j/(numOfPairs[i]-1)) +" "+ xyPairs[i][j]);
                }
                writer.close();
            }
            catch(IOException e)
            {System.out.println(e.getMessage());};
        }
    }
    */
/**
 *   Hauptprogramm:
 *   Simulation eines Wiener Prozesses mit n = 2^6, 2^7 bzw. 2^8 Stützstellen.
 *   Approximation der Ordnung n = 2^6, 2^7 bzw. 2^8.
 */
    public static void main(String[] args) {
	int runs = 3, m=8;
	double[][] xyPairs = new double[runs][2*(int)Math.pow(2,m)+2];
	Random rand_seed = new Random(System.currentTimeMillis());

	/* Algorithmus 1 */
	int[] numOfPairs = new int[]{(int)Math.pow(2,6)+1,(int)Math.pow(2,7)+1,(int)Math.pow(2,8)+1 };
  	/* Aufruf von Random() in Methode algorithmus1() bzw. algorithmus2() 
	 * Unabhängige Simulationen werden über Initialisierung mit verschiedenen random seeds erzielt
 	 * Zur Veranschaulichung auch seed = fester Wert denkbar */
	xyPairs[0] = algorithm1(6,rand_seed.nextLong()); // oder z.B. xyPairs[0] = algorithm1(4,2);
	xyPairs[1] = algorithm1(7,rand_seed.nextLong()); // oder z.B. xyPairs[1] = algorithm1(7,2);
	xyPairs[2] = algorithm1(8,rand_seed.nextLong()); // oder z.B. xyPairs[2] = algorithm1(8,2);

	/* Anzeigen des simulierten Pfades am Bildschrim */
	plot(xyPairs, numOfPairs, "Algorithmus 1");
        //export(xyPairs, numOfPairs, "C:/java/packages/AlgorithmusA");
        
        
	/* Algorithmus 2 */
	xyPairs[0] = algorithm2(6,rand_seed.nextLong()); // oder z.B. xyPairs[0] = algorithm2(4,2);
	xyPairs[1] = algorithm2(7,rand_seed.nextLong()); // oder z.B. xyPairs[1] = algorithm2(7,2);
	xyPairs[2] = algorithm2(8,rand_seed.nextLong()); // oder z.B. xyPairs[2] = algorithm2(8,2);

	/* Anzeigen des simulierten Pfades am Bildschrim */
	plot(xyPairs, numOfPairs, "Algorithmus 2");
	
        
        
    }
    
    
      
}