- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/modules/dampIt.cc File Reference

`#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.

## Functions
void | intODE (GlobExp *ex, Glob *globs, int doDerivs, int doPlot, long expNr)
double | computeRight (Glob *globs, GlobExp *ex)
void | solvLin (Glob *globs, GlobExp *ex, int computeCovar)
double | Norm (Glob *globs, GlobExp *ex, double ***dS, double **dP, double uS)
double | subIt (Glob *globs, GlobExp *ex, double t, double normdx, double **p0, double ***s0, double **dp, double ***ds, double uS)
double | dampIt (Glob *globs, GlobExp *ex, double **P0, double ***S0, double **dP, double ***dS, double *uS)

## Function Documentation

double computeRight | ( | Glob * | globs,
|  | GlobExp * | ex
| ) |

| Definition at line 39 of file computeRight.cc. References GlobExp::h, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, and GlobExp::nvar.
    00040 {
    00041     long i,j,k;
    00042
    00043     // calculate discrepancies, residues and constraints
    00044     long nPoints=ex->nPoints, nvar=ex->nvar;
    00045     long nMeasure=ex->nMeasure, nobs=ex->nobs;
    00046     long nP=globs->npar;
    00047
    00048     // no continuity constraints for first and last mesh point
    00049     for (j=1;j<=nvar;++j) {
    00050         ex->h[1][j]=0;
    00051         ex->h[nPoints][j]=0;
    00052     }
    00053     for (i=2;i<nPoints;++i)
    00054       {
    00055         for (j=1;j<=nvar;++j)
    00056           ex->h[i][j]=ex->yComp[i][j]-ex->yTry[i][j];
    00057       }
    00058     for (i=1;i<=nMeasure;++i)
    00059       {
    00060         for (j=1;j<=nobs;++j)
    00061           ex->residues[i][j]=ex->yPred[i][j]-ex->yMeasure[i][j];
    00062       }
    00063
    00064     //extra equality and inequality constraints
    00065     R2(globs,ex, TRUE);
    00066     R3(globs,ex ,TRUE);
    00067
    00068 #ifdef PRINTDERIVS
    00069     for (i=1; i<=nPoints; ++i) {
    00070         *dbg << "dyds[" << i << "] = (";
    00071         for (j=1; j<=nvar; ++j) {
    00072             for (k=1; k<=nvar; ++k)
    00073                 *dbg << " " << ex->dyds[i][j][k];
    00074             *dbg << ", ";
    00075         }
    00076         *dbg << ")\ndydp[" << i << "] = (";
    00077         for (j=1; j<=nvar; ++j) {
    00078             for (k=1; k<=nP; ++k)
    00079                 *dbg << " " << ex->dydp[i][j][k];
    00080             *dbg << ", ";
    00081         }
    00082         *dbg << ")\n";
    00083     }
    00084     for (i=1; i<=nMeasure; ++i) {
    00085         *dbg << "dmds[" << i << "] = (";
    00086         for (j=1; j<=nobs; ++j) {
    00087             for (k=1; k<=nvar; ++k)
    00088                 *dbg << " " << ex->dmds[i][j][k];
    00089             *dbg << ", ";
    00090         }
    00091         *dbg << ")\ndmdp[" << i << "] = (";
    00092         for (j=1; j<=nobs; ++j) {
    00093             for (k=1; k<=nP; ++k)
    00094                 *dbg << " " << ex->dmdp[i][j][k];
    00095             *dbg << ", ";
    00096         }
    00097         *dbg << ")\n";
    00098     }
    00099     for (i=1; i<=ex->me; ++i) {
    00100         *dbg << "dR2ds[" << i << "] = (";
    00101         for (j=1; j<=nPoints; ++j) {
    00102             for (k=1; k<=nvar; ++k)
    00103                 *dbg << " " << ex->dR2ds[i][j][k];
    00104             *dbg << ", ";
    00105         }
    00106         *dbg << ")\n";
    00107         *dbg << "dR2dp[" << i << "] = (";
    00108         for (j=1; j<=nP; ++j)
    00109             *dbg << " " << ex->dR2dp[i][j];
    00110         *dbg << ")\n";
    00111     }
    00112     for (i=1; i<=ex->mg; ++i) {
    00113         *dbg << "dR3ds[" << i << "] = (";
    00114         for (j=1; j<=nPoints; ++j) {
    00115             for (k=1; k<=nvar; ++k)
    00116                 *dbg << " " << ex->dR3ds[i][j][k];
    00117             *dbg << ", ";
    00118         }
    00119         *dbg << ")\n";
    00120         *dbg << "dR3dp[" << i << "] = (";
    00121         for (j=1; j<=nP; ++j)
    00122             *dbg << " " << ex->dR3dp[i][j];
    00123         *dbg << ")\n";
    00124     }
    00125     dbg->flush();
    00126 #endif
    00127
    00128     return(objectiveF(globs,ex));
    00129 }

double dampIt | ( | Glob * | globs,
|  | GlobExp * | ex,
|  | double ** | P0,
|  | double *** | S0,
|  | double ** | dP,
|  | double *** | dS,
|  | double * | uS
| ) |

| Definition at line 155 of file dampIt.cc. References Glob::chisq, FALSE, Norm(), Glob::silent, subIt(), TRUE, and Glob::wquer.
    00157 {
    00158   int etaOk=FALSE,i=1;
    00159   double tau=0.95;
    00160   double t;
    00161   double tau_min=0.1;
    00162   double eta0=1.;
    00163   double eta2=1.8; //should be parsed
    00164
    00165   double normdx=Norm(globs,ex,dS,dP,*uS);
    00166   // damping algorithm
    00167
    00168   if(globs->wquer<0)
    00169     {
    00170       globs->wquer=1.*eta0/normdx;
    00171     }
    00172
    00173   else
    00174     t=eta0/(globs->wquer*normdx);
    00175
    00176   if (t > tau)
    00177     {
    00178       t=1;
    00179     }
    00180   else if(t <= tau_min)
    00181     {
    00182       t=tau_min;
    00183     }
    00184
    00185   globs->wquer=subIt(globs, ex, t,normdx,P0, S0, dP, dS,*uS);
    00186   if(globs->silent!=TRUE)
    00187     {
    00188       cout << "Start damping:\n";
    00189       cout << "#" << i << " (Predicted)  t = " << t << "\t";
    00190       cout << "chisq=" << globs->chisq;
    00191       cout << endl;
    00192     }
    00193
    00194   etaOk=((globs->wquer*t*normdx)<eta2);
    00195   i++;
    00196
    00197   while (!etaOk)
    00198     {
    00199       t=eta0/(normdx*globs->wquer);
    00200
    00201       if (t > tau)
    00202         {
    00203           t=1;
    00204           etaOk=TRUE;
    00205           continue;
    00206         }
    00207       else if(t < tau_min)
    00208         {
    00209               t=tau_min;
    00210               etaOk=TRUE;
    00211               if(globs->silent!=TRUE)
    00212                 {
    00213                   cout << "#" << i <<  " (Corrected)  t = " << t << "\t";
    00214                   cout << "\nMinimal damping factor tau_min= " << tau_min << " reached.\n";
    00215                 }
    00216               continue;
    00217         }
    00218       globs->wquer=subIt(globs, ex, t, normdx,P0, S0, dP, dS,*uS);
    00219       if(globs->silent!=TRUE)
    00220         {
    00221           cout << "#" << i <<  " (Corrected)  t = " << t << "\t";
    00222           cout << "chisq=" << globs->chisq;
    00223           cout << endl;
    00224         }
    00225       etaOk=((globs->wquer*t*normdx)<eta2);
    00226       i++;
    00227       if(i>=10)
    00228         {
    00229           cerr << "Too many corrections.\n";
    00230           throw 1;
    00231         }
    00232     }
    00233   return(t);
    00234 }

void intODE | ( | GlobExp * | ex,
|  | Glob * | globs,
|  | int | doDerivs,
|  | int | doPlot,
|  | long | expNr
| ) |

| Definition at line 313 of file intODE.cc. References d3tensor(), dmatrix(), GlobExp::dmdp, GlobExp::dmds, dvector(), GlobExp::dydp, GlobExp::dyds, maxTableLen(), GlobExp::mesh, GlobExp::nMeasure, GlobExp::nobs, Glob::npar, GlobExp::nPoints, GlobExp::nvar, and GlobExp::xMeasure.
    00314 {
    00315   long i,j,k,ksav,l,m,n,nok,nbad;
    00316   long nobs=ex->nobs;
    00317   long nPoints=ex->nPoints, nvar=ex->nvar, nMeasure=ex->nMeasure;
    00318   double *mesh=ex->mesh, *xMeasure=ex->xMeasure;
    00319   long tableLen=maxTableLen(mesh, nPoints, xMeasure, nMeasure)+1;
    00320   double temp, h,ymin,ymax;
    00321
    00322   // allocate memory
    00323   double *ys=dvector(1,nvar);
    00324   double *xTab=dvector(1,tableLen);
    00325   double **yTab=dmatrix(1,tableLen,1,nvar);
    00326   double ***dTabds, ***dTabdp;
    00327   char name[50];
    00328   ofstream gnuout;
    00329   ofstream outout;
    00330
    00331
    00332   //gnuout.open("gnuout.dat");
    00333
    00334   // compute derivatives?
    00335   // if not, don't write to graphics either
    00336   // (no genuine integration, but only estimation of omega)
    00337   if (doDerivs)
    00338     {
    00339       dTabds=d3tensor(1,tableLen,1,nvar,1,nvar);
    00340       dTabdp=d3tensor(1,tableLen,1,nvar,1,globs->npar);
    00341     }
    00342   else
    00343     {
    00344       dTabds=NULL;
    00345       dTabdp=NULL;
    00346     }
    00347
    00348
    00349   // first mesh point
    00350   if (doDerivs) {
    00351     for (i=1; i<=nvar; ++i)
    00352       {
    00353         for (j=1; j<=nvar; ++j)
    00354           {
    00355             ex->dyds[1][i][j] = (i==j) ? 1 : 0;
    00356             ex->dmds[1][i][j] = (i==j) ? 1 : 0;
    00357           }
    00358         for (j=1; j<=globs->npar; ++j)
    00359           {
    00360             ex->dydp[1][i][j] = 0;
    00361             ex->dmdp[1][i][j] = 0;
    00362           }
    00363       }
    00364   }
    00365   for (j=1; j<=nvar; ++j)
    00366     ex->yComp[1][j]=ex->yTry[1][j];
    00367
    00368
    00369   if (doDerivs && mesh[1]==xMeasure[1])
    00370     {
    00371       observation(globs,ex,mesh[1],ex->yTry[1],ex->yPred[1],ex->par,ex->dmds[1],ex->dmdp[1]);
    00372     }
    00373
    00374   // integrate
    00375   k=1;
    00376   ksav=1;
    00377   if (mesh[1]==xMeasure[1])
    00378     {
    00379       k=2;
    00380       ksav=2;
    00381     }
    00382   for (i=1; i< nPoints; ++i)
    00383     {   // given starting values
    00384       for (j=1; j<=nvar; ++j)
    00385         ys[j]=ex->yTry[i][j];
    00386       if (xMeasure[k]<mesh[i])
    00387         {
    00388           cerr << "fatal: intODE: xMeasure table corrupted (mesh: "
    00389                << mesh[i] << ", xMeasure: " << xMeasure[k] << ")\n";
    00390           throw 1;
    00391         }
    00392       // build table of measuring points
    00393       n=1;
    00394       while (k <= nMeasure && xMeasure[k]<=mesh[i+1])
    00395         xTab[n++]=xMeasure[k++];
    00396       xTab[n]=mesh[i+1];
    00397 #ifdef PRINTMESH
    00398       *dbg << "Mesh point " << mesh[i] << ": xTab = (";
    00399       for (j=1; j<=n; ++j)
    00400         *dbg << '\t' << xTab[j];
    00401       *dbg << ")\n";
    00402       dbg->flush();
    00403 #endif
    00404       //integration
    00405       tabulateValues(globs,ex,mesh[i],mesh[i+1],xTab,n,ys,yTab, dTabds, dTabdp,ex->dyds[i+1],ex->dydp[i+1]);
    00406
    00407 #ifdef PRINTINTEGRATE
    00408       *dbg << "Integration returned\t  (";
    00409       for (j=1; j<=n; ++j)
    00410         {
    00411           for (l=1; l<=nvar; ++l)
    00412             *dbg << " " << yTab[j][l];
    00413           *dbg << ",";
    00414         }
    00415       *dbg << ")\n";
    00416       dbg->flush();
    00417 #endif
    00418       for (j=1; j<=ex->nvar; ++j)
    00419         ex->yComp[i+1][j]=ys[j];   // store value at next mesh point
    00420
    00421       for (j=ksav; j<k; ++j)          // store value at measuring points
    00422         {
    00423           for (l=1; l<=ex->nobs; ++l)
    00424             {
    00425               ex->yPred[j][l]=yTab[j-ksav+1][l];
    00426               if (doDerivs)
    00427                 {
    00428                   for (m=1; m<=ex->nvar; ++m)
    00429                     ex->dmds[j][l][m]=dTabds[j-ksav+1][l][m];
    00430                   for (m=1; m<=globs->npar; ++m)
    00431                     ex->dmdp[j][l][m]=dTabdp[j-ksav+1][l][m];
    00432                 }
    00433             }
    00434         }
    00435       ksav=k;
    00436     }
    00437
    00438       //GnuPlotting
    00439   long mind=2;
    00440   if(doDerivs && (globs->noGnu==FALSE) && doPlot==TRUE)
    00441     {
    00442       for(l=1;l<=ex->nobs;l++)
    00443         {
    00444           if(globs->gnuindex <= globs->ngnu)
    00445             {
    00446               mind=2;
    00447               ofstream gnuout;
    00448               gnuout.open("gnuout.dat");
    00449               ymax=ex->yPred[1][l];
    00450               ymin=ex->yPred[1][l];
    00451               for(j=1;j<=ex->nMeasure;j++)
    00452                 {
    00453                   gnuout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl;
    00454                   if(mesh[mind] <= ex->xMeasure[j])
    00455                     {
    00456                       observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL);
    00457                       gnuout << mesh[mind] << "\t" << ys[l] << endl << endl;
    00458                       observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL);
    00459                       gnuout << mesh[mind] << "\t" << ys[l] << endl;
    00460                       mind++;
    00461                     }
    00462                   //determine y-range
    00463                   if(ymax <= ex->yPred[j][l])
    00464                     ymax=ex->yPred[j][l];
    00465                   if(ymin >= ex->yPred[j][l])
    00466                     ymin=ex->yPred[j][l];
    00467                   if(ymax <= ex->yMeasure[j][l])
    00468                     ymax=ex->yMeasure[j][l];
    00469                   if(ymin >= ex->yMeasure[j][l])
    00470                     ymin=ex->yMeasure[j][l];
    00471                 }
    00472               gnuout.flush();
    00473               gnuout.close();
    00474
    00475               fprintf(globs->gnuFp[globs->gnuindex],"reset\n");
    00476               fprintf(globs->gnuFp[globs->gnuindex],"set yrange[%f:%f]\n",ymin-0.01*ymin,ymax+0.01*ymax);
    00477               fprintf(globs->gnuFp[globs->gnuindex],"set xrange[%f:%f]\n",ex->xMeasure[1],ex->xMeasure[ex->nMeasure]);
    00478               fflush(globs->gnuFp[globs->gnuindex]);
    00479               fprintf(globs->gnuFp[globs->gnuindex],"plot \"gnuout.dat\" title \"\"w l 3,\"gnuout.dat\" u 1:3 title \"\" 1 \n");
    00480               fflush(globs->gnuFp[globs->gnuindex]);
    00481               //DELAY
    00482               for(long dd=1;dd<=1000000;dd++)
    00483                1+1;
    00484             }
    00485           sprintf(name,"Gnuout_Exp%d_Obs%d.dat",expNr,l);
    00486           outout.open(name);
    00487
    00488           globs->gnuindex++;
    00489           mind=2;
    00490           for(j=1;j<=ex->nMeasure;j++)
    00491             {
    00492
    00493               outout << ex->xMeasure[j] << "\t" << ex->yPred[j][l] << "\t" << ex->yMeasure[j][l] << endl;
    00494               if(mesh[mind] <= ex->xMeasure[j])
    00495                 {
    00496                   observation (globs,ex,mesh[mind],ex->yComp[mind],ys,ex->par,NULL,NULL);
    00497                   outout << mesh[mind] << "\t" << ys[l] << endl << endl;
    00498                   observation (globs,ex,mesh[mind],ex->yTry[mind],ys,ex->par,NULL,NULL);
    00499                   outout << mesh[mind] << "\t" << ys[l] << endl;
    00500                   mind++;
    00501                 }
    00502             }
    00503           outout.flush();
    00504           outout.close();
    00505           for(long dd=1;dd<=1000000;dd++)
    00506             1+1;
    00507         }
    00508     }
    00509
    00510   if (doDerivs) {
    00511     free_d3tensor(dTabdp,1,tableLen,1,nvar,1,globs->npar);
    00512     free_d3tensor(dTabds,1,tableLen,1,nvar,1,nvar);
    00513   }
    00514   free_dmatrix(yTab,1,tableLen,1,nvar);
    00515   free_dvector(xTab,1,tableLen);
    00516   free_dvector(ys,1,nvar);
    00517 }

double Norm | ( | Glob * | globs,
|  | GlobExp * | ex,
|  | double *** | dS,
|  | double ** | dP,
|  | double | uS
| ) |

| Definition at line 29 of file dampIt.cc. References Glob::doP, Glob::npar, GlobExp::nPoints, Glob::nrExp, and GlobExp::nvar. Referenced by dampIt().
    00030 {
    00031   double norm=0;
    00032   if (!dS)
    00033     {
    00034     cerr <<"dS==NULL in Norm\n";
    00035     throw 1;
    00036     }
    00037   long i,ivar,iexp,ip;
    00038   for (iexp=1;iexp<=globs->nrExp;++iexp)
    00039     {
    00040       for (i=1; i< ex[iexp].nPoints; ++i)
    00041         for (ivar=1; ivar<=ex[iexp].nvar; ++ivar)
    00042           norm+=pow(uS*dS[iexp][i][ivar],2);
    00043
    00044       for (ip=1;ip<=globs->npar; ++ip)
    00045         {
    00046           if(globs->doP[ip]=='L')
    00047             norm+=pow(dP[iexp][ip],2);
    00048         }
    00049     } // end of: for (iexp=1;iexp<=globs->nrExp;++iexp)
    00050   for (ip=1;ip<=globs->npar; ++ip)
    00051     {
    00052       if(globs->doP[ip]==TRUE)
    00053         norm+=pow(dP[1][ip],2);
    00054     }
    00055   return sqrt(norm);
    00056
    00057 } // Norm

void solvLin | ( | Glob * | globs,
|  | GlobExp * | ex,
|  | int | computeCovar
| ) |

| 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 }

double subIt | ( | Glob * | globs,
|  | GlobExp * | ex,
|  | double | t,
|  | double | normdx,
|  | double ** | p0,
|  | double *** | s0,
|  | double ** | dp,
|  | double *** | ds,
|  | double | uS
| ) |

| Definition at line 62 of file dampIt.cc. References Glob::npar, and Glob::nrExp. Referenced by dampIt().
    00064 {
    00065   long iexp,i,j;
    00066   long maxPoints=0,maxVar=0;
    00067   long nP=globs->npar;
    00068
    00069   for (i=1;i<=globs->nrExp;++i)
    00070     {
    00071       if (ex[i].nPoints>maxPoints)
    00072         maxPoints=ex[i].nPoints;
    00073       if (ex[i].nvar>maxVar)
    00074         maxVar=ex[i].nvar;
    00075     }
    00076
    00077   double **newdP=dmatrix(1,globs->nrExp,1,nP);                      // parameter updates
    00078   double ***newdS=d3tensor(1,globs->nrExp,1,maxPoints,1,maxVar); // yTry updates
    00079
    00080
    00081   for(iexp=1;iexp<=globs->nrExp;iexp++)
    00082     {
    00083       GlobExp *expi=&ex[iexp];
    00084       for(i=1; i< expi->nPoints; i++)
    00085         {
    00086           for(j=1; j<=expi->nvar; j++)
    00087             {
    00088               expi->yTry[i][j] = s0[iexp][i][j]+t*uS*ds[iexp][i][j];
    00089               // positiveness of initial values --> test
    00090               if(expi->yTry[i][j] < 0.0)
    00091                 expi->yTry[i][j] = 0.0;
    00092               expi->h[i][j]*=uS;
    00093             }
    00094         }
    00095       for(i=1; i<=nP; i++)
    00096         {
    00097           expi->par[i] = p0[iexp][i]+t*dp[iexp][i];
    00098         }
    00099     }
    00100
    00101   //integrate and compute value of objective function
    00102   globs->chisq=0.;
    00103   for (iexp=1;iexp<=globs->nrExp;iexp++)
    00104     {
    00105       intODE(&ex[iexp],globs,FALSE,FALSE,iexp);
    00106       globs->chisq+=computeRight(globs,&ex[iexp]);
    00107
    00108     }
    00109   *dbg  <<"subtotal chisq=" << globs->chisq << endl;
    00110   solvLin(globs,ex,FALSE);
    00111
    00112   for (iexp=1;iexp<=globs->nrExp;iexp++)
    00113     {
    00114       GlobExp *expi=&ex[iexp];
    00115       for (i=1; i< expi->nPoints; i++)
    00116         {
    00117           for (j=1; j<=expi->nvar; j++)
    00118             {
    00119               newdS[iexp][i][j]=ex[iexp].dS[i][j];
    00120             }
    00121         }
    00122       for(i=1; i<=nP; i++)
    00123         {
    00124         newdP[iexp][i]=ex[iexp].dP[i];
    00125         }
    00126     }
    00127   //new increments in newdP and newdS
    00128
    00129   //w-estimator
    00130   double sum=0;
    00131   for (iexp=1;iexp<=globs->nrExp;++iexp)
    00132     {
    00133       GlobExp *expi=&ex[iexp];
    00134       for (i=1; i< expi->nPoints; ++i)
    00135         for (j=1; j<=expi->nvar; ++j)
    00136           sum+=pow(newdS[iexp][i][j]-(1-t)*uS*ds[iexp][i][j],2);
    00137       for (i=1; i<=nP; ++i)
    00138         {
    00139           if(globs->doP[i]=='L')
    00140             sum+=pow(newdP[iexp][i]-(1-t)*dp[iexp][i],2);
    00141         }
    00142     }
    00143   for (i=1; i<=nP; ++i)
    00144     {
    00145       if(globs->doP[i]==TRUE)
    00146         sum+=pow(newdP[1][i]-(1-t)*dp[1][i],2);
    00147     }
    00148
    00149   free_d3tensor(newdS,1,globs->nrExp,1,maxPoints,1,maxVar);
    00150   free_dmatrix(newdP,1,globs->nrExp,1,nP);
    00151   return 2*sqrt(sum)/(pow(t*normdx,2));
    00152 } // subIt