00001
00002
00004
00005 #include "SingleNetworkMinimizable.h"
00006
00008
00010
00011 SingleNetworkMinimizable::SingleNetworkMinimizable()
00012 {
00013
00014 }
00015
00016 SingleNetworkMinimizable::SingleNetworkMinimizable(Experiment *experiment,CReactionMover *reactionMover, bool rateConstantsOptimizable, bool initialChemConcOptimizable, bool logsInObjectiveFunction, double timeSeriesWeight, double rateConstantsWeight, double initialChemConcWeight, int nCells)
00017 {
00018 Initialize(experiment,reactionMover,nCells,rateConstantsOptimizable,initialChemConcOptimizable,logsInObjectiveFunction,timeSeriesWeight,rateConstantsWeight,initialChemConcWeight);
00019 }
00020
00021 SingleNetworkMinimizable::~SingleNetworkMinimizable()
00022 {
00023 delete [] initialNetworkData;
00024 delete [] currentNetworkData;
00025 this->reactionMover->Detach(cellObserver);
00026 delete cellObserver;
00027 for(int i = 0; i < m_pConversionFactors.size(); i++)
00028 {
00029 delete m_pConversionFactors[i];
00030 }
00031 }
00032
00033
00034
00035
00036
00037
00038 void SingleNetworkMinimizable::Initialize(Experiment *experiment, CReactionMover *reactionMover, int nCells, bool rateConstantsOptimizable, bool initialChemConcOptimizable, bool logsInObjectiveFunction, double timeSeriesWeight, double rateConstantsWeight, double initialChemConcWeight)
00039 {
00040 this->experiment = experiment;
00041 this->reactionMover = reactionMover;
00042 this->nCells = nCells;
00043 this->rateConstantsOptimizable = rateConstantsOptimizable;
00044 this->initialChemConcOptimizable = initialChemConcOptimizable;
00045 this->logsInObjectiveFunction = logsInObjectiveFunction;
00046
00047 this->SetTimeSeriesWeight(timeSeriesWeight);
00048 this->SetRateConstantsWeight(rateConstantsWeight);
00049 this->SetInitialChemConcWeight(initialChemConcWeight);
00050
00051
00052 this->experiment->ReadData();
00053
00054 this->experiment->ReadForcingData();
00055
00056 int nChem = this->experiment->GetReactionNetwork()->GetNumberOfChemicals();
00057 int nRateConstants = this->experiment->GetReactionNetwork()->GetNumberOfRateConstants();
00058 int intTime = ((int) this->experiment->GetChemicalTimeSeriesData()->GetLatestDataTime()) + 2;
00059 cellObserver = new CellAverageObserver(nChem,intTime,nCells);
00060
00061
00062 this->reactionMover->Attach(cellObserver);
00063
00064 initialNetworkData = new double[nRateConstants + nChem];
00065 currentNetworkData = new double[nRateConstants + nChem];
00066
00067
00068
00069
00070 for(int l = 0; l < nChem; l++)
00071 {
00072 if((experiment->GetChemicalTimeSeriesData()->GetTimeSeries(l)).size() > 0 )
00073 {
00074 m_pConversionFactors.push_back(new CConversionFactor(true,false,1.0));
00075 }
00076 else
00077 {
00078 m_pConversionFactors.push_back(new CConversionFactor(false));
00079 }
00080 }
00081
00082 m_iNBFactors = 0;
00083
00084
00085
00086 for(int i = 0; i < nRateConstants; i++)
00087 {
00088 initialNetworkData[i] = this->experiment->GetReactionNetwork()->GetRateConstant(i)->GetInitialValue();
00089 currentNetworkData[i] = initialNetworkData[i];
00090 }
00091 for(int i = nRateConstants; i < nRateConstants + nChem; i++)
00092 {
00093 initialNetworkData[i] = this->experiment->GetReactionNetwork()->GetChemical(i-nRateConstants)->GetInitialAmount();
00094 currentNetworkData[i] = initialNetworkData[i];
00095 }
00096
00097 int localResidualsSize = 0;
00098 if(this->timeSeriesWeight > 0.0) {localResidualsSize += this->experiment->GetChemicalTimeSeriesData()->GetNDataPoints();}
00099 if(this->rateConstantsWeight > 0.0) {localResidualsSize += nRateConstants;}
00100 if(this->initialChemConcWeight > 0.0) {localResidualsSize += nChem;}
00101
00102 Allocate(localResidualsSize);
00103
00104 }
00105
00106 double SingleNetworkMinimizable::ComputeES(double *parameters, double T)
00107 {
00108 m_dEnergyLastComputed = this->ComputeResiduals(parameters);
00109 double S = this->EntropyShift(T) - 0.5*m_dEntropyLastComputed;
00110 m_dEntropyLastComputed = S;
00111 return m_dEnergyLastComputed;
00112 }
00113
00114
00115
00116
00117
00118
00119 double SingleNetworkMinimizable::F(double *parameters, double T)
00120 {
00121 m_dEnergyLastComputed = this->ComputeResiduals(parameters);
00122 double S = -0.5*m_dEntropyLastComputed;
00123 m_dEntropyLastComputed = S;
00124 return (m_dEnergyLastComputed - T*m_dEntropyLastComputed);
00125 }
00126
00127 double SingleNetworkMinimizable::F0(double *parameters, double T)
00128 {
00129 m_dEnergyLastComputed = this->ComputeResiduals(parameters);
00130 double S = this->EntropyShift(T) - 0.5*m_dEntropyLastComputed;
00131 m_dEntropyLastComputed = S;
00132 return (m_dEnergyLastComputed - T*m_dEntropyLastComputed);
00133 }
00134
00135 double SingleNetworkMinimizable::ComputeResiduals(double *parameters)
00136 {
00137 m_iNBFactors = 0;
00138 m_dEntropyLastComputed = 0.0;
00139 double partialEntropy = 0.0;
00140
00141 double objFuncVal = 0.0;
00142 int nRateConstants = experiment->GetReactionNetwork()->GetNumberOfRateConstants();
00143 int nChem = experiment->GetReactionNetwork()->GetNumberOfChemicals();
00144
00145
00146
00147 if(rateConstantsOptimizable && initialChemConcOptimizable)
00148 {
00149 for(int i = 0; i < nRateConstants; i++)
00150 {
00151 experiment->GetReactionNetwork()->GetRateConstant(i)->SetRateConstant(parameters[i]);
00152 currentNetworkData[i] = parameters[i];
00153 }
00154 for(int i = nRateConstants; i < nRateConstants + nChem; i++)
00155 {
00156 experiment->GetReactionNetwork()->GetChemical(i - nRateConstants)->SetAmount(parameters[i]);
00157 currentNetworkData[i] = parameters[i];
00158 }
00159 }
00160 else if(rateConstantsOptimizable)
00161 {
00162 for(int i = 0; i < nRateConstants; i++)
00163 {
00164 experiment->GetReactionNetwork()->GetRateConstant(i)->SetRateConstant(parameters[i]);
00165 currentNetworkData[i] = parameters[i];
00166 }
00167 experiment->GetReactionNetwork()->ChemicalReset();
00168 }
00169 else if(initialChemConcOptimizable)
00170 {
00171 for(int i = 0; i < nChem; i++)
00172 {
00173 experiment->GetReactionNetwork()->GetChemical(i)->SetAmount(parameters[i]);
00174 currentNetworkData[i+nRateConstants] = parameters[i];
00175 }
00176 experiment->GetReactionNetwork()->RateConstantReset();
00177 }
00178 else
00179 {
00180
00181
00182
00183
00184
00185 experiment->GetReactionNetwork()->RateConstantReset();
00186 experiment->GetReactionNetwork()->ChemicalReset();
00187 }
00188
00189
00190 cellObserver->ZeroConcentrations();
00191
00192 for(int cellCounter = 1; cellCounter <= nCells; cellCounter++)
00193 {
00194
00195
00196
00197
00198 if(initialChemConcOptimizable)
00199 {
00200 for(int i = nRateConstants; i < nRateConstants + nChem; i++)
00201 {
00202 experiment->GetReactionNetwork()->GetChemical(i - nRateConstants)->SetAmount(parameters[i]);
00203 }
00204 }
00205 else
00206 {
00207 experiment->GetReactionNetwork()->ChemicalReset();
00208 }
00209
00210
00211
00212
00213
00214 reactionMover->ResetStepStatistics(experiment->GetChemicalTimeSeriesData()->GetLatestDataTime());
00215
00216
00217 if(experiment->GetForcingData()->GetForcingVector().size() == 0)
00218 {
00219
00220 reactionMover->Move(0.0,0.0,experiment->GetReactionNetwork());
00221 reactionMover->Move(0.0,experiment->GetChemicalTimeSeriesData()->GetLatestDataTime(),experiment->GetReactionNetwork());
00222 }
00223 else
00224 {
00225 int forceCount = 0;
00226 double time = 0.0;
00227 int lastForce = experiment->GetForcingData()->GetForcingVector().size();
00228 double finaltime = experiment->GetChemicalTimeSeriesData()->GetLatestDataTime();
00229 double finalForceTime = experiment->GetForcingData()->GetLatestForcingTime();
00230 int index;
00231 int start,stop;
00232 double t1,t2,value;
00233 while(time < finalForceTime)
00234 {
00235 start = forceCount;
00236 stop = forceCount;
00237 t1 = experiment->GetForcingData()->GetForcingVector()[start]->GetTime();
00238
00239
00240 for(int i = start; i < lastForce; i++)
00241 {
00242 t2 = experiment->GetForcingData()->GetForcingVector()[i]->GetTime();
00243 if(t1 == t2)
00244 {
00245 }
00246 else
00247 {
00248 stop = i-1;
00249 for(int j = start; j <= stop; j++)
00250 {
00251 index = experiment->GetForcingData()->GetForcingVector()[j]->GetChemNumber();
00252 value = experiment->GetForcingData()->GetForcingVector()[j]->GetValue();
00253 experiment->GetReactionNetwork()->GetChemical(index)->SetAmount(value);
00254 }
00255 break;
00256 }
00257 }
00258 forceCount = stop+1;
00259
00260 reactionMover->Move(time,t1,experiment->GetReactionNetwork());
00261 time = t1;
00262 }
00263
00264
00265 t2 = finalForceTime;
00266 stop = lastForce-1;
00267 for(int i = lastForce-1; i > 0; i--)
00268 {
00269 t1 = experiment->GetForcingData()->GetForcingVector()[i]->GetTime();
00270 if(t1 == t2)
00271 {
00272 }
00273 else
00274 {
00275 start = i+1;
00276 for(int j = start; j <= stop; j++)
00277 {
00278 index = experiment->GetForcingData()->GetForcingVector()[j]->GetChemNumber();
00279 value = experiment->GetForcingData()->GetForcingVector()[j]->GetValue();
00280 experiment->GetReactionNetwork()->GetChemical(index)->SetAmount(value);
00281 }
00282 break;
00283 }
00284 }
00285
00286
00287 if(finaltime > time)
00288 {
00289 reactionMover->Move(time,finaltime,experiment->GetReactionNetwork());
00290 }
00291 time = finaltime;
00292 }
00293 }
00294
00295
00296
00297
00298 for(int chemCount = 0; chemCount < nChem; chemCount++)
00299 {
00300 double num = 0.0;
00301 double den = 0.0;
00302
00303
00304 if(m_pConversionFactors[chemCount]->IsFactorNeeded())
00305 {
00306 if(m_pConversionFactors[chemCount]->IsFactorFixed())
00307 {
00308 num = m_pConversionFactors[chemCount]->GetFactorValue();
00309 den = 1.0;
00310 }
00311 else
00312 {
00313 m_iNBFactors++;
00314 int dataSize = (experiment->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)).size();
00315
00316 for(int nData = 0; nData < dataSize; nData++)
00317 {
00318 double time = experiment->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetTime();
00319 double data = experiment->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetDatum();
00320 double error = experiment->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetError();
00321 double sim = (cellObserver->GetAverageConcentration())[chemCount][int(time)];
00322 num += ((2*sim)/(error*error))*data;
00323 den += ((2*sim)/(error*error))*sim;
00324
00325 partialEntropy += ((sim*sim)/(error*error));
00326 }
00327 m_dEntropyLastComputed += log(partialEntropy);
00328 }
00329 }
00330 else
00331 {
00332 num = 1.0;
00333 den = 1.0;
00334 }
00335
00336 double factor = num/den;
00337
00338 m_pConversionFactors[chemCount]->SetFactorValue(factor);
00339
00340 }
00341
00342
00343 int nextAvailableIndex = 0;
00344 if(timeSeriesWeight > 0.0)
00345 {
00346
00347 for(int chemIt = 0; chemIt < (experiment->GetChemicalTimeSeriesData()->GetTimeVector()).size(); chemIt++)
00348 {
00349
00350 double targetTime = (experiment->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetTime();
00351 int targetChem = (experiment->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetChemNumber();
00352 int vectorIndex = (experiment->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetVectorIndex();
00353
00354 double data = (experiment->GetChemicalTimeSeriesData()->GetTimeSeries(targetChem))[vectorIndex]->GetDatum();
00355 double sigma = (experiment->GetChemicalTimeSeriesData()->GetTimeSeries(targetChem))[vectorIndex]->GetError();
00356 double sim = (cellObserver->GetAverageConcentration())[targetChem][int(targetTime)];
00357
00358 if(m_pConversionFactors[targetChem]->IsFactorNeeded())
00359 {
00360 sim = sim*(m_pConversionFactors[targetChem]->GetFactorValue());
00361 }
00362
00363 residuals[nextAvailableIndex] = (sqrt(timeSeriesWeight/2.0))*((data - sim)/sigma);
00364 objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
00365 nextAvailableIndex++;
00366 }
00367 }
00368 if(rateConstantsWeight > 0.0)
00369 {
00370
00371 if(logsInObjectiveFunction)
00372 {
00373 for(int i = 0; i < nRateConstants; i++)
00374 {
00375 double sigma = experiment->GetReactionNetwork()->GetRateConstant(i)->GetErrorInInitialValue();
00376 double sigmalog = sigma/initialNetworkData[i];
00377 residuals[nextAvailableIndex] = (sqrt(rateConstantsWeight/2.0))*(log(initialNetworkData[i]/currentNetworkData[i])/sigmalog);
00378 objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
00379 nextAvailableIndex++;
00380 }
00381 }
00382 else
00383 {
00384 for(int i = 0; i < nRateConstants; i++)
00385 {
00386 double sigma = experiment->GetReactionNetwork()->GetRateConstant(i)->GetErrorInInitialValue();
00387 residuals[nextAvailableIndex] = (sqrt(rateConstantsWeight/2.0))*((initialNetworkData[i] - currentNetworkData[i])/sigma);
00388 objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
00389 nextAvailableIndex++;
00390 }
00391 }
00392 }
00393 if(initialChemConcWeight > 0.0)
00394 {
00395
00396 for(int i = 0; i < nChem; i++)
00397 {
00398 double sigma = experiment->GetReactionNetwork()->GetChemical(i)->GetErrorInInitialAmount();
00399 residuals[nextAvailableIndex] = (sqrt(initialChemConcWeight/2.0))*((initialNetworkData[i+nRateConstants] - currentNetworkData[i+nRateConstants])/sigma);
00400 objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
00401 nextAvailableIndex++;
00402 }
00403 }
00404
00405
00406 Notify();
00407
00408 return objFuncVal;
00409 }
00410
00411 int SingleNetworkMinimizable::GetNParameters()
00412 {
00413 int nRateConstants = experiment->GetReactionNetwork()->GetNumberOfRateConstants();
00414 int nChem = experiment->GetReactionNetwork()->GetNumberOfChemicals();
00415 int nParameters;
00416
00417 if(rateConstantsOptimizable && initialChemConcOptimizable)
00418 {
00419 nParameters = nRateConstants + nChem;
00420 }
00421 else if(rateConstantsOptimizable)
00422 {
00423 nParameters = nRateConstants;
00424 }
00425 else if(initialChemConcOptimizable)
00426 {
00427 nParameters = nChem;
00428 }
00429 else
00430 {
00431 nParameters = 0;
00432 }
00433
00434 return nParameters;
00435 }
00436
00437 double SingleNetworkMinimizable::GetParameter(int parIndex)
00438 {
00439 int nP = GetNParameters();
00440
00441 if(parIndex > nP)
00442 {
00443 throw std::runtime_error("Parameter index is off parameter list.");
00444 }
00445 else
00446 {
00447 return currentNetworkData[parIndex];
00448 }
00449
00450 return 0.0;
00451 }