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

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

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

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

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

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

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

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

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

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

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;

        }
    }
}
 

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

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

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

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

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

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

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

            }