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

Example 6 with Taxa

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

the class PLCoalescentSimulator method main.

public static void main(String[] arg) throws IOException {
    // READ DEMOGRAPHIC FUNCTION
    String filename = arg[0];
    BufferedReader reader = new BufferedReader(new FileReader(filename));
    double popSizeScale = 1.0;
    double generationTime = 1.0;
    if (arg.length > 2) {
        popSizeScale = Double.parseDouble(arg[2]);
    }
    if (arg.length > 3) {
        generationTime = Double.parseDouble(arg[3]);
    }
    PrintWriter populationFuncLogger = null;
    if (arg.length > 5) {
        String logFileName = arg[5];
        if (logFileName.equals("-")) {
            populationFuncLogger = new PrintWriter(System.out);
        } else {
            populationFuncLogger = new PrintWriter(new FileWriter(logFileName));
        }
    }
    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();
    // READ SAMPLE TIMES
    String samplesFilename = arg[1];
    reader = new BufferedReader(new FileReader(samplesFilename));
    line = reader.readLine();
    Taxa taxa = new Taxa();
    int id = 0;
    while (line != null) {
        if (!line.startsWith("#")) {
            tokens = line.split("[\t ]+");
            if (tokens.length == 4) {
                double t0 = Double.parseDouble(tokens[0]);
                double t1 = Double.parseDouble(tokens[1]);
                double dt = Double.parseDouble(tokens[2]);
                int k = Integer.parseInt(tokens[3]);
                for (double time = t0; time <= t1; time += dt) {
                    double sampleTime = time / generationTime;
                    for (int i = 0; i < k; i++) {
                        Taxon taxon = new Taxon("t" + id);
                        taxon.setAttribute(dr.evolution.util.Date.DATE, new Date(sampleTime, Units.Type.GENERATIONS, true));
                        taxa.addTaxon(taxon);
                        id += 1;
                    }
                }
            } else {
                // sample times are in the same units as simulation
                double sampleTime = Double.parseDouble(tokens[0]) / generationTime;
                int count = Integer.parseInt(tokens[1]);
                for (int i = 0; i < count; i++) {
                    Taxon taxon = new Taxon(id + "");
                    taxon.setAttribute(dr.evolution.util.Date.DATE, new Date(sampleTime, Units.Type.GENERATIONS, true));
                    taxa.addTaxon(taxon);
                    id += 1;
                }
            }
        }
        line = reader.readLine();
    }
    double minTheta = Double.MAX_VALUE;
    double maxTheta = 0.0;
    PrintWriter out;
    if (arg.length > 4) {
        out = new PrintWriter(new FileWriter(arg[4]));
    } else {
        out = new PrintWriter(System.out);
    }
    int pp = 0;
    for (List<Double> popSize : popSizes) {
        double[] thetas = new double[popSize.size()];
        double[] intervals = new double[times.size() - 1];
        if (populationFuncLogger != null) {
            populationFuncLogger.println("# " + pp);
            ++pp;
        }
        // must reverse the direction of the model
        for (int j = intervals.length; j > 0; j--) {
            intervals[intervals.length - j] = times.get(j) - times.get(j - 1);
            final double theta = popSize.get(j) * popSizeScale;
            thetas[intervals.length - j] = theta;
            if (theta < minTheta) {
                minTheta = theta;
            }
            if (theta > maxTheta) {
                maxTheta = theta;
            }
            final double t = times.get(intervals.length) - times.get(j);
            if (populationFuncLogger != null) {
                populationFuncLogger.println(t + "\t" + theta);
            }
        }
        if (debug != null) {
            debug.println("min theta = " + minTheta);
            debug.println("max theta = " + maxTheta);
        }
        PiecewiseLinearPopulation demo = new PiecewiseLinearPopulation(intervals, thetas, Units.Type.GENERATIONS);
        CoalescentSimulator simulator = new CoalescentSimulator();
        Tree tree = simulator.simulateTree(taxa, demo);
        out.println(TreeUtils.newick(tree));
        if (debug != null) {
            debug.println(TreeUtils.newick(tree));
        }
    }
    if (populationFuncLogger != null) {
        populationFuncLogger.flush();
        populationFuncLogger.close();
    }
    out.flush();
    out.close();
}

Example 7 with Taxa

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

the class PlinkImporter method addTaxonAttribute.

public void addTaxonAttribute(TaxonList inTaxonList, String traitName) {
    taxonList = new Taxa();
    for (Taxon taxon : inTaxonList) {
        List<Double> valueList = taxonHash.get(taxon.getId());
        if (valueList == null) {
            Logger.getLogger("dr.evolution").warning("Taxon " + taxon.getId() + " not found in PLINK data");
        } else {
            String string = makeStringAttribute(valueList);
            ((Taxa) taxonList).addTaxon(taxon);
            taxon.setAttribute(traitName, string);
        }
        if (DEBUG) {
            System.err.println("Added trait for " + taxon.getId());
        }
    }
}

Example 8 with Taxa

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

the class MultiTreeIntervalsParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    XMLObject cxo = xo.getChild(TREES);
    List<Tree> trees = new ArrayList<Tree>(cxo.getAllChildren(Tree.class));
    Taxa singletonTaxa = null;
    if (xo.hasChildNamed(SINGLETONS)) {
        singletonTaxa = (Taxa) xo.getElementFirstChild(SINGLETONS);
    }
    boolean includeStems = xo.getBooleanAttribute(INCLUDE_STEMS);
    double cutoffTime = 0.0;
    if (includeStems) {
        if (!xo.hasAttribute(CUTOFF)) {
            throw new XMLParseException("MultiTreeIntervals needs a cutoff time if it is to include stems");
        }
        cutoffTime = xo.getDoubleAttribute(CUTOFF);
    }
    return new MultiTreeIntervals(trees, singletonTaxa, includeStems, cutoffTime);
}

Example 9 with Taxa

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

the class SpeciationLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    XMLObject cxo = xo.getChild(MODEL);
    final SpeciationModel specModel = (SpeciationModel) cxo.getChild(SpeciationModel.class);
    cxo = xo.getChild(TREE);
    final Tree tree = (Tree) cxo.getChild(Tree.class);
    Set<Taxon> excludeTaxa = null;
    if (xo.hasChildNamed(INCLUDE)) {
        excludeTaxa = new HashSet<Taxon>();
        for (int i = 0; i < tree.getTaxonCount(); i++) {
            excludeTaxa.add(tree.getTaxon(i));
        }
        cxo = xo.getChild(INCLUDE);
        for (int i = 0; i < cxo.getChildCount(); i++) {
            TaxonList taxonList = (TaxonList) cxo.getChild(i);
            for (int j = 0; j < taxonList.getTaxonCount(); j++) {
                excludeTaxa.remove(taxonList.getTaxon(j));
            }
        }
    }
    if (xo.hasChildNamed(EXCLUDE)) {
        excludeTaxa = new HashSet<Taxon>();
        cxo = xo.getChild(EXCLUDE);
        for (int i = 0; i < cxo.getChildCount(); i++) {
            TaxonList taxonList = (TaxonList) cxo.getChild(i);
            for (int j = 0; j < taxonList.getTaxonCount(); j++) {
                excludeTaxa.add(taxonList.getTaxon(j));
            }
        }
    }
    if (excludeTaxa != null) {
        Logger.getLogger("dr.evomodel").info("Speciation model excluding " + excludeTaxa.size() + " taxa from prior - " + (tree.getTaxonCount() - excludeTaxa.size()) + " taxa remaining.");
    }
    final XMLObject cal = xo.getChild(CALIBRATION);
    if (cal != null) {
        if (excludeTaxa != null) {
            throw new XMLParseException("Sorry, not implemented: internal calibration prior + excluded taxa");
        }
        List<Distribution> dists = new ArrayList<Distribution>();
        List<Taxa> taxa = new ArrayList<Taxa>();
        List<Boolean> forParent = new ArrayList<Boolean>();
        // (Statistic) cal.getChild(Statistic.class);
        Statistic userPDF = null;
        for (int k = 0; k < cal.getChildCount(); ++k) {
            final Object ck = cal.getChild(k);
            if (DistributionLikelihood.class.isInstance(ck)) {
                dists.add(((DistributionLikelihood) ck).getDistribution());
            } else if (Distribution.class.isInstance(ck)) {
                dists.add((Distribution) ck);
            } else if (Taxa.class.isInstance(ck)) {
                final Taxa tx = (Taxa) ck;
                taxa.add(tx);
                forParent.add(tx.getTaxonCount() == 1);
            } else if (Statistic.class.isInstance(ck)) {
                if (userPDF != null) {
                    throw new XMLParseException("more than one userPDF correction???");
                }
                userPDF = (Statistic) cal.getChild(Statistic.class);
            } else {
                XMLObject cko = (XMLObject) ck;
                assert cko.getChildCount() == 2;
                for (int i = 0; i < 2; ++i) {
                    final Object chi = cko.getChild(i);
                    if (DistributionLikelihood.class.isInstance(chi)) {
                        dists.add(((DistributionLikelihood) chi).getDistribution());
                    } else if (Distribution.class.isInstance(chi)) {
                        dists.add((Distribution) chi);
                    } else if (Taxa.class.isInstance(chi)) {
                        taxa.add((Taxa) chi);
                        boolean fp = ((Taxa) chi).getTaxonCount() == 1;
                        if (cko.hasAttribute(PARENT)) {
                            boolean ufp = cko.getBooleanAttribute(PARENT);
                            if (fp && !ufp) {
                                throw new XMLParseException("forParent==false for a single taxon?? (must be true)");
                            }
                            fp = ufp;
                        }
                        forParent.add(fp);
                    } else {
                        assert false;
                    }
                }
            }
        }
        if (dists.size() != taxa.size()) {
            throw new XMLParseException("Mismatch in number of distributions and taxa specs");
        }
        try {
            final String correction = cal.getAttribute(CORRECTION, EXACT);
            final CalibrationPoints.CorrectionType type = correction.equals(EXACT) ? CalibrationPoints.CorrectionType.EXACT : (correction.equals(APPROX) ? CalibrationPoints.CorrectionType.APPROXIMATED : (correction.equals(NONE) ? CalibrationPoints.CorrectionType.NONE : (correction.equals(PEXACT) ? CalibrationPoints.CorrectionType.PEXACT : null)));
            if (cal.hasAttribute(CORRECTION) && type == null) {
                throw new XMLParseException("correction type == " + correction + "???");
            }
            final CalibrationPoints calib = new CalibrationPoints(tree, specModel.isYule(), dists, taxa, forParent, userPDF, type);
            final SpeciationLikelihood speciationLikelihood = new SpeciationLikelihood(tree, specModel, null, calib);
            return speciationLikelihood;
        } catch (IllegalArgumentException e) {
            throw new XMLParseException(e.getMessage());
        }
    }
    return new SpeciationLikelihood(tree, specModel, excludeTaxa, null);
}

Example 10 with Taxa

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

the class MicrosatellitePatternParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Taxa taxonList = (Taxa) xo.getChild(Taxa.class);
    Microsatellite microsatellite = (Microsatellite) xo.getChild(Microsatellite.class);
    String[] strLengths = ((String) xo.getElementFirstChild(MICROSAT_SEQ)).split(",");
    int[] pattern = new int[strLengths.length];
    try {
        for (int i = 0; i < strLengths.length; i++) {
            pattern[i] = microsatellite.getState(strLengths[i]);
        }
    } catch (NumberFormatException nfe) {
        throw new XMLParseException("Unable to parse microsatellite data: " + nfe.getMessage());
    } catch (IllegalArgumentException iae) {
        throw new XMLParseException("Unable to parse microsatellite data: " + iae.getMessage());
    }
    Patterns microsatPat = new Patterns(microsatellite, taxonList);
    microsatPat.addPattern(pattern);
    microsatPat.setId((String) xo.getAttribute(ID));
    if (xo.getAttribute(PRINT_DETAILS, true)) {
        printDetails(microsatPat);
    }
    if (xo.getAttribute(PRINT_MSAT_PATTERN_CONTENT, true)) {
        printMicrosatContent(microsatPat);
    }
    return microsatPat;
}