C++ DagNode::getLnProbabilityRatio方法代码示例

double CompoundMove::performMove( double &probRatio ) {
    
    if (changed) 
    {
        throw RbException("Trying to execute a Compound moves twice without accept/reject in the meantime.");
    }
    changed = true;
    
    double hr = performCompoundMove();
    
    if ( hr != hr || hr == RbConstants::Double::inf ) 
    {
        return RbConstants::Double::neginf;
    }
    
    std::set<DagNode* > affectedNodes;
    for (std::set<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); it++)
    {
        // touch each node
        (*it)->touch();
    
        // calculate the probability ratio for the node we just changed
        //std::cout << (*it)->getName() << " " << (*it)->getLnProbabilityRatio() << " " << (*it)->getLnProbability() << "\n";
            
        probRatio += (*it)->getLnProbabilityRatio();
 
        if ( probRatio != RbConstants::Double::inf && probRatio != RbConstants::Double::neginf )
        {
            // should contain unique nodes, since it is a set
            (*it)->getAffectedNodes(affectedNodes);
        }
    }
    
    for (std::set<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); it++)
    {
//        if ( nodes.find(oldN) == nodes.end() ) {
//            throw RbException("Cannot replace DAG node in this move because the move doesn't hold this DAG node.");
//        }
        
        affectedNodes.erase( *it );
    }
    
    if ( probRatio != RbConstants::Double::inf && probRatio != RbConstants::Double::neginf )
    {
        for (std::set<DagNode* >::iterator i=affectedNodes.begin(); i!=affectedNodes.end(); ++i)
        {
            DagNode* theAffectedNode = *i;
            
            // do not double-count the prob ratio for any elt in both theNodes and affectedNodes
            //if ( find(theNodes.begin(), theNodes.end(), *i) == theNodes.end() )
            {
                //std::cout << "  " << theAffectedNode->getName() << "  " << theAffectedNode->getLnProbabilityRatio() <<  " " << theAffectedNode->getLnProbability() << "\n";
                probRatio += theAffectedNode->getLnProbabilityRatio();
            }
        }
    }

    
    return hr;
}

double AdmixtureEdgeReplaceResidualsFNPR::performMove( double &probRatio ) {
    
    
    if (changed) {
        throw RbException("Trying to execute a simple moves twice without accept/reject in the meantime.");
    }
    changed = true;
    
    double hr = performSimpleMove();
    
    if ( hr != hr || hr == RbConstants::Double::inf ) {
        return RbConstants::Double::neginf;
    }
    
    // touch the node
    variable->touch();
    probRatio = variable->getLnProbabilityRatio();
    
    if ( probRatio != RbConstants::Double::inf && probRatio != RbConstants::Double::neginf ) {
        
        std::set<DagNode* > affectedNodes;
        variable->getAffectedNodes(affectedNodes);
        for (std::set<DagNode* >::iterator i=affectedNodes.begin(); i!=affectedNodes.end(); ++i) {
            DagNode* theNode = *i;
            probRatio += theNode->getLnProbabilityRatio();
            //std::cout << theNode->getName() << "\t" << probRatio << "\n";
        }
    }
    
    return hr;
}

double AdmixtureNearestNeighborInterchangeAndRateShift::performMove( double &probRatio ) {
    
    //std::cout << "NNI::performMove()\n";
    
    if (changed) {
   //;     throw RbException("Trying to execute a simple moves twice without accept/reject in the meantime.");
    }
    changed = true;
    
    double hr = performSimpleMove();
    
    if ( hr != hr || hr == RbConstants::Double::inf ) {
        return RbConstants::Double::neginf;
    }
    
    // touch the node
    variable->touch();

    // calculate the probability ratio for the node we just changed
    probRatio = variable->getLnProbabilityRatio();
    //std::cout << "n\t" << variable->getName() << "\t" << probRatio << "\n";
 
    if ( probRatio != RbConstants::Double::inf && probRatio != RbConstants::Double::neginf ) {
        
        if (!failed)
        {
            // these three rates should be touched due to setValue()
            branchRates[storedChildRateIndex]->touch();
            probRatio += branchRates[storedChildRateIndex]->getLnProbabilityRatio();
            branchRates[storedBrotherRateIndex]->touch();
            probRatio += branchRates[storedBrotherRateIndex]->getLnProbabilityRatio();
           // probRatio += branchRates[storedNodeRateIndex]->getLnProbabilityRatio();
           // std::cout << probRatio << "\n";
        }
    
        std::set<DagNode* > affectedNodes;
        variable->getAffectedNodes(affectedNodes);
        for (std::set<DagNode* >::iterator i=affectedNodes.begin(); i!=affectedNodes.end(); ++i) {
            DagNode* theNode = *i;
            probRatio += theNode->getLnProbabilityRatio();
            //std::cout << "\tch\t" << theNode->getName() << "\t" << theNode->getLnProbabilityRatio() << "  " <<  probRatio << "\n";
        }
    }
    //std::cout << "pr " << probRatio << "    hr " << hr << "\n";
    
    return hr;
}

double AdmixtureEdgeReplaceResidualWeights::performMove( double &probRatio ) {
    
    //std::cout << "\nAdmix Edge Replace\n";
    
    if (changed) {
        throw RbException("Trying to execute a simple moves twice without accept/reject in the meantime.");
    }
    changed = true;
    
    double hr = performSimpleMove();
    
    if ( hr != hr || hr == RbConstants::Double::inf || hr == RbConstants::Double::neginf ) {
        return RbConstants::Double::neginf;
    }
    
    // touch the node
    variable->touch(); // if previously touched, this will overwrite lnProb??
    probRatio = variable->getLnProbabilityRatio();
    //probRatio = 0.0;
    
    for (std::map<int,double>::iterator it = storedBranchRates.begin(); it != storedBranchRates.end(); it++)
    {
        branchRates[it->first]->touch();
        probRatio += branchRates[it->first]->getLnProbabilityRatio();
        //std::cout << branchRates[it->first]->getLnProbabilityRatio() << "\n";
    }

    if ( probRatio != RbConstants::Double::inf && probRatio != RbConstants::Double::neginf ) {
        
        std::set<DagNode* > affectedNodes;
        variable->getAffectedNodes(affectedNodes);
        for (std::set<DagNode* >::iterator i=affectedNodes.begin(); i!=affectedNodes.end(); ++i) {
            DagNode* theNode = *i;
            probRatio += theNode->getLnProbabilityRatio();
            //std::cout << probRatio << " " << theNode->getName() << "\t" << theNode->getLnProbability() << " " << theNode->getLnProbabilityRatio() << "\n";
        }
    }
   // std::cout << probRatio << "\n";

    
    return hr;
}

double AdmixtureShiftTreeRates::performMove( double &probRatio ) {
    
    if (changed) {
        throw RbException("Trying to execute a simple moves twice without accept/reject in the meantime.");
    }
    changed = true;
    
    double hr = performSimpleMove();
    
    if ( hr != hr || hr == RbConstants::Double::inf ) {
        return RbConstants::Double::neginf;
    }
    
    // touch the node
    treeRate->touch();
    
    // calculate the probability ratio for the node we just changed
    probRatio = treeRate->getLnProbabilityRatio();
    //std::cout << "n\t" << treeRate->getName() << "\t" << probRatio << "\n";
    
    //for (size_t i = 0; i < branchRates.size(); i++)
    for (std::map<int,double>::iterator it = storedRates.begin(); it != storedRates.end(); it++)
    {
        int idx = it->first;
        branchRates[idx]->touch();
        probRatio += branchRates[idx]->getLnProbabilityRatio();
        //std::cout << "n\t" << branchRates[idx]->getName() << " " << branchRates[idx]->getLnProbability() << " " << exp(branchRates[idx]->getLnProbability()) << " " << probRatio << "\n";
    }
    
    if ( probRatio != RbConstants::Double::inf && probRatio != RbConstants::Double::neginf ) {
        
        std::set<DagNode* > affectedNodes;
        treeRate->getAffectedNodes(affectedNodes);
        for (std::set<DagNode* >::iterator i=affectedNodes.begin(); i!=affectedNodes.end(); ++i) {
            DagNode* theNode = *i;
            probRatio += theNode->getLnProbabilityRatio();
           // std::cout << "\tch\t" << theNode->getName() << "\t" << theNode->getLnProbabilityRatio() << "  " <<  probRatio << "\n";
        }
    }
    
    return hr;
}

double FossilSafeScaleMove::performMove( double &probRatio ) {
    
    if (changed)
    {
        throw RbException("Trying to execute a simple move twice without accept/reject in the meantime.");
    }
    changed = true;
    
    double hr = doMove();
    
    if ( hr != hr || hr == RbConstants::Double::inf )
    {
        return RbConstants::Double::neginf;
    }
    
    // touch the node
    scaler->touch();
    tree->touch();
    
    // calculate the probability ratio for the node we just changed
    probRatio = scaler->getLnProbabilityRatio();
    probRatio += tree->getLnProbabilityRatio();
    
    if ( probRatio != RbConstants::Double::inf && probRatio != RbConstants::Double::neginf )
    {
        
        std::set<DagNode* > affectedNodes;
        scaler->getAffectedNodes(affectedNodes);
        tree->getAffectedNodes(affectedNodes);
        for (std::set<DagNode* >::iterator i=affectedNodes.begin(); i!=affectedNodes.end(); ++i)
        {
            DagNode* theAffectedNode = *i;
            if (theAffectedNode != tree && theAffectedNode != scaler) {
                //std::cout << theAffectedNode->getName() << "  " << theAffectedNode->getLnProbabilityRatio() << " " << theAffectedNode->getLnProbability() << "\n";
                probRatio += theAffectedNode->getLnProbabilityRatio();
            }
        }
    }
    
    return hr;
}