Java Code Examples for org.ejml.data.DenseMatrix64F#unsafe_set()

public static void scatterRowsAndColumns(final DenseMatrix64F source, final DenseMatrix64F destination,
                                         final int[] rowIdices, final int[] colIndices, final boolean clear) {
    if (clear) {
        Arrays.fill(destination.getData(), 0.0);
    }

    final int rowLength = rowIdices.length;
    final int colLength = colIndices.length;
    final double[] in = source.getData();

    int index = 0;
    for (int i = 0; i < rowLength; ++i) {
        final int rowIndex = rowIdices[i];
        for (int j = 0; j < colLength; ++j) {
            destination.unsafe_set(rowIndex, colIndices[j], in[index]);
            ++index;
        }
    }
}
 

private DenseMatrix64F getGradientBranchVarianceWrtAttenuationDiagonal(ContinuousDiffusionIntegrator cdi,
                                                                    int nodeIndex, DenseMatrix64F gradient) {

    double[] attenuation = elasticModel.getEigenValuesStrengthOfSelection();
    DenseMatrix64F variance = wrap(
            ((MultivariateIntegrator) cdi).getVariance(getEigenBufferOffsetIndex(0)),
            0, dim, dim);

    double ti = cdi.getBranchLength(getMatrixBufferOffsetIndex(nodeIndex));

    DenseMatrix64F res = new DenseMatrix64F(dim, 1);

    CommonOps.elementMult(variance, gradient);

    for (int k = 0; k < dim; k++) {
        double sum = 0.0;
        for (int l = 0; l < dim; l++) {
            sum -= variance.unsafe_get(k, l) * computeAttenuationFactorActualized(attenuation[k] + attenuation[l], ti);
        }
        res.unsafe_set(k, 0, sum);
    }

    return res;
}
 

public static void copyRowsAndColumns(final DenseMatrix64F source, final DenseMatrix64F destination,
                                      final int[] rowIndices, final int[] colIndices, final boolean clear) {
    if (clear) {
        Arrays.fill(destination.getData(), 0.0);
    }
    for (int row : rowIndices) {
        for (int col : colIndices) {
            destination.unsafe_set(row, col, source.unsafe_get(row, col));
        }
    }
}
 

public void makeDeltaVector(DenseMatrix64F delta, NormalSufficientStatistics jointStatistics,
                            NormalSufficientStatistics above) {
    for (int row = 0; row < dim; ++row) {
        delta.unsafe_set(row, 0,
                jointStatistics.getRawMean().unsafe_get(row, 0) - above.getMean(row)
        );
    }
}
 

private void makeCompletedUnobserved(final DenseMatrix64F matrix, double diagonal) {
    for (int row = 0; row < numFactors; ++row) {
        for (int col = 0; col < numFactors; ++col) {
            double x = (row == col) ? diagonal : 0.0;
            matrix.unsafe_set(row, col, x);
        }
    }
}
 

private void actualizeGradientDiagonal(ContinuousDiffusionIntegrator cdi, int nodeIndex, DenseMatrix64F gradient) {
    double[] attenuation = elasticModel.getEigenValuesStrengthOfSelection();
    double edgeLength = cdi.getBranchLength(getMatrixBufferOffsetIndex(nodeIndex));
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {
            gradient.unsafe_set(i, j,
                    factorFunction(attenuation[i] + attenuation[j], edgeLength) * gradient.unsafe_get(i, j));
        }
    }
}
 

private static void removeMissing(DenseMatrix64F M, int[] missing) {
    for (int m : missing) {
        for (int j = 0; j < M.getNumCols(); j++) {
            M.unsafe_set(m, j, 0.0);
            M.unsafe_set(j, m, 0.0);
        }
    }
}
 

public static void addToDiagonal(DenseMatrix64F source, double increment) {
    final int width = source.getNumRows();
    for (int i = 0; i < width; ++i) {
        source.unsafe_set(i, i, source.unsafe_get(i, i) + increment);
    }
}
 

private static NormalSufficientStatistics computeJointPartiallyMissing(NormalSufficientStatistics below,
                                                                       NormalSufficientStatistics above,
                                                                       PermutationIndices indices,
                                                                       int dim) {

    DenseMatrix64F mean = new DenseMatrix64F(dim, 1);
    DenseMatrix64F precision = new DenseMatrix64F(dim, dim);
    DenseMatrix64F variance = new DenseMatrix64F(dim, dim);

    if (indices.getNumberOfNonZeroFiniteDiagonals() != 0) {
        throw new RuntimeException("Unsure if this works for latent trait below");
        // TODO Probably need to add in below.precision somewhere below
    }

    ConditionalPrecisionAndTransform2 transform = new ConditionalPrecisionAndTransform2(
            above.getRawPrecision(),
            indices.getZeroIndices(),
            indices.getInfiniteIndices()
    );

    double[] result = transform.getConditionalMean(below.getRawMean().getData(), 0,
            above.getRawMean().getData(), 0);

    copyRowsAndColumns(above.getRawPrecision(), precision,
            indices.getZeroIndices(), indices.getZeroIndices(), false);

    scatterRowsAndColumns(transform.getConditionalVariance(), variance,
            indices.getZeroIndices(), indices.getZeroIndices(), false);

    int index = 0;
    for (int zero : indices.getZeroIndices()) {
        mean.unsafe_set(zero, 0, result[index++]);
    }

    for (int infinite : indices.getInfiniteIndices()) {
        mean.unsafe_set(infinite, 0, below.getMean(infinite));
        precision.unsafe_set(infinite, infinite, Double.POSITIVE_INFINITY);
    }

    return new NormalSufficientStatistics(mean, precision, variance);
}
 

private void computePrecisionForTaxon(final DenseMatrix64F precision, final int taxon,
                                      final int numFactors) {

    final double[] observed = observedIndicators[taxon];

    final HashedMissingArray observedArray;
    DenseMatrix64F hashedPrecision;

    if (USE_PRECISION_CACHE) {
        observedArray = new HashedMissingArray(observed);
        hashedPrecision = precisionMatrixMap.get(observedArray);
    }

    // TODO Only need to compute for each unique set of observed[] << numTaxa

    if (!USE_PRECISION_CACHE || hashedPrecision == null) {

        // Compute L D_i \Gamma D_i^t L^t
        for (int row = 0; row < numFactors; ++row) {
            for (int col = row; col < numFactors; ++col) {
                double sum = 0;
                for (int k = 0; k < dimTrait; ++k) {
                    double thisPrecision = (observed[k] == 1.0) ?
                            gamma[k] // traitPrecision.getParameterValue(k)
                            : nuggetPrecision;
                    sum += loadings[row * dimTrait + k] * //loadingsTransposed.getParameterValue(k, row) *
                            thisPrecision *
                            loadings[col * dimTrait + k]; // loadingsTransposed.getParameterValue(k, col);
                }
                precision.unsafe_set(row, col, sum);
                precision.unsafe_set(col, row, sum); // Symmetric matrix
            }
        }

        if (USE_PRECISION_CACHE) {
            precisionMatrixMap.put(observedArray, precision);
        }

    } else {
        System.arraycopy(hashedPrecision.getData(), 0,
                precision.getData(), 0, numFactors * numFactors);
    }
}
 

void getSufficientStatistics(BranchSufficientStatistics statistics) {
    final NormalSufficientStatistics below = statistics.getBelow();
    final NormalSufficientStatistics above = statistics.getAbove();

    // Sampling Parameters: Statistic data model
    DenseMatrix64F samplingVariance = dataModel.getExtensionVariance();

    // One more pre-order step
    // TODO This is just one more pre-order step. Should maybe be moved elsewhere ?
    // TODO Below only works for one data point (no repetition)
    // above data
    DenseMatrix64F aboveDataVariance = new DenseMatrix64F(dim, dim);
    DenseMatrix64F aboveDataPrecision = new DenseMatrix64F(dim, dim);
    CommonOps.add(above.getRawVariance(), samplingVariance, aboveDataVariance);
    safeInvert2(aboveDataVariance, aboveDataPrecision, false);
    final NormalSufficientStatistics aboveData = new NormalSufficientStatistics(
            above.getRawMeanCopy(),
            aboveDataPrecision,
            aboveDataVariance);
    // Below data
    int[] missings = statistics.getMissing();
    DenseMatrix64F precisionData = new DenseMatrix64F(dim, dim);
    for (int i = 0; i < dim; i++) {
        precisionData.unsafe_set(i, i, Double.POSITIVE_INFINITY);
    }
    for (int m : missings) {
        precisionData.unsafe_set(m, m, 0.0);
    }
    final NormalSufficientStatistics belowData = new NormalSufficientStatistics(
            below.getRawMeanCopy(),
            precisionData);

    // Joint data
    NormalSufficientStatistics jointStatisticsData =
            BranchRateGradient.ContinuousTraitGradientForBranch.Default.computeJointStatistics(
                    belowData, aboveData, dim
            );

    // Gradient
    matrixQ = aboveData.getRawPrecision();
    matrixW = aboveData.getRawVariance();
    matrixV = jointStatisticsData.getRawVariance();

    // Delta
    for (int row = 0; row < dim; ++row) {
        matrixDelta.unsafe_set(row, 0,
                jointStatisticsData.getMean(row) - aboveData.getMean(row)
        );
    }
}