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

Example 61 with Parameter

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

the class Tree_Clustering_Shared_Routines method setMembershipTreeToVirusIndexes.

public static int[] setMembershipTreeToVirusIndexes(int numdata, MatrixParameter virusLocations, int numNodes, TreeModel treeModel) {
    //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[] correspondingTreeIndexForVirus = 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);
                correspondingTreeIndexForVirus[i] = j;
                isFound = 1;
                break;
            }
        }
        if (isFound == 0) {
            System.out.println("not found. Exit now.");
            System.exit(0);
        }
    }
    return (correspondingTreeIndexForVirus);
}

Example 62 with Parameter

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

the class TreeClusteringVirusesPrior 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 63 with Parameter

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

the class TreeClusterAlgorithmOperator method Proposal_HotMultistepOnNodeFlipMu.

private double Proposal_HotMultistepOnNodeFlipMu() {
    int rootNum = treeModel.getRoot().getNumber();
    //unlike the old version, self-move isn't allowed.
    //find an on-node
    int originalNode1 = findAnOnNodeRandomly();
    //System.out.print("Try " + originalNode1);
    int[] numStepsFromI_selected = determineTreeNeighborhood(originalNode1, 100000);
    //1. Select an unoccupied site within some steps away from it.	 
    LinkedList<Integer> possibilities1 = new LinkedList<Integer>();
    for (int i = 0; i < numNodes; i++) {
        // System.out.println("#steps from I_selected " + numStepsFromI_selected[i]);
        //make sure no self select
        boolean isIn1 = numStepsFromI_selected[i] <= maxNodeLevel && numStepsFromI_selected[i] != 0 && i != rootNum;
        if (isIn1 && hotNodes[i] == 1) {
            possibilities1.addLast(new Integer(i));
        }
    }
    //end for		
    int numPossibilities1 = possibilities1.size();
    //if there is a single legal configuration to switch to
    if (numPossibilities1 > 0) {
        //choose from possibilities
        int whichMove = (int) (Math.floor(Math.random() * numPossibilities1));
        int site_add1 = possibilities1.get(whichMove).intValue();
        //	 System.out.println(" and select " + site_add1);
        // System.out.println("selected node = " + site_add1 + " that's " + numStepsFromI_selected[site_add1] + " steps from " + originalNode1);
        // the existing indicator is now set to be off
        indicators.setParameterValue(originalNode1, 0);
        //set the new selected index to the new node.
        indicators.setParameterValue(site_add1, 1);
        //Flip mu - so the neighbor that replaces the original node now also inherits the existing node's mu
        //Parameter originalNodeMu = mu.getParameter(originalNode1+1); //offset of 1
        Parameter originalNodeMu = mu.getParameter(originalNode1);
        double[] tmp = originalNodeMu.getParameterValues();
        //Parameter newMu = mu.getParameter(site_add1+1); //offset of 1
        Parameter newMu = mu.getParameter(site_add1);
        double[] tmpNew = newMu.getParameterValues();
        originalNodeMu.setParameterValue(0, tmpNew[0]);
        originalNodeMu.setParameterValue(1, tmpNew[1]);
        newMu.setParameterValue(0, tmp[0]);
        newMu.setParameterValue(1, tmp[1]);
        //backward calculation
        int[] numStepsBackward = determineTreeNeighborhood(site_add1, 100000);
        //1. Select an unoccupied site within some steps away from it.
        LinkedList<Integer> possibilities2 = new LinkedList<Integer>();
        for (int i = 0; i < numNodes; i++) {
            // System.out.println("#steps from I_selected " + numStepsFromI_selected[i]);
            //make sure no self select
            boolean isIn2 = numStepsBackward[i] <= maxNodeLevel && numStepsBackward[i] != 0 && i != rootNum;
            if (isIn2 && hotNodes[i] == 1) {
                possibilities2.addLast(new Integer(i));
            }
        }
        //end for
        int numPossibilities2 = possibilities2.size();
        //	 System.out.println("numPossibilities1=" + numPossibilities1 + " numPossibilities2 = " + numPossibilities2);
        double logHastingRatio = Math.log((1 / (double) numPossibilities2) / (1 / (double) numPossibilities1));
        //System.out.println("logHastingRatio = " + logHastingRatio);
        return logHastingRatio;
    } else {
        return Double.NEGATIVE_INFINITY;
    }
}

Example 64 with Parameter

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

the class TreeClusterAlgorithmOperator method findRestrictedNodeRandomly.

private int findRestrictedNodeRandomly(double neighborhood) {
    double mu1ScaleValue = mu1Scale.getParameterValue(0);
    double mu2ScaleValue = mu2Scale.getParameterValue(0);
    int rootNode = treeModel.getRoot().getNumber();
    int numQualified = 0;
    int[] qualifiedNodes = new int[numNodes];
    for (int i = 0; i < numNodes; i++) {
        if (i != rootNode) {
            Parameter curMu = mu.getParameter(i);
            double coord1 = mu1ScaleValue * curMu.getParameterValue(0);
            double coord2 = mu2ScaleValue * curMu.getParameterValue(1);
            //		double coord1 = curMu.getParameterValue(0);
            //	double coord2 = curMu.getParameterValue(1);
            double dist = Math.sqrt(coord1 * coord1 + coord2 * coord2);
            if (dist < neighborhood) {
                qualifiedNodes[numQualified] = i;
                numQualified++;
            }
        }
    }
    //now draw 
    if (numQualified > 0) {
        int ranSelect = (int) (Math.random() * numQualified);
        int selectedNode = qualifiedNodes[ranSelect];
        return selectedNode;
    }
    // no node qualified, return -1
    return -1;
}

Example 65 with Parameter

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

the class TreeNodeSlide method operateOneNode.

public void operateOneNode(final double factor) {
    //            #print "operate: tree", ut.treerep(t)
    //   if( verbose)  System.out.println("  Mau at start: " + tree.getSimpleTree());
    final int count = tree.getExternalNodeCount();
    assert count == species.nSpecies();
    NodeRef[] order = new NodeRef[2 * count - 1];
    boolean[] swapped = new boolean[count - 1];
    mauCanonical(tree, order, swapped);
    // internal node to change
    // count-1 - number of internal nodes
    int which = MathUtils.nextInt(count - 1);
    //        if( outgroupOnly ) {
    //           if( order[1] == tree.getRoot() ) {
    //               which = 0;
    //           } else if( order[2*count - 3] == tree.getRoot() ) {
    //               which = count - 2;
    //           }
    //        }
    //        if( verbose)  {
    //            System.out.print("order:");
    //            for(int k = 0 ; k < 2*count; k += 2) {
    //                System.out.print(tree.getNodeTaxon(order[k]).getId() + ((k == 2*which) ? "*" : " "));
    //            }
    //            System.out.println();
    //        }
    FixedBitSet left = new FixedBitSet(count);
    FixedBitSet right = new FixedBitSet(count);
    for (int k = 0; k < 2 * which + 1; k += 2) {
        left.set(tree.speciesIndex(order[k]));
    }
    for (int k = 2 * (which + 1); k < 2 * count; k += 2) {
        right.set(tree.speciesIndex(order[k]));
    }
    double newHeight;
    if (factor > 0) {
        newHeight = tree.getNodeHeight(order[2 * which + 1]) * factor;
    } else {
        final double limit = species.speciationUpperBound(left, right);
        newHeight = MathUtils.nextDouble() * limit;
    }
    if (false) {
        double totChange = 0;
        double tot = 0;
        double totChangeUp = 0.0;
        int topChange = 0;
        for (int which1 = 0; which1 < count - 1; which1++) {
            FixedBitSet left1 = new FixedBitSet(count);
            FixedBitSet right1 = new FixedBitSet(count);
            final NodeRef node = order[2 * which1 + 1];
            final double h = tree.getNodeHeight(node);
            for (int k = 0; k < 2 * which1 + 1; k += 2) {
                final NodeRef ref = order[k];
                left1.set(tree.speciesIndex(ref));
            }
            for (int k = 2 * (which1 + 1); k < 2 * count; k += 2) {
                right1.set(tree.speciesIndex(order[k]));
            }
            final double limit = species.speciationUpperBound(left1, right1);
            final NodeRef parent = tree.getParent(node);
            final double ph = parent != null ? tree.getNodeHeight(parent) : limit;
            final double h1 = tree.getNodeHeight(tree.getChild(node, 0));
            final double h2 = tree.getNodeHeight(tree.getChild(node, 1));
            final double chDown = Math.max(h1, h2);
            final double chup = Math.max(limit - ph, 0);
            totChangeUp += chup;
            double change = chDown + chup;
            totChange += change;
            if (change > limit) {
                assert change <= limit;
            }
            tot += limit;
            if (which == which1 && (newHeight < chDown || newHeight > limit - chup)) {
                topChange = 1;
            }
        //               final double h1 = tree.getNodeHeight(order[2 * which1]);
        //               final double h2 = tree.getNodeHeight(order[2 * which1 + 2]);
        //               double lowerLimit = 0;
        //               double upperLimit = limit;
        //               if( h > h1 ) {
        //                   assert tree.getParent(order[2 * which]) == node;
        //                   lowerLimit = h1;
        //               } else {
        //                   assert tree.getParent() == node;
        //                   upperLimit = Math.min(upperLimit, h1);
        //               }
        //               if( h > h2 && h2 > lowerLimit ) {
        //                   assert tree.getParent(order[2 * which+2]) == node;
        //                   lowerLimit = h2;
        //               }
        //
        }
        final double p = totChange / tot;
        System.err.println("top% " + p + " " + totChangeUp / tot + (topChange > 0 ? "++" : ""));
    }
    tree.beginTreeEdit();
    tree.setPreorderIndices(preOrderIndexBefore);
    final NodeRef node = order[2 * which + 1];
    if (false) {
        System.err.println("Before " + node.getNumber() + " " + newHeight);
        ptree(tree, order, preOrderIndexBefore);
        if (newHeight > tree.getNodeHeight(node)) {
            NodeRef parent = tree.getParent(node);
            while (parent != null && tree.getNodeHeight(parent) < newHeight) {
                // swap
                final int pi = getIndexOf(parent, order);
                final int side = pi < which ? 0 : 1;
                // 0 for right, 1 for left
                //final int side = (1+getSide(which, parent, order, swapped))/2;
                swapPops(node, swapped[which] ? 1 - side : side, parent, swapped[pi] ? side : 1 - side, preOrderIndexBefore);
                parent = tree.getParent(parent);
            }
        } else {
            // NodeRef child = tree.getChild(node, 0);
            if (which > 0) {
                NodeRef nb = order[2 * which - 1];
                if (tree.getNodeHeight(node) > tree.getNodeHeight(nb)) {
                    while (nb != node && tree.getNodeHeight(nb) < newHeight) {
                        nb = tree.getParent(nb);
                    }
                    if (nb != node) {
                        final int side = swapped[which] ? 1 : 0;
                        swapPops(node, side, nb, 1 - side, preOrderIndexBefore);
                        NodeRef pnb = tree.getParent(nb);
                        while (pnb != node) {
                            swapPops(nb, 1 - side, pnb, 1 - side, preOrderIndexBefore);
                            nb = pnb;
                            pnb = tree.getParent(nb);
                        }
                    }
                }
            }
            if (which < count - 2) {
                NodeRef nb = order[2 * which + 3];
                if (tree.getNodeHeight(node) > tree.getNodeHeight(nb)) {
                    while (nb != node && tree.getNodeHeight(nb) < newHeight) {
                        nb = tree.getParent(nb);
                    }
                    if (nb != node) {
                        final int side = swapped[which] ? 0 : 1;
                        swapPops(node, side, nb, 1 - side, preOrderIndexBefore);
                        NodeRef pnb = tree.getParent(nb);
                        while (pnb != node) {
                            swapPops(nb, 1 - side, pnb, 1 - side, preOrderIndexBefore);
                            nb = pnb;
                            pnb = tree.getParent(nb);
                        }
                    }
                }
            }
        }
    }
    tree.setNodeHeight(node, newHeight);
    mauReconstruct(tree, order, swapped);
    // restore pre-order of pops -
    {
        tree.setPreorderIndices(preOrderIndexAfter);
        double[] splitPopValues = null;
        for (int k = 0; k < preOrderIndexBefore.length; ++k) {
            final int b = preOrderIndexBefore[k];
            if (b >= 0) {
                final int a = preOrderIndexAfter[k];
                if (a != b) {
                    //if( verbose)  System.out.println("pops: " + a + " <- " + b);
                    final Parameter p1 = tree.sppSplitPopulations;
                    if (splitPopValues == null) {
                        splitPopValues = p1.getParameterValues();
                    }
                    if (tree.constPopulation()) {
                        p1.setParameterValue(count + a, splitPopValues[count + b]);
                    } else {
                        for (int i = 0; i < 2; ++i) {
                            p1.setParameterValue(count + 2 * a + i, splitPopValues[count + 2 * b + i]);
                        }
                    }
                }
            }
        }
    }
    // System.err.println("After");
    //ptree(tree, order, preOrderIndexAfter);
    // if( verbose)  System.out.println("  Mau after: " + tree.getSimpleTree());
    // }
    tree.endTreeEdit();
}