Examples with NodeRef dr.evolution.tree.NodeRef used on opensource projects

Example 11 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class AbstractARGLikelihood method updateNodeAndDescendents.

/**
     * Set update flag for a node and its children
     */
protected void updateNodeAndDescendents(NodeRef node) {
    updateNode[node.getNumber()] = true;
    for (int i = 0; i < treeModel.getChildCount(node); i++) {
        NodeRef child = treeModel.getChild(node, i);
        updateNodeAndDescendents(child);
    }
    likelihoodKnown = false;
}

Example 12 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class simulateClusters method determine_membership_v3.

//traverse down the tree, top down, do calculation
int[] determine_membership_v3(TreeModel treeModel) {
    NodeRef root = treeModel.getRoot();
    int numClusters = 1;
    LinkedList<NodeRef> list = new LinkedList<NodeRef>();
    list.addFirst(root);
    int[] membership = new int[treeModel.getNodeCount()];
    for (int i = 0; i < treeModel.getNodeCount(); i++) {
        membership[i] = -1;
    }
    //root always given the first cluster
    membership[root.getNumber()] = 0;
    while (!list.isEmpty()) {
        //do things with the current object
        NodeRef curElement = list.pop();
        //String content = "node #" + curElement.getNumber() +", taxon=" + treeModel.getNodeTaxon(curElement) + " and parent is = " ;
        String content = "node #" + curElement.getNumber() + ", taxon= ";
        if (treeModel.getNodeTaxon(curElement) == null) {
            content += "internal node\t";
        } else {
            content += treeModel.getNodeTaxon(curElement).getId() + "\t";
        //content += treeModel.getTaxonIndex(treeModel.getNodeTaxon(curElement)) + "\t";
        }
        if (treeModel.getParent(curElement) == null) {
        //content += "no parent";
        } else {
        //content += "parent node#=" + treeModel.getParent(curElement).getNumber();
        }
        //cluster assignment:
        if (!treeModel.isRoot(curElement)) {
            if ((int) indicators.getParameterValue(curElement.getNumber()) == 1) {
                numClusters++;
                membership[curElement.getNumber()] = numClusters - 1;
                //parent cluster label
                int parentClusterLabel = membership[treeModel.getParent(curElement).getNumber()];
                //assign distCluster
                for (int i = 0; i < (numClusters - 1); i++) {
                    distCluster[numClusters - 1][i] = distCluster[parentClusterLabel][i] + 1;
                    distCluster[i][numClusters - 1] = distCluster[i][parentClusterLabel] + 1;
                }
            } else {
                //inherit from parent's cluster assignment
                membership[curElement.getNumber()] = membership[treeModel.getParent(curElement).getNumber()];
            }
        }
        //is not Root
        content += " cluster = " + membership[curElement.getNumber()];
        for (int childNum = 0; childNum < treeModel.getChildCount(curElement); childNum++) {
            list.addFirst(treeModel.getChild(curElement, childNum));
        }
    }
    return (membership);
}

Example 13 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class BeagleTreeLikelihood method traverse.

/**
     * Traverse the tree calculating partial likelihoods.
     *
     * @param tree           tree
     * @param node           node
     * @param operatorNumber operatorNumber
     * @param flip           flip
     * @return boolean
     */
private boolean traverse(Tree tree, NodeRef node, int[] operatorNumber, boolean flip) {
    boolean update = false;
    int nodeNum = node.getNumber();
    NodeRef parent = tree.getParent(node);
    if (operatorNumber != null) {
        operatorNumber[0] = -1;
    }
    // First update the transition probability matrix(ices) for this branch
    if (parent != null && updateNode[nodeNum]) {
        final double branchRate;
        synchronized (branchRateModel) {
            branchRate = branchRateModel.getBranchRate(tree, node);
        }
        final double parentHeight = tree.getNodeHeight(parent);
        final double nodeHeight = tree.getNodeHeight(node);
        // Get the operational time of the branch
        final double branchLength = branchRate * (parentHeight - nodeHeight);
        if (branchLength < 0.0) {
            throw new RuntimeException("Negative branch length: " + branchLength + " (parent: " + parent + "; height: " + parentHeight + " - child: " + node + "height: " + nodeHeight + ")");
        }
        if (flip) {
            substitutionModelDelegate.flipMatrixBuffer(nodeNum);
        }
        branchUpdateIndices[branchUpdateCount] = nodeNum;
        branchLengths[branchUpdateCount] = branchLength;
        branchUpdateCount++;
        update = true;
    }
    // If the node is internal, update the partial likelihoods.
    if (!tree.isExternal(node)) {
        // Traverse down the two child nodes
        NodeRef child1 = tree.getChild(node, 0);
        final int[] op1 = { -1 };
        final boolean update1 = traverse(tree, child1, op1, flip);
        NodeRef child2 = tree.getChild(node, 1);
        final int[] op2 = { -1 };
        final boolean update2 = traverse(tree, child2, op2, flip);
        // If either child node was updated then update this node too
        if (update1 || update2) {
            int x = operationCount[operationListCount] * Beagle.OPERATION_TUPLE_SIZE;
            if (flip) {
                // first flip the partialBufferHelper
                partialBufferHelper.flipOffset(nodeNum);
            }
            final int[] operations = this.operations[operationListCount];
            operations[x] = partialBufferHelper.getOffsetIndex(nodeNum);
            if (useScaleFactors) {
                // get the index of this scaling buffer
                int n = nodeNum - tipCount;
                if (recomputeScaleFactors) {
                    // flip the indicator: can take either n or (internalNodeCount + 1) - n
                    scaleBufferHelper.flipOffset(n);
                    // store the index
                    scaleBufferIndices[n] = scaleBufferHelper.getOffsetIndex(n);
                    // Write new scaleFactor
                    operations[x + 1] = scaleBufferIndices[n];
                    operations[x + 2] = Beagle.NONE;
                } else {
                    operations[x + 1] = Beagle.NONE;
                    // Read existing scaleFactor
                    operations[x + 2] = scaleBufferIndices[n];
                }
            } else {
                if (useAutoScaling) {
                    scaleBufferIndices[nodeNum - tipCount] = partialBufferHelper.getOffsetIndex(nodeNum);
                }
                // Not using scaleFactors
                operations[x + 1] = Beagle.NONE;
                operations[x + 2] = Beagle.NONE;
            }
            // source node 1
            operations[x + 3] = partialBufferHelper.getOffsetIndex(child1.getNumber());
            // source matrix 1
            operations[x + 4] = substitutionModelDelegate.getMatrixIndex(child1.getNumber());
            // source node 2
            operations[x + 5] = partialBufferHelper.getOffsetIndex(child2.getNumber());
            // source matrix 2
            operations[x + 6] = substitutionModelDelegate.getMatrixIndex(child2.getNumber());
            operationCount[operationListCount]++;
            update = true;
            if (hasRestrictedPartials) {
                // Test if this set of partials should be restricted
                if (updateRestrictedNodePartials) {
                    // Recompute map
                    computeNodeToRestrictionMap();
                    updateRestrictedNodePartials = false;
                }
                if (partialsMap[nodeNum] != null) {
                }
            }
        }
    }
    return update;
}

Example 14 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class MarkovJumpsBeagleTreeLikelihood method addRegister.

public void addRegister(Parameter addRegisterParameter, MarkovJumpsType type, boolean scaleByTime) {
    if ((type == MarkovJumpsType.COUNTS && addRegisterParameter.getDimension() != stateCount * stateCount) || (type == MarkovJumpsType.REWARDS && addRegisterParameter.getDimension() != stateCount)) {
        throw new RuntimeException("Register parameter of wrong dimension");
    }
    addVariable(addRegisterParameter);
    final String tag = addRegisterParameter.getId();
    for (int i = 0; i < substitutionModelDelegate.getSubstitutionModelCount(); ++i) {
        registerParameter.add(addRegisterParameter);
        MarkovJumpsSubstitutionModel mjModel;
        SubstitutionModel substitutionModel = substitutionModelDelegate.getSubstitutionModel(i);
        if (useUniformization) {
            mjModel = new UniformizedSubstitutionModel(substitutionModel, type, nSimulants);
        } else {
            if (type == MarkovJumpsType.HISTORY) {
                throw new RuntimeException("Can only report complete history using uniformization");
            }
            mjModel = new MarkovJumpsSubstitutionModel(substitutionModel, type);
        }
        markovjumps.add(mjModel);
        branchModelNumber.add(i);
        addModel(mjModel);
        setupRegistration(numRegisters);
        String traitName;
        if (substitutionModelDelegate.getSubstitutionModelCount() == 1) {
            traitName = tag;
        } else {
            traitName = tag + i;
        }
        jumpTag.add(traitName);
        expectedJumps.add(new double[treeModel.getNodeCount()][patternCount]);
        //        storedExpectedJumps.add(new double[treeModel.getNodeCount()][patternCount]);
        boolean[] oldScaleByTime = this.scaleByTime;
        int oldScaleByTimeLength = (oldScaleByTime == null ? 0 : oldScaleByTime.length);
        this.scaleByTime = new boolean[oldScaleByTimeLength + 1];
        if (oldScaleByTimeLength > 0) {
            System.arraycopy(oldScaleByTime, 0, this.scaleByTime, 0, oldScaleByTimeLength);
        }
        this.scaleByTime[oldScaleByTimeLength] = scaleByTime;
        if (type != MarkovJumpsType.HISTORY) {
            TreeTrait.DA da = new TreeTrait.DA() {

                final int registerNumber = numRegisters;

                public String getTraitName() {
                    return tag;
                }

                public Intent getIntent() {
                    return Intent.BRANCH;
                }

                public double[] getTrait(Tree tree, NodeRef node) {
                    return getMarkovJumpsForNodeAndRegister(tree, node, registerNumber);
                }
            };
            treeTraits.addTrait(traitName + "_base", da);
            treeTraits.addTrait(addRegisterParameter.getId(), new TreeTrait.SumAcrossArrayD(new TreeTrait.SumOverTreeDA(da)));
        } else {
            if (histories == null) {
                histories = new String[treeModel.getNodeCount()][patternCount];
            } else {
                throw new RuntimeException("Only one complete history per markovJumpTreeLikelihood is allowed");
            }
            if (nSimulants > 1) {
                throw new RuntimeException("Only one simulant allowed when saving complete history");
            }
            // Add total number of changes over tree trait
            TreeTrait da = new TreeTrait.DA() {

                final int registerNumber = numRegisters;

                public String getTraitName() {
                    return tag;
                }

                public Intent getIntent() {
                    return Intent.BRANCH;
                }

                public double[] getTrait(Tree tree, NodeRef node) {
                    return getMarkovJumpsForNodeAndRegister(tree, node, registerNumber);
                }
            };
            treeTraits.addTrait(addRegisterParameter.getId(), new TreeTrait.SumOverTreeDA(da));
            // Record the complete history for this register
            historyRegisterNumber = numRegisters;
            ((UniformizedSubstitutionModel) mjModel).setSaveCompleteHistory(true);
            if (useCompactHistory && logHistory) {
                treeTraits.addTrait(ALL_HISTORY, new TreeTrait.SA() {

                    public String getTraitName() {
                        return ALL_HISTORY;
                    }

                    public Intent getIntent() {
                        return Intent.BRANCH;
                    }

                    public boolean getFormatAsArray() {
                        return true;
                    }

                    public String[] getTrait(Tree tree, NodeRef node) {
                        List<String> events = new ArrayList<String>();
                        for (int i = 0; i < patternCount; i++) {
                            String eventString = getHistoryForNode(tree, node, i);
                            if (eventString != null && eventString.compareTo("{}") != 0) {
                                eventString = eventString.substring(1, eventString.length() - 1);
                                if (eventString.contains("},{")) {
                                    // There are multiple events
                                    String[] elements = eventString.split("(?<=\\}),(?=\\{)");
                                    for (String e : elements) {
                                        events.add(e);
                                    }
                                } else {
                                    events.add(eventString);
                                }
                            }
                        }
                        String[] array = new String[events.size()];
                        events.toArray(array);
                        return array;
                    }

                    public boolean getLoggable() {
                        return true;
                    }
                });
            }
            for (int site = 0; site < patternCount; ++site) {
                final String anonName = (patternCount == 1) ? HISTORY : HISTORY + "_" + (site + 1);
                final int anonSite = site;
                treeTraits.addTrait(anonName, new TreeTrait.S() {

                    public String getTraitName() {
                        return anonName;
                    }

                    public Intent getIntent() {
                        return Intent.BRANCH;
                    }

                    public String getTrait(Tree tree, NodeRef node) {
                        String history = getHistoryForNode(tree, node, anonSite);
                        // Return null if empty
                        return (history.compareTo("{}") != 0) ? history : null;
                    }

                    public boolean getLoggable() {
                        return logHistory && !useCompactHistory;
                    }
                });
            }
        }
        numRegisters++;
    }
// End of loop over branch models
}

Example 15 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class SimulationTreeTraversal method traversePreOrder.

/**
     * Traverse the tree in pre order.
     *
     * @param tree tree
     * @param node node
     * @return boolean
     */
private boolean traversePreOrder(final Tree tree, final NodeRef node, final NodeRef parent, final NodeRef sibling) {
    boolean update = (parent == null) ? true : false;
    final int nodeNum = node.getNumber();
    // First update the transition probability matrix(ices) for this branch
    if (parent != null && updateNode[nodeNum]) {
        // @todo - at the moment a matrix is updated even if a branch length doesn't change
        branchNodeOperations.add(new DataLikelihoodDelegate.BranchNodeOperation(nodeNum, parent.getNumber(), computeBranchLength(tree, node)));
        nodeOperations.add(new ProcessOnTreeDelegate.NodeOperation(parent.getNumber(), nodeNum, sibling.getNumber()));
        update = true;
    }
    if (!tree.isExternal(node)) {
        final NodeRef child1 = tree.getChild(node, 0);
        final NodeRef child2 = tree.getChild(node, 1);
        if (update) {
            final boolean update1 = traversePreOrder(tree, child1, node, child2);
            final boolean update2 = traversePreOrder(tree, child2, node, child1);
        }
    }
    return update;
}