00001 #include<iostream>
00002 #include<fstream>
00003 #include<math.h>
00004 #include<string.h>
00005 #include<stdio.h>
00006
00007 #include "../def.h"
00008 #include "../model.h"
00009 #include "../nr.h"
00010
00011 using namespace std;
00012
00013
00014
00015 void setInitialValues(GlobExp *ex, int nP, double *parameters, double *yValues);
00016
00017
00018
00019
00020 #define EPSTOL 1e-8
00021
00022
00023 int initValuesOK(Glob *globs,GlobExp *ex)
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 }
00044
00045
00046 void setupGnuFp(Glob *globs, GlobExp ex[])
00047 {
00048 long gnu_n=20,ncomp,k;
00049 int fullx=700;
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
00062 gnu_n = LMIN (gnu_n, ncomp);
00063 globs->ngnu=gnu_n;
00064
00065
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 }
00094
00095
00096
00097
00098 void initialise(GlobExp ex[],Glob *globs,int simit)
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
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
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
00134 ex[nExp].errP=dvector(1,nP);
00135 ex[nExp].errY0=dvector(1,nvar);
00136
00137 ex[nExp].yTry=dmatrix(1,nPoints,1,nvar);
00138 ex[nExp].yTrySave=dmatrix(1,nPoints,1,nvar);
00139
00140 ex[nExp].yComp=dmatrix(1,nPoints,1,nvar);
00141
00142 ex[nExp].yPred=dmatrix(1,nMeas,1,nobs);
00143
00144 ex[nExp].h=dmatrix(1,nPoints,1,nvar);
00145 ex[nExp].residues=dmatrix(1,nMeas,1,nobs);
00146
00147
00148
00149 ex[nExp].dyds=d3tensor(1,nPoints,1,nvar,1,nvar);
00150
00151 ex[nExp].dydp=d3tensor(1,nPoints,1,nvar,1,nP);
00152
00153 ex[nExp].dmds=d3tensor(1,nMeas,1,nobs,1,nvar);
00154
00155 ex[nExp].dmdp=d3tensor(1,nMeas,1,nobs,1,nP);
00156 if (me>0)
00157 {
00158 ex[nExp].r2=dvector(1,me);
00159
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
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);
00179
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
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
00197 setInitialValues(&ex[nExp],globs->npar,ex[nExp].par,ex[nExp].y0);
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
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
00227
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
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
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
00241 ex[nExp].dS=dmatrix(1,nPoints,1,nvar);
00242 ex[nExp].dP=dvector(1,nP);
00243 }
00244
00245 globs->covar=NULL;
00246 }
00247
00248
00249 if(globs->noGnu==FALSE)
00250 setupGnuFp(globs,ex);
00251
00252 globs->cond=0;
00253 globs->Lambda=1;
00254
00255 char line[1000];
00256
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 }