Examples with Matrix Jama.Matrix used on opensource projects

Example 11 with Matrix

use of Jama.Matrix in project h2o-3 by h2oai.

the class LinearAlgebraUtils method computeR.

/**
   * Get R = L' from Cholesky decomposition Y'Y = LL' (same as R from Y = QR)
   * @param jobKey Job key for Gram calculation
   * @param yinfo DataInfo for Y matrix
   * @param transpose Should result be transposed to get L?
   * @return L or R matrix from Cholesky of Y Gram
   */
public static double[][] computeR(Key<Job> jobKey, DataInfo yinfo, boolean transpose) {
    // Calculate Cholesky of Y Gram to get R' = L matrix
    // Gram is Y'Y/n where n = nrow(Y)
    Gram.GramTask gtsk = new Gram.GramTask(jobKey, yinfo);
    gtsk.doAll(yinfo._adaptedFrame);
    // Gram.Cholesky chol = gtsk._gram.cholesky(null);   // If Y'Y = LL' Cholesky, then R = L'
    Matrix ygram = new Matrix(gtsk._gram.getXX());
    CholeskyDecomposition chol = new CholeskyDecomposition(ygram);
    double[][] L = chol.getL().getArray();
    // Must scale since Cholesky of Y'Y/n where nobs = nrow(Y)
    ArrayUtils.mult(L, Math.sqrt(gtsk._nobs));
    return transpose ? L : ArrayUtils.transpose(L);
}

Example 12 with Matrix

use of Jama.Matrix in project gatk by broadinstitute.

the class MultivariateGaussian method evaluateFinalModelParameters.

public void evaluateFinalModelParameters(final List<VariantDatum> data) {
    sumProb = 0.0;
    zeroOutMu();
    zeroOutSigma();
    int datumIndex = 0;
    for (final VariantDatum datum : data) {
        final double prob = pVarInGaussian[datumIndex++];
        sumProb += prob;
        incrementMu(datum, prob);
    }
    divideEqualsMu(sumProb);
    datumIndex = 0;
    final Matrix pVarSigma = new Matrix(mu.length, mu.length);
    for (final VariantDatum datum : data) {
        final double prob = pVarInGaussian[datumIndex++];
        for (int iii = 0; iii < mu.length; iii++) {
            for (int jjj = 0; jjj < mu.length; jjj++) {
                pVarSigma.set(iii, jjj, prob * (datum.annotations[iii] - mu[iii]) * (datum.annotations[jjj] - mu[jjj]));
            }
        }
        sigma.plusEquals(pVarSigma);
    }
    sigma.timesEquals(1.0 / sumProb);
    // clean up some memory
    resetPVarInGaussian();
}

Example 13 with Matrix

use of Jama.Matrix in project gatk by broadinstitute.

the class MultivariateGaussian method zeroOutSigma.

public void zeroOutSigma() {
    final double[][] zeroSigma = new double[mu.length][mu.length];
    for (final double[] row : zeroSigma) {
        Arrays.fill(row, 0);
    }
    final Matrix tmp = new Matrix(zeroSigma);
    sigma.setMatrix(0, mu.length - 1, 0, mu.length - 1, tmp);
}

Example 14 with Matrix

use of Jama.Matrix in project gatk by broadinstitute.

the class MultivariateGaussian method initializeRandomSigma.

public void initializeRandomSigma(final Random rand) {
    final double[][] randSigma = new double[mu.length][mu.length];
    for (int iii = 0; iii < mu.length; iii++) {
        for (int jjj = iii; jjj < mu.length; jjj++) {
            randSigma[jjj][iii] = 0.55 + 1.25 * rand.nextDouble();
            if (rand.nextBoolean()) {
                randSigma[jjj][iii] *= -1.0;
            }
            // Sigma is a symmetric, positive-definite matrix created by taking a lower diagonal matrix and multiplying it by its transpose
            if (iii != jjj) {
                randSigma[iii][jjj] = 0.0;
            }
        }
    }
    Matrix tmp = new Matrix(randSigma);
    tmp = tmp.times(tmp.transpose());
    sigma.setMatrix(0, mu.length - 1, 0, mu.length - 1, tmp);
}

Example 15 with Matrix

use of Jama.Matrix in project gatk by broadinstitute.

the class MultivariateGaussian method maximizeGaussian.

public void maximizeGaussian(final List<VariantDatum> data, final double[] empiricalMu, final Matrix empiricalSigma, final double SHRINKAGE, final double DIRICHLET_PARAMETER, final double DEGREES_OF_FREEDOM) {
    sumProb = 1E-10;
    final Matrix wishart = new Matrix(mu.length, mu.length);
    zeroOutMu();
    zeroOutSigma();
    int datumIndex = 0;
    for (final VariantDatum datum : data) {
        final double prob = pVarInGaussian[datumIndex++];
        sumProb += prob;
        incrementMu(datum, prob);
    }
    divideEqualsMu(sumProb);
    final double shrinkageFactor = (SHRINKAGE * sumProb) / (SHRINKAGE + sumProb);
    for (int iii = 0; iii < mu.length; iii++) {
        double deltaMu = shrinkageFactor * (mu[iii] - empiricalMu[iii]);
        for (int jjj = 0; jjj < mu.length; jjj++) {
            wishart.set(iii, jjj, deltaMu * (mu[jjj] - empiricalMu[jjj]));
        }
    }
    datumIndex = 0;
    final Matrix pVarSigma = new Matrix(mu.length, mu.length);
    for (final VariantDatum datum : data) {
        final double prob = pVarInGaussian[datumIndex++];
        for (int iii = 0; iii < mu.length; iii++) {
            double deltaMu = prob * (datum.annotations[iii] - mu[iii]);
            for (int jjj = 0; jjj < mu.length; jjj++) {
                pVarSigma.set(iii, jjj, deltaMu * (datum.annotations[jjj] - mu[jjj]));
            }
        }
        sigma.plusEquals(pVarSigma);
    }
    sigma.plusEquals(empiricalSigma);
    sigma.plusEquals(wishart);
    for (int iii = 0; iii < mu.length; iii++) {
        mu[iii] = (sumProb * mu[iii] + SHRINKAGE * empiricalMu[iii]) / (sumProb + SHRINKAGE);
    }
    hyperParameter_a = sumProb + DEGREES_OF_FREEDOM;
    hyperParameter_b = sumProb + SHRINKAGE;
    hyperParameter_lambda = sumProb + DIRICHLET_PARAMETER;
    // clean up some memory
    resetPVarInGaussian();
}