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
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
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
00044 long nPoints=ex->nPoints, nvar=ex->nvar;
00045 long nMeasure=ex->nMeasure, nobs=ex->nobs;
00046 long nP=globs->npar;
00047
00048
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
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 }