- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/globOpt.cc

Go to the documentation of this file.

    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 //Begin definition of module prototypes
    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 //Pointer to the structures ex and globs
    00029 Glob *globsPtr;
    00030 GlobExp *exPtr;
    00031
    00032
    00033 //fitness function for libSRES
    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               // set up internal data (residues, etc.)
    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   //variables for libSRES
    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   //keeping old values
    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   //determine # of decision variables
    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   //bounds ...for test only
    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, &param, 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   //copying old stuff back
    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 }