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