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

Example 6 with SubstitutionModel

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

the class EpochTreeBranchSubstitutionModel method getTransitionProbabilities.

public void getTransitionProbabilities(Tree tree, NodeRef node, int rateCategory, double[] matrix) {
    NodeRef parent = tree.getParent(node);
    final double branchRate = branchModel.getBranchRate(tree, node);
    // Get the operational time of the branch
    final double startTime = tree.getNodeHeight(parent);
    final double endTime = tree.getNodeHeight(node);
    final double branchTime = branchRate * (startTime - endTime);
    if (branchTime < 0.0) {
        throw new RuntimeException("Negative branch length: " + branchTime);
    }
    double distance = siteModel.getRateForCategory(rateCategory) * branchTime;
    int matrixCount = 0;
    boolean oneMatrix = (getEpochWeights(startTime, endTime, weight) == 1);
    for (int m = 0; m < numberModels; m++) {
        if (weight[m] > 0) {
            SubstitutionModel model = modelList.get(m);
            if (matrixCount == 0) {
                if (oneMatrix) {
                    model.getTransitionProbabilities(distance, matrix);
                    break;
                } else
                    model.getTransitionProbabilities(distance * weight[m], resultMatrix);
                matrixCount++;
            } else {
                model.getTransitionProbabilities(distance * weight[m], stepMatrix);
                // Sum over unobserved state
                int index = 0;
                for (int i = 0; i < stateCount; i++) {
                    for (int j = 0; j < stateCount; j++) {
                        productMatrix[index] = 0;
                        for (int k = 0; k < stateCount; k++) {
                            productMatrix[index] += resultMatrix[i * stateCount + k] * stepMatrix[k * stateCount + j];
                        }
                        index++;
                    }
                }
                // Swap pointers
                double[] tmpMatrix = resultMatrix;
                resultMatrix = productMatrix;
                productMatrix = tmpMatrix;
            }
        }
    }
    if (!oneMatrix)
        System.arraycopy(resultMatrix, 0, matrix, 0, stateCount * stateCount);
}

Example 7 with SubstitutionModel

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

the class AlignmentScore method main.

public static void main(String[] args) throws java.io.IOException, Importer.ImportException {
    NexusImporter importer = new NexusImporter(new FileReader(args[0]));
    Alignment alignment = importer.importAlignment();
    ExtractPairs pairs = new ExtractPairs(alignment);
    Parameter muParam = new Parameter.Default(1.0);
    Parameter kappaParam = new Parameter.Default(1.0);
    kappaParam.addBounds(new Parameter.DefaultBounds(100.0, 0.0, 1));
    muParam.addBounds(new Parameter.DefaultBounds(1.0, 1.0, 1));
    Parameter freqParam = new Parameter.Default(alignment.getStateFrequencies());
    FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqParam);
    SubstitutionModel substModel = new HKY(kappaParam, freqModel);
    SiteModel siteModel = new GammaSiteModel(substModel, muParam, null, 1, null);
    ScoreMatrix scoreMatrix = new ScoreMatrix(siteModel, 0.1);
    double threshold = 0.1;
    List<PairDistance> pairDistances = new ArrayList<PairDistance>();
    Set<Integer> sequencesUsed = new HashSet<Integer>();
    List<Integer> allGaps = new ArrayList<Integer>();
    for (int i = 0; i < alignment.getSequenceCount(); i++) {
        for (int j = i + 1; j < alignment.getSequenceCount(); j++) {
            Alignment pairAlignment = pairs.getPairAlignment(i, j);
            if (pairAlignment != null) {
                SitePatterns patterns = new SitePatterns(pairAlignment);
                double distance = getGeneticDistance(scoreMatrix, patterns);
                if (distance < threshold) {
                    List gaps = new ArrayList();
                    GapUtils.getGapSizes(pairAlignment, gaps);
                    pairDistances.add(new PairDistance(i, j, distance, gaps, pairAlignment.getSiteCount()));
                    System.out.print(".");
                } else {
                    System.out.print("*");
                }
            } else {
                System.out.print("x");
            }
        }
        System.out.println();
    }
    Collections.sort(pairDistances);
    int totalLength = 0;
    for (PairDistance pairDistance : pairDistances) {
        Integer x = pairDistance.x;
        Integer y = pairDistance.y;
        if (!sequencesUsed.contains(x) && !sequencesUsed.contains(y)) {
            allGaps.addAll(pairDistance.gaps);
            sequencesUsed.add(x);
            sequencesUsed.add(y);
            System.out.println("Added pair (" + x + "," + y + ") d=" + pairDistance.distance + " L=" + pairDistance.alignmentLength);
            totalLength += pairDistance.alignmentLength;
        }
    }
    printFrequencyTable(allGaps);
    System.out.println("total length=" + totalLength);
}

Example 8 with SubstitutionModel

use of dr.oldevomodel.substmodel.SubstitutionModel 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 9 with SubstitutionModel

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

the class ALSSiteModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    SubstitutionModel substitutionModel = (SubstitutionModel) xo.getElementFirstChild(SUBSTITUTION_MODEL);
    String msg = "";
    Parameter muParam = null;
    if (xo.hasChildNamed(SUBSTITUTION_RATE)) {
        muParam = (Parameter) xo.getElementFirstChild(SUBSTITUTION_RATE);
        msg += "\n  with initial substitution rate = " + muParam.getParameterValue(0);
    } else if (xo.hasChildNamed(MUTATION_RATE)) {
        muParam = (Parameter) xo.getElementFirstChild(MUTATION_RATE);
        msg += "\n  with initial substitution rate = " + muParam.getParameterValue(0);
    } else if (xo.hasChildNamed(RELATIVE_RATE)) {
        muParam = (Parameter) xo.getElementFirstChild(RELATIVE_RATE);
        msg += "\n  with initial relative rate = " + muParam.getParameterValue(0);
    }
    Parameter shapeParam = null;
    int catCount = 4;
    if (xo.hasChildNamed(GAMMA_SHAPE)) {
        final XMLObject cxo = xo.getChild(GAMMA_SHAPE);
        catCount = cxo.getIntegerAttribute(GAMMA_CATEGORIES);
        shapeParam = (Parameter) cxo.getChild(Parameter.class);
        msg += "\n  " + catCount + " category discrete gamma with initial shape = " + shapeParam.getParameterValue(0);
    }
    Parameter invarParam = null;
    if (xo.hasChildNamed(PROPORTION_INVARIANT)) {
        invarParam = (Parameter) xo.getElementFirstChild(PROPORTION_INVARIANT);
        msg += "\n  initial proportion of invariant sites = " + invarParam.getParameterValue(0);
    }
    Logger.getLogger("dr.evomodel").info("Creating site model." + (msg.length() > 0 ? msg : ""));
    return new GammaSiteModel(substitutionModel, muParam, shapeParam, catCount, invarParam);
}

Example 10 with SubstitutionModel

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

the class CategorySiteModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    XMLObject cxo = xo.getChild(SUBSTITUTION_MODEL);
    SubstitutionModel substitutionModel = (SubstitutionModel) cxo.getChild(SubstitutionModel.class);
    cxo = xo.getChild(MUTATION_RATE);
    Parameter muParam = (Parameter) cxo.getChild(Parameter.class);
    cxo = xo.getChild(RATE_PARAMETER);
    Parameter rateParam = null;
    int relativeTo = 0;
    if (cxo != null) {
        rateParam = (Parameter) cxo.getChild(Parameter.class);
        relativeTo = cxo.getIntegerAttribute(RELATIVE_TO);
    }
    cxo = xo.getChild(CATEGORIES);
    String categories = "";
    String states = "";
    if (cxo != null) {
        categories = cxo.getStringAttribute(CATEGORY_STRING);
        states = cxo.getStringAttribute(CATEGORY_STATES);
    }
    return new CategorySiteModel(substitutionModel, muParam, rateParam, categories, states, relativeTo);
}