#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
#include "libSRES/sharefunc.h"
#include "libSRES/ESSRSort.h"
#include "libSRES/ESES.h"
#include "def.h"
#include "model.h"
#include "nr.h"
Go to the source code of this file.
Functions | |
void | intODE (GlobExp *ex, Glob *globs, int doDerivs, int doPlot, long expNr) |
double | computeRight (Glob *globs, GlobExp *ex) |
void | solvLin (Glob *globs, GlobExp *ex, int computeCovar) |
double | dampIt (Glob *globs, GlobExp *ex, double **P0, double ***S0, double **dP, double ***dS, double *uS) |
void | fitness (double *x, double *f, double *g) |
void | globOpt (GlobExp ex[], Glob *globs) |
Variables | |
Glob * | globsPtr |
GlobExp * | exPtr |
|
Definition at line 39 of file computeRight.cc. 00040 { 00041 long i,j,k; 00042 00043 // calculate discrepancies, residues and constraints 00044 long nPoints=ex->nPoints, nvar=ex->nvar; 00045 long nMeasure=ex->nMeasure, nobs=ex->nobs; 00046 long nP=globs->npar; 00047 00048 // no continuity constraints for first and last mesh point 00049 for (j=1;j<=nvar;++j) { 00050 ex->h[1][j]=0; 00051 ex->h[nPoints][j]=0; 00052 } 00053 for (i=2;i<nPoints;++i) 00054 { 00055 for (j=1;j<=nvar;++j) 00056 ex->h[i][j]=ex->yComp[i][j]-ex->yTry[i][j]; 00057 } 00058 for (i=1;i<=nMeasure;++i) 00059 { 00060 for (j=1;j<=nobs;++j) 00061 ex->residues[i][j]=ex->yPred[i][j]-ex->yMeasure[i][j]; 00062 } 00063 00064 //extra equality and inequality constraints 00065 R2(globs,ex, TRUE); 00066 R3(globs,ex ,TRUE); 00067 00068 #ifdef PRINTDERIVS 00069 for (i=1; i<=nPoints; ++i) { 00070 *dbg << "dyds[" << i << "] = ("; 00071 for (j=1; j<=nvar; ++j) { 00072 for (k=1; k<=nvar; ++k) 00073 *dbg << " " << ex->dyds[i][j][k]; 00074 *dbg << ", "; 00075 } 00076 *dbg << ")\ndydp[" << i << "] = ("; 00077 for (j=1; j<=nvar; ++j) { 00078 for (k=1; k<=nP; ++k) 00079 *dbg << " " << ex->dydp[i][j][k]; 00080 *dbg << ", "; 00081 } 00082 *dbg << ")\n"; 00083 } 00084 for (i=1; i<=nMeasure; ++i) { 00085 *dbg << "dmds[" << i << "] = ("; 00086 for (j=1; j<=nobs; ++j) { 00087 for (k=1; k<=nvar; ++k) 00088 *dbg << " " << ex->dmds[i][j][k]; 00089 *dbg << ", "; 00090 } 00091 *dbg << ")\ndmdp[" << i << "] = ("; 00092 for (j=1; j<=nobs; ++j) { 00093 for (k=1; k<=nP; ++k) 00094 *dbg << " " << ex->dmdp[i][j][k]; 00095 *dbg << ", "; 00096 } 00097 *dbg << ")\n"; 00098 } 00099 for (i=1; i<=ex->me; ++i) { 00100 *dbg << "dR2ds[" << i << "] = ("; 00101 for (j=1; j<=nPoints; ++j) { 00102 for (k=1; k<=nvar; ++k) 00103 *dbg << " " << ex->dR2ds[i][j][k]; 00104 *dbg << ", "; 00105 } 00106 *dbg << ")\n"; 00107 *dbg << "dR2dp[" << i << "] = ("; 00108 for (j=1; j<=nP; ++j) 00109 *dbg << " " << ex->dR2dp[i][j]; 00110 *dbg << ")\n"; 00111 } 00112 for (i=1; i<=ex->mg; ++i) { 00113 *dbg << "dR3ds[" << i << "] = ("; 00114 for (j=1; j<=nPoints; ++j) { 00115 for (k=1; k<=nvar; ++k) 00116 *dbg << " " << ex->dR3ds[i][j][k]; 00117 *dbg << ", "; 00118 } 00119 *dbg << ")\n"; 00120 *dbg << "dR3dp[" << i << "] = ("; 00121 for (j=1; j<=nP; ++j) 00122 *dbg << " " << ex->dR3dp[i][j]; 00123 *dbg << ")\n"; 00124 } 00125 dbg->flush(); 00126 #endif 00127 00128 return(objectiveF(globs,ex)); 00129 }
|
|
Definition at line 155 of file dampIt.cc. References Glob::chisq, FALSE, Norm(), Glob::silent, subIt(), TRUE, and Glob::wquer. 00157 { 00158 int etaOk=FALSE,i=1; 00159 double tau=0.95; 00160 double t; 00161 double tau_min=0.1; 00162 double eta0=1.; 00163 double eta2=1.8; //should be parsed 00164 00165 double normdx=Norm(globs,ex,dS,dP,*uS); 00166 // damping algorithm 00167 00168 if(globs->wquer<0) 00169 { 00170 globs->wquer=1.*eta0/normdx; 00171 } 00172 00173 else 00174 t=eta0/(globs->wquer*normdx); 00175 00176 if (t > tau) 00177 { 00178 t=1; 00179 } 00180 else if(t <= tau_min) 00181 { 00182 t=tau_min; 00183 } 00184 00185 globs->wquer=subIt(globs, ex, t,normdx,P0, S0, dP, dS,*uS); 00186 if(globs->silent!=TRUE) 00187 { 00188 cout << "Start damping:\n"; 00189 cout << "#" << i << " (Predicted) t = " << t << "\t"; 00190 cout << "chisq=" << globs->chisq; 00191 cout << endl; 00192 } 00193 00194 etaOk=((globs->wquer*t*normdx)<eta2); 00195 i++; 00196 00197 while (!etaOk) 00198 { 00199 t=eta0/(normdx*globs->wquer); 00200 00201 if (t > tau) 00202 { 00203 t=1; 00204 etaOk=TRUE; 00205 continue; 00206 } 00207 else if(t < tau_min) 00208 { 00209 t=tau_min; 00210 etaOk=TRUE; 00211 if(globs->silent!=TRUE) 00212 { 00213 cout << "#" << i << " (Corrected) t = " << t << "\t"; 00214 cout << "\nMinimal damping factor tau_min= " << tau_min << " reached.\n"; 00215 } 00216 continue; 00217 } 00218 globs->wquer=subIt(globs, ex, t, normdx,P0, S0, dP, dS,*uS); 00219 if(globs->silent!=TRUE) 00220 { 00221 cout << "#" << i << " (Corrected) t = " << t << "\t"; 00222 cout << "chisq=" << globs->chisq; 00223 cout << endl; 00224 } 00225 etaOk=((globs->wquer*t*normdx)<eta2); 00226 i++; 00227 if(i>=10) 00228 { 00229 cerr << "Too many corrections.\n"; 00230 throw 1; 00231 } 00232 } 00233 return(t); 00234 }
|
|
Definition at line 35 of file globOpt.cc. References exPtr, globsPtr, Glob::nrExp, TRUE, Glob::y0fix, and GlobExp::yTry. 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 }
|
|
Definition at line 89 of file globOpt.cc. References exPtr, globsPtr, lvector(), Glob::nIter, GlobExp::nPoints, and Glob::nrExp. 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, ¶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 //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 }
|
|
Definition at line 313 of file intODE.cc. 00314 { 00315 long i,j,k,ksav,l,m,n,nok,nbad; 00316 long nobs=ex->nobs; 00317 long nPoints=ex->nPoints, nvar=ex->nvar, nMeasure=ex->nMeasure; 00318 double *mesh=ex->mesh, *xMeasure=ex->xMeasure; 00319 long tableLen=maxTableLen(mesh, nPoints, xMeasure, nMeasure)+1; 00320 double temp, h,ymin,ymax; 00321 00322 // allocate memory 00323 double *ys=dvector(1,nvar); 00324 double *xTab=dvector(1,tableLen); 00325 double **yTab=dmatrix(1,tableLen,1,nvar); 00326 double ***dTabds, ***dTabdp; 00327 char name[50]; 00328 ofstream gnuout; 00329 ofstream outout; 00330 00331 00332 //gnuout.open("gnuout.dat"); 00333 00334 // compute derivatives? 00335 // if not, don't write to graphics either 00336 // (no genuine integration, but only estimation of omega) 00337 if (doDerivs) 00338 { 00339 dTabds=d3tensor(1,tableLen,1,nvar,1,nvar); 00340 dTabdp=d3tensor(1,tableLen,1,nvar,1,globs->npar); 00341 } 00342 else 00343 { 00344 dTabds=NULL; 00345 dTabdp=NULL; 00346 } 00347 00348 00349 // first mesh point 00350 if (doDerivs) { 00351 for (i=1; i<=nvar; ++i) 00352 { 00353 for (j=1; j<=nvar; ++j) 00354 { 00355 ex->dyds[1][i][j] = (i==j) ? 1 : 0; 00356 ex->dmds[1][i][j] = (i==j) ? 1 : 0; 00357 } 00358 for (j=1; j<=globs->npar; ++j) 00359 { 00360 ex->dydp[1][i][j] = 0; 00361 ex->dmdp[1][i][j] = 0; 00362 } 00363 } 00364 } 00365 for (j=1; j<=nvar; ++j) 00366 ex->yComp[1][j]=ex->yTry[1][j]; 00367 00368 00369 if (doDerivs && mesh[1]==xMeasure[1]) 00370 { 00371 observation(globs,ex,mesh[1],ex->yTry[1],ex->yPred[1],ex->par,ex->dmds[1],ex->dmdp[1]); 00372 } 00373 00374 // integrate 00375 k=1; 00376 ksav=1; 00377 if (mesh[1]==xMeasure[1]) 00378 { 00379 k=2; 00380 ksav=2; 00381 } 00382 for (i=1; i< nPoints; ++i) 00383 { // given starting values 00384 for (j=1; j<=nvar; ++j) 00385 ys[j]=ex->yTry[i][j]; 00386 if (xMeasure[k]<mesh[i]) 00387 { 00388 cerr << "fatal: intODE: xMeasure table corrupted (mesh: " 00389 << mesh[i] << ", xMeasure: " << xMeasure[k] << ")\n"; 00390 throw 1; 00391 } 00392 // build table of measuring points 00393 n=1; 00394 while (k <= nMeasure && xMeasure[k]<=mesh[i+1]) 00395 xTab[n++]=xMeasure[k++]; 00396 xTab[n]=mesh[i+1]; 00397 #ifdef PRINTMESH 00398 *dbg << "Mesh point " << mesh[i] << ": xTab = ("; 00399 for (j=1; j<=n; ++j) 00400 *dbg << '\t' << xTab[j]; 00401 *dbg << ")\n"; 00402 dbg->flush(); 00403 #endif 00404 //integration 00405 tabulateValues(globs,ex,mesh[i],mesh[i+1],xTab,n,ys,yTab, dTabds, dTabdp,ex->dyds[i+1],ex->dydp[i+1]); 00406 00407 #ifdef PRINTINTEGRATE 00408 *dbg << "Integration returned\t ("; 00409 for (j=1; j<=n; ++j) 00410 { 00411 for (l=1; l<=nvar; ++l) 00412 *dbg << " " << yTab[j][l]; 00413 *dbg << ","; 00414 } 00415 *dbg << ")\n"; 00416 dbg->flush(); 00417 #endif 00418 for (j=1; j<=ex->nvar; ++j) 00419 ex->yComp[i+1][j]=ys[j]; // store value at next mesh point 00420 00421 for (j=ksav; j<k; ++j) // store value at measuring points 00422 { 00423 for (l=1; l<=ex->nobs; ++l) 00424 { 00425 ex->yPred[j][l]=yTab[j-ksav+1][l]; 00426 if (doDerivs) 00427 { 00428 for (m=1; m<=ex->nvar; ++m) 00429 ex->dmds[j][l][m]=dTabds[j-ksav+1][l][m]; 00430 for (m=1; m<=globs->npar; ++m) 00431 ex->dmdp[j][l][m]=dTabdp[j-ksav+1][l][m]; 00432 } 00433 } 00434 } 00435 ksav=k; 00436 } 00437 00438 //GnuPlotting 00439 long mind=2; 00440 if(doDerivs && (globs->noGnu==FALSE) && doPlot==TRUE) 00441 { 00442 for(l=1;l<=ex->nobs;l++) 00443 { 00444 if(globs->gnuindex <= globs->ngnu) 00445 { 00446 mind=2; 00447 ofstream gnuout; 00448 gnuout.open("gnuout.dat"); 00449 ymax=ex->yPred[1][l]; 00450 ymin=ex->yPred[1][l]; 00451 for(j=1;j<=ex->nMeasure;j++) 00452 { 00453 gnuout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl; 00454 if(mesh[mind] <= ex->xMeasure[j]) 00455 { 00456 observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL); 00457 gnuout << mesh[mind] << "\t" << ys[l] << endl << endl; 00458 observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL); 00459 gnuout << mesh[mind] << "\t" << ys[l] << endl; 00460 mind++; 00461 } 00462 //determine y-range 00463 if(ymax <= ex->yPred[j][l]) 00464 ymax=ex->yPred[j][l]; 00465 if(ymin >= ex->yPred[j][l]) 00466 ymin=ex->yPred[j][l]; 00467 if(ymax <= ex->yMeasure[j][l]) 00468 ymax=ex->yMeasure[j][l]; 00469 if(ymin >= ex->yMeasure[j][l]) 00470 ymin=ex->yMeasure[j][l]; 00471 } 00472 gnuout.flush(); 00473 gnuout.close(); 00474 00475 fprintf(globs->gnuFp[globs->gnuindex],"reset\n"); 00476 fprintf(globs->gnuFp[globs->gnuindex],"set yrange[%f:%f]\n",ymin-0.01*ymin,ymax+0.01*ymax); 00477 fprintf(globs->gnuFp[globs->gnuindex],"set xrange[%f:%f]\n",ex->xMeasure[1],ex->xMeasure[ex->nMeasure]); 00478 fflush(globs->gnuFp[globs->gnuindex]); 00479 fprintf(globs->gnuFp[globs->gnuindex],"plot \"gnuout.dat\" title \"\"w l 3,\"gnuout.dat\" u 1:3 title \"\" 1 \n"); 00480 fflush(globs->gnuFp[globs->gnuindex]); 00481 //DELAY 00482 for(long dd=1;dd<=1000000;dd++) 00483 1+1; 00484 } 00485 sprintf(name,"Gnuout_Exp%d_Obs%d.dat",expNr,l); 00486 outout.open(name); 00487 00488 globs->gnuindex++; 00489 mind=2; 00490 for(j=1;j<=ex->nMeasure;j++) 00491 { 00492 00493 outout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl; 00494 if(mesh[mind] <= ex->xMeasure[j]) 00495 { 00496 observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL); 00497 outout << mesh[mind] << "\t" << ys[l] << endl << endl; 00498 observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL); 00499 outout << mesh[mind] << "\t" << ys[l] << endl; 00500 mind++; 00501 } 00502 } 00503 outout.flush(); 00504 outout.close(); 00505 for(long dd=1;dd<=1000000;dd++) 00506 1+1; 00507 } 00508 } 00509 00510 if (doDerivs) { 00511 free_d3tensor(dTabdp,1,tableLen,1,nvar,1,globs->npar); 00512 free_d3tensor(dTabds,1,tableLen,1,nvar,1,nvar); 00513 } 00514 free_dmatrix(yTab,1,tableLen,1,nvar); 00515 free_dvector(xTab,1,tableLen); 00516 free_dvector(ys,1,nvar); 00517 }
|
|
Definition at line 313 of file solvLin.cc. 00314 { 00315 00316 long i,j,k,l,xind=0,yind=0,nL; 00317 long nvarFit=0; // # der zu fittenden Anfangswerte 00318 long nLpara=0; // # lokale Parameter 00319 long nparL=0; // nLpara*nrExp 00320 long npar=globs->npar; 00321 long nvar=ex[1].nvar; 00322 long nrExp=globs->nrExp; 00323 long nparG=0; 00324 long fitDim=0; 00325 long aDim=0,eDim=0,gDim=0; 00326 00327 //big system (condensed) 00328 double **Ma,*Ra; 00329 double **Me,*Re; 00330 double **Mg,*Rg; 00331 00332 //solution 00333 double *dX; 00334 //condense experiments 00335 condense(globs,ex); 00336 00337 //determine dimension of big system 00338 for(i=1;i<=nvar;i++) 00339 { 00340 if(globs->y0fix[i]!=FALSE) 00341 nvarFit++; 00342 } 00343 nvarFit*=nrExp; 00344 for(i=1;i<=npar;i++) 00345 { 00346 if(globs->doP[i]=='L') 00347 nparL++; 00348 else if(globs->doP[i]!=FALSE) 00349 nparG++; 00350 } 00351 nLpara=nparL; 00352 nparL*=nrExp; 00353 fitDim=nvarFit+nparL+nparG; 00354 for(i=1;i<=nrExp;i++) 00355 { 00356 aDim+=ex[i].nMeasure*ex[i].nobs; 00357 eDim+=ex[i].me; 00358 gDim+=ex[i].mg; 00359 } 00360 if(fitDim==0) 00361 { 00362 cerr << "Nothing to fit.\n"; 00363 throw 1; 00364 } 00365 #ifdef PRINTDIMENSION 00366 //print dimension of bigsys; 00367 *dbg << "Dimension of Big-System\n"; 00368 *dbg << "-----------------------\n"; 00369 *dbg << "Least-squares:\n"; 00370 *dbg << "#Rows : " << aDim << " #Columns : " << fitDim << endl; 00371 *dbg << "Equality constraints:\n"; 00372 *dbg << "#Rows : " << eDim << " #Columns : " << fitDim << endl; 00373 *dbg << "Inequality constraints:\n"; 00374 *dbg << "#Rows : " << gDim << " #Columns : " << fitDim; 00375 *dbg << "\n\n"; 00376 dbg->flush(); 00377 #endif 00378 00379 00380 long aDimold=aDim; 00381 long eDimold=eDim; 00382 long fitDimold=fitDim; 00383 00384 // Beginn Local_Parameter_Constraints 00385 if(globs->faktorLexist) 00386 { 00387 // Erweiterung der Dimensionen um nLpara*nExp 00388 aDim=aDim+nparL; 00389 eDim=eDim+nparL; 00390 fitDim=fitDim+nparL; 00391 } 00392 // Ende Local_Parameter_Constraints 00393 00394 //allocate memory 00395 Ma=dmatrix(0,aDim,1,fitDim); 00396 Me=dmatrix(0,eDim,1,fitDim); 00397 Mg=dmatrix(0,gDim,1,fitDim); 00398 Ra=dvector(0,aDim); 00399 Re=dvector(0,eDim); 00400 Rg=dvector(0,gDim); 00401 00402 //initialise 00403 for(i=1;i<=aDim;i++) 00404 { 00405 Ra[i]=0.; 00406 for(j=1;j<=fitDim;j++) 00407 Ma[i][j]=0.; 00408 } 00409 for(i=1;i<=eDim;i++) 00410 { 00411 Re[i]=0.; 00412 for(j=1;j<=fitDim;j++) 00413 Me[i][j]=0.; 00414 } 00415 for(i=1;i<=gDim;i++) 00416 { 00417 Rg[i]=0.; 00418 for(j=1;j<=fitDim;j++) 00419 Mg[i][j]=0.; 00420 } 00421 00422 //least-squares 00423 for(i=1;i<=nrExp;i++) 00424 { 00425 for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++) 00426 Ra[yind+k]=-ex[i].ua[k]; 00427 //initial values 00428 for(j=1;j<=nvar;j++) 00429 { 00430 if(globs->y0fix[j]!=FALSE) 00431 { 00432 xind++; 00433 for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++) 00434 Ma[yind+k][xind]=ex[i].Ea[k][j]; 00435 } 00436 00437 } 00438 //parameters 00439 nL=0; 00440 for(j=1;j<=npar;j++) 00441 { 00442 if(globs->doP[j]=='L') 00443 { 00444 xind++; 00445 for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++) 00446 Ma[yind+k][xind]=ex[i].Pa[k][j]; 00447 } 00448 else if(globs->doP[j]!=FALSE) 00449 { 00450 nL++; 00451 for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++) 00452 Ma[yind+k][nvarFit+nparL+nL]=ex[i].Pa[k][j]; 00453 } 00454 } 00455 yind=yind+(ex[i].nMeasure*ex[i].nobs); 00456 }//end loop over experiments 00457 //-------------> end lest squares 00458 00459 //equality constraints 00460 yind=0; 00461 xind=0; 00462 for(i=1;i<=nrExp;i++) 00463 { 00464 for(k=1;k<=ex[i].me;k++) 00465 Re[yind+k]=-ex[i].ue[k]; 00466 //initial values 00467 for(j=1;j<=nvar;j++) 00468 { 00469 if(globs->y0fix[j]!=FALSE) 00470 { 00471 xind++; 00472 for(k=1;k<=ex[i].me;k++) 00473 Me[yind+k][xind]=ex[i].Ee[k][xind]; 00474 } 00475 } 00476 //parameters 00477 nL=0; 00478 for(j=1;j<=npar;j++) 00479 { 00480 if(globs->doP[j]=='L') 00481 { 00482 xind++; 00483 for(k=1;k<=ex[i].me;k++) 00484 Me[yind+k][xind]=ex[i].Pe[k][j]; 00485 } 00486 else if(globs->doP[j]!=FALSE) 00487 { 00488 nL++; 00489 for(k=1;k<=ex[i].me;k++) 00490 Me[yind+k][nvarFit+nparL+nL]=ex[i].Pe[k][j]; 00491 } 00492 } 00493 yind+=ex[i].me; 00494 }//end loop over experiments 00495 //-------------> end equality constraints 00496 00497 //inequality constraints 00498 yind=0; 00499 xind=0; 00500 for(i=1;i<=nrExp;i++) 00501 { 00502 for(k=1;k<=ex[i].mg;k++) 00503 Rg[yind+k]=-ex[i].ug[k]; 00504 //initial values 00505 for(j=1;j<=nvar;j++) 00506 { 00507 if(globs->y0fix[j]!=FALSE) 00508 { 00509 xind++; 00510 for(k=1;k<=ex[i].mg;k++) 00511 Mg[yind+k][xind]=ex[i].Eg[k][xind]; 00512 } 00513 } 00514 //parameters 00515 nL=0; 00516 for(j=1;j<=npar;j++) 00517 { 00518 if(globs->doP[j]=='L') 00519 { 00520 xind++; 00521 for(k=1;k<=ex[i].mg;k++) 00522 Mg[yind+k][xind]=ex[i].Pg[k][j]; 00523 } 00524 else if(globs->doP[j]!=FALSE) 00525 { 00526 nL++; 00527 for(k=1;k<=ex[i].mg;k++) 00528 Mg[yind+k][nvarFit+nparL+nL]=ex[i].Pg[k][j]; 00529 } 00530 } 00531 yind+=ex[i].mg; 00532 }//end loop over experiments 00533 //-------------> end equality constraints 00534 00535 //Big System Ready to Solve 00536 //------------------------- 00537 // 00538 // Ma*dX ~ Ra 00539 // Me*dX = Re 00540 // Mg*dX >= Rg 00541 // 00542 //------------------------- 00543 00544 //**************************************************************************** 00545 // Beginn Local_Parameter_Constraints 00546 // Erweiterung der Systems auf Local_Parameter_Constraints: 00547 // 00548 // pLmean und pLsig bestimmen 00549 if(globs->faktorLexist) 00550 { 00551 double *pLmean, *pLsig; 00552 int index; 00553 int index_local_parameter; 00554 pLmean=dvector(1,nLpara); 00555 pLsig=dvector(1,nLpara); 00556 00557 // for(i=1;i<=nLpara;i++) 00558 // { 00559 // cerr<<"i="<< i<<" faktorL[i]="<<globs->faktorL[i]<<endl; 00560 // if(globs->faktorL[i]<0) 00561 // { 00562 // //i; 00563 // break; 00564 // } 00565 // } 00566 // cerr<<"i="<< i<<" faktorL[i]="<<globs->faktorL[i]<<endl; 00567 // if(i<=nLpara) 00568 // { 00569 // cerr<<"number of local parameter constraints too small"<<endl; 00570 // throw(1); 00571 // } 00572 for(i=1;i<=nLpara;i++) 00573 { 00574 pLmean[i]=0; 00575 index=0; 00576 index_local_parameter=0; 00577 while(index_local_parameter!=i) 00578 { 00579 index++; 00580 // cerr<<"doP[index]="<<globs->doP[index]<<" index="<<index<<" i="<<i<<endl; 00581 if(globs->doP[index]=='L') 00582 { 00583 index_local_parameter++; 00584 // cerr<<"index="<<index<<" i="<<i<<endl; 00585 } 00586 00587 } 00588 for(j=1;j<=nrExp;j++) 00589 { 00590 pLmean[i]=pLmean[i]+ex[j].par[index]; 00591 // cerr<<"ex[j].par[index]="<<ex[j].par[index]<<endl; 00592 } 00593 pLmean[i]=pLmean[i]/double(nrExp); 00594 // cerr<<"pLmean["<<i<<"]="<<pLmean[i]<<" globs->faktorL["<<i<<"]"<<globs->faktorL[i]<<endl; 00595 } 00596 // Berechnung pLsig 00597 for(i=1;i<=nLpara;i++) 00598 { 00599 pLsig[i]=(globs->faktorL[i]*pLmean[i]-pLmean[i])/(1+globs->faktorL[i]); 00600 if(pLsig[i]<0) 00601 { 00602 cerr<<"ERROR - Negative varianz bei den Lokalen Parameter Constraints..."<<endl; 00603 } 00604 } 00605 // Ra,Re 00606 k=1; 00607 for(i=1;i<=nrExp;i++) 00608 { 00609 for(j=1;j<=nLpara;j++) 00610 { 00611 index=0; 00612 index_local_parameter=0; 00613 while(index_local_parameter!=j) 00614 { 00615 index++; 00616 if(globs->doP[index]=='L') 00617 { 00618 index_local_parameter++; 00619 } 00620 00621 } 00622 Ra[aDimold+k]=(ex[i].par[index]-pLmean[j])/pLsig[j]; 00623 Re[eDimold+k]=0; 00624 double pLsig_real=0; 00625 for(int z=1;z<=nrExp;z++) 00626 { 00627 pLsig_real=pLsig_real+(ex[z].par[index]-pLmean[j])*(ex[z].par[index]-pLmean[j]); 00628 } 00629 pLsig_real=sqrt(pLsig_real/nrExp); 00630 //*dbg<<"ex["<<i<<"].par["<<index<<"]="<<ex[i].par[index]<<" pLmean["<<j<<"]="<<pLmean[j]<<" pLsig["<<j<<"]="<<pLsig[j]<<" pLsig_real="<<pLsig_real<<endl; 00631 // Ma 00632 for(int l=1;l<=nrExp;l++) 00633 { 00634 Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]=1/(pLsig[j]*nrExp); 00635 if(i == l) 00636 { 00637 Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]=Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]-1/pLsig[j]; 00638 } 00639 } 00640 k++; 00641 } 00642 } 00643 // Me 00644 for(i=eDimold+1;i<=eDim;i++) 00645 { 00646 for(j=nvarFit+1;j<=fitDim;j++) 00647 { 00648 if(j==fitDimold+i) 00649 { 00650 Me[i][j]=Me[i][j]+1; 00651 } 00652 if(j==nvarFit+i) 00653 { 00654 Me[i][j]=Me[i][j]-1; 00655 } 00656 } 00657 } 00658 } 00659 // Ende Local_Parameter_Constraints 00660 //**************************************************************************** 00661 00662 //regularising 00663 if(globs->reg==TRUE) 00664 { 00665 double *diag,**V,**MMa; 00666 double cond,maxx; 00667 double reg=0.; 00668 00669 V=dmatrix(1,fitDim,1,fitDim); 00670 MMa=dmatrix(1,aDim,1,fitDim); 00671 diag=dvector(1,fitDim); 00672 00673 for(i=1;i<=aDim;i++) 00674 { 00675 for(j=1;j<=fitDim;j++) 00676 { 00677 MMa[i][j]=Ma[i][j]; 00678 } 00679 } 00680 00681 //SVD 00682 svdcmp(MMa,aDim,fitDim,diag,V); 00683 00684 //determine condition number 00685 *dbg << "\ncondition number : "; 00686 maxx=diag[1]; 00687 cond=diag[1]; 00688 for(i=1;i<=fitDim;i++) 00689 { 00690 if(diag[i]>maxx) 00691 maxx=diag[i]; 00692 if(diag[i]<cond) 00693 cond=diag[i]; 00694 } 00695 00696 //Is the matrix approximately singular? 00697 for(i=1;i<=fitDim;i++) 00698 { 00699 if(fabs(diag[i]/maxx)<=globs->epsilon) 00700 { 00701 diag[i]=0.; 00702 reg=globs->lambda*maxx; 00703 } 00704 } 00705 00706 cond*=1./maxx; 00707 *dbg << cond << endl; 00708 globs->cond=cond; 00709 00710 //writing regularised system back 00711 for(i=1;i<=aDim;i++) 00712 { 00713 for(j=1;j<=fitDim;j++) 00714 { 00715 Ma[i][j]=0.; 00716 for(k=1;k<=fitDim;k++) 00717 { 00718 Ma[i][j]+=MMa[i][k]*diag[k]*V[j][k]; 00719 } 00720 if(i <=fitDim) 00721 { 00722 if(diag[i]==0.) 00723 Ma[i][j]+=reg*V[j][i]; 00724 } 00725 } 00726 } 00727 if(reg!=0.) 00728 { 00729 *dbg << "regularisation active, linear dependency of parameters:\n"; 00730 for(i=1;i<=fitDim;i++) 00731 { 00732 if(diag[i]==0.) 00733 *dbg << "line " << i << ": "; 00734 for(j=1;j<=fitDim;j++) 00735 { 00736 if(diag[i]==0.) 00737 *dbg << V[j][i] << " "; 00738 } 00739 if(diag[i]==0.) 00740 *dbg << endl; 00741 } 00742 } 00743 00744 free_dmatrix(MMa,1,aDim,1,fitDim); 00745 free_dmatrix(V,1,fitDim,1,fitDim); 00746 free_dvector(diag,1,fitDim); 00747 }//end if(globs->reg==TRUE) 00748 00749 #ifdef PRINTBIGSYS 00750 *dbg << "\nLeast-squares :\n"; 00751 for(i=1;i<=aDim;i++) 00752 { 00753 for(j=1;j<=fitDim;j++) 00754 { 00755 *dbg << Ma[i][j] << " "; 00756 } 00757 *dbg << "| " << Ra[i] << endl; 00758 } 00759 *dbg << "\nEquality constraints :\n"; 00760 for(i=1;i<=eDim;i++) 00761 { 00762 for(j=1;j<=fitDim;j++) 00763 { 00764 *dbg << Me[i][j] << " "; 00765 } 00766 *dbg << "| " << Re[i] << endl; 00767 } 00768 *dbg << "\nInquality constraints :\n"; 00769 for(i=1;i<=gDim;i++) 00770 { 00771 for(j=1;j<=fitDim;j++) 00772 { 00773 *dbg << Mg[i][j] << " "; 00774 } 00775 *dbg << "| " << Rg[i] << endl; 00776 } 00777 dbg->flush(); 00778 #endif 00779 //allocate solution vector 00780 dX=(double*)malloc((fitDim+10)*sizeof(double)); 00781 //allocate covariance matrix 00782 globs->fitdim=fitDim; 00783 if(globs->covar==NULL) 00784 globs->covar=dmatrix(1,fitDim,1,fitDim); 00785 00786 if(globs->minimiser==1) 00787 { 00788 // ############### 00789 // # LSEI # 00790 // ############### 00791 00792 //allocation/initialisation for LSEI 00793 double *w; 00794 double *ws; 00795 double *prgopt; 00796 double rnorme,rnorml; 00797 long mdw=aDim+eDim+gDim; 00798 long mode; 00799 long ipstor=mdw+npar+2*fitDim+2; 00800 long wsstor=3*(mdw+npar)+(mdw+npar+2)*(fitDim+8)+fitDim+npar; 00801 long *ip; 00802 00803 w=(double*)malloc(((fitDim+1)*mdw)*sizeof(double)); 00804 ws=(double*)malloc((wsstor+1)*sizeof(double)); 00805 ip=(long*)malloc((ipstor+10)*sizeof(long)); 00806 prgopt=(double*)malloc(30*sizeof(double)); 00807 00808 xind=0; 00809 yind=0; 00810 00811 for(i=1;i<=fitDim;i++) 00812 { 00813 for(j=1;j<=eDim;j++) 00814 { 00815 w[xind*mdw+yind]=Me[j][i]; 00816 yind++; 00817 } 00818 for(j=1;j<=aDim;j++) 00819 { 00820 w[xind*mdw+yind]=Ma[j][i]; 00821 yind++; 00822 } 00823 for(j=1;j<=gDim;j++) 00824 { 00825 w[xind*mdw+yind]=Mg[j][i]; 00826 yind++; 00827 } 00828 yind=0; 00829 xind++; 00830 } 00831 00832 for(j=1;j<=eDim;j++) 00833 { 00834 w[(fitDim)*mdw+yind]=Re[j]; 00835 yind++; 00836 } 00837 for(j=1;j<=aDim;j++) 00838 { 00839 w[(fitDim)*mdw+yind]=Ra[j]; 00840 yind++; 00841 } 00842 for(j=1;j<=gDim;j++) 00843 { 00844 w[(fitDim)*mdw+yind]=Rg[j]; 00845 yind++; 00846 } 00847 00848 ip[0]=wsstor; 00849 ip[1]=ipstor; 00850 00851 if(computeCovar==TRUE) 00852 { 00853 prgopt[0]=4.0; 00854 prgopt[1]=1.0; 00855 prgopt[2]=1.0; 00856 prgopt[3]=7.0; 00857 prgopt[4]=10.0; 00858 prgopt[5]=1.0; 00859 } 00860 else 00861 { 00862 prgopt[0]=4.0; 00863 prgopt[1]=0.0; 00864 prgopt[2]=1.0; 00865 prgopt[3]=7.0; 00866 prgopt[4]=10.0; 00867 prgopt[5]=0.0; 00868 } 00869 prgopt[6]=1.0; 00870 00871 //call LSEI 00872 lsei_(w,mdw,eDim,aDim,gDim,fitDim,prgopt,dX,&rnorme,&rnorml,&mode,ws,ip); 00873 00874 if(mode>=4) 00875 { 00876 cerr << "Minimisation (LSEI) failed.\n"; 00877 throw 1; 00878 } 00879 00880 if(computeCovar==TRUE) 00881 { 00882 for (i=1; i<=fitDim; ++i) 00883 { 00884 for (j=1;j<=fitDim;++j) 00885 globs->covar[i][j] = w[(j-1)*mdw+i-1]; 00886 } 00887 } 00888 00889 free(w); 00890 free(ws); 00891 free(ip); 00892 free(prgopt); 00893 00894 // ####### END OF LSEI ###### 00895 } 00896 else if(globs->minimiser==2) 00897 { 00898 //to solve -----> 00899 //------------------------- 00900 // 00901 // Ma*dX ~ Ra 00902 // Me*dX = Re 00903 // Mg*dX >= Rg 00904 // 00905 //------------------------- 00906 00907 00908 //initialisation of E04NCF 00909 long NCLIN=eDim+gDim; 00910 double *C=(double*)malloc(((eDim+gDim+1)*fitDim)*sizeof(double)); 00911 double *BL=(double*)malloc((fitDim+eDim+gDim+1)*sizeof(double)); 00912 double *BU=(double*)malloc((fitDim+eDim+gDim+1)*sizeof(double)); 00913 double *CVEC=(double*)malloc((fitDim+1)*sizeof(double)); 00914 long *ISTATE=(long*)malloc((fitDim+eDim+gDim+1)*sizeof(long)); 00915 long *KX=(long*)malloc((fitDim+1)*sizeof(long)); 00916 double *X=(double*)malloc((1+fitDim)*sizeof(double)); 00917 double *A=(double*)malloc((1+fitDim*aDim)*sizeof(double)); 00918 double *B=(double*)malloc((1+aDim)*sizeof(double)); 00919 long ITER; 00920 double OBJ; 00921 double *CLAMDA=(double*)malloc((1+fitDim+eDim+gDim)*sizeof(double)); 00922 long *IWORK=(long*)malloc((1+fitDim)*sizeof(long)); 00923 long LWORK=2*fitDim*fitDim+9*fitDim+6*(eDim+gDim); 00924 double *WORK=(double*)malloc((1+LWORK)*sizeof(double)); 00925 long IFAIL=-1; 00926 double **CTMP=dmatrix(1,1+eDim+gDim,1,1+fitDim); 00927 00928 //right-hand side 00929 for(i=1;i<=aDim;i++) 00930 B[i-1]=Ra[i]; 00931 00932 //bounds 00933 for(i=1;i<=fitDim;i++) 00934 { 00935 BL[i-1]=-1e25; 00936 BU[i-1]=1e25; 00937 X[i-1]=0.; //initialisation of X 00938 } 00939 for(i=fitDim+1;i<=fitDim+eDim;i++) 00940 { 00941 BL[i-1]=Re[i-fitDim]; 00942 BU[i-1]=Re[i-fitDim]; 00943 } 00944 for(i=fitDim+eDim+1;i<=fitDim+eDim+gDim;i++) 00945 { 00946 BL[i-1]=Rg[i-fitDim-eDim]; 00947 BU[i-1]=1e25; 00948 } 00949 00950 for(i=1;i<=fitDim;i++) 00951 { 00952 for(j=1;j<=eDim+gDim;j++) 00953 { 00954 if(j<=eDim) 00955 CTMP[j][i]=Me[j][i]; 00956 else 00957 CTMP[j][i]=Mg[j-eDim][i]; 00958 } 00959 } 00960 00961 for(i=1;i<=fitDim+eDim+gDim;i++) 00962 { 00963 ISTATE[i]=0; 00964 } 00965 00966 cmat2fvec(Ma,aDim,fitDim,A); 00967 cmat2fvec(CTMP,eDim+gDim,fitDim,C); 00968 00969 setopte04ncf_(); 00970 e04ncf_(aDim,fitDim,NCLIN,NCLIN,aDim,C,BL,BU,CVEC,ISTATE, 00971 KX,X,A,B,ITER,OBJ,CLAMDA,IWORK,fitDim,WORK,LWORK,IFAIL); 00972 for(i=0;i< fitDim;i++) 00973 { 00974 dX[i]=X[i]; 00975 } 00976 00977 //free memory 00978 free(C); 00979 free(BL); 00980 free(BU); 00981 free(CVEC); 00982 free(ISTATE); 00983 free(KX); 00984 free(X); 00985 free(A); 00986 free(B); 00987 free(CLAMDA); 00988 free(IWORK); 00989 free(WORK); 00990 free_dmatrix(CTMP,1,eDim+gDim,1,fitDim); 00991 } 00992 00993 //decondense experiments 00994 decondense(globs,ex,dX); 00995 00996 #ifdef PRINTSTEP 00997 for(i=1;i<=nrExp;i++) 00998 { 00999 *dbg << "\nExperiment #" << i << ":\n\n"; 01000 for(k=1;k<=ex[i].nPoints-1;k++) 01001 { 01002 *dbg << "dS(msP=" << k << ") = "; 01003 for(j=1;j<=nvar;j++) 01004 *dbg << ex[i].dS[k][j] << " "; 01005 *dbg << endl; 01006 } 01007 *dbg << "dP = "; 01008 for(j=1;j<=npar;j++) 01009 *dbg << ex[i].dP[j] << " "; 01010 } 01011 *dbg << endl; 01012 dbg->flush(); 01013 #endif 01014 01015 01016 //free memory 01017 free_dmatrix(Ma,0,aDim,1,fitDim); 01018 free_dmatrix(Me,0,eDim,1,fitDim); 01019 free_dmatrix(Mg,0,gDim,1,fitDim); 01020 free_dvector(Ra,0,aDim); 01021 free_dvector(Re,0,eDim); 01022 free_dvector(Rg,0,gDim); 01023 01024 dX[0]; 01025 free(dX); 01026 }
|
|
Definition at line 30 of file globOpt.cc. |
|
Definition at line 29 of file globOpt.cc. |