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

The following examples show how to use org.ejml.data.DenseMatrix64F#unsafe_get() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.

Example 1

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

6 votes

public static void weightedSumActualized(final double[] ipartial,
                                         final int ibo,
                                         final DenseMatrix64F Pi,
                                         final double[] iactualization,
                                         final int ido,
                                         final double[] jpartial,
                                         final int jbo,
                                         final DenseMatrix64F Pj,
                                         final double[] jactualization,
                                         final int jdo,
                                         final int dimTrait,
                                         final double[] out) {
    for (int g = 0; g < dimTrait; ++g) {
        double sum = 0.0;
        for (int h = 0; h < dimTrait; ++h) {
            sum += iactualization[ido + g] * Pi.unsafe_get(g, h) * ipartial[ibo + h];
            sum += jactualization[jdo + g] * Pj.unsafe_get(g, h) * jpartial[jbo + h];
        }
        out[g] = sum;
    }
}

Example 2

Source File: OUDiffusionModelDelegate.java From beast-mcmc with GNU Lesser General Public License v2.1

6 votes

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;
}

Example 3

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

6 votes

public static double weightedInnerProductOfDifferences(final double[] source1,
                                                       final int source1Offset,
                                                       final double[] source2,
                                                       final int source2Offset,
                                                       final DenseMatrix64F P,
                                                       final int dimTrait) {
    double SS = 0;
    for (int g = 0; g < dimTrait; ++g) {
        final double gDifference = source1[source1Offset + g] - source2[source2Offset + g];

        for (int h = 0; h < dimTrait; ++h) {
            final double hDifference = source1[source1Offset + h] - source2[source2Offset + h];

            SS += gDifference * P.unsafe_get(g, h) * hDifference;
        }
    }

    return SS;
}

Example 4

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

6 votes

public static double weightedInnerProduct(final double[] partials,
                                          final int bo,
                                          final DenseMatrix64F P,
                                          final int dimTrait) {
    double SS = 0;

    // vector-matrix-vector
    for (int g = 0; g < dimTrait; ++g) {
        final double ig = partials[bo + g];
        for (int h = 0; h < dimTrait; ++h) {
            final double ih = partials[bo + h];
            SS += ig * P.unsafe_get(g, h) * ih;
        }
    }

    return SS;
}

Example 5

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

6 votes

public static void gatherRowsAndColumns(final DenseMatrix64F source, final DenseMatrix64F destination,
                                        final int[] rowIndices, final int[] colIndices) {

    final int rowLength = rowIndices.length;
    final int colLength = colIndices.length;
    final double[] out = destination.getData();

    int index = 0;
    for (int i = 0; i < rowLength; ++i) {
        final int rowIndex = rowIndices[i];
        for (int j = 0; j < colLength; ++j) {
            out[index] = source.unsafe_get(rowIndex, colIndices[j]);
            ++index;
        }
    }
}

Example 6

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

6 votes

public static void weightedSum(final double[] ipartial,
                               final int ibo,
                               final DenseMatrix64F Pi,
                               final double[] jpartial,
                               final int jbo,
                               final DenseMatrix64F Pj,
                               final int dimTrait,
                               final double[] out) {
    for (int g = 0; g < dimTrait; ++g) {
        double sum = 0.0;
        for (int h = 0; h < dimTrait; ++h) {
            sum += Pi.unsafe_get(g, h) * ipartial[ibo + h];
            sum += Pj.unsafe_get(g, h) * jpartial[jbo + h];
        }
        out[g] = sum;
    }
}

Example 7

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

6 votes

public static void matrixVectorMultiple(final DenseMatrix64F A,
                                        final WrappedVector x,
                                        final WrappedVector y,
                                        final int dim) {
    if (buffer.length < dim) {
        buffer = new double[dim];
    }

    for (int row = 0; row < dim; ++row) {
        double sum = 0.0;
        for (int col = 0; col < dim; ++col) {
            sum += A.unsafe_get(row, col) * x.get(col);
        }
        buffer[row] = sum;
    }

    for (int col = 0; col < dim; ++col) {
        y.set(col, buffer[col]);
    }
}

Example 8

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static void getFiniteNonZeroDiagonalIndices(final DenseMatrix64F source, final int[] indices) {
    final int length = source.getNumCols();

    int index = 0;
    for (int i = 0; i < length; ++i) {
        final double d = source.unsafe_get(i, i);
        if (!Double.isInfinite(d) && d != 0.0) {
            indices[index] = i;
            ++index;
        }
    }
}

Example 9

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static int countFiniteNonZeroDiagonals(DenseMatrix64F source) {
    final int length = source.getNumCols();

    int count = 0;
    for (int i = 0; i < length; ++i) {
        final double d = source.unsafe_get(i, i);
        if (!Double.isInfinite(d) && d != 0.0) {
            ++count;
        }
    }
    return count;
}

Example 10

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static void getFiniteDiagonalIndices(final DenseMatrix64F source, final int[] indices) {
    final int length = source.getNumCols();

    int index = 0;
    for (int i = 0; i < length; ++i) {
        final double d = source.unsafe_get(i, i);
        if (!Double.isInfinite(d)) {
            indices[index] = i;
            ++index;
        }
    }
}

Example 11

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static void weightedAverage(final double[] ipartial,
                                   final int ibo,
                                   final DenseMatrix64F Pi,
                                   final double[] jpartial,
                                   final int jbo,
                                   final DenseMatrix64F Pj,
                                   final double[] kpartial,
                                   final int kbo,
                                   final DenseMatrix64F Vk,
                                   final int dimTrait,
                                   final double[] tmp) {
    weightedSum(ipartial, ibo, Pi, jpartial, jbo, Pj, dimTrait, tmp);
    for (int g = 0; g < dimTrait; ++g) {
        if (!Double.isInfinite(Vk.unsafe_get(g, g))) {
            double sum = 0.0;
            for (int h = 0; h < dimTrait; ++h) {
                if (!Double.isInfinite(Vk.unsafe_get(h, h))) {

                    sum += Vk.unsafe_get(g, h) * tmp[h];

                }
            }
            kpartial[kbo + g] = sum;

        }
    }
}

Example 12

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static boolean allZeroDiagonals(DenseMatrix64F source) {
    final int length = source.getNumCols();

    for (int i = 0; i < length; ++i) {
        if (source.unsafe_get(i, i) != 0.0) {
            return false;
        }
    }
    return true;
}

Example 13

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static int countZeroDiagonals(DenseMatrix64F source) {
    final int length = source.getNumCols();

    int count = 0;
    for (int i = 0; i < length; ++i) {
        final double d = source.unsafe_get(i, i);
        if (d == 0.0) {
            ++count;
        }
    }
    return count;
}

Example 14

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static double weightedThreeInnerProduct(final double[] ipartials,
                                               final int ibo,
                                               final DenseMatrix64F Pip,
                                               final double[] jpartials,
                                               final int jbo,
                                               final DenseMatrix64F Pjp,
                                               final double[] kpartials,
                                               final int kbo,
                                               final DenseMatrix64F Pk,
                                               final int dimTrait) {

    // TODO Is it better to split into 3 separate calls to weightedInnerProduct?

    double SSi = 0;
    double SSj = 0;
    double SSk = 0;


    // vector-matrix-vector TODO in parallel
    for (int g = 0; g < dimTrait; ++g) {
        final double ig = ipartials[ibo + g];
        final double jg = jpartials[jbo + g];
        final double kg = kpartials[kbo + g];

        for (int h = 0; h < dimTrait; ++h) {
            final double ih = ipartials[ibo + h];
            final double jh = jpartials[jbo + h];
            final double kh = kpartials[kbo + h];


            SSi += ig * Pip.unsafe_get(g, h) * ih;
            SSj += jg * Pjp.unsafe_get(g, h) * jh;
            SSk += kg * Pk.unsafe_get(g, h) * kh;

        }
    }

    return SSi + SSj - SSk;
}

Example 15

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static double weightedThreeInnerProductNormalized(final double[] ipartials,
                                                         final int ibo,
                                                         final DenseMatrix64F Pip,
                                                         final double[] jpartials,
                                                         final int jbo,
                                                         final DenseMatrix64F Pjp,
                                                         final double[] kpartials,
                                                         final int kbo,
                                                         final double[] kpartialsBis,
                                                         final int kboBis,
                                                         final int dimTrait) {

    double SSi = 0;
    double SSj = 0;
    double SSk = 0;

    // vector-matrix-vector TODO in parallel
    for (int g = 0; g < dimTrait; ++g) {
        final double ig = ipartials[ibo + g];
        final double jg = jpartials[jbo + g];
        final double kg = kpartials[kbo + g];
        final double kgBis = kpartialsBis[kboBis + g];

        for (int h = 0; h < dimTrait; ++h) {
            final double ih = ipartials[ibo + h];
            final double jh = jpartials[jbo + h];

            SSi += ig * Pip.unsafe_get(g, h) * ih;
            SSj += jg * Pjp.unsafe_get(g, h) * jh;
        }
        SSk += kg * kgBis;
    }

    return SSi + SSj - SSk;
}

Example 16

Source File: IntegratedFactorAnalysisLikelihood.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

private double computeFactorInnerProduct(final WrappedVector mean, final DenseMatrix64F precision) {
    // Compute \mu_i^t P_i \mu^t
    double sum = 0;
    for (int row = 0; row < numFactors; ++row) {
        for (int col = 0; col < numFactors; ++col) {
            sum += mean.get(row) * precision.unsafe_get(row, col) * mean.get(col);
        }
    }
    return sum;
}

Example 17

Source File: MissingOps.java From beast-mcmc with GNU Lesser General Public License v2.1

5 votes

public static int countFiniteDiagonals(DenseMatrix64F source) {
    final int length = source.getNumCols();

    int count = 0;
    for (int i = 0; i < length; ++i) {
        final double d = source.unsafe_get(i, i);
        if (!Double.isInfinite(d)) {
            ++count;
        }
    }
    return count;
}

Example 18

Source File: BranchRateGradient.java From beast-mcmc with GNU Lesser General Public License v2.1

4 votes

@Override
            public double getGradientForBranch(BranchSufficientStatistics statistics, double differentialScaling) {

                final NormalSufficientStatistics below = statistics.getBelow();
                final NormalSufficientStatistics branch = statistics.getBranch();
                final NormalSufficientStatistics above = statistics.getAbove();

                if (DEBUG) {
                    System.err.println("B = " + statistics.toVectorizedString());
                }

//                if (DEBUG) {
//                    System.err.println("\tQi = " + NormalSufficientStatistics.toVectorizedString(Qi));
//                    System.err.println("\tSi = " + NormalSufficientStatistics.toVectorizedString(Si));
//                }

                DenseMatrix64F Qi = above.getRawPrecision();

                DenseMatrix64F logDetComponent = matrix0;
                DenseMatrix64F quadraticComponent = matrix1;
                makeGradientMatrices0(quadraticComponent, logDetComponent, statistics, differentialScaling);

                double grad1 = 0.0;
                for (int row = 0; row < dim; ++row) {
                    grad1 -= 0.5 * logDetComponent.unsafe_get(row, row);
                }

                if (DEBUG) {
                    System.err.println("grad1 = " + grad1);
                }


//                if (DEBUG) {
//                    System.err.println("\tQi = " + NormalSufficientStatistics.toVectorizedString(Qi));
//                    System.err.println("\tQQ = " + NormalSufficientStatistics.toVectorizedString(quadraticComponent));
//                }

                NormalSufficientStatistics jointStatistics = computeJointStatistics(below, above, dim);

                DenseMatrix64F additionalVariance = matrix0; //new DenseMatrix64F(dim, dim);
                makeGradientMatrices1(additionalVariance, quadraticComponent,
                        jointStatistics);


                DenseMatrix64F delta = vector0;
                makeDeltaVector(delta, jointStatistics, above);

                double grad2 = 0.0;
                for (int row = 0; row < dim; ++row) {
                    for (int col = 0; col < dim; ++col) {

                        grad2 += 0.5 * delta.unsafe_get(row, 0)
                                * quadraticComponent.unsafe_get(row, col)
                                * delta.unsafe_get(col, 0);

                    }

                    grad2 += 0.5 * additionalVariance.unsafe_get(row, row);
                }

                if (DEBUG) {
                    System.err.println("\tjoint mean  = " + NormalSufficientStatistics.toVectorizedString(jointStatistics.getRawMean()));
                    System.err.println("\tabove mean = " + NormalSufficientStatistics.toVectorizedString(above.getRawMean()));
                    System.err.println("\tbelow mean  = " + NormalSufficientStatistics.toVectorizedString(below.getRawMean()));
                    System.err.println("\tjoint precision  = " + NormalSufficientStatistics.toVectorizedString(jointStatistics.getRawPrecision()));
                    System.err.println("\tabove precision = " + NormalSufficientStatistics.toVectorizedString(above.getRawPrecision()));
                    System.err.println("\tbelow precision  = " + NormalSufficientStatistics.toVectorizedString(below.getRawPrecision()));
                    System.err.println("\tbelow variance   = " + NormalSufficientStatistics.toVectorizedString(below.getRawVariance()));
                    System.err.println("\tquadratic      = " + NormalSufficientStatistics.toVectorizedString(quadraticComponent));
                    System.err.println("\tadditional     = " + NormalSufficientStatistics.toVectorizedString(additionalVariance));
                    System.err.println("delta: " + NormalSufficientStatistics.toVectorizedString(delta));
                    System.err.println("grad2 = " + grad2);
                }

                // W.r.t. drift
                // TODO: Fix delegate to (possibly) un-link drift from arbitrary rate
                DenseMatrix64F Di = new DenseMatrix64F(dim, 1);
                CommonOps.scale(differentialScaling, branch.getRawMean(), Di);

                double grad3 = 0.0;
                for (int row = 0; row < dim; ++row) {
                    for (int col = 0; col < dim; ++col) {

                        grad3 += delta.unsafe_get(row, 0)
                                * Qi.unsafe_get(row, col)
                                * Di.unsafe_get(col, 0);

                    }
                }

                if (DEBUG) {
                    System.err.println("\tDi     = " + NormalSufficientStatistics.toVectorizedString(branch.getRawMean()));
                    System.err.println("grad3 = " + grad3);
                }

                return grad1 + grad2 + grad3;

            }