/*******************************
This file gives a simple implementation for the Gillespie Algorithm introduced
in Project 1. It demonstrates how the reactions are simulated, how the trajectory
and ensemble data are collected, and how the stastistics are collected, based on
a simple reaction. The system you will handle in project 1 will be much more
complicated and special data structures will be needed for efficiency.
********************************/

import java.io.*;

public class SimpleReactionExample {

  static final double T = 10; // final time

  // calculate the offset of the time for the next reaction to fire
  public static double nexttau(double propensity) {
      if (propensity == 0) return T+1; // A big time
      double r = Math.random();   // Generate a random number 
      return (- Math.log(r)/propensity);
  }

  public static void main(String[] args) {
    
      // This example shows how to simulate a simple reaction S1 -> 0
      FileOutputStream file_trajectory, file_ensemble; 
      PrintStream p_trajectory, p_ensemble; 	  
	  
      double k = 0.1; // reaction rate for the decaying process
      int x, initial_x = 5000; // initial state for S1
      double t=0, initial_t = 0; // initial time 
      double a, tau, next_t; 
      double[] record_t = new double[initial_x+1]; 
      int[] record_x = new int[initial_x+1]; 
      int record_i, j, total_run = 10000;
      int[] ensemble_x = new int[total_run]; 
      
      record_t[0] = t; 
      record_x[0] = initial_x; 
      record_i = 1; 
	  
      // Do total_run simulation runs
      for (j = 0; j < total_run; j++) { 
	  // Initialize the run.
	  t = initial_t; 
	  // You will need to initialize propensities and popluations.
	  // But that is pretty simple for this example.
	  x = initial_x; 
	  tau = nexttau(k*x);

	  // Now run a simulation loop. We stop when we reach the end time.
	  while (t <= T) {
	      // Determine the reaction with the smallest nexttau increment.
	      // That is the next reaction that will fire.
	      // In this example, there is only one, so we can skip that step.

	      next_t = t + tau; // This is the time of the next event
		  
	      x--; // When this reaction fires, x decreases by 1. 
	      t = next_t; // Update the simulation time
	      if (j == total_run - 1) { // Save the trajectory of only the last run
		  record_t[record_i] = t; 
		  record_x[record_i] = x; 
		  record_i++; // You should check if the record is out of bound.
		              // We skipped this step for this simple code. 
	      }
	      // Now update the propensities for all affected reactions.
	      // Again, this example has only one reaction.
	      a = k*x; // Calculate propensity function for this simple reaction
	      tau = nexttau(a); // Generate the random time for next firing of this reaction
	  }
	  ensemble_x[j] = x; // Add final state to ensemble when run is finished. 
      } 
	
      try  {
	  file_trajectory = new FileOutputStream("SimpleReaction.trajectory");
	  file_ensemble  = new FileOutputStream("SimpleReaction.ensemble");

	  p_trajectory = new PrintStream(file_trajectory);
	  p_ensemble = new PrintStream(file_ensemble); 

	  for (int i = 0; i < record_i; i ++) {
	      p_trajectory.println(record_t[i] + "\t" + record_x[i]); 
	  }
					  
	  p_trajectory.close();
		  
	  for (int i = 0; i < total_run; i++) { 
	      p_ensemble.println(ensemble_x[i]); 
	  } 
	  p_ensemble.close(); 
      } catch (Exception e) {
	  System.err.println ("Error writing to file");
      }
	  
      double sum_x = 0; 			  
      for (j = 0; j < total_run; j++) {
	  sum_x += ensemble_x[j]; 
      }
      double mean_x = sum_x/total_run; 
	  
      double var_x = 0; 
      for (j = 0; j < total_run; j++) {
	  var_x += (ensemble_x[j] - mean_x)*(ensemble_x[j] - mean_x); 
      }
      var_x /= (total_run - 1); 
      double std_x = Math.sqrt(var_x); 
				  

      System.out.println("The mean of the final state of x at time " + T + " is: " +   mean_x);
      System.out.println("The variance of the final state of x at time " + T + " is: " + var_x);
      System.out.println("The standard deviation of the final state of x at time " + T + " is: " + std_x);
	  
  }

}