---

EnsembleCombinationDirector.cpp

Go to the documentation of this file.

// EnsembleCombinationDirector.cpp: implementation of the CEnsembleCombinationDirector class.
//

#include "EnsembleCombinationDirector.h"

// Construction/Destruction

CEnsembleCombinationDirector::CEnsembleCombinationDirector(int ensembleSize)
{
        m_iEnsembleSize = ensembleSize;
}

CEnsembleCombinationDirector::~CEnsembleCombinationDirector()
{
        int i;
        // lots to delete
        for(i = 0; i < m_vpdEnsembleParameters.size(); i++)
        {
                delete [] m_vpdEnsembleParameters[i];
        }
        for(i = 0; i < m_vpdNoDataSumC.size(); i++)
        {
                delete [] m_vpdNoDataSumC[i][0];
                delete [] m_vpdNoDataSumC[i];
                delete [] m_vpdNoDataSumCSqr[i][0];
                delete [] m_vpdNoDataSumCSqr[i];
        }
        for(i = 0; i < m_vpdDataSumC.size(); i++)
        {
                delete [] m_vpdDataSumC[i][0];
                delete [] m_vpdDataSumC[i];
                delete [] m_vpdDataSumCSqr[i][0];
                delete [] m_vpdDataSumCSqr[i];
        }
        for(i = 0; i < m_vpEnsembleDataPlotters.size(); i++)
        {
                delete m_vpEnsembleDataPlotters[i];
        }
        for(i = 0; i < m_vpEnsembleRunPlotters.size(); i++)
        {
                delete m_vpEnsembleRunPlotters[i];
        }
        for(i = 0; i < m_vpMinimizables.size(); i++)
        {
                delete m_vpMinimizables[i];
        }
        for(i = 0; i < m_vpNetworks.size(); i++)
        {
                delete m_vpNetworks[i];
        }
        for(i = 0; i < m_vpExperiments.size(); i++)
        {
                delete m_vpExperiments[i];
        }
        for(i = 0; i < m_vpMovers.size(); i++)
        {
                delete m_vpMovers[i];
        }

}

void CEnsembleCombinationDirector::SumAllocate()
{
        int i;
        // allocate for the runnables
        for(i = 0; i < m_vpRunnables.size(); i++)
        {
                int nSteps = m_vpRunnables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                int nChem = m_vpRunnables[i]->GetReactionNetwork()->GetNumberOfChemicals();
                m_vpdNoDataSumC.push_back(new double *[nChem]);
                m_vpdNoDataSumCSqr.push_back(new double *[nChem]);
                m_vpdNoDataSumC[i][0] = new double[nChem*nSteps];
                m_vpdNoDataSumCSqr[i][0] = new double[nChem*nSteps];
                for(int chem = 1; chem < nChem; chem++)
                {
                        m_vpdNoDataSumC[i][chem] = &(m_vpdNoDataSumC[i][0][chem*nSteps]);
                        m_vpdNoDataSumCSqr[i][chem] = &(m_vpdNoDataSumCSqr[i][0][chem*nSteps]);
                }
                // zero the array just created
                for(int l = 0; l < nChem; l++)
                {
                        for(int m = 0; m < nSteps; m++)
                        {
                                m_vpdNoDataSumC[i][l][m] = 0;
                                m_vpdNoDataSumCSqr[i][l][m] = 0;
                        }
                }
        }
        // do the same thing for the minimizables
        for(i = 0; i < m_vpMinimizables.size(); i++)
        {
                int nSteps = m_vpMinimizables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                int nChem = m_vpMinimizables[i]->GetExperiment()->GetChemicalTimeSeriesData()->GetNChemicals();
                m_vpdDataSumC.push_back(new double *[nChem]);
                m_vpdDataSumCSqr.push_back(new double *[nChem]);
                m_vpdDataSumC[i][0] = new double[nChem*nSteps];
                m_vpdDataSumCSqr[i][0] = new double[nChem*nSteps];
                for(int chem = 1; chem < nChem; chem++)
                {
                        m_vpdDataSumC[i][chem] = &(m_vpdDataSumC[i][0][chem*nSteps]);
                        m_vpdDataSumCSqr[i][chem] = &(m_vpdDataSumCSqr[i][0][chem*nSteps]);
                }
                // zero the array just created
                for(int l = 0; l < nChem; l++)
                {
                        for(int m = 0; m < nSteps; m++)
                        {
                                m_vpdDataSumC[i][l][m] = 0;
                                m_vpdDataSumCSqr[i][l][m] = 0;
                        }
                }
        }
}

void CEnsembleCombinationDirector::LoadEnsembleInfo(std::string parameterFile)
{
        // read in the vectors of parameters - this assumes absolutely everybody has
        // exactly the same number of parameters; all runs are just different facets
        // of a single network
        int nRC = 0;
        if(m_vpRunnables.size() > 0)
        {
                nRC = m_vpRunnables[0]->GetReactionNetwork()->GetNumberOfRateConstants();
        }
        else
        {
                nRC = m_vpMinimizables[0]->GetNParameters();
        }
        ParameterReader *pEnsReader = new ParameterReader(parameterFile.c_str());
        for(int i = 0; i < m_iEnsembleSize; i++)
        {
                m_vpdEnsembleParameters.push_back(new double[nRC]);
                for(int j = 0; j < nRC; j++)
                {
                        double parameter = pEnsReader->ReadParameter();
                        m_vpdEnsembleParameters[i][j] = parameter;
                }
        }
        delete pEnsReader;

        return;
}

void CEnsembleCombinationDirector::Execute()
{
        int i;
        // run over all runnables and all minimizables for each ensemble member
        for(int eCount = 0; eCount < m_iEnsembleSize; eCount++)
        {
                // this is for all the runnables
                for(i = 0; i < m_vpRunnables.size(); i++)
                {
                        // set the parameters to the correct ensemble parameters
                        int nRC = m_vpRunnables[i]->GetReactionNetwork()->GetNumberOfRateConstants();
                        for(int j = 0; j < nRC; j++)
                        {
                                double rate = m_vpdEnsembleParameters[eCount][j];
                                m_vpRunnables[i]->GetReactionNetwork()->GetRateConstant(j)->SetRateConstant(rate);
                        }
                
                        // run the simulation
                        m_vpRunnables[i]->Run();
                        
                        // now increment the statistical arrays; each curve is normalized so 
                        // so that C(tmax) = 1.0
                        int nChem = m_vpRunnables[i]->GetReactionNetwork()->GetNumberOfChemicals();
                        int nSteps = m_vpRunnables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                        for(int cCount = 0; cCount < nChem; cCount++)
                        {
                                //double norm = m_vpRunnables[i]->GetCellObserver()->GetAverageConcentration()[cCount][nSteps-1];
                                double norm = 1.0;
                                for(int sCount = 0; sCount < nSteps; sCount++)
                                {
                                        double conc = m_vpRunnables[i]->GetCellObserver()->GetAverageConcentration()[cCount][sCount];
                                        m_vpdNoDataSumC[i][cCount][sCount] += conc/norm;
                                        m_vpdNoDataSumCSqr[i][cCount][sCount] += (conc/norm)*(conc/norm);
                                }
                        }
                }
                // this is for all the minimizables
                for(i = 0; i < m_vpMinimizables.size(); i++)
                {
                        // set the parameters to the correct ensemble parameters
                        int nRC = m_vpMinimizables[i]->GetNParameters();
                        double *parameters = new double[nRC];
                        for(int j = 0; j < nRC; j++)
                        {
                                double rate = m_vpdEnsembleParameters[eCount][j];
                                parameters[j] = rate;
                        }
                
                        // run the simulation
                        double cost = m_vpMinimizables[i]->ObjectiveFunction(parameters);
                        
                        // now increment the statistical arrays; multiply by the conversion 
                        // factor before incrementing the array
                        int nChem = m_vpMinimizables[i]->GetExperiment()->GetChemicalTimeSeriesData()->GetNChemicals();
                        int nSteps = m_vpMinimizables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                        for(int cCount = 0; cCount < nChem; cCount++)
                        {
                                for(int sCount = 0; sCount < nSteps; sCount++)
                                {
                                        double conc = m_vpMinimizables[i]->GetCellObserver()->GetAverageConcentration()[cCount][sCount];
                                        double scale = m_vpMinimizables[i]->GetConversionFactor(cCount)->GetFactorValue();
                                        m_vpdDataSumC[i][cCount][sCount] += conc*scale;
                                        m_vpdDataSumCSqr[i][cCount][sCount] += (conc*scale)*(conc*scale);
                                }
                        }
                        delete [] parameters;
                }
                // XXX debug
                cout << "Ensemble member " << eCount << " run." << endl;
        }

        // convert the sums/squared sums to average/std. dev.
        ConvertToStatistics();
        // Notify observers (like plotters)
        Notify();
        // Now dump all the data to a file, one file for each minimizable/runnable
        // Format should be:  t   avg   avg+std   avg-std
        // Separator is \n\n#, with the name and chemical number of the 
        //    measurement on the second line so they can later be identified
        
        // runnables
        char filename[200];
        for(i = 0; i < m_vpRunnables.size(); i++)
        {
                // set the name of the file
                sprintf(filename,"runnable%i.dat",i);
                // open the file
                FILE *pfData = fopen(filename,"w");
                // write a short header
                fprintf(pfData,"##################################################################\n");
                fprintf(pfData,"# ENSEMBLE CHEMICAL DATA FOR RUNNABLE %i\n",i);
                fprintf(pfData,"##################################################################\n");
                // number of chemicals/number of data points
                int nChem = m_vpRunnables[i]->GetReactionNetwork()->GetNumberOfChemicals();
                int nSteps = m_vpRunnables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                // loop over chemicals and print the necessary info
                for(int cCount = 0; cCount < nChem; cCount++)
                {
                        // write chemical number/name on comment line
                        fprintf(pfData,"# %i %s\n",cCount,m_vpRunnables[i]->GetReactionNetwork()->GetChemical(cCount)->GetName().c_str());
                        // loop over number of measurements and write
                        for(int sCount = 0; sCount < nSteps; sCount++)
                        {
                                double t = (double) sCount;
                                double avg = m_vpdNoDataSumC[i][cCount][sCount];
                                double sig = m_vpdNoDataSumCSqr[i][cCount][sCount];
                                fprintf(pfData,"%16.6e %16.6e %16.6e %16.6e\n",t,avg,avg+sig,avg-sig);
                        }
                        fprintf(pfData,"\n\n");
                }
                fclose(pfData);
        }
        // minimizables
        for(i = 0; i < m_vpMinimizables.size(); i++)
        {
                // set the name of the file
                sprintf(filename,"minimizable%i.dat",i);
                // open the file for reading - use std c functions because C++
                // has awful numerical formatting
                FILE *pfData = fopen(filename,"w");
                fprintf(pfData,"##################################################################\n");
                fprintf(pfData,"# ENSEMBLE CHEMICAL DATA FOR MINIMIZABLE %i\n",i);
                fprintf(pfData,"##################################################################\n");
                // number of chemicals/number of data points
                int nChem = m_vpMinimizables[i]->GetExperiment()->GetChemicalTimeSeriesData()->GetNChemicals();
                int nSteps = m_vpMinimizables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                // loop over chemicals and print the necessary info
                for(int cCount = 0; cCount < nChem; cCount++)
                {
                        // write chemical number/name on a comment line
                        fprintf(pfData,"# %i %s\n",cCount,m_vpMinimizables[i]->GetExperiment()->GetReactionNetwork()->GetChemical(cCount)->GetName().c_str());
                        // loop over number of measurements and write
                        for(int sCount = 0; sCount < nSteps; sCount++)
                        {
                                double t = (double) sCount;
                                double avg = m_vpdDataSumC[i][cCount][sCount];
                                double sig = m_vpdDataSumCSqr[i][cCount][sCount];
                                fprintf(pfData,"%16.6e %16.6e %16.6e %16.6e\n",t,avg,avg+sig,avg-sig);
                        }
                        // double newlines to next data set
                        fprintf(pfData,"\n\n");
                }               
                fclose(pfData);
        }
        
        cout << "Data dumped to files." << endl;

        return;
}


void CEnsembleCombinationDirector::ConvertToStatistics()
{
        int i;
        // now convert all the arrays
        for(i = 0; i < m_vpRunnables.size(); i++)
        {
                int nChem = m_vpRunnables[i]->GetReactionNetwork()->GetNumberOfChemicals();
                int nSteps = m_vpRunnables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                for(int chem = 0; chem < nChem; chem++)
                {
                        for(int step = 0; step < nSteps; step++)
                        {
                                // find averages
                                m_vpdNoDataSumC[i][chem][step] = m_vpdNoDataSumC[i][chem][step]/m_iEnsembleSize;
                                m_vpdNoDataSumCSqr[i][chem][step] = m_vpdNoDataSumCSqr[i][chem][step]/m_iEnsembleSize;
                                // convert sum(C^2) to sigma (fabs ensures numerical solvency for constants)
                                m_vpdNoDataSumCSqr[i][chem][step] -= m_vpdNoDataSumC[i][chem][step]*m_vpdNoDataSumC[i][chem][step];
                                m_vpdNoDataSumCSqr[i][chem][step] = sqrt(fabs(m_vpdNoDataSumCSqr[i][chem][step]));
                        }
                }
        }
        for(i = 0; i < m_vpMinimizables.size(); i++)
        {
                int nChem = m_vpMinimizables[i]->GetExperiment()->GetChemicalTimeSeriesData()->GetNChemicals();
                int nSteps = m_vpMinimizables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                for(int chem = 0; chem < nChem; chem++)
                {
                        for(int step = 0; step < nSteps; step++)
                        {
                                // find averages
                                m_vpdDataSumC[i][chem][step] = m_vpdDataSumC[i][chem][step]/m_iEnsembleSize;
                                m_vpdDataSumCSqr[i][chem][step] = m_vpdDataSumCSqr[i][chem][step]/m_iEnsembleSize;
                                // convert sum(C^2) to sigma (fabs ensures numerical solvency for constants)
                                m_vpdDataSumCSqr[i][chem][step] -= m_vpdDataSumC[i][chem][step]*m_vpdDataSumC[i][chem][step];
                                m_vpdDataSumCSqr[i][chem][step] = sqrt(fabs(m_vpdDataSumCSqr[i][chem][step]));
                        }
                }
        }
        
        cout << "Sum arrays converted to <C> and sigma(C)" << endl;
}

---