/*********************************************************************
 ** Stochastic Ranking Evolution Strategy                           **
 ** example with SRES for hiv model                                 **
 **                                                                 **
 ** 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 28, 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. Petr Kuzmic. 1996. Program DYNAFIT for the analysis of     **
 **      enzyme kinetic data: Application to HIV proteinase.        **
 **      Analytic Biochemistry, 237:260-273.                        **
 **   4. Pedro Mendes and Douglas B. Kell. 1998. Non-linear         **
 **      optimization of biochemical pathways: applications to      **
 **      metabolic engineering and parameter estimation.            **
 **   5. 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)

/*********************************************************************
 ** hiv model                                                       **
 **-----------------------------------------------------------------**
 **    M + M  <=>  E       kmd, kdm                                 **
 **    S + E  <=>  ES      kon, ks                                  **
 **    ES     ==>  P + E   kcat                                     **
 **    P + E  <=>  EP      kon, kp                                  **
 **    I + E  <=>  EI      kon, ki                                  **
 **    EI     ==>  EJ      kde                                      **
 **-----------------------------------------------------------------**
 **    dM/dt  = -2*kmd*M*M + 2*kdm*E                                **
 **    dP/dt  = kcat*ES - kon*P*E + kp*EP                           **
 **    dS/dt  = -kon*S*E + ks*ES                                    **
 **    dI/dt  = -kon*I*E + ki*EI                                    **
 **    dES/dt = kon*S*E - ks*ES - kcat*ES                           **
 **    dEP/dt = kon*P*E - kp*EP                                     **
 **    dE/dt  = kmd*M*M - kdm*E - kon*S*E + ks*ES + kcat*ES         **
 **                     - kon*P*E +kp*EP - kon*I*E + ki*EI          **
 **    dEI/dt = kon*I*E - ki*EI - kde*EI                            **
 **    dEJ/dt = kde*EI                                              **
 **-----------------------------------------------------------------**
 **    one solution:                                                **
 **    k1=kmd=0.1; k2=kdm=0.001; k3=kr=kcat=9.46; k4=kon=100;       **
 **    k5=kp=1117; k6=ks=179.7; k7=ki=0.0831; k8=kde=0.1224;        **
 **    k1=0.1; k2=0.001; k3=9.46; k4=100; k5=1117;                  **
 **    k6=179.7; k7=0.0831; k8=0.1224;                              **
 **    y0 = [0; 0; S; I; 0; 0; E; 0; 0];                            **
 **    signal = 0.024*P + offset                                    **
 *********************************************************************/
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                                       **
 ** t=0:0.4:3600, SEP = 30                                          **
 *********************************************************************/
double T0, T1, Tm;
int SEP;
/*********************************************************************
 ** number of experiments                                           **
 ** with different cI                                               **
 *********************************************************************/
int NEPR;
static int nepr;
/*********************************************************************
 ** initializing values, experimental data,                         **
 ** function to read experimental data                              **
 ** Y0: 0->M, 1->P, 2->S, 3->I, 4->ES, 5->EP, 6->E, 7->EI, 8->EJ    **
 ** signal = slope * P + Offset                                     **
 ** kmd = 0.1, kdm = 0.001, kon = 100, fixed                        **
 *********************************************************************/
double ***DataEpr, ***DataPrd;
double *cI0;
double kmd, kdm, kon;
double ks, kcat, kp, ki, kde;
double *cS0, *cE0, *cOffset;
double slope;
double M0, P0, S0, I0, ES0, EP0, E0, EI0, EJ0, Offset;

/*********************************************************************
 ** 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;
  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 the 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 = 20;
  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<5; i++)
  {
    lb[i] = 1e-6;
    ub[i] = 1e+6;
  }
  for(i=5; i<10; i++)
  {
    lb[i] = 10;
    ub[i] = 40;
  }
  for(i=10; i<15; i++)
  {
    lb[i] = 0.002;
    ub[i] = 0.005;
  }
  for(i=15; i<20; i++)
  {
    lb[i] = -0.2;
    ub[i] = 0.2;
  }

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

/*********************************************************************
 ** set ODE parameters                                              **
 *********************************************************************/
  NEQ = 9;
  RTOL = 1e-6;
  ATOL = 1e-6;
  T0 = 0.0;
  T1 = 0.4;
  Tm = 3600;
  SEP = 30;
  NEPR = 5;

  M0 = 0;
  P0 = 0;
  cS0 = NULL;
  cS0 = ShareMallocM1d(NEPR);
  cI0 = NULL;
  cI0 = ShareMallocM1d(NEPR);
  cI0[0] = 0; cI0[1] = 0.0015; cI0[2] = 0.003;
  cI0[3] = 0.004; cI0[4] = 0.004;
  ES0 = 0;
  EP0 = 0;
  cE0 = NULL;
  cE0 = ShareMallocM1d(NEPR);
  EI0 = 0;
  EJ0 = 0;
  slope = 0.024;
  kmd = 0.1;
  kdm = 0.001;
  kon = 100;
  cOffset = NULL;
  cOffset = ShareMallocM1d(NEPR);

  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);
  printf("\n======\nseed=%u\n======\n", param->seed);

  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(cS0);
  cS0 = NULL;
  ShareFreeM1d(cI0);
  cI0 = NULL;
  ShareFreeM1d(cE0);
  cE0 = NULL;
  ShareFreeM1d(cOffset);
  cOffset = 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,k;
  int t1;
  double value = 0.0;
  double sum;

  ks = (trsfm[0])(x[0]);   kcat = (trsfm[1])(x[1]);   kp = (trsfm[2])(x[2]);
  ki = (trsfm[3])(x[3]);   kde = (trsfm[4])(x[4]);
  cS0[0] = (trsfm[5])(x[5]); cS0[1] = (trsfm[6])(x[6]);
  cS0[2] = (trsfm[7])(x[7]); cS0[3] = (trsfm[8])(x[8]);
  cS0[4] = (trsfm[9])(x[9]);
  cE0[0] = (trsfm[10])(x[10]); cE0[1] = (trsfm[11])(x[11]);
  cE0[2] = (trsfm[12])(x[12]); cE0[3] = (trsfm[13])(x[13]);
  cE0[4] = (trsfm[14])(x[14]);
  cOffset[0] = (trsfm[15])(x[15]); cOffset[1] = (trsfm[16])(x[16]);
  cOffset[2] = (trsfm[17])(x[17]); cOffset[3] = (trsfm[18])(x[18]);
  cOffset[4] = (trsfm[19])(x[19]);
  
  for(i=0; i<NEPR; i++)
  {
    nepr = i;
    I0 = cI0[i];
    S0 = cS0[i];
    E0 = cE0[i];
    Offset = cOffset[i];

    y=N_VNew(NEQ, NULL);
    Ith(y,1) = M0;
    Ith(y,2) = P0;
    Ith(y,3) = S0;
    Ith(y,4) = I0;
    Ith(y,5) = ES0;
    Ith(y,6) = EP0;
    Ith(y,7) = E0;
    Ith(y,8) = EI0;
    Ith(y,9) = EJ0;
    
    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][1][k] = slope*Ith(y,2)+Offset;
      }
    }

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

  t1 = tn;
  for(i=0;i<NEPR;i++)
  {
    nepr = i;
    for(k=0; k<t1; k++)
    {
      sum = (DataEpr[nepr][1][k]-DataPrd[nepr][1][k])  \
            *(DataEpr[nepr][1][k]-DataPrd[nepr][1][k]);
      value += sum;
    }
  }
                                                                                
  (*f) = value;
  g[0] = 0.0;

  return;
}

/*********************************************************************
 ** read from pseudoexperimental file                               **
 **      to DataEpr[epr][chemical][time]                            **
 *********************************************************************/
void funReadEpr(void)
{
  char file[] = "irrevers.txt";
  char buf[shareDefMaxLine];
  char **s1;
  FILE *fp;
  int i,j;
  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!=2)
      {
        printf("line failed: %s\n", buf);
        exit(1);
      }
      DataEpr[i][1][j] = atof(s1[1]);
      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 M, P, S, I, ES, EP, E, EI, EJ;
  real dM, dP, dS, dI, dES, dEP, dE, dEI, dEJ;
                                                                                
  M = Ith(y,1); P = Ith(y,2);  S = Ith(y,3);
  I = Ith(y,4); ES = Ith(y,5); EP = Ith(y,6);
  E = Ith(y,7); EI = Ith(y,8); EJ = Ith(y,9);

  dM  = -2*kmd*M*M + 2*kdm*E;
  dP  = kcat*ES - kon*P*E + kp*EP;
  dS  = -kon*S*E + ks*ES;
  dI  = -kon*I*E + ki*EI;
  dES = kon*S*E - ks*ES - kcat*ES;
  dEP = kon*P*E - kp*EP;
  dE  = kmd*M*M - kdm*E - kon*S*E + ks*ES + kcat*ES  \
                   - kon*P*E +kp*EP - kon*I*E + ki*EI;
  dEI = kon*I*E - ki*EI - kde*EI;
  dEJ = kde*EI;
                                                                                
  Ith(ydot,1) = dM;  Ith(ydot,2) = dP;  Ith(ydot,3) = dS;
  Ith(ydot,4) = dI;  Ith(ydot,5) = dES; Ith(ydot,6) = dEP;
  Ith(ydot,7) = dE;  Ith(ydot,8) = dEI; Ith(ydot,9) = dEJ;
                                                                                
  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 = x;
                                                                                
  return y;
}