Examples with SubstitutionModel dr.evomodel.substmodel.SubstitutionModel used on opensource projects

Example 1 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel 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 SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class BranchAssignmentModel method getRootSubstitutionModel.

// END: getSubstitutionModels
@Override
public SubstitutionModel getRootSubstitutionModel() {
    Object nodeAttribute = treeModel.getNodeAttribute(treeModel.getRoot(), BranchAssignmentModelParser.ANNOTATION_VALUE);
    SubstitutionModel model = null;
    if (nodeAttribute == null) {
        model = this.baseModel;
    } else {
        Integer modelIndex = (Integer) nodeAttribute;
        model = this.modelIndexMap.get(modelIndex);
    }
    return model;
}

Example 3 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class BranchAssignmentModel method setup.

// END: Constructor
private void setup() {
    for (NodeRef node : this.treeModel.getNodes()) {
        if (!treeModel.isRoot(node)) {
            Integer modelIndex = Integer.MAX_VALUE;
            SubstitutionModel model = null;
            Object nodeAttribute = treeModel.getNodeAttribute(node, annotation);
            if (nodeAttribute == null) {
                System.out.println("Attribute " + annotation + " missing from node. Using base model as branch model.");
                modelIndex = this.baseModelIndex;
                model = this.baseModel;
            } else {
                modelIndex = (Integer) nodeAttribute;
                model = this.modelIndexMap.get(modelIndex);
            }
            branchAssignmentMap.put(node, modelIndex);
            // if (substitutionModels.get(modelIndex) == null) {
            // substitutionModels.set(modelIndex, model);
            // }
            substitutionModels.add(model);
        }
    // END: root check
    }
// END: nodes loop
}

Example 4 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class LineageSpecificBranchModelParser method parseXMLObject.

@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    FrequencyModel rootFrequencyModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
    XMLObject cxo = xo.getChild(MODELS);
    for (int i = 0; i < cxo.getChildCount(); i++) {
        SubstitutionModel substModel = (SubstitutionModel) cxo.getChild(i);
        substitutionModels.add(substModel);
    }
    // END: models loop
    // TODO: check if categories numbering starts from zero
    Parameter categories = (Parameter) xo.getElementFirstChild(CATEGORIES);
    return new // provider,
    LineageSpecificBranchModel(// provider,
    treeModel, // provider,
    rootFrequencyModel, // provider,
    substitutionModels, categories);
}

Example 5 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class DiscreteTraitBranchRateDelegate method cacheDifferentialMassMatrix.

protected void cacheDifferentialMassMatrix(Tree tree, boolean cacheSquaredMatrix) {
    for (int i = 0; i < evolutionaryProcessDelegate.getSubstitutionModelCount(); i++) {
        double[] infinitesimalMatrix = new double[stateCount * stateCount];
        SubstitutionModel substitutionModel = evolutionaryProcessDelegate.getSubstitutionModel(i);
        substitutionModel.getInfinitesimalMatrix(infinitesimalMatrix);
        // if (stateCount > 4) {
        // double[] transpose = new double[stateCount * stateCount];
        // for (int row = 0; row < stateCount; ++row) {
        // for (int col = 0; col < stateCount; ++col) {
        // transpose[col * stateCount + row] = infinitesimalMatrix[row * stateCount + col];
        // }
        // }
        // 
        // infinitesimalMatrix = transpose;
        // }
        double[] scaledInfinitesimalMatrix = scaleInfinitesimalMatrixByRates(infinitesimalMatrix, DifferentialChoice.GRADIENT, siteRateModel);
        evolutionaryProcessDelegate.cacheInfinitesimalMatrix(beagle, i, scaledInfinitesimalMatrix);
        if (cacheSquaredMatrix) {
            double[] infinitesimalMatrixSquared = new double[stateCount * stateCount];
            for (int l = 0; l < stateCount; l++) {
                for (int j = 0; j < stateCount; j++) {
                    double sumOverState = 0.0;
                    for (int k = 0; k < stateCount; k++) {
                        sumOverState += infinitesimalMatrix[l * stateCount + k] * infinitesimalMatrix[k * stateCount + j];
                    }
                    infinitesimalMatrixSquared[l * stateCount + j] = sumOverState;
                }
            }
            double[] scaledInfinitesimalMatrixSquared = scaleInfinitesimalMatrixByRates(infinitesimalMatrixSquared, DifferentialChoice.HESSIAN, siteRateModel);
            evolutionaryProcessDelegate.cacheInfinitesimalSquaredMatrix(beagle, i, scaledInfinitesimalMatrixSquared);
        }
    }
}