Examples with Taxa dr.evolution.util.Taxa used on opensource projects

Example 66 with Taxa

use of dr.evolution.util.Taxa in project beast-mcmc by beast-dev.

the class VariableCoalescentSimulator method readSampleFile.

private static Taxa readSampleFile(String fileName, double generationTime) throws IOException {
    BufferedReader reader = new BufferedReader(new FileReader(fileName));
    String line = reader.readLine();
    Taxa taxa = new Taxa();
    int id = 0;
    while (line != null) {
        if (!line.startsWith("#")) {
            String[] tokens = line.split("[\t ]+");
            // sample times are in the same units as simulation
            double sampleTime = Double.parseDouble(tokens[0]) / generationTime;
            int count = Integer.parseInt(tokens[1]);
            for (int i = 0; i < count; i++) {
                Taxon taxon = new Taxon(id + "");
                taxon.setAttribute(dr.evolution.util.Date.DATE, new Date(sampleTime, Units.Type.GENERATIONS, true));
                taxa.addTaxon(taxon);
                id += 1;
            }
        }
        line = reader.readLine();
    }
    return taxa;
}

Example 67 with Taxa

use of dr.evolution.util.Taxa in project beast-mcmc by beast-dev.

the class VariableCoalescentSimulator method main.

public static void main(String[] arg) throws IOException {
    long startTime = System.currentTimeMillis();
    Options options = new Options(arg, 0, 7);
    options.getSet().addOption("g", Options.Separator.EQUALS, Options.Multiplicity.ZERO_OR_ONE);
    options.getSet().addOption("n", Options.Separator.EQUALS, Options.Multiplicity.ZERO_OR_ONE);
    options.getSet().addOption("p", Options.Separator.EQUALS, Options.Multiplicity.ZERO_OR_ONE);
    options.getSet().addOption("se", Options.Separator.EQUALS, Options.Multiplicity.ZERO_OR_ONE);
    options.getSet().addOption("ss", Options.Separator.EQUALS, Options.Multiplicity.ZERO_OR_ONE);
    options.getSet().addOption("reps", Options.Separator.EQUALS, Options.Multiplicity.ZERO_OR_ONE);
    options.getSet().addOption("f", Options.Multiplicity.ZERO_OR_ONE);
    if (!options.check()) {
        System.out.println(options.getCheckErrors());
        System.out.println();
        printUsage();
        System.exit(1);
    }
    double generationTime = 1.0;
    double popSizeScale = 1.0;
    int n = 50;
    double ss = -1;
    double se = -1;
    int reps = 1;
    boolean timeForward = options.getSet().isSet("f");
    if (options.getSet().isSet("g")) {
        String g = options.getSet().getOption("g").getResultValue(0);
        generationTime = Double.parseDouble(g);
        System.out.println("generation time = " + g);
    }
    if (options.getSet().isSet("n")) {
        String sampleSize = options.getSet().getOption("n").getResultValue(0);
        n = Integer.parseInt(sampleSize);
        System.out.println("sample size = " + n);
    }
    if (options.getSet().isSet("p")) {
        String p = options.getSet().getOption("p").getResultValue(0);
        popSizeScale = Double.parseDouble(p);
        System.out.println("population size scale factor = " + p);
    }
    if (options.getSet().isSet("ss")) {
        String sampleStart = options.getSet().getOption("ss").getResultValue(0);
        ss = Double.parseDouble(sampleStart);
        System.out.println("sample start time = " + ss);
    }
    if (options.getSet().isSet("se")) {
        String sampleEnd = options.getSet().getOption("se").getResultValue(0);
        se = Double.parseDouble(sampleEnd);
        System.out.println("sample end time = " + se);
    }
    if (options.getSet().isSet("reps")) {
        String replicates = options.getSet().getOption("reps").getResultValue(0);
        reps = Integer.parseInt(replicates);
        System.out.println("replicates = " + reps);
    }
    String filename = options.getSet().getData().get(0);
    String outfile = options.getSet().getData().get(1);
    // READ DEMOGRAPHIC FUNCTION
    BufferedReader reader = new BufferedReader(new FileReader(filename));
    List<Double> times = new ArrayList<Double>();
    String line = reader.readLine();
    String[] tokens = line.trim().split("[\t ]+");
    if (tokens.length < 2)
        throw new RuntimeException();
    ArrayList<ArrayList> popSizes = new ArrayList<ArrayList>();
    while (line != null) {
        double time = Double.parseDouble(tokens[0]) / generationTime;
        times.add(time);
        for (int i = 1; i < tokens.length; i++) {
            popSizes.add(new ArrayList<Double>());
            popSizes.get(i - 1).add(Double.parseDouble(tokens[i]));
        }
        line = reader.readLine();
        if (line != null) {
            tokens = line.trim().split("[\t ]+");
            if (tokens.length != popSizes.size() + 1)
                throw new RuntimeException();
        }
    }
    reader.close();
    // GENERATE SAMPLE
    double lastTime = times.get(times.size() - 1);
    if (ss == -1) {
        ss = timeForward ? lastTime : times.get(0);
    }
    if (se == -1) {
        se = timeForward ? lastTime : times.get(0);
    }
    double dt = (se - ss) / ((double) n - 1.0);
    double time = ss;
    Taxa taxa = new Taxa();
    for (int i = 0; i < n; i++) {
        double sampleTime;
        if (timeForward) {
            sampleTime = (lastTime - time) / generationTime;
        } else
            sampleTime = time / generationTime;
        Taxon taxon = new Taxon(i + "");
        taxon.setAttribute(dr.evolution.util.Date.DATE, new Date(sampleTime, Units.Type.GENERATIONS, true));
        taxa.addTaxon(taxon);
        time += dt;
    }
    double minTheta = Double.MAX_VALUE;
    double maxTheta = 0.0;
    PrintWriter out = new PrintWriter(new FileWriter(outfile));
    int popHistory = 0;
    PiecewiseLinearPopulation[] demography = new PiecewiseLinearPopulation[popSizes.size()];
    for (List<Double> popSize : popSizes) {
        double[] thetas = new double[popSize.size()];
        double[] intervals = new double[times.size() - 1];
        for (int i = intervals.length; i >= 0; i--) {
            int j = timeForward ? intervals.length - i : i - 1;
            int k = timeForward ? i : intervals.length - i + 1;
            if (i != 0)
                intervals[j] = times.get(k) - times.get(k - 1);
            double theta = popSize.get(k) * popSizeScale;
            thetas[j] = theta;
            if (theta < minTheta) {
                minTheta = theta;
            }
            if (theta > maxTheta) {
                maxTheta = theta;
            }
        // System.out.println(t + "\t" + theta);
        }
        System.out.println("N" + popHistory + "(t) range = [" + minTheta + ", " + maxTheta + "]");
        demography[popHistory] = new PiecewiseLinearPopulation(intervals, thetas, Units.Type.GENERATIONS);
        popHistory += 1;
    }
    CoalescentSimulator simulator = new CoalescentSimulator();
    for (int i = 0; i < reps; i++) {
        out.println("#rep " + i);
        for (int j = 0; j < demography.length; j++) {
            Tree tree = simulator.simulateTree(taxa, demography[j]);
            out.println(TreeUtils.newick(tree));
        // System.err.println(Tree.Utils.newick(tree));
        }
    }
    out.flush();
    out.close();
    long stopTime = System.currentTimeMillis();
    System.out.println("Took " + (stopTime - startTime) / 1000.0 + " seconds");
}

Example 68 with Taxa

use of dr.evolution.util.Taxa in project beast-mcmc by beast-dev.

the class GibbsIndependentCoalescentOperator method doOperation.

/**
 * change the parameter and return the hastings ratio.
 */
public double doOperation() {
    CoalescentSimulator simulator = new CoalescentSimulator();
    // store the generated tree
    Tree tree;
    List<TaxonList> taxonLists = new ArrayList<TaxonList>();
    List<Tree> subtrees = new ArrayList<Tree>();
    // should have one child that is node
    for (int i = 0; i < xo.getChildCount(); i++) {
        final Object child = xo.getChild(i);
        // AER - swapped the order of these round because Trees are TaxonLists...
        if (child instanceof TreeModel) {
        // do nothing
        } else if (child instanceof Tree) {
            // don't do this as all subtrees will remain identical
            // subtrees.add((Tree) child);
            Tree currentTree = (Tree) child;
            Tree[] newTrees = new Tree[currentTree.getTaxonCount()];
            for (int j = 0; j < currentTree.getTaxonCount(); j++) {
                Taxa tip = new Taxa();
                tip.addTaxon(currentTree.getTaxon(j));
                newTrees[j] = simulator.simulateTree(tip, demoModel);
            }
            Tree newTree = simulator.simulateTree(newTrees, demoModel, rootHeight, newTrees.length != 1);
            subtrees.add(newTree);
        } else if (child instanceof TaxonList) {
            taxonLists.add((TaxonList) child);
        }
    }
    if (taxonLists.size() == 0) {
        if (subtrees.size() == 1) {
            // System.out.println(" *** 1 *** ");
            tree = subtrees.get(0);
            // this would be the normal way to do it
            treeModel.beginTreeEdit();
            // now it's allowed to adjust the tree structure
            treeModel.adoptTreeStructure(tree);
            // endTreeEdit() would then fire the events
            treeModel.endTreeEdit();
            return 0;
        }
        throw new RuntimeException("Expected at least one taxonList or two subtrees in " + getOperatorName() + " XML specification.");
    }
    Taxa remainingTaxa = new Taxa();
    for (int i = 0; i < taxonLists.size(); i++) {
        remainingTaxa.addTaxa(taxonLists.get(i));
    }
    for (int i = 0; i < subtrees.size(); i++) {
        remainingTaxa.removeTaxa(subtrees.get(i));
    }
    Tree[] trees = new Tree[subtrees.size() + remainingTaxa.getTaxonCount()];
    // add the preset trees
    for (int i = 0; i < subtrees.size(); i++) {
        trees[i] = subtrees.get(i);
    }
    // add all the remaining taxa in as single tip trees...
    for (int i = 0; i < remainingTaxa.getTaxonCount(); i++) {
        Taxa tip = new Taxa();
        tip.addTaxon(remainingTaxa.getTaxon(i));
        trees[i + subtrees.size()] = simulator.simulateTree(tip, demoModel);
    }
    tree = simulator.simulateTree(trees, demoModel, rootHeight, trees.length != 1);
    // this would be the normal way to do it
    treeModel.beginTreeEdit();
    // now it's allowed to adjust the tree structure
    treeModel.adoptTreeStructure(tree);
    // endTreeEdit() would then fire the events
    treeModel.endTreeEdit();
    return 0;
}

Example 69 with Taxa

use of dr.evolution.util.Taxa in project beast-mcmc by beast-dev.

the class IndependentCoalescentSampler method doOperation.

/**
 * change the parameter and return the hastings ratio.
 */
public double doOperation() {
    CoalescentSimulator simulator = new CoalescentSimulator();
    List<TaxonList> taxonLists = new ArrayList<TaxonList>();
    double rootHeight = -1.0;
    double oldLikelihood = 0.0;
    double newLikelihood = 0.0;
    // should have one child that is node
    for (int i = 0; i < xo.getChildCount(); i++) {
        final Object child = xo.getChild(i);
        // careful: Trees are TaxonLists ... (AER); see OldCoalescentSimulatorParser
        if (child instanceof Tree) {
        // do nothing
        } else if (child instanceof TaxonList) {
            // taxonLists.add((TaxonList) child);
            taxonLists.add((Taxa) child);
            // taxa added
            break;
        }
    }
    try {
        Tree[] trees = new Tree[taxonLists.size()];
        // simulate each taxonList separately
        for (int i = 0; i < taxonLists.size(); i++) {
            trees[i] = simulator.simulateTree(taxonLists.get(i), demoModel);
        }
        oldLikelihood = coalescent.getLogLikelihood();
        SimpleTree simTree = simulator.simulateTree(trees, demoModel, rootHeight, trees.length != 1);
        // this would be the normal way to do it
        treeModel.beginTreeEdit();
        // now it's allowed to adjust the tree structure
        treeModel.adoptTreeStructure(simTree);
        // endTreeEdit() would then fire the events
        treeModel.endTreeEdit();
        newLikelihood = coalescent.getLogLikelihood();
    } catch (IllegalArgumentException iae) {
        try {
            throw new XMLParseException(iae.getMessage());
        } catch (XMLParseException e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
    }
    return oldLikelihood - newLikelihood;
}

Example 70 with Taxa

use of dr.evolution.util.Taxa in project beast-mcmc by beast-dev.

the class CoalescentSimulatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    CoalescentSimulator simulator = new CoalescentSimulator();
    DemographicModel demoModel = (DemographicModel) xo.getChild(DemographicModel.class);
    List<TaxonList> taxonLists = new ArrayList<TaxonList>();
    List<Tree> subtrees = new ArrayList<Tree>();
    double height = xo.getAttribute(HEIGHT, Double.NaN);
    Tree constraintsTree = null;
    if (xo.hasChildNamed((CONSTRAINTS_TREE))) {
        XMLObject cxo = xo.getChild(CONSTRAINTS_TREE);
        constraintsTree = (Tree) cxo.getChild(Tree.class);
    }
    // should have one child that is node
    for (int i = 0; i < xo.getChildCount(); i++) {
        final Object child = xo.getChild(i);
        if (child != constraintsTree) {
            // AER - swapped the order of these round because Trees are TaxonLists...
            if (child instanceof Tree) {
                subtrees.add((Tree) child);
            } else if (child instanceof TaxonList) {
                taxonLists.add((TaxonList) child);
            }
        }
    }
    if (taxonLists.size() == 0) {
        if (subtrees.size() == 1) {
            return subtrees.get(0);
        }
        if (constraintsTree == null) {
            throw new XMLParseException("Expected at least one taxonList or two subtrees or a constraints tree in " + getParserName() + " element.");
        }
    }
    Taxa remainingTaxa = new Taxa();
    for (int i = 0; i < taxonLists.size(); i++) {
        remainingTaxa.addTaxa(taxonLists.get(i));
    }
    for (int i = 0; i < subtrees.size(); i++) {
        remainingTaxa.removeTaxa(subtrees.get(i));
    }
    if (constraintsTree != null) {
        remainingTaxa.removeTaxa(constraintsTree);
    }
    try {
        if (constraintsTree != null) {
            subtrees.add(simulator.simulateTree(constraintsTree, constraintsTree.getRoot(), demoModel));
        }
        Tree[] trees = new Tree[subtrees.size() + remainingTaxa.getTaxonCount()];
        // add the preset trees
        for (int i = 0; i < subtrees.size(); i++) {
            trees[i] = subtrees.get(i);
        }
        // add all the remaining taxa in as single tip trees...
        for (int i = 0; i < remainingTaxa.getTaxonCount(); i++) {
            Taxa tip = new Taxa();
            tip.addTaxon(remainingTaxa.getTaxon(i));
            trees[i + subtrees.size()] = simulator.simulateTree(tip, demoModel);
        }
        return simulator.simulateTree(trees, demoModel, height, trees.length != 1);
    } catch (IllegalArgumentException iae) {
        throw new XMLParseException(iae.getMessage());
    }
}