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

Example 51 with Taxa

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

the class BeagleSequenceSimulatorConsoleApp method simulate.

// END: constructor
public void simulate(String[] args) {
    try {
        // /////////////////////////////////
        // ---SPLIT PARTITION ARGUMENTS---//
        // /////////////////////////////////
        int from = 0;
        int to = 0;
        ArrayList<String[]> argsList = new ArrayList<String[]>();
        for (String arg : args) {
            if (arg.equalsIgnoreCase(SPLIT_PARTITION)) {
                argsList.add(Arrays.copyOfRange(args, from, to));
                from = to + 1;
            }
            // END: split check
            to++;
        }
        if (args[0].contains(HELP)) {
            gracefullyExit(null);
        } else if (argsList.size() == 0) {
            gracefullyExit("Empty or incorrect arguments list.");
        }
        // END: failed split check
        String[] leftoverArguments = Arrays.copyOfRange(args, from, args.length);
        if (leftoverArguments.length > 2) {
            gracefullyExit("Unrecognized option " + leftoverArguments[2]);
        }
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        for (String[] partitionArgs : argsList) {
            // /////////////
            // ---PARSE---//
            // /////////////
            arguments.parseArguments(partitionArgs);
            // ///////////////////
            // ---INTERROGATE---//
            // ///////////////////
            PartitionData data = new PartitionData();
            String option = null;
            double[] values = null;
            if (partitionArgs.length == 0 || arguments.hasOption(HELP)) {
                gracefullyExit(null);
            }
            // Tree / Taxa
            if (arguments.hasOption(TREE_FILE)) {
                File treeFile = new File(arguments.getStringOption(TREE_FILE));
                Tree tree = Utils.importTreeFromFile(treeFile);
                data.record = new TreesTableRecord(treeFile.getName(), tree);
            } else if (arguments.hasOption(TAXA_SET)) {
                File taxaFile = new File(arguments.getStringOption(TAXA_SET));
                Taxa taxa = Utils.importTaxaFromFile(taxaFile);
                data.record = new TreesTableRecord(taxaFile.getName(), taxa);
            } else {
                throw new RuntimeException("Tree file / Taxa set not specified.");
            }
            // Demographic Model
            if (arguments.hasOption(DEMOGRAPHIC_MODEL)) {
                option = arguments.getStringOption(DEMOGRAPHIC_MODEL);
                if (option.equalsIgnoreCase(NO_DEMOGRAPHIC_MODEL)) {
                    int index = 0;
                    data.demographicModelIndex = index;
                } else if (option.equalsIgnoreCase(CONSTANT_POPULATION)) {
                    int index = 1;
                    data.demographicModelIndex = index;
                    if (arguments.hasOption(CONSTANT_POPULATION_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(CONSTANT_POPULATION_PARAMETER_VALUES);
                        parseDemographicValues(index, values, data);
                    }
                } else if (option.equalsIgnoreCase(EXPONENTIAL_GROWTH_RATE)) {
                    int index = 2;
                    data.demographicModelIndex = index;
                    if (arguments.hasOption(EXPONENTIAL_GROWTH_RATE_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(EXPONENTIAL_GROWTH_RATE_PARAMETER_VALUES);
                        parseDemographicValues(index, values, data);
                    }
                } else if (option.equalsIgnoreCase(EXPONENTIAL_DOUBLING_TIME)) {
                    int index = 3;
                    data.demographicModelIndex = index;
                    if (arguments.hasOption(EXPONENTIAL_GROWTH_DOUBLING_TIME_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(EXPONENTIAL_GROWTH_DOUBLING_TIME_PARAMETER_VALUES);
                        parseDemographicValues(index, values, data);
                    }
                } else {
                    gracefullyExit("Unrecognized option.");
                }
            }
            // Branch Substitution Model
            if (arguments.hasOption(BRANCH_SUBSTITUTION_MODEL)) {
                option = arguments.getStringOption(BRANCH_SUBSTITUTION_MODEL);
                if (option.equalsIgnoreCase(HKY)) {
                    int index = 0;
                    data.substitutionModelIndex = index;
                    if (arguments.hasOption(HKY_SUBSTITUTION_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(HKY_SUBSTITUTION_PARAMETER_VALUES);
                        parseSubstitutionValues(index, values, data);
                    }
                } else if (option.equalsIgnoreCase(GTR)) {
                    int index = 1;
                    data.substitutionModelIndex = index;
                    if (arguments.hasOption(GTR_SUBSTITUTION_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(GTR_SUBSTITUTION_PARAMETER_VALUES);
                        parseSubstitutionValues(index, values, data);
                    }
                } else if (option.equalsIgnoreCase(TN93)) {
                    int index = 2;
                    data.substitutionModelIndex = index;
                    if (arguments.hasOption(TN93_SUBSTITUTION_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(TN93_SUBSTITUTION_PARAMETER_VALUES);
                        parseSubstitutionValues(index, values, data);
                    }
                } else if (option.equalsIgnoreCase(GY94_CODON_MODEL)) {
                    int index = 3;
                    data.substitutionModelIndex = index;
                    if (arguments.hasOption(GY94_SUBSTITUTION_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(GY94_SUBSTITUTION_PARAMETER_VALUES);
                        parseSubstitutionValues(index, values, data);
                    }
                } else {
                    gracefullyExit("Unrecognized option.");
                }
            }
            // Site Rate Model
            if (arguments.hasOption(SITE_RATE_MODEL)) {
                option = arguments.getStringOption(SITE_RATE_MODEL);
                if (option.equalsIgnoreCase(NO_SITE_RATE_MODEL)) {
                    int index = 0;
                    data.siteRateModelIndex = index;
                } else if (option.equalsIgnoreCase(GAMMA_SITE_RATE_MODEL)) {
                    int index = 1;
                    data.siteRateModelIndex = index;
                    if (arguments.hasOption(GAMMA_SITE_RATE_MODEL_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(GAMMA_SITE_RATE_MODEL_PARAMETER_VALUES);
                        parseSiteRateValues(index, values, data);
                    }
                } else {
                    gracefullyExit("Unrecognized option.");
                }
            }
            // Clock Rate Model
            if (arguments.hasOption(CLOCK_RATE_MODEL)) {
                option = arguments.getStringOption(CLOCK_RATE_MODEL);
                if (option.equalsIgnoreCase(STRICT_CLOCK)) {
                    int index = 0;
                    data.clockModelIndex = index;
                    if (arguments.hasOption(STRICT_CLOCK_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(STRICT_CLOCK_PARAMETER_VALUES);
                        parseClockValues(index, values, data);
                    }
                } else if (option.equalsIgnoreCase(LOGNORMAL_RELAXED_CLOCK)) {
                    int index = 1;
                    data.clockModelIndex = index;
                    if (arguments.hasOption(LOGNORMAL_RELAXED_CLOCK_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(LOGNORMAL_RELAXED_CLOCK_PARAMETER_VALUES);
                        parseClockValues(index, values, data);
                    }
                    // set to true/false
                    if (arguments.hasOption(LRC_PARAMETERS_IN_REAL_SPACE)) {
                        boolean lrcParametersInRealSpace = (arguments.getStringOption(LRC_PARAMETERS_IN_REAL_SPACE).equalsIgnoreCase(TRUE) ? true : false);
                        data.lrcParametersInRealSpace = lrcParametersInRealSpace;
                    }
                } else if (option.equalsIgnoreCase(EXPONENTIAL_RELAXED_CLOCK)) {
                    int index = 2;
                    data.clockModelIndex = index;
                    if (arguments.hasOption(EXPONENTIAL_RELAXED_CLOCK_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(EXPONENTIAL_RELAXED_CLOCK_PARAMETER_VALUES);
                        parseClockValues(index, values, data);
                    }
                } else if (option.equalsIgnoreCase(INVERSE_GAUSSIAN_RELAXED_CLOCK)) {
                    int index = 3;
                    data.clockModelIndex = index;
                    if (arguments.hasOption(INVERSE_GAUSSIAN_RELAXED_CLOCK_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(INVERSE_GAUSSIAN_RELAXED_CLOCK_PARAMETER_VALUES);
                        parseClockValues(index, values, data);
                    }
                } else {
                    gracefullyExit("Unrecognized option.");
                }
            }
            // Frequency Model
            if (arguments.hasOption(BASE_FREQUENCIES)) {
                option = arguments.getStringOption(BASE_FREQUENCIES);
                if (option.equalsIgnoreCase(NUCLEOTIDE_FREQUENCIES)) {
                    int index = 0;
                    data.frequencyModelIndex = index;
                    if (arguments.hasOption(NUCLEOTIDE_FREQUENCY_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(NUCLEOTIDE_FREQUENCY_PARAMETER_VALUES);
                        parseFrequencyValues(index, values, data);
                    }
                } else if (option.equalsIgnoreCase(CODON_FREQUENCIES)) {
                    int index = 1;
                    data.frequencyModelIndex = index;
                    if (arguments.hasOption(CODON_FREQUENCY_PARAMETER_VALUES)) {
                        values = arguments.getRealArrayOption(CODON_FREQUENCY_PARAMETER_VALUES);
                        parseFrequencyValues(index, values, data);
                    }
                } else {
                    gracefullyExit("Unrecognized option.");
                }
            }
            if (arguments.hasOption(FROM)) {
                data.from = arguments.getIntegerOption(FROM);
            }
            if (arguments.hasOption(TO)) {
                data.to = arguments.getIntegerOption(TO);
            }
            if (arguments.hasOption(EVERY)) {
                data.every = arguments.getIntegerOption(EVERY);
            }
            // END: EVERY option check
            // create partition
            Partition partition = new Partition(//
            data.createTreeModel(), //
            data.createBranchModel(), //
            data.createSiteRateModel(), //
            data.createClockRateModel(), //
            data.createFrequencyModel(), // from
            data.from - 1, // to
            data.to - 1, // every
            data.every);
            if (arguments.hasOption(ROOT_SEQUENCE)) {
                data.ancestralSequenceString = arguments.getStringOption(ROOT_SEQUENCE);
                partition.setRootSequence(data.createAncestralSequence());
            }
            // END: ANCESTRAL_SEQUENCE option check
            partitionsList.add(partition);
            //                System.err.println(data.from);
            dataList.add(data);
        }
        if (this.VERBOSE) {
            //        		System.out.println(dataList.get(0).from + " " + dataList.get(1).from);
            Utils.printPartitionDataList(dataList);
            System.out.println();
        }
        SimpleAlignment alignment = new SimpleAlignment();
        String outputFile = null;
        if (leftoverArguments.length > 0) {
            outputFile = leftoverArguments[0];
            String[] file = outputFile.split("\\.", 2);
            if (file.length > 1) {
                String extension = outputFile.split("\\.", 2)[1];
                // TODO Delegate here to enum-class; switches are not generic
                if (extension.equalsIgnoreCase(SimpleAlignment.OutputType.FASTA.getText()) || extension.equalsIgnoreCase("fst")) {
                    dataList.outputFormat = SimpleAlignment.OutputType.FASTA;
                } else if (extension.equalsIgnoreCase(SimpleAlignment.OutputType.NEXUS.getText()) || extension.equalsIgnoreCase("nxs")) {
                    dataList.outputFormat = SimpleAlignment.OutputType.NEXUS;
                } else if (extension.equalsIgnoreCase(SimpleAlignment.OutputType.XML.getText())) {
                    dataList.outputFormat = SimpleAlignment.OutputType.XML;
                } else {
                    dataList.outputFormat = SimpleAlignment.OutputType.FASTA;
                }
            //END: extension check
            } else {
                outputFile = file[0] + "." + SimpleAlignment.OutputType.FASTA.toString().toLowerCase();
                dataList.outputFormat = SimpleAlignment.OutputType.FASTA;
            }
        //END: 
        }
        if (leftoverArguments.length > 1) {
            dataList.startingSeed = Long.parseLong(leftoverArguments[1]);
            dataList.setSeed = true;
        }
        if (dataList.setSeed) {
            MathUtils.setSeed(dataList.startingSeed);
        }
        if (leftoverArguments.length > 2) {
            dataList.outputAncestralSequences = Boolean.parseBoolean(leftoverArguments[2]);
        }
        if (leftoverArguments.length > 3) {
            dataList.useParallel = Boolean.parseBoolean(leftoverArguments[3]);
        }
        BeagleSequenceSimulator beagleSequenceSimulator = new BeagleSequenceSimulator(partitionsList);
        alignment = beagleSequenceSimulator.simulate(dataList.useParallel, dataList.outputAncestralSequences);
        alignment.setOutputType(dataList.outputFormat);
        PrintWriter writer = new PrintWriter(new FileWriter(outputFile));
        writer.println(alignment.toString());
        writer.close();
    } catch (ArgumentException e) {
        System.out.println();
        printUsage(arguments);
        System.out.println();
        System.out.println(e.getMessage());
        e.printStackTrace();
        System.exit(1);
    } catch (IOException e) {
        System.out.println();
        printUsage(arguments);
        System.out.println();
        System.out.println(e.getMessage());
        e.printStackTrace();
        System.exit(1);
    } catch (ImportException e) {
        System.out.println();
        printUsage(arguments);
        System.out.println();
        System.out.println(e.getMessage());
        e.printStackTrace();
        System.exit(1);
    }
// END: try-catch block
}

Example 52 with Taxa

use of dr.evolution.util.Taxa 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 53 with Taxa

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

the class TaxaEditor method doOk.

// END: ListenOk
private void doOk() {
    frame.setBusy();
    SwingWorker<Void, Void> worker = new SwingWorker<Void, Void>() {

        // Executed in background thread
        public Void doInBackground() {
            try {
                // delete taxa connected to this row
                String value = dataList.recordsList.get(row).getName();
                Utils.removeTaxaWithAttributeValue(dataList, Utils.TREE_FILENAME, value);
                String name = String.valueOf("TaxaSet").concat(String.valueOf(row + 1));
                Taxa taxa = taxaEditorTableModel.getTaxaSet();
                TreesTableRecord record = new TreesTableRecord(name, taxa);
                dataList.recordsList.set(row, record);
                dataList.allTaxa.addTaxa(taxa);
            // treesTableModel.setRow(row, record);
            } catch (Exception e) {
                Utils.handleException(e);
            }
            return null;
        }

        // END: doInBackground()
        // Executed in event dispatch thread
        public void done() {
            frame.setIdle();
            frame.fireTaxaChanged();
            window.setVisible(false);
        }
    };
    worker.execute();
}

Example 54 with Taxa

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

the class TaxaEditor method loadTaxaFile.

// END: ListenLoadTaxaFile
private void loadTaxaFile(File file) {
    try {
        Taxa taxa = new Taxa();
        String[] lines = Utils.loadStrings(file.getAbsolutePath());
        Taxon taxon;
        for (int i = 0; i < lines.length; i++) {
            String[] line = lines[i].split("\\s+");
            taxon = new Taxon(line[TaxaEditorTableModel.NAME_INDEX]);
            taxon.setAttribute(Utils.ABSOLUTE_HEIGHT, Double.valueOf(line[TaxaEditorTableModel.HEIGHT_INDEX]));
            taxa.addTaxon(taxon);
        }
        // END: i loop
        updateTable(taxa);
    } catch (Exception e) {
        Utils.handleException(e);
    }
// END: try-catch block
}

Example 55 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");
}