Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Compound Members   File Members  

SingleNetworkMinimizable.cpp

Go to the documentation of this file.
00001 // SingleNetworkMinimizable.cpp: implementation of the SingleNetworkMinimizable class.
00002 //
00004 
00005 #include "SingleNetworkMinimizable.h"
00006 
00008 // Construction/Destruction
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 /*Experiment *SingleNetworkMinimizable::GetExperiment()
00034 {
00035         return this->experiment;
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         // Read in the data
00052         this->experiment->ReadData();
00053         // Read in forcing data (if any)
00054         this->experiment->ReadForcingData();
00055         // Initialize variables dependent on the Experiment class
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         // attach the observer
00062         this->reactionMover->Attach(cellObserver);
00063 
00064         initialNetworkData = new double[nRateConstants + nChem];
00065         currentNetworkData = new double[nRateConstants + nChem];
00066 
00067         // load up the vector of conversion factors with dummies assuming all data is 
00068         // expressed in terms of number of molecules; if we don't have data we don't 
00069         // need conversion factors
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         // store beginning network data - this allows us to restore the 
00085         // initial state of the network to recover from the beginning
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 double SingleNetworkMinimizable::EntropyShift(double T)
00115 {
00116         return 0.5*m_iNBFactors*log(2*3.14159265358979323846*T);
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         // put the parameters into the network
00146         // rate constants first, then initial chemical concentrations
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                 // do one run at the existing parameter values
00181                 // this ensures that multiple calls with the same parameters
00182                 // give the same cost, but it also ensures that only time
00183                 // series contribute to the cost - this is probably not what
00184                 // we want?
00185                 experiment->GetReactionNetwork()->RateConstantReset();
00186                 experiment->GetReactionNetwork()->ChemicalReset();
00187         }
00188 
00189         // zero the observer
00190         cellObserver->ZeroConcentrations();
00191 
00192         for(int cellCounter = 1; cellCounter <= nCells; cellCounter++)
00193         {
00194                 // make sure the concentration in the network is re-initialized, which
00195                 // is done differently whether you allow initial conc. optimization
00196                 // or not
00197                 //experiment->GetReactionNetwork()->ChemicalReset();
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                 // if there's no forcing data at all, we can just run from (0,tfinal). 
00210                 // Otherwise we have to stop running in between forcing changes.  This
00211                 // block might be able to be rearranged to be much faster. KSB 4/16/02
00212 
00213                 // reset step size statistics in either case
00214                 reactionMover->ResetStepStatistics(experiment->GetChemicalTimeSeriesData()->GetLatestDataTime());
00215                 
00216                 // now move through the rest
00217                 if(experiment->GetForcingData()->GetForcingVector().size() == 0)
00218                 {
00219                         // take data at time 0
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                                 // find the starting and stopping points
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                                 // move for the correct amount of time
00260                                 reactionMover->Move(time,t1,experiment->GetReactionNetwork());
00261                                 time = t1;
00262                         }
00263                         // take care of the last force points
00264                         // find the starting and stopping points
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                         // now run to the end, unless we're already there
00287                         if(finaltime > time)
00288                         {
00289                                 reactionMover->Move(time,finaltime,experiment->GetReactionNetwork());
00290                         }
00291                         time = finaltime;
00292                 }
00293         }
00294 
00295         // compute the scaling factors; successive computation of the cost is somewhat of a
00296         // repeat of this operation, so combining the two calculations in the future may 
00297         // speed things up
00298         for(int chemCount = 0; chemCount < nChem; chemCount++)
00299         {
00300                 double num = 0.0;
00301                 double den = 0.0;
00302                 //ChemicalTimeSeriesData *cTSD = experiment->GetChemicalTimeSeriesData();
00303                 // do we need a scaling factor AND is it allowed to vary?
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                                 // loop over the number of data points and compute
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 //              cerr << "chemCount " << chemCount << " has B factor " << factor << endl;
00338                 m_pConversionFactors[chemCount]->SetFactorValue(factor);
00339 
00340         }
00341                 
00342         // compute the cost
00343         int nextAvailableIndex = 0;
00344         if(timeSeriesWeight > 0.0)
00345         {
00346                 // load the time series residuals in the master list
00347                 for(int chemIt = 0; chemIt < (experiment->GetChemicalTimeSeriesData()->GetTimeVector()).size(); chemIt++)
00348                 {
00349                         // get the keys
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                         // use the keys
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                         // rescale sim if there is a conversion factor
00358                         if(m_pConversionFactors[targetChem]->IsFactorNeeded())
00359                         {
00360                                 sim = sim*(m_pConversionFactors[targetChem]->GetFactorValue());
00361                         }
00362                         // compute quantities of interest
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                 // load in the rate constants residuals in the master list
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                 // load in the initial chemical concentration residuals in the master list
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         // notify observers
00406         Notify();
00407 //      cerr << "Leaving SNM::CompRes\n";
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 }

Generated on Mon Nov 3 09:38:15 2003 by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002