#include "runsim.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;	
	float 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], my_snap[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;
      double ai_sum;
      double tnext;
      int rxnnum;

	FILE *snap;
	FILE *rundata;
	FILE *f_evlog;


	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");
  
	name = scan_line( infile );
    strcpy(snap_file,output_dir);
    strcat(snap_file, name);

    name = NULL;
    
	name = scan_line( infile );
    strcpy(rundata_file,output_dir);
    strcat(rundata_file, name);    
      name = NULL;
      name = scan_line( infile );
      strcpy(evlog_file,output_dir);
      strcat(evlog_file, name);    

     free(name);
     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);

//        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(rundata_file, "wt");
    if (rundata == NULL)
    { 
        printf("Can't open %s file. Exit.\n",rundata_file);
        return(1);
    }
    if (r_c.evlog_flag == 1)
        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 verification
  */                                                                        
/*		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];

        t_points = (int) floor(r_c.time_end / r_c.snapshot_interval) +1; 
        run_id = 0;
//loop for param sets
	for (iparam = 0; iparam < r_c.num_sets; iparam++ ) 
        {

		num_repeats =  r_c.num_repeats[iparam];

                 if (iparam > 0)
                {
//    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;

                 }
 
		//loop for repeats
		for (irepeat = 0; irepeat < num_repeats; irepeat++) 
        {
            run_id ++;
			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);       

           find_next_reaction_erg(&tnext, &rxnnum, &ai_sum, rxns, r_c.stoch_algor, 1, &arsz, 0.0);
            r_c.tnext = tnext;
            r_c.rxnnum = rxnnum;
            r_c.ai_sum = ai_sum;
    
		fprintf(rundata, "grp 1 rxns %d\n", arsz.num_rxns);
		fprintf(rundata, "set %d run %d\n", iparam+1, run_id);
        
        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",snap_file);
            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)   /*  Save every reaction event   */
            { 
                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);
                    }
                }
			}     /* end individual run    */
            
				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]);
			}
/*			fprintf(rundata, "total number of events = %d\n\n", event_tally);   */
            if (r_c.evlog_flag == 1)
            {    
                fwrite(evlog[0], sizeof(double), event_num, f_evlog);
                fwrite(evlog[1], sizeof(double), event_num, f_evlog);
                fclose(f_evlog); 
            }
		}
	}


	fclose(rundata);

        free(r_c.num_repeats);
        free(times_executed);
        free(arsz.num_subs);
        free(arsz.num_prods);
        free(arsz.num_dependent);
        free(output_dir);
     
        free (sets);
        free (StoP);	
 
      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;
}

