Examples with Vector dr.math.matrixAlgebra.Vector used on opensource projects

Example 11 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class KDEDistributionTest method testNormalKDE.

public void testNormalKDE() {
    int length = 100;
    Double[] sample = new Double[length];
    for (int i = 0; i < length; i++) {
        sample[i] = (double) (i + 1);
    }
    //        final int gridSize = 256;
    NormalKDEDistribution kde = new NormalKDEDistribution(sample, null, null, null);
    System.out.println("bw = " + kde.getBandWidth());
    assertEquals(rBandWidth[0], kde.getBandWidth(), tolerance);
    final int gridSize = NormalKDEDistribution.MINIMUM_GRID_SIZE;
    double[] testPoints = new double[gridSize];
    double x = kde.getFromPoint();
    double delta = (kde.getToPoint() - kde.getFromPoint()) / (gridSize - 1);
    for (int i = 0; i < gridSize; ++i) {
        testPoints[i] = x;
        x += delta;
    }
    System.err.println("Eval @ " + new Vector(testPoints));
    double[] testDensity = new double[gridSize];
    for (int i = 0; i < gridSize; ++i) {
        testDensity[i] = kde.pdf(testPoints[i]);
    }
    System.err.println("Den    " + new Vector(testDensity));
    System.err.println("den[0] = " + testDensity[0]);
    System.err.println("den[N] = " + testDensity[NormalKDEDistribution.MINIMUM_GRID_SIZE - 1]);
//  System.exit(-1);
}

Example 12 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class MultiEpochExponential method getAnalyticIntegral.

public double getAnalyticIntegral(double start, double finish) {
    if (start == finish) {
        return 0.0;
    }
    double integral = 0.0;
    double logDemographic = 0.0;
    double lastTransitionTime = 0.0;
    int currentEpoch = 0;
    // Account for all epochs before start
    while (currentEpoch < transitionTime.length && start > transitionTime[currentEpoch]) {
        logDemographic += -rate[currentEpoch] * (transitionTime[currentEpoch] - lastTransitionTime);
        lastTransitionTime = transitionTime[currentEpoch];
        ++currentEpoch;
    }
    // Account for all epochs before finish
    while (currentEpoch < transitionTime.length && finish > transitionTime[currentEpoch]) {
        // Add to integral
        double incr = 0.0;
        if (rate[currentEpoch] == 0.0) {
            integral += incr = integrateConstant(start, transitionTime[currentEpoch], logDemographic);
        } else {
            integral += incr = integrateExponential(start - lastTransitionTime, transitionTime[currentEpoch] - lastTransitionTime, logDemographic, rate[currentEpoch]);
        }
        //            System.err.println("begin incr = " + incr + " for " + start + " -> " + transitionTime[currentEpoch] + " or " +
        //                    (start - lastTransitionTime) + " -> " + (transitionTime[currentEpoch] - lastTransitionTime) + " @ " + rate[currentEpoch] + " & " + Math.exp(logDemographic));
        // Update demographic function
        logDemographic += -rate[currentEpoch] * (transitionTime[currentEpoch] - lastTransitionTime);
        lastTransitionTime = transitionTime[currentEpoch];
        start = lastTransitionTime;
        ++currentEpoch;
    }
    // End of integral
    double incr = 0.0;
    if (rate[currentEpoch] == 0.0) {
        integral += incr = integrateConstant(start, finish, logDemographic);
    } else {
        integral += incr = integrateExponential(start - lastTransitionTime, finish - lastTransitionTime, logDemographic, rate[currentEpoch]);
    }
    if (Double.isNaN(integral) || Double.isInfinite(integral)) {
        System.err.println(integral + " " + start + " " + finish + new Vector(rate) + "\n");
    }
    return integral / getN0();
}

Example 13 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class StateHistoryTest method testFreqDistribution.

public void testFreqDistribution() {
    System.out.println("Start of FreqDistribution test");
    int startingState = 0;
    // 10 expected substitutions is close to \infty
    double duration = 10;
    double[] freq = new double[stateCount];
    for (int i = 0; i < N; i++) {
        StateHistory simultant = StateHistory.simulateUnconditionalOnEndingState(0.0, startingState, duration, lambda, stateCount);
        freq[simultant.getEndingState()]++;
    }
    for (int i = 0; i < stateCount; i++) {
        freq[i] /= N;
    }
    System.out.println("freq = " + new Vector(freq));
    assertEquals(freq, freqModel.getFrequencies(), 1E-3);
    System.out.println("End of FreqDistribution test\n");
}

Example 14 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class StateHistoryTest method setUp.

public void setUp() {
    MathUtils.setSeed(666);
    freqModel = new FrequencyModel(Nucleotides.INSTANCE, new double[] { 0.45, 0.25, 0.05, 0.25 });
    baseModel = new HKY(2.0, freqModel);
    stateCount = baseModel.getDataType().getStateCount();
    lambda = new double[stateCount * stateCount];
    baseModel.getInfinitesimalMatrix(lambda);
    System.out.println("lambda = " + new Vector(lambda));
    markovjumps = new MarkovJumpsSubstitutionModel(baseModel);
}

Example 15 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class StateHistoryTest method testCounts.

public void testCounts() {
    System.out.println("State of Counts test");
    int startingState = 2;
    double duration = 0.5;
    int[] counts = new int[stateCount * stateCount];
    double[] expectedCounts = new double[stateCount * stateCount];
    for (int i = 0; i < N; i++) {
        StateHistory simultant = StateHistory.simulateUnconditionalOnEndingState(0.0, startingState, duration, lambda, stateCount);
        simultant.accumulateSufficientStatistics(counts, null);
    }
    for (int i = 0; i < stateCount * stateCount; i++) {
        expectedCounts[i] = (double) counts[i] / (double) N;
    }
    double[] r = new double[stateCount * stateCount];
    double[] joint = new double[stateCount * stateCount];
    double[] analytic = new double[stateCount * stateCount];
    for (int from = 0; from < stateCount; from++) {
        for (int to = 0; to < stateCount; to++) {
            double marginal = 0;
            if (from != to) {
                MarkovJumpsCore.fillRegistrationMatrix(r, from, to, stateCount);
                markovjumps.setRegistration(r);
                markovjumps.computeJointStatMarkovJumps(duration, joint);
                for (int j = 0; j < stateCount; j++) {
                    // Marginalize out ending state
                    marginal += joint[startingState * stateCount + j];
                }
            }
            analytic[from * stateCount + to] = marginal;
        }
    }
    System.out.println("unconditional expected counts = " + new Vector(expectedCounts));
    System.out.println("analytic               counts = " + new Vector(analytic));
    assertEquals(expectedCounts, analytic, 1E-3);
    System.out.println("End of Counts test\n");
}