Examples with Parameter dr.inference.model.Parameter used on opensource projects

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

Example 2 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);
        }
    }
}

Example 3 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());
    }
}

Example 4 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);
        }
    }
}

Example 5 with Parameter

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

the class TreeClusterSequentialSampling method setVirusLocationAndOffsets.

private void setVirusLocationAndOffsets() {
    //change the mu in the toBin and fromBIn
    //borrow from getLogLikelihood:
    double[] meanYear = new double[binSize];
    double[] groupCount = new double[binSize];
    for (int i = 0; i < numdata; i++) {
        int label = (int) clusterLabels.getParameterValue(i);
        double year = 0;
        if (virusOffsetsParameter != null) {
            //	System.out.print("virus Offeset Parameter present"+ ": ");
            //	System.out.print( virusOffsetsParameter.getParameterValue(i) + " ");
            //	System.out.print(" drift= " + drift + " ");
            //just want year[i]
            year = virusOffsetsParameter.getParameterValue(i);
        //make sure that it is equivalent to double offset  = year[virusIndex] - firstYear;
        } else {
            System.out.println("virus Offeset Parameter NOT present. We expect one though. Something is wrong.");
        }
        meanYear[label] = meanYear[label] + year;
        groupCount[label] = groupCount[label] + 1;
    }
    for (int i = 0; i < binSize; i++) {
        if (groupCount[i] > 0) {
            meanYear[i] = meanYear[i] / groupCount[i];
        }
    //System.out.println(meanYear[i]);
    }
    mu0_offset = new double[binSize];
    //now, change the mu..
    for (int i = 0; i < binSize; i++) {
        //System.out.println(meanYear[i]*beta);
        mu0_offset[i] = meanYear[i];
    //System.out.println("group " + i + "\t" + mu0_offset[i]);
    }
    //virus in the same cluster has the same position
    for (int i = 0; i < numdata; i++) {
        int label = (int) clusterLabels.getParameterValue(i);
        Parameter vLoc = virusLocations.getParameter(i);
        //setting the virus locs to be equal to the corresponding mu
        double muValue = mu.getParameter(label).getParameterValue(0);
        vLoc.setParameterValue(0, muValue);
        double muValue2 = mu.getParameter(label).getParameterValue(1);
        vLoc.setParameterValue(1, muValue2);
    //System.out.println("vloc="+ muValue + "," + muValue2);
    }
    for (int i = 0; i < numdata; i++) {
        int label = (int) clusterLabels.getParameterValue(i);
        //if we want to apply the mean year virus cluster offset to the cluster
        if (clusterOffsetsParameter != null) {
            //setting the clusterOffsets to be equal to the mean year of the virus cluster
            // by doing this, the virus changes cluster AND updates the offset simultaneously
            clusterOffsetsParameter.setParameterValue(i, mu0_offset[label]);
        }
    //		System.out.println("mu0_offset[label]=" + mu0_offset[label]);
    //		System.out.println("clusterOffsets " +  i +" now becomes =" + clusterOffsetsParameter.getParameterValue(i) );   			
    }
//    	System.out.println("===The on nodes===");
//    	for(int i=0; i < binSize; i++){	    
//    		if((int) excisionPoints.getParameterValue(i) == 1){
//    			System.out.println("Cluster node " + i + " = " + (int) indicators.getParameterValue(i) + "\tstatus=" + (int) excisionPoints.getParameterValue(i));
//    		}
//    	}
}