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

Example 11 with Matrix

use of dr.math.matrixAlgebra.Matrix 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;
}

Example 12 with Matrix

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

the class MultivariateTraitDebugUtilities method getJointVarianceFactor.

public static Matrix getJointVarianceFactor(final double priorSampleSize, double[][] treeVariance, double[][] treeSharedLengths, double[][] loadingsVariance, double[][] diffusionVariance, DiffusionProcessDelegate diffusionProcessDelegate, Matrix Lt) {
    if (!diffusionProcessDelegate.hasActualization()) {
        double[][] jointVariance = diffusionProcessDelegate.getJointVariance(priorSampleSize, treeVariance, treeVariance, loadingsVariance);
        Matrix loadingsFactorsVariance = new Matrix(jointVariance);
        return loadingsFactorsVariance;
    } else {
        double[][] jointVariance = diffusionProcessDelegate.getJointVariance(priorSampleSize, treeVariance, treeSharedLengths, diffusionVariance);
        Matrix jointVarianceMatrix = new Matrix(jointVariance);
        double[][] identity = KroneckerOperation.makeIdentityMatrixArray(treeSharedLengths[0].length);
        Matrix loadingsTree = new Matrix(KroneckerOperation.product(identity, Lt.toComponents()));
        Matrix loadingsFactorsVariance = null;
        try {
            loadingsFactorsVariance = loadingsTree.product(jointVarianceMatrix.product(loadingsTree.transpose()));
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
        return loadingsFactorsVariance;
    }
}

Example 13 with Matrix

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

the class ContinuousDataLikelihoodDelegate method getReport.

@Override
public String getReport() {
    StringBuilder sb = new StringBuilder();
    sb.append("Tree:\n");
    sb.append(callbackLikelihood.getId()).append("\t");
    sb.append(cdi.getReport());
    final Tree tree = callbackLikelihood.getTree();
    sb.append(tree.toString());
    sb.append("\n\n");
    final double normalization = rateTransformation.getNormalization();
    final double priorSampleSize = rootProcessDelegate.getPseudoObservations();
    double[][] treeStructure = MultivariateTraitDebugUtilities.getTreeVariance(tree, callbackLikelihood.getBranchRateModel(), 1.0, Double.POSITIVE_INFINITY);
    sb.append("Tree structure:\n");
    sb.append(new Matrix(treeStructure));
    sb.append("\n\n");
    double[][] treeSharedLengths = MultivariateTraitDebugUtilities.getTreeVariance(tree, callbackLikelihood.getBranchRateModel(), rateTransformation.getNormalization(), Double.POSITIVE_INFINITY);
    double[][] treeVariance = getTreeVariance();
    double[][] traitPrecision = getDiffusionModel().getPrecisionmatrix();
    Matrix traitVariance = new Matrix(traitPrecision).inverse();
    double[][] jointVariance = diffusionProcessDelegate.getJointVariance(priorSampleSize, treeVariance, treeSharedLengths, traitVariance.toComponents());
    if (dataModel instanceof RepeatedMeasuresTraitDataModel) {
        for (int tip = 0; tip < tipCount; ++tip) {
            double[] partial = dataModel.getTipPartial(tip, false);
            WrappedMatrix tipVariance = new WrappedMatrix.Raw(partial, dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
            for (int row = 0; row < dimTrait; ++row) {
                for (int col = 0; col < dimTrait; ++col) {
                    jointVariance[tip * dimTrait + row][tip * dimTrait + col] += tipVariance.get(row, col);
                }
            }
        }
    }
    Matrix treeV = new Matrix(treeVariance);
    Matrix treeP = treeV.inverse();
    sb.append("Tree variance:\n");
    sb.append(treeV);
    sb.append("Tree precision:\n");
    sb.append(treeP);
    sb.append("\n\n");
    sb.append("Trait variance:\n");
    sb.append(traitVariance);
    sb.append("\n\n");
    sb.append("Joint variance:\n");
    sb.append(new Matrix(jointVariance));
    sb.append("\n\n");
    double[] priorMean = rootPrior.getMean();
    sb.append("prior mean: ").append(new dr.math.matrixAlgebra.Vector(priorMean));
    sb.append("\n\n");
    sb.append("Joint variance:\n");
    sb.append(new Matrix(jointVariance));
    sb.append("\n\n");
    sb.append("Joint precision:\n");
    sb.append(new Matrix(getTreeTraitPrecision()));
    sb.append("\n\n");
    double[][] treeDrift = MultivariateTraitDebugUtilities.getTreeDrift(tree, priorMean, cdi, diffusionProcessDelegate);
    if (diffusionProcessDelegate.hasDrift()) {
        sb.append("Tree drift (including root mean):\n");
        sb.append(new Matrix(treeDrift));
        sb.append("\n\n");
    }
    double[] drift = KroneckerOperation.vectorize(treeDrift);
    final int datumLength = tipCount * dimProcess;
    sb.append("Tree dim : ").append(treeStructure.length).append("\n");
    sb.append("dimTrait : ").append(dimTrait).append("\n");
    sb.append("numTraits: ").append(numTraits).append("\n");
    sb.append("jVar dim : ").append(jointVariance.length).append("\n");
    sb.append("datum dim: ").append(datumLength);
    sb.append("\n\n");
    double[] data = dataModel.getParameter().getParameterValues();
    if (dataModel instanceof ContinuousTraitDataModel) {
        for (int tip = 0; tip < tipCount; ++tip) {
            double[] tipData = ((ContinuousTraitDataModel) dataModel).getTipObservation(tip, precisionType);
            System.arraycopy(tipData, 0, data, tip * numTraits * dimProcess, numTraits * dimProcess);
        }
    }
    sb.append("data: ").append(new dr.math.matrixAlgebra.Vector(data));
    sb.append("\n\n");
    double[][] graphStructure = MultivariateTraitDebugUtilities.getGraphVariance(tree, callbackLikelihood.getBranchRateModel(), normalization, priorSampleSize);
    double[][] jointGraphVariance = KroneckerOperation.product(graphStructure, traitVariance.toComponents());
    sb.append("graph structure:\n");
    sb.append(new Matrix(graphStructure));
    sb.append("\n\n");
    for (int trait = 0; trait < numTraits; ++trait) {
        sb.append("Trait #").append(trait).append("\n");
        double[] rawDatum = new double[datumLength];
        List<Integer> missing = new ArrayList<Integer>();
        int index = 0;
        for (int tip = 0; tip < tipCount; ++tip) {
            for (int dim = 0; dim < dimProcess; ++dim) {
                double d = data[tip * dimProcess * numTraits + trait * dimProcess + dim];
                rawDatum[index] = d;
                if (Double.isNaN(d)) {
                    missing.add(index);
                }
                ++index;
            }
        }
        double[][] varianceDatum = jointVariance;
        double[] datum = rawDatum;
        double[] driftDatum = drift;
        int[] notMissingIndices;
        notMissingIndices = new int[datumLength - missing.size()];
        int offsetNotMissing = 0;
        for (int i = 0; i < datumLength; ++i) {
            if (!missing.contains(i)) {
                notMissingIndices[offsetNotMissing] = i;
                ++offsetNotMissing;
            }
        }
        datum = Matrix.gatherEntries(rawDatum, notMissingIndices);
        varianceDatum = Matrix.gatherRowsAndColumns(jointVariance, notMissingIndices, notMissingIndices);
        driftDatum = Matrix.gatherEntries(drift, notMissingIndices);
        sb.append("datum : ").append(new dr.math.matrixAlgebra.Vector(datum)).append("\n");
        sb.append("drift : ").append(new dr.math.matrixAlgebra.Vector(driftDatum)).append("\n");
        sb.append("variance:\n");
        sb.append(new Matrix(varianceDatum));
        MultivariateNormalDistribution mvn = new MultivariateNormalDistribution(driftDatum, new Matrix(varianceDatum).inverse().toComponents());
        double logDensity = mvn.logPdf(datum);
        sb.append("\n\n");
        sb.append("logDatumLikelihood: ").append(logDensity).append("\n\n");
        // Compute joint for internal nodes
        int[] cNotMissingJoint = new int[dimProcess * tipCount];
        int[] cMissingJoint = new int[dimProcess * (tipCount - 1)];
        // External nodes
        for (int tipTrait = 0; tipTrait < dimProcess * tipCount; ++tipTrait) {
            cNotMissingJoint[tipTrait] = tipTrait;
        }
        // Internal nodes
        for (int tipTrait = dimProcess * tipCount; tipTrait < dimProcess * (2 * tipCount - 1); ++tipTrait) {
            cMissingJoint[tipTrait - dimProcess * tipCount] = tipTrait;
        }
        double[] rawDatumJoint = new double[dimProcess * (2 * tipCount - 1)];
        System.arraycopy(rawDatum, 0, rawDatumJoint, 0, rawDatum.length);
        double[][] driftJointMatrix = MultivariateTraitDebugUtilities.getGraphDrift(tree, cdi, diffusionProcessDelegate);
        double[] driftJoint = KroneckerOperation.vectorize(driftJointMatrix);
        for (int idx = 0; idx < driftJoint.length / dimProcess; ++idx) {
            for (int dim = 0; dim < dimProcess; ++dim) {
                driftJoint[idx * dimProcess + dim] += priorMean[dim];
            }
        }
        ConditionalVarianceAndTransform cVarianceJoint = new ConditionalVarianceAndTransform(new Matrix(jointGraphVariance), cMissingJoint, cNotMissingJoint);
        double[] cMeanJoint = cVarianceJoint.getConditionalMean(rawDatumJoint, 0, driftJoint, 0);
        sb.append("cDriftJoint: ").append(new dr.math.matrixAlgebra.Vector(driftJoint)).append("\n\n");
        sb.append("cMeanInternalJoint: ").append(new dr.math.matrixAlgebra.Vector(cMeanJoint)).append("\n\n");
        // Compute full conditional distributions
        sb.append("Full conditional distributions:\n");
        int offsetNotMissing2 = 0;
        for (int tip = 0; tip < tipCount; ++tip) {
            int offset = tip * dimProcess;
            int dimTip = 0;
            for (int cTrait = 0; cTrait < dimProcess; cTrait++) {
                if ((offsetNotMissing2 + cTrait < notMissingIndices.length) && notMissingIndices[offsetNotMissing2 + cTrait] < offset + dimProcess) {
                    dimTip++;
                }
            }
            int[] cMissing = new int[dimProcess];
            int[] cNotMissing = new int[notMissingIndices.length - dimTip];
            for (int cTrait = 0; cTrait < dimProcess; ++cTrait) {
                cMissing[cTrait] = offset + cTrait;
            }
            for (int m = 0; m < offsetNotMissing2; ++m) {
                cNotMissing[m] = notMissingIndices[m];
            }
            offsetNotMissing2 += dimTip;
            for (int m = offsetNotMissing2; m < notMissingIndices.length; ++m) {
                cNotMissing[m - dimTip] = notMissingIndices[m];
            }
            ConditionalVarianceAndTransform cVariance = new ConditionalVarianceAndTransform(new Matrix(jointVariance), cMissing, cNotMissing);
            double[] cMean = cVariance.getConditionalMean(rawDatum, 0, drift, 0);
            Matrix cVar = cVariance.getConditionalVariance();
            sb.append("cMean #").append(tip).append(" ").append(new dr.math.matrixAlgebra.Vector(cMean)).append(" cVar [").append(cVar).append("]\n");
        }
    }
    return sb.toString();
}

Example 14 with Matrix

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

the class IntegratedFactorAnalysisLikelihood method getReport.

@Override
public String getReport() {
    StringBuilder sb = new StringBuilder();
    double logComponents = 0;
    if (delegate != null) {
        double logInc = delegate.getCallbackLikelihood().getLogLikelihood();
        final Tree tree = delegate.getCallbackLikelihood().getTree();
        final BranchRates branchRates = delegate.getCallbackLikelihood().getBranchRateModel();
        sb.append(tree.toString());
        sb.append("\n\n");
        final double normalization = delegate.getRateTransformation().getNormalization();
        final double priorSampleSize = delegate.getRootProcessDelegate().getPseudoObservations();
        double[][] treeStructure = MultivariateTraitDebugUtilities.getTreeVariance(tree, branchRates, 1.0, Double.POSITIVE_INFINITY);
        sb.append("Tree structure:\n");
        sb.append(new Matrix(treeStructure));
        sb.append("\n\n");
        double[][] treeSharedLengths = MultivariateTraitDebugUtilities.getTreeVariance(tree, branchRates, normalization, Double.POSITIVE_INFINITY);
        double[][] treeVariance = MultivariateTraitDebugUtilities.getTreeVariance(tree, branchRates, normalization, priorSampleSize);
        Matrix treeV = new Matrix(treeVariance);
        Matrix treeP = treeV.inverse();
        sb.append("Tree variance:\n");
        sb.append(treeV);
        sb.append("Tree precision:\n");
        sb.append(treeP);
        sb.append("\n\n");
        Matrix Lt = new Matrix(loadingsTransposed.getParameterAsMatrix());
        sb.append("Loadings:\n");
        sb.append(Lt.transpose());
        sb.append("\n\n");
        double[][] diffusionPrecision = delegate.getDiffusionModel().getPrecisionmatrix();
        Matrix diffusionVariance = new Matrix(diffusionPrecision).inverse();
        Matrix loadingsVariance = null;
        try {
            loadingsVariance = Lt.product(diffusionVariance.product(Lt.transpose()));
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
        sb.append("Loadings variance:\n");
        sb.append(loadingsVariance);
        sb.append("\n\n");
        assert (loadingsVariance != null);
        Matrix loadingsFactorsVariance = MultivariateTraitDebugUtilities.getJointVarianceFactor(priorSampleSize, treeVariance, treeSharedLengths, loadingsVariance.toComponents(), diffusionVariance.toComponents(), delegate.getDiffusionProcessDelegate(), Lt);
        Matrix gamma = buildDiagonalMatrix(traitPrecision.getParameterValues());
        sb.append("Trait precision:\n");
        sb.append(gamma);
        sb.append("\n\n");
        Matrix gammaVariance = gamma.inverse();
        double[] tmp = new double[tree.getExternalNodeCount()];
        Arrays.fill(tmp, 1.0);
        Matrix identity = buildDiagonalMatrix(tmp);
        Matrix errorVariance = new Matrix(KroneckerOperation.product(identity.toComponents(), gammaVariance.toComponents()));
        sb.append("Loadings-factors variance:\n");
        sb.append(loadingsFactorsVariance);
        sb.append("\n\n");
        sb.append("Error variance\n");
        sb.append(errorVariance);
        sb.append("\n\n");
        Matrix totalVariance = null;
        try {
            totalVariance = loadingsFactorsVariance.add(errorVariance);
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
        double[] allData = getParameter().getParameterValues();
        List<Integer> notMissing = new ArrayList<>();
        for (int taxon = 0; taxon < numTaxa; ++taxon) {
            double[] observed = observedIndicators[taxon];
            for (int trait = 0; trait < dimTrait; ++trait) {
                if (observed[trait] == 0.0) {
                    System.err.println("Missing taxon " + taxon + " trait " + trait);
                } else {
                    notMissing.add(taxon * dimTrait + trait);
                }
            }
        }
        double[] priorMean = delegate.getRootPrior().getMean();
        Matrix treeDrift = new Matrix(MultivariateTraitDebugUtilities.getTreeDrift(tree, priorMean, delegate.getIntegrator(), delegate.getDiffusionProcessDelegate()));
        if (delegate.getDiffusionProcessDelegate().hasDrift()) {
            sb.append("Tree drift (including root mean):\n");
            sb.append(new Matrix(treeDrift.toComponents()));
            sb.append("\n\n");
        }
        try {
            loadingsFactorsVariance = treeDrift.product(Lt.transpose());
        } catch (IllegalDimension illegalDimension) {
            illegalDimension.printStackTrace();
        }
        double[] drift = KroneckerOperation.vectorize(loadingsFactorsVariance.toComponents());
        int[] notMissingIndices = new int[notMissing.size()];
        double[] data = new double[notMissing.size()];
        for (int i = 0; i < notMissing.size(); ++i) {
            notMissingIndices[i] = notMissing.get(i);
            data[i] = allData[notMissing.get(i)];
        }
        if (totalVariance != null) {
            totalVariance = new Matrix(Matrix.gatherRowsAndColumns(totalVariance.toComponents(), notMissingIndices, notMissingIndices));
        }
        Matrix totalPrecision = null;
        if (totalVariance != null) {
            totalPrecision = totalVariance.inverse();
        }
        drift = Matrix.gatherEntries(drift, notMissingIndices);
        sb.append("Total variance:\n");
        sb.append(totalVariance);
        sb.append("\n\n");
        sb.append("Total precision:\n");
        sb.append(totalPrecision);
        sb.append("\n\n");
        sb.append("Data:\n");
        sb.append(new Vector(data));
        sb.append("\n\n");
        sb.append("Expectations:\n");
        sb.append(new Vector(drift));
        sb.append("\n\n");
        MultivariateNormalDistribution mvn = null;
        if (totalPrecision != null) {
            mvn = new MultivariateNormalDistribution(drift, totalPrecision.toComponents());
        }
        double logDensity = 0;
        if (mvn != null) {
            logDensity = mvn.logPdf(data);
        }
        sb.append("logMultiVariateNormalDensity = ").append(logDensity).append("\n\n");
        sb.append("traitDataLikelihood = ").append(logInc).append("\n");
        logComponents += logInc;
    }
    sb.append("logLikelihood = ").append(getLogLikelihood()).append("\n");
    if (logComponents != 0.0) {
        sb.append("total likelihood = ").append((getLogLikelihood() + logComponents)).append("\n");
    }
    return sb.toString();
}

Example 15 with Matrix

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

the class RepeatedMeasuresTraitDataModel method restoreState.

@Override
protected void restoreState() {
    Matrix tmp = samplingPrecision;
    samplingPrecision = storedSamplingPrecision;
    storedSamplingPrecision = tmp;
    tmp = samplingVariance;
    samplingVariance = storedSamplingVariance;
    storedSamplingVariance = tmp;
    varianceKnown = storedVarianceKnown;
    variableChanged = storedVariableChanged;
}