Java Randomizer.nextExponential方法代码示例

import beast.util.Randomizer; //导入方法依赖的package包/类
public Sequence mutateLang(Sequence l, double T) throws Exception {
	Sequence newLang = new Sequence("", l.getData());
	for (int i = 0; i < newLang.getData().length(); i++) {
		int currentTrait = (int) newLang.getData().charAt(i);
		// Mutations are exponentially distributed.
		double t = Randomizer.nextExponential(rate);
		String newSeq;
		while (t < T) {
			currentTrait = Character.getNumericValue(newLang.getData().charAt(i));
			// In binary model, a mutation switches trait.
			// If death: check NoEmptyTrait.
			if (1 - currentTrait == 0) {
				if (noEmptyTraitCheck(newLang)) {
					newSeq = replaceCharAt(newLang.getData(), i, Integer.toString((1 - currentTrait)));
				} else {
					newSeq = newLang.getData();
				}
			} else {
				newSeq = replaceCharAt(newLang.getData(), i, Integer.toString((1 - currentTrait)));
			}
			newLang.dataInput.setValue(newSeq, newLang);
			t += Randomizer.nextExponential(rate);
		}
	}
	return newLang;
}
 

import beast.util.Randomizer; //导入方法依赖的package包/类
static Tree randomYuleTree (int nodes, double l) throws Exception {
	Tree tr = new Tree();
	ArrayList<Node> nodeList = new ArrayList<Node>();
	double t = 0.0;
	int label = 1;
	for (int i = 0; i < nodes; i++) {
		Node n = new Node();
		n.setHeight(t);
		n.setNr(label);
		label++;
		nodeList.add(n);
	}

	while (nodeList.size() > 1) {
		t += Randomizer.nextExponential(nodeList.size()*l);

		int p_1_index = Randomizer.nextInt(nodeList.size());
		Node p1 = nodeList.remove(p_1_index);
		int p_2_index = Randomizer.nextInt(nodeList.size());
		Node p2 = nodeList.remove(p_2_index);

		Node parent = new Node();
		parent.setHeight(t);

		parent.setNr(label);
		label++;
		p1.setParent(parent);
		parent.addChild(p1);
		p2.setParent(parent);
		parent.addChild(p2);

		nodeList.add(parent);
	}
	tr.setRoot(nodeList.get(0));
	return new Tree(tr.toString());

}
 

import beast.util.Randomizer; //导入方法依赖的package包/类
private static Tree randomYuleTree (int nodes, double l) throws Exception {
	Tree tr = new Tree();
	ArrayList<Node> nodeList = new ArrayList<Node>();
	double t = 0.0;
	int label = 1;
	for (int i = 0; i < nodes; i++) {
		Node n = new Node();
		n.setHeight(t);
		n.setNr(label);
		label++;
		nodeList.add(n);
	}

	while (nodeList.size() > 1) {
		t += Randomizer.nextExponential(nodeList.size()*l);

		int p_1_index = Randomizer.nextInt(nodeList.size());
		Node p1 = nodeList.remove(p_1_index);
		int p_2_index = Randomizer.nextInt(nodeList.size());
		Node p2 = nodeList.remove(p_2_index);

		Node parent = new Node();
		parent.setHeight(t);

		parent.setNr(label);
		label++;
		p1.setParent(parent);
		parent.addChild(p1);
		p2.setParent(parent);
		parent.addChild(p2);

		nodeList.add(parent);
	}
	tr.setRoot(nodeList.get(0));
	return new Tree(tr.toString());

}
 

import beast.util.Randomizer; //导入方法依赖的package包/类
/**
 * override this for proposals,
 *
 * @return log of Hastings Ratio, or Double.NEGATIVE_INFINITY if proposal should not be accepted *
 */
@Override
public double proposal() {
    // always operate on the root node
    final Node speciesTreeRoot = speciesTree.getRoot();
    
    // don't bother if root node is a sampled ancestor
    // TODO make it work
    if (speciesTreeRoot.isFake()) return Double.NEGATIVE_INFINITY;

    final double currentRootHeight = speciesTreeRoot.getHeight();
    final double leftChildHeight = speciesTreeRoot.getLeft().getHeight();
    final double rightChildHeight = speciesTreeRoot.getRight().getHeight();

    final MinimumDouble tipwardFreedom = new MinimumDouble();
    final SetMultimap<Integer, Node> connectingNodes = getConnectingNodes(speciesTreeRoot, tipwardFreedom);

    tipwardFreedom.set(currentRootHeight - leftChildHeight);
    tipwardFreedom.set(currentRootHeight - rightChildHeight);
    waitingTime = tipwardFreedom.get() * waitingTimeScale;

    // the youngest age the species tree root node can be (preserving topologies)
    final double uniformShift = Randomizer.nextExponential(lambda) - tipwardFreedom.get();

    speciesTreeRoot.setHeight(currentRootHeight + uniformShift);
    for (Node geneTreeNode: connectingNodes.values()) {
        geneTreeNode.setHeight(geneTreeNode.getHeight() + uniformShift);
    }

    // the log ratio of the density of the proposed over the current species tree root heights
    final double fLogHastingsRatio = lambda * uniformShift;

    return fLogHastingsRatio;
}
 

import beast.util.Randomizer; //导入方法依赖的package包/类
public Tree mutateOverTreeBorrowing(Tree base) throws Exception {
	Double[] events = getEvents(base);
	setSubTreeLanguages(base.getRoot(), (Sequence) base.getRoot().getMetaData("lang"));
	// Get root node.
	ArrayList<Node> aliveNodes = new ArrayList<Node>();
	String[] stringAliveNodes = {};
	int[] traits = {};
	int numberOfLangs = 0;
	Double totalRate = null;
	int idx, ind;
	double[] probs = new double[3];
	for (int i = 0; i < events.length - 1; i++) {
		if (events[i] == 0.0) {
			break;
		}
		aliveNodes = getAliveNodes(base, events[i+1]);
		stringAliveNodes = getSequences(aliveNodes);
		numberOfLangs = stringAliveNodes.length;
		traits = getBirths(stringAliveNodes, numberOfLangs);
		totalRate = totalRate(stringAliveNodes, traits, numberOfLangs);
		Double t = events[i] - Randomizer.nextExponential(totalRate);
		//System.out.println();
		//System.out.println("On branch event: " + (i+1)+ " out of " + (events.length/2) + ". Next event at " + events[i+1]);
		while (t > events[i+1]) {
			//System.out.print("\r"+t);
			probs = BorrowingProbs(stringAliveNodes, totalRate,traits, numberOfLangs);
			Integer choice = Randomizer.randomChoicePDF(probs);
			// Birth.
			if (choice == 0) {
				idx = Randomizer.nextInt(stringAliveNodes.length);
				stringAliveNodes[idx] = stringAliveNodes[idx] + "1";
				stringAliveNodes = addEmptyTrait(stringAliveNodes, idx);
				traits[idx]++;
				// Death.
			} else if (choice == 1) {
				idx = getDeathNode(stringAliveNodes, traits, numberOfLangs);
				// Find random alive trait, and kill it.
				ind = getRandomBirthIndex(stringAliveNodes[idx], traits[idx]);
				if (noEmptyTraitCheck(stringAliveNodes[idx]) && ind > -1) {
					stringAliveNodes[idx] = replaceCharAt(stringAliveNodes[idx], ind, Integer.toString(0));
					traits[idx]--;
				}
				// Borrowing.
			} else if (choice == 2) {
				if (aliveNodes.size() > 1) {
					// Pick two distinct languages at random.
					int[] bN = getBorrowingNodes(stringAliveNodes, traits, numberOfLangs);

					if (localDist(aliveNodes.get(bN[0]), aliveNodes.get(bN[1]))) {
						// Randomly iterate through language and find a 1.
						ind = getRandomBirthIndex(stringAliveNodes[bN[0]], traits[bN[0]]);
						// If recieving language is going 0 -> 1.
						if (ind > 1 && stringAliveNodes[bN[1]].charAt(ind) == '0') {
							traits[bN[1]]++;
						}
						// Give the 1 to the receiving language.
						stringAliveNodes[bN[1]] = replaceCharAt(stringAliveNodes[bN[1]], ind, Integer.toString(1));
					}
				}
			}
			totalRate = totalRate(stringAliveNodes, traits, numberOfLangs);
			t -= Randomizer.nextExponential(totalRate);
		}
		setLangs(aliveNodes, stringAliveNodes);
	}
	return base;
}
 

import beast.util.Randomizer; //导入方法依赖的package包/类
public Tree mutateOverTreeBorrowing(Tree base) throws Exception {
	Double[] events = getEvents(base);
	setSubTreeLanguages(base.getRoot(), (Sequence) base.getRoot().getMetaData("lang"));
	// Get root node.
	ArrayList<Node> aliveNodes = new ArrayList<Node>();
	String[] stringAliveNodes = {};
	int[] traits = {};
	int numberOfLangs = 0;
	// Get first event.
	Double totalRate = null;
	//Double t = treeHeight - Randomizer.nextExponential(totalRate);
	// Variable declarations.
	int idx;
	double[] probs;
	for (int i = 0; i < events.length - 1; i++) {
		if (events[i] == 0.0) {
			break;
		}
		aliveNodes = getAliveNodes(base, events[i+1]);
		stringAliveNodes = getSequences(aliveNodes);
		numberOfLangs = stringAliveNodes.length;
		traits = getBirths(stringAliveNodes, numberOfLangs);
		totalRate = totalRate(stringAliveNodes, traits, numberOfLangs);
		Double t = events[i] - Randomizer.nextExponential(totalRate);
		//System.out.println();
		//System.out.println("On branch event: " + (i+1)+ " out of " + (events.length/2) + ". Next event at " + events[i+1]);
		while (t > events[i+1]) {
			//System.out.print("\r"+t);
			// Return array of event probabilities and pick one.
			probs = BorrowingProbs(stringAliveNodes, totalRate, traits, numberOfLangs);
			Integer choice = Randomizer.randomChoicePDF(probs);
			// Mutate.
			if (choice == 0) {
				// Pick a random node at time t.
				idx = Randomizer.nextInt(stringAliveNodes.length);
				// Pick a random position in language.
				int pos = Randomizer.nextInt(stringAliveNodes[idx].length());
				int currentTrait = Character.getNumericValue(stringAliveNodes[idx].charAt(pos));
				// If death and noEmptyTraitCheck fails.
				if (currentTrait == 1 && (!noEmptyTraitCheck(stringAliveNodes[idx]))) {
					stringAliveNodes[idx] = stringAliveNodes[idx];
				} else {
					stringAliveNodes[idx] = replaceCharAt(stringAliveNodes[idx], pos, Integer.toString((1 - currentTrait)));
				} 
				// Change trait structure
				if (stringAliveNodes[idx].charAt(pos) == '1') {
					traits[idx]++;
				} else {
					traits[idx]--;
				}

				// Borrow.
			} else if (choice == 1) {
				if (aliveNodes.size() > 1) {
					// Pick two distinct languages at random.
					int[] bN = getBorrowingNodes(stringAliveNodes, traits, numberOfLangs);

					if (localDist(aliveNodes.get(bN[0]), aliveNodes.get(bN[1]))) {
						// Randomly iterate through language and find a 1.
						int ind = getRandomBirthIndex(stringAliveNodes[bN[0]], traits[bN[0]]);
						// If recieving language is going 0 -> 1.
						if (ind > -1 && stringAliveNodes[bN[1]].charAt(ind) == '0') {
							traits[bN[1]]++;
						}
						// Give the 1 to the receiving language.
						stringAliveNodes[bN[1]] = replaceCharAt(stringAliveNodes[bN[1]], ind, Integer.toString(1));
					}
				}
			}
			totalRate = totalRate(stringAliveNodes, traits, numberOfLangs);
			t -= Randomizer.nextExponential(totalRate);
		}
		setLangs(aliveNodes, stringAliveNodes);
	}
	return base;
}
 

import beast.util.Randomizer; //导入方法依赖的package包/类
/**
 * override this for proposals,
 * returns log of hastingRatio, or Double.NEGATIVE_INFINITY if proposal should not be accepted *
 */
@Override
public double proposal(Evaluator E) {

    int m = 100;

    RealParameter X = parameterInput.get();

    // Find the density at the current point
    Double gx0 = evaluate(E);

    //	System.err.println("gx0 = " + gx0);
    // Get the 1st element
    Double x0 = X.getValue(0);
    //	System.err.println("x0 = " + x0);

    // Determine the slice level, in log terms.
    double logy = gx0 - Randomizer.nextExponential(1);

    // Find the initial interval to sample from.
    Double[] range = find_slice_boundaries_stepping_out(E, X, logy, windowSize, m);
    Double L = range[0];
    Double R = range[1];

    // Sample from the interval, shrinking it on each rejection
    double x_new = search_interval(E, x0, X, L, R, logy);
    X.setValue(x_new);

    if (n_learning_iterations > 0) {
        n_learning_iterations--;

        totalDelta += Math.abs(x_new - x0);
        totalNumber++;

        double W_predicted = totalDelta / totalNumber * 4.0;
        if (totalNumber > 3) {
            W = 0.95 * W + 0.05 * W_predicted;
            windowSize = W;
        }
        //	    System.err.println("W = " + W);
    }

    return Double.POSITIVE_INFINITY;
}
 

import beast.util.Randomizer; //导入方法依赖的package包/类
/**
 * Use coalescent model to simulate clonal frame.
 */
private void simulateClonalFrame() {

    // Initialize leaf nodes
    List<Node> leafNodes = new ArrayList<>();
    for (int i=0; i<m_taxonset.get().getTaxonCount(); i++) {
        Node leaf = new Node();
        leaf.setNr(i);
        leaf.setID(m_taxonset.get().getTaxonId(i));
                    
        if (hasDateTrait())
            leaf.setHeight(getDateTrait().getValue(leaf.getID()));
        else
            leaf.setHeight(0.0);
        
        leafNodes.add(leaf);
    }
    
    // Create and sort list of inactive nodes
    List<Node> inactiveNodes = new ArrayList<>(leafNodes);
    Collections.sort(inactiveNodes, (Node n1, Node n2) -> {
        if (n1.getHeight()<n2.getHeight())
            return -1;
        
        if (n1.getHeight()>n2.getHeight())
            return 1;
        
        return 0;
    });
    
    List<Node> activeNodes = new ArrayList<>();
    
    double tau = 0.0;
    int nextNr = leafNodes.size();
    while (true) {
        
        // Calculate coalescence propensity
        int k = activeNodes.size();
        double chi = 0.5*k*(k-1);
        
        // Draw scaled coalescent time
        if (chi>0.0)
            tau += Randomizer.nextExponential(chi);
        else
            tau = Double.POSITIVE_INFINITY;
        
        // Convert to real time
        double t = popFunc.getInverseIntensity(tau);
        
        // If new time takes us past next sample time, insert that sample
        if (!inactiveNodes.isEmpty() && t>inactiveNodes.get(0).getHeight()) {
            Node nextActive = inactiveNodes.remove(0);
            activeNodes.add(nextActive);
            tau = popFunc.getIntensity(nextActive.getHeight());
            continue;
        }
        
        // Coalesce random pair of active nodes.
        Node node1 = activeNodes.remove(Randomizer.nextInt(k));
        Node node2 = activeNodes.remove(Randomizer.nextInt(k-1));
        
        Node parent = new Node();
        parent.addChild(node1);
        parent.addChild(node2);
        parent.setHeight(t);
        parent.setNr(nextNr++);
        
        activeNodes.add(parent);
        
        if (inactiveNodes.isEmpty() && activeNodes.size()<2)
            break;
    }
    
    // Remaining active node is root
    setRoot(activeNodes.get(0));
}
 

import beast.util.Randomizer; //导入方法依赖的package包/类
/**
 * Take a recombination with an existing departure point and determine
 * the arrival point by allowing it to coalesce with the clonal frame.
 * 
 * @param conv recombination to modify
 * @return log probability density of coalescent point chosen.
 */
public double coalesceEdge(Conversion conv) {
    double logP = 0.0;
    
    List<CFEventList.Event> events = acg.getCFEvents();
    
    // Locate event immediately below departure point
    int startIdx = 0;
    while (events.get(startIdx+1).getHeight()<conv.getHeight1())
        startIdx += 1;
            
    // Choose edge length in dimensionless time.
    double u = Randomizer.nextExponential(1.0);
    
    // Determine arrival point in real time
    for (int i=startIdx; i<events.size(); i++) {
        
        CFEventList.Event event = events.get(i);
        
        double t = Math.max(conv.getHeight1(), event.getHeight());
    
        // Determine length of interval in dimensionless time
        double intervalArea;
        if (i<events.size()-1)
            intervalArea = popFunc.getIntegral(t, events.get(i+1).getHeight());
        else
            intervalArea = Double.POSITIVE_INFINITY;
        
        // Check whether arrival falls within this interval
        if (u<intervalArea*event.getLineageCount()) {
            
            // Set arrival point in real time
            conv.setHeight2(popFunc.getInverseIntensity(
                    popFunc.getIntensity(t) + u/event.getLineageCount()));
            
            // Attach to random clonal frame lineage extant at this time
            int z = Randomizer.nextInt(event.getLineageCount());
            for (Node node : acg.getNodesAsArray()) {
                if (conv.getHeight2()>node.getHeight() &&
                        (node.isRoot() || conv.getHeight2()<node.getParent().getHeight())) {
                    if (z==0) {
                        conv.setNode2(node);
                        break;
                    } else
                        z -= 1;
                }
            }

            logP += -u + Math.log(1.0/popFunc.getPopSize(conv.getHeight2()));
            break;
            
        } else {
            u -= intervalArea*event.getLineageCount();
            logP += -intervalArea*event.getLineageCount();
        }
    }
    
    return logP;
}