Search in sources :

Example 36 with Parameter

use of dr.inference.model.Parameter in project beast-mcmc by beast-dev.

the class NodeHeightsStatisticParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    String name = xo.getAttribute(Statistic.NAME, xo.getId());
    Tree tree = (Tree) xo.getChild(Tree.class);
    Parameter groupSizes = (Parameter) xo.getChild(Parameter.class);
    return new NodeHeightsStatistic(name, tree, groupSizes);
}
Also used : NodeHeightsStatistic(dr.evomodel.tree.NodeHeightsStatistic) Tree(dr.evolution.tree.Tree) Parameter(dr.inference.model.Parameter)

Example 37 with Parameter

use of dr.inference.model.Parameter in project beast-mcmc by beast-dev.

the class ARGAddRemoveEventOperator method AddOperation.

private double AddOperation() {
    double logHastings = 0;
    double treeHeight = arg.getNodeHeight(arg.getRoot());
    double newBifurcationHeight = Double.POSITIVE_INFINITY;
    double newReassortmentHeight = Double.POSITIVE_INFINITY;
    double theta = probBelowRoot / treeHeight;
    while (newBifurcationHeight > treeHeight && newReassortmentHeight > treeHeight) {
        newBifurcationHeight = MathUtils.nextExponential(theta);
        newReassortmentHeight = MathUtils.nextExponential(theta);
    }
    logHastings += theta * (newBifurcationHeight + newReassortmentHeight) - LOG_TWO - 2.0 * Math.log(theta) + Math.log(1 - Math.exp(-2.0 * treeHeight * theta));
    if (newBifurcationHeight < newReassortmentHeight) {
        double temp = newBifurcationHeight;
        newBifurcationHeight = newReassortmentHeight;
        newReassortmentHeight = temp;
    }
    //2. Find the possible re-assortment and bifurcation points.
    ArrayList<NodeRef> potentialBifurcationChildren = new ArrayList<NodeRef>();
    ArrayList<NodeRef> potentialReassortmentChildren = new ArrayList<NodeRef>();
    int totalPotentialBifurcationChildren = findPotentialAttachmentPoints(newBifurcationHeight, potentialBifurcationChildren);
    int totalPotentialReassortmentChildren = findPotentialAttachmentPoints(newReassortmentHeight, potentialReassortmentChildren);
    assert totalPotentialBifurcationChildren > 0;
    assert totalPotentialReassortmentChildren > 0;
    logHastings += Math.log((double) potentialBifurcationChildren.size() * potentialReassortmentChildren.size());
    //3.  Choose your re-assortment location.
    Node reassortChild = (Node) potentialReassortmentChildren.get(MathUtils.nextInt(totalPotentialReassortmentChildren));
    Node reassortLeftParent = reassortChild.leftParent;
    Node reassortRightParent = reassortChild.rightParent;
    Node reassortSplitParent = reassortChild.leftParent;
    boolean splitReassortLeftParent = true;
    if (!reassortChild.bifurcation) {
        boolean[] tester = { arg.getNodeHeight(reassortLeftParent) > newReassortmentHeight, arg.getNodeHeight(reassortRightParent) > newReassortmentHeight };
        if (tester[0] && tester[1]) {
            if (MathUtils.nextBoolean()) {
                splitReassortLeftParent = false;
                reassortSplitParent = reassortRightParent;
            }
        //        		 logHastings += LOG_TWO;
        } else if (tester[0]) {
        //nothing to do, setup above
        } else if (tester[1]) {
            splitReassortLeftParent = false;
            reassortSplitParent = reassortRightParent;
        } else {
            assert false;
        }
    }
    //        if (recParentL != recParentR) {
    //            boolean[] tester = {arg.getNodeHeight(recParentL) > newReassortmentHeight,
    //                    arg.getNodeHeight(recParentR) > newReassortmentHeight};
    //
    //            if (tester[0] && tester[1]) {
    //                if (MathUtils.nextBoolean()) {
    //                    recParent = recParentR;
    //                }
    //
    //                logHastings += LOG_TWO;
    //            } else if (tester[0]) {
    //                recParent = recParentL;
    //            } else {
    //                recParent = recParentR;
    //            }
    //        }
    //4. Choose your bifurcation location.
    Node bifurcationChild = (Node) potentialBifurcationChildren.get(MathUtils.nextInt(potentialBifurcationChildren.size()));
    Node bifurcationLeftParent = bifurcationChild.leftParent;
    Node bifurcationRightParent = bifurcationChild.rightParent;
    Node bifurcationSplitParent = bifurcationLeftParent;
    boolean splitBifurcationLeftParent = true;
    if (!bifurcationChild.bifurcation) {
        boolean[] tester = { arg.getNodeHeight(bifurcationLeftParent) > newBifurcationHeight, arg.getNodeHeight(bifurcationRightParent) > newBifurcationHeight };
        if (tester[0] && tester[1]) {
            if (MathUtils.nextBoolean()) {
                splitBifurcationLeftParent = false;
                bifurcationSplitParent = bifurcationRightParent;
            }
        //    		  logHastings += LOG_TWO;
        } else if (tester[0]) {
        //nothing to do
        } else if (tester[1]) {
            splitBifurcationLeftParent = false;
            bifurcationSplitParent = bifurcationRightParent;
        } else {
            assert false;
        }
    }
    boolean attachNewReassortNewBifurcationThroughLeft = MathUtils.nextBoolean();
    //During the delete step, this guy gets cancelled out!
    logHastings += LOG_TWO;
    //        if (sisParentL != sisParentR) {
    //            boolean[] tester = {arg.getNodeHeight(sisParentL) > newBifurcationHeight,
    //                    arg.getNodeHeight(sisParentR) > newBifurcationHeight};
    //
    //            if (tester[0] && tester[1]) {
    //                if (MathUtils.nextBoolean()) {
    //                    sisParent = sisParentR;
    //                }
    //                logHastings += LOG_TWO;
    //            } else if (tester[0]) {
    //                sisParent = sisParentL;
    //            } else {
    //                sisParent = sisParentR;
    //            }
    //        }
    //5. Make the new nodes.
    //Note: The height stuff is taken care of below.
    Node newReassortment = arg.new Node();
    newReassortment.bifurcation = false;
    double[] reassortmentValues = { 1.0 };
    if (relaxed) {
        reassortmentValues = ratePrior.generateValues();
    }
    //        logHastings += ratePrior.getAddHastingsRatio(reassortmentValues);
    newReassortment.rateParameter = new Parameter.Default(reassortmentValues);
    newReassortment.rateParameter.addBounds(new Parameter.DefaultBounds(Double.POSITIVE_INFINITY, 0, reassortmentValues.length));
    newReassortment.number = arg.getNodeCount() + 1;
    Node newBifurcation = arg.new Node();
    double[] bifurcationValues = { 1.0 };
    if (relaxed) {
        bifurcationValues = ratePrior.generateValues();
        logHastings += ratePrior.getAddHastingsRatio(bifurcationValues);
    }
    newBifurcation.rateParameter = new Parameter.Default(bifurcationValues);
    newBifurcation.rateParameter.addBounds(new Parameter.DefaultBounds(Double.POSITIVE_INFINITY, 0, bifurcationValues.length));
    newBifurcation.number = arg.getNodeCount();
    //6. Begin editing the tree.
    arg.beginTreeEdit();
    adjustRandomPartitioning();
    if (newBifurcationHeight < treeHeight) {
        if (bifurcationChild == reassortChild) {
            if (bifurcationChild.bifurcation) {
                bifurcationChild.leftParent = bifurcationChild.rightParent = newReassortment;
                newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                newReassortment.leftParent = newReassortment.rightParent = newBifurcation;
                newBifurcation.leftChild = newBifurcation.rightChild = newReassortment;
                newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                if (bifurcationSplitParent.bifurcation) {
                    if (bifurcationSplitParent.leftChild == bifurcationChild) {
                        bifurcationSplitParent.leftChild = newBifurcation;
                    } else {
                        bifurcationSplitParent.rightChild = newBifurcation;
                    }
                } else {
                    bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                }
                logHastings -= LOG_TWO;
            } else {
                if (splitBifurcationLeftParent && splitReassortLeftParent) {
                    bifurcationChild.leftParent = newReassortment;
                    newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                    newReassortment.leftParent = newReassortment.rightParent = newBifurcation;
                    newBifurcation.leftChild = newBifurcation.rightChild = newReassortment;
                    newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                    if (bifurcationSplitParent.bifurcation) {
                        if (bifurcationSplitParent.leftChild == bifurcationChild) {
                            bifurcationSplitParent.leftChild = newBifurcation;
                        } else {
                            bifurcationSplitParent.rightChild = newBifurcation;
                        }
                    } else {
                        bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                    }
                } else if (splitBifurcationLeftParent) {
                    //bifurcation on left, reassortment on right
                    bifurcationChild.leftParent = newBifurcation;
                    bifurcationChild.rightParent = newReassortment;
                    newBifurcation.leftChild = bifurcationChild;
                    newBifurcation.rightChild = newReassortment;
                    newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                    if (bifurcationSplitParent.bifurcation) {
                        if (bifurcationSplitParent.leftChild == bifurcationChild) {
                            bifurcationSplitParent.leftChild = newBifurcation;
                        } else {
                            bifurcationSplitParent.rightChild = newBifurcation;
                        }
                    } else {
                        bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                    }
                    newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                    if (attachNewReassortNewBifurcationThroughLeft) {
                        newReassortment.leftParent = newBifurcation;
                        newReassortment.rightParent = reassortSplitParent;
                    } else {
                        newReassortment.rightParent = newBifurcation;
                        newReassortment.leftParent = reassortSplitParent;
                    }
                    if (reassortSplitParent.bifurcation) {
                        if (reassortSplitParent.leftChild == reassortChild) {
                            reassortSplitParent.leftChild = newReassortment;
                        } else {
                            reassortSplitParent.rightChild = newReassortment;
                        }
                    } else {
                        reassortSplitParent.leftChild = reassortSplitParent.rightChild = newReassortment;
                    }
                } else if (splitReassortLeftParent) {
                    //bifurcation on right, reassortment on left
                    bifurcationChild.rightParent = newBifurcation;
                    bifurcationChild.leftParent = newReassortment;
                    newBifurcation.leftChild = bifurcationChild;
                    newBifurcation.rightChild = newReassortment;
                    newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                    if (bifurcationSplitParent.bifurcation) {
                        if (bifurcationSplitParent.leftChild == bifurcationChild) {
                            bifurcationSplitParent.leftChild = newBifurcation;
                        } else {
                            bifurcationSplitParent.rightChild = newBifurcation;
                        }
                    } else {
                        bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                    }
                    newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                    if (attachNewReassortNewBifurcationThroughLeft) {
                        newReassortment.leftParent = newBifurcation;
                        newReassortment.rightParent = reassortSplitParent;
                    } else {
                        newReassortment.rightParent = newBifurcation;
                        newReassortment.leftParent = reassortSplitParent;
                    }
                    if (reassortSplitParent.bifurcation) {
                        if (reassortSplitParent.leftChild == reassortChild) {
                            reassortSplitParent.leftChild = newReassortment;
                        } else {
                            reassortSplitParent.rightChild = newReassortment;
                        }
                    } else {
                        reassortSplitParent.leftChild = reassortSplitParent.rightChild = newReassortment;
                    }
                } else {
                    bifurcationChild.rightParent = newReassortment;
                    newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                    newReassortment.leftParent = newReassortment.rightParent = newBifurcation;
                    newBifurcation.leftChild = newBifurcation.rightChild = newReassortment;
                    newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                    if (bifurcationSplitParent.bifurcation) {
                        if (bifurcationSplitParent.leftChild == bifurcationChild) {
                            bifurcationSplitParent.leftChild = newBifurcation;
                        } else {
                            bifurcationSplitParent.rightChild = newBifurcation;
                        }
                    } else {
                        bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                    }
                    logHastings -= LOG_TWO;
                }
            }
        } else {
            newReassortment.leftChild = newReassortment.rightChild = reassortChild;
            newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
            newBifurcation.leftChild = newReassortment;
            newBifurcation.rightChild = bifurcationChild;
            if (attachNewReassortNewBifurcationThroughLeft) {
                newReassortment.leftParent = newBifurcation;
                newReassortment.rightParent = reassortSplitParent;
            } else {
                newReassortment.rightParent = newBifurcation;
                newReassortment.leftParent = reassortSplitParent;
            }
            if (reassortChild.bifurcation) {
                reassortChild.leftParent = reassortChild.rightParent = newReassortment;
            } else if (splitReassortLeftParent) {
                reassortChild.leftParent = newReassortment;
            } else {
                reassortChild.rightParent = newReassortment;
            }
            if (reassortSplitParent.bifurcation) {
                if (reassortSplitParent.leftChild == reassortChild) {
                    reassortSplitParent.leftChild = newReassortment;
                } else {
                    reassortSplitParent.rightChild = newReassortment;
                }
            } else {
                reassortSplitParent.leftChild = reassortSplitParent.rightChild = newReassortment;
            }
            if (bifurcationChild.bifurcation) {
                bifurcationChild.leftParent = bifurcationChild.rightParent = newBifurcation;
            } else if (splitBifurcationLeftParent) {
                bifurcationChild.leftParent = newBifurcation;
            } else {
                bifurcationChild.rightParent = newBifurcation;
            }
            if (bifurcationSplitParent.bifurcation) {
                if (bifurcationSplitParent.leftChild == bifurcationChild) {
                    bifurcationSplitParent.leftChild = newBifurcation;
                } else {
                    bifurcationSplitParent.rightChild = newBifurcation;
                }
            } else {
                bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
            }
        }
        Parameter partition = new Parameter.Default(arg.getNumberOfPartitions());
        drawRandomPartitioning(partition);
        newReassortment.partitioning = partition;
        newBifurcation.heightParameter = new Parameter.Default(newBifurcationHeight);
        newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
        newBifurcation.setupHeightBounds();
        newReassortment.setupHeightBounds();
        arg.expandARG(newBifurcation, newReassortment, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
        //                     nodeRates);
        assert nodeCheck() : arg.toARGSummary();
    } else {
        assert newReassortmentHeight < treeHeight;
        //New bifurcation takes the place of the old root.
        //Much easier to program.
        newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
        newReassortment.setupHeightBounds();
        bifurcationChild = newBifurcation;
        if (arg.isRoot(reassortSplitParent))
            reassortSplitParent = newBifurcation;
        Node root = (Node) arg.getRoot();
        Node rootLeftChild = root.leftChild;
        Node rootRightChild = root.rightChild;
        arg.singleRemoveChild(root, rootLeftChild);
        arg.singleRemoveChild(root, rootRightChild);
        arg.singleAddChild(newBifurcation, rootLeftChild);
        arg.singleAddChild(newBifurcation, rootRightChild);
        if (reassortSplitParent.isBifurcation())
            arg.singleRemoveChild(reassortSplitParent, reassortChild);
        else
            arg.doubleRemoveChild(reassortSplitParent, reassortChild);
        arg.doubleAddChild(newReassortment, reassortChild);
        arg.singleAddChild(root, newBifurcation);
        Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
        drawRandomPartitioning(partitioning);
        arg.addChildAsRecombinant(root, reassortSplitParent, newReassortment, partitioning);
        if (attachNewReassortNewBifurcationThroughLeft) {
            newReassortment.leftParent = root;
            newReassortment.rightParent = reassortSplitParent;
        } else {
            newReassortment.leftParent = reassortSplitParent;
            newReassortment.rightParent = root;
        }
        newBifurcation.heightParameter = new Parameter.Default(root.getHeight());
        newBifurcation.setupHeightBounds();
        root.heightParameter.setParameterValue(0, newBifurcationHeight);
        arg.expandARG(newBifurcation, newReassortment, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
        assert nodeCheck();
    }
    //6a. This is when we do not create a new root.
    //        if (newBifurcationHeight < treeHeight) {
    //            newBifurcation.heightParameter = new Parameter.Default(newBifurcationHeight);
    //            newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
    //            newBifurcation.setupHeightBounds();
    //            newReassortment.setupHeightBounds();
    //
    //            if (sisParent.bifurcation)
    //                arg.singleRemoveChild(sisParent, bifurcationChild);
    //            else
    //                arg.doubleRemoveChild(sisParent, bifurcationChild);
    //            if (bifurcationChild != reassortChild) {
    //                if (recParent.bifurcation)
    //                    arg.singleRemoveChild(recParent, reassortChild);
    //                else
    //                    arg.doubleRemoveChild(recParent, reassortChild);
    //            }
    //            if (sisParent.bifurcation)
    //                arg.singleAddChild(sisParent, newBifurcation);
    //            else
    //                arg.doubleAddChild(sisParent, newBifurcation);
    //            if (bifurcationChild != reassortChild)
    //                arg.singleAddChild(newBifurcation, bifurcationChild);
    //            arg.doubleAddChild(newReassortment, reassortChild);
    //
    //            Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
    //            drawRandomPartitioning(partitioning);
    //
    //            if (bifurcationChild != reassortChild) {
    //                arg.addChildAsRecombinant(newBifurcation, recParent,
    //                        newReassortment, partitioning);
    //            } else {
    //                arg.addChildAsRecombinant(newBifurcation, newBifurcation,
    //                        newReassortment, partitioning);
    //            }
    //            arg.expandARGWithRecombinant(newBifurcation, newReassortment,
    //                    internalNodeParameters,
    //                    internalAndRootNodeParameters,
    //                    nodeRates);
    //            assert nodeCheck();
    //
    //            //6b. But here we do.
    //        } else if (newReassortmentHeight < treeHeight) {
    //
    //            newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
    //            newReassortment.setupHeightBounds();
    //
    //            bifurcationChild = newBifurcation;
    //            if (arg.isRoot(recParent))
    //                recParent = newBifurcation;
    //
    //
    //            Node root = (Node) arg.getRoot();
    //            Node rootLeftChild = root.leftChild;
    //            Node rootRightChild = root.rightChild;
    //
    //            arg.singleRemoveChild(root, rootLeftChild);
    //            arg.singleRemoveChild(root, rootRightChild);
    //            arg.singleAddChild(newBifurcation, rootLeftChild);
    //            arg.singleAddChild(newBifurcation, rootRightChild);
    //
    //            if (recParent.isBifurcation())
    //                arg.singleRemoveChild(recParent, reassortChild);
    //            else
    //                arg.doubleRemoveChild(recParent, reassortChild);
    //
    //            arg.doubleAddChild(newReassortment, reassortChild);
    //            arg.singleAddChild(root, newBifurcation);
    //
    //            Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
    //            drawRandomPartitioning(partitioning);
    //
    //
    //            arg.addChildAsRecombinant(root, recParent, newReassortment, partitioning);
    //
    //            newBifurcation.heightParameter = new Parameter.Default(root.getHeight());
    //
    //            newBifurcation.setupHeightBounds();
    //            root.heightParameter.setParameterValue(0, newBifurcationHeight);
    //
    //
    //            arg.expandARGWithRecombinant(newBifurcation, newReassortment,
    //                    internalNodeParameters, internalAndRootNodeParameters,
    //                    nodeRates);
    //
    //            assert nodeCheck();
    //
    //        } else {
    //
    //            Node root = (Node) arg.getRoot();
    //            Node rootLeftChild = root.leftChild;
    //            Node rootRightChild = root.rightChild;
    //
    //            arg.singleRemoveChild(root, rootLeftChild);
    //            arg.singleRemoveChild(root, rootRightChild);
    //            arg.singleAddChild(newBifurcation, rootLeftChild);
    //            arg.singleAddChild(newBifurcation, rootRightChild);
    //
    //            arg.doubleAddChild(newReassortment, newBifurcation);
    //            arg.doubleAddChild(root, newReassortment);
    //
    //            Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
    //            drawRandomPartitioning(partitioning);
    //
    //            newReassortment.partitioning = partitioning;
    //
    //            newBifurcation.heightParameter = new Parameter.Default(arg.getNodeHeight(root));
    //            newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
    //            root.heightParameter.setParameterValueQuietly(0, newBifurcationHeight);
    //
    //            newBifurcation.setupHeightBounds();
    //            newReassortment.setupHeightBounds();
    //
    //            arg.expandARGWithRecombinant(newBifurcation, newReassortment,
    //                    internalNodeParameters, internalAndRootNodeParameters,
    //                    nodeRates);
    //
    //            assert nodeCheck();
    //
    //        }
    arg.pushTreeSizeIncreasedEvent();
    arg.endTreeEdit();
    try {
        arg.checkTreeIsValid();
    } catch (MutableTree.InvalidTreeException ite) {
        throw new RuntimeException(ite.toString() + "\n" + arg.toString() + "\n" + TreeUtils.uniqueNewick(arg, arg.getRoot()));
    }
    assert nodeCheck();
    logHastings -= Math.log((double) findPotentialNodesToRemove(null));
    assert nodeCheck();
    assert !Double.isNaN(logHastings) && !Double.isInfinite(logHastings);
    if (newReassortment.leftParent.bifurcation && newReassortment.rightParent.bifurcation && newReassortment.leftParent != newReassortment.rightParent) {
        logHastings -= LOG_TWO;
    }
    //You're done, return the hastings ratio!
    //		System.out.println(logHastings);
    logHastings += getPartitionAddHastingsRatio(newReassortment.partitioning.getParameterValues());
    return logHastings;
}
Also used : Node(dr.evomodel.arg.ARGModel.Node) ArrayList(java.util.ArrayList) MutableTree(dr.evolution.tree.MutableTree) NodeRef(dr.evolution.tree.NodeRef) CompoundParameter(dr.inference.model.CompoundParameter) Parameter(dr.inference.model.Parameter)

Example 38 with Parameter

use of dr.inference.model.Parameter in project beast-mcmc by beast-dev.

the class ARGAddRemoveEventOperator method adjustRandomPartitioning.

private void adjustRandomPartitioning() {
    if (tossSize < 1) {
        return;
    }
    if (arg.getReassortmentNodeCount() > 0) {
        int total = arg.getReassortmentNodeCount();
        Parameter xyz = arg.getPartitioningParameters().getParameter(MathUtils.nextInt(total));
        if (arg.isRecombinationPartitionType()) {
            adjustRecombinationPartition(xyz);
        } else {
            adjustReassortmentPartition(xyz);
        }
    }
}
Also used : CompoundParameter(dr.inference.model.CompoundParameter) Parameter(dr.inference.model.Parameter)

Example 39 with Parameter

use of dr.inference.model.Parameter in project beast-mcmc by beast-dev.

the class TreeClusterGibbsOperator method setClusterLabels.

private void setClusterLabels(int K_int) {
    int numNodes = treeModel.getNodeCount();
    int[] cutNodes = new int[K_int];
    int cutNum = 0;
    String content = "";
    for (int i = 0; i < binSize; i++) {
        if ((int) excisionPoints.getParameterValue(i) == 1) {
            cutNodes[cutNum] = (int) indicators.getParameterValue(i);
            content += (int) indicators.getParameterValue(i) + ",";
            cutNum++;
        }
    }
    if (cutNum != K_int) {
        System.out.println("cutNum != K_int. we got a problem");
    }
    //    for(int i=0; i < K_int; i++){
    // 	   System.out.println(cutNodes[i]);
    //   }
    //int []membership = determine_membership(treeModel, cutNodes, K_int-1);
    int[] membership = determine_membership(treeModel, cutNodes, K_int);
    double uniqueCode = 0;
    for (int i = 0; i < numNodes; i++) {
        uniqueCode += membership[i] * i;
    }
    //  System.out.println(" sum = " + uniqueCode);
    //   System.out.println("number of nodes = " + treeModel.getNodeCount());
    //  for(int i=0; i < treeModel.getNodeCount(); i++){
    //	   System.out.println(membership[i]);
    //  }
    //System.out.println("Done");
    //  for(int i=0; i < numdata; i++){
    //	   Parameter v = virusLocations.getParameter(i);
    //	   String curName = v.getParameterName();
    //	   System.out.println("i=" + i + " = " + curName);       
    //	}       
    //  for(int j=0; j < numdata; j++){
    //	   System.out.println("j=" + j + " = " + treeModel.getTaxonId(j));
    //  }
    //   Parameter vv = virusLocations.getParameter(0);
    //  String curNamev = vv.getParameterName();
    //  System.out.println(curNamev + " and " +treeModel.getTaxonId(392) );
    //System.out.println(  curNamev.equals(treeModel.getTaxonId(392) )  );
    //System.exit(0);
    // System.out.println("numNodes=" + numNodes);
    // System.exit(0);
    //create dictionary:
    //I suspect this is an expensive operation, so I don't want to do it many times,
    //which is also unnecessary  - MAY have to update whenever a different tree is used.
    int[] membershipToClusterLabelIndexes = new int[numdata];
    for (int i = 0; i < numdata; i++) {
        Parameter v = virusLocations.getParameter(i);
        String curName = v.getParameterName();
        // System.out.println(curName);
        int isFound = 0;
        for (int j = 0; j < numNodes; j++) {
            String treeId = treeModel.getTaxonId(j);
            if (curName.equals(treeId)) {
                //	   System.out.println("  isFound at j=" + j);
                membershipToClusterLabelIndexes[i] = j;
                isFound = 1;
                break;
            }
        }
        if (isFound == 0) {
            System.out.println("not found. Exit now.");
            System.exit(0);
        }
    }
    for (int i = 0; i < numdata; i++) {
        //The assumption that the first nodes being external node corresponding to the cluster labels IS FALSE
        //so I have to search for the matching indexes
        Parameter vloc = virusLocations.getParameter(i);
        //must uncomment out because this sets the new partitioning ... now i am doing code testing.     	   
        clusterLabels.setParameterValue(i, membership[membershipToClusterLabelIndexes[i]]);
    //System.out.println(vloc.getParameterName() + " i="+ i + " membership=" + (int) clusterLabels.getParameterValue(i));
    //  Parameter v = virusLocations.getParameter(i);
    // System.out.println(v.getParameterName());
    }
}
Also used : Parameter(dr.inference.model.Parameter) MatrixParameter(dr.inference.model.MatrixParameter)

Example 40 with Parameter

use of dr.inference.model.Parameter in project beast-mcmc by beast-dev.

the class TreeClusterSequentialSampling method setMembershipToClusterLabelIndexes.

private void setMembershipToClusterLabelIndexes() {
    int numNodes = treeModel.getNodeCount();
    //I suspect this is an expensive operation, so I don't want to do it many times,
    //which is also unnecessary  - MAY have to update whenever a different tree is used.
    membershipToClusterLabelIndexes = new int[numdata];
    for (int i = 0; i < numdata; i++) {
        Parameter v = virusLocations.getParameter(i);
        String curName = v.getParameterName();
        // System.out.println(curName);
        int isFound = 0;
        for (int j = 0; j < numNodes; j++) {
            String treeId = treeModel.getTaxonId(j);
            if (curName.equals(treeId)) {
                //	   System.out.println("  isFound at j=" + j);
                membershipToClusterLabelIndexes[i] = j;
                isFound = 1;
                break;
            }
        }
        if (isFound == 0) {
            System.out.println("not found. Exit now.");
            System.exit(0);
        }
    }
}
Also used : Parameter(dr.inference.model.Parameter) MatrixParameter(dr.inference.model.MatrixParameter)

Aggregations

Parameter (dr.inference.model.Parameter)397 TreeModel (dr.evomodel.tree.TreeModel)62 MatrixParameter (dr.inference.model.MatrixParameter)46 ArrayList (java.util.ArrayList)44 FrequencyModel (dr.oldevomodel.substmodel.FrequencyModel)43 FrequencyModel (dr.evomodel.substmodel.FrequencyModel)41 Units (dr.evolution.util.Units)36 XMLUnits (dr.evoxml.util.XMLUnits)36 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)30 Tree (dr.evolution.tree.Tree)25 DataType (dr.evolution.datatype.DataType)24 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)23 CompoundParameter (dr.inference.model.CompoundParameter)23 GammaSiteModel (dr.oldevomodel.sitemodel.GammaSiteModel)21 SitePatterns (dr.evolution.alignment.SitePatterns)20 HKY (dr.evomodel.substmodel.nucleotide.HKY)17 Likelihood (dr.inference.model.Likelihood)17 HomogeneousBranchModel (dr.evomodel.branchmodel.HomogeneousBranchModel)16 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)16 ParametricDistributionModel (dr.inference.distribution.ParametricDistributionModel)16