/home/peifer/diffit/modules/solvLin.cc

Go to the documentation of this file.
00001 #include<iostream>
00002 #include<fstream>
00003 #include<math.h>
00004 #include<stdlib.h>
00005 
00006 //#define PRINTDIMENSION
00007 #include "../def.h"
00008 #include "../model.h"
00009 #include "../nr.h"
00010 
00011 using namespace std;
00012 
00013 //#define PRINTDIMENSION
00014 #define PRINTBIGSYS
00015 //#define PRINTSTEP
00016 
00017 //Sub-module minimiser
00018 
00019 /* Black box routine for linear least squares problem   */
00020 /* with linear equality and inequality constraints      */
00021 /* ACM TOMS 8,3 (1982), 323-333                         */
00022 /* edited for IBM RS/6000 compatibility by TM, 96/02/22 */
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   //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 }
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   //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 }
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;  // # 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 }
01027 
01028 /******************************************************************************
01029  * error routines called from Fortran modules constr.f and odessa.f
01030  */
01031 
01032 // for constr.f
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 /* for odessa.f
01044  * nerr ignored
01045  * iert==2 checked, otherwise ignored
01046  * e.g.: CALL XERR ('AT T=R1, MXSTEP=I1 STEPS WERE TAKEN... ',
01047      1   201, 2, 1, MXSTEP, 0, 1, TN, ZERO)
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 }

Generated on Mon Jan 29 17:09:13 2007 for Diffit by  doxygen 1.4.6