Examples with Locus bacter.Locus used on opensource projects

Example 6 with Locus

use of bacter.Locus 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);
}

Example 7 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class AddRemoveConversion method main.

public static void main(String[] args) throws Exception {
    ConversionGraph acg = new ConversionGraph();
    ConstantPopulation popFunc = new ConstantPopulation();
    AddRemoveConversion operator = new AddRemoveConversion();
    operator.initByName("weight", 1.0, "acg", acg, "populationModel", popFunc, "rho", new RealParameter(Double.toString(1.0 / 10000.0)), "delta", new RealParameter("50.0"));
    popFunc.initByName("popSize", new RealParameter("1.0"));
    TaxonSet taxonSet = new TaxonSet();
    taxonSet.taxonsetInput.setValue(new Taxon("t1"), taxonSet);
    taxonSet.taxonsetInput.setValue(new Taxon("t2"), taxonSet);
    Locus locus = new Locus("locus", 10000);
    try (PrintStream ps = new PrintStream("out.txt")) {
        for (int i = 0; i < 100000; i++) {
            acg.initByName("locus", locus, "taxonset", taxonSet, "fromString", "(0:1.0,1:1.0)2:0.0;");
            operator.drawNewConversion();
            ps.println(acg.getConversions(locus).get(0).getStartSite() + " " + acg.getConversions(locus).get(0).getEndSite());
        }
    }
}

Example 8 with Locus

use of bacter.Locus 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;
}

Example 9 with Locus

use of bacter.Locus 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;
}

Example 10 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class COACGLogFileReader method iterator.

@Override
public Iterator<ConversionGraph> iterator() {
    try {
        reset();
        skipBurnin();
    } catch (XMLStreamException | IOException e) {
        throw new IllegalStateException(e.getMessage());
    }
    return new Iterator<ConversionGraph>() {

        int current = 0;

        @Override
        public boolean hasNext() {
            return current < nACGs;
        }

        @Override
        public ConversionGraph next() {
            current += 1;
            ConversionGraph acg = new ConversionGraph();
            for (Locus locus : loci) acg.lociInput.setValue(locus, acg);
            try {
                acg.initAndValidate();
            } catch (Exception e) {
                throw new IllegalStateException(e.getMessage());
            }
            String newick = null;
            List<Integer> rStarts = new ArrayList<>();
            List<Integer> rEnds = new ArrayList<>();
            List<Integer> rEFroms = new ArrayList<>();
            List<Integer> rETos = new ArrayList<>();
            List<Double> rAFroms = new ArrayList<>();
            List<Double> rATos = new ArrayList<>();
            try {
                skipUntil("iteration");
                while (xmlStreamReader.hasNext()) {
                    int type = xmlStreamReader.next();
                    if (type == XMLStreamConstants.END_ELEMENT && xmlStreamReader.getLocalName().toLowerCase().equals("iteration"))
                        break;
                    if (type != XMLStreamConstants.START_ELEMENT)
                        continue;
                    switch(xmlStreamReader.getLocalName().toLowerCase()) {
                        case "tree":
                            xmlStreamReader.next();
                            newick = xmlStreamReader.getText().trim();
                            break;
                        case "recedge":
                            while ((type = xmlStreamReader.next()) != XMLStreamConstants.END_ELEMENT || !xmlStreamReader.getLocalName().toLowerCase().equals("recedge")) {
                                if (type != XMLStreamConstants.START_ELEMENT)
                                    continue;
                                switch(xmlStreamReader.getLocalName().toLowerCase()) {
                                    case "start":
                                        xmlStreamReader.next();
                                        rStarts.add(Integer.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    case "end":
                                        xmlStreamReader.next();
                                        rEnds.add(Integer.valueOf(xmlStreamReader.getText().trim()) - 1);
                                        break;
                                    case "efrom":
                                        xmlStreamReader.next();
                                        rEFroms.add(Integer.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    case "eto":
                                        xmlStreamReader.next();
                                        rETos.add(Integer.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    case "afrom":
                                        xmlStreamReader.next();
                                        rAFroms.add(Double.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    case "ato":
                                        xmlStreamReader.next();
                                        rATos.add(Double.valueOf(xmlStreamReader.getText().trim()));
                                        break;
                                    default:
                                }
                            }
                            break;
                        default:
                    }
                }
            } catch (XMLStreamException e) {
                throw new IllegalStateException(e.getMessage());
            }
            acg.fromExtendedNewick(newick, true, 0);
            for (int i = 0; i < rStarts.size(); i++) {
                Node fromNode = acg.getNode(rEFroms.get(i));
                Node toNode = acg.getNode(rETos.get(i));
                Conversion conv = new Conversion(toNode, toNode.getHeight() + rATos.get(i), fromNode, fromNode.getHeight() + rAFroms.get(i), rStarts.get(i), rEnds.get(i), acg, loci.get(0));
                acg.addConversion(conv);
            }
            current += 1;
            return acg;
        }
    };
}