/*=================================================================
 *
 * gateway.c	gateway .MEX file called from coordinate_erg.m 
 *	        erg stands for execute reaction group
 * The calling syntax is:
 *
 *	[cmpdnums_m_out, rxnai_m_out, rxntime_m_out, times_executed_m, eventlog_time_m, eventlog_rxn_m] =
 * 	gateway_erg(rxn_m, rxn_depend_m, cmpdnums_m_in, globtime_m, timestep_m, volume_m, stochalgor_m, update)
 *
 * Primary function is to transfor matlab style inputs into c style inputs into erg.c.  
 * Then to transform returned variables
 * from erg.c into matlab style outputs.
 *
 *   Modifiyed by Slava Ch. to work with stand-alone BNS C-functions
 *                                                         12/06    
 *=================================================================*/

#include "mex.h"
#include "mat.h"
#include "runsim.h"

/* Input Arguments */

#define	rxn_m			prhs[0]
#define	rxn_depend_m		prhs[1]
#define	cmpdnums_m_in	 	prhs[2]
#define	run_constants_m		prhs[3]
#define  monitor_comps_m	prhs[4]

/* Output Arguments */

#define	cmpdnums_m_out		plhs[0]
#define	rxnai_m_out		plhs[1]
#define	rxntime_m_out		plhs[2]
#define	times_executed_m	plhs[3]
#define	eventlog_time_m		plhs[4]
#define	eventlog_rxn_m		plhs[5]
#define	ssdata_m		plhs[6]
#define	agevent_m		plhs[7]


#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;
}


/*
 *  Create rnsx structure in the same way as stand-alone
 *            Slava Ch.   12/06
 */
void get_reaction_mat(int i, mxArray *rxn_m_i, struct rxn *rxn_ptr,  struct array_sizes arsz)
{
    mxArray *mxa_tmp;
    double *temp_array;
    
    int index;
    int strlen, status;
    
        mxa_tmp = mxGetField(rxn_m_i,0,"type");
	    strlen = (mxGetN(mxa_tmp)) + 1;

	    status = mxGetString(mxa_tmp, rxn_ptr[i].type, strlen);
        rxn_ptr[i].rxn_index = assign_rxn_index(rxn_ptr[i].type);

  		if (status != 0) 
    		mexWarnMsgTxt("Not enough space. String is truncated.");
/*  Get substrates   */ 		
	    mxa_tmp = mxGetField(rxn_m_i,0,"substrates");
	    arsz.num_subs[i] = mxGetM(mxa_tmp);
        
        rxn_ptr[i].substrates = (int *)calloc(arsz.num_subs[i], sizeof(int));
        rxn_ptr[i].numsubs_used = (int *)calloc(arsz.num_subs[i], sizeof(int));
                  
		temp_array = mxGetData(mxa_tmp);   
        switch (rxn_ptr[i].rxn_index) {
			case BINDING:
			case UNIQUE_BI_BI:
			case UNBINDING:
			case UNI_UNI:
			case UNI_MULTI:
			case HOMO_LTOR:
			{
    		    for (index = 0; index < arsz.num_subs[i]; index++) 
			 		   rxn_ptr[i].substrates[index] = (int) temp_array[index];
 	   	  	
				break;
			}

			case TRANSV1:
			case GENM:
            case UNDEFINED:
			{
    			for (index = 0; index < arsz.num_subs[i]; index++) 
                {
			  rxn_ptr[i].substrates[index]  = (int) temp_array[index];
                    rxn_ptr[i].numsubs_used[index]= (int) temp_array[index + arsz.num_subs[i]];
 	   	  		}
				break;
			}

			case NTH_ORDER:
			{
              rxn_ptr[i].misc_subval = (double *)calloc(arsz.num_subs[i], sizeof(double));
    				for (index = 0; index < arsz.num_subs[i]; index++) {
					rxn_ptr[i].substrates[index] = (int) temp_array[index];
                    rxn_ptr[i].numsubs_used[index] = (int) temp_array[index + arsz.num_subs[i]];
                    rxn_ptr[i].misc_subval[index] = (int) temp_array[index + 2*arsz.num_subs[i]];
 	   	  	}
				break;
			}
                                         default:
                 mexPrintf("\tA specified reaction type is not defined \n");
        }
/*  Get products   */                    
        mxa_tmp = mxGetField(rxn_m_i,0,"products");
		if (mxa_tmp != NULL) 
        {
			arsz.num_prods[i] = mxGetM(mxa_tmp);
			temp_array = mxGetData(mxa_tmp);
            
            rxn_ptr[i].products = (int *)calloc(arsz.num_prods[i], sizeof(int));
            rxn_ptr[i].numprods_made = (int *)calloc(arsz.num_prods[i], sizeof(int));
			switch (rxn_ptr[i].rxn_index) 
            {
				case BINDING:
				case UNIQUE_BI_BI:
				case UNBINDING:
				case UNI_UNI:
				case HOMO_LTOR:
				{
  	  				for (index = 0; index < arsz.num_prods[i]; index++) 
						rxn_ptr[i].products[index] = (int)temp_array[index];
	
                    mxa_tmp = mxGetField(rxn_m_i,0,"StoP");
                    if (mxa_tmp != NULL) 
                    {
                        temp_array = mxGetData(mxa_tmp);
                        rxn_ptr[i].num_rates = mxGetM(mxa_tmp);
                        rxn_ptr[i].rates = mxCalloc(rxn_ptr[i].num_rates , sizeof(struct rate_data));
                        for (index=0; index<rxn_ptr[i].num_rates; index++)
                                    rxn_ptr[i].rates[index].value = temp_array[index];
                    }
                    
                   break;
				}

				case UNI_MULTI:
				case TRANSV1:
                case GENM: 
                case NTH_ORDER:    
				{
                    for (index = 0; index < arsz.num_prods[i]; index++) 
                    {
                        rxn_ptr[i].products[index] = (int)temp_array[index];
                        rxn_ptr[i].numprods_made[index] = (int)temp_array[index + arsz.num_prods[i]];
                    }                  
                    mxa_tmp = mxGetField(rxn_m_i,0,"StoP");
                    if (mxa_tmp != NULL) 
                    {
                        temp_array = mxGetData(mxa_tmp);
                        rxn_ptr[i].num_rates = mxGetN(mxa_tmp);                  
                        rxn_ptr[i].rates = mxCalloc(rxn_ptr[i].num_rates , sizeof(struct rate_data));
                        for (index=0; index<rxn_ptr[i].num_rates; index++)
                                rxn_ptr[i].rates[index].value = temp_array[index];                        
             
                     }
                    
                        
					break;
				}

				case UNDEFINED:
				{
                    for (index = 0; index < arsz.num_prods[i]; index++) 
                    {
                        rxn_ptr[i].products[index] = (int)temp_array[index];
                        rxn_ptr[i].numprods_made[index] = (int)temp_array[index + arsz.num_prods[i]];
                    }
                          
                    mxa_tmp = mxGetField(rxn_m_i,0,"formula");
                    strlen = (mxGetN(mxa_tmp)) + 1;

                    status = mxGetString(mxa_tmp, rxn_ptr[i].kineticlaw, strlen);                 
                    
                    mxa_tmp = mxGetField(rxn_m_i,0,"StoP");
                    if (mxa_tmp != NULL) 
                    {
                        temp_array = mxGetData(mxa_tmp);
                        rxn_ptr[i].num_rates = mxGetN(mxa_tmp);
                        rxn_ptr[i].rates = mxCalloc(rxn_ptr[i].num_rates , sizeof(struct rate_data));
                        for (index=0; index<rxn_ptr[i].num_rates; index++)
                                    rxn_ptr[i].rates[index].value = temp_array[index];
                      }
                 break;
				}
                 default:
                 mexPrintf("a specified reaction type is not defined \n");
			}
		}      
		else 
        {
			arsz.num_prods[i] = 0;
		}
		
        rxn_ptr[i].id = i;
        
/* assign initial ai and time values to output vectors, since these will be modified during execution */
        mxa_tmp = mxGetField(rxn_m_i,0,"ai");
		if (mxa_tmp != NULL)
        	rxn_ptr[i].ai = mxGetScalar(mxa_tmp);


        mxa_tmp = mxGetField(rxn_m_i,0,"time");
		if (mxa_tmp != NULL) 
	        rxn_ptr[i].rxn_time = mxGetScalar(mxa_tmp);
           
      return;
}

static double  ai_sum=0;
static double tnext=0;
static int rxnnum=0;

void mexFunction( int nlhs, mxArray *plhs[], 
		  int nrhs, const mxArray*prhs[] )
     
{ 
    double *cmpd_nums_out, *rxn_ai_out, *rxn_time_out; 
    int *itimes_executed;
    int *itimes_old;
    double *times_executed;
    double *eventlog_time, eventlog_maxelements ;
    double *eventlog_rxn;
    double *ssdata;
    double *agevent;
    int ssdata_num;
    struct rxn *rxns;
    int *cmpd_init;
    double *monitor_comps_tmp;
    struct array_sizes arsz;
    struct run_constants r_c;
    struct cmpd_data *cmpds;
    double *evlog[2];
    
    int ** stoich_mat = NULL;
    double rand_numb, d_rxnnum;
      

	int buflen, status;
    int event_num;
    int i, j; 
    unsigned int ran_seed;
    double sim_time;
    double update;
    mxArray *mxa_tmp, *rxn_m_i;
    double *temp_array;
    
    /* Check for proper number of arguments */

    if (nrhs != 5) { 
	mexErrMsgTxt("Five input arguments required."); 
    } else if (nlhs != 8) {
	mexErrMsgTxt("Eight output arguments required."); 
    } 
   
    arsz.num_rxns = mxGetNumberOfElements(rxn_m);
    arsz.num_cmpds = mxGetM(cmpdnums_m_in); 
    arsz.num_monitor = mxGetN(monitor_comps_m);  	
/********************************************************
* process the input data
*********************************************************/
    
    /* The .rxn structure.  Assign its values to rxns c structure */ 
    
    rxns = (struct rxn *) mxCalloc(arsz.num_rxns,sizeof(struct rxn));
    arsz.num_subs = (int *) mxCalloc(arsz.num_rxns,sizeof(int));
    arsz.num_prods = (int *) mxCalloc(arsz.num_rxns,sizeof(int));
    itimes_executed = (int *) mxCalloc(arsz.num_rxns,sizeof(int));
    itimes_old = (int *) mxCalloc(arsz.num_rxns,sizeof(int));
 
    for (i = 0; i < arsz.num_rxns; i++) {

   	    rxn_m_i = mxGetCell(rxn_m,i);	    
      get_reaction_mat(i, rxn_m_i, rxns, arsz);
           itimes_executed[i] = 0;      
    }
 
/* The .depend input.  Assign its values to rxn_depend c array of pointers to double arrays */ 
	arsz.num_dependent = (int *) mxCalloc(arsz.num_rxns+1,sizeof(int));

     calc_dependency(&arsz, rxns);
     
    /* The .cmpd_nums input.  Assign its values to cmpd_nums_in c structure */ 
       cmpds = mxCalloc( arsz.num_cmpds, sizeof(struct cmpd_data));
       cmpd_init = (int *) mxCalloc(arsz.num_cmpds,sizeof(int));
       temp_array = mxGetData(cmpdnums_m_in);
       for (j=0; j<arsz.num_cmpds; j++)
       {
                cmpds[j].value = (int)temp_array[j];
                cmpd_init[j] = temp_array[j];
       }
    /* the monitor_compounds id values */
    
    if (arsz.num_monitor != 0)
    {
    	monitor_comps_tmp = mxGetPr(monitor_comps_m);
    	r_c.comps_to_mon = (int *) mxCalloc(arsz.num_monitor,sizeof(int));
        for(i = 0; i < arsz.num_monitor; i++) 
		    r_c.comps_to_mon[i] = (int)monitor_comps_tmp[i]; 
    }
    else 
    {
       	r_c.comps_to_mon = (int *) mxCalloc(1,sizeof(int));
        r_c.comps_to_mon[0] = -1;
   	}

    /* global time, timestep, volume, stoch_algor, update */
    mxa_tmp = mxGetField(run_constants_m,0,"start_time");
    r_c.globtime = *(double *) mxGetData(mxa_tmp);
    mxa_tmp = mxGetField(run_constants_m,0,"timestep");
    r_c.timestep = *(double *) mxGetData(mxa_tmp);
    mxa_tmp = mxGetField(run_constants_m,0,"volume"); 
    r_c.volume = *(double *) mxGetData(mxa_tmp);
    mxa_tmp = mxGetField(run_constants_m,0,"stoch_algor");
    r_c.stoch_algor = *(double *) mxGetData(mxa_tmp);
    mxa_tmp = mxGetField(run_constants_m,0,"randseed");
    rand_numb = *(double *) mxGetData(mxa_tmp);
    r_c.randseed = (int ) rand_numb;
    mxa_tmp = mxGetField(run_constants_m,0,"snapshot_interval");
    r_c.snapshot_interval = *(double *) mxGetData(mxa_tmp);
    mxa_tmp = mxGetField(run_constants_m,0,"evlog_flag");
    r_c.evlog_flag = *(double *) mxGetData(mxa_tmp);
    mxa_tmp = mxGetField(run_constants_m,0,"eventlog_arraysize");     
    eventlog_maxelements = *(double *)mxGetData(mxa_tmp);
    r_c.eventlog_maxelements = (int) eventlog_maxelements;
    mxa_tmp = mxGetField(run_constants_m,0,"update");
    update = *(double *) mxGetData(mxa_tmp);
    if (r_c.stoch_algor == 3)
    {
        mxa_tmp = mxGetField(run_constants_m,0,"tau_error");
        r_c.tau_error = *(double *) mxGetData(mxa_tmp);	
    }

/*    mxa_tmp = mxGetField(run_constants_m,0,"tnext");
    tnext = *(double *) mxGetData(mxa_tmp);
    mxa_tmp = mxGetField(run_constants_m,0,"rxnnum");
    d_rxnnum = *(double *) mxGetData(mxa_tmp);
    rxnnum = (int ) d_rxnnum;
    mxa_tmp = mxGetField(run_constants_m,0,"ai_sum");
    ai_sum = *(double *) mxGetData(mxa_tmp);
 */



    
 /********************************************************
* setup the output variables, except for event log stuff
*********************************************************/
            
    /* Create matrices for the return arguments */ 
    cmpdnums_m_out = mxCreateDoubleMatrix(arsz.num_cmpds, 1, mxREAL); 
    rxnai_m_out = mxCreateDoubleMatrix(arsz.num_rxns,1,mxREAL);
    rxntime_m_out = mxCreateDoubleMatrix(arsz.num_rxns,1,mxREAL);
    times_executed_m = mxCreateDoubleMatrix(arsz.num_rxns,1,mxREAL);
    eventlog_time_m = mxCreateDoubleMatrix(eventlog_maxelements,1,mxREAL);	
    eventlog_rxn_m = mxCreateDoubleMatrix(eventlog_maxelements,1,mxREAL);	

	
    arsz.max_ssdata_rows = (int)floor(r_c.timestep/r_c.snapshot_interval);
   
    if (arsz.max_ssdata_rows == 0) 
	    arsz.max_ssdata_rows = 1;
	ssdata_m = mxCreateDoubleMatrix(arsz.num_monitor+1, arsz.max_ssdata_rows,  mxREAL);
    agevent_m = mxCreateDoubleMatrix(arsz.num_rxns+1, arsz.max_ssdata_rows,  mxREAL);
    
    /* Assign pointers for the output parameters */ 
    cmpd_nums_out = mxGetPr(cmpdnums_m_out);
    rxn_ai_out = mxGetPr(rxnai_m_out);
    rxn_time_out = mxGetPr(rxntime_m_out);
    times_executed = mxGetPr(times_executed_m);
    evlog[0] = mxGetPr(eventlog_time_m);
    evlog[1] = mxGetPr(eventlog_rxn_m);
    ssdata = mxGetPr(ssdata_m);
    agevent =  mxGetPr(agevent_m);
   
    
/********************************************************************        
    /* Do the actual computations in the same way as stand-alone
                                Slava Ch 12/06                 */


  if (update != 0) 
    {
		if (update == 2) 
            {
            setall(r_c.randseed, r_c.randseed);
            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;  
           
		}
           else if (update == 1)
           {
       		for (i = 0; i < arsz.num_rxns; i++) 
                  {
			   itimes_executed[i] = 0;
                     itimes_old[i] = 0;
                  }
                 update_rxn_prob_all(rxns, cmpds, &r_c, &arsz, 0);
          }
         else
          {
                  for (i = 0; i < arsz.num_rxns; i++) 
                  {
			 itimes_executed[i] = 0;
                   itimes_old[i] = 0;
                   }
            
            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;  

        }

    }
  
    else 
    {
         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 stoichiometry matrix         */
         create_stoich_matrix(&arsz, rxns,  stoich_mat);

        sim_time = r_c.globtime;
        ssdata_num = 0;
        event_num =0;
        
        r_c.tnext = tnext;
        r_c.rxnnum = rxnnum;
        r_c.ai_sum = ai_sum;  

        while (sim_time < (r_c.globtime + r_c.timestep)) 
        {     
                 if (r_c.stoch_algor == 3) 
                    erg_tau(rxns, &arsz, cmpds, stoich_mat, &r_c, evlog, &event_num, itimes_executed, &sim_time);
                else
                    erg(rxns, &arsz, cmpds, stoich_mat, &r_c, evlog, &event_num, itimes_executed, &sim_time);

            ssdata[ssdata_num*(arsz.num_monitor+1)] = sim_time;
            agevent[ssdata_num*(arsz.num_rxns+1)] = sim_time;
                for (i = 0; i < arsz.num_monitor; i++) 
                {
                j =  r_c.comps_to_mon[i];
                ssdata[ssdata_num*(arsz.num_monitor+1) + i + 1] = - cmpd_init[j] + cmpds[j].value; 
                }
               
                    for (i = 0; i < arsz.num_rxns; i++) 
                    {
                        agevent[ssdata_num*(arsz.num_rxns+1) + i + 1] = 
                                (double)(itimes_executed[i] - itimes_old[i]);
                        itimes_old[i] = itimes_executed[i];
                    }  
            ssdata_num = ssdata_num + 1;
        }
    }  

 /* Create output arrays          */
        for (j=0; j<arsz.num_cmpds; j++)
                cmpd_nums_out[j] = - cmpd_init[j] + cmpds[j].value; 

        for (j=0; j<arsz.num_rxns; j++)
        {
           rxn_ai_out[j] = rxns[j].ai;
           rxn_time_out[j] = rxns[j].rxn_time;
           times_executed[j] = (double) itimes_executed[j];
       }

tnext = r_c.tnext;
ai_sum = r_c.ai_sum;
rxnnum = r_c.rxnnum;

	return;
    
}
