Examples with NewickImporter dr.evolution.io.NewickImporter used on opensource projects

Example 1 with NewickImporter

use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.

the class BeagleTreeLikelihood method main.

public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 1: simulateOnePartition");
        int sequenceLength = 1000;
        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 branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        Parameter kappa2 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        HKY hky2 = new HKY(kappa2, freqModel);
        HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        substitutionModels.add(hky2);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        Parameter epochTimes = new Parameter.Default(1, 20);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 0.001);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogenousBranchSubstitutionModel, //
        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);
        BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
        nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 2 with NewickImporter

use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.

the class TreeTraceAnalysisParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    try {
        Reader reader;
        String fileName = xo.getStringAttribute(FILE_NAME);
        String name;
        try {
            File file = new File(fileName);
            name = file.getName();
            String parent = file.getParent();
            if (!file.isAbsolute()) {
                parent = System.getProperty("user.dir");
            }
            //					System.out.println("Writing log file to "+parent+System.getProperty("path.separator")+name);
            reader = new FileReader(new File(parent, name));
        } catch (FileNotFoundException fnfe) {
            throw new XMLParseException("File '" + fileName + "' can not be opened for " + getParserName() + " element.");
        }
        // leaving the burnin attribute off will result in 10% being used
        int burnin = xo.getAttribute(BURN_IN, -1);
        double minCladeProbability = xo.getAttribute(MIN_CLADE_PROBABILITY, 0.5);
        double credSetProbability = xo.getAttribute(CRED_SET_PROBABILITY, 0.95);
        Tree referenceTree = null;
        Reader refReader;
        if (xo.hasAttribute(REFERENCE_TREE)) {
            String referenceName = xo.getStringAttribute(REFERENCE_TREE);
            try {
                File refFile = new File(referenceName);
                String refName = refFile.getName();
                String parent = refFile.getParent();
                if (!refFile.isAbsolute()) {
                    parent = System.getProperty("user.dir");
                }
                refReader = new FileReader(new File(parent, refName));
            } catch (FileNotFoundException fnfe) {
                throw new XMLParseException("File '" + fileName + "' can not be opened for " + getParserName() + " element.");
            }
            try {
                NewickImporter importTree = new NewickImporter(refReader);
                if (importTree.hasTree()) {
                    referenceTree = importTree.importNextTree();
                }
            } catch (Importer.ImportException iee) {
                throw new XMLParseException("Reference file '" + referenceName + "' is empty.");
            }
        }
        boolean shortReport = xo.getAttribute(SHORT_REPORT, false);
        TreeTraceAnalysis analysis = TreeTraceAnalysis.analyzeLogFile(new Reader[] { reader }, burnin, true);
        if (shortReport) {
            analysis.shortReport(name, referenceTree, true, credSetProbability);
        } else {
            analysis.report(minCladeProbability, credSetProbability, 0);
        }
        System.out.println();
        System.out.flush();
        return analysis;
    } catch (java.io.IOException ioe) {
        throw new XMLParseException(ioe.getMessage());
    }
}

Example 3 with NewickImporter

use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.

the class Branch2dRateToGrid method readAndAnalyzeTrees.

private void readAndAnalyzeTrees(String treeFileName, int burnin, int skipEvery, String traitString, boolean traitNoise, boolean rateNoise, double maxPathLength, Normalization normalize, boolean getStdev) throws IOException, Importer.ImportException {
    int totalTrees = 10000;
    int totalStars = 0;
    if (!printedBar) {
        progressStream.println("Reading and analyzing trees (bar assumes 10,000 trees)...");
        progressStream.println("0              25             50             75            100");
        progressStream.println("|--------------|--------------|--------------|--------------|");
        printedBar = true;
    }
    if (getStdev) {
        progressStream.println("summarizing standard deviations");
    }
    int stepSize = totalTrees / 60;
    if (stepSize < 1)
        stepSize = 1;
    BufferedReader reader1 = new BufferedReader(new FileReader(treeFileName));
    String line1 = reader1.readLine();
    TreeImporter importer1;
    if (line1.toUpperCase().startsWith("#NEXUS")) {
        importer1 = new NexusImporter(new FileReader(treeFileName));
    } else {
        importer1 = new NewickImporter(new FileReader(treeFileName));
    }
    totalTrees = 0;
    while (importer1.hasTree()) {
        Tree treeTime = importer1.importNextTree();
        if (totalTrees % skipEvery == 0) {
            treesRead++;
            if (totalTrees >= burnin) {
                analyzeTree(treeTime, traitString, traitNoise, rateNoise, maxPathLength, normalize, getStdev);
            }
        }
        if (totalTrees > 0 && totalTrees % stepSize == 0) {
            progressStream.print("*");
            totalStars++;
            if (totalStars % 61 == 0)
                progressStream.print("\n");
            progressStream.flush();
        }
        totalTrees++;
    }
    progressStream.print("\n");
}

Example 4 with NewickImporter

use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.

the class BranchRatePlotter method main.

public static void main(String[] args) throws java.io.IOException, Importer.ImportException {
    String controlFile = args[0];
    //String treeFile1 = args[0];
    //String treeFile2 = args[1];
    String targetTreeFile = args[1];
    int burnin = 0;
    if (args.length > 2) {
        burnin = Integer.parseInt(args[2]);
    }
    System.out.println("Ignoring first " + burnin + " trees as burnin.");
    BufferedReader readerTarget = new BufferedReader(new FileReader(targetTreeFile));
    String lineTarget = readerTarget.readLine();
    readerTarget.close();
    TreeImporter targetImporter;
    if (lineTarget.toUpperCase().startsWith("#NEXUS")) {
        targetImporter = new NexusImporter(new FileReader(targetTreeFile));
    } else {
        targetImporter = new NewickImporter(new FileReader(targetTreeFile));
    }
    MutableTree targetTree = new FlexibleTree(targetImporter.importNextTree());
    targetTree = TreeUtils.rotateTreeByComparator(targetTree, TreeUtils.createNodeDensityComparator(targetTree));
    BufferedReader reader = new BufferedReader(new FileReader(controlFile));
    String line = reader.readLine();
    int totalTrees = 0;
    int totalTreesUsed = 0;
    while (line != null) {
        StringTokenizer tokens = new StringTokenizer(line);
        NexusImporter importer1 = new NexusImporter(new FileReader(tokens.nextToken()));
        NexusImporter importer2 = new NexusImporter(new FileReader(tokens.nextToken()));
        int fileTotalTrees = 0;
        while (importer1.hasTree()) {
            Tree timeTree = importer1.importNextTree();
            Tree mutationTree = importer2.importNextTree();
            if (fileTotalTrees >= burnin) {
                annotateRates(targetTree, targetTree.getRoot(), timeTree, mutationTree);
                totalTreesUsed += 1;
            }
            totalTrees += 1;
            fileTotalTrees += 1;
        }
        line = reader.readLine();
    }
    System.out.println("Total trees read: " + totalTrees);
    System.out.println("Total trees summarized: " + totalTreesUsed);
    // collect all rates
    double mutations = 0.0;
    double time = 0.0;
    double[] rates = new double[targetTree.getNodeCount() - 1];
    int index = 0;
    for (int i = 0; i < targetTree.getNodeCount(); i++) {
        NodeRef node = targetTree.getNode(i);
        if (!targetTree.isRoot(node)) {
            Integer count = ((Integer) targetTree.getNodeAttribute(node, "count"));
            if (count == null) {
                throw new RuntimeException("Count missing from node in target tree");
            }
            if (!targetTree.isExternal(node)) {
                double prob = (double) (int) count / (double) (totalTreesUsed);
                if (prob >= 0.5) {
                    String label = "" + (Math.round(prob * 100) / 100.0);
                    targetTree.setNodeAttribute(node, "label", label);
                }
            }
            Number totalMutations = (Number) targetTree.getNodeAttribute(node, "totalMutations");
            Number totalTime = (Number) targetTree.getNodeAttribute(node, "totalTime");
            mutations += totalMutations.doubleValue();
            time += totalTime.doubleValue();
            rates[index] = totalMutations.doubleValue() / totalTime.doubleValue();
            System.out.println(totalMutations.doubleValue() + " / " + totalTime.doubleValue() + " = " + rates[index]);
            targetTree.setNodeRate(node, rates[index]);
            index += 1;
        }
    }
    double minRate = DiscreteStatistics.min(rates);
    double maxRate = DiscreteStatistics.max(rates);
    double medianRate = DiscreteStatistics.median(rates);
    //double topThird = DiscreteStatistics.quantile(2.0/3.0,rates);
    //double bottomThird = DiscreteStatistics.quantile(1.0/3.0,rates);
    //double unweightedMeanRate = DiscreteStatistics.mean(rates);
    double meanRate = mutations / time;
    System.out.println(minRate + "\t" + maxRate + "\t" + medianRate + "\t" + meanRate);
    for (int i = 0; i < targetTree.getNodeCount(); i++) {
        NodeRef node = targetTree.getNode(i);
        if (!targetTree.isRoot(node)) {
            double rate = targetTree.getNodeRate(node);
            //double branchTime = ((Number)targetTree.getNodeAttribute(node, "totalTime")).doubleValue();
            //double branchMutations = ((Number)targetTree.getNodeAttribute(node, "totalMutations")).doubleValue();
            float relativeRate = (float) (rate / maxRate);
            float radius = (float) Math.sqrt(relativeRate * 36.0);
            if (rate > meanRate) {
                targetTree.setNodeAttribute(node, "color", new Color(1.0f, 0.5f, 0.5f));
            } else {
                targetTree.setNodeAttribute(node, "color", new Color(0.5f, 0.5f, 1.0f));
            }
            //targetTree.setNodeAttribute(node, "color", new Color(red, green, blue));
            targetTree.setNodeAttribute(node, "line", new BasicStroke(1.0f));
            targetTree.setNodeAttribute(node, "shape", new java.awt.geom.Ellipse2D.Double(0, 0, radius * 2.0, radius * 2.0));
        }
        java.util.List heightList = (java.util.List) targetTree.getNodeAttribute(node, "heightList");
        if (heightList != null) {
            double[] heights = new double[heightList.size()];
            for (int j = 0; j < heights.length; j++) {
                heights[j] = (Double) heightList.get(j);
            }
            targetTree.setNodeHeight(node, DiscreteStatistics.mean(heights));
            //if (heights.length >= (totalTreesUsed/2)) {
            targetTree.setNodeAttribute(node, "nodeHeight.mean", DiscreteStatistics.mean(heights));
            targetTree.setNodeAttribute(node, "nodeHeight.hpdUpper", DiscreteStatistics.quantile(0.975, heights));
            targetTree.setNodeAttribute(node, "nodeHeight.hpdLower", DiscreteStatistics.quantile(0.025, heights));
        //targetTree.setNodeAttribute(node, "nodeHeight.max", new Double(DiscreteStatistics.max(heights)));
        //targetTree.setNodeAttribute(node, "nodeHeight.min", new Double(DiscreteStatistics.min(heights)));
        //}
        }
    }
    StringBuffer buffer = new StringBuffer();
    writeTree(targetTree, targetTree.getRoot(), buffer, true, false);
    buffer.append(";\n");
    writeTree(targetTree, targetTree.getRoot(), buffer, false, true);
    buffer.append(";\n");
    System.out.println(buffer.toString());
    SquareTreePainter treePainter = new SquareTreePainter();
    treePainter.setColorAttribute("color");
    treePainter.setLineAttribute("line");
    //        treePainter.setShapeAttribute("shape");
    //        treePainter.setLabelAttribute("label");
    JTreeDisplay treeDisplay = new JTreeDisplay(treePainter, targetTree);
    JTreePanel treePanel = new JTreePanel(treeDisplay);
    JFrame frame = new JFrame();
    frame.setSize(800, 600);
    frame.getContentPane().setLayout(new BorderLayout());
    frame.getContentPane().add(treePanel);
    frame.setVisible(true);
    PrinterJob printJob = PrinterJob.getPrinterJob();
    printJob.setPrintable(treeDisplay);
    if (printJob.printDialog()) {
        try {
            printJob.print();
        } catch (Exception ex) {
            throw new RuntimeException(ex);
        }
    }
}

Example 5 with NewickImporter

use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.

the class DataLikelihoodTester2 method createSpecifiedTree.

private static TreeModel createSpecifiedTree(String t) throws Exception {
    NewickImporter importer = new NewickImporter(t);
    Tree tree = importer.importTree(null);
    //treeModel
    return new TreeModel(tree);
}