00001 #include<iostream>
00002 #include<fstream>
00003 #include<math.h>
00004 #include<stdlib.h>
00005
00006
00007 #include "../def.h"
00008 #include "../model.h"
00009 #include "../nr.h"
00010
00011 using namespace std;
00012
00013
00014 #define PRINTBIGSYS
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 extern "C" int lsei_(double *w, long &mdw, long &me, long &ma, long &mg,
00025 long &n, double *prgopt, double *x, double *rnorme,
00026 double *rnorml,
00027 long *mode, double *ws, long *ip);
00028
00029
00030 extern "C" int e04ncf_(long &M,long &N,long &NCLIN,long &LDC,long &LDA,
00031 double *C,double *BL,double *BU,double *CVEC,
00032 long *ISTATE,long *KX,double *X,double *A,double *B,
00033 long &ITER,double &OBJ,double *CLAMBDA,long *IWORK,long &LIWORK,
00034 double *WORK,long &LWORK,long &IFAIL);
00035
00036 extern "C" int setopte04ncf_();
00037
00038
00039
00040 void condense(Glob *globs,GlobExp *ex)
00042 {
00043 long nExp,msP,i,j,k,l;
00044 long nMeas,nvar,nP;
00045 long me,mg,nobs,nms;
00046
00047
00048 double elastic=globs->elastic;
00049 double **Ea,**Ee,**Eg;
00050
00051 nP=globs->npar;
00052 for(nExp=1;nExp<=globs->nrExp;++nExp)
00053 {
00054 nMeas=ex[nExp].nMeasure;
00055 nvar=ex[nExp].nvar;
00056 me=ex[nExp].me;
00057 mg=ex[nExp].mg;
00058 nobs=ex[nExp].nobs;
00059 nms=ex[nExp].nPoints;
00060 nms--;
00061 Ea=dmatrix(1,nMeas*nobs,1,nvar);
00062 Ee=dmatrix(1,me,1,nvar);
00063 Eg=dmatrix(1,mg,1,nvar);
00064
00065
00066
00067
00068 for(i=1;i<=nMeas;i++)
00069 {
00070 for(j=1;j<=nobs;j++)
00071 {
00072
00073 ex[nExp].ua[(j-1)*nMeas+i]=ex[nExp].residues[i][j]/ex[nExp].sigma[i][j];
00074
00075 for(k=1;k<=nvar;k++)
00076 {
00077
00078
00079 if(ex[nExp].xMeasure[i] >= ex[nExp].mesh[nms] &&
00080 ex[nExp].xMeasure[i] <= ex[nExp].mesh[nms+1])
00081 {
00082 ex[nExp].Ea[(j-1)*nMeas+i][k]=ex[nExp].dmds[i][j][k]/ex[nExp].sigma[i][j];
00083 }
00084 else
00085 {
00086 ex[nExp].Ea[(j-1)*nMeas+i][k]=0.0;
00087 }
00088 Ea[(j-1)*nMeas+i][k]=ex[nExp].Ea[(j-1)*nMeas+i][k];
00089 }
00090
00091 for(k=1;k<=nP;k++)
00092 {
00093 ex[nExp].Pa[(j-1)*nMeas+i][k]=ex[nExp].dmdp[i][j][k]/ex[nExp].sigma[i][j];
00094 }
00095 }
00096 }
00097
00098
00099 for(i=1;i<=me;i++)
00100 {
00101
00102 ex[nExp].ue[i]=ex[nExp].r2[i];
00103
00104 for(j=1;j<=nvar;j++)
00105 {
00106 ex[nExp].Ee[i][j]=ex[nExp].dR2ds[i][nms][j];
00107 Ee[i][j]=ex[nExp].Ee[i][j];
00108 }
00109
00110 for(j=1;j<=nP;j++)
00111 ex[nExp].Pe[i][j]=ex[nExp].dR2dp[i][j];
00112 }
00113
00114
00115 for(i=1;i<=mg;i++)
00116 {
00117
00118 ex[nExp].ug[i]=ex[nExp].r3[i];
00119
00120 for(j=1;j<=nvar;j++)
00121 {
00122 ex[nExp].Eg[i][j]=ex[nExp].dR3ds[i][nms][j];
00123 Eg[i][j]=ex[nExp].Eg[i][j];
00124 }
00125
00126 for(j=1;j<=nP;j++)
00127 ex[nExp].Pg[i][j]=ex[nExp].dR3dp[i][j];
00128 }
00129
00130
00131
00132
00133
00134
00135
00136 for(msP=nms;msP>=2;msP--)
00137 {
00138
00139
00140 for(i=1;i<=nMeas;i++)
00141 {
00142 for(j=1;j<=nobs;j++)
00143 {
00144
00145 for(k=1;k<=nvar;k++)
00146 ex[nExp].ua[(j-1)*nMeas+i]+=ex[nExp].Ea[(j-1)*nMeas+i][k]*ex[nExp].h[msP][k]*elastic;
00147
00148 for(k=1;k<=nvar;k++)
00149 {
00150
00151
00152 if(ex[nExp].xMeasure[i] >= ex[nExp].mesh[msP-1] &&
00153 ex[nExp].xMeasure[i] <= ex[nExp].mesh[msP])
00154 {
00155 ex[nExp].Ea[(j-1)*nMeas+i][k]=ex[nExp].dmds[i][j][k]/ex[nExp].sigma[i][j];
00156 }
00157 else
00158 {
00159 ex[nExp].Ea[(j-1)*nMeas+i][k]=0.0;
00160 }
00161 for(l=1;l<=nvar;l++)
00162 ex[nExp].Ea[(j-1)*nMeas+i][k]+=Ea[(j-1)*nMeas+i][l]*ex[nExp].dyds[msP][l][k];
00163 }
00164
00165 for(k=1;k<=nP;k++)
00166 {
00167 for(l=1;l<=nvar;l++)
00168 {
00169 ex[nExp].Pa[(j-1)*nMeas+i][k]+=Ea[(j-1)*nMeas+i][l]*ex[nExp].dydp[msP][l][k];
00170 }
00171 }
00172 }
00173 }
00174
00175
00176 for(i=1;i<=me;i++)
00177 {
00178
00179 for(k=1;k<=nvar;k++)
00180 ex[nExp].ue[i]+=Ee[i][k]*ex[nExp].h[msP][k]*elastic;
00181
00182 for(j=1;j<=nvar;j++)
00183 {
00184 ex[nExp].Ee[i][j]=ex[nExp].dR2ds[i][msP-1][j];
00185 for(l=1;l<=nvar;l++)
00186 ex[nExp].Ee[i][j]+=Ee[i][l]*ex[nExp].dyds[msP][l][j];
00187 }
00188
00189 for(j=1;j<=nP;j++)
00190 {
00191 for(l=1;l<=nvar;l++)
00192 ex[nExp].Pe[i][j]+=Ee[i][l]*ex[nExp].dydp[msP][l][j];
00193 }
00194 }
00195
00196
00197 for(i=1;i<=mg;i++)
00198 {
00199
00200 for(k=1;k<=nvar;k++)
00201 ex[nExp].ug[i]+=Eg[i][k]*ex[nExp].h[msP][k]*elastic;
00202
00203 for(j=1;j<=nvar;j++)
00204 {
00205 ex[nExp].Eg[i][j]=ex[nExp].dR3ds[i][msP-1][j];
00206 for(l=1;l<=nvar;l++)
00207 ex[nExp].Eg[i][j]+=Eg[i][l]*ex[nExp].dyds[msP][l][j];
00208 }
00209
00210 for(j=1;j<=nP;j++)
00211 {
00212 for(l=1;l<=nvar;l++)
00213 ex[nExp].Pg[i][j]+=Eg[i][l]*ex[nExp].dydp[msP][l][j];
00214 }
00215 }
00216
00217
00218 for(i=1;i<=nvar;i++)
00219 {
00220 for(j=1;j<=nMeas*nobs;j++)
00221 Ea[j][i]=ex[nExp].Ea[j][i];
00222 for(j=1;j<=me;j++)
00223 Ee[j][i]=ex[nExp].Ee[j][i];
00224 for(j=1;j<=mg;j++)
00225 Eg[j][i]=ex[nExp].Eg[j][i];
00226 }
00227
00228 }
00229
00230
00231 free_dmatrix(Ea,1,nMeas*nobs,1,nvar);
00232 free_dmatrix(Ee,1,me,1,nvar);
00233 free_dmatrix(Eg,1,mg,1,nvar);
00234 }
00235 }
00236
00237
00238 void decondense(Glob *globs,GlobExp *ex,double *dX)
00240 {
00241 long nExp,msP,i,j,k,l;
00242 long nvar,npar;
00243 long nms,ind=0;
00244
00245
00246 double elastic=globs->elastic;
00247
00248 npar=globs->npar;
00249 nvar=ex[1].nvar;
00250
00251
00252 for(nExp=1;nExp<=globs->nrExp;++nExp)
00253 {
00254 for(i=1;i<=ex[nExp].nPoints;i++)
00255 for(j=1;j<=nvar;j++)
00256 ex[nExp].dS[i][j]=0.;
00257 for(i=1;i<=npar;i++)
00258 ex[nExp].dP[i]=0.;
00259
00260 for(i=1;i<=nvar;i++)
00261 {
00262 if(globs->y0fix[i]!=FALSE)
00263 {
00264 ex[nExp].dS[1][i]=dX[ind];
00265 ind++;
00266 }
00267 }
00268 for(i=1;i<=npar;i++)
00269 {
00270 if(globs->doP[i]=='L')
00271 {
00272 ex[nExp].dP[i]=dX[ind];
00273 ind++;
00274 }
00275 }
00276 }
00277
00278 for(i=1;i<=npar;i++)
00279 {
00280 if(globs->doP[i]==TRUE)
00281 {
00282 for(nExp=1;nExp<=globs->nrExp;++nExp)
00283 {
00284 ex[nExp].dP[i]=dX[ind];
00285 }
00286 ind++;
00287 }
00288 }
00289
00290
00291
00292
00293 for(nExp=1;nExp<=globs->nrExp;++nExp)
00294 {
00295 nms=ex[nExp].nPoints;
00296 nms--;
00297
00298 for(msP=2;msP<=nms;msP++)
00299 {
00300 for(i=1;i<=nvar;i++)
00301 {
00302 ex[nExp].dS[msP][i]=elastic*ex[nExp].h[msP][i];
00303 for(j=1;j<=nvar;j++)
00304 ex[nExp].dS[msP][i]+=ex[nExp].dyds[msP][i][j]*ex[nExp].dS[msP-1][j];
00305 for(j=1;j<=npar;j++)
00306 ex[nExp].dS[msP][i]+=ex[nExp].dydp[msP][i][j]*ex[nExp].dP[j];
00307 }
00308 }
00309
00310 }
00311 }
00312
00313 void solvLin(Glob *globs,GlobExp *ex,int computeCovar)
00314 {
00315
00316 long i,j,k,l,xind=0,yind=0,nL;
00317 long nvarFit=0;
00318 long nLpara=0;
00319 long nparL=0;
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
00328 double **Ma,*Ra;
00329 double **Me,*Re;
00330 double **Mg,*Rg;
00331
00332
00333 double *dX;
00334
00335 condense(globs,ex);
00336
00337
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
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
00385 if(globs->faktorLexist)
00386 {
00387
00388 aDim=aDim+nparL;
00389 eDim=eDim+nparL;
00390 fitDim=fitDim+nparL;
00391 }
00392
00393
00394
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
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
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
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
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 }
00457
00458
00459
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
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
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 }
00495
00496
00497
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
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
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 }
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
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
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
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
00581 if(globs->doP[index]=='L')
00582 {
00583 index_local_parameter++;
00584
00585 }
00586
00587 }
00588 for(j=1;j<=nrExp;j++)
00589 {
00590 pLmean[i]=pLmean[i]+ex[j].par[index];
00591
00592 }
00593 pLmean[i]=pLmean[i]/double(nrExp);
00594
00595 }
00596
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
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
00631
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
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
00660
00661
00662
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
00682 svdcmp(MMa,aDim,fitDim,diag,V);
00683
00684
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
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
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 }
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
00780 dX=(double*)malloc((fitDim+10)*sizeof(double));
00781
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
00790
00791
00792
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
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
00895 }
00896 else if(globs->minimiser==2)
00897 {
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
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
00929 for(i=1;i<=aDim;i++)
00930 B[i-1]=Ra[i];
00931
00932
00933 for(i=1;i<=fitDim;i++)
00934 {
00935 BL[i-1]=-1e25;
00936 BU[i-1]=1e25;
00937 X[i-1]=0.;
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
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
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
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 }
01027
01028
01029
01030
01031
01032
01033 extern "C" void
01034 myerr_ (int *messg, int &nmessg)
01035 {
01036 char *s = (char *) messg;
01037 cerr << endl;
01038 for (int i = 0; i < nmessg; i++)
01039 cerr << s[i];
01040 cerr << endl;
01041 }
01042
01043
01044
01045
01046
01047
01048
01049 extern "C" void
01050 xerr_ (char *msg, int &nerr, int &iert, int &ni, int &i1,
01051 int &i2, int &nr, double &r1, double &r2)
01052 {
01053 cerr << msg << endl;
01054 if (ni > 0)
01055 cerr << i1 << " ";
01056 if (ni > 1)
01057 cerr << i2 << " ";
01058 if (nr > 0)
01059 cerr << r1 << " ";
01060 if (nr > 1)
01061 cerr << r2 << " ";
01062 if (ni > 0 || nr > 0)
01063 cerr << endl;
01064 if (iert >= 2)
01065 cerr << endl;
01066 }