Examples with DataType dr.evolution.datatype.DataType used on opensource projects

Example 61 with DataType

use of dr.evolution.datatype.DataType in project beast-mcmc by beast-dev.

the class IstvanOperator method doOperation.

public double doOperation() {
    Tree tree = likelihood.getTreeModel();
    alignment = likelihood.getAlignment();
    // System.out.println("Incoming alignment");
    // System.out.println(alignment);
    // System.out.println();
    SubstitutionModel substModel = likelihood.getSiteModel().getSubstitutionModel();
    // initialize the iParent and iTau arrays based on the given tree.
    initTree(tree, likelihood.getSiteModel().getMu());
    int[] treeIndex = new int[tree.getTaxonCount()];
    for (int i = 0; i < treeIndex.length; i++) {
        treeIndex[i] = tree.getTaxonIndex(alignment.getTaxonId(i));
    }
    // initialize the iAlignment array from the given alignment.
    initAlignment(alignment, treeIndex);
    // initialize the iProbs array from the substitution model -- must be called after populating tree!
    initSubstitutionModel(substModel);
    DataType dataType = substModel.getDataType();
    proposal.setGapSymbol(dataType.getGapState());
    int[][] returnedAlignment = new int[iAlignment.length][];
    // System.out.println("Initialization done, starting proposal proper...");
    double logq = proposal.propose(iAlignment, iProbs, iBaseFreqs, iParent, iTau, returnedAlignment, tuning, exponent, gapPenalty);
    // System.out.println("Proposal finished, logq=" + logq);
    // create new alignment object
    SimpleAlignment newAlignment = new SimpleAlignment();
    for (int i = 0; i < alignment.getTaxonCount(); i++) {
        StringBuffer seqBuffer = new StringBuffer();
        for (int j = 0; j < returnedAlignment[i].length; j++) {
            seqBuffer.append(dataType.getChar(returnedAlignment[treeIndex[i]][j]));
        }
        // add sequences in order of tree
        String seqString = seqBuffer.toString();
        Sequence sequence = new Sequence(alignment.getTaxon(i), seqString);
        newAlignment.addSequence(sequence);
        String oldunaligned = alignment.getUnalignedSequenceString(i);
        String unaligned = newAlignment.getUnalignedSequenceString(i);
        if (!unaligned.equals(oldunaligned)) {
            System.err.println("Sequence changed from:");
            System.err.println("old:'" + oldunaligned + "'");
            System.err.println("new:'" + unaligned + "'");
            throw new RuntimeException();
        }
    }
    // System.out.println("Outgoing alignment");
    // System.out.println(newAlignment);
    // System.out.println();
    likelihood.setAlignment(newAlignment);
    return logq;
}

Example 62 with DataType

use of dr.evolution.datatype.DataType in project beast-mcmc by beast-dev.

the class ComplexSubstitutionModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    DataType dataType = DataTypeUtils.getDataType(xo);
    if (dataType == null)
        dataType = (DataType) xo.getChild(DataType.class);
    XMLObject cxo = xo.getChild(RATES);
    Parameter ratesParameter = (Parameter) cxo.getChild(Parameter.class);
    int rateCount = (dataType.getStateCount() - 1) * dataType.getStateCount();
    if (ratesParameter.getDimension() != rateCount) {
        throw new XMLParseException("Rates parameter in " + getParserName() + " element should have " + (rateCount) + " dimensions.  However parameter dimension is " + ratesParameter.getDimension());
    }
    cxo = xo.getChild(ROOT_FREQUENCIES);
    FrequencyModel rootFreq = (FrequencyModel) cxo.getChild(FrequencyModel.class);
    if (dataType != rootFreq.getDataType()) {
        throw new XMLParseException("Data type of " + getParserName() + " element does not match that of its rootFrequencyModel.");
    }
    Parameter indicators = null;
    if (xo.hasChildNamed(INDICATOR)) {
        indicators = (Parameter) ((XMLObject) xo.getChild(INDICATOR)).getChild(Parameter.class);
        if (ratesParameter.getDimension() != indicators.getDimension())
            throw new XMLParseException("Rate parameter dimension must match indicator parameter dimension");
    }
    StringBuffer sb = new StringBuffer().append("Constructing a complex substitution model using\n").append("\tRate parameters: ").append(ratesParameter.getId()).append("\n").append("\tRoot frequency model: ").append(rootFreq.getId()).append("\n");
    ComplexSubstitutionModel model;
    if (indicators == null)
        model = new ComplexSubstitutionModel(xo.getId(), dataType, rootFreq, ratesParameter);
    else {
        boolean randomize = xo.getAttribute(RANDOMIZE, false);
        boolean connected = xo.getAttribute(CONNECTED, false);
        model = new SVSComplexSubstitutionModel(xo.getId(), dataType, rootFreq, ratesParameter, indicators);
        if (randomize) {
            BayesianStochasticSearchVariableSelection.Utils.randomize(indicators, dataType.getStateCount(), false);
            boolean valid = !Double.isInfinite(model.getLogLikelihood());
            if (!valid) {
                throw new XMLParseException("Poor tolerance in complex substitution model.  Please retry analysis using BEAGLE");
            }
        }
        sb.append("\tBSSVS indicators: ").append(indicators.getId()).append("\n");
        sb.append("\tGraph must be connected: ").append(connected).append("\n");
    }
    boolean doNormalization = xo.getAttribute(NORMALIZATION, true);
    model.setNormalization(doNormalization);
    sb.append("\tNormalized: ").append(doNormalization).append("\n");
    boolean checkConditioning = xo.getAttribute(CHECK_CONDITIONING, true);
    model.setCheckConditioning(checkConditioning);
    if (checkConditioning) {
        double maxConditionNumber = xo.getAttribute(MAX_CONDITION_NUMBER, 1000);
        model.setMaxConditionNumber(maxConditionNumber);
        sb.append("\tMax. condition number: ").append(maxConditionNumber).append("\n");
    }
    int maxIterations = xo.getAttribute(MAX_ITERATIONS, 1000);
    model.setMaxIterations(maxIterations);
    sb.append("\tMax iterations: ").append(maxIterations).append("\n");
    sb.append("\t\tPlease cite Edwards, Suchard et al. (2011)\n");
    Logger.getLogger("dr.evomodel.substmodel").info(sb.toString());
    return model;
}

Example 63 with DataType

use of dr.evolution.datatype.DataType in project beast-mcmc by beast-dev.

the class GLMSubstitutionModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    DataType dataType = DataTypeUtils.getDataType(xo);
    if (dataType == null)
        dataType = (DataType) xo.getChild(DataType.class);
    int rateCount = (dataType.getStateCount() - 1) * dataType.getStateCount();
    LogLinearModel glm = (LogLinearModel) xo.getChild(GeneralizedLinearModel.class);
    int length = glm.getXBeta().length;
    if (length != rateCount) {
        throw new XMLParseException("Rates parameter in " + getParserName() + " element should have " + (rateCount) + " dimensions.  However GLM dimension is " + length);
    }
    XMLObject cxo = xo.getChild(ComplexSubstitutionModelParser.ROOT_FREQUENCIES);
    FrequencyModel rootFreq = (FrequencyModel) cxo.getChild(FrequencyModel.class);
    if (dataType != rootFreq.getDataType()) {
        throw new XMLParseException("Data type of " + getParserName() + " element does not match that of its rootFrequencyModel.");
    }
    return new GLMSubstitutionModel(xo.getId(), dataType, rootFreq, glm);
}

Example 64 with DataType

use of dr.evolution.datatype.DataType in project beast-mcmc by beast-dev.

the class GeneralSubstitutionModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Parameter ratesParameter = null;
    FrequencyModel freqModel = null;
    if (xo.hasChildNamed(FREQUENCIES)) {
        XMLObject cxo = xo.getChild(FREQUENCIES);
        freqModel = (FrequencyModel) cxo.getChild(FrequencyModel.class);
    }
    DataType dataType = DataTypeUtils.getDataType(xo);
    if (dataType == null)
        dataType = (DataType) xo.getChild(DataType.class);
    if (dataType == null)
        dataType = freqModel.getDataType();
    if (dataType != freqModel.getDataType()) {
        throw new XMLParseException("Data type of " + getParserName() + " element does not match that of its frequencyModel.");
    }
    XMLObject cxo = xo.getChild(RATES);
    ratesParameter = (Parameter) cxo.getChild(Parameter.class);
    int states = dataType.getStateCount();
    Logger.getLogger("dr.evomodel").info("  General Substitution Model (stateCount=" + states + ")");
    boolean hasRelativeRates = cxo.hasChildNamed(RELATIVE_TO) || (cxo.hasAttribute(RELATIVE_TO) && cxo.getIntegerAttribute(RELATIVE_TO) > 0);
    int nonReversibleRateCount = ((dataType.getStateCount() - 1) * dataType.getStateCount());
    int reversibleRateCount = (nonReversibleRateCount / 2);
    boolean isNonReversible = ratesParameter.getDimension() == nonReversibleRateCount;
    boolean hasIndicator = xo.hasChildNamed(INDICATOR);
    if (!hasRelativeRates) {
        Parameter indicatorParameter = null;
        if (ratesParameter.getDimension() != reversibleRateCount && ratesParameter.getDimension() != nonReversibleRateCount) {
            throw new XMLParseException("Rates parameter in " + getParserName() + " element should have " + (reversibleRateCount) + " dimensions for reversible model or " + nonReversibleRateCount + " dimensions for non-reversible. " + "However parameter dimension is " + ratesParameter.getDimension());
        }
        if (hasIndicator) {
            // this is using BSSVS
            cxo = xo.getChild(INDICATOR);
            indicatorParameter = (Parameter) cxo.getChild(Parameter.class);
            if (indicatorParameter.getDimension() != ratesParameter.getDimension()) {
                throw new XMLParseException("Rates and indicator parameters in " + getParserName() + " element must be the same dimension.");
            }
            boolean randomize = xo.getAttribute(ComplexSubstitutionModelParser.RANDOMIZE, false);
            if (randomize) {
                BayesianStochasticSearchVariableSelection.Utils.randomize(indicatorParameter, dataType.getStateCount(), !isNonReversible);
            }
        }
        if (isNonReversible) {
            Logger.getLogger("dr.evomodel").info("  Using BSSVS Complex Substitution Model");
            return new SVSComplexSubstitutionModel(getParserName(), dataType, freqModel, ratesParameter, indicatorParameter);
        } else {
            Logger.getLogger("dr.evomodel").info("  Using BSSVS General Substitution Model");
            return new SVSGeneralSubstitutionModel(dataType, freqModel, ratesParameter, indicatorParameter);
        }
    } else {
        if (ratesParameter.getDimension() != reversibleRateCount - 1) {
            throw new XMLParseException("Rates parameter in " + getParserName() + " element should have " + (reversibleRateCount - 1) + " dimensions. However parameter dimension is " + ratesParameter.getDimension());
        }
        int relativeTo = 0;
        if (hasRelativeRates) {
            relativeTo = cxo.getIntegerAttribute(RELATIVE_TO) - 1;
        }
        if (relativeTo < 0 || relativeTo >= reversibleRateCount) {
            throw new XMLParseException(RELATIVE_TO + " must be 1 or greater");
        } else {
            int t = relativeTo;
            int s = states - 1;
            int row = 0;
            while (t >= s) {
                t -= s;
                s -= 1;
                row += 1;
            }
            int col = t + row + 1;
            Logger.getLogger("dr.evomodel").info("  Rates relative to " + dataType.getCode(row) + "<->" + dataType.getCode(col));
        }
        if (ratesParameter == null) {
            if (reversibleRateCount == 1) {
            // simplest model for binary traits...
            } else {
                throw new XMLParseException("No rates parameter found in " + getParserName());
            }
        }
        return new GeneralSubstitutionModel(dataType, freqModel, ratesParameter, relativeTo);
    }
}

Example 65 with DataType

use of dr.evolution.datatype.DataType in project beast-mcmc by beast-dev.

the class TimeIrreversibleTest method testSVSComplexSubstitutionModel.

private double[] testSVSComplexSubstitutionModel(Original test, double[] rates) {
    System.out.println("\n*** SVS Complex Substitution Model Test: " + test + " ***");
    double[] indicators = test.getIndicators();
    Parameter ratesP = new Parameter.Default(rates);
    Parameter indicatorsP = new Parameter.Default(indicators);
    DataType dataType = test.getDataType();
    FrequencyModel freqModel = new FrequencyModel(dataType, new Parameter.Default(test.getFrequencies()));
    SVSComplexSubstitutionModel substModel = new SVSComplexSubstitutionModel("SVS Complex Substitution Model Test", dataType, freqModel, ratesP, indicatorsP);
    double logL = substModel.getLogLikelihood();
    System.out.println("Prior = " + logL);
    double[] finiteTimeProbs = null;
    if (!Double.isInfinite(logL)) {
        finiteTimeProbs = new double[substModel.getDataType().getStateCount() * substModel.getDataType().getStateCount()];
        substModel.getTransitionProbabilities(time, finiteTimeProbs);
        System.out.println("Probs = ");
        printRateMatrix(finiteTimeProbs, substModel.getDataType().getStateCount());
    }
    // assertEquals(1, 1, 1e-10);
    return finiteTimeProbs;
}