---

EnsembleRunDirector.cpp

Go to the documentation of this file.

// EnsembleRunDirector.cpp: implementation of the CEnsembleRunDirector class.
//

#include "EnsembleRunDirector.h"

// Construction/Destruction

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

CEnsembleRunDirector::~CEnsembleRunDirector()
{
        int i;
        for(i = 0; i < m_vpdEnsembleParameters.size(); i++)
        {
                delete [] m_vpdEnsembleParameters[i];
        }
        for(i = 0; i < m_vpdSumC.size(); i++)
        {
                delete [] m_vpdSumC[i][0];
                delete [] m_vpdSumC[i];
                delete [] m_vpdSumCSqr[i][0];
                delete [] m_vpdSumCSqr[i];
        }
        for(i = 0; i < m_vpEnsemblePlotters.size(); i++)
        {
                delete m_vpEnsemblePlotters[i];
        }
}

void CEnsembleRunDirector::SumAllocate()
{
        // runnables
        for(int i = 0; i < m_vpRunnables.size(); i++)
        {
                int nSteps = m_vpRunnables[i]->GetCellObserver()->GetNumberOfTimeSteps();
                int nChem = m_vpRunnables[i]->GetReactionNetwork()->GetNumberOfChemicals();
                m_vpdSumC.push_back(new double *[nChem]);
                m_vpdSumCSqr.push_back(new double *[nChem]);
                m_vpdSumC[i][0] = new double[nChem*nSteps];
                m_vpdSumCSqr[i][0] = new double[nChem*nSteps];
                for(int chem = 1; chem < nChem; chem++)
                {
                        m_vpdSumC[i][chem] = &(m_vpdSumC[i][0][chem*nSteps]);
                        m_vpdSumCSqr[i][chem] = &(m_vpdSumCSqr[i][0][chem*nSteps]);
                }
                // zero the array just created
                for(int l = 0; l < nChem; l++)
                {
                        for(int m = 0; m < nSteps; m++)
                        {
                                m_vpdSumC[i][l][m] = 0;
                                m_vpdSumCSqr[i][l][m] = 0;
                        }
                }
        }
}

void CEnsembleRunDirector::LoadEnsembleInfo(std::string parameterFile)
{
        // read in the vectors of parameters
        int nRC = m_vpRunnables[0]->GetReactionNetwork()->GetNumberOfRateConstants();
        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 CEnsembleRunDirector::Execute()
{
        // run over all runnables for each ensemble member
        for(int eCount = 0; eCount < m_iEnsembleSize; eCount++)
        {
                // run all the runnables
                for(int 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];
                                for(int sCount = 0; sCount < nSteps; sCount++)
                                {
                                        double conc = m_vpRunnables[i]->GetCellObserver()->GetAverageConcentration()[cCount][sCount];
                                        m_vpdSumC[i][cCount][sCount] += conc/norm;
                                        m_vpdSumCSqr[i][cCount][sCount] += (conc/norm)*(conc/norm);
                                }
                        }
                }
                // XXX debug
                cout << "Ensemble member " << eCount << " run." << endl;
        }

        // convert the sums/squared sums to average/std. dev.
        ConvertToStatistics();
        // Notify observers
        Notify();

        return;
}

void CEnsembleRunDirector::ConvertToStatistics()
{
        // now convert all the arrays
        for(int 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_vpdSumC[i][chem][step] = m_vpdSumC[i][chem][step]/m_iEnsembleSize;
                                m_vpdSumCSqr[i][chem][step] = m_vpdSumCSqr[i][chem][step]/m_iEnsembleSize;
                                // convert sum(C^2) to sigma
                                m_vpdSumCSqr[i][chem][step] -= m_vpdSumC[i][chem][step]*m_vpdSumC[i][chem][step];
                                m_vpdSumCSqr[i][chem][step] = sqrt(m_vpdSumCSqr[i][chem][step]);
                        }
                }
        }
        
        cout << "Sum arrays converted to <C> and sigma(C)" << endl;
}

---