/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 }

Generated on Mon Jan 29 17:09:12 2007 for Diffit by  doxygen 1.4.6