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

Example 1 with SymmetricMatrix

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

the class FactorGibbsOperator method doOperation.

@Override
public double doOperation() {
    setupParameters();
    getPrecision(precision);
    double[][] variance = (new SymmetricMatrix(precision)).inverse().toComponents();
    if (randomScan) {
        int i = MathUtils.nextInt(LFM.getFactors().getColumnDimension());
        randomDraw(i, variance);
    }
    for (int i = 0; i < LFM.getFactors().getColumnDimension(); i++) {
        randomDraw(i, variance);
    }
    LFM.getFactors().fireParameterChangedEvent();
    //To change body of implemented methods use File | Settings | File Templates.
    return 0;
}

Example 2 with SymmetricMatrix

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

the class FactorOperator method doOperation.

@Override
public double doOperation() {
    setupParameters();
    getPrecision(precision);
    double[][] variance = (new SymmetricMatrix(precision)).inverse().toComponents();
    if (randomScan) {
        int i = MathUtils.nextInt(LFM.getFactors().getColumnDimension());
        randomDraw(i, variance);
    }
    for (int i = 0; i < LFM.getFactors().getColumnDimension(); i++) {
        randomDraw(i, variance);
    }
    LFM.getFactors().fireParameterChangedEvent();
    //To change body of implemented methods use File | Settings | File Templates.
    return 0;
}

Example 3 with SymmetricMatrix

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

the class PrecisionMatrixGibbsOperator method getOperationScaleMatrixAndSetObservationCount.

private double[][] getOperationScaleMatrixAndSetObservationCount() {
    // calculate sum-of-the-weighted-squares matrix over tree
    double[][] S = new double[dim][dim];
    SymmetricMatrix S2;
    SymmetricMatrix inverseS2 = null;
    // Need to reset, as incrementOuterProduct can be recursive
    numberObservations = 0;
    if (isSampledTraitLikelihood) {
        incrementOuterProduct(S, treeModel.getRoot());
    } else {
        // IntegratedTraitLikelihood
        if (traitModel != null) {
            // is a tree
            // TODO deprecate usage
            incrementOuterProduct(S, (ConjugateWishartStatisticsProvider) traitModel);
        } else if (conjugateWishartProvider != null) {
            incrementOuterProduct(S, conjugateWishartProvider);
        } else {
            // is a normal-normal-wishart model
            incrementOuterProduct(S, multivariateLikelihood);
        }
    }
    try {
        S2 = new SymmetricMatrix(S);
        if (pathWeight != 1.0) {
            S2 = (SymmetricMatrix) S2.product(pathWeight);
        }
        if (priorInverseScaleMatrix != null)
            S2 = priorInverseScaleMatrix.add(S2);
        inverseS2 = (SymmetricMatrix) S2.inverse();
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    assert inverseS2 != null;
    return inverseS2.toComponents();
}

Example 4 with SymmetricMatrix

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

the class TraitGibbsOperator method operateRoot.

private MeanPrecision operateRoot(NodeRef node) {
    double[] trait;
    double weightTotal = 0.0;
    double[] weightedAverage = new double[dim];
    double[][] precision = precisionMatrixParameter.getParameterAsMatrix();
    for (int k = 0; k < treeModel.getChildCount(node); k++) {
        NodeRef child = treeModel.getChild(node, k);
        trait = treeModel.getMultivariateNodeTrait(child, traitName);
        final double weight = 1.0 / traitModel.getRescaledBranchLengthForPrecision(child);
        for (int i = 0; i < dim; i++) {
            for (int j = 0; j < dim; j++) weightedAverage[i] += precision[i][j] * weight * trait[j];
        }
        weightTotal += weight;
    }
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {
            weightedAverage[i] += rootPriorPrecision[i][j] * rootPriorMean[j];
            precision[i][j] = precision[i][j] * weightTotal + rootPriorPrecision[i][j];
        }
    }
    double[][] variance = new SymmetricMatrix(precision).inverse().toComponents();
    trait = new double[dim];
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) trait[i] += variance[i][j] * weightedAverage[j];
    }
    return new MeanPrecision(trait, precision);
}

Example 5 with SymmetricMatrix

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

the class LatentLiabilityGibbs method getComponentConditionalMean.

private double[] getComponentConditionalMean(double[][] thisP, double[] oldValue, double[] mean, double[][] condP) {
    double[] condMean = new double[numUpdate];
    double[][] H = new double[numUpdate][numFixed];
    Matrix prod = new Matrix(numUpdate, numFixed);
    Vector dif = new Vector(numUpdate);
    double[] contMeans = new double[numFixed];
    for (int i = 0; i < numUpdate; i++) {
        for (int j = 0; j < numFixed; j++) {
            H[i][j] = thisP[doUpdate[i]][dontUpdate[j]];
        }
    }
    for (int i = 0; i < numFixed; i++) {
        contMeans[i] = oldValue[dontUpdate[i]] - mean[dontUpdate[i]];
    }
    Matrix invK = new SymmetricMatrix(condP).inverse();
    Matrix HH = new Matrix(H);
    try {
        prod = invK.product(HH);
        dif = prod.product(new Vector(contMeans));
    } catch (IllegalDimension illegalDimension) {
        illegalDimension.printStackTrace();
    }
    for (int i = 0; i < numUpdate; i++) {
        condMean[i] = mean[doUpdate[i]] - dif.component(i);
    }
    return condMean;
}