simulation_clib.c

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<stdbool.h>
#include<string.h>
#include "pcg_basic.h"
#include "simulation_clib.h"

const double ZERO=1.0e-10;

const size_t LINE_SIZE=300; 
const char CONFIG_FILE_DELIM[5]=":";

//SUBSET OF CONFIG PARAMETERS:
struct config_dict {
       int NUM_RANDOM_ICS;
       long double CONVERGENCE_PROXIMITY;
       double SAMEPOINT_PROXIMITY;
       int ITER_FOR_ODE;
       int ITER_FOR_RELAXATION;
       int ITER_FOR_LIMIT_CYCLE;
       double EULER_SIM_TIME;
       double EULER_SIM_TIME_LAST_ANN_STEP;
       double STARTING_NOISE_ANN;
       double STARTING_NOISE_CONSTANT;
       double NOISE_SCALING_FACTOR_ANN;
       double LIMIT_CYCLE_SIM_TIME;
       double EULER_SIM_STEP_SIZE;
       double LIMIT_CYCLE_SIM_STEP_SIZE;
       int MAX_STABLE_STATES;
       int MAX_LIMIT_CYCLES;
       int MAX_ALLOWED_PERIODS;
       int NO_OF_SAMPLED_PERIODS;
       int ALLOWED_ERROR_IN_PERIODS;
}; 

//DEFAULT VALUES:
struct config_dict sim_config={
       .NUM_RANDOM_ICS=1000, 
       .CONVERGENCE_PROXIMITY=1.0e-12,
       .SAMEPOINT_PROXIMITY=0.1,
       .ITER_FOR_ODE=200,
       .ITER_FOR_RELAXATION=10,
       .ITER_FOR_LIMIT_CYCLE=20,
       .EULER_SIM_TIME=10,
       .EULER_SIM_TIME_LAST_ANN_STEP=2000,
       .STARTING_NOISE_ANN=20,
       .STARTING_NOISE_CONSTANT=20,
       .NOISE_SCALING_FACTOR_ANN=0.95,
       .LIMIT_CYCLE_SIM_TIME=10,
       .EULER_SIM_STEP_SIZE=0.01,
       .LIMIT_CYCLE_SIM_STEP_SIZE=0.01,
       .MAX_STABLE_STATES=20,
       .MAX_LIMIT_CYCLES=10,
       .MAX_ALLOWED_PERIODS=100,
       .NO_OF_SAMPLED_PERIODS=3,
       .ALLOWED_ERROR_IN_PERIODS=3
};


/**********************************************************************/
void find_solutions(char *WORK_DIR, \
                    char *FNAME_STATES, \
                    char *FNAME_LIMITCYCLES, \
                    char *FNAME_SUMMARY, \
                    const int MODEL_NO,\
                    int NUM_NODES,int NUM_EDGES,\
                    int NUM_RANDOM_ICS,int ITER_FOR_ODE,\
                    double EULER_SIM_TIME,\
                    double EULER_STEP_SIZE,\
                    long double CONVERGENCE_PROXIMITY,\
                    double TRANS_RATE_FACTOR,\
                    const double *MPR_arrv,\
                    const double *DNR_arrv,\
                    const int *node_type_arrv,\
                    const int *edge_source_arrv,\
                    const int *edge_target_arrv,\
                    const int *edge_type_arrv,\
                    const double *TSH_arrv,\
                    const int *HCO_arrv,\
                    const double *FCH_arrv,\
                    double *EXP_dict_arrv) 
/*------------------------------------------------------------------
It carries out two main operations 
   (1) For each of the initial conditions, it estimates stable state 
       of the network (an n dimensional point) by repeatedly invoking 
       estimate_stable_expression function. 
   (2) After obtaining stable states for all initial conditions, 
       it clusters those stable state solutions by invoking 
       cluster_solutions function.

Definitions of the main data structures: 
   MPR_arrv: maximum production rates for each node
   DNR_arrv: degradation rates for each node  
   node_type_arrv: type for each node 

   edge_source_arrv: source id of each edge. Id of a source is the 
      corresponding index of the arrays MPR_arrv, DNR_arrv, and 
      node_type_arrv
   edge_target_arrv: target id of each edge
   edge_type_arrv: type of each edge 
   TSH_arrv: threshold of each edge
   HCO_arrv: hill coefficient for each edge
   FCH_arrv: fold change for each edge

   EXP_dict_arrv: stable expression level of all the nodes
--------------------------------------------------------------------*/
{
    //CONSTANTS:
    const int TOTAL_EULER_STEPS=(int) (EULER_SIM_TIME/EULER_STEP_SIZE);

    //VARIABLES:
    double *EXP_dict_arr=(double *) EXP_dict_arrv;
    double IC_arr[NUM_NODES]; //saves ICs generated by set_ICs()
    double exp_arr[NUM_NODES]; //saves exps from estimate_stable_expression()

    double minEXP_arr[NUM_NODES]; //minimum expression for all nodes 
    double maxEXP_arr[NUM_NODES]; //maximum expression for all nodes
   
    double EXP_arr[NUM_RANDOM_ICS][NUM_NODES];
    double fX_arr_norm[NUM_RANDOM_ICS];
    double fX_arr[NUM_NODES];
    double zero_arr[NUM_NODES];

    //INDEX VARIABLES:
    int count_pos=0; //keeps track of position in EXP_dict_arr 
    int count_iteration; //counts iteration over IC initialization
    int i=0; //general purpose 

    //read the simulation configations:
    char fname_cfg[30]="racipe.cfg";

    read_config(WORK_DIR,fname_cfg);

    // calculate minimum and maximum allowable expression for each node:
    cal_EXPrange(NUM_NODES,NUM_EDGES,MPR_arrv,DNR_arrv,\
                 edge_source_arrv, edge_target_arrv,\
                 edge_type_arrv, FCH_arrv, minEXP_arr, maxEXP_arr);

    //set exp_arr to zeros: 
    for (i=0;i<NUM_NODES;++i) {
        exp_arr[i]=0.0;
        fX_arr[i]=0.0;
        zero_arr[i]=0.0;
    }
    count_iteration=0;
    while(count_iteration<NUM_RANDOM_ICS)
    { 
        // randomly select initial condition:
        set_ICs(NUM_NODES,minEXP_arr,maxEXP_arr,IC_arr);

        // estimate the stable state for the current initial condition:
        estimate_stable_expression(NUM_NODES,NUM_EDGES,ITER_FOR_ODE,\
                                   TOTAL_EULER_STEPS,EULER_STEP_SIZE,\
				                       CONVERGENCE_PROXIMITY,\
                                   TRANS_RATE_FACTOR,\
                                   MPR_arrv,DNR_arrv,node_type_arrv,
                                   edge_source_arrv,\
                                   edge_target_arrv,edge_type_arrv,\
                                   TSH_arrv,HCO_arrv,FCH_arrv,\
				                       IC_arr,exp_arr,fX_arr);

        // calculate norm of fX:
        fX_arr_norm[count_iteration]=cal_norm(fX_arr, NUM_NODES);

        //save the exp_arr to EXP_dict_arr: 
        //if any value is found <= zero, assign it min value
        for(i=0;i<NUM_NODES;i++){
            EXP_dict_arr[count_pos+i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            EXP_arr[count_iteration][i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
        }
        count_pos+=NUM_NODES;
        count_iteration+=1;
    } //end of while(count_iteration<NUM_RANDOM_ICS) 

    
    // cluster the stable state solutions:
    cluster_solutions(WORK_DIR, \
                      FNAME_STATES, FNAME_LIMITCYCLES, FNAME_SUMMARY,\
                      MODEL_NO, NUM_NODES,NUM_EDGES,\
                      NUM_RANDOM_ICS,ITER_FOR_ODE,\
                      EULER_SIM_TIME,\
                      EULER_STEP_SIZE,\
                      CONVERGENCE_PROXIMITY,\
                      TRANS_RATE_FACTOR,\
                      MPR_arrv,\
                      DNR_arrv,\
                      node_type_arrv,\
                      edge_source_arrv,\
                      edge_target_arrv,\
                      edge_type_arrv,\
                      TSH_arrv,\
                      HCO_arrv,\
                      FCH_arrv,\
                      fX_arr_norm,\
                      EXP_arr); 
    return;
}

/*-----------------------------------------------------------------------*/
void cluster_solutions(char *WORK_DIR, \
                       char *FNAME_STATES, char *FNAME_LIMITCYCLES, \
                       char *FNAME_SUMMARY,\
                       const int MODEL_NO,\
                       int NUM_NODES,int NUM_EDGES,\
                       int NUM_RANDOM_ICS,int ITER_FOR_ODE,\
                       double EULER_SIM_TIME,\
                       double EULER_STEP_SIZE,\
                       long double CONVERGENCE_PROXIMITY,\
                       double TRANS_RATE_FACTOR,\
                       const double *MPR_arrv,\
                       const double *DNR_arrv,\
                       const int *node_type_arrv,\
                       const int *edge_source_arrv,\
                       const int *edge_target_arrv,\
                       const int *edge_type_arrv,\
                       const double *TSH_arrv,\
                       const int *HCO_arrv,\
                       const double *FCH_arrv,\
                       double *fX_arr_norm,\
                       double EXP_arr[][NUM_NODES])
/*------------------------------------------------------------------
This function clusters the solutions obtained from all initial 
conditions into states and limit cycles.

It invokes two functions to carry out its operation:  
  (1) find_states: clusters all the solutions into states
  (2) find_limitcycles: for each solution not used in state calculation,
      this function runs further simulation to obtain the limit cycle
      for that solution.
--------------------------------------------------------------------*/
{
    //LOCAL VARIABLES:
    //variables for limit cycle (LC) calculation:
    bool fX_arr_bool[NUM_RANDOM_ICS];

    //COUNTING VARIABLES:
    int state_count;
    int LC_count;
    int i;

    clock_t time_begin, time_end; 
    double time_spent; 

    //DEFINE FILE NAME:
    //char fname_summary[300]="\0";
    //strcat(fname_summary,WORK_DIR);
    //strcat(fname_summary,"/");
    //strcat(fname_summary,"cellcycle.summary.txt");
    //strcat(fname_summary,"TS.summary.txt");

    state_count=0; 
    LC_count=0; 

    state_count=find_states(WORK_DIR, FNAME_STATES, MODEL_NO,\
                            NUM_NODES, NUM_RANDOM_ICS,\
                            CONVERGENCE_PROXIMITY,\
                            fX_arr_norm,\
                            EXP_arr,\
                            fX_arr_bool);

    time_begin=clock(); 
    LC_count=find_limitcycles(WORK_DIR, FNAME_LIMITCYCLES, MODEL_NO,\
                              NUM_NODES,NUM_EDGES,\
                              NUM_RANDOM_ICS,ITER_FOR_ODE,\
                              EULER_SIM_TIME,\
                              EULER_STEP_SIZE,\
                              CONVERGENCE_PROXIMITY,\
                              TRANS_RATE_FACTOR,\
                              MPR_arrv,\
                              DNR_arrv,\
                              node_type_arrv,\
                              edge_source_arrv,\
                              edge_target_arrv,\
                              edge_type_arrv,\
                              TSH_arrv,\
                              HCO_arrv,\
                              FCH_arrv,\
                              fX_arr_norm,\
                              EXP_arr,\
                              fX_arr_bool);

    time_end=clock();
    time_spent=(double) (time_end-time_begin)/CLOCKS_PER_SEC;
    //write summary:
    FILE *fh_summary;
    //fh_summary=fopen(fname_summary,"a");
    fh_summary=fopen(FNAME_SUMMARY,"a");
    fprintf(fh_summary,"%d\t%d\t%d\t%8.6f\n",MODEL_NO, state_count,\
                                              LC_count,time_spent);
    fflush(fh_summary);
    fclose(fh_summary); 
    return;
}

/*-----------------------------------------------------------------------*/
int find_states(char *WORK_DIR, char *FNAME_STATES, \
                const int MODEL_NO,\
                int NUM_NODES, int NUM_RANDOM_ICS,\
                long double CONVERGENCE_PROXIMITY,\
                double *fX_arr_norm,\
                double EXP_arr[][NUM_NODES],\
                bool fX_arr_bool[NUM_RANDOM_ICS])
{
    //LOCAL VARIABLES:
    int MAX_STABLE_STATES=sim_config.MAX_STABLE_STATES;
    //solution_count_arr stores the number of solutions in each state:
    int solution_count_arr[MAX_STABLE_STATES];
    //solution_arr stores the average expression for all states: 
    double solution_arr[MAX_STABLE_STATES][NUM_NODES];
    //solution_arr_avg stores the average solution for any state:
    double solution_arr_avg[NUM_NODES];

    double testDelta=0.0;
    bool found;

    //COUNTING VARIABLES:
    int count_state;
    int i;

    //DEFINE FILE NAME:
    //char fname_states[300]="\0";
    //strcat(fname_states,WORK_DIR);
    //strcat(fname_states,"/");
    //strcat(fname_states,"cellcycle.states.txt");
    //strcat(fname_states,"TS.states.txt");

    //initialize the solution count for each state: 
    for (i=0;i<MAX_STABLE_STATES;i++) solution_count_arr[i]=0; 

    //SEPARATE ALL THE STABLE STATES INTO CLUSTERS:
    count_state=0; //initialize state count variable
    for (i=0;i<NUM_RANDOM_ICS;i++)
    {
        if (fX_arr_norm[i]>CONVERGENCE_PROXIMITY)
        {
             fX_arr_bool[i]=true; 
             continue;
        }
        fX_arr_bool[i]=false; 
        found = false;
        int sol_no;
        for (sol_no=0; sol_no<count_state;sol_no++)
        {
            testDelta=sum_delta(solution_arr[sol_no],EXP_arr[i],NUM_NODES);
            if (testDelta<=CONVERGENCE_PROXIMITY)
            {
		         int j;
                for(j=0;j<NUM_NODES;j++)
                {
                    solution_arr[sol_no][j]=(solution_arr[sol_no][j]*\
                                             solution_count_arr[sol_no]+\
                                             EXP_arr[i][j])/\
                                             (solution_count_arr[sol_no]+1);
                }
		          solution_count_arr[sol_no]++;
		          found=true;
                break;
             }
         }
         //new solution found:
         if(!found){
        		  int j;
              for (j=0;j<NUM_NODES;j++) {
                  solution_arr[count_state][j]=EXP_arr[i][j]; 
              }
              //printf("fX_arr_norm[i]: %d\t%f\n",count_state,fX_arr_norm[i]);
	           solution_count_arr[count_state]=1;
	           count_state+=1; 
              if(count_state>=MAX_STABLE_STATES) break; 
         }
    }

    //SAVE SOLUTIONS: solution_arr
    FILE *fh_states;
    fh_states=fopen(FNAME_STATES,"a");
    int j;
    int k;

    //write a new line to move the file pointer to the next line 
    //fprintf(fh_states,"\n");  
    for (j=0; j<count_state;j++){
         fprintf(fh_states,"%d\t",MODEL_NO);  //write model index
         fprintf(fh_states,"%d\t",count_state); //write total solution
         fprintf(fh_states,"%d\t",(j+1));     //write solution index
         for(k=0;k<NUM_NODES;k++){            //write gene levels 
             fprintf(fh_states,"%10.6f\t",log(solution_arr[j][k])/log(2));
         }
         fprintf(fh_states,"\n");
    }

    fflush(fh_states);
    fclose(fh_states);
    return count_state;
}

/*-----------------------------------------------------------------------*/
int find_limitcycles(char *WORK_DIR, char *FNAME_LIMITCYCLES, \
                     const int MODEL_NO,\
                     int NUM_NODES, int NUM_EDGES,\
                     int NUM_RANDOM_ICS, int ITER_FOR_ODE,\
                     double EULER_SIM_TIME,\
                     double EULER_STEP_SIZE,\
                     long double CONVERGENCE_PROXIMITY,\
                     double TRANS_RATE_FACTOR,\
                     const double *MPR_arrv,\
                     const double *DNR_arrv,\
                     const int *node_type_arrv,\
                     const int *edge_source_arrv,\
                     const int *edge_target_arrv,\
                     const int *edge_type_arrv,\
                     const double *TSH_arrv,\
                     const int *HCO_arrv,\
                     const double *FCH_arrv,\
                     double *fX_arr_norm,\
                     double EXP_arr[][NUM_NODES],\
                     bool fX_arr_bool[NUM_RANDOM_ICS])
/*------------------------------------------------------------------
This function finds and saves limit cycles. It calls other functions
to carry these out:
   (1) detect_limitcycle: this returns the estimated size of 
      the limit cycle.
   (2) cal_limitcycle: this calculates expression levels along the 
       limit cycle.
    return value: number of limit cycle found for a specific model.
--------------------------------------------------------------------*/
{
    //CONSTANTS:
    const double LIMIT_CYCLE_SIM_TIME=sim_config.LIMIT_CYCLE_SIM_TIME;
    const double LIMIT_CYCLE_SIM_STEP_SIZE=sim_config.LIMIT_CYCLE_SIM_STEP_SIZE;
    const int LIMIT_CYCLE_SIM_STEPS=(int) (LIMIT_CYCLE_SIM_TIME/\
                                           LIMIT_CYCLE_SIM_STEP_SIZE);
    const int MAX_LIMIT_CYCLES=sim_config.MAX_LIMIT_CYCLES;
    const int ALLOWED_ERROR_IN_PERIODS=sim_config.ALLOWED_ERROR_IN_PERIODS;

    //LOCAL VARIABLES:
    int LC_period_arr[MAX_LIMIT_CYCLES];
    double **LC_exp_arr;

    double start_exp_arr[NUM_NODES]; //start point for calculating limit cycle
    double LC_start_exp_arr[NUM_NODES];

    int period=0; //period of the limit cycle

    //COUNTING VARIABLE:
    int count_LC;
    int i;

    //DEFINE FILE HEADER:
    FILE *fh_LCs;
    //fh_LCs=fopen(fname_LCs,"a");
    fh_LCs=fopen(FNAME_LIMITCYCLES,"a");

    //write a new line to move the file pointer to the next line 
    //fprintf(fh_LCs,"\n");  

    //initialize LC_period_arr with zeros:
    for (i=0;i<MAX_LIMIT_CYCLES;i++){
        LC_period_arr[i]=0; 
    }

    //initialize LC_start_exp_arr
    //it will be defined in detect_limitcycle()
    for (i=0;i<NUM_NODES;i++){
       LC_start_exp_arr[i]=0.0;
    }

    //initialize limit cycle count variable:
    count_LC=0;
    for (i=0;i<NUM_RANDOM_ICS;i++)
    {
	      if (fX_arr_bool[i]==false) continue;

         //calculate limit cycle:
	      fX_arr_bool[i]=false;
         period=detect_limitcycle(NUM_NODES,NUM_EDGES,\
                                  LIMIT_CYCLE_SIM_STEPS,\
                                  LIMIT_CYCLE_SIM_STEP_SIZE,\
                                  CONVERGENCE_PROXIMITY,\
                                  TRANS_RATE_FACTOR,\
                                  MPR_arrv,DNR_arrv,\
                                  node_type_arrv,\
                                  edge_source_arrv,\
                                  edge_target_arrv,\
                                  edge_type_arrv,\
                                  TSH_arrv,\
                                  HCO_arrv,\
                                  FCH_arrv,\
                                  EXP_arr[i],\
                                  LC_start_exp_arr);

        if(period>0){
             //check for repeated limit cycle:
             int j;
             int k;
             //NEW_SIM_STEPS: simulation steps for calculating expression 
             //levels along limit cycle:
             int NEW_SIM_STEPS;

             //MAX_DIST: used as an allowed error between the converged expression 
             //levels (based on Initial Conditions) and the expression 
             //level along the limit cycle:
             double MAX_DIST; 
             double tmp_norm=0.0;
             bool is_same;

             is_same=false;
             for(k=0;k<count_LC;k++){
                  if(abs(period-LC_period_arr[k])<=ALLOWED_ERROR_IN_PERIODS){
                     is_same=true;
                     break;
                  }
             }
             //if it is the same limit cycle, skip subsequent calculation:
             if(is_same) continue;

             LC_period_arr[count_LC]=period;
             count_LC++;
             LC_exp_arr=allocate_LCmemory(period+2,NUM_NODES);
             NEW_SIM_STEPS=period+1;

             MAX_DIST=cal_limitcycle(NUM_NODES,NUM_EDGES,\
                                     NEW_SIM_STEPS,\
                                     LIMIT_CYCLE_SIM_STEP_SIZE,\
                                     CONVERGENCE_PROXIMITY,\
                                     TRANS_RATE_FACTOR,\
                                     MPR_arrv,DNR_arrv,\
                                     node_type_arrv,\
                                     edge_source_arrv,\
                                     edge_target_arrv,\
                                     edge_type_arrv,\
                                     TSH_arrv,\
                                     HCO_arrv,\
                                     FCH_arrv,\
                                     LC_start_exp_arr,\
                                     LC_exp_arr);

             //relax MAX_DIST by factoring it by 10:
             MAX_DIST *= 10.0;

             //change fX_arr_bool status to false for all the positions that 
             //all give to this limit cycle:
             for (j=0;j<NUM_RANDOM_ICS;j++){
	             if (fX_arr_bool[j]==false) continue;
                for(k=0;k<period+2;k++){
                     tmp_norm=sum_delta(EXP_arr[j],LC_exp_arr[k],NUM_NODES);
                     if(tmp_norm<=MAX_DIST) {
	                     fX_arr_bool[j]=false;
                        break;
                     }
                }
             }
            
             //write LC expressions to file: 
             write_limitcycle(fh_LCs, MODEL_NO, NUM_NODES,\
                              count_LC, period, NEW_SIM_STEPS,\
                              LC_exp_arr);
             //note: NEW_SIM_STEPS is larger than period by 2

             //free the memory allocated for the current LC: 
             free_LCmemory(LC_exp_arr,period+2,NUM_NODES);
             if(count_LC>=MAX_LIMIT_CYCLES) break;
        }
    }
    fflush(fh_LCs);
    fclose(fh_LCs); 
    return count_LC;
}

/*-----------------------------------------------------------------------*/
int detect_limitcycle(int NUM_NODES,int NUM_EDGES,\
                      int SIM_STEPS,double SIM_STEP_SIZE,\
                      long double CONVERGENCE_PROXIMITY,\
                      double TRANS_RATE_FACTOR,\
                      const double *MPR_arrv,const double *DNR_arrv,\
                      const int *node_type_arrv,\
                      const int *edge_source_arrv,\
                      const int *edge_target_arrv,\
                      const int *edge_type_arrv,\
                      const double *TSH_arrv,\
                      const int *HCO_arrv,\
                      const double *FCH_arrv,\
                      double *start_exp_arr,\
                      double *LC_start_exp_arr)
{
    //CONSTANTS:
    const int ITER_FOR_LIMIT_CYCLE=sim_config.ITER_FOR_LIMIT_CYCLE;
    const int MAX_PERIODS=sim_config.MAX_ALLOWED_PERIODS;

    //variables:
    double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double prev_exp_arr[NUM_NODES]; 
    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    int i;

    //storage for expressions at each peak:
    double max_exp_arr[MAX_PERIODS][NUM_NODES];
    //storage for expressions at each valley:
    double min_exp_arr[MAX_PERIODS][NUM_NODES];
    //storage for index at each peak:
    int max_idx_arr[MAX_PERIODS];
    //storage for index at each valley:
    int min_idx_arr[MAX_PERIODS];
    //store for distances from the start to each valley:
    double min_dist_arr[MAX_PERIODS];

    //storage for differences between min_dist and prev_dist:
    //double diff_dist_arr[MAX_PERIODS];
    
    //counts for peaks and valleys:
    int count_max_exp, count_min_exp;
    
    double prev_dist, curr_dist; //keeps track of distance from the start
    bool moving_uphill; //keeps track of simulation direction

    int count_exp; //global count for simulation steps
    int period; //period for the LC
    int count_iter; //count for the outer loop

    //copy the start expressions to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=start_exp_arr[i]; 
        prev_exp_arr[i]=curr_exp_arr[i]; 
        //printf("%2.6f\t",curr_exp_arr[i]);
    }
    
    //initialize the control variables:
    prev_dist=curr_dist=0.0;
    moving_uphill=true;
    count_exp=0;
    count_max_exp=count_min_exp=0;
    period=0;

    //calculate ODE-outer loop:
    for(count_iter=0; count_iter<ITER_FOR_LIMIT_CYCLE;count_iter++) 
    {
        //calculate ODE-inner loop:
        int count_step=0;
        while(count_step<SIM_STEPS)
        {
            cal_fX(NUM_NODES,NUM_EDGES,SIM_STEP_SIZE,\
                   TRANS_RATE_FACTOR,\
                   MPR_arrv,DNR_arrv,\
                   node_type_arrv,edge_source_arrv,\
                   edge_target_arrv,edge_type_arrv,\
                   TSH_arrv,HCO_arrv,FCH_arrv,\
                   curr_exp_arr,fX_arr);
            //calculate next expressions:
            for (i=0;i<NUM_NODES;++i) {
                 next_exp_arr[i]=curr_exp_arr[i]+fX_arr[i]*SIM_STEP_SIZE;
            }
 
            curr_dist=sum_delta(start_exp_arr,next_exp_arr,NUM_NODES);
            if(moving_uphill){ //while moving uphill:
                if (curr_dist<prev_dist){
                    //hit the peak already. now, switch to downhill:
                    moving_uphill=false;
                }
            }
            else { //while moving downhill:
                if (curr_dist>prev_dist){
                    //hit the bottom already. now, switch to uphill:
                    moving_uphill=true;
                    //save information about this valley:
                    for (i=0;i<NUM_NODES;++i){
                        min_exp_arr[count_min_exp][i]=curr_exp_arr[i];
                    }
                    min_idx_arr[count_min_exp]=count_exp;
                    min_dist_arr[count_min_exp]=prev_dist;
                    //diff_dist_arr[count_min_exp]=sum_delta(curr_exp_arr,\
                    //                                       next_exp_arr,\
                    //                                       NUM_NODES);
                    count_min_exp++;
                } //(curr_dist>prev_dist)
            } //end of the block for moving downhill

            if((count_max_exp>=MAX_PERIODS) ||
               (count_min_exp>=MAX_PERIODS)) break;
             
            prev_dist=curr_dist; //uphill motion
            //copy next exp to curr exp and save curr exp
            for (i=0;i<NUM_NODES;++i){
                curr_exp_arr[i]=next_exp_arr[i]>=0?next_exp_arr[i]:0;
            }
            count_exp++;
            count_step+=1;
        } //end of while(count_step<SIM_STEPS)

        double fX_norm=cal_norm(fX_arr,NUM_NODES);
        //if a stable state is found or a boundary situation is met,
        //then return -1:
        if(fX_norm<=CONVERGENCE_PROXIMITY) return -1;

        period=cal_period(NUM_NODES,\
                          MAX_PERIODS,\
                          count_min_exp,\
                          min_idx_arr,\
                          min_exp_arr,\
                          LC_start_exp_arr);
        //if nonzero period found, then break of the loop:
        if (period!=0) break;
    }//for(count_iter=0; count_iter<ITER_FOR_LIMIT_CYCLE;count_iter++)

    return period;
} 

/*-----------------------------------------------------------------------*/
int cal_period(int NUM_NODES,\
               int MAX_PERIODS,\
               int count_min_exp,\
               int *min_idx_arr,\
               double min_exp_arr[][NUM_NODES],\
               double *LC_start_exp_arr)
{
    //CALCULATE PERIODS from all valleys 
    //CONSTANTS:
    const int NO_OF_SAMPLED_PERIODS=sim_config.NO_OF_SAMPLED_PERIODS;
    const int ALLOWED_ERROR_IN_PERIODS=sim_config.ALLOWED_ERROR_IN_PERIODS;
    const double SAMEPOINT_PROXIMITY=sim_config.SAMEPOINT_PROXIMITY;

    //LOCAL VARIABLES:
    int period;
    int count_same_min_dist;
    int count_period; //counts number of periods found
    int period_T[MAX_PERIODS]; //saves the periods
    int idx_last_period; //keeps track of the index of the last period
    double diff_exp; //saves diff between expr levels at two time points
    bool found = false; //marks whether period is found

    int i;

    //initialize all variables:
    period=0;
    count_period=0; 
    diff_exp=0.0; 
    count_same_min_dist=0;   
    //set a marker at the index of the position of the last valley:
    idx_last_period=count_min_exp-1; 

    //check in backward direction starting from
    //the last element in min_dist_arr:
    for (i=count_min_exp-2;i>=0;i--){
       diff_exp=sum_delta(min_exp_arr[count_min_exp-1],\
                          min_exp_arr[i],NUM_NODES);
       if(diff_exp<=SAMEPOINT_PROXIMITY){
          count_same_min_dist++;
          period_T[count_period]=min_idx_arr[idx_last_period]-min_idx_arr[i];
          //save position of the valley where the last period was found:
          idx_last_period=i;
          count_period++; //increment count of the periods found
       }
    }
    if(count_period>=NO_OF_SAMPLED_PERIODS){
        found = true;
        //check whether each of these periods are 
        //within the allowed error limit:
        for (i=count_period-1;i>=1;i--){
            if(abs(period_T[i]-period_T[i-1])>ALLOWED_ERROR_IN_PERIODS){
                found = false;
                break;
            }
        }
    }
    if(found){
        //copy the last sampled period to period varialbe for returning:
        period=period_T[count_period-1];
        //copy the last minimum expression to the 
        //LC start variable for returning:
        for(i=0;i<NUM_NODES;i++){
           LC_start_exp_arr[i]=min_exp_arr[count_min_exp-1][i];
        }
    }
    return period;
}

/*-----------------------------------------------------------------------*/
double cal_limitcycle(int NUM_NODES,int NUM_EDGES,\
                      int SIM_STEPS,double SIM_STEP_SIZE,\
                      long double CONVERGENCE_PROXIMITY,\
                      double TRANS_RATE_FACTOR,\
                      const double *MPR_arrv,const double *DNR_arrv,\
                      const int *node_type_arrv,\
                      const int *edge_source_arrv,\
                      const int *edge_target_arrv,\
                      const int *edge_type_arrv,\
                      const double *TSH_arrv,\
                      const int *HCO_arrv,\
                      const double *FCH_arrv,\
                      double *start_exp_arr,\
                      double **LC_exp_arr)
{

    double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    int i;
    double max_dist;
    //copy the start expressions to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=start_exp_arr[i]; 
        LC_exp_arr[0][i]=curr_exp_arr[i];
    }

    int count_step=0;
    max_dist=0.0;
    while(count_step<SIM_STEPS)
    {
        cal_fX(NUM_NODES,NUM_EDGES,SIM_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,fX_arr);
        //calculate next expressions:
        for (i=0;i<NUM_NODES;++i) {
             next_exp_arr[i]=curr_exp_arr[i]+fX_arr[i]*SIM_STEP_SIZE;
        }

        double tmp_dist;
        tmp_dist=sum_delta(curr_exp_arr,next_exp_arr,NUM_NODES);
        if(tmp_dist>max_dist) max_dist=tmp_dist;

        count_step+=1;
        //copy next exp to curr exp and save curr exp
        for (i=0;i<NUM_NODES;++i){
            curr_exp_arr[i]=next_exp_arr[i]>=0?next_exp_arr[i]:0;
            LC_exp_arr[count_step][i]=curr_exp_arr[i];
        }
   }
   //printf("%f\n",max_dist);
   return max_dist;
}

/*-----------------------------------------------------------------------*/
void write_limitcycle(FILE *fh_LCs, int MODEL_NO, int NUM_NODES,\
                      int count_LC, int SIZE_OF_PERIOD, int SIZE_OF_LIMIT_CYCLE,\
                      double **LC_exp_arr)
{ 
    //This function writes the expressions of the network states 
    //along the limit cycle.
    int j; //outer loop counter
    int k; //inner loop coutner

    //for (j=0;j<NEW_SIM_STEPS+1;j++){
    for (j=0;j<SIZE_OF_LIMIT_CYCLE+1;j++){
       fprintf(fh_LCs,"%d\t",MODEL_NO);
       //write limit cycle number:
       fprintf(fh_LCs,"%d\t",count_LC);
       //write period of LC:
       fprintf(fh_LCs,"%d\t",SIZE_OF_PERIOD);
         for(k=0;k<NUM_NODES;k++){ 
            //save log2(gene expression):
            fprintf(fh_LCs,"%10.6f\t",log(LC_exp_arr[j][k])/log(2));
         }
         //new line to separate the expressions from 
         //those of the next time step: 
         fprintf(fh_LCs,"%s","\n");
    } 
    //new line after writing the whole LC
    //this separates expressions between limit cycles:
    fprintf(fh_LCs,"%s","\n");
    return;
}

/*-----------------------------------------------------------------------*/
void cal_fX(int NUM_NODES,int NUM_EDGES,\
            double EULER_STEP_SIZE,\
            double TRANS_RATE_FACTOR,\
            const double *MPR_arrv,
            const double *DNR_arrv,\
            const int *node_type_arrv,\
            const int *edge_source_arrv,\
            const int *edge_target_arrv,\
            const int *edge_type_arrv,\
            const double *TSH_arrv,\
            const int *HCO_arrv,\
            const double *FCH_arrv,\
            double const *curr_exp_arr,\
            double *fX_arr)
{
    //double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part 
    int i;

    //fX=Gx:
    //fX2=kX*X:
    for (i=0;i<NUM_NODES;++i) {
        fX_arr[i]=MPR_arrv[i]; 
        fX2_arr[i]=DNR_arrv[i]*curr_exp_arr[i]; 
    } 
    //update fX with regulations:
    for (i=0;i<NUM_EDGES;i++){
        if (edge_type_arrv[i]==1 || edge_type_arrv[i]==5) {
           fX_arr[edge_target_arrv[i]] *=\
	          eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                             	 TSH_arrv[i],HCO_arrv[i],\
                           		 FCH_arrv[i])/FCH_arrv[i]; 
        }
        else if (edge_type_arrv[i]==2 || edge_type_arrv[i]==6) {
           fX_arr[edge_target_arrv[i]] *=\
                  eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                           		 TSH_arrv[i],HCO_arrv[i],\
 				         FCH_arrv[i]);
        }
        else if (edge_type_arrv[i]==3 || edge_type_arrv[i]==4) {
           fX2_arr[edge_target_arrv[i]] *=\
                  eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                           		 TSH_arrv[i],HCO_arrv[i],\
 				         FCH_arrv[i]); 
        }
    }
    //adjust fX by kX*X:
    for (i=0;i<NUM_NODES;++i) {
        fX_arr[i]-=fX2_arr[i];
        if (node_type_arrv[i]==2){
           fX_arr[i]*=TRANS_RATE_FACTOR;
        } 
    }
    return;
}

/*-----------------------------------------------------------------------*/
void estimate_stable_expression(int NUM_NODES,int NUM_EDGES,\
                                int ITER_FOR_ODE,int TOTAL_EULER_STEPS,\
                                double EULER_STEP_SIZE,\
                                long double CONVERGENCE_PROXIMITY,\
                                double TRANS_RATE_FACTOR,\
                                const double *MPR_arrv,\
                                const double *DNR_arrv,\
                                const int *node_type_arrv,
                    		        const int *edge_source_arrv,\
                    		        const int *edge_target_arrv,\
                    		        const int *edge_type_arrv,\
                                const double *TSH_arrv,\
                    		        const int *HCO_arrv,\
                    		        const double *FCH_arrv,\
                                const double *IC_arr, double *exp_arr,\
                                double *fX_arr)
{
    //CONSTANTS:
    const int ITER_FOR_RELAXATION=sim_config.ITER_FOR_RELAXATION;

    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    double test_delta=0.0;
    int iter_count=0;
    int i=0;

    double fX_norm; 
    double fX_norm_now;

    //copy ICs to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=IC_arr[i];
    }

    while(iter_count<ITER_FOR_ODE)
    {
        cal_euler_approximation(NUM_NODES,NUM_EDGES,TOTAL_EULER_STEPS,\
                                EULER_STEP_SIZE,TRANS_RATE_FACTOR,\
                                MPR_arrv,DNR_arrv,node_type_arrv,\
                                edge_source_arrv,edge_target_arrv,\
                                edge_type_arrv,TSH_arrv,HCO_arrv,FCH_arrv,\
    			                    curr_exp_arr,next_exp_arr,fX_arr);

        test_delta=sum_delta(curr_exp_arr,next_exp_arr,NUM_NODES);
        
        //check for convergence to a state:
        //if the simulation converges to a state, then break out of the loop:
        if (test_delta<=CONVERGENCE_PROXIMITY){
            break;
        }
        fX_norm_now=cal_norm(fX_arr,NUM_NODES);
        //check for convergence to a limit cycle:
        //after ITER_RELAXATION times of iteration, check whether limit
        //cycle situation arises. If it does, break outf of the loop:
        if (iter_count>ITER_FOR_RELAXATION) 
        {
            if(fX_norm_now>fX_norm) break;
        }
        fX_norm=fX_norm_now;
       
        //update curr_exp_arr:
        for (i=0;i<NUM_NODES;++i) {
            curr_exp_arr[i]=next_exp_arr[i]; 
        }
        iter_count+=1;
    } 
    //copy the expression to exp_arr for returning:
    for (i=0;i<NUM_NODES;i++){
        exp_arr[i]=next_exp_arr[i];
    }
    //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]); 
    //printf("\n"); 
    return;
}

/*-----------------------------------------------------------------------*/
void cal_euler_approximation(int NUM_NODES,int NUM_EDGES,\
                             int TOTAL_EULER_STEPS,double EULER_STEP_SIZE,\
                             double TRANS_RATE_FACTOR,\
                             const double *MPR_arrv,const double *DNR_arrv,\
                             const int *node_type_arrv,\
                             const int *edge_source_arrv,\
                             const int *edge_target_arrv,\
                             const int *edge_type_arrv,\
                             const double *TSH_arrv,\
                             const int *HCO_arrv,\
                             const double *FCH_arrv,\
                             double const *curr_exp_arrv,\
                             double *next_exp_arr,double *fX_arr)
{

   /* 
     Given the current value of an n dimensional point, this function 
     calculates the final value of the point after TOTAL_EULER_STEPS, 
     by using the Euler approximation method.
    */
    //double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part 
    double curr_exp_arr[NUM_NODES]; //start expression (n dimensional point)
    double hill_impact=0.0;
    int fX_POINTS=(int) (TOTAL_EULER_STEPS/100);
    int i;
    //EULER_STEP_SIZE*=0.1;
//    const double FAST_EULER_STEP_SIZE=EULER_STEP_SIZE*TRANS_RATE_FACTOR; 
    double fX_norm=0.0; 

    int step_count=0; 
 
    //copy start expressions to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=curr_exp_arrv[i]; 
    }
    //fX_norm=0.0;
    while(step_count<TOTAL_EULER_STEPS)
    {
        //calculate the slope/rate at the current point: 
        cal_fX(NUM_NODES,NUM_EDGES,EULER_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,fX_arr);
        //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]); 
        //printf("\n");  
        //exit(0);
        //calculate next expressions: y=y+hF
        for (i=0;i<NUM_NODES;++i) {
             next_exp_arr[i]=curr_exp_arr[i]+fX_arr[i]*EULER_STEP_SIZE;
        }
        //update curr_exp_arr:
        //prevent curr_exp_arr value from getting -ve value:
        //set it to zero, if negative exp value is encountered:
        for (i=0;i<NUM_NODES;++i) {
            //curr_exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i]>=0?next_exp_arr[i]:0; 
        }
        step_count+=1;
        
        /*
	if ((step_count%fX_POINTS)==0)
        {
            fX_norm+=cal_norm(fX_arr,NUM_NODES);
        }
        */
    }
    //printf("%f\n",fX_norm);
    return;
}

/* ---- ODE simulation  ----- */
void cal_ode(int NUM_NODES,int NUM_EDGES,\
			    int TOTAL_EULER_STEPS,double EULER_STEP_SIZE,\
             double TRANS_RATE_FACTOR,\
             const double *MPR_arrv,const double *DNR_arrv,\
             const int *node_type_arrv,\
             const int *edge_source_arrv,\
             const int *edge_target_arrv,\
             const int *edge_type_arrv,\
             const double *TSH_arrv,\
             const int *HCO_arrv,\
             const double *FCH_arrv,\
             double const *curr_exp_arrv)
{

    double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    //double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part 
    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    double hill_impact=0.0; 
    int fX_POINTS=(int) (TOTAL_EULER_STEPS/100);
    int i;
    //EULER_STEP_SIZE*=0.1;
    //const double FAST_EULER_STEP_SIZE=EULER_STEP_SIZE*TRANS_RATE_FACTOR; 
    double fX_norm; 
    

    int step_count=0; 
    FILE *fh_out=NULL; 
    char fname[100]="out.txt";

    fh_out=fopen(fname,"w");
    if (fh_out==NULL) { 
        printf("error opening file\n"); 
        exit(0);
    } 
    //copy start expressions to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=curr_exp_arrv[i]; 
    }
    fX_norm=0.0;
    while(step_count<TOTAL_EULER_STEPS)
    {
        cal_fX(NUM_NODES,NUM_EDGES,EULER_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,fX_arr);
        //calculate next expressions:
        for (i=0;i<NUM_NODES;++i) {
             next_exp_arr[i]=curr_exp_arr[i]+fX_arr[i]*EULER_STEP_SIZE;
        }
        //update curr_exp_arr:
        //prevent curr_exp_arr value from getting -ve value:
        //set it to zero, if negative exp value is encountered:
        step_count+=1;
	fprintf(fh_out,"%d\t",step_count);
        for (i=0;i<NUM_NODES;++i) {
            //curr_exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i]>=0?next_exp_arr[i]:0; 
	    fprintf(fh_out,"%f\t", curr_exp_arr[i]);
        }
	fprintf(fh_out,"\n");
    }
    fclose(fh_out);
    return;
}

/*-----------------------------------------------------------------------*/
void cal_EXPrange(int NUM_NODES,int NUM_EDGES,const double *MPR_arrv,\
             const double *DNR_arrv,const int *edge_source_arrv,\
             const int *edge_target_arrv,const int *edge_type_arrv,\
             const double *FCH_arrv,double *minEXP_arr, double *maxEXP_arr)
/*------------------------------------------------------------------
 It calculates gene expression range to be used in IC estimation 
considering the edge type and FCH.

It calculates the allowable minimum and maximum expression for each 
node, based on the following rules:
   (1) Initially, set the minimum and maximum expressions of each 
       node based on MPR and DNR:
       minEXP(X) <- MPR(X)/DNR(X)
       maxEXP(X) <- minEXP(X)

   (2) Then, update these values based on regulator type
       for each regulator Y (edge Y->X):
         if edge type 1 or 5: 
            minEXP(X) <- minEXP(X)/FCH(Y->X)
         if edge type 2 or 6: 
            minEXP(X) <- minEXP(X)*FCH(Y->X)
         if edge type 4: 
            minEXP(X) <- FCH(Y->X)
         if edge type 3: 
            minEXP(X) <- FCH(Y->X)
--------------------------------------------------------------------*/
{
    int i;
    
    //set minEXP_arr and maxEXP_arr based on MPR and DNR:
    for (i=0;i<NUM_NODES;++i) {
        minEXP_arr[i]=MPR_arrv[i]/DNR_arrv[i]; 
        maxEXP_arr[i]=minEXP_arr[i];
    }

    //update minEXP_arr based on type of regulators:
    for (i=0;i<NUM_EDGES;i++){
        //if (edge_type_arrv[i]==REG_EXCITATORY) {
        if (edge_type_arrv[i]==1 || edge_type_arrv[i]==5) {
           minEXP_arr[edge_target_arrv[i]] /= FCH_arrv[i];
        }
        //else if (edge_type_arrv[i]==REG_INHIBITORY) {
        else if (edge_type_arrv[i]==2 || edge_type_arrv[i]==6) {
           minEXP_arr[edge_target_arrv[i]] *= FCH_arrv[i];
        }
        else if (edge_type_arrv[i]==4) {
           minEXP_arr[edge_target_arrv[i]] /= FCH_arrv[i];
        }
        else if (edge_type_arrv[i]==3) {
           maxEXP_arr[edge_target_arrv[i]] /= FCH_arrv[i];
        }
    }
    return;
}

/*-----------------------------------------------------------------------*/
void set_ICs(int NUM_NODES, const double *minEXP_arr, \
             const double *maxEXP_arr, double *IC_arr)
{

    //double maxEXP_arr[NUM_NODES],minEXP_arr[NUM_NODES];
    int i;

    //generate random ICs within a range:
    for (i=0;i<NUM_NODES;++i) {
        IC_arr[i]=exp2(randu(log2(minEXP_arr[i]),log2(maxEXP_arr[i])));
    }
    return;
}

/*-----------------------------------------------------------------------*/
double eval_shiftedHill_fn (double x, double x0, double nx, double lamda)
{
  return (lamda+(1.0-lamda)*(1.0/(1.0+pow((x/x0),nx))));
}

/*-----------------------------------------------------------------------*/
double sum_delta (double *curr_exp, double *next_exp, int NUM_NODES)
{
    int i = 0;
    double ssq = 0.0;
    
    for (i=0;i<NUM_NODES;i++){
        ssq+=pow(curr_exp[i]-next_exp[i],2);
    }
    return ssq;
}

/*-----------------------------------------------------------------------*/
double cal_norm (double *curr_exp, int NUM_NODES)
{
    int i = 0;
    double ssq = 0.0;
    
    for (i=0;i<NUM_NODES;i++){
        ssq+=pow(curr_exp[i],2);
    }
    return ssq;
}

/*-----------------------------------------------------------------------*/
bool is_element(int key, const int *arr, int size) 
{ 
    int i;
    for (i=0;i<size;i++) 
    {
        if (key== arr[i]) 
            return true; 
    }
    return false; 
}

/*-----------------------------------------------------------------------*/
double **allocate_LCmemory(int LC_SIZE, int NUM_NODES)
{
    //allocates memory for limit cycle:
    double **EXP_arr; 
    int i,j; 
    EXP_arr=(double **) calloc(LC_SIZE,sizeof(double));
    for (i=0; i<LC_SIZE;i++){ 
       EXP_arr[i]=(double *) calloc(NUM_NODES,sizeof(double));
    }

    for (i=0; i<LC_SIZE;i++){ 
        for(j=0;j<NUM_NODES;j++){
           EXP_arr[i][j]=0.0;
        }
    }
    return EXP_arr;
}

/*-----------------------------------------------------------------------*/
void free_LCmemory(double **EXP_arr, int LC_SIZE, int NUM_NODES)
{
    //free the allocated memory:
    int i;
    //first free allocated memory for gene expression for each node:
    for (i=0; i<LC_SIZE;i++){
       free(EXP_arr[i]);
    }
    //second, free memory allocated for holding the memory address:
    free(EXP_arr);
}

/*-----------------------------------------------------------------------*/
void read_config(char *WORK_DIR, char* fname_cfg)
{ 
    //reads the deafult values from the config file 

    char fname[300]="\0";
    strcat(fname,WORK_DIR);
    strcat(fname,"/");
    strcat(fname,fname_cfg);

    FILE *fh_config;
    fh_config=fopen(fname,"r");

    char* line=malloc(LINE_SIZE); 
    char delim[5];
    char *token; 
    strcpy(delim,CONFIG_FILE_DELIM);
    //printf("%s\n",CONFIG_FILE_DELIM);
    //printf("within read config\n");

    while(fgets(line,LINE_SIZE,fh_config)!=NULL) {
         token=strtok(line,delim);
         if(strcmp(token,"NUM_RANDOM_ICS")==0){
               token=strtok(NULL,"\n"); 
               sim_config.NUM_RANDOM_ICS=strtol(token,NULL,10);
         }
         else if(strcmp(token,"CONVERGENCE_PROXIMITY")==0){
               token=strtok(NULL,"\n"); 
               sim_config.CONVERGENCE_PROXIMITY=strtod(token,NULL);
         }
         else if(strcmp(token,"SAMEPOINT_PROXIMITY")==0){
               token=strtok(NULL,"\n"); 
               sim_config.SAMEPOINT_PROXIMITY=strtod(token,NULL);
         }
         else if(strcmp(token,"ITER_FOR_ODE")==0){
               token=strtok(NULL,"\n"); 
               sim_config.ITER_FOR_ODE=strtol(token,NULL,10);
         }
         else if(strcmp(token,"ITER_FOR_RELAXATION")==0){
               token=strtok(NULL,"\n"); 
               sim_config.ITER_FOR_RELAXATION=strtol(token,NULL,10);
         }
         else if(strcmp(token,"ITER_FOR_LIMIT_CYCLE")==0){
               token=strtok(NULL,"\n"); 
               sim_config.ITER_FOR_LIMIT_CYCLE=strtol(token,NULL,10);
         }
         else if(strcmp(token,"EULER_SIM_TIME")==0){
               token=strtok(NULL,"\n"); 
               sim_config.EULER_SIM_TIME=strtod(token,NULL);
         }
         else if(strcmp(token,"EULER_SIM_TIME_LAST_ANN_STEP")==0){
               token=strtok(NULL,"\n"); 
               sim_config.EULER_SIM_TIME_LAST_ANN_STEP=strtod(token,NULL);
         }
         else if(strcmp(token,"STARTING_NOISE_ANN")==0){
               token=strtok(NULL,"\n"); 
               sim_config.STARTING_NOISE_ANN=strtod(token,NULL);
         }
         else if(strcmp(token,"STARTING_NOISE_CONSTANT")==0){
               token=strtok(NULL,"\n"); 
               sim_config.STARTING_NOISE_CONSTANT=strtod(token,NULL);
         }
         else if(strcmp(token,"NOISE_SCALING_FACTOR_ANN")==0){
               token=strtok(NULL,"\n"); 
               sim_config.NOISE_SCALING_FACTOR_ANN=strtod(token,NULL);
         }
         else if(strcmp(token,"LIMIT_CYCLE_SIM_TIME")==0){
               token=strtok(NULL,"\n"); 
               sim_config.LIMIT_CYCLE_SIM_TIME=strtod(token,NULL);
         }
         else if(strcmp(token,"EULER_SIM_STEP_SIZE")==0){
               token=strtok(NULL,"\n"); 
               sim_config.EULER_SIM_STEP_SIZE=strtod(token,NULL);
         }
         else if(strcmp(token,"LIMIT_CYCLE_SIM_STEP_SIZE")==0){
               token=strtok(NULL,"\n"); 
               sim_config.LIMIT_CYCLE_SIM_STEP_SIZE=strtod(token,NULL);
         }
         else if(strcmp(token,"MAX_STABLE_STATES")==0){
               token=strtok(NULL,"\n"); 
               sim_config.MAX_STABLE_STATES=strtol(token,NULL,10);
         }
         else if(strcmp(token,"MAX_LIMIT_CYCLES")==0){
               token=strtok(NULL,"\n"); 
               sim_config.MAX_LIMIT_CYCLES=strtol(token,NULL,10);
         }
         else if(strcmp(token,"MAX_ALLOWED_PERIODS")==0){
               token=strtok(NULL,"\n"); 
               sim_config.MAX_ALLOWED_PERIODS=strtol(token,NULL,10);
         }
         else if(strcmp(token,"NO_OF_SAMPLED_PERIODS")==0){
               token=strtok(NULL,"\n"); 
               sim_config.NO_OF_SAMPLED_PERIODS=strtol(token,NULL,10);
         }
         else if(strcmp(token,"ALLOWED_ERROR_IN_PERIODS")==0){
               token=strtok(NULL,"\n"); 
               sim_config.ALLOWED_ERROR_IN_PERIODS=strtol(token,NULL,10);
         }
    }
    fclose(fh_config);
    free(line);
    return;
} 

/*-----------------------------------------------------------------------*/
void print_config()
{ 
    printf("NUM_RANDOM_ICS: %d\n",sim_config.NUM_RANDOM_ICS);
    printf("CONVERGENCE_PROXIMITY: %LE\n",sim_config.CONVERGENCE_PROXIMITY);
    printf("SAMEPOINT_PROXIMITY: %f\n",sim_config.SAMEPOINT_PROXIMITY);
    printf("ITER_FOR_ODE: %d\n",sim_config.ITER_FOR_ODE);
    printf("ITER_FOR_RELAXATION: %d\n",sim_config.ITER_FOR_RELAXATION);
    printf("ITER_FOR_LIMIT_CYCLE: %d\n",sim_config.ITER_FOR_LIMIT_CYCLE);
    printf("EULER_SIM_TIME: %f\n",sim_config.EULER_SIM_TIME);
    printf("LIMIT_CYCLE_SIM_TIME: %f\n",sim_config.LIMIT_CYCLE_SIM_TIME);
    printf("EULER_SIM_STEP_SIZE: %f\n",sim_config.EULER_SIM_STEP_SIZE);
    printf("LIMIT_CYCLE_SIM_STEP_SIZE: %f\n",sim_config.LIMIT_CYCLE_SIM_STEP_SIZE);
    printf("MAX_STABLE_STATES: %d\n",sim_config.MAX_STABLE_STATES);
    printf("MAX_LIMIT_CYCLES: %d\n",sim_config.MAX_LIMIT_CYCLES);
    printf("MAX_ALLOWED_PERIODS: %d\n",sim_config.MAX_ALLOWED_PERIODS);
    printf("NO_OF_SAMPLED_PERIODS: %d\n",sim_config.NO_OF_SAMPLED_PERIODS);
    printf("ALLOWED_ERROR_IN_PERIODS: %d\n",sim_config.ALLOWED_ERROR_IN_PERIODS);
    return;
}


/********* ------- COMMON FUNCTIONS for stochastic simulation ------- ********/
/*-----------------------------   START   -----------------------------------*/
void write_state_expressions(FILE *fh_states, int MODEL_NO, \
                             int NUM_NODES, double *exp_arr)
{
    const int STATE_NO=1; 
    const int NUM_STATES=1;
   //write model index, number of states, and state number:
    fprintf(fh_states,"%d\t%d\t%d\t",MODEL_NO, NUM_STATES, STATE_NO);  

    int i;
    for (i=0;i<NUM_NODES;i++){
        fprintf(fh_states,"%10.6f\t",log(exp_arr[i])/log(2));
    }
    fprintf(fh_states,"\n");
    return;
} 

void write_state_expressions_annealed(FILE *fh_states, int MODEL_NO, \
                                     double NOISE_LEVEL, int NUM_NODES, \
                                     double *exp_arr)
{
    const int STATE_NO=1; 
    const int NUM_STATES=1;
   //write model index, number of states, and state number:
    fprintf(fh_states,"%d\t%10.6f\t%d\t",MODEL_NO, NOISE_LEVEL, STATE_NO);  

    int i;
    for (i=0;i<NUM_NODES;i++){
        fprintf(fh_states,"%10.6f\t",log(exp_arr[i])/log(2));
    }
    fprintf(fh_states,"\n");
    return;
}


/*************************** ----- COMMON FUNCTIONS ----- ********************/
/*---------------------------------   END   ---------------------------------*/



/* probing models using stochastic simulation with GNW constant noise */
/*---------------------------- START ------------------------------------*/
void simulate_network_cnoise_gnw(char *WORK_DIR,\
                             char *FNAME_STATES, \
                             char *FNAME_LIMITCYCLES, \
                             char *FNAME_SUMMARY, \
                             const int MODEL_NO,\
                             int NUM_NODES,\
                             int NUM_EDGES,\
                             int NUM_RANDOM_ICS,\
                             int ITER_FOR_ODE,\
                             double EULER_SIM_TIME,\
                             double EULER_STEP_SIZE,\
                             long double CONVERGENCE_PROXIMITY,\
                             double TRANS_RATE_FACTOR,\
                             const double *MPR_arrv,\
                             const double *DNR_arrv,\
                             const int *node_type_arrv,\
                             const int *edge_source_arrv,\
                             const int *edge_target_arrv,\
                             const int *edge_type_arrv,\
                             const double *TSH_arrv,\
                             const int *HCO_arrv,\
                             const double *FCH_arrv,\
                             double *EXP_dict_arrv)
{
    //printf("Stochastic model simulation!\n");
    //printf("%f\n",NOISE);
    //CONSTANTS:
    const int TOTAL_EULER_STEPS=(int) (EULER_SIM_TIME/EULER_STEP_SIZE);
    
    //VARIABLES:
    double NOISE=0.0;
    double *EXP_dict_arr=(double *) EXP_dict_arrv;
    double IC_arr[NUM_NODES]; //saves ICs generated by set_ICs()
    double exp_arr[NUM_NODES]; //saves exps from estimate_stable_expression()
    
    double minEXP_arr[NUM_NODES]; //minimum expression for all nodes
    double maxEXP_arr[NUM_NODES]; //maximum expression for all nodes
    
    double EXP_arr[NUM_RANDOM_ICS][NUM_NODES];
    double fX_arr_norm[NUM_RANDOM_ICS];
    double fX_arr[NUM_NODES];
    double zero_arr[NUM_NODES];
    
    //INDEX VARIABLES:
    int count_pos=0; //keeps track of position in EXP_dict_arr
    int count_iteration; //counts iteration over IC initialization
    int i=0; //general purpose
    
    //read the simulation configations:
    char fname_cfg[30]="racipe.cfg";
    
    read_config(WORK_DIR,fname_cfg);
    //printf("%d\n",sim_config.NUM_RANDOM_ICS);
    //printf("%f\n",sim_config.EULER_SIM_TIME);
    //printf("%f\n",sim_config.STARTING_NOISE_CONSTANT);
    
    NOISE = sim_config.STARTING_NOISE_CONSTANT;

    cal_EXPrange(NUM_NODES,NUM_EDGES,MPR_arrv,DNR_arrv,\
                 edge_source_arrv, edge_target_arrv,\
                 edge_type_arrv,FCH_arrv,minEXP_arr, maxEXP_arr);
    
    //set exp_arr to zeros:
    for (i=0;i<NUM_NODES;++i) {
        exp_arr[i]=0.0;
        fX_arr[i]=0.0;
        zero_arr[i]=0.0;
    }
    
    //SAVE SOLUTIONS: solution_arr
    FILE *fh_states;
    //char str_noise[5];
    //snprintf(str_noise, 5, "%lf", NOISE);

    //open file for appending expressions:
    fh_states=fopen(FNAME_STATES, "a");
    //fprintf(fh_states,"%d\t",MODEL_NO);  //write model index
    
    count_iteration=0;
    while(count_iteration<NUM_RANDOM_ICS)
    {
        set_ICs(NUM_NODES, minEXP_arr, maxEXP_arr, IC_arr);
     
        estimate_stable_expression_cnoise_gnw(NUM_NODES, NUM_EDGES, ITER_FOR_ODE,\
                                   TOTAL_EULER_STEPS, EULER_STEP_SIZE,\
                                   CONVERGENCE_PROXIMITY,\
                                   TRANS_RATE_FACTOR,\
                                   MPR_arrv, DNR_arrv,\
                                   node_type_arrv,\
                                   edge_source_arrv,\
                                   edge_target_arrv,\
                                   edge_type_arrv,\
                                   TSH_arrv, HCO_arrv, FCH_arrv,\
                                   IC_arr, exp_arr, fX_arr, \
                                   NOISE);
       
        /*
        estimate_stable_expression(NUM_NODES, NUM_EDGES, ITER_FOR_ODE,\
                                   TOTAL_EULER_STEPS, EULER_STEP_SIZE,\
                                   CONVERGENCE_PROXIMITY,\
                                   TRANS_RATE_FACTOR,\
                                   MPR_arrv, DNR_arrv,\
                                   node_type_arrv,\
                                   edge_source_arrv,\
                                   edge_target_arrv,\
                                   edge_type_arrv,\
                                   TSH_arrv, HCO_arrv, FCH_arrv,\
                                   IC_arr, exp_arr, fX_arr);
      */ 


       //write_state_expressions(fh_states, MODEL_NO, NUM_NODES, exp_arr);

        //fX_arr_norm[count_iteration]=cal_norm(fX_arr,NUM_NODES);
        
        //write model index, number of states, state number
        fprintf(fh_states,"%d\t%d\t%d\t",MODEL_NO, 1, 1);  

        //copy the exp_arr to EXP_dict_arr:
        //if any value is found <= zero, assign it min value
        for(i=0;i<NUM_NODES;i++){
            EXP_dict_arr[count_pos+i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            EXP_arr[count_iteration][i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            fprintf(fh_states,"%10.6f\t",log(EXP_arr[count_iteration][i])/log(2));
        }

        count_pos+=NUM_NODES;
        count_iteration+=1;
        fprintf(fh_states,"\n");
    } //end of while(count_iteration<NUM_RANDOM_ICS)

    fflush(fh_states);
    fclose(fh_states);
    
    return;
}

void estimate_stable_expression_cnoise_gnw(int NUM_NODES,int NUM_EDGES,\
                                           int ITER_FOR_ODE,int TOTAL_EULER_STEPS,\
                                           double EULER_STEP_SIZE,\
                                           long double CONVERGENCE_PROXIMITY,\
                                           double TRANS_RATE_FACTOR,\
                                           const double *MPR_arrv,\
                                           const double *DNR_arrv,\
                                           const int *node_type_arrv,
                                           const int *edge_source_arrv,\
                                           const int *edge_target_arrv,\
                                           const int *edge_type_arrv,\
                                           const double *TSH_arrv,\
                                           const int *HCO_arrv,\
                                           const double *FCH_arrv,\
                                           const double *IC_arr, double *exp_arr,\
                                           double *fX_arr,\
                                           double NOISE)
{
    //CONSTANTS:
    const int ITER_FOR_RELAXATION=sim_config.ITER_FOR_RELAXATION;
    
    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    double test_delta=0.0;
    int iter_count=0;
    int i=0;
    
    double fX_norm;
    double fX_norm_now;
    //printf("%f\n",NOISE);
    //copy ICs to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=IC_arr[i];
    }

    cal_euler_approximation_cnoise_gnw(NUM_NODES, NUM_EDGES, TOTAL_EULER_STEPS,\
                                       EULER_STEP_SIZE,TRANS_RATE_FACTOR,\
                                       MPR_arrv,DNR_arrv,node_type_arrv,\
                                       edge_source_arrv,edge_target_arrv,\
                                       edge_type_arrv,TSH_arrv,HCO_arrv,FCH_arrv,\
                                       curr_exp_arr,next_exp_arr,fX_arr,\
                                       NOISE);

    //update curr_exp_arr:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=next_exp_arr[i];
    }
    iter_count+=1;

    //copy the expression to exp_arr for returning:
    for (i=0;i<NUM_NODES;i++){
        exp_arr[i]=next_exp_arr[i];
    }

    //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
    //printf("\n"); 
    return;
}

void cal_euler_approximation_cnoise_gnw(int NUM_NODES,int NUM_EDGES,\
                             int TOTAL_EULER_STEPS,double EULER_STEP_SIZE,\
                             double TRANS_RATE_FACTOR,\
                             const double *MPR_arrv,\
                            const double *DNR_arrv,\
                             const int *node_type_arrv,\
                             const int *edge_source_arrv,\
                             const int *edge_target_arrv,\
                             const int *edge_type_arrv,\
                             const double *TSH_arrv,\
                             const int *HCO_arrv,\
                             const double *FCH_arrv,\
                             double const *curr_exp_arrv,\
                             double *next_exp_arr,\
                            double *fX_arr,\
                            double NOISE)
{
    //double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part
    double curr_exp_arr[NUM_NODES];
    double hill_impact=0.0;
    int fX_POINTS=(int) (TOTAL_EULER_STEPS/100);
    int i;
    double sqrt_EULER_STEP_SIZE=sqrt(EULER_STEP_SIZE);
    //EULER_STEP_SIZE*=0.1;
    //    const double FAST_EULER_STEP_SIZE=EULER_STEP_SIZE*TRANS_RATE_FACTOR;
    double fX_norm;
    //printf("%s\n", "NOISE");
    //printf("%f\n",NOISE);
    
    int step_count=0;
    
    //copy start expressions to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=curr_exp_arrv[i];
    }
    //fX_norm=0.0;
    while(step_count<TOTAL_EULER_STEPS)
    {
        cal_fX_gnw(NUM_NODES,NUM_EDGES,EULER_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr, fX_arr, fX2_arr);
        //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
        //printf("\n");
        //exit(0);
        //calculate next expressions:
        for (i=0;i<NUM_NODES;++i) {
            next_exp_arr[i]=curr_exp_arr[i]+fX_arr[i]*EULER_STEP_SIZE + \
                            NOISE*(sqrt(fX_arr[i]+fX2_arr[i])*pcg32_gasdev() + \
                                  sqrt(fX2_arr[i])*pcg32_gasdev())* \
                            sqrt_EULER_STEP_SIZE;
        }
        //update curr_exp_arr:
        //prevent curr_exp_arr value from getting -ve value:
        //set it to zero, if negative exp value is encountered:
        for (i=0;i<NUM_NODES;++i) {
            //curr_exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i]>=0?next_exp_arr[i]:0;
        }
        step_count+=1;
    }
    return;
}

/* probing models using stochastic simulation with GNW constant noise */
/*---------------------------- END ------------------------------------*/


/****** ----- stochastic simulation with constant noise ----- ************/
/*---------------------------- START ------------------------------------*/
void find_solutions_stochastic(char *WORK_DIR, const int MODEL_NO,\
                    int NUM_NODES,int NUM_EDGES,\
                    int NUM_RANDOM_ICS,int ITER_FOR_ODE,\
                    double EULER_SIM_TIME,\
                    double EULER_STEP_SIZE,\
                    long double CONVERGENCE_PROXIMITY,\
                    double TRANS_RATE_FACTOR,\
                    const double *MPR_arrv,\
                    const double *DNR_arrv,\
                    const int *node_type_arrv,\
                    const int *edge_source_arrv,\
                    const int *edge_target_arrv,\
                    const int *edge_type_arrv,\
                    const double *TSH_arrv,\
                    const int *HCO_arrv,\
                    const double *FCH_arrv,\
                    double *EXP_dict_arrv,\
                    double *NOISE_strength,\
                    double *NOISE_strength_shot,\
                    double NOISE, double NOISE_SHOT)
{
    //printf("Stochastic model simulation!\n");
    //printf("%f\n",NOISE);
    //CONSTANTS:
    const int TOTAL_EULER_STEPS=(int) (EULER_SIM_TIME/EULER_STEP_SIZE);
    
    //VARIABLES:
    double *EXP_dict_arr=(double *) EXP_dict_arrv;
    double IC_arr[NUM_NODES]; //saves ICs generated by set_ICs()
    double exp_arr[NUM_NODES]; //saves exps from estimate_stable_expression()
    
    double minEXP_arr[NUM_NODES]; //minimum expression for all nodes
    double maxEXP_arr[NUM_NODES]; //maximum expression for all nodes
    
    double EXP_arr[NUM_RANDOM_ICS][NUM_NODES];
    double fX_arr_norm[NUM_RANDOM_ICS];
    double fX_arr[NUM_NODES];
    double zero_arr[NUM_NODES];
    
    //INDEX VARIABLES:
    int count_pos=0; //keeps track of position in EXP_dict_arr
    int count_iteration; //counts iteration over IC initialization
    int i=0; //general purpose
    
    //read the simulation configations:
    char fname_cfg[30]="racipe.cfg";
    
    read_config(WORK_DIR,fname_cfg);
    
    cal_EXPrange(NUM_NODES,NUM_EDGES,MPR_arrv,DNR_arrv,\
                 edge_source_arrv, edge_target_arrv,\
                 edge_type_arrv,FCH_arrv,minEXP_arr, maxEXP_arr);
    
    //set exp_arr to zeros:
    for (i=0;i<NUM_NODES;++i) {
        exp_arr[i]=0.0;
        fX_arr[i]=0.0;
        zero_arr[i]=0.0;
    }
    
    //SAVE SOLUTIONS: solution_arr
    FILE *fh_states;
    char str_noise[5];
    snprintf(str_noise, 5, "%f", NOISE);
    char str[80]="Stochastic_NOISE";
    strcat(str,str_noise);
    strcat(str,".txt");
    fh_states=fopen(str,"a");
    fprintf(fh_states,"%d\t",MODEL_NO);  //write model index
    
        
    count_iteration=0;
    while(count_iteration<NUM_RANDOM_ICS)
    {
        set_ICs(NUM_NODES,minEXP_arr,maxEXP_arr,IC_arr);
        estimate_stable_expression_stochastic(NUM_NODES,NUM_EDGES,ITER_FOR_ODE,\
                                   TOTAL_EULER_STEPS,EULER_STEP_SIZE,\
                                   CONVERGENCE_PROXIMITY,\
                                   TRANS_RATE_FACTOR,\
                                   MPR_arrv,DNR_arrv,node_type_arrv,
                                   edge_source_arrv,\
                                   edge_target_arrv,edge_type_arrv,\
                                   TSH_arrv,HCO_arrv,FCH_arrv,\
                                   IC_arr,exp_arr,fX_arr, NOISE_strength,NOISE_strength_shot,\
                                              NOISE,NOISE_SHOT);
        
        //fX_arr_norm[count_iteration]=cal_norm(fX_arr,NUM_NODES);
        
        //copy the exp_arr to EXP_dict_arr:
        //if any value is found <= zero, assign it min value
        for(i=0;i<NUM_NODES;i++){
            EXP_dict_arr[count_pos+i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            EXP_arr[count_iteration][i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            fprintf(fh_states,"%10.6f\t",log(EXP_arr[count_iteration][i])/log(2));
            
                
        }
        count_pos+=NUM_NODES;
        count_iteration+=1;
        fprintf(fh_states,"\n");
    } //end of while(count_iteration<NUM_RANDOM_ICS)

    
    fflush(fh_states);
    fclose(fh_states);
    
    return;
}

/**********************************************************************/
void estimate_stable_expression_stochastic(int NUM_NODES,int NUM_EDGES,\
                                int ITER_FOR_ODE,int TOTAL_EULER_STEPS,\
                                double EULER_STEP_SIZE,\
                                long double CONVERGENCE_PROXIMITY,\
                                double TRANS_RATE_FACTOR,\
                                const double *MPR_arrv,\
                                const double *DNR_arrv,\
                                const int *node_type_arrv,
                                const int *edge_source_arrv,\
                                const int *edge_target_arrv,\
                                const int *edge_type_arrv,\
                                const double *TSH_arrv,\
                                const int *HCO_arrv,\
                                const double *FCH_arrv,\
                                const double *IC_arr, double *exp_arr,\
                                double *fX_arr,\
                                double *NOISE_strength,\
                                double *NOISE_strength_shot,\
                                           double NOISE,\
                                           double NOISE_SHOT)
{
    //CONSTANTS:
    const int ITER_FOR_RELAXATION=sim_config.ITER_FOR_RELAXATION;
    
    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    double test_delta=0.0;
    int iter_count=0;
    int i=0;
    
    double fX_norm;
    double fX_norm_now;
    //printf("%f\n",NOISE);
    //copy ICs to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=IC_arr[i];
    }

    cal_euler_approximation_stochastic(NUM_NODES, NUM_EDGES, TOTAL_EULER_STEPS,\
                                       EULER_STEP_SIZE, TRANS_RATE_FACTOR,\
                                       MPR_arrv, DNR_arrv, node_type_arrv,\
                                       edge_source_arrv, edge_target_arrv,\
                                       edge_type_arrv, TSH_arrv, HCO_arrv, FCH_arrv,\
                                       curr_exp_arr, next_exp_arr, fX_arr, NOISE_strength,\
                                       NOISE_strength_shot, NOISE, NOISE_SHOT);


    //update curr_exp_arr:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=next_exp_arr[i];
    }
    iter_count+=1;

    //copy the expression to exp_arr for returning:
    for (i=0;i<NUM_NODES;i++){
        exp_arr[i]=next_exp_arr[i];
    }

    //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
    //printf("\n"); 
    return;
}

/****** ----- stochastic simulation with constant noise ----- ************/
/*---------------------------- END ------------------------------------*/






/**********----- stochastic annealing using old method ----- ***********/
/*---------------------------- START ------------------------------------*/
void find_solutions_stochastic_annealing(char *WORK_DIR, const int MODEL_NO,\
                               int NUM_NODES,int NUM_EDGES,\
                               int NUM_RANDOM_ICS,int ITER_FOR_ODE,\
                               double EULER_SIM_TIME,\
                               double EULER_STEP_SIZE,\
                               long double CONVERGENCE_PROXIMITY,\
                               double TRANS_RATE_FACTOR,\
                               const double *MPR_arrv,\
                               const double *DNR_arrv,\
                               const int *node_type_arrv,\
                               const int *edge_source_arrv,\
                               const int *edge_target_arrv,\
                               const int *edge_type_arrv,\
                               const double *TSH_arrv,\
                               const int *HCO_arrv,\
                               const double *FCH_arrv,\
                               double *EXP_dict_arrv,\
                               double *NOISE_strength,\
                               double *NOISE_strength_shot,\
                               double NOISE, double NOISE_SHOT, double SCALING)
{
    //printf("Stochastic model annealing simulation!\n");
    //CONSTANTS:
    const int TOTAL_EULER_STEPS=(int) (EULER_SIM_TIME/EULER_STEP_SIZE);
    
    //VARIABLES:
    double *EXP_dict_arr=(double *) EXP_dict_arrv;
    double IC_arr[NUM_NODES]; //saves ICs generated by set_ICs()
    double exp_arr[NUM_NODES]; //saves exps from estimate_stable_expression()
    
    double minEXP_arr[NUM_NODES]; //minimum expression for all nodes
    double maxEXP_arr[NUM_NODES]; //maximum expression for all nodes
    
    double EXP_arr[NUM_RANDOM_ICS][NUM_NODES];
    double fX_arr_norm[NUM_RANDOM_ICS];
    double fX_arr[NUM_NODES];
    double zero_arr[NUM_NODES];
    
    //INDEX VARIABLES:
    int count_pos=0; //keeps track of position in EXP_dict_arr
    int count_iteration; //counts iteration over IC initialization
    int i=0; //general purpose
    
    //read the simulation configations:
    char fname_cfg[30]="racipe.cfg";
    
    read_config(WORK_DIR,fname_cfg);
    
    cal_EXPrange(NUM_NODES,NUM_EDGES,MPR_arrv,DNR_arrv,\
                 edge_source_arrv, edge_target_arrv,\
                 edge_type_arrv,FCH_arrv,minEXP_arr, maxEXP_arr);
        
    //set exp_arr to zeros:
    for (i=0;i<NUM_NODES;++i) {
        exp_arr[i]=0.0;
        fX_arr[i]=0.0;
        zero_arr[i]=0.0;
    }
    
    //SAVE SOLUTIONS: solution_arr
    //FILE *fh_states;
    //fh_states=fopen("Stochastic.txt","a");
    //fprintf(fh_states,"%d\t",MODEL_NO);  //write model index
    
    count_iteration=0;
    while(count_iteration<NUM_RANDOM_ICS)
    {
        set_ICs(NUM_NODES,minEXP_arr,maxEXP_arr,IC_arr);
        estimate_stable_expression_stochastic_annealing(MODEL_NO,NUM_NODES,\
                                                        NUM_EDGES,ITER_FOR_ODE,\
                                              TOTAL_EULER_STEPS,EULER_STEP_SIZE,\
                                              CONVERGENCE_PROXIMITY,\
                                              TRANS_RATE_FACTOR,\
                                              MPR_arrv,DNR_arrv,node_type_arrv,
                                              edge_source_arrv,\
                                              edge_target_arrv,edge_type_arrv,\
                                              TSH_arrv,HCO_arrv,FCH_arrv,\
                                              IC_arr,exp_arr,fX_arr, NOISE_strength,NOISE_strength_shot,\
                                              NOISE,NOISE_SHOT,SCALING);
        
        //fX_arr_norm[count_iteration]=cal_norm(fX_arr,NUM_NODES);
        
        //copy the exp_arr to EXP_dict_arr:
        //if any value is found <= zero, assign it min value
        for(i=0;i<NUM_NODES;i++){
            EXP_dict_arr[count_pos+i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            EXP_arr[count_iteration][i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            //fprintf(fh_states,"%10.6f\t",log(EXP_arr[count_iteration][i])/log(2));
        }
        count_pos+=NUM_NODES;
        count_iteration+=1;
        //fprintf(fh_states,"\n");
    } //end of while(count_iteration<NUM_RANDOM_ICS)
    
    
    //fflush(fh_states);
    //fclose(fh_states);
    
    return;
}

/*----------------------stochastic annealing-----------------------*/
void estimate_stable_expression_stochastic_annealing(const int MODEL_NO,\
                                                     int NUM_NODES,int NUM_EDGES,\
                                           int ITER_FOR_ODE,int TOTAL_EULER_STEPS,\
                                           double EULER_STEP_SIZE,\
                                           long double CONVERGENCE_PROXIMITY,\
                                           double TRANS_RATE_FACTOR,\
                                           const double *MPR_arrv,\
                                           const double *DNR_arrv,\
                                           const int *node_type_arrv,
                                           const int *edge_source_arrv,\
                                           const int *edge_target_arrv,\
                                           const int *edge_type_arrv,\
                                           const double *TSH_arrv,\
                                           const int *HCO_arrv,\
                                           const double *FCH_arrv,\
                                           const double *IC_arr, double *exp_arr,\
                                           double *fX_arr,\
                                           double *NOISE_strength,\
                                           double *NOISE_strength_shot,\
                                           double NOISE, double NOISE_SHOT,\
                                           double SCALING)
{
    //CONSTANTS:
    const int ITER_FOR_RELAXATION=sim_config.ITER_FOR_RELAXATION;
    
    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    double test_delta=0.0;
    int iter_count=0;
    int i=0;
    
    double fX_norm;
    double fX_norm_now;
    
    //copy ICs to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=IC_arr[i];
    }
    FILE *fh_states;
    fh_states=fopen("Stochastic_annealing.txt","a");
    
    while(iter_count<ITER_FOR_ODE)
    {
        double NOISE_annealing=pow(SCALING, iter_count)*NOISE;
        double NOISE_SHOT_annealing=pow(SCALING, iter_count)*NOISE_SHOT;
        if(iter_count==ITER_FOR_ODE-1)
        {NOISE_annealing=0; NOISE_SHOT_annealing=0;}
        cal_euler_approximation_stochastic(NUM_NODES,NUM_EDGES,TOTAL_EULER_STEPS,\
                                           EULER_STEP_SIZE,TRANS_RATE_FACTOR,\
                                           MPR_arrv,DNR_arrv,node_type_arrv,\
                                           edge_source_arrv,edge_target_arrv,\
                                           edge_type_arrv,TSH_arrv,HCO_arrv,FCH_arrv,\
                                           curr_exp_arr,next_exp_arr,fX_arr,NOISE_strength,\
                                           NOISE_strength_shot,NOISE_annealing,NOISE_SHOT_annealing);
        
        /*
         test_delta=sum_delta(curr_exp_arr,next_exp_arr,NUM_NODES);
         */
        //check for convergence to a state:
        //if the simulation converges to a state, then break out of the loop:
        /*
         if (test_delta<=CONVERGENCE_PROXIMITY){
         break;
         }
         */
        /*
         fX_norm_now=cal_norm(fX_arr,NUM_NODES);
         //check for convergence to a limit cycle:
         //after ITER_RELAXATION times of iteration, check whether limit
         //cycle situation arises. If it does, break outf of the loop:
         if (iter_count>ITER_FOR_RELAXATION)
         {
         if(fX_norm_now>fX_norm) break;
         }
         fX_norm=fX_norm_now;
         */
        //update curr_exp_arr:
        fprintf(fh_states,"%f\t",NOISE_annealing);
        fprintf(fh_states,"%d\t",MODEL_NO);  //write model index

        for (i=0;i<NUM_NODES;i++){
            exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i];

            fprintf(fh_states,"%10.6f\t",log(exp_arr[i])/log(2));
        }

       fprintf(fh_states,"\n");
       
        iter_count+=1;
    }
    //copy the expression to exp_arr for returning:
    fflush(fh_states);
    fclose(fh_states);
    
    //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
    //printf("\n"); 
    return;
}










/*------------------ stochastic annealing of the network ----------------------*/
void stochastic_anneal_network(char *WORK_DIR, 
                               char *FNAME_STATES, \
                               char *FNAME_LIMITCYCLES, \
                               char *FNAME_SUMMARY, \
                               const int MODEL_NO,\
                               int NUM_NODES,int NUM_EDGES,\
                               int NUM_RANDOM_ICS,int ITER_FOR_ODE,\
                               double EULER_SIM_TIME,\
                               double EULER_STEP_SIZE,\
                               long double CONVERGENCE_PROXIMITY,\
                               double TRANS_RATE_FACTOR,\
                               const double *MPR_arrv,\
                               const double *DNR_arrv,\
                               const int *node_type_arrv,\
                               const int *edge_source_arrv,\
                               const int *edge_target_arrv,\
                               const int *edge_type_arrv,\
                               const double *TSH_arrv,\
                               const int *HCO_arrv,\
                               const double *FCH_arrv,\
                               double *EXP_dict_arrv,\
                               double *NOISE_strength,\
                               double *NOISE_strength_shot,\
                               double NOISE, double NOISE_SHOT, double SCALING)
{

    //printf("Stochastic model annealing simulation!\n");
    //CONSTANTS:
    const int TOTAL_EULER_STEPS=(int) (EULER_SIM_TIME/EULER_STEP_SIZE);
    
    //VARIABLES:
    double *EXP_dict_arr=(double *) EXP_dict_arrv;
    double IC_arr[NUM_NODES]; //saves ICs generated by set_ICs()
    double exp_arr[NUM_NODES]; //saves exps from estimate_stable_expression()
    
    double minEXP_arr[NUM_NODES]; //minimum expression for all nodes
    double maxEXP_arr[NUM_NODES]; //maximum expression for all nodes
    
    double EXP_arr[NUM_RANDOM_ICS][NUM_NODES];
    double fX_arr_norm[NUM_RANDOM_ICS];
    double fX_arr[NUM_NODES];
    double zero_arr[NUM_NODES];
    
    //INDEX VARIABLES:
    int count_pos=0; //keeps track of position in EXP_dict_arr
    int count_iteration; //counts iteration over IC initialization
    int i=0; //general purpose
    
    //read the simulation configations:
    char fname_cfg[30]="racipe.cfg";
    
    FILE *fh_states;
    fh_states=fopen(FNAME_STATES, "a");

    //printf("%d\n",sim_config.ITER_FOR_ODE);
    //printf("%f\n",sim_config.EULER_SIM_TIME);
    //printf("%f\n",sim_config.EULER_SIM_TIME_LAST_ANN_STEP);
    //printf("%f\n",sim_config.STARTING_NOISE_ANN);
    
    //printf("%d\n",sim_config.NUM_RANDOM_ICS);
    read_config(WORK_DIR,fname_cfg);
    //printf("%d\n",sim_config.NUM_RANDOM_ICS);

    //printf("%d\n",sim_config.ITER_FOR_ODE);
    //printf("%f\n",sim_config.EULER_SIM_TIME);
    //printf("%f\n",sim_config.EULER_SIM_TIME_LAST_ANN_STEP);
    //printf("%f\n",sim_config.STARTING_NOISE_ANN);
    //exit(0);

    cal_EXPrange(NUM_NODES,NUM_EDGES,MPR_arrv,DNR_arrv,\
                 edge_source_arrv, edge_target_arrv,\
                 edge_type_arrv,FCH_arrv,minEXP_arr, maxEXP_arr);
        
    //set exp_arr to zeros:
    for (i=0;i<NUM_NODES;++i) {
        exp_arr[i]=0.0;
        fX_arr[i]=0.0;
        zero_arr[i]=0.0;
    }
    
    //SAVE SOLUTIONS: solution_arr
    //FILE *fh_states;
    //fh_states=fopen("Stochastic.txt","a");
    //fprintf(fh_states,"%d\t",MODEL_NO);  //write model index
    
    count_iteration=0;
    while(count_iteration<NUM_RANDOM_ICS)
    {
        set_ICs(NUM_NODES,minEXP_arr,maxEXP_arr,IC_arr);
        estimate_stable_expression_stochastic_anneal_network(\
                                              FNAME_STATES,\
                                              MODEL_NO,\
                                              NUM_NODES,\
                                              NUM_EDGES,ITER_FOR_ODE,\
                                              TOTAL_EULER_STEPS,EULER_STEP_SIZE,\
                                              CONVERGENCE_PROXIMITY,\
                                              TRANS_RATE_FACTOR,\
                                              MPR_arrv,DNR_arrv,node_type_arrv,
                                              edge_source_arrv,\
                                              edge_target_arrv,edge_type_arrv,\
                                              TSH_arrv,HCO_arrv,FCH_arrv,\
                                              IC_arr,exp_arr,fX_arr, NOISE_strength,NOISE_strength_shot,\
                                              NOISE,NOISE_SHOT,SCALING);
        
        write_state_expressions(fh_states, MODEL_NO, NUM_NODES, exp_arr);
        //fX_arr_norm[count_iteration]=cal_norm(fX_arr,NUM_NODES);
        
        //copy the exp_arr to EXP_dict_arr:
        //if any value is found <= zero, assign it min value
        for(i=0;i<NUM_NODES;i++){
            EXP_dict_arr[count_pos+i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            EXP_arr[count_iteration][i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            //fprintf(fh_states,"%10.6f\t",log(EXP_arr[count_iteration][i])/log(2));
            
            
        }
        count_pos+=NUM_NODES;
        count_iteration+=1;
        //fprintf(fh_states,"\n");
    } //end of while(count_iteration<NUM_RANDOM_ICS)
    
    
    fflush(fh_states);
    fclose(fh_states);
    
    return;
}

/*----------------------stochastic annealing-----------------------*/
void estimate_stable_expression_stochastic_anneal_network(\
                                           char *FNAME_STATES, \
                                           const int MODEL_NO,\
                                           int NUM_NODES,int NUM_EDGES,\
                                           int ITER_FOR_ODE,int TOTAL_EULER_STEPS,\
                                           double EULER_STEP_SIZE,\
                                           long double CONVERGENCE_PROXIMITY,\
                                           double TRANS_RATE_FACTOR,\
                                           const double *MPR_arrv,\
                                           const double *DNR_arrv,\
                                           const int *node_type_arrv,
                                           const int *edge_source_arrv,\
                                           const int *edge_target_arrv,\
                                           const int *edge_type_arrv,\
                                           const double *TSH_arrv,\
                                           const int *HCO_arrv,\
                                           const double *FCH_arrv,\
                                           const double *IC_arr, double *exp_arr,\
                                           double *fX_arr,\
                                           double *NOISE_strength,\
                                           double *NOISE_strength_shot,\
                                           double NOISE, double NOISE_SHOT,\
                                           double SCALING)
{
    //CONSTANTS:
    const int ITER_FOR_RELAXATION=sim_config.ITER_FOR_RELAXATION;
    
    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    double test_delta=0.0;
    int iter_count=0;
    int i=0;
    
    double fX_norm;
    double fX_norm_now;
    
    //NOISE parameters for annealing:
    double NOISE_annealing=0.0;
    double NOISE_SHOT_annealing=0.0;

    //copy ICs to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=IC_arr[i];
    }
    //FILE *fh_states;
    //fh_states=fopen("Stochastic_annealing.txt","a");
    //fh_states=fopen(FNAME_STATES, "a");
    
    //while(iter_count<ITER_FOR_ODE)
    while(iter_count<ITER_FOR_ODE-1)
    {
        //change noise level for annealing:
        NOISE_annealing=pow(SCALING, iter_count)*NOISE;
        NOISE_SHOT_annealing=pow(SCALING, iter_count)*NOISE_SHOT;

        cal_euler_approximation_stochastic(NUM_NODES, NUM_EDGES, TOTAL_EULER_STEPS,\
                                           EULER_STEP_SIZE,TRANS_RATE_FACTOR,\
                                           MPR_arrv,DNR_arrv,node_type_arrv,\
                                           edge_source_arrv,edge_target_arrv,\
                                           edge_type_arrv,TSH_arrv,HCO_arrv,FCH_arrv,\
                                           curr_exp_arr,next_exp_arr,fX_arr,\
                                           NOISE_strength,NOISE_strength_shot,\
                                           NOISE_annealing,NOISE_SHOT_annealing);
        
        for (i=0;i<NUM_NODES;i++){
            exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i];
            //fprintf(fh_states,"%10.6f\t",log(exp_arr[i])/log(2));
        }

        iter_count+=1;
    } // end of while loop

    
    //at last iteration, force set NOISE parameters to be zero:
    NOISE_annealing=0; 
    NOISE_SHOT_annealing=0;
   
    //printf("%d\n", TOTAL_EULER_STEPS);
    //set new value for TOTAL_EULER_STEPS:

    //printf("sim_config.EULER_SIM_TIME_LAST_ANN_STEP: %f\n",\
    //      sim_config.EULER_SIM_TIME_LAST_ANN_STEP);

    TOTAL_EULER_STEPS=(int) (sim_config.EULER_SIM_TIME_LAST_ANN_STEP/EULER_STEP_SIZE);

    //printf("%d\n", TOTAL_EULER_STEPS);
    //simulate at Zero noise level: 
    cal_euler_approximation_stochastic(NUM_NODES, NUM_EDGES, TOTAL_EULER_STEPS,\
                                       EULER_STEP_SIZE,TRANS_RATE_FACTOR,\
                                       MPR_arrv,DNR_arrv,node_type_arrv,\
                                       edge_source_arrv,edge_target_arrv,\
                                       edge_type_arrv,TSH_arrv,HCO_arrv,FCH_arrv,\
                                       curr_exp_arr,next_exp_arr,fX_arr,\
                                       NOISE_strength,NOISE_strength_shot,\
                                       NOISE_annealing,NOISE_SHOT_annealing);
    for (i=0;i<NUM_NODES;i++){
         exp_arr[i]=next_exp_arr[i];
         curr_exp_arr[i]=next_exp_arr[i];
         //fprintf(fh_states,"%10.6f\t",log(exp_arr[i])/log(2));
    }

    //write_state_expressions(fh_states, MODEL_NO, NUM_NODES, exp_arr);

    //copy the expression to exp_arr for returning:
        
        //fflush(fh_states);
        //fclose(fh_states);
    
    //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
    //printf("\n"); 
    return;
}

/*-----------------------------------------------------------------------*/
void cal_rk_approximation(int NUM_NODES,int NUM_EDGES,\
                          double RK_TOLERANCE,\
                           double *ptr_TOTAL_TIME,\
                            double *ptr_RK_STEP_SIZE,\
                             double TRANS_RATE_FACTOR,\
                             const double *MPR_arrv,const double *DNR_arrv,\
                             const int *node_type_arrv,\
                             const int *edge_source_arrv,\
                             const int *edge_target_arrv,\
                             const int *edge_type_arrv,\
                             const double *TSH_arrv,\
                             const int *HCO_arrv,\
                             const double *FCH_arrv,\
                             double const *curr_exp_arrv,\
                             double *next_exp_arr,double *fX_arr)
{
        //double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part
    double rk_arr1[NUM_NODES], rk_arr2[NUM_NODES], rk_arr3[NUM_NODES]; 
    double rk_arr4[NUM_NODES], rk_arr5[NUM_NODES], rk_arr6[NUM_NODES];
    double rk_arr7[NUM_NODES];

    double curr_exp_arr[NUM_NODES];
    double hill_impact=0.0;
    double RK_STEP_SIZE=*ptr_RK_STEP_SIZE;
    double tot_time=*ptr_TOTAL_TIME;
    
    //int fX_POINTS=(int) (TOTAL_EULER_STEPS/100);
    int i;
    //EULER_STEP_SIZE*=0.1;
    //    const double FAST_EULER_STEP_SIZE=EULER_STEP_SIZE*TRANS_RATE_FACTOR;
    double fX_norm;
    
    double step_count=0;
    
    //copy start expressions to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=curr_exp_arrv[i];
    }
    //fX_norm=0.0;
    while(step_count<tot_time)
    {
        cal_fX(NUM_NODES,NUM_EDGES,RK_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,rk_arr1);
        for (i=0;i<NUM_NODES;++i) {
            curr_exp_arr[i]=curr_exp_arrv[i]+rk_arr1[i]*RK_STEP_SIZE/5;
        }
        cal_fX(NUM_NODES,NUM_EDGES,RK_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,rk_arr2);
        for (i=0;i<NUM_NODES;++i) {
            curr_exp_arr[i]=curr_exp_arrv[i]+\
                            ((1/4)*rk_arr1[i]+(3/4)*rk_arr2[i])*\
                            RK_STEP_SIZE*3/10;
        }
        cal_fX(NUM_NODES,NUM_EDGES,RK_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,rk_arr3);
        for (i=0;i<NUM_NODES;++i) {
            curr_exp_arr[i]=curr_exp_arrv[i]+\
                            ((11/9)*rk_arr1[i]+(-14/3)*rk_arr2[i]+(40/9)*rk_arr3[i])*\
                            RK_STEP_SIZE*4/5;
        }
        
        cal_fX(NUM_NODES,NUM_EDGES,RK_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,rk_arr4);
        for (i=0;i<NUM_NODES;++i) {
            curr_exp_arr[i]=curr_exp_arrv[i]+\
                            ((4843/1458)*rk_arr1[i]+(-3170/243)*rk_arr2[i]+(8056/729)*rk_arr3[i]+\
                             (-53/162)*rk_arr4[i])*\
                            RK_STEP_SIZE*8/9;
        }
        
        cal_fX(NUM_NODES,NUM_EDGES,RK_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,rk_arr5);
        for (i=0;i<NUM_NODES;++i) {
            curr_exp_arr[i]=curr_exp_arrv[i]+\
                            ((9017/3168)*rk_arr1[i]+(-355/33)*rk_arr2[i]+(46732/5247)*rk_arr3[i]+\
                             (49/176)*rk_arr4[i]+(-5103/18656)*rk_arr5[i])*\
                            RK_STEP_SIZE;
        }
        
        cal_fX(NUM_NODES,NUM_EDGES,RK_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,rk_arr6);
        for (i=0;i<NUM_NODES;++i) {
            curr_exp_arr[i]=curr_exp_arrv[i]+\
                            ((35/384)*rk_arr1[i]+(500/113)*rk_arr3[i]+\
                             (125/192)*rk_arr4[i]+(-2187/6784)*rk_arr5[i]+(11/84)*rk_arr6[i])*\
                            RK_STEP_SIZE;
        }
        cal_fX(NUM_NODES,NUM_EDGES,RK_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,rk_arr7);
        
        //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
        //printf("\n");
        //exit(0);
        //calculate next expressions:
        double max_diff_o4_o5;
        for (i=0;i<NUM_NODES;++i) {
            curr_exp_arr[i]=curr_exp_arrv[i]+\
                            ((5179/57600)*rk_arr1[i]+(7571/16695)*rk_arr3[i]+\
                             (393/640)*rk_arr4[i]+(-92097/339200)*rk_arr5[i]+\
                             (187/2100)*rk_arr6[i]+(1/40)*rk_arr7[i])*\
                            RK_STEP_SIZE;
            
            next_exp_arr[i]=curr_exp_arrv[i]+\
                            ((35/384)*rk_arr1[i]+(500/1113)*rk_arr3[i]+\
                             (125/192)*rk_arr4[i]+(-2187/6784)*rk_arr5[i]+(11/84)*rk_arr6[i])*\
                            RK_STEP_SIZE;
            
            double diff_o4_o5=next_exp_arr[i]-curr_exp_arr[i];
            diff_o4_o5=diff_o4_o5>=0?diff_o4_o5:-diff_o4_o5;
            max_diff_o4_o5=max_diff_o4_o5>diff_o4_o5?max_diff_o4_o5:diff_o4_o5;
        }
        double s_rk=RK_STEP_SIZE*RK_TOLERANCE/(2*(tot_time)*max_diff_o4_o5);
        s_rk=pow(s_rk,0.25);
        
        if(s_rk<1) RK_STEP_SIZE=0.5*RK_STEP_SIZE;
        else {
            step_count+=RK_STEP_SIZE;
            if(s_rk>2) {RK_STEP_SIZE=2*RK_STEP_SIZE;}
            
        }
        //update curr_exp_arr:
        //prevent curr_exp_arr value from getting -ve value:
        //set it to zero, if negative exp value is encountered:
        for (i=0;i<NUM_NODES;++i) {
            //curr_exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i]>=0?next_exp_arr[i]:0;
        }
        //step_count+=1;
        
        /*
         if ((step_count%fX_POINTS)==0)
         {
         fX_norm+=cal_norm(fX_arr,NUM_NODES);
         }
         */
    }
    *ptr_RK_STEP_SIZE=RK_STEP_SIZE;
    *ptr_TOTAL_TIME=step_count;
    //printf("%f\n",fX_norm);
    return;
}

/*-----------------------------Euler Maruyama approximtion-------------------------------------*/
void cal_euler_approximation_stochastic (int NUM_NODES,int NUM_EDGES,\
                             int TOTAL_EULER_STEPS,double EULER_STEP_SIZE,\
                             double TRANS_RATE_FACTOR,\
                             const double *MPR_arrv,\
                            const double *DNR_arrv,\
                             const int *node_type_arrv,\
                             const int *edge_source_arrv,\
                             const int *edge_target_arrv,\
                             const int *edge_type_arrv,\
                             const double *TSH_arrv,\
                             const int *HCO_arrv,\
                             const double *FCH_arrv,\
                             double const *curr_exp_arrv,\
                             double *next_exp_arr,\
                            double *fX_arr,\
                            double *NOISE_strength,\
                            double *NOISE_strength_shot,\
                            double NOISE, double NOISE_SHOT )
{
    
    //double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part
    double curr_exp_arr[NUM_NODES];
    double hill_impact=0.0;
    int fX_POINTS=(int) (TOTAL_EULER_STEPS/100);
    int i;
    double sqrt_EULER_STEP_SIZE=sqrt(EULER_STEP_SIZE);
    //EULER_STEP_SIZE*=0.1;
    //    const double FAST_EULER_STEP_SIZE=EULER_STEP_SIZE*TRANS_RATE_FACTOR;
    double fX_norm;
    //printf("%s\n", "NOISE");
    //printf("%f\n",NOISE);
    
    int step_count=0;
    
    //copy start expressions to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=curr_exp_arrv[i];
    }
    //fX_norm=0.0;
    while(step_count<TOTAL_EULER_STEPS)
    {
        cal_fX(NUM_NODES,NUM_EDGES,EULER_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr,fX_arr);
        //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
        //printf("\n");
        //exit(0);
        //calculate next expressions:
        for (i=0;i<NUM_NODES;++i) {
            next_exp_arr[i]=curr_exp_arr[i]+fX_arr[i]*EULER_STEP_SIZE
            +NOISE*NOISE_strength[i]*sqrt_EULER_STEP_SIZE*pcg32_gasdev()
            +NOISE_SHOT*curr_exp_arr[i]*NOISE_strength_shot[i]*sqrt_EULER_STEP_SIZE*pcg32_gasdev();
            //printf("%f\n",NOISE_strength[i]);
        }
        //update curr_exp_arr:
        //prevent curr_exp_arr value from getting -ve value:
        //set it to zero, if negative exp value is encountered:
        for (i=0;i<NUM_NODES;++i) {
            //curr_exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i]>=0?next_exp_arr[i]:0;
        }
        step_count+=1;
        
        /*
         if ((step_count%fX_POINTS)==0)
         {
         fX_norm+=cal_norm(fX_arr,NUM_NODES);
         }
         */
    }
    //printf("%f\n",fX_norm);
    return;
}


/**********----- stochastic annealing using old method ----- ***********/
/*---------------------------- END -----------------------------------*/


/*--------------------- stochastic annealing using GNW -----------------*/
/*************------------START-------------*****************************/
void stochastic_anneal_network_gnw(char *WORK_DIR, 
                               char *FNAME_STATES, \
                               char *FNAME_LIMITCYCLES, \
                               char *FNAME_SUMMARY, \
                               const int MODEL_NO,\
                               int NUM_NODES,int NUM_EDGES,\
                               int NUM_RANDOM_ICS,int ITER_FOR_ODE,\
                               double EULER_SIM_TIME,\
                               double EULER_STEP_SIZE,\
                               long double CONVERGENCE_PROXIMITY,\
                               double TRANS_RATE_FACTOR,\
                               const double *MPR_arrv,\
                               const double *DNR_arrv,\
                               const int *node_type_arrv,\
                               const int *edge_source_arrv,\
                               const int *edge_target_arrv,\
                               const int *edge_type_arrv,\
                               const double *TSH_arrv,\
                               const int *HCO_arrv,\
                               const double *FCH_arrv,\
                               double *EXP_dict_arrv,\
                               double NOISE, \
                               double SCALING)
{
    //printf("Stochastic model annealing simulation!\n");
    //CONSTANTS:
    const int TOTAL_EULER_STEPS=(int) (EULER_SIM_TIME/EULER_STEP_SIZE);
    
    //VARIABLES:
    double *EXP_dict_arr=(double *) EXP_dict_arrv;
    double IC_arr[NUM_NODES]; //saves ICs generated by set_ICs()
    double exp_arr[NUM_NODES]; //saves exps from estimate_stable_expression()
    
    double minEXP_arr[NUM_NODES]; //minimum expression for all nodes
    double maxEXP_arr[NUM_NODES]; //maximum expression for all nodes
    
    double EXP_arr[NUM_RANDOM_ICS][NUM_NODES];
    double fX_arr_norm[NUM_RANDOM_ICS];
    double fX_arr[NUM_NODES];
    double zero_arr[NUM_NODES];
    
    //INDEX VARIABLES:
    int count_pos=0; //keeps track of position in EXP_dict_arr
    int count_iteration; //counts iteration over IC initialization
    int i=0; //general purpose
    
    //read the simulation configations:
    char fname_cfg[30]="racipe.cfg";
    
    //FILE *fh_states;
    //fh_states=fopen(FNAME_STATES, "a");

    read_config(WORK_DIR,fname_cfg);

    cal_EXPrange(NUM_NODES,NUM_EDGES,MPR_arrv,DNR_arrv,\
                 edge_source_arrv, edge_target_arrv,\
                 edge_type_arrv,FCH_arrv,minEXP_arr, maxEXP_arr);
        
    //set exp_arr to zeros:
    for (i=0;i<NUM_NODES;++i) {
        exp_arr[i]=0.0;
        fX_arr[i]=0.0;
        zero_arr[i]=0.0;
    }
    
    //SAVE SOLUTIONS: solution_arr
    //FILE *fh_states;
    //fh_states=fopen("Stochastic.txt","a");
    //fprintf(fh_states,"%d\t",MODEL_NO);  //write model index
    
    count_iteration=0;
    while(count_iteration<NUM_RANDOM_ICS)
    {
        set_ICs(NUM_NODES,minEXP_arr,maxEXP_arr,IC_arr);

        estimate_stable_expression_stochastic_anneal_network_gnw(\
                                              FNAME_STATES,\
                                              MODEL_NO,\
                                              NUM_NODES,\
                                              NUM_EDGES,ITER_FOR_ODE,\
                                              TOTAL_EULER_STEPS,EULER_STEP_SIZE,\
                                              CONVERGENCE_PROXIMITY,\
                                              TRANS_RATE_FACTOR,\
                                              MPR_arrv,DNR_arrv,node_type_arrv,
                                              edge_source_arrv,\
                                              edge_target_arrv,edge_type_arrv,\
                                              TSH_arrv,HCO_arrv,FCH_arrv,\
                                              IC_arr,exp_arr,fX_arr, \
                                              NOISE, SCALING);
        
        //write_state_expressions(fh_states, MODEL_NO, NUM_NODES, exp_arr);
        //fX_arr_norm[count_iteration]=cal_norm(fX_arr,NUM_NODES);
        
        //copy the exp_arr to EXP_dict_arr:
        //if any value is found <= zero, assign it min value
        for(i=0;i<NUM_NODES;i++){
            EXP_dict_arr[count_pos+i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            EXP_arr[count_iteration][i]=exp_arr[i]>0?exp_arr[i]:minEXP_arr[i];
            //fprintf(fh_states,"%10.6f\t",log(EXP_arr[count_iteration][i])/log(2));
            
            
        }
        count_pos+=NUM_NODES;
        count_iteration+=1;
        //fprintf(fh_states,"\n");
    } //end of while(count_iteration<NUM_RANDOM_ICS)
    
   // fflush(fh_states);
   // fclose(fh_states);
    
    return;
}

void estimate_stable_expression_stochastic_anneal_network_gnw(\
                                           char *FNAME_STATES, \
                                           const int MODEL_NO,\
                                           int NUM_NODES, int NUM_EDGES,\
                                           int ITER_FOR_ODE, int TOTAL_EULER_STEPS,\
                                           double EULER_STEP_SIZE,\
                                           long double CONVERGENCE_PROXIMITY,\
                                           double TRANS_RATE_FACTOR,\
                                           const double *MPR_arrv,\
                                           const double *DNR_arrv,\
                                           const int *node_type_arrv,\
                                           const int *edge_source_arrv,\
                                           const int *edge_target_arrv,\
                                           const int *edge_type_arrv,\
                                           const double *TSH_arrv,\
                                           const int *HCO_arrv,\
                                           const double *FCH_arrv,\
                                           const double *IC_arr,\
                                           double *exp_arr,\
                                           double *fX_arr,\
                                           double NOISE,\
                                           double SCALING)
{
    //CONSTANTS:
    const int ITER_FOR_RELAXATION=sim_config.ITER_FOR_RELAXATION;
    
    double curr_exp_arr[NUM_NODES];
    double next_exp_arr[NUM_NODES];
    double test_delta=0.0;
    int iter_count=0;
    int i=0;
    
    double fX_norm;
    double fX_norm_now;
    
    FILE *fh_states;
    fh_states=fopen(FNAME_STATES, "a");


    //NOISE parameters for annealing:
    double NOISE_annealing=0.0;

    //copy ICs to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=IC_arr[i];
    }
    //FILE *fh_states;
    //fh_states=fopen("Stochastic_annealing.txt","a");
    //fh_states=fopen(FNAME_STATES, "a");
   
    //while(iter_count<ITER_FOR_ODE)
    while(iter_count<ITER_FOR_ODE-1)
    {
        //change noise level for annealing:
        NOISE_annealing=pow(SCALING, iter_count)*NOISE;
       
        cal_euler_approximation_stochastic_gnw(NUM_NODES, NUM_EDGES, TOTAL_EULER_STEPS,\
                                               EULER_STEP_SIZE,TRANS_RATE_FACTOR,\
                                               MPR_arrv,DNR_arrv,node_type_arrv,\
                                               edge_source_arrv,edge_target_arrv,\
                                               edge_type_arrv,TSH_arrv,HCO_arrv,FCH_arrv,\
                                               curr_exp_arr,next_exp_arr,fX_arr,\
                                               NOISE_annealing);
       
        for (i=0;i<NUM_NODES;i++){
            exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i];
            //fprintf(fh_states,"%10.6f\t",log(exp_arr[i])/log(2));
        }
        //write_state_expressions_annealed(fh_states, MODEL_NO, NOISE_annealing,\
        //                                 NUM_NODES, exp_arr);

        iter_count+=1;
    } // end of while loop
    

    // at last iteration, force set NOISE parameters to be zero:
    NOISE_annealing=0; 
   
    TOTAL_EULER_STEPS=(int) (sim_config.EULER_SIM_TIME_LAST_ANN_STEP/EULER_STEP_SIZE);

    // simulate at Zero noise level: 
    
    cal_euler_approximation_stochastic_gnw(NUM_NODES, NUM_EDGES, TOTAL_EULER_STEPS,\
                                       EULER_STEP_SIZE,TRANS_RATE_FACTOR,\
                                       MPR_arrv,DNR_arrv,node_type_arrv,\
                                       edge_source_arrv,edge_target_arrv,\
                                       edge_type_arrv,TSH_arrv,HCO_arrv,FCH_arrv,\
                                       curr_exp_arr,next_exp_arr,fX_arr,\
                                       NOISE_annealing); 
    for (i=0;i<NUM_NODES;i++){
         exp_arr[i]=next_exp_arr[i];
         curr_exp_arr[i]=next_exp_arr[i];
         //fprintf(fh_states,"%10.6f\t",log(exp_arr[i])/log(2));
    }

    //write_state_expressions(fh_states, MODEL_NO, NUM_NODES, exp_arr);
    write_state_expressions_annealed(fh_states, MODEL_NO, NOISE_annealing, \
                                     NUM_NODES, exp_arr);
    //copy the expression to exp_arr for returning:
        //fflush(fh_states);
        //fclose(fh_states);
    fflush(fh_states);
    fclose(fh_states);
 
    //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
    //printf("\n"); 
    return;
}

void cal_euler_approximation_stochastic_gnw(int NUM_NODES,int NUM_EDGES,\
                             int TOTAL_EULER_STEPS,double EULER_STEP_SIZE,\
                             double TRANS_RATE_FACTOR,\
                             const double *MPR_arrv,\
                            const double *DNR_arrv,\
                             const int *node_type_arrv,\
                             const int *edge_source_arrv,\
                             const int *edge_target_arrv,\
                             const int *edge_type_arrv,\
                             const double *TSH_arrv,\
                             const int *HCO_arrv,\
                             const double *FCH_arrv,\
                             double const *curr_exp_arrv,\
                             double *next_exp_arr,\
                            double *fX_arr,\
                            double NOISE)
{
    //double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part
    double curr_exp_arr[NUM_NODES];
    double hill_impact=0.0;
    int fX_POINTS=(int) (TOTAL_EULER_STEPS/100);
    int i;
    double sqrt_EULER_STEP_SIZE=sqrt(EULER_STEP_SIZE);
    //EULER_STEP_SIZE*=0.1;
    //    const double FAST_EULER_STEP_SIZE=EULER_STEP_SIZE*TRANS_RATE_FACTOR;
    double fX_norm;
    //printf("%s\n", "NOISE");
    //printf("%f\n",NOISE);
    
    int step_count=0;
    
    //copy start expressions to a local variable:
    for (i=0;i<NUM_NODES;++i) {
        curr_exp_arr[i]=curr_exp_arrv[i];
    }
    //fX_norm=0.0;
    while(step_count<TOTAL_EULER_STEPS)
    {
        cal_fX_gnw(NUM_NODES,NUM_EDGES,EULER_STEP_SIZE,\
               TRANS_RATE_FACTOR,\
               MPR_arrv,DNR_arrv,\
               node_type_arrv,edge_source_arrv,\
               edge_target_arrv,edge_type_arrv,\
               TSH_arrv,HCO_arrv,FCH_arrv,\
               curr_exp_arr, fX_arr, fX2_arr);
        //for (int j=0;j<NUM_NODES;j++) printf("%f\t",fX_arr[j]);
        //printf("\n");
        //exit(0);
        //calculate next expressions:
        for (i=0;i<NUM_NODES;++i) {
            next_exp_arr[i]=curr_exp_arr[i]+fX_arr[i]*EULER_STEP_SIZE + \
                                            NOISE*(sqrt(fX_arr[i]+fX2_arr[i])*pcg32_gasdev() + \
                                                   sqrt(fX2_arr[i])*pcg32_gasdev())*sqrt_EULER_STEP_SIZE;
        }
        //update curr_exp_arr:
        //prevent curr_exp_arr value from getting -ve value:
        //set it to zero, if negative exp value is encountered:
        for (i=0;i<NUM_NODES;++i) {
            //curr_exp_arr[i]=next_exp_arr[i];
            curr_exp_arr[i]=next_exp_arr[i]>=0?next_exp_arr[i]:0;
        }
        step_count+=1;
    }
    return;
}

void cal_fX_gnw(int NUM_NODES,int NUM_EDGES,\
            double EULER_STEP_SIZE,\
            double TRANS_RATE_FACTOR,\
            const double *MPR_arrv,
            const double *DNR_arrv,\
            const int *node_type_arrv,\
            const int *edge_source_arrv,\
            const int *edge_target_arrv,\
            const int *edge_type_arrv,\
            const double *TSH_arrv,\
            const int *HCO_arrv,\
            const double *FCH_arrv,\
            double const *curr_exp_arr,\
            double *fX_arr,\
            double *fX2_arr)
{
    //double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    //double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part 
    int i;

    //fX=Gx:
    //fX2=kX*X:
    for (i=0;i<NUM_NODES;++i) {
        fX_arr[i]=MPR_arrv[i]; 
        fX2_arr[i]=DNR_arrv[i]*curr_exp_arr[i]; 
    } 
    //update fX with regulations:
    for (i=0;i<NUM_EDGES;i++){
        if (edge_type_arrv[i]==1 || edge_type_arrv[i]==5) {
           fX_arr[edge_target_arrv[i]] *=\
	          eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                             	 TSH_arrv[i],HCO_arrv[i],\
                           		 FCH_arrv[i])/FCH_arrv[i]; 
        }
        else if (edge_type_arrv[i]==2 || edge_type_arrv[i]==6) {
           fX_arr[edge_target_arrv[i]] *=\
                  eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                           		 TSH_arrv[i],HCO_arrv[i],\
 				         FCH_arrv[i]);
        }
        else if (edge_type_arrv[i]==3 || edge_type_arrv[i]==4) {
           fX2_arr[edge_target_arrv[i]] *=\
                  eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                           		 TSH_arrv[i],HCO_arrv[i],\
 				         FCH_arrv[i]); 
        }
    }
    //adjust fX by kX*X:
    for (i=0;i<NUM_NODES;++i) {
        fX_arr[i]-=fX2_arr[i];
        if (node_type_arrv[i]==2){
           fX_arr[i]*=TRANS_RATE_FACTOR;
           fX2_arr[i]*=TRANS_RATE_FACTOR;
        }
    }
    return;
}

/*******----- functions for stochastic annealing using GNW ---- *******/
/*************------------END-------------*****************************/


/*************----- cal_expressions ----- *****************************/
/*************------------START-------------****************************/
void cal_expressions(char *WORK_DIR, \
                    char *FNAME_STATES, \
                    char *FNAME_LIMITCYCLES, \
                    char *FNAME_SUMMARY, \
                    const int MODEL_NO,\
                    int NUM_NODES,int NUM_EDGES,\
                    int NUM_RANDOM_ICS,int ITER_FOR_ODE,\
                    double EULER_SIM_TIME,\
                    double EULER_STEP_SIZE,\
                    long double CONVERGENCE_PROXIMITY,\
                    double TRANS_RATE_FACTOR,\
                    const double *MPR_arrv,\
                    const double *DNR_arrv,\
                    const int *node_type_arrv,\
                    const int *edge_source_arrv,\
                    const int *edge_target_arrv,\
                    const int *edge_type_arrv,\
                    const double *TSH_arrv,\
                    const int *HCO_arrv,\
                    const double *FCH_arrv,\
                    double *EXP_dict_arrv) 
/*------------------------------------------------------------------
It carries out two main operations 
   (1) For each of the initial conditions, it estimates stable state 
       of the network (an n dimensional point) by repeatedly invoking 
       estimate_stable_expression function. 
   (2) After obtaining stable states for all initial conditions, 
       it clusters those stable state solutions by invoking 
       cluster_solutions function.

Definitions of the main data structures: 
   MPR_arrv: maximum production rates for each node
   DNR_arrv: degradation rates for each node  
   node_type_arrv: type for each node 

   edge_source_arrv: source id of each edge. Id of a source is the 
      corresponding index of the arrays MPR_arrv, DNR_arrv, and 
      node_type_arrv
   edge_target_arrv: target id of each edge
   edge_type_arrv: type of each edge 
   TSH_arrv: threshold of each edge
   HCO_arrv: hill coefficient for each edge
   FCH_arrv: fold change for each edge

   EXP_dict_arrv: stable expression level of all the nodes
--------------------------------------------------------------------*/
{
    //CONSTANTS:
    const int TOTAL_EULER_STEPS=(int) (EULER_SIM_TIME/EULER_STEP_SIZE);

    //VARIABLES:
    double *EXP_dict_arr=(double *) EXP_dict_arrv;

    double IC_arr[NUM_NODES]; //saves ICs generated by set_ICs()
    double exp_arr[NUM_NODES]; //saves exps from estimate_stable_expression()

    double minEXP_arr[NUM_NODES]; //minimum expression for all nodes 
    double maxEXP_arr[NUM_NODES]; //maximum expression for all nodes
   
    double EXP_arr[NUM_RANDOM_ICS][NUM_NODES];
    double fX_arr_norm[NUM_RANDOM_ICS];
    double fX_arr[NUM_NODES];
    double zero_arr[NUM_NODES];

    //INDEX VARIABLES:
    int count_pos=0; //keeps track of position in EXP_dict_arr 
    int count_iteration; //counts iteration over IC initialization
    int i=0; //general purpose

    //read the simulation configations:
    char fname_cfg[30]="racipe.cfg";

    double *curr_exp_arr=(double *) EXP_dict_arrv;
   
    for(i=0; i<NUM_NODES; i++) {
       curr_exp_arr[i] = pow(2, curr_exp_arr[i]);
    }

    //DEFINE FILE HEADER:
    FILE *fh_states_exp;
    fh_states_exp=fopen(FNAME_STATES,"a");

    //printf("%s\n","in c");
    //exit(0);

    read_config(WORK_DIR,fname_cfg);

    //printf("%s\n",FNAME_STATES);

    cal_fX2(NUM_NODES, NUM_EDGES, EULER_STEP_SIZE,\
           TRANS_RATE_FACTOR,\
           MPR_arrv, DNR_arrv,\
           node_type_arrv, edge_source_arrv,\
           edge_target_arrv, edge_type_arrv,\
           TSH_arrv, HCO_arrv, FCH_arrv,\
           curr_exp_arr, fX_arr);

    //printf("%s\n","cal_fX");
    //printf("%d\t", MODEL_NO); 

    /*
    for (i=0; i<NUM_NODES; i++){
      printf("%8.6f\t", fX_arr[i]); 
    } */

    //write MODEL_NO:
    fprintf(fh_states_exp,"%d\t%d\t%d\t",MODEL_NO, 1, 1);
    int k;
    //write the expressions to the file:
    //for(int k=0;k<NUM_NODES;k++){
    for(k=0;k<NUM_NODES;k++){
        //fprintf(fh_states_exp,"%10.6f\t", fX_arr[k]);
        fprintf(fh_states_exp,"%10.6f\t", log2(fX_arr[k]));
    }

    //new line to separate the expressions from 
    //those of the next time step:
    fprintf(fh_states_exp,"%s","\n");
    fclose(fh_states_exp);
    return;
}

/*-----------------------------------------------------------------------*/
void cal_fX2(int NUM_NODES,int NUM_EDGES,\
            double EULER_STEP_SIZE,\
            double TRANS_RATE_FACTOR,\
            const double *MPR_arrv,
            const double *DNR_arrv,\
            const int *node_type_arrv,\
            const int *edge_source_arrv,\
            const int *edge_target_arrv,\
            const int *edge_type_arrv,\
            const double *TSH_arrv,\
            const int *HCO_arrv,\
            const double *FCH_arrv,\
            double const *curr_exp_arr,\
            double *fX_arr)
{
    //double fX_arr[NUM_NODES]; //fX_arr saves fX values in each Euler step
    double fX2_arr[NUM_NODES]; //fX2_arr saves fX values for degradation part 
    int i;

    //fX=Gx:
    //fX2=kX*X:
    for (i=0;i<NUM_NODES;++i) {
        fX_arr[i]=MPR_arrv[i]; 
        fX2_arr[i]=DNR_arrv[i]*curr_exp_arr[i]; 
    } 
    //update fX with regulations:
    for (i=0;i<NUM_EDGES;i++){
        if(edge_type_arrv[i]==5)
           continue;
        if (edge_type_arrv[i]==1 || edge_type_arrv[i]==5) {
           fX_arr[edge_target_arrv[i]] *=\
	          eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                             	 TSH_arrv[i],HCO_arrv[i],\
                           		 FCH_arrv[i])/FCH_arrv[i];
        }
        else if (edge_type_arrv[i]==2 || edge_type_arrv[i]==6) {
           fX_arr[edge_target_arrv[i]] *=\
                  eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                           		 TSH_arrv[i],HCO_arrv[i],\
 				         FCH_arrv[i]);
        }
        else if (edge_type_arrv[i]==3 || edge_type_arrv[i]==4) {
           fX2_arr[edge_target_arrv[i]] *=\
                  eval_shiftedHill_fn(curr_exp_arr[edge_source_arrv[i]],\
                           		 TSH_arrv[i],HCO_arrv[i],\
 				         FCH_arrv[i]);
        }
    }
    //adjust fX by kX*X:
    for (i=0;i<NUM_NODES;++i) {
        fX_arr[i]=fX_arr[i]/DNR_arrv[i];
    }
    return;
}

/*************----- cal_expressions ----- *****************************/
/*************------------END-------------*****************************/