---

SingleNetworkMinimizable.cpp

Go to the documentation of this file.

// SingleNetworkMinimizable.cpp: implementation of the SingleNetworkMinimizable class.
//

#include "SingleNetworkMinimizable.h"

// Construction/Destruction

SingleNetworkMinimizable::SingleNetworkMinimizable()
{

}

SingleNetworkMinimizable::SingleNetworkMinimizable(Experiment *experiment,CReactionMover *reactionMover, bool rateConstantsOptimizable, bool initialChemConcOptimizable, bool logsInObjectiveFunction, double timeSeriesWeight, double rateConstantsWeight, double initialChemConcWeight, int nCells)
{
        Initialize(experiment,reactionMover,nCells,rateConstantsOptimizable,initialChemConcOptimizable,logsInObjectiveFunction,timeSeriesWeight,rateConstantsWeight,initialChemConcWeight);
}

SingleNetworkMinimizable::~SingleNetworkMinimizable()
{
        delete [] initialNetworkData;
        delete [] currentNetworkData;
        this->reactionMover->Detach(cellObserver);
        delete cellObserver;
        for(int i = 0; i < m_pConversionFactors.size(); i++)
        {
                delete m_pConversionFactors[i];
        }
}

/*Experiment *SingleNetworkMinimizable::GetExperiment()
{
        return this->experiment;
}
*/
void SingleNetworkMinimizable::Initialize(Experiment *experiment, CReactionMover *reactionMover, int nCells, bool rateConstantsOptimizable, bool initialChemConcOptimizable, bool logsInObjectiveFunction, double timeSeriesWeight, double rateConstantsWeight, double initialChemConcWeight)
{
        this->experiment = experiment;
        this->reactionMover = reactionMover;
        this->nCells = nCells;
        this->rateConstantsOptimizable = rateConstantsOptimizable;
        this->initialChemConcOptimizable = initialChemConcOptimizable;
        this->logsInObjectiveFunction = logsInObjectiveFunction;

        this->SetTimeSeriesWeight(timeSeriesWeight);
        this->SetRateConstantsWeight(rateConstantsWeight);
        this->SetInitialChemConcWeight(initialChemConcWeight);

        // Read in the data
        this->experiment->ReadData();
        // Read in forcing data (if any)
        this->experiment->ReadForcingData();
        // Initialize variables dependent on the Experiment class
        int nChem = this->experiment->GetReactionNetwork()->GetNumberOfChemicals();
        int nRateConstants = this->experiment->GetReactionNetwork()->GetNumberOfRateConstants();
        int intTime = ((int) this->experiment->GetChemicalTimeSeriesData()->GetLatestDataTime()) + 2;
        cellObserver = new CellAverageObserver(nChem,intTime,nCells);

        // attach the observer
        this->reactionMover->Attach(cellObserver);

        initialNetworkData = new double[nRateConstants + nChem];
        currentNetworkData = new double[nRateConstants + nChem];

        // load up the vector of conversion factors with dummies assuming all data is 
        // expressed in terms of number of molecules; if we don't have data we don't 
        // need conversion factors
        for(int l = 0; l < nChem; l++)
        {
                if((experiment->GetChemicalTimeSeriesData()->GetTimeSeries(l)).size() > 0 )
                {
                        m_pConversionFactors.push_back(new CConversionFactor(true,false,1.0));
                }
                else
                {
                        m_pConversionFactors.push_back(new CConversionFactor(false));
                }
        }

        m_iNBFactors = 0;

        // store beginning network data - this allows us to restore the 
        // initial state of the network to recover from the beginning
        for(int i = 0; i < nRateConstants; i++)
        {
                initialNetworkData[i] = this->experiment->GetReactionNetwork()->GetRateConstant(i)->GetInitialValue();
                currentNetworkData[i] = initialNetworkData[i];
        }
        for(int i = nRateConstants; i < nRateConstants + nChem; i++)
        {
                initialNetworkData[i] = this->experiment->GetReactionNetwork()->GetChemical(i-nRateConstants)->GetInitialAmount();
                currentNetworkData[i] = initialNetworkData[i];
        }

        int localResidualsSize = 0;
        if(this->timeSeriesWeight > 0.0) {localResidualsSize += this->experiment->GetChemicalTimeSeriesData()->GetNDataPoints();}
        if(this->rateConstantsWeight > 0.0) {localResidualsSize += nRateConstants;}
        if(this->initialChemConcWeight > 0.0) {localResidualsSize       += nChem;}

        Allocate(localResidualsSize);

}

double SingleNetworkMinimizable::ComputeES(double *parameters, double T)
{
        m_dEnergyLastComputed = this->ComputeResiduals(parameters);
        double S = this->EntropyShift(T) - 0.5*m_dEntropyLastComputed;
        m_dEntropyLastComputed = S;
        return m_dEnergyLastComputed;
}
/*
double SingleNetworkMinimizable::EntropyShift(double T)
{
        return 0.5*m_iNBFactors*log(2*3.14159265358979323846*T);
}
*/
double SingleNetworkMinimizable::F(double *parameters, double T)
{
        m_dEnergyLastComputed = this->ComputeResiduals(parameters);
        double S = -0.5*m_dEntropyLastComputed;
        m_dEntropyLastComputed = S;
        return (m_dEnergyLastComputed - T*m_dEntropyLastComputed);
}

double SingleNetworkMinimizable::F0(double *parameters, double T)
{
        m_dEnergyLastComputed = this->ComputeResiduals(parameters);
        double S = this->EntropyShift(T) - 0.5*m_dEntropyLastComputed;
        m_dEntropyLastComputed = S;
        return (m_dEnergyLastComputed - T*m_dEntropyLastComputed);
}

double SingleNetworkMinimizable::ComputeResiduals(double *parameters)
{
        m_iNBFactors = 0;
        m_dEntropyLastComputed = 0.0;
        double partialEntropy = 0.0;
        
        double objFuncVal = 0.0;
        int nRateConstants = experiment->GetReactionNetwork()->GetNumberOfRateConstants();
        int nChem = experiment->GetReactionNetwork()->GetNumberOfChemicals();
        
        // put the parameters into the network
        // rate constants first, then initial chemical concentrations
        if(rateConstantsOptimizable && initialChemConcOptimizable)
        {
                for(int i = 0; i < nRateConstants; i++)
                {
                        experiment->GetReactionNetwork()->GetRateConstant(i)->SetRateConstant(parameters[i]);
                        currentNetworkData[i] = parameters[i];
                }
                for(int i = nRateConstants; i < nRateConstants + nChem; i++)
                {
                        experiment->GetReactionNetwork()->GetChemical(i - nRateConstants)->SetAmount(parameters[i]);
                        currentNetworkData[i] = parameters[i];
                }
        }
        else if(rateConstantsOptimizable)
        {
                for(int i = 0; i < nRateConstants; i++)
                {
                        experiment->GetReactionNetwork()->GetRateConstant(i)->SetRateConstant(parameters[i]);
                        currentNetworkData[i] = parameters[i];
                }
                experiment->GetReactionNetwork()->ChemicalReset();
        }
        else if(initialChemConcOptimizable)
        {
                for(int i = 0; i < nChem; i++)
                {
                        experiment->GetReactionNetwork()->GetChemical(i)->SetAmount(parameters[i]);
                        currentNetworkData[i+nRateConstants] = parameters[i];
                }
                experiment->GetReactionNetwork()->RateConstantReset();
        }
        else
        {
                // do one run at the existing parameter values
                // this ensures that multiple calls with the same parameters
                // give the same cost, but it also ensures that only time
                // series contribute to the cost - this is probably not what
                // we want?
                experiment->GetReactionNetwork()->RateConstantReset();
                experiment->GetReactionNetwork()->ChemicalReset();
        }

        // zero the observer
        cellObserver->ZeroConcentrations();

        for(int cellCounter = 1; cellCounter <= nCells; cellCounter++)
        {
                // make sure the concentration in the network is re-initialized, which
                // is done differently whether you allow initial conc. optimization
                // or not
                //experiment->GetReactionNetwork()->ChemicalReset();
                if(initialChemConcOptimizable)
                {
                        for(int i = nRateConstants; i < nRateConstants + nChem; i++)
                        {
                                experiment->GetReactionNetwork()->GetChemical(i - nRateConstants)->SetAmount(parameters[i]);
                        }
                }
                else
                {
                        experiment->GetReactionNetwork()->ChemicalReset();
                }
                // if there's no forcing data at all, we can just run from (0,tfinal). 
                // Otherwise we have to stop running in between forcing changes.  This
                // block might be able to be rearranged to be much faster. KSB 4/16/02

                // reset step size statistics in either case
                reactionMover->ResetStepStatistics(experiment->GetChemicalTimeSeriesData()->GetLatestDataTime());
                
                // now move through the rest
                if(experiment->GetForcingData()->GetForcingVector().size() == 0)
                {
                        // take data at time 0
                        reactionMover->Move(0.0,0.0,experiment->GetReactionNetwork());
                        reactionMover->Move(0.0,experiment->GetChemicalTimeSeriesData()->GetLatestDataTime(),experiment->GetReactionNetwork());
                }
                else
                {
                        int forceCount = 0;
                        double time = 0.0;
                        int lastForce = experiment->GetForcingData()->GetForcingVector().size();
                        double finaltime = experiment->GetChemicalTimeSeriesData()->GetLatestDataTime();
                        double finalForceTime = experiment->GetForcingData()->GetLatestForcingTime();
                        int index;
                        int start,stop;
                        double t1,t2,value;
                        while(time < finalForceTime)
                        {
                                start = forceCount;
                                stop = forceCount;
                                t1 = experiment->GetForcingData()->GetForcingVector()[start]->GetTime();

                                // find the starting and stopping points
                                for(int i = start; i < lastForce; i++)
                                {
                                        t2 = experiment->GetForcingData()->GetForcingVector()[i]->GetTime();
                                        if(t1 == t2)
                                        {
                                        }
                                        else
                                        {
                                                stop = i-1;
                                                for(int j = start; j <= stop; j++)
                                                {
                                                        index = experiment->GetForcingData()->GetForcingVector()[j]->GetChemNumber();
                                                        value = experiment->GetForcingData()->GetForcingVector()[j]->GetValue();
                                                        experiment->GetReactionNetwork()->GetChemical(index)->SetAmount(value);
                                                }
                                                break;
                                        }
                                }
                                forceCount = stop+1;
                                // move for the correct amount of time
                                reactionMover->Move(time,t1,experiment->GetReactionNetwork());
                                time = t1;
                        }
                        // take care of the last force points
                        // find the starting and stopping points
                        t2 = finalForceTime;
                        stop = lastForce-1;
                        for(int i = lastForce-1; i > 0; i--)
                        {
                                        t1 = experiment->GetForcingData()->GetForcingVector()[i]->GetTime();
                                        if(t1 == t2)
                                        {
                                        }
                                        else
                                        {
                                                start = i+1;
                                                for(int j = start; j <= stop; j++)
                                                {
                                                        index = experiment->GetForcingData()->GetForcingVector()[j]->GetChemNumber();
                                                        value = experiment->GetForcingData()->GetForcingVector()[j]->GetValue();
                                                        experiment->GetReactionNetwork()->GetChemical(index)->SetAmount(value);
                                                }
                                                break;
                                        }
                        }
                        
                        // now run to the end, unless we're already there
                        if(finaltime > time)
                        {
                                reactionMover->Move(time,finaltime,experiment->GetReactionNetwork());
                        }
                        time = finaltime;
                }
        }

        // compute the scaling factors; successive computation of the cost is somewhat of a
        // repeat of this operation, so combining the two calculations in the future may 
        // speed things up
        for(int chemCount = 0; chemCount < nChem; chemCount++)
        {
                double num = 0.0;
                double den = 0.0;
                //ChemicalTimeSeriesData *cTSD = experiment->GetChemicalTimeSeriesData();
                // do we need a scaling factor AND is it allowed to vary?
                if(m_pConversionFactors[chemCount]->IsFactorNeeded())
                {
                        if(m_pConversionFactors[chemCount]->IsFactorFixed())
                        {
                                num = m_pConversionFactors[chemCount]->GetFactorValue();
                                den = 1.0;
                        }
                        else
                        {
                                m_iNBFactors++;
                                int dataSize = (experiment->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)).size();
                                // loop over the number of data points and compute
                                for(int nData = 0; nData < dataSize; nData++)
                                {
                                        double time = experiment->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetTime();
                                        double data = experiment->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetDatum();
                                        double error = experiment->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetError();
                                        double sim = (cellObserver->GetAverageConcentration())[chemCount][int(time)];
                                        num += ((2*sim)/(error*error))*data;
                                        den += ((2*sim)/(error*error))*sim;

                                        partialEntropy += ((sim*sim)/(error*error));
                                }
                                m_dEntropyLastComputed += log(partialEntropy);
                        }
                }
                else
                {
                        num = 1.0;
                        den = 1.0;
                }

                double factor = num/den;
//              cerr << "chemCount " << chemCount << " has B factor " << factor << endl;
                m_pConversionFactors[chemCount]->SetFactorValue(factor);

        }
                
        // compute the cost
        int nextAvailableIndex = 0;
        if(timeSeriesWeight > 0.0)
        {
                // load the time series residuals in the master list
                for(int chemIt = 0; chemIt < (experiment->GetChemicalTimeSeriesData()->GetTimeVector()).size(); chemIt++)
                {
                        // get the keys
                        double targetTime = (experiment->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetTime();
                        int targetChem = (experiment->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetChemNumber();
                        int vectorIndex = (experiment->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetVectorIndex();
                        // use the keys
                        double data = (experiment->GetChemicalTimeSeriesData()->GetTimeSeries(targetChem))[vectorIndex]->GetDatum();
                        double sigma = (experiment->GetChemicalTimeSeriesData()->GetTimeSeries(targetChem))[vectorIndex]->GetError();
                        double sim = (cellObserver->GetAverageConcentration())[targetChem][int(targetTime)];
                        // rescale sim if there is a conversion factor
                        if(m_pConversionFactors[targetChem]->IsFactorNeeded())
                        {
                                sim = sim*(m_pConversionFactors[targetChem]->GetFactorValue());
                        }
                        // compute quantities of interest
                        residuals[nextAvailableIndex] = (sqrt(timeSeriesWeight/2.0))*((data - sim)/sigma);
                        objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
                        nextAvailableIndex++;
                }
        }
        if(rateConstantsWeight > 0.0)
        {
                // load in the rate constants residuals in the master list
                if(logsInObjectiveFunction)
                {
                        for(int i = 0; i < nRateConstants; i++)
                        {
                                double sigma = experiment->GetReactionNetwork()->GetRateConstant(i)->GetErrorInInitialValue();
                                double sigmalog = sigma/initialNetworkData[i];
                                residuals[nextAvailableIndex] = (sqrt(rateConstantsWeight/2.0))*(log(initialNetworkData[i]/currentNetworkData[i])/sigmalog);
                                objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
                                nextAvailableIndex++;
                        }
                }
                else
                {
                        for(int i = 0; i < nRateConstants; i++)
                        {
                                double sigma = experiment->GetReactionNetwork()->GetRateConstant(i)->GetErrorInInitialValue();
                                residuals[nextAvailableIndex] = (sqrt(rateConstantsWeight/2.0))*((initialNetworkData[i] - currentNetworkData[i])/sigma);
                                objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
                                nextAvailableIndex++;
                        }
                }
        }
        if(initialChemConcWeight > 0.0)
        {
                // load in the initial chemical concentration residuals in the master list
                for(int i = 0; i < nChem; i++)
                {       
                        double sigma = experiment->GetReactionNetwork()->GetChemical(i)->GetErrorInInitialAmount();
                        residuals[nextAvailableIndex] = (sqrt(initialChemConcWeight/2.0))*((initialNetworkData[i+nRateConstants] - currentNetworkData[i+nRateConstants])/sigma);
                        objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
                        nextAvailableIndex++;
                }
        }
        
        // notify observers
        Notify();
//      cerr << "Leaving SNM::CompRes\n";
        return objFuncVal;
}

int SingleNetworkMinimizable::GetNParameters()
{
        int nRateConstants = experiment->GetReactionNetwork()->GetNumberOfRateConstants();
        int nChem = experiment->GetReactionNetwork()->GetNumberOfChemicals();
        int nParameters;

        if(rateConstantsOptimizable && initialChemConcOptimizable)
        {
                nParameters = nRateConstants + nChem;
        }
        else if(rateConstantsOptimizable)
        {
                nParameters = nRateConstants;
        }
        else if(initialChemConcOptimizable)
        {
                nParameters = nChem;
        }
        else
        {
                nParameters = 0;
        }

        return nParameters;
}

double SingleNetworkMinimizable::GetParameter(int parIndex)
{
        int nP = GetNParameters();
        
        if(parIndex > nP)
        {
                throw std::runtime_error("Parameter index is off parameter list.");
        }
        else
        {
                return currentNetworkData[parIndex];
        }

        return 0.0;
}

---