Search in sources :

Example 1 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class ErfGaussian2DFunctionTest method functionComputesExtendedGradientForEach.

@Test
public void functionComputesExtendedGradientForEach() {
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) this.f1;
    final int nparams = f1.getNumberOfGradients();
    final int[] gradientIndices = f1.gradientIndices();
    final double[] du_da = new double[f1.getNumberOfGradients()];
    final double[] du_db = new double[f1.getNumberOfGradients()];
    final ErfGaussian2DFunction[] fHigh = new ErfGaussian2DFunction[nparams];
    final ErfGaussian2DFunction[] fLow = new ErfGaussian2DFunction[nparams];
    final double[] delta = new double[nparams];
    for (int j = 0; j < nparams; j++) {
        fHigh[j] = f1.copy();
        fLow[j] = f1.copy();
    }
    for (double background : testbackground) // Peak 1
    for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
        double[] a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        f1.initialiseExtended2(a);
        // Create a set of functions initialised +/- delta in each parameter
        for (int j = 0; j < nparams; j++) {
            int targetParameter = gradientIndices[j];
            // Numerically solve gradient. 
            // Calculate the step size h to be an exact numerical representation
            final double xx = a[targetParameter];
            // Get h to minimise roundoff error
            double h = Precision.representableDelta(xx, h_);
            // Evaluate at (x+h) and (x-h)
            a[targetParameter] = xx + h;
            fHigh[j].initialise1(a);
            a[targetParameter] = xx - h;
            fLow[j].initialise1(a);
            a[targetParameter] = xx;
            delta[j] = 2 * h;
        }
        f1.forEach(new ExtendedGradient2Procedure() {

            int i = -1;

            public void executeExtended(double value, double[] dy_da, double[] d2y_dadb) {
                i++;
                DenseMatrix64F m = DenseMatrix64F.wrap(nparams, nparams, d2y_dadb);
                for (int j = 0; j < nparams; j++) {
                    // Evaluate the function +/- delta for parameter j
                    fHigh[j].eval(i, du_da);
                    fLow[j].eval(i, du_db);
                    // Check the gradient with respect to parameter k
                    for (int k = 0; k < nparams; k++) {
                        double gradient = (du_da[k] - du_db[k]) / delta[j];
                        boolean ok = eq.almostEqualRelativeOrAbsolute(gradient, m.get(j, k));
                        if (!ok) {
                            System.out.printf("%d [%d,%d] %f ?= %f\n", i, j, k, gradient, m.get(j, k));
                            Assert.fail(String.format("%d [%d,%d] %f != %f", i, j, k, gradient, m.get(j, k)));
                        }
                    }
                }
            }
        });
    }
}
Also used : ExtendedGradient2Procedure(gdsc.smlm.function.ExtendedGradient2Procedure) DenseMatrix64F(org.ejml.data.DenseMatrix64F) Gaussian2DFunctionTest(gdsc.smlm.function.gaussian.Gaussian2DFunctionTest) Test(org.junit.Test)

Example 2 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class EJMLLinearSolver method invertPseudoInverse.

/**
	 * Invert symmetric positive definite matrix A. On output a replaced by A^-1.
	 * 
	 * @param a
	 *            the matrix a
	 * @return False if the matrix is singular (no solution)
	 */
public boolean invertPseudoInverse(double[][] a) {
    DenseMatrix64F A = toA(a);
    if (!invertPseudoInverse(A))
        return false;
    toSquareData(A, a);
    return true;
}
Also used : DenseMatrix64F(org.ejml.data.DenseMatrix64F)

Example 3 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class EJMLLinearSolver method invertLinear.

/**
	 * Invert symmetric positive definite matrix A. On output a replaced by A^-1.
	 * 
	 * @param a
	 *            the matrix a
	 * @return False if the matrix is singular (no solution)
	 */
public boolean invertLinear(double[][] a) {
    DenseMatrix64F A = toA(a);
    if (!invertLinear(A))
        return false;
    toSquareData(A, a);
    return true;
}
Also used : DenseMatrix64F(org.ejml.data.DenseMatrix64F)

Example 4 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class EJMLLinearSolver method invert.

/**
	 * Computes the inverse of the symmetric positive definite matrix. On output a is replaced by the inverse of a.
	 * <p>
	 * Note: If the matrix is singular then a pseudo inverse will be computed.
	 *
	 * @param a
	 *            the matrix a
	 * @return False if there is no solution
	 */
public boolean invert(double[][] a) {
    DenseMatrix64F A = toA(a);
    if (!invert(A))
        return false;
    toSquareData(A, a);
    return true;
}
Also used : DenseMatrix64F(org.ejml.data.DenseMatrix64F)

Example 5 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class EJMLLinearSolver method invertCholeskyLDLT.

/**
	 * Invert symmetric positive definite matrix A. On output a replaced by A^-1.
	 * 
	 * @param a
	 *            the matrix a
	 * @return False if the matrix is singular (no solution)
	 */
public boolean invertCholeskyLDLT(double[][] a) {
    DenseMatrix64F A = toA(a);
    if (!invertCholeskyLDLT(A))
        return false;
    toSquareData(A, a);
    return true;
}
Also used : DenseMatrix64F(org.ejml.data.DenseMatrix64F)

Aggregations

DenseMatrix64F (org.ejml.data.DenseMatrix64F)20 RandomDataGenerator (org.apache.commons.math3.random.RandomDataGenerator)6 Well19937c (org.apache.commons.math3.random.Well19937c)6 ArrayList (java.util.ArrayList)4 FakeGradientFunction (gdsc.smlm.function.FakeGradientFunction)3 Test (org.junit.Test)3 TimingService (gdsc.core.test.TimingService)2 TurboList (gdsc.core.utils.TurboList)2 GradientCalculator (gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator)2 ExtendedGradient2Procedure (gdsc.smlm.function.ExtendedGradient2Procedure)2 ValueProcedure (gdsc.smlm.function.ValueProcedure)2 Gaussian2DFunction (gdsc.smlm.function.gaussian.Gaussian2DFunction)2 Gaussian2DFunctionTest (gdsc.smlm.function.gaussian.Gaussian2DFunctionTest)2 DoubleEquality (gdsc.core.utils.DoubleEquality)1 ChiSquareTest (org.apache.commons.math3.stat.inference.ChiSquareTest)1