#include "runsim.h"
#include <mpi.h>

#define DRAND_RANGE (1.0/(((double)RAND_MAX + 1.0)*((double)RAND_MAX + 1.0)*((double)RAND_MAX + 1.0)))
double drand (void) {
    double d;
    do {
       d = ((double)rand () * (((double)RAND_MAX + 1.0) * ((double)RAND_MAX + 1.0)) 
       		+ (double)rand () * ((double)RAND_MAX + 1.0) + (double)rand ()) * DRAND_RANGE;
    } while (d >= 1 || d == 0); /* avoid Round off, exclude 0 */
    return d;
}


static char*
get_line( FILE *fp );
static char*
scan_line ( FILE *fp);

int  main(int argc, char *argv[])
{
    	int i, iparam, irepeat, num_repeats, event_tally;
        int j, num_rates, run_id;	
	double rand_numb;
   	double sim_time;
	int *cmpd_nums;
	int *cmpd_init;
      char *output_dir; 
      char *name;
      char *pch;
	char snap_file[264], evlog_file[264];
	char constants_file[264], rundata_file[264];
    char  my_evlog[264];
	FILE *infile;
	struct array_sizes arsz;
	struct run_constants r_c;
	time_t start_time;
	struct rxn *rxns;
        struct cmpd_data *cmpds;
	struct set_data *sets;
        double *evlog[2];
	int *times_executed;
    int *times_exec_old;
        int ** stoich_mat = NULL;
        double **StoP;
    unsigned int ran_seed;
    int event_num, num_logs, t_points;
    int log_maxelements = 1000000;
    
	FILE *snap;
	FILE *rundata;
	FILE *f_evlog;

/*    MPI part       Slava Ch. 01/07          */
        
        int my_id, numprocs;
        char my_snap[256], my_rundata[256];
        int *param_sets[2];
        int tot_runs, mat_start;
        int node_runs, remain;
        int iparam_old;


           MPI_Init(&argc,&argv);
           MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
           MPI_Comm_rank(MPI_COMM_WORLD,&my_id);

	time(&start_time);
     
    r_c.eventlog_maxelements = log_maxelements;
    
      infile = fopen("constants.txt", "rt");
      output_dir = get_line(infile);

  /* Add slash to the end of the string 
   It depends on the system (Windows vs. UNIX)       */
      pch = strchr(output_dir, '/');
      if (pch != NULL)
       strcat(output_dir, "/");
      else
       strcat(output_dir, "\\");          
      name = scan_line( infile );
      fclose(infile);
/* Create constants_file name     */
      strcpy(constants_file,output_dir);
      strcat(constants_file, name);
     name = NULL;
/* Read constants_file and create output file_names  */ 
	infile = fopen(constants_file, "rt");

        fscanf(infile, "%s\n", snap_file);
        fscanf(infile, "%s\n", rundata_file);
        fscanf(infile, "%s\n\n", evlog_file);
  
     free(name);

        sprintf(my_rundata, "%s.%d", rundata_file,my_id);


     fscanf(infile,"\n");
  
	fscanf(infile, "num_sets = %d\n", &r_c.num_sets);
	r_c.num_repeats = (int *)calloc(r_c.num_sets,sizeof(int));
    
	fscanf(infile, "num_rxns = %d\n", &arsz.num_rxns);
	rxns = calloc(arsz.num_rxns, sizeof(struct rxn));
        times_executed = (int *)calloc(arsz.num_rxns,sizeof(int));
        times_exec_old = (int *)calloc(arsz.num_rxns,sizeof(int));
        StoP = (double **) calloc(arsz.num_rxns,sizeof(double *));

	arsz.num_subs = (int *)calloc(arsz.num_rxns,sizeof(int));
	arsz.num_prods = (int *)calloc(arsz.num_rxns,sizeof(int));
	arsz.num_dependent = (int *)calloc(arsz.num_rxns+1,sizeof(int));

	sets  = calloc(r_c.num_sets, sizeof(struct set_data));

	fscanf(infile,"num_cmpds = %d\n\n", &arsz.num_cmpds);
    
	cmpd_nums = (int *)calloc(arsz.num_cmpds, sizeof(int));
	cmpd_init = (int *)calloc(arsz.num_cmpds, sizeof(int));
	fclose(infile);
        cmpds = calloc(arsz.num_cmpds, sizeof(struct cmpd_data)); 

	read_constants(rxns, &arsz, cmpds, &r_c, sets, constants_file);

        r_c.randseed += my_id;

/*        Save initial number of molecules              */
		for (i = 0; i < arsz.num_cmpds; i++)
			cmpd_init[i] = cmpds[i].value;
/*        Save initial kinetic constants         */
		for (i = 0; i < arsz.num_rxns; i++)
                {       
                        num_rates = rxns[i].num_rates;
                        StoP[i] = (double *) calloc(num_rates,sizeof(double ));
                        for (j=0; j < num_rates; j++)    
			StoP[i][j] = rxns[i].rates[j].value;
                }
    
        calc_dependency(&arsz, rxns);

	rundata = fopen(my_rundata, "wt");
    if (rundata == NULL)
    { 
        printf("Can't open %s file. Exit.\n",my_rundata);
        return(1);
    }

           tot_runs = 0;
/*     Calculate ithe total number of runs and number of runs per processor   */        
	  for (iparam = 0; iparam < r_c.num_sets; iparam++ ) 
		                tot_runs +=  r_c.num_repeats[iparam];
               if ( numprocs > tot_runs) 
                                  numprocs = tot_runs;
               
               node_runs = floor(tot_runs / numprocs);
               remain = tot_runs %  numprocs;
         if (my_id < remain)
                 node_runs += 1;
 
          mat_start = my_id*node_runs;
         if (my_id >= remain)
            mat_start = my_id*node_runs +  remain;
                        
/*    Allocate memory for param_sets matrix.
      It will have parameters sets/run_number info for each of the nodes   */ 
          for (i=0; i<2; i++)
                param_sets[i] = (int *)calloc(tot_runs, sizeof(int)); 

/*     Fill up param_sets  matrix       */
         
             j= 0; 
            for (iparam = 0; iparam < r_c.num_sets; iparam++ )
             for (i=0; i <  r_c.num_repeats[iparam]; i++ )
             {
                  param_sets[0][j] = iparam;
                  param_sets[1][j] = i+1;   /* Run number for this parameters set  */
                  j ++;           
             }


        for (i=0; i<2; i++)
            evlog[i] = (double *)calloc(log_maxelements, sizeof(double));
        
        stoich_mat = (int **) calloc(arsz.num_rxns + 1, sizeof(int *));
         for (i=0; i< arsz.num_rxns +1; i++)
            stoich_mat[i] = (int *) calloc(arsz.num_cmpds, sizeof(int ));
     
         create_stoich_matrix(&arsz, rxns,  stoich_mat);

/*	rand31pmc_seedi(r_c.randseed);
 Use original rand() for validation against bns2p5p4.
  */    
/*		ran_seed = (unsigned int)floor(r_c.randseed * (double)UINT_MAX);
        srand(ran_seed);	*/
        setall(r_c.randseed, r_c.randseed);
        
/*    Initialize working number of molecules         */
		for (i = 0; i < arsz.num_cmpds; i++)
			cmpd_nums[i] = cmpd_init[i];

/*     Run simulations                */
       
        iparam_old = 0; 
		t_points = (int) floor(r_c.time_end / r_c.snapshot_interval) +1; 

	for (irepeat = 0; irepeat < node_runs; irepeat++ ) 
        {
                iparam = param_sets[0][mat_start + irepeat]; 
                run_id = mat_start + irepeat+1;

               fprintf(rundata, "grp 1 rxns %d\n", arsz.num_rxns);
               fprintf(rundata, "set %d run %d\n", iparam+1, run_id);

                 if (iparam > 0 || (iparam != iparam_old) )
                {
/*    Initialize number of molecules       */
		for (i = 0; i < arsz.num_cmpds; i++)
			cmpd_nums[i] = cmpd_init[i];

/*    Initialize kinetic constants   */
		for (i = 0; i < arsz.num_rxns; i++)
                {       
                 num_rates = rxns[i].num_rates;
                        for (j=0; j < num_rates; j++)    
				rxns[i].rates[j].value = StoP[i][j];
                }

/*       Modify number of molecules           */
                  for (i=0; i< sets[iparam].num_cmpds; i++)
                      cmpd_nums[sets[iparam].cmpds[i].id] = sets[iparam].cmpds[i].value;
 
/*        Modify kinetic constants           */
                  for (i=0; i< sets[iparam].num_rates; i++)
               rxns[sets[iparam].rxn_id[i]].rates[sets[iparam].rates[i].id].value = sets[iparam].rates[i].value;
               iparam_old = iparam;

                 }
 

			sim_time = r_c.time_start;
/*       Assign number of molecules    */
 			for (i = 0; i < arsz.num_cmpds; i++)
                                 cmpds[i].value = cmpd_nums[i];
			for (i = 0; i < arsz.num_rxns; i++)
            {
			            	times_executed[i] = 0;
                            times_exec_old[i] = times_executed[i];
            }
            update_rxn_prob_all(rxns, cmpds, &r_c, &arsz, 0);
            

         sprintf(my_snap, "%s_%d", snap_file,run_id);
	     snap = fopen(my_snap, "wt");
    
    if (snap == NULL)
    { 
        printf("Can't open %s- file. Exit.\n",my_snap);
        return(1);
    }
        	fprintf(snap, "run %d points %d n_mon %d n_rxn %d\n", irepeat+1, t_points, arsz.num_monitor, arsz.num_rxns);
            if (r_c.evlog_flag == 1)
            { 
                num_logs = 0;
                sprintf(my_evlog,"%s_%d_%d",evlog_file,run_id,num_logs+1);
                f_evlog = fopen(my_evlog, "wb");
                if (f_evlog == NULL)
                { 
                    printf("Can't open %s file. Exit.\n",evlog_file);
                    return(1);
                }
            }   
            
            event_num = 0; 
   	  while (sim_time < r_c.time_end) 
            {  
                if (fmod(sim_time, r_c.snapshot_interval) == 0)
                {
                    fprintf(snap, "%g", sim_time);    
                    for (i = 0; i < arsz.num_monitor; i++) 
                        fprintf(snap,"\t%d",cmpds[r_c.comps_to_mon[i]].value);
				    
                    for (i = 0; i < arsz.num_rxns; i++)
                       { 
                           fprintf(snap,"\t%d",times_executed[i]- times_exec_old[i]);
                           times_exec_old[i] = times_executed[i];
                       }  

                    fprintf(snap,"\n");
                }
                if (r_c.stoch_algor == 3) 
                    erg_tau(rxns, &arsz, cmpds, stoich_mat, &r_c, evlog, &event_num, times_executed, &sim_time);
                else
                    erg(rxns, &arsz, cmpds, stoich_mat, &r_c, evlog, &event_num, times_executed, &sim_time);
                
                if (event_num == log_maxelements)
                {
                   fwrite(evlog[0], sizeof(double), event_num, f_evlog);
                   fwrite(evlog[1], sizeof(double), event_num, f_evlog);
                   fclose(f_evlog);
                   num_logs ++;
                   event_num = 0;
                    sprintf(my_evlog,"%s_%d_%d",evlog_file,run_id,num_logs+1);
                    f_evlog = fopen(my_evlog, "wb");
                    if (f_evlog == NULL)
                    { 
                        printf("Can't open %s file. Exit.\n",evlog_file);
                        return(1);
                    }
                }
            }
            
				fprintf(snap, "%g", sim_time);    

				for (i = 0; i < arsz.num_monitor; i++) 
					fprintf(snap,"\t%d",cmpds[r_c.comps_to_mon[i]].value);
			    for (i = 0; i < arsz.num_rxns; i++)
                       { 
                           fprintf(snap,"\t%d",times_executed[i]-times_exec_old[i]);
                           times_exec_old[i] = times_executed[i];
                       } 
                fclose(snap);
				
			event_tally = 0;
            
			for(i = 0; i < arsz.num_rxns; i++) 
                       {
				event_tally += times_executed[i];
				fprintf(rundata, "%d %d\n", i+1, times_executed[i]);
			}
            if (r_c.evlog_flag == 1)
            {    
                fwrite(evlog[0], sizeof(double), event_num, f_evlog);
                fwrite(evlog[1], sizeof(double), event_num, f_evlog);
            }
	        }



        if (r_c.evlog_flag == 1)
	                    fclose(f_evlog); 
	fclose(rundata);

        free(r_c.num_repeats);
        free(times_executed);
        free(times_exec_old);
        free(arsz.num_subs);
        free(arsz.num_prods);
        free(arsz.num_dependent);
        free(output_dir);
     
        free (sets);
        free (StoP);
      MPI_Finalize();	
 
      return(0);
}


/**
 * The function get_line reads a line from a file and
 * returns it as a string.  
 */
static char*
scan_line( FILE *fp )
{
  /* reads lines of arbitrary length from fp */
  char s[264], *line;
  int len,i;
  
  line = NULL;
  fscanf(fp,"%s",s);
    len = strlen(s);
    
      line = calloc( len+1, sizeof(char) );
      for (i=0; i<len; i++)
      {
          line[i] = s[i];
      }
  return line;
}

/**
 * The function get_line reads a line from a file (in this case "stdin" and
 * returns it as a string.  It is taken from the utilities library of the
 * VIENNA RNA PACKAGE ( http://www.tbi.univie.ac.at/~ivo/RNA/ )
 */
static char*
get_line( FILE *fp )
{
  /* reads lines of arbitrary length from fp */
  char s[264], *line, *cp;
  
  line = NULL;
  do
  {
    if ( fgets( s, 512, fp ) == NULL ) break;
    cp = strchr( s, '\n' );
    if ( cp != NULL ) *cp = '\0';
    if ( line == NULL )
      line = calloc( 1+strlen(s), sizeof(char) );
    else
      line = (char *)realloc( line, 1+strlen( s )+strlen( line ) );
    strcat( line, s );
  } while ( cp == NULL );

  return line;
}

