/*********************************************************************
 ** Stochastic Ranking Evolution Strategy                           **
 ** example with SRES for threestep model                           **
 **   to fit                                                        **
 **                                                                 **
 ** For ACADEMIC RESEARCH, this is licensed with GPL license        **
 ** For COMMERCIAL ACTIVITIES, please contact the authors           **
 **                                                                 **
 ** Copyright (C) 2005 Xinglai Ji (jix1@ornl.gov)                   **
 **                                                                 **
 ** This program is distributed in the hope that it will be useful, **
 ** but WITHOUT ANY WARRANTY; without even the implied warranty of  **
 ** MERCHANTABILITY of FITNESS FOR A PARTICULAR PURPOSE. See the    **
 ** GNU General Public License for more details.                    **
 **                                                                 **
 ** You should have received a copy of the GNU General Public       **
 ** License along with is program; if not, write to the Free        **
 ** Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, **
 ** MA 02111-1307, USA.                                             **
 **                                                                 **
 ** Author: Xinglai Ji (jix1@ornl.gov)                              **
 ** Date:   Mar 11, 2005;  Mar 14, 2005;                            **
 **                                                                 **
 ** Organization: Oak Ridge National Laboratory                     **
 ** References:                                                     **
 **   1. Thomas P. Runarsson and Xin Yao. 2000. Stochastic Ranking  **
 **      for Constrained Evolutionary Optimization. 4(3):284-294.   **
 **      http://cerium.raunvis.hi.is/~tpr/software/sres/            **
 **   2. Thomas Philip Runarsson and Xin Yao. 2005. Search Biases   **
 **      in Constrained Evolutionary Optimization. IEEE             **
 **      Transactions on Systems, Man and Cybernetics -- Part C:    **
 **      Applications and Reviews. 35(2):233-243.                   **
 **   3. Carmen G. Moles, Pedro Mendes and Julio R. Banga. 2003.    **
 **      Parameter estimation in biochemical pathways: A comparison **
 **      of global optimization methods. Genome Research,           **
 **      13:2467-2474.                                              **
 **   4. Scott D. Cohen, Alan C. Hindmarsh. CVODE User Guide.       **
 **      Numerical Mathematics Group, Lawrence Livermore National   **
 **      Laboratory.                                                **
 *********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*********************************************************************
 ** header files for SRES                                           **
 *********************************************************************/
#include "sharefunc.h"
#include "ESSRSort.h"
#include "ESES.h"

/*********************************************************************
 ** header files and macros for CVODE                               **
 *********************************************************************/
#include "llnltyps.h"
#include "cvode.h"
#include "cvdense.h"
#include "nvector.h"
#include "dense.h"
#define Ith(v,i) N_VIth(v,i-1)
#define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1)
/*********************************************************************
 ** threestep model:                                                **
 ** dG1/dt = V1/(1+pow(cP/Ki1,ni1)+pow(Ka1/cS, na1))-k1*G1          **
 ** dG2/dt = V2/(1+pow(cP/Ki2,ni2)+pow(Ka2/M1, na2))-k2*G2          **
 ** dG3/dt = V3/(1+pow(cP/Ki3,ni3)+pow(Ka3/M2, na3))-k3*G3          **
 ** dE1/dt = V4*G1/(K4+G1) - k4*E1                                  **
 ** dE2/dt = V5*G2/(K5+G2) - k5*E2                                  **
 ** dE3/dt = V6*G3/(K6+G3) - k6*E3                                  **
 ** dM1/dt =  kcat1*E1*(1/Km1)*(cS-M1)/(1+cS/Km1+M1/Km2)  \         **
 **          -kcat2*E2*(1/Km3)*(M1-M2)/(1+M1/Km3+M2/Km4)            **
 ** dM2/dt =  kcat2*E2*(1/Km3)*(M1-M2)/(1+M1/Km3+M2/Km4)  \         **
 **          -kcat3*E3*(1/Km5)*(M2-cP)/(1+M2/Km5+cP/Km6)            **
 ** initialization:                                                 **
 **   G1=0.66667, G2=0.57254, G3=0.41758,                           **
 **   E1=0.4, E2=0.36409, E3=0.29457,                               **
 **   M1=1.419, M2=0.93464,                                         **
 *********************************************************************/
/*********************************************************************
 ** NEQ differential equations                                      ** 
 *********************************************************************/
int NEQ;
/*********************************************************************
 ** the objective value is about order 1                            **
 ** set RTOL = ATOL = 1e-6                                          **
 ** error = ATOL + RTOL * value                                     **
 ** set itol = SS                                                   **
 *********************************************************************/
double RTOL, ATOL;
/*********************************************************************
 ** for(t=T0; t<=Tm; t+=T1)                                         **
 ** use every SEP time points                                       **
 *********************************************************************/
double T0, T1, Tm;
int SEP;
/*********************************************************************
 ** number of experiments                                           **
 ** with different cS and cP which keep fixed                       **
 *********************************************************************/
int NEPR;
int NEPRcomb;
static int nepr;
/*********************************************************************
 ** initializing values, pseudoexperimental data,                   **
 ** function to read experimental data                              **
 ** Y0: 0->G1, 1->G2, 2->G3, 3->E1, 4->E2, 5->E3, 6->M1, 7->M2      **
 *********************************************************************/
double ***DataEpr, ***DataPrd;
double *p,*s, cS, cP;
double V1,Ki1,ni1,Ka1,na1,k1,V2,Ki2,ni2,Ka2,na2,k2,V3,Ki3,ni3,Ka3,na3,k3,  \
       V4,K4,k4,V5,K5,k5,V6,K6,k6,kcat1,Km1,Km2,kcat2,Km3,Km4,kcat3,Km5,Km6;
double G10, G20, G30, E10, E20, E30, M10, M20;

/*********************************************************************
 ** it's so strange that tn=(int)((Tm-T0)/(T1*SEP)) failed          **
 ** so it has to be set tn=(int)((Tm-T0)/(T1*SEP)+0.0001)           **
 *********************************************************************/
int tn;

/*********************************************************************
 ** fitness(double *x, double *f, double *g)                        **
 ** change the function name as you want                            **
 ** x[dim], g[constraint]                                           **
 *********************************************************************/
void fitness(double *, double *, double *);

ESfcnTrsfm *trsfm;

/*********************************************************************
 ** double transform(double x)                                      **
 ** transform x: for coding                                         **
 *********************************************************************/
double transform(double);

/*********************************************************************
 ** read from pseudoexperimental file                               **
 **      to DataEpr[epr][chemical][time]                            **
 *********************************************************************/
void funReadEpr(void);
/*********************************************************************
 ** differential equations                                          **
 *********************************************************************/
static void difeq(integer N, real t, N_Vector y, N_Vector ydot, void *f_data);

static void PrintFinalStats(long int iopt[]);

/*********************************************************************
 ** main function                                                   **
 *********************************************************************/
int main(int argc, char ** argv)
{
  int i;
/*********************************************************************
 ** variables /parameters for SRES                                  **
 *********************************************************************/
  ESParameter *param;
  ESPopulation *population;
  ESStatistics *stats;
  unsigned int seed;
  int es;
  int constraint, dim;
  double *ub, *lb;
  int miu, lambda, gen;
  double gamma, alpha, varphi;
  int retry;
  double pf;

/*********************************************************************
 ** change SRES parameters here as you want                         **
 ** random seed, gamma, alpha, varphi, retry number, pf,            **
 ** if you dont know how to set, keep default settings              **
 **                                                                 **
 ** constraint number, x dimension, miu:parent number,              **
 ** lambda:offspring number, generation number,                     **
 ** up and low bounds on x,                                         **
 *********************************************************************/
  seed = shareDefSeed;
  gamma = esDefGamma;
  alpha = esDefAlpha;
  varphi = esDefVarphi;
  retry = esDefRetry;

  pf = essrDefPf;

  es = esDefESSlash;
  constraint = 1;
  dim = 36;
  miu = 30;
  lambda = 350;
  gen = 100000;
  ub = NULL;
  lb = NULL;
  ub = ShareMallocM1d(dim);
  lb = ShareMallocM1d(dim);

/*********************************************************************
 ** set the up and low bounds on x here                             **
 ** lb[dim] and ub[dim]                                             **
 *********************************************************************/
  for(i=0; i<dim; i++)
  {
    lb[i] = 12;
    ub[i] = 6;
  }
  lb[2]=-1; ub[2]=+1;  
  lb[4]=-1; ub[4]=+1;  
  lb[8]=-1; ub[8]=+1;  
  lb[10]=-1; ub[10]=+1;  
  lb[14]=-1; ub[14]=+1;  
  lb[16]=-1; ub[16]=+1;  

  trsfm = (ESfcnTrsfm *)ShareMallocM1c(dim*sizeof(ESfcnTrsfm));
  for(i=0;i<dim;i++)
    trsfm[i] = transform;

/*********************************************************************
 ** set ODE parameters                                              **
 *********************************************************************/
  NEQ = 8;
  RTOL = 1e-6; 
  ATOL = 1e-6;
  T0 = 0.0;
  T1 = 0.2;
  Tm = 120;
  SEP = 10;
  NEPR = 16;
  NEPRcomb = 4;

  G10 = 0.66667;
  G20 = 0.57254;
  G30 = 0.41758;
  E10 = 0.4;
  E20 = 0.36409;
  E30 = 0.29457;
  M10 = 1.419;
  M20 = 0.93464;
  p = ShareMallocM1d(NEPRcomb);
  s = ShareMallocM1d(NEPRcomb);
  p[0] = 0.05; p[1]=0.13572; p[2]=0.36840;p[3]=1.0;
  s[0]=0.1;  s[1]=0.46416; s[2]=2.1544; s[3]=10;
  tn = (int)((Tm-T0)/(T1*SEP)+0.00001);

  DataEpr = ShareMallocM3d(NEPR,NEQ,tn);
  DataPrd = ShareMallocM3d(NEPR,NEQ,tn);

/********************************************************************
 ** read from pseudoexperiment file                                **
 ********************************************************************/
  funReadEpr();

/********************************************************************
 ** end of user parameter setting                                  **
 ** get started                                                    **
 ********************************************************************/
  ESInitial(seed, &param, trsfm,fitness, es, constraint, dim, ub, lb,  \
            miu, lambda, gen, gamma, alpha, varphi,  \
            retry, &population, &stats);

  while(stats->curgen < param->gen)
    ESStep(population, param, stats, pf);


  ESDeInitial(param, population, stats);

  ShareFreeM1d(ub);
  ub = NULL;
  ShareFreeM1d(lb);
  lb = NULL;
  ShareFreeM1c((char*)trsfm);
  trsfm = NULL;
  ShareFreeM1d(p);
  p = NULL;
  ShareFreeM1d(s);
  s = NULL;
  ShareFreeM3d(DataEpr, NEPR, NEQ);
  ShareFreeM3d(DataPrd, NEPR, NEQ);

  return 0;
}


/*********************************************************************
 ** set the fitness function here                                   **
 ** x[dim], g[constraint]                                           **
 *********************************************************************/
void fitness(double *x, double *f, double *g)
{
  real ropt[OPT_SIZE], reltol, t, tout;
  long int iopt[OPT_SIZE];
  N_Vector y;
  real abstol;
  void *cvode_mem;
  int iout, flag;
  int i,j,k,l;
  int t1;
  double value = 0.0;
  double w,sum;

  V1 = (trsfm[0])(x[0]); Ki1 = (trsfm[1])(x[1]); ni1 = (trsfm[2])(x[2]); 
  Ka1 = (trsfm[3])(x[3]); na1 = (trsfm[4])(x[4]); k1 = (trsfm[5])(x[5]);
  V2 = (trsfm[6])(x[6]); Ki2 = (trsfm[7])(x[7]); ni2 = (trsfm[8])(x[8]); 
  Ka2 = (trsfm[9])(x[9]); na2 = (trsfm[10])(x[10]); k2 = (trsfm[11])(x[11]);
  V3 = (trsfm[12])(x[12]); Ki3 = (trsfm[13])(x[13]); ni3 = (trsfm[14])(x[14]); 
  Ka3 = (trsfm[15])(x[15]); na3 = (trsfm[16])(x[16]); k3 = (trsfm[17])(x[17]);
  V4 = (trsfm[18])(x[18]); K4 = (trsfm[19])(x[19]); k4 = (trsfm[20])(x[20]); 
  V5 = (trsfm[21])(x[21]); K5 = (trsfm[22])(x[22]); k5 = (trsfm[23])(x[23]);
  V6 = (trsfm[24])(x[24]); K6 = (trsfm[25])(x[25]); k6 = (trsfm[26])(x[26]);
  kcat1 = (trsfm[27])(x[27]);Km1 = (trsfm[28])(x[28]); Km2 = (trsfm[29])(x[29]);
  kcat2 = (trsfm[30])(x[30]);Km3 = (trsfm[31])(x[31]); Km4 = (trsfm[32])(x[32]);
  kcat3 = (trsfm[33])(x[33]);Km5 = (trsfm[34])(x[34]); Km6 = (trsfm[35])(x[35]);

  for(i=0; i<NEPRcomb; i++)
    for(j=0; j<NEPRcomb; j++)
    {
      cP = p[i];
      cS = s[j];
      nepr = NEPRcomb*i+j;

      y=N_VNew(NEQ, NULL);
      Ith(y,1) = G10;
      Ith(y,2) = G20;
      Ith(y,3) = G30;
      Ith(y,4) = E10;
      Ith(y,5) = E20;
      Ith(y,6) = E30;
      Ith(y,7) = M10;
      Ith(y,8) = M20;

      reltol = RTOL;
      abstol = ATOL;

      cvode_mem = CVodeMalloc(NEQ, difeq, T0, y, BDF, NEWTON, SS, &reltol,   \
                              &abstol, NULL, NULL, FALSE, iopt, ropt, NULL);
      if(cvode_mem == NULL)
      {
        printf("CVodeMalloc failed.\n"); 
        exit(1);
      }
      CVDense(cvode_mem, NULL, NULL);

      for (iout=1, tout=T1; tout <= Tm; iout++, tout = iout*T1) 
      {
        flag = CVode(cvode_mem, tout, y, &t, NORMAL);
        if (flag != SUCCESS)
        {  
          printf("CVode failed, flag=%d.\n", flag); 
          (*f) = 80000;
          g[0] = 0.0;
          return;
        }
        if(iout%SEP==0)
        {
          k = iout/SEP-1;
          DataPrd[nepr][0][k] = Ith(y,1);
          DataPrd[nepr][1][k] = Ith(y,2);
          DataPrd[nepr][2][k] = Ith(y,3);
          DataPrd[nepr][3][k] = Ith(y,4);
          DataPrd[nepr][4][k] = Ith(y,5);
          DataPrd[nepr][5][k] = Ith(y,6);
          DataPrd[nepr][6][k] = Ith(y,7);
          DataPrd[nepr][7][k] = Ith(y,8);
        }
      }

      N_VFree(y);
      CVodeFree(cvode_mem);
      PrintFinalStats(iopt);
    }

  
  t1 = tn;
  for(i=0;i<NEPRcomb;i++)
    for(j=0;j<NEPRcomb;j++)
    {
      nepr = NEPRcomb*i+j;
      for(k=0; k<t1; k++)
      {
//        w = DataEpr[nepr][0][k];
//        for(l=1; l<NEQ; l++)
//          if(DataEpr[nepr][l][k] > w)
//            w = DataEpr[nepr][l][k];
//        w = (1/w)*(1/w);
        w = 1;
        for(l=0,sum=0.0; l<NEQ; l++)
          sum += (DataEpr[nepr][l][k]-DataPrd[nepr][l][k])  \
                 *(DataEpr[nepr][l][k]-DataPrd[nepr][l][k]);
        value += sum * w;        
      }
    }

  (*f) = value;
  g[0] = 0.0;

  return;
}

/*********************************************************************
 ** read from pseudoexperimental file                               **
 **      to DataEpr[epr][chemical][time]                            **
 *********************************************************************/
void funReadEpr(void)
{
  char file[] = "pseudo.epr";
  char buf[shareDefMaxLine];
  char **s1;
  FILE *fp;
  int i,j,k;
  int t,n;

  if((fp = fopen(file, "r")) == NULL)
  {
    printf("fopen %s failed!\n", file);
    exit(1);
  }

  t = tn;
  for(i=0; i<NEPR; i++)
    for(j=0; j<t; j++)
    {
      fgets(buf, shareDefMaxLine, fp);
      ShareChop(buf);
      s1 = ShareSplitStr(buf, "\t", &n, shareDefNullNo);
      if(n!=NEQ)
      {
        printf("line failed: %s\n", buf);
        exit(1);
      }
      for(k=0;k<n;k++)
        DataEpr[i][k][j] = atof(s1[k]);
      ShareFreeM2c(s1, n);
    }

  fclose(fp);

  return;
}


/*********************************************************************
 ** differential equations                                          **
 *********************************************************************/
static void difeq(integer N, real t, N_Vector y, N_Vector ydot, void *f_data)
{
  real G1, G2, G3, E1, E2, E3, M1, M2;
  real dG1, dG2, dG3, dE1, dE2, dE3, dM1, dM2;

  G1 = Ith(y,1); G2 = Ith(y,2); G3 = Ith(y,3);
  E1 = Ith(y,4); E2 = Ith(y,5); E3 = Ith(y,6);
  M1 = Ith(y,7); M2 = Ith(y,8);

  dG1 = V1/(1+pow(cP/Ki1,ni1)+pow(Ka1/cS, na1))-k1*G1;          
  dG2 = V2/(1+pow(cP/Ki2,ni2)+pow(Ka2/M1, na2))-k2*G2;         
  dG3 = V3/(1+pow(cP/Ki3,ni3)+pow(Ka3/M2, na3))-k3*G3;         
  dE1 = V4*G1/(K4+G1) - k4*E1;         
  dE2 = V5*G2/(K5+G2) - k5*E2;                                 
  dE3 = V6*G3/(K6+G3) - k6*E3;                                
  dM1 =  kcat1*E1*(1/Km1)*(cS-M1)/(1+cS/Km1+M1/Km2)  \
           -kcat2*E2*(1/Km3)*(M1-M2)/(1+M1/Km3+M2/Km4);
  dM2 =  kcat2*E2*(1/Km3)*(M1-M2)/(1+M1/Km3+M2/Km4)  \
           -kcat3*E3*(1/Km5)*(M2-cP)/(1+M2/Km5+cP/Km6);           

  Ith(ydot,1) = dG1; Ith(ydot,2) = dG2; Ith(ydot,3) = dG3;
  Ith(ydot,4) = dE1; Ith(ydot,5) = dE2; Ith(ydot,6) = dE3;
  Ith(ydot,7) = dM1; Ith(ydot,8) = dM2;

  return;
}

static void PrintFinalStats(long int iopt[])
{
/*  printf("\nFinal Statistics.. \n\n");
  printf("nst = %-6ld nfe  = %-6ld nsetups = %-6ld nje = %ld\n",
         iopt[NST], iopt[NFE], iopt[NSETUPS], iopt[DENSE_NJE]);
  printf("nni = %-6ld ncfn = %-6ld netf = %ld\n \n",
         iopt[NNI], iopt[NCFN], iopt[NETF]);*/
}

/*********************************************************************
 ** double transform(double x)                                      **
 ** transform x: for coding                                         **
 *********************************************************************/
double transform(double x)
{
  double y;
                                                                                
  y = pow(10, x);
                                                                                
  return y;
}