Examples with MultiTypeNode beast.evolution.tree.MultiTypeNode used on opensource projects

Example 11 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class BeerliFelsenstein method getReverseMoveProb.

/**
 * Calculate probability with which the current state is proposed from
 * the new state.
 * @param eventList List of events on partial genealogy
 * @param migProp Pre-calculated Migration propensities
 * @param node Node below edge which is modified during proposal
 * @param nodeSis Sister of node selected for proposal
 * @param oldCoalTime Original coalescence time of selected edge
 * @param newRootHeight Height of tree following proposal
 * @return log of proposal density
 */
private double getReverseMoveProb(List<Event> eventList, double[] migProp, MultiTypeNode node, MultiTypeNode nodeSis, double oldCoalTime, double newRootHeight) {
    double logP = 0.0;
    MultiTypeNode mtNode = (MultiTypeNode) node;
    // Number of lineages in each deme - needed for coalescence propensity
    int[] lineageCounts = new int[migModel.getNTypes()];
    // Next change index
    int changeIdx = 0;
    // Current deme
    int deme = mtNode.getNodeType();
    // Flag to indicate this interval will terminate early due to
    // the start of a two-lineage simulation proposal phase
    boolean switchPhase = false;
    for (int eidx = 0; (eidx < eventList.size() - 1 && !switchPhase); eidx++) {
        Event event = eventList.get(eidx);
        double intervalEndTime = eventList.get(eidx + 1).time;
        if (intervalEndTime > newRootHeight) {
            intervalEndTime = newRootHeight;
            switchPhase = true;
        }
        switch(event.type) {
            case COALESCENCE:
                lineageCounts[event.thisDeme] -= 1;
                break;
            case SAMPLE:
                lineageCounts[event.thisDeme] += 1;
                break;
            case MIGRATION:
                lineageCounts[event.prevDeme] -= 1;
                lineageCounts[event.thisDeme] += 1;
                break;
        }
        if (node.getHeight() > intervalEndTime)
            continue;
        double t = Math.max(event.time, node.getHeight());
        // Loop over changes within this interval
        while (true) {
            // Calculate coalescence propensities
            double coalProp = lineageCounts[deme] / migModelSC.getPopSize(deme);
            double nextTime;
            if (changeIdx < mtNode.getChangeCount())
                nextTime = mtNode.getChangeTime(eidx);
            else
                nextTime = oldCoalTime;
            double dt = Math.min(nextTime, intervalEndTime) - t;
            logP += -(coalProp + migProp[deme]) * dt;
            t += dt;
            if (t > intervalEndTime)
                break;
            if (changeIdx < mtNode.getChangeCount()) {
                // Migration
                int toDeme = mtNode.getChangeType(changeIdx);
                logP += Math.log(migModel.getBackwardRate(deme, toDeme));
                deme = toDeme;
            } else {
                // Coalescence
                logP += Math.log(1.0 / migModelSC.getPopSize(deme));
                return logP;
            }
        }
    }
    double t = newRootHeight;
    int sisChangeIdx = 0;
    int demeSis = nodeSis.getNodeType();
    while (sisChangeIdx < nodeSis.getChangeCount() && nodeSis.getChangeTime(sisChangeIdx) < newRootHeight) {
        demeSis = nodeSis.getChangeType(sisChangeIdx);
        sisChangeIdx += 1;
    }
    while (true) {
        // Calculate propensities
        double coalProp;
        if (deme == demeSis)
            coalProp = 1.0 / migModelSC.getPopSize(deme);
        else
            coalProp = 0.0;
    }
// return logP;
}

Example 12 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class MultiTypeTreeScale method proposal.

@Override
public double proposal() {
    // Choose scale factor:
    double u = Randomizer.nextDouble();
    double f = u * scaleFactorInput.get() + (1.0 - u) / scaleFactorInput.get();
    // Keep track of Hastings ratio:
    double logf = Math.log(f);
    double logHR = -2 * logf;
    // Scale colour change times on external branches:
    if (!useOldTreeScalerInput.get()) {
        for (Node leaf : mtTree.getExternalNodes()) {
            double lold = leaf.getParent().getHeight() - leaf.getHeight();
            double lnew = f * leaf.getParent().getHeight() - leaf.getHeight();
            for (int c = 0; c < ((MultiTypeNode) leaf).getChangeCount(); c++) {
                double oldTime = ((MultiTypeNode) leaf).getChangeTime(c);
                double newTime = leaf.getHeight() + (oldTime - leaf.getHeight()) * lnew / lold;
                ((MultiTypeNode) leaf).setChangeTime(c, newTime);
            }
            logHR += ((MultiTypeNode) leaf).getChangeCount() * Math.log(lnew / lold);
        }
    } else {
        for (Node leaf : mtTree.getExternalNodes()) {
            for (int c = 0; c < ((MultiTypeNode) leaf).getChangeCount(); c++) {
                double oldTime = ((MultiTypeNode) leaf).getChangeTime(c);
                ((MultiTypeNode) leaf).setChangeTime(c, f * oldTime);
                logHR += logf;
            }
        }
    }
    // Scale internal node heights and colour change times:
    for (Node node : mtTree.getInternalNodes()) {
        node.setHeight(node.getHeight() * f);
        logHR += logf;
        for (int c = 0; c < ((MultiTypeNode) node).getChangeCount(); c++) {
            double oldTime = ((MultiTypeNode) node).getChangeTime(c);
            ((MultiTypeNode) node).setChangeTime(c, f * oldTime);
            logHR += logf;
        }
    }
    // Reject invalid tree scalings:
    if (f < 1.0) {
        for (Node leaf : mtTree.getExternalNodes()) {
            if (leaf.getParent().getHeight() < leaf.getHeight())
                return Double.NEGATIVE_INFINITY;
            if (((MultiTypeNode) leaf).getChangeCount() > 0 && ((MultiTypeNode) leaf).getChangeTime(0) < leaf.getHeight())
                if (!useOldTreeScalerInput.get())
                    throw new IllegalStateException("Scaled colour change time " + "has dipped below age of leaf - this should never " + "happen!");
                else
                    return Double.NEGATIVE_INFINITY;
        }
    }
    // Scale parameters:
    for (int pidx = 0; pidx < parametersInput.get().size(); pidx++) {
        RealParameter param = parametersInput.get().get(pidx);
        for (int i = 0; i < param.getDimension(); i++) {
            if (!indicatorsUsed || indicatorsInput.get().get(pidx).getValue(i)) {
                double oldValue = param.getValue(i);
                double newValue = oldValue * f;
                if (newValue < param.getLower() || newValue > param.getUpper())
                    return Double.NEGATIVE_INFINITY;
                param.setValue(i, newValue);
                logHR += logf;
            }
        }
    }
    // Scale parameters inversely:
    for (int pidx = 0; pidx < parametersInverseInput.get().size(); pidx++) {
        RealParameter param = parametersInverseInput.get().get(pidx);
        for (int i = 0; i < param.getDimension(); i++) {
            if (!indicatorsInverseUsed || indicatorsInverseInput.get().get(pidx).getValue(i)) {
                double oldValue = param.getValue(i);
                double newValue = oldValue / f;
                if (newValue < param.getLower() || newValue > param.getUpper())
                    return Double.NEGATIVE_INFINITY;
                param.setValue(i, newValue);
                logHR -= logf;
            }
        }
    }
    // Return Hastings ratio:
    return logHR;
}

Example 13 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class NodeRetype method proposal.

@Override
public double proposal() {
    double logHR = 0.0;
    // Select node:
    Node node = mtTree.getNode(mtTree.getLeafNodeCount() + Randomizer.nextInt(mtTree.getInternalNodeCount()));
    // Record probability of current types along attached branches:
    if (!node.isRoot())
        logHR += getBranchTypeProb(node);
    logHR += getBranchTypeProb(node.getLeft()) + getBranchTypeProb(node.getRight());
    // Select new node type:
    ((MultiTypeNode) node).setNodeType(Randomizer.nextInt(migModel.getNTypes()));
    // Retype attached branches:
    try {
        if (!node.isRoot())
            logHR -= retypeBranch(node);
        logHR -= retypeBranch(node.getLeft()) + retypeBranch(node.getRight());
    } catch (NoValidPathException e) {
        return Double.NEGATIVE_INFINITY;
    }
    return logHR;
}

Example 14 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class StructuredCoalescentTreeDensity method updateEventSequence.

/**
 * Determines the sequence of migration, coalescence and sampling events
 * which make up the coloured tree.
 */
protected void updateEventSequence() {
    // Clean up previous list:
    eventList.clear();
    lineageCountList.clear();
    Node rootNode = mtTree.getRoot();
    // Initialise map of active nodes to active change indices:
    Map<Node, Integer> changeIdx = new HashMap<>();
    changeIdx.put(rootNode, -1);
    // Initialise lineage count per colour array:
    Integer[] lineageCount = new Integer[migrationModel.getNTypes()];
    for (int c = 0; c < migrationModel.getNTypes(); c++) if (c == ((MultiTypeNode) rootNode).getNodeType())
        lineageCount[c] = 1;
    else
        lineageCount[c] = 0;
    // Calculate event sequence:
    while (!changeIdx.isEmpty()) {
        SCEvent nextEvent = new SCEvent();
        nextEvent.time = Double.NEGATIVE_INFINITY;
        // Initial assignment not significant
        nextEvent.node = rootNode;
        // Determine next event
        for (Node node : changeIdx.keySet()) if (changeIdx.get(node) < 0) {
            if (node.isLeaf()) {
                // Next event is a sample
                if (node.getHeight() > nextEvent.time) {
                    nextEvent.time = node.getHeight();
                    nextEvent.kind = SCEventKind.SAMPLE;
                    nextEvent.type = ((MultiTypeNode) node).getNodeType();
                    nextEvent.node = node;
                }
            } else {
                // Next event is a coalescence
                if (node.getHeight() > nextEvent.time) {
                    nextEvent.time = node.getHeight();
                    nextEvent.kind = SCEventKind.COALESCE;
                    nextEvent.type = ((MultiTypeNode) node).getNodeType();
                    nextEvent.node = node;
                }
            }
        } else {
            // Next event is a migration
            double thisChangeTime = ((MultiTypeNode) node).getChangeTime(changeIdx.get(node));
            if (thisChangeTime > nextEvent.time) {
                nextEvent.time = thisChangeTime;
                nextEvent.kind = SCEventKind.MIGRATE;
                nextEvent.destType = ((MultiTypeNode) node).getChangeType(changeIdx.get(node));
                if (changeIdx.get(node) > 0)
                    nextEvent.type = ((MultiTypeNode) node).getChangeType(changeIdx.get(node) - 1);
                else
                    nextEvent.type = ((MultiTypeNode) node).getNodeType();
                nextEvent.node = node;
            }
        }
        // Update active node list (changeIdx) and lineage count appropriately:
        switch(nextEvent.kind) {
            case COALESCE:
                Node leftChild = nextEvent.node.getLeft();
                Node rightChild = nextEvent.node.getRight();
                changeIdx.remove(nextEvent.node);
                changeIdx.put(leftChild, ((MultiTypeNode) leftChild).getChangeCount() - 1);
                changeIdx.put(rightChild, ((MultiTypeNode) rightChild).getChangeCount() - 1);
                lineageCount[nextEvent.type]++;
                break;
            case SAMPLE:
                changeIdx.remove(nextEvent.node);
                lineageCount[nextEvent.type]--;
                break;
            case MIGRATE:
                lineageCount[nextEvent.destType]--;
                lineageCount[nextEvent.type]++;
                int oldIdx = changeIdx.get(nextEvent.node);
                changeIdx.put(nextEvent.node, oldIdx - 1);
                break;
        }
        // Add event to list:
        eventList.add(nextEvent);
        lineageCountList.add(Arrays.copyOf(lineageCount, lineageCount.length));
    }
    // Reverse event and lineage count lists (order them from tips to root):
    eventList = Lists.reverse(eventList);
    lineageCountList = Lists.reverse(lineageCountList);
}

Example 15 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class MultiTypeUniform method proposal.

/**
 * Change the node height and return the hastings ratio.
 *
 * @return log of Hastings Ratio
 */
@Override
public double proposal() {
    // Randomly select event on tree:
    int event = Randomizer.nextInt(mtTree.getInternalNodeCount() + mtTree.getTotalNumberOfChanges());
    MultiTypeNode node = null;
    int changeIdx = -1;
    if (event < mtTree.getInternalNodeCount()) {
        node = (MultiTypeNode) mtTree.getNode(mtTree.getLeafNodeCount() + event);
        if (!includeRootInput.get() && node.isRoot())
            return Double.NEGATIVE_INFINITY;
    } else {
        event -= mtTree.getInternalNodeCount();
        for (Node thisNode : mtTree.getNodesAsArray()) {
            if (thisNode.isRoot())
                continue;
            if (event < ((MultiTypeNode) thisNode).getChangeCount()) {
                node = (MultiTypeNode) thisNode;
                changeIdx = event;
                break;
            }
            event -= ((MultiTypeNode) thisNode).getChangeCount();
        }
    }
    if (node == null)
        throw new IllegalStateException("Event selection loop fell through!");
    if (changeIdx == -1) {
        if (node.isRoot()) {
            // Scale distance root and closest event
            double tmin = Math.max(((MultiTypeNode) node.getLeft()).getFinalChangeTime(), ((MultiTypeNode) node.getRight()).getFinalChangeTime());
            double u = Randomizer.nextDouble();
            double f = u * rootScaleFactorInput.get() + (1 - u) / rootScaleFactorInput.get();
            double tnew = tmin + f * (node.getHeight() - tmin);
            node.setHeight(tnew);
            return -Math.log(f);
        } else {
            // Reposition node randomly between closest events
            double tmin = Math.max(((MultiTypeNode) node.getLeft()).getFinalChangeTime(), ((MultiTypeNode) node.getRight()).getFinalChangeTime());
            double tmax = node.getChangeCount() > 0 ? node.getChangeTime(0) : node.getParent().getHeight();
            double u = Randomizer.nextDouble();
            double tnew = u * tmin + (1.0 - u) * tmax;
            node.setHeight(tnew);
            return 0.0;
        }
    } else {
        double tmin, tmax;
        if (changeIdx + 1 < node.getChangeCount())
            tmax = node.getChangeTime(changeIdx + 1);
        else
            tmax = node.getParent().getHeight();
        if (changeIdx - 1 < 0)
            tmin = node.getHeight();
        else
            tmin = node.getChangeTime(changeIdx - 1);
        double u = Randomizer.nextDouble();
        double tnew = u * tmin + (1 - u) * tmax;
        node.setChangeTime(changeIdx, tnew);
        return 0.0;
    }
}