Examples with DummyLatentTruncationProvider dr.evomodel.continuous.DummyLatentTruncationProvider used on opensource projects

Example 1 with DummyLatentTruncationProvider

use of dr.evomodel.continuous.DummyLatentTruncationProvider in project beast-mcmc by beast-dev.

the class NewLatentLiabilityGibbs method sampleNode.

private double sampleNode(NodeRef node, WrappedNormalSufficientStatistics statistics) {
    final int thisNumber = node.getNumber();
    final int[] obsInds = maskDelegate.getObservedIndices(node);
    final int obsDim = obsInds.length;
    if (obsDim == dim) {
        return 0;
    }
    ReadableVector mean = statistics.getMean();
    ReadableMatrix thisP = statistics.getPrecision();
    double precisionScalar = statistics.getPrecisionScalar();
    for (int i = 0; i < mean.getDim(); ++i) {
        fcdMean[i] = mean.get(i);
    }
    for (int i = 0; i < mean.getDim(); ++i) {
        for (int j = 0; j < mean.getDim(); ++j) {
            fcdPrecision[i][j] = thisP.get(i, j) * precisionScalar;
        }
    }
    if (repeatedMeasuresModel != null) {
        // TODO: preallocate memory for all these matrices/vectors
        // storing original fcd precision
        DenseMatrix64F Q = new DenseMatrix64F(fcdPrecision);
        double[] tipPartial = repeatedMeasuresModel.getTipPartial(thisNumber, false);
        int offset = dim;
        DenseMatrix64F P = MissingOps.wrap(tipPartial, offset, dim, dim);
        DenseMatrix64F preOrderMean = new DenseMatrix64F(dim, 1);
        for (int i = 0; i < dim; i++) {
            preOrderMean.set(i, 0, fcdMean[i]);
        }
        DenseMatrix64F dataMean = new DenseMatrix64F(dim, 1);
        for (int i = 0; i < dim; i++) {
            dataMean.set(i, 0, tipPartial[i]);
        }
        D1Matrix64F bufferMean = new DenseMatrix64F(dim, 1);
        // bufferMean = Q * preOrderMean
        mult(Q, preOrderMean, bufferMean);
        // bufferMean = Q * preOderMean + P * dataMean
        multAdd(P, dataMean, bufferMean);
        // P = P + Q
        CommonOps.addEquals(P, Q);
        DenseMatrix64F V = new DenseMatrix64F(dim, dim);
        // V = inv(P + Q)
        CommonOps.invert(P, V);
        // dataMean = inv(P + Q) * (Q * preOderMean + P * dataMean)
        mult(V, bufferMean, dataMean);
        for (int i = 0; i < dim; i++) {
            fcdMean[i] = dataMean.get(i);
            for (int j = 0; j < dim; j++) {
                fcdPrecision[i][j] = P.get(i, j);
            }
        }
    }
    // TODO: should this not be declared until 'else' statement?
    MultivariateNormalDistribution fullDistribution = new MultivariateNormalDistribution(fcdMean, fcdPrecision);
    MultivariateNormalDistribution drawDistribution;
    if (mask != null && obsDim > 0) {
        addMaskOnContiuousTraitsPrecisionSpace(thisNumber);
        drawDistribution = new MultivariateNormalDistribution(maskedMean, maskedPrecision);
    } else {
        drawDistribution = fullDistribution;
    }
    double[] oldValue = getNodeTrait(node);
    int attempt = 0;
    boolean validTip = false;
    while (!validTip & attempt < maxAttempts) {
        setNodeTrait(node, drawDistribution.nextMultivariateNormal());
        if (latentLiability.validTraitForTip(thisNumber)) {
            validTip = true;
        }
        attempt++;
    }
    if (attempt == maxAttempts) {
        return Double.NEGATIVE_INFINITY;
    }
    double[] newValue = getNodeTrait(node);
    if (latentLiability instanceof DummyLatentTruncationProvider) {
        return Double.POSITIVE_INFINITY;
    } else {
        return fullDistribution.logPdf(oldValue) - fullDistribution.logPdf(newValue);
    }
}