#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include "../def.h"
#include "../model.h"
#include "../nr.h"
Go to the source code of this file.
Defines | |
#define | PRINTBIGSYS |
Functions | |
int | lsei_ (double *w, long &mdw, long &me, long &ma, long &mg, long &n, double *prgopt, double *x, double *rnorme, double *rnorml, long *mode, double *ws, long *ip) |
int | e04ncf_ (long &M, long &N, long &NCLIN, long &LDC, long &LDA, double *C, double *BL, double *BU, double *CVEC, long *ISTATE, long *KX, double *X, double *A, double *B, long &ITER, double &OBJ, double *CLAMBDA, long *IWORK, long &LIWORK, double *WORK, long &LWORK, long &IFAIL) |
int | setopte04ncf_ () |
void | condense (Glob *globs, GlobExp *ex) |
Condenses each experiment. | |
void | decondense (Glob *globs, GlobExp *ex, double *dX) |
Decondenses each experiment. | |
void | solvLin (Glob *globs, GlobExp *ex, int computeCovar) |
void | myerr_ (int *messg, int &nmessg) |
void | xerr_ (char *msg, int &nerr, int &iert, int &ni, int &i1, int &i2, int &nr, double &r1, double &r2) |
|
Definition at line 14 of file solvLin.cc. |
|
Condenses each experiment.
Definition at line 40 of file solvLin.cc. References dmatrix(), Glob::elastic, GlobExp::me, GlobExp::mg, GlobExp::nobs, Glob::npar, GlobExp::nPoints, Glob::nrExp, and GlobExp::nvar. Referenced by solvLin(). 00042 { 00043 long nExp,msP,i,j,k,l; 00044 long nMeas,nvar,nP; 00045 long me,mg,nobs,nms; 00046 00047 //weakens continuty constraints 00048 double elastic=globs->elastic; 00049 double **Ea,**Ee,**Eg; //save the previous E-matrices 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 //initialisation of backward recursion 00066 //------------------------------------ 00067 //least squares: 00068 for(i=1;i<=nMeas;i++) 00069 { 00070 for(j=1;j<=nobs;j++) 00071 { 00072 //ua 00073 ex[nExp].ua[(j-1)*nMeas+i]=ex[nExp].residues[i][j]/ex[nExp].sigma[i][j]; 00074 //Ea 00075 for(k=1;k<=nvar;k++) 00076 { 00077 //find out which datapoints are falling in the last MS-interval 00078 //the other derivatives are zero 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 //Pa 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 }// end of least squares init. 00097 00098 //equality constraints 00099 for(i=1;i<=me;i++) 00100 { 00101 //ue 00102 ex[nExp].ue[i]=ex[nExp].r2[i]; 00103 //Ee 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 //Pe 00110 for(j=1;j<=nP;j++) 00111 ex[nExp].Pe[i][j]=ex[nExp].dR2dp[i][j]; 00112 }//end of equality constraints 00113 00114 //inequality constraints 00115 for(i=1;i<=mg;i++) 00116 { 00117 //ug 00118 ex[nExp].ug[i]=ex[nExp].r3[i]; 00119 //Ee 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 //Pe 00126 for(j=1;j<=nP;j++) 00127 ex[nExp].Pg[i][j]=ex[nExp].dR3dp[i][j]; 00128 }//end of inequality constraints 00129 //initialisation of backward recursion complete 00130 00131 00132 //backward recursion 00133 //----------------- 00134 00135 //loop over all multiple shooting intervals 00136 for(msP=nms;msP>=2;msP--) 00137 { 00138 //least squares 00139 00140 for(i=1;i<=nMeas;i++) 00141 { 00142 for(j=1;j<=nobs;j++) 00143 { 00144 //ua 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 //Ea 00148 for(k=1;k<=nvar;k++) 00149 { 00150 //find out which datapoints are falling in the MS-interval 00151 //the other derivatives are zero 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 //Pa 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 } //end of least squares 00174 00175 //equality constraints 00176 for(i=1;i<=me;i++) 00177 { 00178 //ue 00179 for(k=1;k<=nvar;k++) 00180 ex[nExp].ue[i]+=Ee[i][k]*ex[nExp].h[msP][k]*elastic; 00181 //Ee 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 //Pe 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 }//end of equality constraints 00195 00196 //inequality constraints 00197 for(i=1;i<=mg;i++) 00198 { 00199 //ue 00200 for(k=1;k<=nvar;k++) 00201 ex[nExp].ug[i]+=Eg[i][k]*ex[nExp].h[msP][k]*elastic; 00202 //Ee 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 //Pe 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 }//end of inequality constraints 00216 00217 //updating new Ea,Ee,Ea 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 }//end loop over multiple shooting intervals 00229 00230 //clearing memory 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 }//end loop over all experiments 00235 }
|
|
Decondenses each experiment.
Definition at line 238 of file solvLin.cc. References GlobExp::dS, Glob::elastic, FALSE, Glob::npar, GlobExp::nPoints, Glob::nrExp, GlobExp::nvar, and Glob::y0fix. 00240 { 00241 long nExp,msP,i,j,k,l; 00242 long nvar,npar; 00243 long nms,ind=0; 00244 00245 //weakens continuty constraints 00246 double elastic=globs->elastic; 00247 00248 npar=globs->npar; 00249 nvar=ex[1].nvar; 00250 00251 //initialise 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 //end initialise 00291 00292 //forward iteration 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 }//end loop over experiments 00311 }
|
|
|
|
|
|
Definition at line 1034 of file solvLin.cc. 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 }
|
|
|
|
Definition at line 313 of file solvLin.cc. References condense(), FALSE, Glob::npar, Glob::nrExp, GlobExp::nvar, and Glob::y0fix. 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 1050 of file solvLin.cc. 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 }
|