Examples with IntervalList dr.evolution.coalescent.IntervalList used on opensource projects

Example 11 with IntervalList

use of dr.evolution.coalescent.IntervalList in project beast-mcmc by beast-dev.

the class BigFastTreeTreeIntervalsTest method testCompareIntervals.

public void testCompareIntervals() throws TreeUtils.MissingTaxonException, IOException, Importer.ImportException {
    NewickImporter importer = new NewickImporter("(Lishui/LS557/2020:0,((Netherlands/Utrecht_10015/2020:0.00006795400000000001,USA/IL-NM073/2020:0.00006799599999999999):0.000033976,England/LOND-D604F/2020:0.000101963):0.000033968,Guangdong/2020XN4459-P0041/2020:0.000000005,(Portugal/PT0063/2020:0,(Spain/Zaragoza2486/2020:0.000102605,Scotland/CVR746/2020:0.000000005,Spain/COV000882/2020:0.000067956,Colombia/INS-79253/2020:0.000101944,Uruguay/UY-4/2020:0.000031515):0.000033979,(Spain/CastillaLaMancha201329/2020:0.000000005,Netherlands/NoordHolland_10011/2020:0.000033987):0.00006799,England/LIVE-9CE87/2020:0.00013727299999999998,Spain/Granada-COV002916/2020:0.000033979999999999997):0.000033968,((USA/VI-CDC-3705/2020:0.000000005,Australia/VIC229/2020:0,USA/MA-MGH-00063/2020:0,(USA/WA-S41/2020:0.000068895,USA/WA-UW114/2020:0.000067978,USA/WA-UW17/2020:0.000000005,(USA/WA-S582/2020:0,USA/WA-UW-1682/2020:0.000000005,USA/WA-S994/2020:0.000101934):0.000033955,USA/WA-S121/2020:0.000000005,USA/WA-S154/2020:0.000067982,USA/WA-UW37/2020:0,USA/WA-S321/2020:0,USA/WA-S445/2020:0,USA/WA-S512/2020:0,USA/WA-S33/2020:0.000033979,Canada/BC_6981299/2020:0.000033972,USA/WA-UW-1294/2020:0.000033972,USA/WA-UW-2247/2020:0.000033988,Australia/VIC140/2020:0.000033984,USA/WA-UW61/2020:0.000033972,Canada/BC_8606204/2020:0.000166157,(USA/WA-S734/2020:0,USA/WA-S844/2020:0.000033983):0.000067965,(USA/WA-S1191/2020:0.000067947,USA/WA-S951/2020:0.000101914):0.000095803,Australia/NSW99/2020:0.000101953,(USA/WA-S317/2020:0.000000005,USA/WA-S721/2020:0.00003397):0.00010195700000000001,USA/WA-UW139/2020:0.000135916,USA/WA-S572/2020:0.000033979999999999997,USA/WA-S279/2020:0.000033972,USA/WA-UW28/2020:0.000034002,USA/WA-S114/2020:0.000033969,(USA/WA-S852/2020:0.000203899,(USA/WA-S568/2020:0,USA/WA-S791/2020:0.00006794599999999999):0.000033983):0.000101964,USA/WA-S842/2020:0.000067951):0.000033986,Singapore/302/2020:0.000101947):0.00016677,(((USA/IL-NM0112/2020:0.00003397,USA/IL-NM053/2020:0.000034229,USA/IL-NM059/2020:0.000101967):0.00003397,USA/WI-UW-218/2020:0.000033995):0.000030539,(USA/UT-QDX-63/2020:0,USA/CA-QDX-111/2020:0,USA/TX-HMH0427/2020:0.000203861):0.000101955):0.00023787300000000002):0.000033959,(((Scotland/CVR3203/2020:0.000000005,Scotland/CVR2246/2020:0.000000005,Scotland/GCVR-1714B2/2020:0.000033975999999999995,Scotland/CVR3514/2020:0.000068628):0.000067954,Australia/NT08/2020:0.000034000999999999995):0.00003397,Spain/COV001440/2020:0,Spain/Alcaniz2449/2020:0.000068985,Spain/COV001548/2020:0,USA/WI-WSLH-200057/2020:0.000000005,Spain/Valencia6/2020:0.0000343,Spain/Granada-COV002944/2020:0.000000005,Spain/COV001929/2020:0.000000005,Spain/COV002049/2020:0.000000005,(Spain/Valencia59/2020:0,Spain/Valencia306/2020:0.000000005):0.000033996,(Spain/COV001117/2020:0.00010265,Spain/COV002055/2020:0.000000005,England/20126000104/2020:0.00006758400000000001):0.000033997,Spain/COV001576/2020:0.000000005,Chile/Santiago-1/2020:0.000000005,Spain/COV000721/2020:0.000000005,(Spain/COV001575/2020:0,Spain/COV001505/2020:0):0.000067968,Spain/Madrid_H12_28/2020:0.000067957,Spain/COV001568/2020:0.000033975,England/CAMB-83357/2020:0.000068619,(Spain/Almeria-COV002842/2020:0.000000005,Spain/Malaga-COV002841/2020:0.000000005):0.000067851):0.000169854,Spain/Madrid_LP16_6193/2020:0.00006795299999999999,Singapore/51/2020:0.000044697,(Thailand/Nonthaburi_193/2020:0,Thailand/Bangkok_237/2020:0,Thailand/Bangkok_238/2020:0,((Thailand/Bangkok-0034/2020:0.000000005,Thailand/Bangkok_2295/2020:0,Thailand/Bangkok-0065/2020:0.000047826,Thailand/Bangkok-CONI-0147/2020:0.000033997):0.000033983,Thailand/SI202769-NT/2020:0.000203899):0.000101951):0.00006797500000000001,Shenzhen/SZTH-002/2020:0.000033999);");
    // NewickImporter importer = new NewickImporter("(((0:0.5,(1:1.0,2:1.0)n6:1.0)n7:1.0,3:1.5)n8:1.0,(4:2.0,5:1.51)n9:1.5)n10;");
    // NewickImporter constraintsImporter = new NewickImporter("(((0:0.5,(1:1.0,2:1.0)n6:1.0)n7:1.0,3:1.5)n8:1.0,(4:2.0,5:1.51)n9:1.5)n10;");
    tree = new DefaultTreeModel(importer.importTree(null));
    IntervalList intervals = new TreeIntervals(tree, null, null);
    BigFastTreeIntervals bigFastTreeIntervals = new BigFastTreeIntervals(tree);
    SubtreeLeapOperator op = new SubtreeLeapOperator(tree, 1, 0.0001, SubtreeLeapOperator.DistanceKernelType.NORMAL, AdaptationMode.ADAPTATION_OFF, 0.2);
    NodeHeightOperator nh = new NodeHeightOperator(tree, 1, 1, NodeHeightOperator.OperatorType.UNIFORM, AdaptationMode.ADAPTATION_OFF, 0.25);
    NodeHeightOperator root = new NodeHeightOperator(tree, 1, 0.75, NodeHeightOperator.OperatorType.SCALEROOT, AdaptationMode.ADAPTATION_OFF, 0.25);
    boolean pass = true;
    MathUtils.setSeed(2);
    for (int i = 0; i < 100000; i++) {
        op.doOperation();
        intervals.calculateIntervals();
        // bigFastIntervals.makeDirty();
        bigFastTreeIntervals.calculateIntervals();
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getInterval(j) != bigFastTreeIntervals.getInterval(j)) {
                System.out.println(i);
                System.out.println("interval wrong");
                pass = false;
                break;
            }
        }
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getLineageCount(j) != bigFastTreeIntervals.getLineageCount(j)) {
                System.out.println(i);
                System.out.println("lineage Counts wrong: " + j);
                System.out.println("expected: " + intervals.getLineageCount(j));
                System.out.println("got " + bigFastTreeIntervals.getLineageCount(j));
                pass = false;
                break;
            }
        }
        for (int j = 0; j < bigFastTreeIntervals.getIntervalCount(); j++) {
            if (intervals.getIntervalTime(j) != bigFastTreeIntervals.getIntervalTime(j)) {
                System.out.println(i);
                System.out.println("times wrong");
                pass = false;
                break;
            }
        }
        if (!pass) {
            break;
        }
    }
    assertTrue(pass);
}

Example 12 with IntervalList

use of dr.evolution.coalescent.IntervalList in project beast-mcmc by beast-dev.

the class GMRFSkyrideLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    XMLObject cxo = xo.getChild(POPULATION_PARAMETER);
    Parameter popParameter = (Parameter) cxo.getChild(Parameter.class);
    cxo = xo.getChild(PRECISION_PARAMETER);
    Parameter precParameter = (Parameter) cxo.getChild(Parameter.class);
    boolean buildIntervalNodeMapping = xo.getAttribute(BUILD_MAPPING, false);
    List<IntervalList> intervalsList = new ArrayList<IntervalList>();
    List<Tree> treeList = new ArrayList<Tree>();
    if (xo.getChild(POPULATION_TREE) != null) {
        cxo = xo.getChild(POPULATION_TREE);
        for (int i = 0; i < cxo.getChildCount(); i++) {
            Object testObject = cxo.getChild(i);
            if (testObject instanceof Tree) {
                treeList.add((TreeModel) testObject);
            // TreeIntervals treeIntervals;
            // try {
            // treeIntervals = new TreeIntervals((Tree) testObject, null, null);
            // } catch (TreeUtils.MissingTaxonException mte) {
            // throw new XMLParseException("Taxon, " + mte + ", in " + getParserName() + " was not found in the tree.");
            // }
            // intervalsList.add(treeIntervals);
            }
        }
    }
    if (xo.getChild(INTERVALS) != null) {
        cxo = xo.getChild(INTERVALS);
        intervalsList = new ArrayList<IntervalList>();
        for (int i = 0; i < cxo.getChildCount(); i++) {
            Object testObject = cxo.getChild(i);
            if (testObject instanceof IntervalList) {
                intervalsList.add((IntervalList) testObject);
            }
        }
    }
    cxo = xo.getChild(GROUP_SIZES);
    Parameter groupParameter = null;
    if (cxo != null) {
        groupParameter = (Parameter) cxo.getChild(Parameter.class);
        if (popParameter.getDimension() != groupParameter.getDimension())
            throw new XMLParseException("Population and group size parameters must have the same length");
    }
    Parameter lambda;
    if (xo.getChild(LAMBDA_PARAMETER) != null) {
        cxo = xo.getChild(LAMBDA_PARAMETER);
        lambda = (Parameter) cxo.getChild(Parameter.class);
    } else {
        lambda = new Parameter.Default(LAMBDA_PARAMETER, 1.0);
    }
    Parameter gridPoints = null;
    if (xo.getChild(GRID_POINTS) != null) {
        cxo = xo.getChild(GRID_POINTS);
        gridPoints = (Parameter) cxo.getChild(Parameter.class);
    }
    Parameter numGridPoints = null;
    if (xo.getChild(NUM_GRID_POINTS) != null) {
        cxo = xo.getChild(NUM_GRID_POINTS);
        numGridPoints = (Parameter) cxo.getChild(Parameter.class);
    }
    Parameter cutOff = null;
    if (xo.getChild(CUT_OFF) != null) {
        cxo = xo.getChild(CUT_OFF);
        cutOff = (Parameter) cxo.getChild(Parameter.class);
    }
    Parameter phi = null;
    if (xo.getChild(PHI_PARAMETER) != null) {
        cxo = xo.getChild(PHI_PARAMETER);
        phi = (Parameter) cxo.getChild(Parameter.class);
    }
    List<Parameter> firstObservedIndex = null;
    if (xo.hasChildNamed(FIRST_OBSERVED_INDEX)) {
        firstObservedIndex = new ArrayList<Parameter>();
        cxo = xo.getChild(FIRST_OBSERVED_INDEX);
        final int numInd = cxo.getChildCount();
        for (int i = 0; i < numInd; ++i) {
            firstObservedIndex.add((Parameter) cxo.getChild(i));
        }
    }
    List<Parameter> lastObservedIndex = null;
    if (xo.hasChildNamed(LAST_OBSERVED_INDEX)) {
        lastObservedIndex = new ArrayList<Parameter>();
        cxo = xo.getChild(LAST_OBSERVED_INDEX);
        final int numObsInd = cxo.getChildCount();
        for (int i = 0; i < numObsInd; ++i) {
            lastObservedIndex.add((Parameter) cxo.getChild(i));
        }
    }
    Parameter ploidyFactors = null;
    if (xo.getChild(PLOIDY) != null) {
        cxo = xo.getChild(PLOIDY);
        ploidyFactors = (Parameter) cxo.getChild(Parameter.class);
    } else {
        if (intervalsList.size() != 0) {
            ploidyFactors = new Parameter.Default(PLOIDY, intervalsList.size());
            for (int i = 0; i < intervalsList.size(); i++) {
                ploidyFactors.setParameterValue(i, 1.0);
            }
        } else {
            ploidyFactors = new Parameter.Default(PLOIDY, treeList.size());
            for (int i = 0; i < treeList.size(); i++) {
                ploidyFactors.setParameterValue(i, 1.0);
            }
        }
    }
    Parameter betaParameter = null;
    if (xo.hasChildNamed(SINGLE_BETA)) {
        betaParameter = (Parameter) xo.getElementFirstChild(SINGLE_BETA);
    }
    List<Parameter> betaList = null;
    if (xo.getChild(BETA_PARAMETER) != null) {
        betaList = new ArrayList<Parameter>();
        cxo = xo.getChild(BETA_PARAMETER);
        final int numBeta = cxo.getChildCount();
        for (int i = 0; i < numBeta; ++i) {
            betaList.add((Parameter) cxo.getChild(i));
        }
    }
    List<Parameter> deltaList = new ArrayList<Parameter>();
    if (xo.getChild(DELTA_PARAMETER) != null) {
        cxo = xo.getChild(DELTA_PARAMETER);
        final int numDelta = cxo.getChildCount();
        if (numDelta != betaList.size()) {
            throw new XMLParseException("Cannot have different number of delta and beta parameters");
        }
        for (int i = 0; i < numDelta; ++i) {
            deltaList.add((Parameter) cxo.getChild(i));
        }
    } else {
        deltaList = null;
    }
    MatrixParameter dMatrix = null;
    if (xo.getChild(COVARIATE_MATRIX) != null) {
        cxo = xo.getChild(COVARIATE_MATRIX);
        dMatrix = (MatrixParameter) cxo.getChild(MatrixParameter.class);
    }
    boolean timeAwareSmoothing = GMRFSkyrideLikelihood.TIME_AWARE_IS_ON_BY_DEFAULT;
    if (xo.hasAttribute(TIME_AWARE_SMOOTHING)) {
        timeAwareSmoothing = xo.getBooleanAttribute(TIME_AWARE_SMOOTHING);
    }
    if (dMatrix != null) {
        if (dMatrix.getRowDimension() != popParameter.getDimension())
            throw new XMLParseException("Design matrix row dimension must equal the population parameter length.");
        if (dMatrix.getColumnDimension() != betaParameter.getDimension())
            throw new XMLParseException("Design matrix column dimension must equal the regression coefficient length.");
    }
    List<Parameter> covPrecParamRecent = null;
    List<Parameter> covPrecParamDistant = null;
    if (xo.hasChildNamed(COV_PREC_REC)) {
        covPrecParamRecent = new ArrayList<Parameter>();
        cxo = xo.getChild(COV_PREC_REC);
        for (int i = 0; i < cxo.getChildCount(); ++i) {
            covPrecParamRecent.add((Parameter) cxo.getChild(i));
        }
    }
    if (xo.hasChildNamed(COV_PREC_DIST)) {
        covPrecParamDistant = new ArrayList<Parameter>();
        cxo = xo.getChild(COV_PREC_DIST);
        for (int i = 0; i < cxo.getChildCount(); ++i) {
            covPrecParamDistant.add((Parameter) cxo.getChild(i));
        }
    }
    if (xo.hasChildNamed(COV_PREC_PARAM)) {
        if (firstObservedIndex != null) {
            covPrecParamRecent = new ArrayList<Parameter>();
        }
        if (lastObservedIndex != null) {
            covPrecParamDistant = new ArrayList<Parameter>();
        }
        cxo = xo.getChild(COV_PREC_PARAM);
        for (int i = 0; i < cxo.getChildCount(); ++i) {
            if (firstObservedIndex != null) {
                covPrecParamRecent.add((Parameter) cxo.getChild(i));
            }
            if (lastObservedIndex != null) {
                covPrecParamDistant.add((Parameter) cxo.getChild(i));
            }
        }
    }
    if ((covPrecParamDistant == null && lastObservedIndex != null) || (covPrecParamDistant != null && lastObservedIndex == null)) {
        throw new XMLParseException("Must specify both lastObservedIndex and covariatePrecision");
    }
    if ((covPrecParamRecent == null && firstObservedIndex != null) || (covPrecParamRecent != null && firstObservedIndex == null)) {
        throw new XMLParseException("Must specify both firstObservedIndex and covariatePrecision");
    }
    Parameter recentIndices = null;
    if (xo.getChild(REC_INDICES) != null) {
        cxo = xo.getChild(REC_INDICES);
        recentIndices = (Parameter) cxo.getChild(Parameter.class);
    }
    if (firstObservedIndex == null && recentIndices != null) {
        throw new XMLParseException("Cannot specify covIndicesMissingRecent without specifying firstObservedIndex");
    }
    Parameter distantIndices = null;
    if (xo.getChild(DIST_INDICES) != null) {
        cxo = xo.getChild(DIST_INDICES);
        distantIndices = (Parameter) cxo.getChild(Parameter.class);
    }
    if (lastObservedIndex == null && distantIndices != null) {
        throw new XMLParseException("Cannot specify covIndicesMissingDistant without specifying lastObservedIndex");
    }
    List<MatrixParameter> covariates = null;
    if (xo.hasChildNamed(COVARIATES)) {
        covariates = new ArrayList<MatrixParameter>();
        cxo = xo.getChild(COVARIATES);
        final int numCov = cxo.getChildCount();
        for (int i = 0; i < numCov; ++i) {
            covariates.add((MatrixParameter) cxo.getChild(i));
        }
    }
    if ((covariates != null && betaList == null) || (covariates == null && betaList != null))
        throw new XMLParseException("Must specify both a set of regression coefficients and a design matrix.");
    boolean useGlmModel = xo.getAttribute(USE_GLM_MODEL, false);
    if (useGlmModel) {
        GeneralizedLinearModel glm = (GeneralizedLinearModel) xo.getChild(GeneralizedLinearModel.class);
        covariates = new ArrayList<MatrixParameter>();
        betaList = new ArrayList<Parameter>();
        List<DesignMatrix> designMat = glm.getDesignMatrix();
        List<Parameter> indepParam = glm.getIndependentParameter();
        List<Parameter> indepParamDelta = glm.getIndependentParameterDelta();
        deltaList = new ArrayList<Parameter>();
        for (int i = 0; i < indepParam.get(0).getSize(); i++) {
            MatrixParameter matParam = new MatrixParameter("covariate values", 1, designMat.get(0).getRowDimension());
            for (int j = 0; j < matParam.getRowDimension(); j++) {
                matParam.setParameterValue(0, j, designMat.get(0).getParameterValue(0, j));
            }
            covariates.add(matParam);
            Parameter betaParam = new Parameter.Default(1);
            betaParam.setParameterValue(0, indepParam.get(0).getParameterValue(i));
            betaList.add(betaParam);
            if (indepParamDelta != null) {
                Parameter deltaParam = new Parameter.Default(1);
                deltaParam.setParameterValue(0, indepParamDelta.get(0).getParameterValue(i));
                deltaList.add(deltaParam);
            }
        }
    }
    /*
        if (xo.getAttribute(RANDOMIZE_TREE, false)) {
            for (Tree tree : treeList) {
                if (tree instanceof TreeModel) {
                    GMRFSkyrideLikelihood.checkTree((TreeModel) tree);
                } else {
                    throw new XMLParseException("Can not randomize a fixed tree");
                }
            }
        }*/
    boolean rescaleByRootHeight = xo.getAttribute(RESCALE_BY_ROOT_ISSUE, true);
    Logger.getLogger("dr.evomodel").info("The " + SKYLINE_LIKELIHOOD + " has " + (timeAwareSmoothing ? "time aware smoothing" : "uniform smoothing"));
    if (xo.getAttribute(OLD_SKYRIDE, true) && xo.getName().compareTo(SKYGRID_LIKELIHOOD) != 0) {
        return new OldGMRFSkyrideLikelihood(treeList, popParameter, groupParameter, precParameter, lambda, betaParameter, dMatrix, timeAwareSmoothing, rescaleByRootHeight, buildIntervalNodeMapping);
    } else {
        if (intervalsList.size() > 0) {
            if (xo.getChild(GRID_POINTS) != null) {
                return new GMRFSkygridLikelihood(intervalsList, popParameter, groupParameter, precParameter, lambda, betaParameter, dMatrix, timeAwareSmoothing, gridPoints, covariates, ploidyFactors, firstObservedIndex, lastObservedIndex, covPrecParamRecent, covPrecParamDistant, recentIndices, distantIndices, betaList);
            } else {
                return new GMRFSkygridLikelihood(intervalsList, popParameter, groupParameter, precParameter, lambda, betaParameter, dMatrix, timeAwareSmoothing, cutOff.getParameterValue(0), (int) numGridPoints.getParameterValue(0), phi, ploidyFactors);
            }
        } else {
            if (xo.getChild(GRID_POINTS) != null) {
                return new GMRFMultilocusSkyrideLikelihood(treeList, popParameter, groupParameter, precParameter, lambda, betaParameter, dMatrix, timeAwareSmoothing, gridPoints, covariates, ploidyFactors, firstObservedIndex, lastObservedIndex, covPrecParamRecent, covPrecParamDistant, recentIndices, distantIndices, betaList, deltaList);
            } else {
                return new GMRFMultilocusSkyrideLikelihood(treeList, popParameter, groupParameter, precParameter, lambda, betaParameter, dMatrix, timeAwareSmoothing, cutOff.getParameterValue(0), (int) numGridPoints.getParameterValue(0), phi, ploidyFactors);
            }
        }
    }
}

Example 13 with IntervalList

use of dr.evolution.coalescent.IntervalList in project beast-mcmc by beast-dev.

the class SmoothSkygridLikelihoodParser method parseXMLObject.

@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    List<IntervalList> intervalList = new ArrayList<>();
    List<TreeIntervals> debugIntervalList = new ArrayList<>();
    if (xo.hasChildNamed(INTERVALS)) {
        XMLObject cxo = xo.getChild(INTERVALS);
        for (int i = 0; i < cxo.getChildCount(); ++i) {
            intervalList.add((IntervalList) cxo.getChild(i));
        }
    } else {
        XMLObject cxo = xo.getChild(POPULATION_TREE);
        for (int i = 0; i < cxo.getChildCount(); ++i) {
            TreeModel tree = (TreeModel) cxo.getChild(i);
            intervalList.add(new BigFastTreeIntervals(tree));
            debugIntervalList.add(new TreeIntervals(tree));
        }
    }
    Parameter logPopSizes = (Parameter) xo.getElementFirstChild(POPULATION_PARAMETER);
    Parameter gridPoints = (Parameter) xo.getElementFirstChild(GRID_POINTS);
    if (!SmoothSkygridLikelihood.checkValidParameters(logPopSizes, gridPoints)) {
        throw new XMLParseException("Invalid initial parameters");
    }
    SmoothSkygridLikelihood likelihood = new SmoothSkygridLikelihood(xo.getId(), intervalList, logPopSizes, gridPoints);
    likelihood.setDebugIntervalList(debugIntervalList);
    return likelihood;
}