00001 #include<iostream>
00002 #include<fstream>
00003 #include<stdlib.h>
00004 #include<stdio.h>
00005
00006 #include "libSRES/sharefunc.h"
00007 #include "libSRES/ESSRSort.h"
00008 #include "libSRES/ESES.h"
00009
00010 #include "def.h"
00011 #include "model.h"
00012 #include "nr.h"
00013
00014 using namespace std;
00015
00016
00017
00018 void intODE(GlobExp *ex,Glob *globs,int doDerivs,int doPlot,long expNr);
00019
00020 double computeRight(Glob *globs,GlobExp *ex);
00021
00022 void solvLin(Glob *globs,GlobExp *ex,int computeCovar);
00023
00024 double dampIt(Glob *globs,GlobExp *ex,double **P0, double ***S0,
00025 double **dP, double ***dS,double *uS);
00026
00027
00028
00029 Glob *globsPtr;
00030 GlobExp *exPtr;
00031
00032
00033
00034
00035 void fitness(double *x,double *f,double *g)
00036 {
00037 long nExp,j;
00038 long pos=0;
00039 double sum=0;
00040
00041 for (nExp=1;nExp<=globsPtr->nrExp;++nExp)
00042 {
00043 for(j=1;j<=NEQNS;j++)
00044 {
00045 if(globsPtr->y0fix[j]==TRUE)
00046 {
00047 exPtr[nExp].yTry[1][j]=x[pos];
00048 pos++;
00049 }
00050 }
00051 for(j=1;j<=globsPtr->npar;j++)
00052 {
00053 if(globsPtr->doP[j]=='L')
00054 {
00055 exPtr[nExp].par[j]=x[pos];
00056 pos++;
00057 }
00058 }
00059 }
00060 for(j=1;j<=globsPtr->npar;j++)
00061 {
00062 if(globsPtr->doP[j]==TRUE)
00063 {
00064 for (nExp=1;nExp<=globsPtr->nrExp;++nExp)
00065 exPtr[nExp].par[j]=x[pos];
00066 pos++;
00067 }
00068 }
00069
00070
00071 for (nExp=1;nExp<=globsPtr->nrExp;++nExp)
00072 {
00073 try
00074 {
00075 intODE(&exPtr[nExp],globsPtr,FALSE,FALSE,nExp);
00076
00077 sum+=computeRight(globsPtr,&exPtr[nExp]);
00078 }
00079 catch(int i)
00080 {
00081 sum+=1e10;
00082 }
00083 }
00084 globsPtr->chisq=sum;
00085 (*f)=sum;
00086 }
00087
00088
00089 void globOpt(GlobExp ex[],Glob *globs)
00090 {
00091 globsPtr=globs;
00092 exPtr=ex;
00093
00094
00095 globsPtr->nIter=0;
00096 long i,j;
00097
00098 ESParameter *param;
00099 ESPopulation *population;
00100 ESStatistics *stats;
00101 ESfcnTrsfm *trsfm;
00102 unsigned int seed;
00103 int es;
00104 int dim=0;
00105 double *ub, *lb;
00106 int miu, lambda, gen;
00107 double gamma, alpha, varphi;
00108 int retry;
00109 double pf;
00110
00111
00112 long unsigned *nms=lvector(1,globsPtr->nrExp);
00113 for(i=1;i<=globsPtr->nrExp;i++)
00114 {
00115 nms[i]=exPtr[i].nPoints;
00116 exPtr[i].nPoints=2;
00117 }
00118
00119 seed = shareDefSeed;
00120 gamma = esDefGamma;
00121 alpha = esDefAlpha;
00122 varphi = esDefVarphi;
00123 retry = esDefRetry;
00124 pf = essrDefPf;
00125 es = esDefESSlash;
00126 miu = 30;
00127 lambda = 200;
00128 gen = 1750;
00129 trsfm = NULL;
00130
00131
00132 for(i=1;i<=globsPtr->nrExp;i++)
00133 {
00134 for(j=1;j<=NEQNS;j++)
00135 {
00136 if(globsPtr->y0fix[j]==TRUE)
00137 dim++;
00138 }
00139 for(j=1;j<=globsPtr->npar;j++)
00140 {
00141 if(globsPtr->doP[j]=='L')
00142 dim++;
00143 }
00144 }
00145 for(j=1;j<=globsPtr->npar;j++)
00146 {
00147 if(globsPtr->doP[j]==TRUE)
00148 dim++;
00149 }
00150
00151
00152
00153
00154 lb=dvector(0,dim);
00155 ub=dvector(0,dim);
00156
00157 for(i=0;i<dim;i++)
00158 {
00159 lb[i]=0;
00160 ub[i]=2;
00161 }
00162
00163 ESInitial(seed, ¶m, trsfm, fitness, es,0, dim, ub, lb, \
00164 miu, lambda, gen, gamma, alpha, varphi, \
00165 retry, &population, &stats);
00166
00167 for(int i=1;i<=globsPtr->maxit;i++)
00168 {
00169 globsPtr->nIter++;
00170 ESStep(population, param, stats, pf);
00171 if(stats->curgen >= param->gen)
00172 break;
00173 }
00174
00175
00176 for(i=1;i<=globsPtr->nrExp;i++)
00177 {
00178 exPtr[i].nPoints=nms[i];
00179 }
00180
00181 cout << endl;
00182
00183 free_dvector(lb,0,dim);
00184 free_dvector(ub,0,dim);
00185 free_lvector(nms,1,globsPtr->nrExp);
00186 }