Examples with CoalescentSimulator dr.evolution.coalescent.CoalescentSimulator used on opensource projects

Example 1 with CoalescentSimulator

use of dr.evolution.coalescent.CoalescentSimulator in project beast-mcmc by beast-dev.

the class TransmissionTreeToVirusTree method simulateCoalescent.

private FlexibleNode simulateCoalescent(ArrayList<SimpleNode> nodes, DemographicFunction demogFunct, double maxHeight) {
    double earliestNodeHeight = Double.NEGATIVE_INFINITY;
    for (SimpleNode node : nodes) {
        if (node.getHeight() > earliestNodeHeight) {
            earliestNodeHeight = node.getHeight();
        }
    }
    double maxLastInterval = earliestNodeHeight;
    double probNoCoalesenceInTime = Math.exp(demogFunct.getIntensity(maxLastInterval));
    coalescentProbability *= (1 - probNoCoalesenceInTime);
    CoalescentSimulator simulator = new CoalescentSimulator();
    SimpleNode root;
    SimpleNode[] simResults;
    int failCount = 0;
    do {
        simResults = simulator.simulateCoalescent(nodes.toArray(new SimpleNode[nodes.size()]), demogFunct, -maxHeight, 0, true);
        if (simResults.length > 1) {
            failCount++;
            System.out.println("Failed to coalesce lineages: " + failCount);
        }
    } while (simResults.length != 1);
    root = simResults[0];
    SimpleTree simpleTreelet = new SimpleTree(root);
    for (int i = 0; i < simpleTreelet.getNodeCount(); i++) {
        SimpleNode node = (SimpleNode) simpleTreelet.getNode(i);
        node.setHeight(node.getHeight() + maxHeight);
    }
    return new FlexibleNode(simpleTreelet, root, true);
}

Example 2 with CoalescentSimulator

use of dr.evolution.coalescent.CoalescentSimulator in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateTopology.

// END: annotateTree
static void simulateTopology() {
    try {
        System.out.println("Test case 1: simulateTopology");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        // set demographic function
        ExponentialGrowth exponentialGrowth = new ExponentialGrowth(Units.Type.YEARS);
        exponentialGrowth.setN0(10);
        exponentialGrowth.setGrowthRate(0.5);
        Taxa taxa = new Taxa();
        for (Taxon taxon : tree.asList()) {
            double absoluteHeight = Utils.getAbsoluteTaxonHeight(taxon, tree);
            taxon.setAttribute(Utils.ABSOLUTE_HEIGHT, absoluteHeight);
            // taxon.setAttribute("date", new Date(absoluteHeight,
            // Units.Type.YEARS, true));
            taxa.addTaxon(taxon);
        }
        // END: taxon loop
        CoalescentSimulator topologySimulator = new CoalescentSimulator();
        TreeModel treeModel = new TreeModel(topologySimulator.simulateTree(taxa, exponentialGrowth));
        System.out.println(treeModel.toString());
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create substitution model
        Parameter kappa = new Parameter.Default(1, 10);
        HKY hky = new HKY(kappa, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(hky);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        System.out.println(simulator.simulate(simulateInPar, false).toString());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 3 with CoalescentSimulator

use of dr.evolution.coalescent.CoalescentSimulator in project beast-mcmc by beast-dev.

the class PartitionData method createTreeModel.

public TreeModel createTreeModel() {
    TreeModel treeModel = null;
    if (this.demographicModelIndex == 0 && this.record.isTreeSet()) {
        treeModel = new TreeModel(this.record.getTree());
    } else if ((this.demographicModelIndex > 0 && this.demographicModelIndex <= lastImplementedIndex) && this.record.isTreeSet()) {
        Taxa taxa = new Taxa(this.record.getTree().asList());
        CoalescentSimulator topologySimulator = new CoalescentSimulator();
        treeModel = new TreeModel(topologySimulator.simulateTree(taxa, createDemographicFunction()));
    } else if ((this.demographicModelIndex > 0 && this.demographicModelIndex <= lastImplementedIndex) && this.record.isTaxaSet()) {
        Taxa taxa = this.record.getTaxa();
        CoalescentSimulator topologySimulator = new CoalescentSimulator();
        treeModel = new TreeModel(topologySimulator.simulateTree(taxa, createDemographicFunction()));
    //			} else if (this.demographicModelIndex == 0 && this.record.taxaSet) { 
    //			throw new RuntimeException("Data and demographic model incompatible for partition ");	
    } else {
        throw new RuntimeException("Data and demographic model incompatible.");
    }
    return treeModel;
}

Example 4 with CoalescentSimulator

use of dr.evolution.coalescent.CoalescentSimulator 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 5 with CoalescentSimulator

use of dr.evolution.coalescent.CoalescentSimulator in project beast-mcmc by beast-dev.

the class TraceCorrelationAssert method createTreeModel.

private void createTreeModel(ConstantPopulation constant) {
    CoalescentSimulator simulator = new CoalescentSimulator();
    Tree tree = simulator.simulateTree(alignment, constant);
    //treeModel
    treeModel = new TreeModel(tree);
}