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

Example 1 with WrappedMatrix

use of dr.math.matrixAlgebra.WrappedMatrix 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 2 with WrappedMatrix

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

the class DifferentiableSubstitutionModelUtil method getDifferentialMassMatrix.

public static double[] getDifferentialMassMatrix(double time, int stateCount, WrappedMatrix differentialMassMatrix, EigenDecomposition eigenDecomposition) {
    double[] eigenValues = eigenDecomposition.getEigenValues();
    WrappedMatrix eigenVectors = new WrappedMatrix.Raw(eigenDecomposition.getEigenVectors(), 0, stateCount, stateCount);
    WrappedMatrix inverseEigenVectors = new WrappedMatrix.Raw(eigenDecomposition.getInverseEigenVectors(), 0, stateCount, stateCount);
    getTripleMatrixMultiplication(stateCount, inverseEigenVectors, differentialMassMatrix, eigenVectors);
    for (int i = 0; i < stateCount; i++) {
        for (int j = 0; j < stateCount; j++) {
            if (i == j || eigenValues[i] == eigenValues[j]) {
                differentialMassMatrix.set(i, j, differentialMassMatrix.get(i, j) * time);
            } else {
                differentialMassMatrix.set(i, j, differentialMassMatrix.get(i, j) * (1.0 - Math.exp((eigenValues[j] - eigenValues[i]) * time)) / (eigenValues[i] - eigenValues[j]));
            }
        }
    }
    getTripleMatrixMultiplication(stateCount, eigenVectors, differentialMassMatrix, inverseEigenVectors);
    double[] outputArray = new double[stateCount * stateCount];
    for (int i = 0, length = stateCount * stateCount; i < length; ++i) {
        outputArray[i] = differentialMassMatrix.get(i);
    }
    return outputArray;
}

Example 3 with WrappedMatrix

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

the class MultivariateConditionalOnTipsRealizedDelegate method simulateTraitForInternalNode.

// private ReadableVector getMeanWithDrift(final ReadableVector mean,
// final ReadableVector drift) {
// return new ReadableVector.Sum(mean, drift);
// }
// private ReadableVector getMeanWithDrift(double[] mean, int offsetMean, double[] drift, int dim) {
// for (int i = 0;i < dim; ++i) {
// tmpDrift[i] = mean[offsetMean + i] + drift[i];
// }
// return new WrappedVector.Raw(tmpDrift, 0, dimTrait);
// }
private void simulateTraitForInternalNode(final int offsetSample, final int offsetParent, final int offsetPartial, final double branchPrecision) {
    if (!Double.isInfinite(branchPrecision)) {
        // Here we simulate X_j | X_pa(j), Y
        final WrappedVector M0 = new WrappedVector.Raw(partialNodeBuffer, offsetPartial, dimTrait);
        final DenseMatrix64F P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
        // final ReadableVector parentSample = new WrappedVector.Raw(sample, offsetParent, dimTrait);
        final ReadableVector M1;
        final DenseMatrix64F P1;
        M1 = getMeanBranch(offsetParent);
        P1 = getPrecisionBranch(branchPrecision);
        // if (hasNoDrift) {
        // M1 = parentSample; // new WrappedVector.Raw(sample, offsetParent, dimTrait);
        // P1 = new DenseMatrix64F(dimTrait, dimTrait);
        // CommonOps.scale(branchPrecision, Pd, P1);
        // } else {
        // M1 = getMeanWithDrift(parentSample,
        // new WrappedVector.Raw(displacementBuffer, 0, dimTrait)); //getMeanWithDrift(sample, offsetParent, displacementBuffer, dimTrait);
        // P1 = DenseMatrix64F.wrap(dimTrait, dimTrait, precisionBuffer);
        // }
        final WrappedVector M2 = new WrappedVector.Raw(tmpMean, 0, dimTrait);
        final DenseMatrix64F P2 = new DenseMatrix64F(dimTrait, dimTrait);
        final DenseMatrix64F V2 = new DenseMatrix64F(dimTrait, dimTrait);
        CommonOps.add(P0, P1, P2);
        safeInvert2(P2, V2, false);
        weightedAverage(M0, P0, M1, P1, M2, V2, dimTrait);
        final WrappedMatrix C2;
        if (NEW_CHOLESKY) {
            DenseMatrix64F tC2 = getCholeskyOfVariance(V2, dimTrait);
            C2 = new WrappedMatrix.WrappedDenseMatrix(tC2);
            MultivariateNormalDistribution.nextMultivariateNormalCholesky(// input mean
            M2, // input variance
            C2, // input variance
            1.0, new WrappedVector.Raw(sample, offsetSample, dimTrait), tmpEpsilon);
        } else {
            double[][] tC2 = getCholeskyOfVariance(V2.getData(), dimTrait);
            C2 = new WrappedMatrix.ArrayOfArray(tC2);
            MultivariateNormalDistribution.nextMultivariateNormalCholesky(// input mean
            M2, // input variance
            C2, // input variance
            1.0, new WrappedVector.Raw(sample, offsetSample, dimTrait), tmpEpsilon);
        }
        if (DEBUG) {
            System.err.println("sT F I N");
            System.err.println("M0: " + M0);
            System.err.println("P0: " + P0);
            System.err.println("M1: " + M1);
            System.err.println("P1: " + P1);
            System.err.println("M2: " + M2);
            System.err.println("V2: " + V2);
            System.err.println("C2: " + C2);
            System.err.println("SS: " + new WrappedVector.Raw(sample, offsetSample, dimTrait));
            System.err.println("");
            if (!check(M2)) {
                System.exit(-1);
            }
        }
    } else {
        System.arraycopy(sample, offsetParent, sample, offsetSample, dimTrait);
        if (DEBUG) {
            System.err.println("sT F I N infinite branch precision");
            System.err.println("SS: " + new WrappedVector.Raw(sample, offsetSample, dimTrait));
        }
    }
}

Example 4 with WrappedMatrix

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

the class CachedMatrixInverse method restoreValues.

@Override
protected void restoreValues() {
    super.restoreValues();
    inverseKnown = savedInverseKnown;
    WrappedMatrix tmp = inverse;
    inverse = savedInverse;
    savedInverse = tmp;
}