#include <iostream>
#include <fstream>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "../def.h"
#include "../model.h"
#include "../nr.h"
Go to the source code of this file.
Defines | |
#define | EPSTOL 1e-8 |
Functions | |
void | setInitialValues (GlobExp *ex, int nP, double *parameters, double *yValues) |
int | initValuesOK (Glob *globs, GlobExp *ex) |
void | setupGnuFp (Glob *globs, GlobExp ex[]) |
void | initialise (GlobExp ex[], Glob *globs, int simit) |
|
Definition at line 20 of file initialise.cc. Referenced by initValuesOK(). |
|
Definition at line 98 of file initialise.cc. References d3tensor(), dbg, dmatrix(), GlobExp::dmdp, GlobExp::dmds, Glob::doP, GlobExp::dR2dp, GlobExp::dR2ds, GlobExp::dR3dp, GlobExp::dR3ds, dvector(), GlobExp::dydp, GlobExp::dyds, GlobExp::errP, GlobExp::errY0, FALSE, GlobExp::h, GlobExp::me, GlobExp::mg, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, Glob::nrExp, GlobExp::nvar, GlobExp::par, GlobExp::r2, GlobExp::residues, TRUE, GlobExp::yComp, GlobExp::yPred, GlobExp::yTry, and GlobExp::yTrySave. 00099 { 00100 long i,j,k,nExp; 00101 long nrExp=globs->nrExp; 00102 00103 ifstream inSplines; 00104 00105 if(simit==FALSE) 00106 { 00107 for (nExp=1;nExp<=nrExp;++nExp) 00108 { 00109 00110 //non-local parameters are equal for each experiment 00111 for(i=1;i<=globs->npar;i++) 00112 { 00113 if(globs->doP[i]==TRUE) 00114 ex[nExp].par[i]=ex[1].par[i]; 00115 } 00116 00117 // set ex->me and ex->mg given in model.cc 00118 setNrConstraints(&ex[nExp],globs->npar,ex[nExp].par); 00119 00120 long nPoints=ex[nExp].nPoints; 00121 long nP=globs->npar; 00122 long nvar=ex[nExp].nvar, nobs=ex[nExp].nobs; 00123 long nMeas=ex[nExp].nMeasure; 00124 long me=ex[nExp].me, mg=ex[nExp].mg; 00125 00126 #ifdef PRINTINITVALUES 00127 *dbg << "initialize: exp. #" << nExp << ", nPoints=" << nPoints 00128 << ", nvar=" << nvar << ", nobs=" << nobs << "\nnMeas=" 00129 << nMeas << ", nVal=" << nVal << ", me=" << me 00130 << ", mg=" << mg << ", nP=" << nP << "\n"; 00131 dbg->flush(); 00132 #endif 00133 //std. error -> output 00134 ex[nExp].errP=dvector(1,nP); 00135 ex[nExp].errY0=dvector(1,nvar); 00136 // initial guesses at mesh points 00137 ex[nExp].yTry=dmatrix(1,nPoints,1,nvar); 00138 ex[nExp].yTrySave=dmatrix(1,nPoints,1,nvar); 00139 // computed values at mesh points 00140 ex[nExp].yComp=dmatrix(1,nPoints,1,nvar); 00141 // computed values at measuring points 00142 ex[nExp].yPred=dmatrix(1,nMeas,1,nobs); 00143 00144 ex[nExp].h=dmatrix(1,nPoints,1,nvar); // discrepancies 00145 ex[nExp].residues=dmatrix(1,nMeas,1,nobs); // residues 00146 00147 // derivatives 00148 // d(yComp[i][j]) / d(yTry[i-1][k]) 00149 ex[nExp].dyds=d3tensor(1,nPoints,1,nvar,1,nvar); 00150 // d(yComp[i][j]) / d(p[k]) 00151 ex[nExp].dydp=d3tensor(1,nPoints,1,nvar,1,nP); 00152 // d(yPred[i][j]) / d(yTry[**][k]) 00153 ex[nExp].dmds=d3tensor(1,nMeas,1,nobs,1,nvar); 00154 // d(yPred[i][j]) / d(p[k]) 00155 ex[nExp].dmdp=d3tensor(1,nMeas,1,nobs,1,nP); 00156 if (me>0) 00157 { 00158 ex[nExp].r2=dvector(1,me); // equal. constr. 00159 // d(yR2[i]) / d(yTry[j][k]) 00160 ex[nExp].dR2ds=d3tensor(1,me,1,nPoints,1,nvar); 00161 for(i=1;i<=me;i++) 00162 for(j=1;j<=nPoints;j++) 00163 for(k=1;k<=nvar;k++) 00164 ex[nExp].dR2ds[i][j][k]=0.; 00165 00166 // d(yR2[i]) / d(p[j]) 00167 ex[nExp].dR2dp=dmatrix(1,me,1,nP); 00168 for(i=1;i<=me;i++) 00169 for(j=1;j<=nP;j++) 00170 ex[nExp].dR2dp[i][j]=0.; 00171 } 00172 else 00173 { 00174 ex[nExp].r2=NULL; ex[nExp].dR2ds=NULL; ex[nExp].dR2dp=NULL; 00175 } 00176 if (mg>0) 00177 { 00178 ex[nExp].r3=dvector(1,mg); // ineq. constr. 00179 // d(yR3[i]) / d(yTry[j][k]) 00180 ex[nExp].dR3ds=d3tensor(1,mg,1,nPoints,1,nvar); 00181 for(i=1;i<=mg;i++) 00182 for(j=1;j<=nPoints;j++) 00183 for(k=1;k<=nvar;k++) 00184 ex[nExp].dR3ds[i][j][k]=0.; 00185 // d(yR3[i]) / d(p[j]) 00186 ex[nExp].dR3dp=dmatrix(1,mg,1,nP); 00187 for(i=1;i<=mg;i++) 00188 for(j=1;j<=nP;j++) 00189 ex[nExp].dR3dp[i][j]=0.; 00190 } 00191 else 00192 { 00193 ex[nExp].r3=NULL; ex[nExp].dR3ds=NULL; ex[nExp].dR3dp=NULL; 00194 } 00195 00196 // compute yTry from ex (all ex read from file still present) 00197 setInitialValues(&ex[nExp],globs->npar,ex[nExp].par,ex[nExp].y0); 00198 // check if initial values are compatible with constraints 00199 // CURRENTLY DISABLED (1. Dez. 2004) 00200 // if (!initValuesOK(globs,&ex[nExp])) 00201 // { 00202 // cerr << "Experiment #" << nExp 00203 // << ": initial values are not compatible with constraints\n"; 00204 // cerr << "R2:"; 00205 // for (i=1;i<=ex[nExp].me;++i) 00206 // cerr << "\t" << ex[nExp].r2[i]; 00207 // cerr << "\nR3:"; 00208 // for (i=1;i<=ex[nExp].mg;++i) 00209 // cerr << "\t" << ex[nExp].r3[i]; 00210 // cerr << "\n"; 00211 // exit(1); 00212 // } 00213 00214 #ifdef PRINTEX 00215 *dbg << "Ex used for fitting: \n\n"; 00216 for (i=1; i<=ex[nExp].nMeasure; ++i) 00217 { 00218 *dbg << ex[nExp].xMeasure[i]; 00219 for (j=1; j<=ex[nExp].nobs; ++j) 00220 *dbg << "\t" << ex[nExp].yMeasure[i][j] << "\t" << ex[nExp].sigma[i][j]; 00221 *dbg << '\n'; 00222 } 00223 *dbg << '\n'; 00224 dbg->flush(); 00225 #endif 00226 //Allocate memory for condensation 00227 //least squares 00228 ex[nExp].ua=dvector(1,nMeas*nobs); 00229 ex[nExp].Ea=dmatrix(1,nMeas*nobs,1,nvar); 00230 ex[nExp].Pa=dmatrix(1,nMeas*nobs,1,nP); 00231 //equality constraints 00232 ex[nExp].ue=dvector(1,me); 00233 ex[nExp].Ee=dmatrix(1,me,1,nvar); 00234 ex[nExp].Pe=dmatrix(1,me,1,nP); 00235 //inequality constraints 00236 ex[nExp].ug=dvector(1,mg); 00237 ex[nExp].Eg=dmatrix(1,mg,1,nvar); 00238 ex[nExp].Pg=dmatrix(1,mg,1,nP); 00239 00240 //update steps 00241 ex[nExp].dS=dmatrix(1,nPoints,1,nvar); 00242 ex[nExp].dP=dvector(1,nP); 00243 }// matches: for(nExp=1;nExp<=nrExp;++nExp) loop over experiments 00244 //covar in solvLin allocated 00245 globs->covar=NULL; 00246 } // matches with if(simit==FALSE)... 00247 00248 //Initialise GNUplot 00249 if(globs->noGnu==FALSE) 00250 setupGnuFp(globs,ex); 00251 00252 globs->cond=0; 00253 globs->Lambda=1; 00254 00255 char line[1000]; 00256 //Initialise Splines 00257 if(globs->initSpline==TRUE) 00258 { 00259 for (nExp=1;nExp<=nrExp;++nExp) 00260 { 00261 ex[nExp].splineNodes=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00262 ex[nExp].splineY=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00263 ex[nExp].splineGam=(double **) malloc((size_t)(NSPLINES+1)*sizeof(double*)); 00264 ex[nExp].nNodes=(long unsigned*) malloc((size_t)(NSPLINES+1)*sizeof(long unsigned*)); 00265 00266 for(i=1;i<=NSPLINES;i++) 00267 { 00268 j=0; 00269 while(ex[nExp].splineFile[i][j]!='\0') 00270 { 00271 line[j]=ex[nExp].splineFile[i][j]; 00272 j++; 00273 } 00274 line[j]='\0'; 00275 inSplines.open(line); 00276 if(inSplines==NULL) 00277 { 00278 cerr << "Cannot open " << ex[nExp].splineFile[i] << ".\n"; 00279 exit(1); 00280 } 00281 inSplines >> ex[nExp].nNodes[i]; 00282 ex[nExp].splineNodes[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00283 ex[nExp].splineY[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00284 ex[nExp].splineGam[i]=(double *) malloc((size_t)(ex[nExp].nNodes[i]+1)*sizeof(double)); 00285 for(j=1;j<= ex[nExp].nNodes[i];j++) 00286 { 00287 inSplines >> ex[nExp].splineNodes[i][j]; 00288 inSplines >> ex[nExp].splineY[i][j]; 00289 inSplines >> ex[nExp].splineGam[i][j]; 00290 } 00291 inSplines.close(); 00292 } 00293 } 00294 } 00295 }
|
|
Definition at line 23 of file initialise.cc. References EPSTOL, FALSE, GlobExp::me, GlobExp::mg, GlobExp::r2, GlobExp::r3, and TRUE. 00024 { 00025 int i,good=TRUE; 00026 00027 if (ex->me > 0) 00028 { 00029 R2(globs,ex, FALSE); 00030 for (i=1; i<=ex->me; ++i) 00031 if (good) 00032 good = fabs(ex->r2[i]) < EPSTOL; 00033 } 00034 00035 if (ex->mg>0) 00036 { 00037 R3(globs,ex, FALSE); 00038 for (i=1; i<=ex->mg; ++i) 00039 good = good && (ex->r3[i]>=0); 00040 } 00041 00042 return(good); 00043 }
|
|
Definition at line 8 of file setInitialValues.cc. References GlobExp::mesh, GlobExp::nMeasure, GlobExp::nobs, GlobExp::nPoints, GlobExp::nvar, GlobExp::xMeasure, GlobExp::yMeasure, and GlobExp::yTry. 00009 { 00010 long i,j,k; 00011 double xi, d, w1, w2; 00012 long nobs=ex->nobs, nvar=ex->nvar, nMeasure=ex->nMeasure; 00013 long nPoints=ex->nPoints; 00014 00015 /* set initial values */ 00016 if (yValues) 00017 { 00018 for (j=1; j<=nvar; ++j) 00019 ex->yTry[1][j]=yValues[j]; 00020 i=2; 00021 } 00022 else 00023 i=1; 00024 k=1; 00025 00026 while (i<=nPoints) 00027 { 00028 xi=ex->mesh[i]; 00029 // make educated guess for yTry(xi) 00030 00031 // could be improved (polynomial interpolation/extrapolation etc.) 00032 00033 if (ex->xMeasure[k]>=xi) 00034 { 00035 // before first measurement: take first measured value 00036 for (j=1; j<=nvar; ++j) 00037 if (j<=nobs) 00038 ex->yTry[i][j]=ex->yMeasure[k][j]; 00039 else 00040 ex->yTry[i][j]=1; 00041 } else if (k==nMeasure) 00042 { 00043 // after last measurement: take last measured value 00044 for (j=1; j<=nvar; ++j) 00045 if (j<=nobs) 00046 ex->yTry[i][j]=ex->yMeasure[k][j]; 00047 else 00048 ex->yTry[i][j]=1; 00049 } 00050 else 00051 { 00052 // interpolate linearly 00053 while (k<nMeasure && ex->xMeasure[k]<xi) 00054 ++k; 00055 d=ex->xMeasure[k]-ex->xMeasure[k-1]; 00056 w1=(ex->xMeasure[k]-xi)/d; 00057 w2=(xi-ex->xMeasure[k-1])/d; 00058 for (j=1; j<=nvar; ++j) 00059 if (j<=nobs) 00060 ex->yTry[i][j]=w1*ex->yMeasure[k-1][j] + w2*ex->yMeasure[k][j]; 00061 else 00062 ex->yTry[i][j]=1; 00063 } 00064 ++i; 00065 } 00066 }
|
|
Definition at line 46 of file initialise.cc. References Glob::nrExp. 00047 { 00048 long gnu_n=20,ncomp,k; 00049 int fullx=700; //max. geometry 00050 int fully=500; 00051 int nwi,dx=fullx/2,dy,xpos,ypos; 00052 char outstr[100]; 00053 00054 00055 ncomp=0; 00056 for(k=1;k<=globs->nrExp;k++) 00057 { 00058 ncomp+=ex[k].nobs; 00059 } 00060 00061 // determine nr of windows in y-direction 00062 gnu_n = LMIN (gnu_n, ncomp); 00063 globs->ngnu=gnu_n; 00064 00065 // allocate streams 00066 globs->gnuFp = new FILE *[gnu_n+1]; 00067 00068 if((gnu_n % 2)==1) 00069 nwi=(gnu_n+1)/2; 00070 else 00071 nwi=gnu_n/2; 00072 dy=fully/nwi; 00073 00074 for(k=1;k<=gnu_n;k++) 00075 { 00076 if(k==1) 00077 { 00078 xpos=0; 00079 ypos=0; 00080 } 00081 else if((k % 2)==0) 00082 { 00083 xpos=dx+10; 00084 } 00085 else 00086 { 00087 xpos=0; 00088 ypos+=dy+25; 00089 } 00090 sprintf(outstr,"gnuplot -noraise -geometry %ix%i+%i+%i -title \"%i\"\0",dx,dy,xpos,ypos,k); 00091 globs->gnuFp[k]=popen(outstr,"w"); 00092 } 00093 } // setupGnuFp
|