Java SitePatterns类代码示例

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public double doOperation() {

        // Select boundary to update, 0 => btw partition 0 and 1,  1 => btw partition 1 and 2, etc.
        int whichBoundary = MathUtils.nextInt(breakPoints.getDimension());

        int cBreakPt = (int) breakPoints.getParameterValue(whichBoundary);
        SitePatterns left = partitions.get(whichBoundary);
        SitePatterns right = partitions.get(whichBoundary + 1);
        int min = left.getFrom();
        int max = right.getTo();

        int pBreakPt = min;
        while (pBreakPt <= min || pBreakPt >= max) {        // cBreakPt + [windowSize-1, windowSize+1] (and not 0)
            if (MathUtils.nextBoolean())
                pBreakPt = cBreakPt + MathUtils.nextInt(windowSize) + 1;
            else
                pBreakPt = cBreakPt - MathUtils.nextInt(windowSize) - 1;
        }


        return 0;
    }
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
protected double[] computeSitePatternLikelihoods(SitePatterns patterns) {
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 29.739445, 0, 100);

    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);

    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);

    //treeLikelihood
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null,
            false, false, true, false, false);

    return treeLikelihood.getPatternLogLikelihoods();
}
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public static double[] getAlignmentScore(ScoreMatrix scoreMatrix, SitePatterns sitePatterns) {

        AlignmentScore alignmentScore = new AlignmentScore(scoreMatrix, sitePatterns);

        double[] params = new double[alignmentScore.getNumArguments()];
        DifferentialEvolution de = new DifferentialEvolution(params.length,params.length*10);
        for (int i = 0; i < params.length; i++) params[i] = 0.5;

        de.optimize(alignmentScore, params,1e-6, 1e-6);
        double score = alignmentScore.evaluate(params);

        double[] results = new double[params.length + 1];
        System.arraycopy(params, 0, results, 0, params.length);
        results[params.length] = score;

        return results;
    }
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public final double getScore(SitePatterns patterns) {

        if (patterns.getPatternLength() != 2) throw new IllegalArgumentException();

        double logL = 0.0;
        for (int i = 0; i < patterns.getPatternCount(); i++) {
            double weight = patterns.getPatternWeight(i);
            int[] pattern = patterns.getPattern(i);

            int x = pattern[0];
            int y = pattern[1];

            if (x < stateCount && y < stateCount) {
                logL += getScore(x,y) * weight;
            }

        }
        return logL;
    }
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public SlidingPatternsOperator(List<SitePatterns> list, Parameter breakPoints, int windowSize, int weight, CoercionMode mode) {
        super(mode);
        this.partitions = list;
        this.windowSize = windowSize;
//		this.weight = weight;
//		this.mode = mode;
        this.breakPoints = breakPoints;
//        System.err.println("Starting values: "+currentBreakPointsString() );
//        System.exit(-1) ;
    }
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public static boolean arePartitionsContiguous(List<SitePatterns> list) {
        int current = -1;
        int index = 0;
        for (SitePatterns patterns : list) {
            int start = patterns.getFrom();
            int end = patterns.getTo();
            /* System.err.println(start+" -> "+end+" : "+patterns.getSiteCount()+" "+patterns.getPatternCount());
                           int[] data = patterns.getSitePattern(0);
                           System.err.print("Data 0:");
                           for (int i : data)
                               System.err.print(" "+i);
                           System.err.println("");*/
//            if (current == -1)
//                current = end;
            if (current != -1 && start != (current + 1))
//                throw new NonContiguousPartitionsException("Partition #"+0+" does not start contiguously");
                return false;
            current = end;
        }

        return true;
    }
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public Object parseXMLObject(XMLObject xo) throws XMLParseException {

        SitePatterns patternList = (SitePatterns) xo.getChild(PatternList.class);
        TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
        BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
        Alignment alignment = (Alignment) xo.getChild(Alignment.class);
        GammaSiteRateModel substitutionModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);

        PrintWriter branch = null, operation = null;

        if (xo.hasAttribute(BRANCH_FILE_NAME)) {
            branch = XMLParser.getFilePrintWriter(xo, OPERATION_REPORT, BRANCH_FILE_NAME);
        }
        if (xo.hasAttribute(OPERATION_FILE_NAME)) {
            operation = XMLParser.getFilePrintWriter(xo, OPERATION_REPORT, OPERATION_FILE_NAME);
        }

        return new BeagleOperationReport(treeModel, patternList, rateModel, substitutionModel, alignment, branch, operation);
    }
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
private void tryAllPossibleAlignments(int length, SitePatterns patterns) {
    System.out.println("Trying all possible " + length + " site alignments");
    double[] patternLogLikeihoods = computeSitePatternLikelihoods(patterns);

    System.out.println("Site logLikelihoods: " + new Vector(patternLogLikeihoods));

    double total = 0.0;

    List<Double> allAlignmentLogProbabilities = new ArrayList<Double>();
    recursivelyComputeAlignmentLikelihood(allAlignmentLogProbabilities, patternLogLikeihoods, length, 0, 0);

    System.out.println("Total possible alignments: " + allAlignmentLogProbabilities.size());

    for (Double x : allAlignmentLogProbabilities) {
        total += Math.exp(x);
    }
    System.out.println("Total probability = " + total);
    assertEquals("uncorrected", 1.0, total, tolerance);
}
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public void testLikelihoodJC69() {
    System.out.println("\nTest Likelihood using JC69:");
    // Sub model
    Parameter freqs = new Parameter.Default(new double[]{0.25, 0.25, 0.25, 0.25});
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, 100);

    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);

    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);

    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);

    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null,
            false, false, true, false, false);

    assertEquals("treeLikelihoodJC69", format.format(-1992.20564), format.format(treeLikelihood.getLogLikelihood()));
}
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public void testLikelihoodK80() {
    System.out.println("\nTest Likelihood using K80:");
    // Sub model
    Parameter freqs = new Parameter.Default(new double[]{0.25, 0.25, 0.25, 0.25});
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 27.402591, 0, 100);

    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);

    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);

    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);

    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null,
            false, false, true, false, false);

    assertEquals("treeLikelihoodK80", format.format(-1856.30305), format.format(treeLikelihood.getLogLikelihood()));
}
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public void testLikelihoodHKY85() {
    System.out.println("\nTest Likelihood using HKY85:");
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 29.739445, 0, 100);

    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);

    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);

    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);

    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null,
            false, false, true, false, false);          

    assertEquals("treeLikelihoodHKY85", format.format(-1825.21317), format.format(treeLikelihood.getLogLikelihood()));
}
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public void testLikelihoodHKY85G() {
    System.out.println("\nTest Likelihood using HKY85G:");
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 38.829740, 0, 100);

    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);

    //siteModel
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    Parameter shape = new Parameter.Default(GammaSiteModelParser.GAMMA_SHAPE, 0.137064, 0, 1000.0);

    GammaSiteModel siteModel = new GammaSiteModel(hky, mu, shape, 4, null);

    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);

    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null,
            false, false, true, false, false);

    assertEquals("treeLikelihoodHKY85G", format.format(-1789.75936), format.format(treeLikelihood.getLogLikelihood()));
}
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public void testLikelihoodHKY85I() {
    System.out.println("\nTest Likelihood using HKY85I:");
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 38.564672, 0, 100);

    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);

    //siteModel
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    Parameter invar = new Parameter.Default(GammaSiteModelParser.PROPORTION_INVARIANT, 0.701211, 0, 1.0);

    GammaSiteModel siteModel = new GammaSiteModel(hky, mu, null, 4, invar);

    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);

    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null,
            false, false, true, false, false);

    assertEquals("treeLikelihoodHKY85I", format.format(-1789.91240), format.format(treeLikelihood.getLogLikelihood()));
}
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public void testLikelihoodHKY85GI() {
    System.out.println("\nTest Likelihood using HKY85GI:");
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 39.464538, 0, 100);

    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);

    //siteModel
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    Parameter shape = new Parameter.Default(GammaSiteModelParser.GAMMA_SHAPE, 0.587649, 0, 1000.0);
    Parameter invar = new Parameter.Default(GammaSiteModelParser.PROPORTION_INVARIANT, 0.486548, 0, 1.0);

    GammaSiteModel siteModel = new GammaSiteModel(hky, mu, shape, 4, invar);

    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);

    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null,
            false, false, true, false, false);

    assertEquals("treeLikelihoodHKY85GI", format.format(-1789.63923), format.format(treeLikelihood.getLogLikelihood()));
}
 

import dr.evolution.alignment.SitePatterns; //导入依赖的package包/类
public double doOperation() throws OperatorFailedException {

        // Select boundary to update, 0 => btw partition 0 and 1,  1 => btw partition 1 and 2, etc.
        int whichBoundary = MathUtils.nextInt(breakPoints.getDimension());

        int cBreakPt = (int) breakPoints.getParameterValue(whichBoundary);
        SitePatterns left = partitions.get(whichBoundary);
        SitePatterns right = partitions.get(whichBoundary + 1);
        int min = left.getFrom();
        int max = right.getTo();

        int pBreakPt = min;
        while (pBreakPt <= min || pBreakPt >= max) {        // cBreakPt + [windowSize-1, windowSize+1] (and not 0)
            if (MathUtils.nextBoolean())
                pBreakPt = cBreakPt + MathUtils.nextInt(windowSize) + 1;
            else
                pBreakPt = cBreakPt - MathUtils.nextInt(windowSize) - 1;
        }


        return 0;
    }