---

NetworkMinimizableDirector.cpp

Go to the documentation of this file.

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

#include "NetworkMinimizableDirector.h"

// Construction/Destruction

NetworkMinimizableDirector::NetworkMinimizableDirector()
{
 
}

NetworkMinimizableDirector::~NetworkMinimizableDirector()
{
        int i;
        for(i = 0; i < _minimizableList.size(); i++)
        {
                delete _minimizableList[i];
        }
        for(i = 0; i < _networkList.size(); i++)
        {
                delete _networkList[i];
        }
        for(i = 0; i < _experimentList.size(); i++)
        {
                delete _experimentList[i];
        }
        for(i = 0; i < _moverList.size(); i++)
        {
                delete _moverList[i];
        }
}

void NetworkMinimizableDirector::DumpResidualInfo()
{
        fstream *pFile = new fstream("resinfo.par",ios::out);
        *pFile << "##################################################################" << endl;
        *pFile << "# RESIDUALS PRINTED IN THE ORDER THEY ARE STORED INTERNALLY" << endl;
        *pFile << "# C1 IS MINIMIZABLE NUM., C2 IS CHEMICAL NUM., C3 IS TIME" << endl;
        *pFile << "##################################################################" << endl;

        for(int i = 0; i < _experimentList.size(); i++)
        {
                int tvecsize = _experimentList[i]->GetChemicalTimeSeriesData()->GetTimeVector().size();
                for(int j = 0; j < tvecsize; j++)
                {
                        int cnum = _experimentList[i]->GetChemicalTimeSeriesData()->GetTimeVector()[j]->GetChemNumber();
                        double time = _experimentList[i]->GetChemicalTimeSeriesData()->GetTimeVector()[j]->GetTime();
                        *pFile << i << "\t\t" << cnum << "\t\t" << time << endl;
                }
        }
        delete pFile;
}

double NetworkMinimizableDirector::ObjectiveFunction(double *parameters)
{
        double summedCost = 0.0;
        summedCost = this->ComputeResiduals(parameters);
        return summedCost;
}

double NetworkMinimizableDirector::ComputeES(double *parameters, double T)
{
        m_dEnergyLastComputed = this->ComputeResiduals(parameters);
        double S = 0;
        for(int i = 0; i < _minimizableList.size(); i++)
        {
                S += _minimizableList[i]->EntropyShift(T) - 0.5*(_minimizableList[i]->GetEntropyLastComputed());
        }
        m_dEntropyLastComputed = S;
        return m_dEnergyLastComputed;
}

double NetworkMinimizableDirector::EntropyShift(double T)
{
        // total up the net entropy shift
        double S = 0;
        for(int i = 0; i < _minimizableList.size(); i++)
        {
                S += 0.5*(_minimizableList[i]->GetNBFactors())*log(2*3.14159265358979323846*T);
        }
        return S;
}

double NetworkMinimizableDirector::F(double *parameters, double T)
{
        m_dEnergyLastComputed = this->ComputeResiduals(parameters);
        double S = 0;
        for(int i = 0; i < _minimizableList.size(); i++)
        {
                S += -0.5*(_minimizableList[i]->GetEntropyLastComputed());
        }
        m_dEntropyLastComputed = S;
        return (m_dEnergyLastComputed - T*m_dEntropyLastComputed);
}

double NetworkMinimizableDirector::F0(double *parameters, double T)
{
        m_dEnergyLastComputed = this->ComputeResiduals(parameters);
        double S = 0;
        for(int i = 0; i < _minimizableList.size(); i++)
        {
                S += _minimizableList[i]->EntropyShift(T) - 0.5*(_minimizableList[i]->GetEntropyLastComputed());
        }
        m_dEntropyLastComputed = S;
        return (m_dEnergyLastComputed - T*m_dEntropyLastComputed);
}

double NetworkMinimizableDirector::ComputeResiduals(double *parameters)
{
        int i;
        // just sum up the costs from each director and then add the 
        // integration penalty, = (gamma^2/2)*(sum_ti)^2
        double summedCost = 0.0;
        double timeCost = 0.0;
        double pcost = 0.0;
        
        // compute the residuals/cost from each minimizable
        for(i = 0; i < _minimizableList.size(); i++)
        {
                summedCost += _minimizableList[i]->ObjectiveFunction(parameters);
        }
        
        // load up the master list with each set of residuals
        int nextResidual = 0;
        for(i = 0; i < _minimizableList.size(); i++)
        {
                int nRes = _minimizableList[i]->GetNResiduals();
                for(int j = 0; j < nRes; j++)
                {
                        residuals[nextResidual+j] = (_minimizableList[i]->GetResiduals())[j];
                }
                nextResidual += nRes;
        }
        
        // now compute the cost for deviation from guessed values, if any
        for(i = 0; i < m_iPriorList.size(); i++)
        {
                int which = m_iPriorList[i];
                //double presid = parameters[which] - _networkList[0]->GetRateConstant(which)->GetInitialValue();
                double presid = _networkList[0]->GetRateConstant(which)->GetValue() - _networkList[0]->GetRateConstant(which)->GetInitialValue();
                presid = presid/(sqrt(2.0)*_networkList[0]->GetRateConstant(which)->GetErrorInInitialValue());
                residuals[nextResidual + i] = presid;
                pcost += presid*presid;
        }
        nextResidual += m_iPriorList.size();
        
        // if the prefactor for the integration cost is > 0, compute its residual (otherwise
        // that residual shouldn't exist)
        if(m_dGammaSquared > 0.0)
        {
                double totalSteps = 0.0;
                for(int l = 0; l < _moverList.size(); l++)
                {
                        totalSteps += _moverList[l]->GetTotalStepsTaken();
                }
                residuals[nResiduals - 1] = sqrt(m_dGammaSquared/2.0)*totalSteps;
                timeCost = residuals[nResiduals - 1]*residuals[nResiduals - 1];
        }
        return summedCost + pcost + timeCost;
}

---