---

SharedBFactorNetworkMinimizable.cpp

Go to the documentation of this file.

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

#include "SharedBFactorNetworkMinimizable.h"

// Construction/Destruction

SharedBFactorNetworkMinimizable::SharedBFactorNetworkMinimizable()
{

}

SharedBFactorNetworkMinimizable::SharedBFactorNetworkMinimizable(std::vector<Experiment *> experiments,std::vector<CReactionMover *> reactionMovers, bool rateConstantsOptimizable, bool initialChemConcOptimizable, bool logsInObjectiveFunction, double timeSeriesWeight, double rateConstantsWeight, double initialChemConcWeight, int nCells)
{
        Initialize(experiments,reactionMovers,nCells,rateConstantsOptimizable,initialChemConcOptimizable,logsInObjectiveFunction,timeSeriesWeight,rateConstantsWeight,initialChemConcWeight);

}

SharedBFactorNetworkMinimizable::~SharedBFactorNetworkMinimizable()
{
//      delete [] initialNetworkData;
//      delete [] currentNetworkData;
        for (int expCounter=0; expCounter < nExp; expCounter++) {
                this->reactionMovers[expCounter]->Detach(cellObservers[expCounter]);
                delete cellObservers[expCounter];
        }
        for(int i = 0; i < m_pConversionFactors.size(); i++)
        {
                delete m_pConversionFactors[i];
        }
}

void SharedBFactorNetworkMinimizable::Initialize(std::vector<Experiment *> experiments, std::vector<CReactionMover *> reactionMovers, int nCells, bool rateConstantsOptimizable, bool initialChemConcOptimizable, bool logsInObjectiveFunction, double timeSeriesWeight, double rateConstantsWeight, double initialChemConcWeight)
{
        this->nExp = experiments.size();
        this->experiments = experiments;
        this->reactionMovers = reactionMovers;
        this->nCells = nCells;
        this->rateConstantsOptimizable = rateConstantsOptimizable;
        this->initialChemConcOptimizable = initialChemConcOptimizable;
        this->logsInObjectiveFunction = logsInObjectiveFunction;

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

        for (int i=0; i < nExp; i++) {
                // Read in the data
                this->experiments[i]->ReadData();
                // Read in forcing data (if any)
                this->experiments[i]->ReadForcingData();
        }

        // Initialize variables dependent on the Experiment class
        // Should be safe using nChem and nRateConstants from the first experiment because I don't think this can be changed. -JJW 9-18-03
        int nChem = this->experiments[0]->GetReactionNetwork()->GetNumberOfChemicals();
        int nRateConstants = this->experiments[0]->GetReactionNetwork()->GetNumberOfRateConstants();
//      int intTime = ((int) this->experiments[0]->GetChemicalTimeSeriesData()->GetLatestDataTime()) + 2;

        for(int expCounter=0; expCounter < nExp; expCounter++) {
                assert(nChem == this->experiments[expCounter]->GetReactionNetwork()->GetNumberOfChemicals());
                assert(nRateConstants == this->experiments[expCounter]->GetReactionNetwork()->GetNumberOfRateConstants() );
                int intTime = ((int) this->experiments[expCounter]->GetChemicalTimeSeriesData()->GetLatestDataTime()) + 2;

                this->cellObservers.push_back(new CellAverageObserver(nChem,intTime,nCells));
                // attach the observer
                this->reactionMovers[expCounter]->Attach(this->cellObservers[expCounter]);
        }

        initialNetworkData.reserve(nExp);
        currentNetworkData.reserve(nExp);

        for (int expCounter=0; expCounter<nExp; expCounter++) {
                currentNetworkData.push_back(new std::vector<double *>);
                (*currentNetworkData[expCounter]).reserve(nRateConstants + nChem);
                initialNetworkData.push_back(new std::vector<double *>);
                (*initialNetworkData[expCounter]).reserve(nRateConstants + nChem);
        }

//      initialNetworkData = new double[nExp*(nRateConstants + nChem)];
//      currentNetworkData = new double[nExp*(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
        // need one conversion factor for every chemical which has any data on it in any
        // of the experiments but need same conversionfactor for same chemical in 
        // different experiments
        for(int l = 0; l < nChem; l++)
        {
                bool needFactorI = false;
                for (int i = 0; i < nExp; i++) {
                        if((experiments[i]->GetChemicalTimeSeriesData()->GetTimeSeries(l)).size() > 0 )
                        {
                                needFactorI = true;
                                break;
                        }
                }
                if (needFactorI) {
                        m_pConversionFactors.push_back(new CConversionFactor(true,false,1.0));
                }
                else
                {
                        m_pConversionFactors.push_back(new CConversionFactor(false));
                }
        }
        // I don't know what this is for -JJW 9-12-03
        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 expCounter=0; expCounter < nExp; expCounter++) {
                for(int i = 0; i < nRateConstants; i++)
                {
                        (*initialNetworkData[expCounter]).push_back(new double);
                        (*(*initialNetworkData[expCounter])[i]) = this->experiments[expCounter]->GetReactionNetwork()->GetRateConstant(i)->GetInitialValue();
                        (*currentNetworkData[expCounter]).push_back(new double);
                        (*(*currentNetworkData[expCounter])[i]) = (*(*initialNetworkData[expCounter])[i]);
                }
                for(int i = nRateConstants; i < nRateConstants + nChem; i++)
                {
                        (*initialNetworkData[expCounter]).push_back(new double);
                        (*(*initialNetworkData[expCounter])[i]) = this->experiments[expCounter]->GetReactionNetwork()->GetChemical(i-nRateConstants)->GetInitialAmount();
                        (*currentNetworkData[expCounter]).push_back(new double);
                        (*(*currentNetworkData[expCounter])[i]) = (*(*initialNetworkData[expCounter])[i]);

                }
        }

        int localResidualsSize = 0;
        if(this->timeSeriesWeight > 0.0) {
                for (int expCounter=0; expCounter < nExp; expCounter++) {
                        localResidualsSize += this->experiments[expCounter]->GetChemicalTimeSeriesData()->GetNDataPoints();
                }
        }
        if(this->rateConstantsWeight > 0.0) {localResidualsSize += nExp*nRateConstants;}
        if(this->initialChemConcWeight > 0.0) {localResidualsSize       += nExp*nChem;}

        Allocate(localResidualsSize);

}

double SharedBFactorNetworkMinimizable::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 SharedBFactorNetworkMinimizable::EntropyShift(double T)
{
        return 0.5*m_iNBFactors*log(2*3.14159265358979323846*T);
}
*/
double SharedBFactorNetworkMinimizable::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 SharedBFactorNetworkMinimizable::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 SharedBFactorNetworkMinimizable::ComputeResiduals(double *parameters)
{
        m_iNBFactors = 0;
        m_dEntropyLastComputed = 0.0;
        double partialEntropy = 0.0;
        
        double objFuncVal = 0.0;
        //should be safe since all the experiments should have the same number of rateconstants and chemicals (maybe) -JJW
        int nRateConstants = experiments[0]->GetReactionNetwork()->GetNumberOfRateConstants();
        int nChem = experiments[0]->GetReactionNetwork()->GetNumberOfChemicals();
        // put the parameters into the network
        // rate constants first, then initial chemical concentrations
//              int base = j*(nRateConstants+nChem);
//              for(int i = base+0; i < base+nRateConstants; i++)
//                      initialNetworkData[i] = this->experiments[j]->GetReactionNetwork()->GetRateConstant(i-base)->GetInitialValue();
//              for(int i = base+nRateConstants; i < base+nRateConstants + nChem; i++)
//                      initialNetworkData[i] = this->experiments[j]->GetReactionNetwork()->GetChemical(i-base-nRateConstants)->GetInitialAmount();

        if(rateConstantsOptimizable && initialChemConcOptimizable)
        {
                for (int expCounter = 0; expCounter < nExp; expCounter++) {
                        for(int i = 0; i < nRateConstants; i++)
                        {
                                experiments[expCounter]->GetReactionNetwork()->GetRateConstant(i)->SetRateConstant(parameters[i]);
                                (*(*currentNetworkData[expCounter])[i]) = parameters[i];
                        }
                        for(int i = nRateConstants; i < nRateConstants + nChem; i++)
                        {
                                experiments[expCounter]->GetReactionNetwork()->GetChemical(i - nRateConstants)->SetAmount(parameters[i]);
                                (*(*currentNetworkData[expCounter])[i]) = parameters[i];
                        }
                }
        }
        else if(rateConstantsOptimizable)
        {
                for (int expCounter = 0; expCounter < nExp; expCounter++) {
                        for(int i = 0; i < nRateConstants; i++)
                        {
                                experiments[expCounter]->GetReactionNetwork()->GetRateConstant(i)->SetRateConstant(parameters[i]);
                                (*(*currentNetworkData[expCounter])[i]) = parameters[i];
                        }
                        experiments[expCounter]->GetReactionNetwork()->ChemicalReset();
                }
        }
        else if(initialChemConcOptimizable)
        {
                for (int expCounter = 0; expCounter < nExp; expCounter++) {
                        for(int i = 0; i < nChem; i++)
                        {
                                experiments[expCounter]->GetReactionNetwork()->GetChemical(i)->SetAmount(parameters[i]);
                                (*(*currentNetworkData[expCounter])[i+nRateConstants]) = parameters[i];
                        }
                        experiments[expCounter]->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?
                for (int expCounter = 0; expCounter < nExp; expCounter++) {
                        experiments[expCounter]->GetReactionNetwork()->RateConstantReset();
                        experiments[expCounter]->GetReactionNetwork()->ChemicalReset();
                }
        }

        // zero the observer
        for (int expCounter = 0; expCounter < nExp; expCounter++) {
                cellObservers[expCounter]->ZeroConcentrations();
        }
        for (int expCounter = 0; expCounter < nExp; expCounter++) {

                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++)
                                {
                                        experiments[expCounter]->GetReactionNetwork()->GetChemical(i - nRateConstants)->SetAmount(parameters[i]);
                                }
                        }
                        else
                        {
                                experiments[expCounter]->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
                        reactionMovers[expCounter]->ResetStepStatistics(experiments[expCounter]->GetChemicalTimeSeriesData()->GetLatestDataTime());

                        // now move through the rest
                        if(experiments[expCounter]->GetForcingData()->GetForcingVector().size() == 0)
                        {
                                // take data at time 0
                                reactionMovers[expCounter]->Move(0.0,0.0,experiments[expCounter]->GetReactionNetwork());
                                reactionMovers[expCounter]->Move(0.0,experiments[expCounter]->GetChemicalTimeSeriesData()->GetLatestDataTime(),experiments[expCounter]->GetReactionNetwork());
                        }
                        else
                        {
                                int forceCount = 0;
                                double time = 0.0;
                                int lastForce = experiments[expCounter]->GetForcingData()->GetForcingVector().size();
                                double finaltime = experiments[expCounter]->GetChemicalTimeSeriesData()->GetLatestDataTime();
                                double finalForceTime = experiments[expCounter]->GetForcingData()->GetLatestForcingTime();
                                int index;
                                int start,stop;
                                double t1,t2,value;
                                while(time < finalForceTime)
                                {
                                        start = forceCount;
                                        stop = forceCount;
                                        t1 = experiments[expCounter]->GetForcingData()->GetForcingVector()[start]->GetTime();
                                        // find the starting and stopping points
                                        for(int i = start; i < lastForce; i++)
                                        {
                                                t2 = experiments[expCounter]->GetForcingData()->GetForcingVector()[i]->GetTime();
                                                if(t1 == t2)
                                                {
                                                }
                                                else
                                                {
                                                        stop = i-1;
                                                        for(int k = start; k <= stop; k++)
                                                        {
                                                                index = experiments[expCounter]->GetForcingData()->GetForcingVector()[k]->GetChemNumber();
                                                                value = experiments[expCounter]->GetForcingData()->GetForcingVector()[k]->GetValue();
                                                                experiments[expCounter]->GetReactionNetwork()->GetChemical(index)->SetAmount(value);
                                                        }
                                                        break;
                                                }
                                        }
                                        forceCount = stop+1;
                                        // move for the correct amount of time
                                        reactionMovers[expCounter]->Move(time,t1,experiments[expCounter]->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 = experiments[expCounter]->GetForcingData()->GetForcingVector()[i]->GetTime();
                                        if(t1 == t2)
                                        {
                                        }
                                        else
                                        {
                                                start = i+1;
                                                for(int k = start; k <= stop; k++)
                                                {
                                                        index = experiments[expCounter]->GetForcingData()->GetForcingVector()[k]->GetChemNumber();
                                                        value = experiments[expCounter]->GetForcingData()->GetForcingVector()[k]->GetValue();
                                                        experiments[expCounter]->GetReactionNetwork()->GetChemical(index)->SetAmount(value);
                                                }
                                                break;
                                        }
                                }
                                // now run to the end, unless we're already there
                                if(finaltime > time)
                                {
                                        reactionMovers[expCounter]->Move(time,finaltime,experiments[expCounter]->GetReactionNetwork());
                                }
                                time = finaltime;
                        }
                } //end cellCount
        } //end nExp

        // 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++;
                                for (int expCounter=0; expCounter < nExp; expCounter++) {
                                        int dataSize = (experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)).size();
                                        // loop over the number of data points and compute
                                        for(int nData = 0; nData < dataSize; nData++)
                                        {
                                                double time = experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetTime();
                                                double data = experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetDatum();
                                                double error = experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeSeries(chemCount)[nData]->GetError();
                                                double sim = (cellObservers[expCounter]->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);
                                } // end expCounter
                        }
                }
                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)
        {
                for (int expCounter=0; expCounter<nExp; expCounter++) {
                        // load the time series residuals in the master list
                        for(int chemIt = 0; chemIt < (experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeVector()).size(); chemIt++)
                        {
                                // get the keys
                                double targetTime = (experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetTime();
                                int targetChem = (experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetChemNumber();
                                int vectorIndex = (experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeVector())[chemIt]->GetVectorIndex();
                                // use the keys
                                double data = (experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeSeries(targetChem))[vectorIndex]->GetDatum();
                                double sigma = (experiments[expCounter]->GetChemicalTimeSeriesData()->GetTimeSeries(targetChem))[vectorIndex]->GetError();
                                double sim = (cellObservers[expCounter]->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++;
                        }
                } // end expCounter
        }

        if(rateConstantsWeight > 0.0)
        {
                for (int expCounter=0; expCounter<nExp; expCounter++) {
                        // load in the rate constants residuals in the master list
                        if(logsInObjectiveFunction)
                        {
                                for(int i = 0; i < nRateConstants; i++)
                                {
                                        double sigma = experiments[expCounter]->GetReactionNetwork()->GetRateConstant(i)->GetErrorInInitialValue();
                                        double sigmalog = sigma/(*(*initialNetworkData[expCounter])[i]);
                                        residuals[nextAvailableIndex] = (sqrt(rateConstantsWeight/2.0))*(log((*(*initialNetworkData[expCounter])[i])/(*(*currentNetworkData[expCounter])[i]))/sigmalog);
                                        objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
                                        nextAvailableIndex++;
                                }
                        }
                        else
                        {
                                for(int i = 0; i < nRateConstants; i++)
                                {
                                        double sigma = experiments[expCounter]->GetReactionNetwork()->GetRateConstant(i)->GetErrorInInitialValue();
                                        residuals[nextAvailableIndex] = (sqrt(rateConstantsWeight/2.0))*((initialNetworkData[i+expCounter*(nRateConstants+nChem)] - currentNetworkData[i+expCounter*(nRateConstants+nChem)])/sigma);
                                        objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
                                        nextAvailableIndex++;
                                }
                        }
                } // end expCounter
        }

        if(initialChemConcWeight > 0.0)
        {
                for (int expCounter=0; expCounter<nExp; expCounter++) {
                        // load in the initial chemical concentration residuals in the master list
                        for(int i = 0; i < nChem; i++)
                        {       
                                double sigma = experiments[expCounter]->GetReactionNetwork()->GetChemical(i)->GetErrorInInitialAmount();
                                residuals[nextAvailableIndex] = (sqrt(initialChemConcWeight/2.0))*((initialNetworkData[i+nRateConstants+expCounter*(nRateConstants+nChem)] - currentNetworkData[i+nRateConstants+expCounter*(nRateConstants+nChem)])/sigma);
                                objFuncVal += residuals[nextAvailableIndex]*residuals[nextAvailableIndex];
                                nextAvailableIndex++;
                        }
                } //end expCounter
        }

        // notify observers
        Notify();

        return objFuncVal;
}

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

        int nRateConstants = experiments[0]->GetReactionNetwork()->GetNumberOfRateConstants();
        int nChem = experiments[0]->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 SharedBFactorNetworkMinimizable::GetParameter(int parIndex)
{
        int nP = GetNParameters();
        
        if(parIndex > nP)
        {
                throw std::runtime_error("Parameter index is off parameter list.");
        }
        else
        {
                return (*(*currentNetworkData[0])[parIndex]);
        }

        return 0.0;
}

---