Examples with TreeModel dr.evomodel.tree.TreeModel used on opensource projects

Example 21 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateAminoAcid.

// END: simulateThreePartitions
static void simulateAminoAcid() {
    try {
        System.out.println("Test case 4: simulateAminoAcid");
        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);
        TreeModel treeModel = new TreeModel(tree);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05 });
        FrequencyModel freqModel = new FrequencyModel(AminoAcids.INSTANCE, freqs);
        // create substitution model
        EmpiricalRateMatrix rateMatrix = Blosum62.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
        // 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
}

Example 22 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateRandomBranchAssignment.

// END: main
static void simulateRandomBranchAssignment() {
    try {
        System.out.println("Test case I dunno which: simulate random branch assignments");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        File treeFile = new File("/home/filip/Dropbox/BeagleSequenceSimulator/SimTree/SimTree.figtree");
        Tree tree = Utils.importTreeFromFile(treeFile);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(Utils.UNIFORM_CODON_FREQUENCIES);
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create base subst model
        Parameter omegaParameter = new Parameter.Default("omega", 1, 1.0);
        Parameter kappaParameter = new Parameter.Default("kappa", 1, 1.0);
        GY94CodonModel baseSubModel = new GY94CodonModel(Codons.UNIVERSAL, omegaParameter, kappaParameter, freqModel);
        RandomBranchModel substitutionModel = new RandomBranchModel(treeModel, baseSubModel, 0.25, false, -1);
        // 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, true);
        // alignment.setOutputType(SimpleAlignment.OutputType.NEXUS);
        alignment.setOutputType(SimpleAlignment.OutputType.XML);
        System.out.println(alignment.toString());
    } catch (Exception e) {
        e.printStackTrace();
    }
// END: try-catch
}

Example 23 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class BeastCheckpointer method readStateFromFile.

private long readStateFromFile(File file, MarkovChain markovChain, double[] lnL) {
    OperatorSchedule operatorSchedule = markovChain.getSchedule();
    long state = -1;
    ArrayList<TreeParameterModel> traitModels = new ArrayList<TreeParameterModel>();
    try {
        FileReader fileIn = new FileReader(file);
        BufferedReader in = new BufferedReader(fileIn);
        int[] rngState = null;
        String line = in.readLine();
        String[] fields = line.split("\t");
        if (fields[0].equals("rng")) {
            // if there is a random number generator state present then load it...
            try {
                rngState = new int[fields.length - 1];
                for (int i = 0; i < rngState.length; i++) {
                    rngState[i] = Integer.parseInt(fields[i + 1]);
                }
            } catch (NumberFormatException nfe) {
                throw new RuntimeException("Unable to read state number from state file");
            }
            line = in.readLine();
            fields = line.split("\t");
        }
        try {
            if (!fields[0].equals("state")) {
                throw new RuntimeException("Unable to read state number from state file");
            }
            state = Long.parseLong(fields[1]);
        } catch (NumberFormatException nfe) {
            throw new RuntimeException("Unable to read state number from state file");
        }
        line = in.readLine();
        fields = line.split("\t");
        try {
            if (!fields[0].equals("lnL")) {
                throw new RuntimeException("Unable to read lnL from state file");
            }
            if (lnL != null) {
                lnL[0] = Double.parseDouble(fields[1]);
            }
        } catch (NumberFormatException nfe) {
            throw new RuntimeException("Unable to read lnL from state file");
        }
        for (Parameter parameter : Parameter.CONNECTED_PARAMETER_SET) {
            line = in.readLine();
            fields = line.split("\t");
            //if (!fields[0].equals(parameter.getParameterName())) {
            //  System.err.println("Unable to match state parameter: " + fields[0] + ", expecting " + parameter.getParameterName());
            //}
            int dimension = Integer.parseInt(fields[2]);
            if (dimension != parameter.getDimension()) {
                System.err.println("Unable to match state parameter dimension: " + dimension + ", expecting " + parameter.getDimension() + " for parameter: " + parameter.getParameterName());
                System.err.print("Read from file: ");
                for (int i = 0; i < fields.length; i++) {
                    System.err.print(fields[i] + "\t");
                }
                System.err.println();
            }
            if (fields[1].equals("branchRates.categories.rootNodeNumber")) {
                // System.out.println("eek");
                double value = Double.parseDouble(fields[3]);
                parameter.setParameterValue(0, value);
                if (DEBUG) {
                    System.out.println("restoring " + fields[1] + " with value " + value);
                }
            } else {
                if (DEBUG) {
                    System.out.print("restoring " + fields[1] + " with values ");
                }
                for (int dim = 0; dim < parameter.getDimension(); dim++) {
                    parameter.setParameterValue(dim, Double.parseDouble(fields[dim + 3]));
                    if (DEBUG) {
                        System.out.print(Double.parseDouble(fields[dim + 3]) + " ");
                    }
                }
                if (DEBUG) {
                    System.out.println();
                }
            }
        }
        for (int i = 0; i < operatorSchedule.getOperatorCount(); i++) {
            MCMCOperator operator = operatorSchedule.getOperator(i);
            line = in.readLine();
            fields = line.split("\t");
            if (!fields[1].equals(operator.getOperatorName())) {
                throw new RuntimeException("Unable to match operator: " + fields[1]);
            }
            if (fields.length < 4) {
                throw new RuntimeException("Operator missing values: " + fields[1]);
            }
            operator.setAcceptCount(Integer.parseInt(fields[2]));
            operator.setRejectCount(Integer.parseInt(fields[3]));
            if (operator instanceof CoercableMCMCOperator) {
                if (fields.length != 5) {
                    throw new RuntimeException("Coercable operator missing parameter: " + fields[1]);
                }
                ((CoercableMCMCOperator) operator).setCoercableParameter(Double.parseDouble(fields[4]));
            }
        }
        // load the tree models last as we get the node heights from the tree (not the parameters which
        // which may not be associated with the right node
        Set<String> expectedTreeModelNames = new HashSet<String>();
        for (Model model : Model.CONNECTED_MODEL_SET) {
            if (model instanceof TreeModel) {
                if (DEBUG) {
                    System.out.println("model " + model.getModelName());
                }
                expectedTreeModelNames.add(model.getModelName());
                if (DEBUG) {
                    for (String s : expectedTreeModelNames) {
                        System.out.println(s);
                    }
                }
            }
            if (model instanceof TreeParameterModel) {
                traitModels.add((TreeParameterModel) model);
            }
        }
        line = in.readLine();
        fields = line.split("\t");
        // Read in all (possibly more than one) trees
        while (fields[0].equals("tree")) {
            if (DEBUG) {
                System.out.println("tree: " + fields[1]);
            }
            for (Model model : Model.CONNECTED_MODEL_SET) {
                if (model instanceof TreeModel && fields[1].equals(model.getModelName())) {
                    line = in.readLine();
                    line = in.readLine();
                    fields = line.split("\t");
                    //read number of nodes
                    int nodeCount = Integer.parseInt(fields[0]);
                    double[] nodeHeights = new double[nodeCount];
                    for (int i = 0; i < nodeCount; i++) {
                        line = in.readLine();
                        fields = line.split("\t");
                        nodeHeights[i] = Double.parseDouble(fields[1]);
                    }
                    //on to reading edge information
                    line = in.readLine();
                    line = in.readLine();
                    line = in.readLine();
                    fields = line.split("\t");
                    int edgeCount = Integer.parseInt(fields[0]);
                    //create data matrix of doubles to store information from list of TreeParameterModels
                    double[][] traitValues = new double[traitModels.size()][edgeCount];
                    //create array to store whether a node is left or right child of its parent
                    //can be important for certain tree transition kernels
                    int[] childOrder = new int[edgeCount];
                    for (int i = 0; i < childOrder.length; i++) {
                        childOrder[i] = -1;
                    }
                    int[] parents = new int[edgeCount];
                    for (int i = 0; i < edgeCount; i++) {
                        parents[i] = -1;
                    }
                    for (int i = 0; i < edgeCount; i++) {
                        line = in.readLine();
                        if (line != null) {
                            fields = line.split("\t");
                            parents[Integer.parseInt(fields[0])] = Integer.parseInt(fields[1]);
                            childOrder[i] = Integer.parseInt(fields[2]);
                            for (int j = 0; j < traitModels.size(); j++) {
                                traitValues[j][i] = Double.parseDouble(fields[3 + j]);
                            }
                        }
                    }
                    //perform magic with the acquired information
                    if (DEBUG) {
                        System.out.println("adopting tree structure");
                    }
                    //adopt the loaded tree structure; this does not yet copy the traits on the branches
                    ((TreeModel) model).beginTreeEdit();
                    ((TreeModel) model).adoptTreeStructure(parents, nodeHeights, childOrder);
                    ((TreeModel) model).endTreeEdit();
                    expectedTreeModelNames.remove(model.getModelName());
                }
            }
            line = in.readLine();
            if (line != null) {
                fields = line.split("\t");
            }
        }
        if (expectedTreeModelNames.size() > 0) {
            StringBuilder sb = new StringBuilder();
            for (String notFoundName : expectedTreeModelNames) {
                sb.append("Expecting, but unable to match state parameter:" + notFoundName + "\n");
            }
            throw new RuntimeException(sb.toString());
        }
        if (DEBUG) {
            System.out.println("\nDouble checking:");
            for (Parameter parameter : Parameter.CONNECTED_PARAMETER_SET) {
                if (parameter.getParameterName().equals("branchRates.categories.rootNodeNumber")) {
                    System.out.println(parameter.getParameterName() + ": " + parameter.getParameterValue(0));
                }
            }
        }
        if (rngState != null) {
            MathUtils.setRandomState(rngState);
        }
        in.close();
        fileIn.close();
    // This shouldn't be necessary and if it is then it might be hiding a bug...
    //            for (Likelihood likelihood : Likelihood.CONNECTED_LIKELIHOOD_SET) {
    //                likelihood.makeDirty();
    //            }
    } catch (IOException ioe) {
        throw new RuntimeException("Unable to read file: " + ioe.getMessage());
    }
    return state;
}

Example 24 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class Utils method partitionDataListToString.

// END: partitionDataToString
public static //
String partitionDataListToString(//
PartitionDataList dataList, ArrayList<TreeModel> simulatedTreeModelList) // ,LinkedHashMap<Integer,LinkedHashMap<NodeRef, int[]>>
// partitionSequencesMap
{
    String string = "";
    TreeModel simulatedTreeModel;
    // LinkedHashMap<NodeRef, int[]> sequencesMap;
    string += ("Site count: " + getSiteCount(dataList)) + ("\n");
    if (dataList.setSeed) {
        string += ("Starting seed: " + dataList.startingSeed) + ("\n");
    }
    int row = 0;
    for (PartitionData data : dataList) {
        simulatedTreeModel = simulatedTreeModelList.get(row);
        // sequencesMap = partitionSequencesMap.get(row);
        string += ("Partition: " + (row + 1)) + ("\n");
        string += partitionDataToString(data, simulatedTreeModel);
        string += ("\n");
        row++;
    }
    return string;
}

Example 25 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class RandomWalkIntegerSetSizeWeightedOperator method computeSampleWeights.

private void computeSampleWeights() {
    TreeModel tree = msatSampleTreeModel.getTreeModel();
    int intNodeCount = tree.getInternalNodeCount();
    int extNodeCount = tree.getExternalNodeCount();
    weights = new double[intNodeCount];
    for (int i = 0; i < intNodeCount; i++) {
        NodeRef node = tree.getNode(i + extNodeCount);
        int lcState = msatSampleTreeModel.getNodeValue(tree.getChild(node, 0));
        int rcState = msatSampleTreeModel.getNodeValue(tree.getChild(node, 1));
        weights[i] = Math.abs(lcState - rcState) + baseSetSize;
    }
}