---

ReplicaMinimizableDirector.cpp

Go to the documentation of this file.

// ReplicaMinimizableDirector.cpp: implementation of the CReplicaMinimizableDirector class.
//

#include "ReplicaMinimizableDirector.h"

// Construction/Destruction

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

CReplicaMinimizableDirector::~CReplicaMinimizableDirector()
{

}

int CReplicaMinimizableDirector::GetNParameters()
{
        int summedNP = 0;
        for(int i = 0; i < _minimizableList.size(); i++)
        {
                summedNP += _minimizableList[i]->GetNParameters();
        }
        return summedNP;
}

double CReplicaMinimizableDirector::GetParameter(int parIndex)
{
        // it's in the first list
        int firstNP = _minimizableList[0]->GetNParameters();
        if(parIndex < firstNP)
        {
                return _minimizableList[0]->GetParameter(parIndex);
        }
        else // it's in the second list
        {
                return _minimizableList[1]->GetParameter(parIndex - firstNP); 
        }
}

void CReplicaMinimizableDirector::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 < _networkList[0]->GetNumberOfRateConstants(); i++)
        {
                std::string firstName = _networkList[0]->GetRateConstant(i)->GetName();
                for(int j = 0; j < _networkList[1]->GetNumberOfRateConstants(); j++)
                {
                        if(firstName == _networkList[1]->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);
                        }
                }
        }
}

double CReplicaMinimizableDirector::ComputeResiduals(double *parameters)
{
        int i;
        // divide the residuals by 1/sqrt(N_i) and the summed cost by 1/N_i, where
        // N_i is the number of residuals in each replica (replicas are separated in 
        // this calculation in a way that sub-experiments in minimizables are not)
        double summedCost = 0.0;
        double psiSquared = 0.0;
        double psiFourth = 0.0;
        double Pcost = 0.0;
        double gammaCost = 0.0;
        int subResSize = 0;
        if(m_dGammaSquared > 0.0)
        {
                subResSize = this->nResiduals - 1;
        }
        else
        {
                subResSize = this->nResiduals;
        }
        // compute the residuals/cost from each replica
        for(i = 0; i < _minimizableList.size(); i++)
        {
                int startingIndex = 0;
                for(int j = 0; j < i; j++)
                {
                        startingIndex += _minimizableList[j]->GetNParameters();
                }
                double tempCost = _minimizableList[i]->ObjectiveFunction(&parameters[startingIndex]);
                summedCost += tempCost/(_minimizableList[i]->GetNResiduals());
        }
        
        // 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])/sqrt((double)nRes);
                }
                nextResidual += nRes;
        }
        // the intersection list is the same size for both replicas, so it doesn't matter
        // which one you ask for the size
        double rcCost = 0.0;
        int intersectionSize = m_viIntersectOne.size();
        int n1 = _minimizableList[0]->GetNParameters();
        for(int 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/intersectionSize)*psi;
                rcCost += residuals[nextResidual+j]*residuals[nextResidual+j];
        }
        // if the integration penalty is present, the P ratio is the next-to-last
        // residual.  Otherwise it's the last one
        if(m_dGammaSquared > 0.0)
        {
                // take care of P ratio
                residuals[this->nResiduals - 2] = sqrt(m_dDeltaSquared)*(2.0*log(psiSquared) - log(psiFourth));
                Pcost = residuals[this->nResiduals-2]*residuals[this->nResiduals-2];    
                // now the time cost
                double totalSteps = 0.0;
                for(int l = 0; l < _moverList.size(); l++)
                {
                        totalSteps += _moverList[l]->GetTotalStepsTaken();
                }
                residuals[this->nResiduals - 1] = sqrt(m_dGammaSquared)*(totalSteps/_moverList.size());
                gammaCost = residuals[this->nResiduals - 1]*residuals[this->nResiduals - 1];
        }
        else
        {
                // fill in the next to last residual, which is the participation ratio
                residuals[this->nResiduals - 1] = sqrt(m_dDeltaSquared)*(2.0*log(psiSquared) - log(psiFourth));
                Pcost = residuals[this->nResiduals - 1]*residuals[this->nResiduals - 1]; 
                // no time cost
                gammaCost = 0.0;
        }
        // XXX debug
        cout << "Summed cost = " << summedCost << ", RC cost = " << rcCost << endl;
        return summedCost + rcCost + Pcost + gammaCost;
}

---