package WT;

/*******************************************************************
 *       Simulation eines zeitstetigen Markov-Prozesses            *
 *       mit Hilfe von exponentialverteilten Zwischenankunfts-     *
 *       zeiten und der Simulation der eingebetteten zeit-         *
 *       diskreten Markov-Kette.                                   *
 *                                                                 *
 *   siehe Abschnitt 2.3.4 der Vorlesung sowie Uebungsaufgabe 3.7. *
 *******************************************************************/

import java.util.Random;
import java.text.DecimalFormat;

public class Markov{

    /* Formatierung der Ausgabe */
    static DecimalFormat df = new DecimalFormat ("#0.000");
    static public String out(double y){	return df.format(y);  }

   /*---------------------------------------------------------------------------*
                        Variablen Deklaration
    *---------------------------------------------------------------------------*/

    /* kumulierte Uebergangswahrscheinlichkeiten plus Extra-0-Spalte als Hilfe  */
    static private double[][] p_add = null;

     /*  Intensitaeten q_i = -q_ii  */
    static private double[] q = null;   

    /* kumulierte Anfangsverteilung plus a[-1] = 0.0 als Hilfsgroesse */
    static private double[] a_add = null;   

    /* Initialisiert ein neues Markov-Objekt mit     
     * @param intMatrix_Q Intensitaetsmatrix Q 
     * @param a           Anfangsverteilung     */
    public Markov(double[][] intMatrix_Q, double[] a){
	int dim = a.length;

	/* Dimension der Anfangsverteilung und der Intensitaetsmartix 
	 * muessen uebereinstimmen */
	if(intMatrix_Q == null || 
	   intMatrix_Q.length != dim || intMatrix_Q[0].length != dim){
	    System.out.println("Fehler in der Intensitaetsmatrix Q");
	    System.exit(0);
	}

	/* Definition der Intensitaeten q[i], aufsummierten Werte der 
	 * Anfangsverteilung und der Uebergangswahrscheinlichkeiten 
	 * der eingebetteten Markov-Kette */
	this.a_add = new double[dim+1];
	this.q = new double[dim];
	this.p_add = new double[dim+1][dim+1];
	for(int i = 0; i<dim; i++){
	    this.q[i] = -intMatrix_Q[i][i];

	    this.a_add[i+1]  = this.a_add[i]+ a[i];

	    for(int j = 0; j< dim; j++){
		if(i != j && q[i] > 0){   
		    this.p_add[i][j+1] = this.p_add[i][j] 
			                 + intMatrix_Q[i][j]/this.q[i];
		}
		else{ 
		    this.p_add[i][j+1] = this.p_add[i][j];
		}
	    }
	}

    }


   /*------------------------- Transformation: ---------------------------------*
    *	Bestimmung des neuen Zustandes "new_state" der Markov-Kette 
    *   gegeben X_n = "current_state"
    *
    * phi(i,z) = sum_{k ¤ E} k * 1(sum_{j=1...k-1} p_ij < z <= sum_{j=1...k} p_ij) 
    ----------------------------------------------------------------------------*/
    static private int transform(int current_state, double z){
	int new_state = 0; 
	int i = current_state-1; /* "-1" nur wegen Indizierung des Array! */

	for(int k=1; k <= 3; k++){
	    if(p_add[i][k-1] < z && z <= p_add[i][k]){  
		new_state = k;  break;  
	    }
	}
	
	return new_state;
    }


   /*---------------------------------------------------------------------------*
                            Hauptprogramm                      
    *---------------------------------------------------------------------------*/
    static public void main(String[] args){
	
	if(args.length < 2){
	    System.out.println("Aufruf mit Stichprobenumfang Maximalzeit"); 
	    System.exit(0);	}

	/* args[0] = Anzahl der Simulationen, x = aktueller Zustand der MK */
	int sample_size = Integer.parseInt(args[0]), x = 0;
	/* args[1] = Zeit bis Abbruch einer einzelnen Simulation */
	double max_time = Double.parseDouble(args[1]);


	/* Hilfsvariablen:  z = Zustand,  t = Zeitpunkt */
	double z = 0, t = 0;

	/* Zur Bestimmung der relativen Häufigkeiten der Zustände 
	   bei Simulationsabbruch */
	double[] freq = new double[3];
	double add = 1.0/(double)sample_size;

	/* Intensitaetsmatrix Q */           
        double[][] intMartix_Q = new double[][]{
	    new double[]{-1.0/3.0, 1.0/9.0, 2.0/9.0},
	    new double[]{1.0/9.0, -1.0/3.0, 2.0/9.0},
	    new double[]{1.0/6.0,  0.0,    -1.0/6.0}};
	    
	/* Anfangsverteilung a */           
        double[] a = new double[]{1.0/3.0,1.0/3.0,1.0/3.0};
	    
	/* Initialisierung einer Zufallsgroesse 
	   (zufaelliger Startwert = aktuelle Zeit) */
	Random random = new Random(System.currentTimeMillis());

	Markov sim = new Markov(intMartix_Q, a);
	/* Simulationsschleife */
	for(int n = 0; n < sample_size; n++){
	/**
	 * Bestimmung des Startzustandes x_0
	 * Entscheidungskriterium:
	 * x = sum_{k ¤ E} k * 1(sum_{i=1...k-1} a_i < z <= sum_{i=1...k}) a_i  
	 **/
	    z = random.nextDouble(); /* erzeuge U(0,1)-Pseodozufallszahl Z_0 */
	    for(int k=1; k <= 3; k++){
		if(sim.a_add[k-1] < z && z <= sim.a_add[k]){ 
		    x = k;  break;  }
	    }

	    /*  Verweildauer T_0 ~ Exp(q[X_0]), wobei hier 
		x-1 eingesetzt wird, da 1. Eintrag im Array den Index 0 hat */
	    t = -Math.log(random.nextDouble())/sim.q[x-1];

	    while(t < max_time){
                /* erzeuge U(0,1)-Pseodozufallszahl Z_n */
		z = random.nextDouble();  
		/*  X_{n+1} = phi(X_{n-1},Z_n)   */
		x = sim.transform(x,z);
  	        /*  Verweildauer T_n ~ Exp(q[X_n])  */
		t = t-Math.log(random.nextDouble())/sim.q[x-1];
	    }

	    /* relative Häufigkeiten der Zustände zur Zeit max_time */
	    freq[x-1] += add;
	}

	/* Ausgabe des Ergebnisses */
	 System.out.println("Simulationsergebnis fuer n = "+sample_size+" und t = "+max_time);
	for(int i = 0; i < 3; i++){
	 System.out.println("relative Häufigkeit von Zustand "+(i+1)+" "+out(freq[i]));
	}
    }
}

/*Ergebnisse für 1000 Iterationen und t=10
 *relative Häufigkeit von Zustand 1: 0.321
 *relative Häufigkeit von Zustand 2: 0.104
 *relative Häufigkeit von Zustand 3: 0.575
 *
 *Ergebnisse für 100 Iterationen und t=100
 *relative Häufigkeit von Zustand 1: 0.334
 *relative Häufigkeit von Zustand 2: 0.113
 *relative Häufigkeit von Zustand 3: 0.553
 *
 *Ergebnisse für 100 Iterationen und t=1000
 *relative Häufigkeit von Zustand 1: 0.339
 *relative Häufigkeit von Zustand 2: 0.107
 *relative Häufigkeit von Zustand 3: 0.554
 */