C++ DagNode类代码示例

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;
}

/** 
 * Monitor value at generation gen 
 */
void AbstractFileMonitor::monitorVariables(unsigned long gen)
{
        
    for (std::vector<DagNode*>::iterator i = nodes.begin(); i != nodes.end(); ++i)
    {
        // add a separator before every new element
        out_stream << separator;
            
        // get the node
        DagNode *node = *i;
            
        // print the value
        node->printValueElements(out_stream, separator, -1, true, flatten);
    }
    
}

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 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 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;
}

/**
 * Basic destructor doing nothing.
 */
MetropolisHastingsMove::~MetropolisHastingsMove( void )
{
    for (std::set<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it)
    {
        // get the pointer to the current node
        DagNode* theNode = *it;

        // add myself to the set of moves
        theNode->removeMove( this );

        // decrease the DAG node reference count because we also have a pointer to it
        if ( theNode->decrementReferenceCount() == 0 )
        {
            delete theNode;
        }

    }

    delete proposal;
}

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;
}

/**
 * Constructor
 *
 * Here we simply allocate and initialize the move object.
 *
 * \param[in]    w   The weight how often the proposal will be used (per iteration).
 * \param[in]    t   If auto tuning should be used.
 */
AbstractMove::AbstractMove( const std::vector<DagNode*> &n, double w, bool t ) :
    nodes( n ),
    affectedNodes( ),
    weight( w ),
    autoTuning( t )
{
    
    for (std::vector<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it)
    {
        // get the pointer to the current node
        DagNode* theNode = *it;
        
        // add myself to the set of moves
        theNode->addMove( this );
        
        // increase the DAG node reference count because we also have a pointer to it
        theNode->incrementReferenceCount();
        
    }
    
    
    // remove all "core" nodes from affectedNodes so their probabilities are not double-counted
    for (size_t i = 0; i < affectedNodes.size(); ++i)
    {
        RbOrderedSet<DagNode*>::iterator it = affectedNodes.begin();
        std::advance(it, i);
        
        for (size_t j = 0; j < nodes.size(); ++j)
        {
            if ( nodes[j] == *it )
            {
                affectedNodes.erase(*it);
                --i;
                break;
            }
            
        }
        
    }
    
}

/**
 * Copy constructor.
 * We need to create a deep copy of the proposal here.
 *
 * \param[in]   m   The object to copy.
 *
 */
MetropolisHastingsMove::MetropolisHastingsMove(const MetropolisHastingsMove &m) : AbstractMove(m),
    affectedNodes( m.affectedNodes ),
    nodes( m.nodes ),
    numAccepted( m.numAccepted ),
    proposal( m.proposal->clone() )
{

    for (std::set<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it)
    {
        // get the pointer to the current node
        DagNode* theNode = *it;

        // add myself to the set of moves
        theNode->addMove( this );

        // increase the DAG node reference count because we also have a pointer to it
        theNode->incrementReferenceCount();

    }

}

/**
 * Basic destructor doing nothing.
 */
AbstractMove::~AbstractMove( void )
{
    
    // clean up my pointers to the nodes
    for (size_t i = 0; i < nodes.size(); ++i)
    {
        // get the pointer to the current node
        DagNode* theNode = nodes[i];
        
        // add myself to the set of moves
        theNode->removeMove( this );
        
        // decrease the DAG node reference count because we also have a pointer to it
        if ( theNode->decrementReferenceCount() == 0 )
        {
            delete theNode;
        }
        
    }
    
}

AbstractMove::AbstractMove( const AbstractMove &m ) : Move( m ),
    nodes( m.nodes ),
    affectedNodes( m.affectedNodes ),
    weight( m.weight ),
    autoTuning( m.autoTuning  )
{
    
    
    for (size_t i = 0; i < nodes.size(); ++i)
    {
        // get the pointer to the current node
        DagNode* theNode = nodes[i];
        
        // add myself to the set of moves
        theNode->addMove( this );
        
        // increase the DAG node reference count because we also have a pointer to it
        theNode->incrementReferenceCount();
        
    }
    
}

/**
 * Reset the currently monitored DAG nodes by extracting the DAG nodes from the model again 
 * and store this in the set of DAG nodes.
 */
void ModelMonitor::resetDagNodes( void )
{
    
    // for savety we empty our dag nodes
    nodes.clear();
    
    if ( model != NULL )
    {
        // we only want to have each nodes once
        // this should by default happen by here we check again
        std::set<std::string> varNames;
        
        const std::vector<DagNode*> &n = model->getDagNodes();
        for (std::vector<DagNode*>::const_iterator it = n.begin(); it != n.end(); ++it) 
        {
            
            DagNode *theNode = *it;
            
            // only simple numeric variable can be monitored (i.e. only integer and real numbers)
            if ( theNode->isSimpleNumeric() && !theNode->isClamped())
            {
                if ( (!stochasticNodesOnly && !theNode->isConstant() && theNode->getName() != "" && !theNode->isHidden() && !theNode->isElementVariable() ) || ( theNode->isStochastic() && !theNode->isClamped() && theNode->isHidden() == false  && theNode->isElementVariable() == false ) )
                {
                    const std::string &name = theNode->getName();
                    if ( varNames.find( name ) == varNames.end() )
                    {
                        nodes.push_back( theNode );
                        varNames.insert( name );
                    }
                    else
                    {
                        /*
                        std::cerr << "Trying to add variable with name '" << name << "' twice." << std::endl;
                         */
                    }
                }
            }
        
        }
    }
    
}

/** Monitor value at generation gen */
void ScreenMonitor::monitor(unsigned long gen) {
    
    // get the printing frequency
    unsigned long samplingFrequency = printgen;
    
    if (gen % samplingFrequency == 0) {
        // print the iteration number first
        std::cout << gen;
        
        if ( posterior ) {
            // add a separator before every new element
            std::cout << separator;
            
            const std::vector<DagNode*> &n = model->getDagNodes();
            double pp = 0.0;
            for (std::vector<DagNode*>::const_iterator it = n.begin(); it != n.end(); ++it) {
                pp += (*it)->getLnProbability();
            }
            std::cout << pp;
        }
        
        if ( likelihood ) {
            // add a separator before every new element
            std::cout << separator;
            
            const std::vector<DagNode*> &n = model->getDagNodes();
            double pp = 0.0;
            for (std::vector<DagNode*>::const_iterator it = n.begin(); it != n.end(); ++it) {
                if ( (*it)->isClamped() ) {
                    pp += (*it)->getLnProbability();
                }
            }
            std::cout << pp;
        }
        
        if ( prior ) {
            // add a separator before every new element
            std::cout << separator;
            
            const std::vector<DagNode*> &n = model->getDagNodes();
            double pp = 0.0;
            for (std::vector<DagNode*>::const_iterator it = n.begin(); it != n.end(); ++it) {
                if ( !(*it)->isClamped() ) {
                    pp += (*it)->getLnProbability();
                }
            }
            std::cout << pp;
        }
        
        for (std::vector<DagNode*>::iterator i = nodes.begin(); i != nodes.end(); ++i) {
            // add a separator before every new element
            std::cout << separator;
            
            // get the node
            DagNode *node = *i;
            
            // print the value
            node->printValue(std::cout, separator);
        }
        
        std::cout << std::endl;
        
    }
}

/**
 * Swap the current variable for a new one.
 *
 * \param[in]     oldN     The old variable that needs to be replaced.
 * \param[in]     newN     The new variable.
 */
void AbstractMove::swapNode(DagNode *oldN, DagNode *newN)
{
    
    // find the old node
    
    for (size_t i = 0; i < nodes.size(); ++i)
    {
        // get the pointer to the current node
        DagNode* theNode = nodes[i];
        if ( theNode == oldN )
        {
            nodes[i] = newN;
        }
        
    }
    
    // remove myself from the old node and add myself to the new node
    oldN->removeMove( this );
    newN->addMove( this );
    
    // increment and decrement the reference counts
    newN->incrementReferenceCount();
    if ( oldN->decrementReferenceCount() == 0 )
    {
        throw RbException("Memory leak in Metropolis-Hastings move. Please report this bug to Sebastian.");
    }
    
    affectedNodes.clear();
    
    for (std::vector<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it)
    {
        // get the pointer to the current node
        DagNode* theNode = *it;
        
        // get the affected nodes if we would update this node
        // then we don't need to get the affected nodes every time again
        theNode->getAffectedNodes( affectedNodes );
    }
    
    // remove all "core" nodes from affectedNodes so their probabilities are not double-counted
    for (size_t i = 0; i < affectedNodes.size(); ++i)
    {
        RbOrderedSet<DagNode*>::iterator it = affectedNodes.begin();
        std::advance(it, i);
        
        for (size_t j = 0; j < nodes.size(); ++j)
        {
            if ( nodes[j] == *it )
            {
                affectedNodes.erase(*it);
                --i;
                break;
            }
            
        }
        
    }
    
    swapNodeInternal(oldN, newN);
    
}