Examples with Partition dr.app.beagle.tools.Partition used on opensource projects

Example 11 with Partition

use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateOnePartition.

// END: simulate topology
static void simulateOnePartition() {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 2: simulateOnePartition");
        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);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        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);
        Sequence ancestralSequence = new Sequence();
        ancestralSequence.appendSequenceString("TCAAGTGAGG");
        partition1.setRootSequence(ancestralSequence);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        SimpleAlignment alignment = simulator.simulate(simulateInPar, false);
        // alignment.setOutputType(SimpleAlignment.OutputType.NEXUS);
        alignment.setOutputType(SimpleAlignment.OutputType.XML);
        System.out.println(alignment.toString());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 12 with Partition

use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.

the class BeagleBranchLikelihood method main.

// END: LikelihoodColumn class
// ////////////
// ---TEST---//
// ////////////
public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("((SimSeq1:22.0,SimSeq2:22.0):12.0,(SimSeq3:23.1,SimSeq4:23.1):10.899999999999999);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 1.000);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogeneousBranchModel, //
        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);
        Alignment alignment = simulator.simulate(false, false);
        System.out.println(alignment);
        BeagleTreeLikelihood btl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("BTL(homogeneous) = " + btl.getLogLikelihood());
        BeagleBranchLikelihood bbl = new BeagleBranchLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, freqModel, branchRateModel);
        int branchIndex = 4;
        System.out.println(bbl.getBranchLogLikelihood(branchIndex));
        bbl.finalizeBeagle();
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    } catch (Throwable e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 13 with Partition

use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.

the class LineageSpecificBranchModel method main.

// END: acceptState
public static void main(String[] args) {
    try {
        // the seed of the BEAST
        MathUtils.setSeed(666);
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        TreeModel tree = new TreeModel(importer.importTree(null));
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { //
        0.0163936, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344 });
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create substitution model
        Parameter alpha = new Parameter.Default(1, 10);
        Parameter beta = new Parameter.Default(1, 5);
        MG94CodonModel mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(mg94);
        // create partition
        Partition partition1 = new //
        Partition(//
        tree, //
        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);
        Alignment alignment = simulator.simulate(false, false);
        ConvertAlignment convert = new ConvertAlignment(Nucleotides.INSTANCE, GeneticCode.UNIVERSAL, alignment);
        List<SubstitutionModel> substModels = new ArrayList<SubstitutionModel>();
        for (int i = 0; i < 2; i++) {
            //				alpha = new Parameter.Default(1, 10 );
            //				beta = new Parameter.Default(1, 5 );
            //				mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta,
            //						freqModel);
            substModels.add(mg94);
        }
        Parameter uCategories = new Parameter.Default(2, 0);
        //            CountableBranchCategoryProvider provider = new CountableBranchCategoryProvider.IndependentBranchCategoryModel(tree, uCategories);
        LineageSpecificBranchModel branchSpecific = new //provider, 
        LineageSpecificBranchModel(//provider, 
        tree, //provider, 
        freqModel, //provider, 
        substModels, uCategories);
        BeagleTreeLikelihood like = new //
        BeagleTreeLikelihood(//
        convert, //
        tree, //
        branchSpecific, //
        siteRateModel, //
        branchRateModel, //
        null, //
        false, PartialsRescalingScheme.DEFAULT, true);
        BeagleTreeLikelihood gold = new //
        BeagleTreeLikelihood(//
        convert, //
        tree, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        null, //
        false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("likelihood (gold) = " + gold.getLogLikelihood());
        System.out.println("likelihood = " + like.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
    }
}

Example 14 with Partition

use of dr.app.beagle.tools.Partition 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
}