/home/peifer/diffit/modules/computeRight.cc

Go to the documentation of this file.
00001 #include<iostream>
00002 #include<math.h>
00003 #include<stdlib.h>
00004 
00005 #include "../def.h"
00006 #include "../model.h"
00007 #include "../nr.h"
00008 
00009 using namespace std;
00010 
00011 // objective Function
00012 double objectiveF(Glob *globs,GlobExp *ex)
00013 {
00014   int i,j;
00015   double df,f;
00016   
00017   f=0;
00018   for (i=1; i<=ex->nMeasure; ++i)
00019     {
00020       for (j=1; j<=ex->nobs; ++j) 
00021         {
00022           //residues calculated in function computerRight
00023           df = (ex->residues[i][j]*ex->residues[i][j])/
00024             (ex->sigma[i][j]!=0 ? (ex->sigma[i][j]*ex->sigma[i][j]) : 1.0);
00025 #ifdef DEBUGF
00026           if (df>0.1) 
00027             *dbg << "f increment " << df << " for i = " << i << ", j = " << j << "\n";
00028 #endif
00029           f += df;
00030         }
00031     }
00032   ex->objF=f;
00033 #ifdef DEBUGF
00034   *dbg << "nMeasure = " << ex->nMeasure << ", nobs = " << ex->nobs << ", f = " << f << "\n";
00035 #endif
00036   return(f);
00037 }
00038 
00039 double computeRight(Glob *globs,GlobExp *ex)
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 }

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