Search in sources :

Example 6 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ACGScaler method proposal.

@Override
public double proposal() {
    // Keep track of number of positively scaled elements minus
    // negatively scaled elements.
    int count = 0;
    // Choose scaling factor:
    double f = scaleParam + Randomizer.nextDouble() * (1.0 / scaleParam - scaleParam);
    // Scale clonal frame:
    if (rootOnly) {
        acg.getRoot().setHeight(acg.getRoot().getHeight() * f);
        count += 1;
    } else {
        for (Node node : acg.getInternalNodes()) {
            node.setHeight(node.getHeight() * f);
            count += 1;
        }
    }
    // Scale recombinant edges:
    for (Locus locus : acg.getLoci()) {
        for (Conversion recomb : acg.getConversions(locus)) {
            if (!rootOnly || recomb.getNode1().getParent().isRoot()) {
                recomb.setHeight1(recomb.getHeight1() * f);
                count += 1;
            }
            if (!rootOnly || recomb.getNode2().isRoot() || recomb.getNode2().getParent().isRoot()) {
                recomb.setHeight2(recomb.getHeight2() * f);
                count += 1;
            }
            if (recomb.getHeight1() < recomb.getNode1().getHeight() || recomb.getHeight2() < recomb.getNode2().getHeight() || recomb.getHeight1() > recomb.getHeight2())
                return Double.NEGATIVE_INFINITY;
        }
    }
    // Check for illegal node heights:
    if (rootOnly) {
        for (Node node : acg.getRoot().getChildren()) {
            if (node.getHeight() > node.getParent().getHeight())
                return Double.NEGATIVE_INFINITY;
        }
    } else {
        for (Node node : acg.getExternalNodes()) {
            if (node.getHeight() > node.getParent().getHeight())
                return Double.NEGATIVE_INFINITY;
        }
    }
    for (RealParameter param : parametersInput.get()) {
        try {
            param.startEditing(null);
            param.scale(f);
        } catch (Exception e) {
            // bounds.  Needs to change!
            return Double.NEGATIVE_INFINITY;
        }
        count += param.getDimension();
    }
    for (RealParameter paramInv : parametersInverseInput.get()) {
        try {
            paramInv.startEditing(null);
            paramInv.scale(1.0 / f);
        } catch (Exception e) {
            // bounds.  Needs to change!
            return Double.NEGATIVE_INFINITY;
        }
        count -= paramInv.getDimension();
    }
    // Return log of Hastings ratio:
    return (count - 2) * Math.log(f);
}
Also used : Node(beast.evolution.tree.Node) RealParameter(beast.core.parameter.RealParameter) Locus(bacter.Locus) Conversion(bacter.Conversion)

Example 7 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class AddRemoveConversion method drawNewConversion.

/**
 * Add new conversion to ACG, returning the probability density of the
 * new edge and converted region location.
 *
 * @return log of proposal density
 */
public double drawNewConversion() {
    Conversion newConversion = new Conversion();
    double logP = attachEdge(newConversion) + drawAffectedRegion(newConversion);
    acg.addConversion(newConversion);
    return logP;
}
Also used : Conversion(bacter.Conversion)

Example 8 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class AddRemoveDetour method removeDetour.

/**
 * Detour deletion variant of move.
 *
 * @return log HGF
 */
private double removeDetour() {
    double logHGF = 0.0;
    // Choose non-root detour edge at random
    Node detour = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
    logHGF -= Math.log(1.0 / (acg.getNodeCount() - 1));
    List<Conversion> convApotentials = new ArrayList<>();
    List<Conversion> convBpotentials = new ArrayList<>();
    for (Locus locus : acg.getLoci()) {
        for (Conversion conv : acg.getConversions(locus)) {
            if (conv.getNode2() == detour && conv.getNode1() != detour)
                convApotentials.add(conv);
            if (conv.getNode1() == detour && conv.getNode2() != detour)
                convBpotentials.add(conv);
        }
    }
    if (convApotentials.isEmpty() || convBpotentials.isEmpty())
        return Double.NEGATIVE_INFINITY;
    Conversion convA = convApotentials.get(Randomizer.nextInt(convApotentials.size()));
    Conversion convB = convBpotentials.get(Randomizer.nextInt(convBpotentials.size()));
    if (convA.getHeight2() > convB.getHeight1())
        return Double.NEGATIVE_INFINITY;
    logHGF -= Math.log(1.0 / (convApotentials.size() * convBpotentials.size()));
    double tLowerBound = convA.getHeight1();
    double tUpperBound = convB.getHeight2();
    Conversion conv = new Conversion();
    conv.setNode1(convA.getNode1());
    conv.setHeight1(convA.getHeight1());
    conv.setNode2(convB.getNode2());
    conv.setHeight2(convB.getHeight2());
    conv.setStartSite(convA.getStartSite());
    conv.setEndSite(convA.getEndSite());
    conv.setLocus(convA.getLocus());
    acg.deleteConversion(convA);
    acg.deleteConversion(convB);
    acg.addConversion(conv);
    logHGF += Math.log(1.0 / acg.getTotalConvCount()) + Math.log(1.0 / (acg.getNodeCount() - 1)) + Math.log(2.0) - 2.0 * Math.log(tUpperBound - tLowerBound) + getAffectedRegionProb(convB);
    return logHGF;
}
Also used : Node(beast.evolution.tree.Node) ArrayList(java.util.ArrayList) Locus(bacter.Locus) Conversion(bacter.Conversion)

Example 9 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class AddRemoveRedundantConversion method getRedundantConversions.

/**
 * Obtain list of redundant conversions.
 *
 * @param cfNode node identifying edge on CF
 * @return conversion list
 */
private List<Conversion> getRedundantConversions(Node cfNode) {
    Node cfParent = cfNode.getParent();
    List<Conversion> redundantConversions = new ArrayList<>();
    double maxL = acg.getRoot().getHeight() * apertureInput.get();
    for (Locus locus : acg.getLoci()) {
        for (Conversion conv : acg.getConversions(locus)) {
            if (((conv.getNode1() == cfNode || conv.getNode1().getParent() == cfNode) && Math.abs(conv.getHeight1() - cfNode.getHeight()) < maxL) && (conv.getNode2() == cfParent || conv.getNode2().getParent() == cfParent) && Math.abs(conv.getHeight2() - cfParent.getHeight()) < maxL) {
                redundantConversions.add(conv);
            }
        }
    }
    return redundantConversions;
}
Also used : Node(beast.evolution.tree.Node) ArrayList(java.util.ArrayList) Locus(bacter.Locus) Conversion(bacter.Conversion)

Example 10 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class AddRemoveRedundantConversion method proposal.

@Override
public double proposal() {
    double logHGF = 0;
    // Select non-root CF node
    Node cfNode = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
    Node cfParent = cfNode.getParent();
    double maxL = apertureInput.get() * acg.getRoot().getHeight();
    if (Randomizer.nextBoolean()) {
        // Add redundant conversion
        Conversion newConv = new Conversion();
        double L = Math.min(getEdgeLength(cfNode), maxL);
        for (Node child : cfNode.getChildren()) L += Math.min(getEdgeLength(child), maxL);
        logHGF -= Math.log(1.0 / L);
        double fromPoint = L * Randomizer.nextDouble();
        if (fromPoint < Math.min(getEdgeLength(cfNode), maxL)) {
            newConv.setNode1(cfNode);
            newConv.setHeight1(cfNode.getHeight() + fromPoint);
        } else {
            fromPoint -= Math.min(getEdgeLength(cfNode), maxL);
            for (Node child : cfNode.getChildren()) {
                if (fromPoint < Math.min(getEdgeLength(child), maxL)) {
                    newConv.setNode1(child);
                    newConv.setHeight1(cfNode.getHeight() - fromPoint);
                    break;
                }
                fromPoint -= Math.min(getEdgeLength(child), maxL);
            }
        }
        L = Math.min(getEdgeLength(cfParent), maxL);
        for (Node child : cfParent.getChildren()) L += Math.min(getEdgeLength(child), maxL);
        logHGF -= Math.log(1.0 / L);
        double toPoint = L * Randomizer.nextDouble();
        if (toPoint < Math.min(getEdgeLength(cfParent), maxL)) {
            newConv.setNode2(cfParent);
            newConv.setHeight2(cfParent.getHeight() + toPoint);
        } else {
            toPoint -= Math.min(getEdgeLength(cfParent), maxL);
            for (Node child : cfParent.getChildren()) {
                if (toPoint < Math.min(getEdgeLength(child), maxL)) {
                    newConv.setNode2(child);
                    newConv.setHeight2(cfParent.getHeight() - toPoint);
                    break;
                }
                toPoint -= Math.min(getEdgeLength(child), maxL);
            }
        }
        if (newConv.getHeight1() > newConv.getHeight2())
            return Double.NEGATIVE_INFINITY;
        logHGF -= drawAffectedRegion(newConv);
        // Add conversion
        acg.addConversion(newConv);
        // Add probability of reverse move deleting this conversion
        // to HGF:
        logHGF += Math.log(1.0 / getRedundantConversions(cfNode).size());
    } else {
        // Remove
        // Identify redundant conversions
        List<Conversion> redundantConversions = getRedundantConversions(cfNode);
        if (redundantConversions.size() == 0)
            return Double.NEGATIVE_INFINITY;
        // Choose conversion to remove
        Conversion conv = redundantConversions.get(Randomizer.nextInt(redundantConversions.size()));
        logHGF -= Math.log(1.0 / redundantConversions.size());
        // Add probability of reverse move generating this conversion
        // to HGF:
        double L = Math.min(getEdgeLength(cfNode), maxL);
        for (Node child : cfNode.getChildren()) L += Math.min(getEdgeLength(child), maxL);
        logHGF += Math.log(1.0 / L);
        L = Math.min(getEdgeLength(cfParent), maxL);
        for (Node child : cfParent.getChildren()) L += Math.min(getEdgeLength(child), maxL);
        logHGF += Math.log(1.0 / L);
        logHGF += getAffectedRegionProb(conv);
        // Remove conversion
        acg.deleteConversion(conv);
    }
    return logHGF;
}
Also used : Node(beast.evolution.tree.Node) Conversion(bacter.Conversion)

Aggregations

Conversion (bacter.Conversion)49 Locus (bacter.Locus)27 Node (beast.evolution.tree.Node)24 ConversionGraph (bacter.ConversionGraph)7 ArrayList (java.util.ArrayList)7 Test (org.junit.Test)6 RealParameter (beast.core.parameter.RealParameter)4 SiteModel (beast.evolution.sitemodel.SiteModel)4 JukesCantor (beast.evolution.substitutionmodel.JukesCantor)4 TaxonSet (beast.evolution.alignment.TaxonSet)3 ConstantPopulation (beast.evolution.tree.coalescent.ConstantPopulation)3 ClusterTree (beast.util.ClusterTree)3 PrintStream (java.io.PrintStream)3 SimulatedACG (bacter.model.SimulatedACG)2 CFEventList (bacter.CFEventList)1 ACGCoalescent (bacter.model.ACGCoalescent)1 ACGLogReader (bacter.util.ACGLogReader)1 BacterACGLogReader (bacter.util.BacterACGLogReader)1 COACGLogFileReader (bacter.util.COACGLogFileReader)1 State (beast.core.State)1