---

CoupledReplicasMinimizableDirector.cpp

Go to the documentation of this file.

// CoupledReplicasMinimizableDirector.cpp: implementation of the CCoupledReplicasMinimizableDirector class.
//

#include "CoupledReplicasMinimizableDirector.h"

// Construction/Destruction

CCoupledReplicasMinimizableDirector::CCoupledReplicasMinimizableDirector(double epsSq, double deltaSq)
{
        m_dEpsSquared = epsSq;
        m_dDeltaSquared = deltaSq;
}

CCoupledReplicasMinimizableDirector::~CCoupledReplicasMinimizableDirector()
{
        delete m_pReplicaOne;
        delete m_pReplicaTwo;
}

int CCoupledReplicasMinimizableDirector::GetNParameters()
{
        return (m_pReplicaOne->GetNParameters() + m_pReplicaTwo->GetNParameters());
}

double CCoupledReplicasMinimizableDirector::GetParameter(int parIndex)
{
        // it's in the first list
        int firstNP = m_pReplicaOne->GetNParameters();
        if(parIndex < firstNP)
        {
                return m_pReplicaOne->GetParameter(parIndex);
        }
        else // it's in the second list
        {
                return m_pReplicaTwo->GetParameter(parIndex - firstNP); 
        }

}

void CCoupledReplicasMinimizableDirector::SetIntersectionLists()
{
        // look at all the parameters in the first minimizable/replica and determine
        // if they have counterparts in the second network.  When they are found, push
        // them onto the intersection lists
        for(int i = 0; i < m_pReplicaOne->GetReactionNetwork(0)->GetNumberOfRateConstants(); i++)
        {
                std::string firstName = m_pReplicaOne->GetReactionNetwork(0)->GetRateConstant(i)->GetName();
                for(int j = 0; j < m_pReplicaTwo->GetReactionNetwork(0)->GetNumberOfRateConstants(); j++)
                {
                        if(firstName == m_pReplicaTwo->GetReactionNetwork(0)->GetRateConstant(j)->GetName())
                        {
                                // put i on the first intersection list and j on the second
                                // intersection list
                                m_viIntersectOne.push_back(i);
                                m_viIntersectTwo.push_back(j);
                        }
                }
        }

        // before returning, if the intersection list is null print a warning
        if(m_viIntersectOne.size() == 0)
        {
                cout << "WARNING : Coupled models have no shared parameters." << endl;
        }
}

double CCoupledReplicasMinimizableDirector::ComputeResiduals(double *parameters)
{
        int j;
        double summedCost = 0.0;
        double psiSquared = 0.0;
        double psiFourth = 0.0;
        double PCost = 0.0;
        double gammaCost = 0.0;
        int subResSize = 0;
        
        // compute the residuals/cost from each replica
        int startingIndex = m_pReplicaOne->GetNParameters();
        summedCost = m_pReplicaOne->ObjectiveFunction(&parameters[0]) + m_pReplicaTwo->ObjectiveFunction(&parameters[startingIndex]);
        
        int ResOne = m_pReplicaOne->GetNResiduals();
        int ResTwo = m_pReplicaTwo->GetNResiduals();

        for(j = 0; j < ResOne; j++)
        {
                residuals[j] = (m_pReplicaOne->GetResiduals())[j];
        }
        for(j = ResOne; j < ResOne + ResTwo; j++)
        {
                residuals[j] = (m_pReplicaTwo->GetResiduals())[j-ResOne];
        }
        int nextResidual = ResOne + ResTwo;
        
        double rcCost = 0.0;
        int intersectionSize = m_viIntersectOne.size();
        int n1 = m_pReplicaOne->GetNParameters();
        for(j = 0; j < intersectionSize; j++)
        {
                double psi = log(parameters[m_viIntersectOne[j]]) - log(parameters[n1 + m_viIntersectTwo[j]]);
                psiSquared += psi*psi;
                psiFourth += psi*psi*psi*psi;
                residuals[nextResidual+j] = sqrt(m_dEpsSquared/2.0)*psi;
                rcCost += residuals[nextResidual+j]*residuals[nextResidual+j];
        }
        // last residual is the participation ratio
        residuals[nResiduals-1] = sqrt(m_dDeltaSquared/2.0)*(2*log(psiSquared) - log(psiFourth));
        PCost = residuals[nResiduals-1]*residuals[nResiduals-1];
        
        return summedCost + rcCost + PCost;
}

---