Examples with DenseMatrix64F org.ejml.data.DenseMatrix64F used on opensource projects

Example 16 with DenseMatrix64F

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

the class EJMLLinearSolver method invertSafe.

private boolean invertSafe(LinearSolver<DenseMatrix64F> solver, DenseMatrix64F A, boolean pseudoInverse) {
    DenseMatrix64F Ain = (solver.modifiesA() || isInversionTolerance()) ? A.copy() : A;
    if (!initialiseSolver(solver, Ain))
        return false;
    solver.invert(getA_inv());
    // Check for NaN or Infinity
    double[] a_inv = A_inv.data;
    for (int i = a_inv.length; i-- > 0; ) if (!Maths.isFinite(a_inv[i]))
        return false;
    if (isInversionTolerance() && invalidInversion(A, pseudoInverse))
        return false;
    System.arraycopy(a_inv, 0, A.data, 0, a_inv.length);
    return true;
}

Example 17 with DenseMatrix64F

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

the class FastMLESteppingFunctionSolver method computeStep.

/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.smlm.fitting.nonlinear.SteppingFunctionSolver#computeStep(double[])
	 */
@Override
protected void computeStep(double[] step) {
    final double[] d1 = gradientProcedure.d1;
    final double[] d2 = gradientProcedure.d2;
    if (solver != null) {
        // Should the target for A x = b be created differently?
        for (int i = 0; i < step.length; i++) step[i] = -d1[i];
        jacobianGradientProcedure.getJacobianLinear(jacobian);
        DenseMatrix64F m = DenseMatrix64F.wrap(d1.length, d1.length, jacobian);
        System.out.println(m.toString());
        if (solver.solve(jacobian, step)) {
            // XXX - debug the difference
            double[] step2 = new double[d1.length];
            for (int i = 0; i < step.length; i++) step2[i] = -d1[i] / d2[i];
            System.out.printf("[%d] Jacobian Step %s vs %s\n", tc.getIterations(), Arrays.toString(step), Arrays.toString(step2));
            return;
        }
    }
    // => new parameter = parameter - delta  
    for (int i = 0; i < step.length; i++) step[i] = -d1[i] / d2[i];
}

Example 18 with DenseMatrix64F

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

the class EJMLLinearSolver method invalidInversion.

private boolean invalidInversion(DenseMatrix64F A, boolean pseudoInverse) {
    // Check for the identity matrix:
    // Compute A A_inv = I
    final int n = A.numCols;
    DenseMatrix64F I = new DenseMatrix64F(n, n);
    CommonOps.mult(A, A_inv, I);
    if (pseudoInverse) {
        for (int i = n, index = I.data.length; i-- > 0; ) {
            for (int j = n; j-- > 0; ) {
                if (j == i) {
                    --index;
                    // If using the pseudo inverse then the diagonal can be zero or 1
                    if (invalid(I.data[index], 1) && invalid(I.data[index], 0))
                        return true;
                } else {
                    if (invalid(I.data[--index], 0))
                        return true;
                }
            }
        }
    } else {
        for (int i = n, index = I.data.length; i-- > 0; ) {
            for (int j = n; j-- > 0; ) {
                if (invalid(I.data[--index], (j == i) ? 1 : 0))
                    return true;
            }
        }
    }
    return false;
}

Example 19 with DenseMatrix64F

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

the class EJMLLinearSolver method invertDirectInversion.

/**
	 * 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 invertDirectInversion(double[][] a) {
    DenseMatrix64F A = toA(a);
    if (!invertDirectInversion(A))
        return false;
    toSquareData(A, a);
    return true;
}

Example 20 with DenseMatrix64F

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

the class EJMLLinearSolverTest method runSolverSpeedTest.

private void runSolverSpeedTest(int flags) {
    final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
    int n = f0.size();
    final double[] y = new double[n];
    final TurboList<DenseMatrix64F> aList = new TurboList<DenseMatrix64F>();
    final TurboList<DenseMatrix64F> bList = new TurboList<DenseMatrix64F>();
    double[] testbackground = new double[] { 0.2, 0.7 };
    double[] testsignal1 = new double[] { 30, 100, 300 };
    double[] testcx1 = new double[] { 4.9, 5.3 };
    double[] testcy1 = new double[] { 4.8, 5.2 };
    double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
    int np = f0.getNumberOfGradients();
    GradientCalculator calc = GradientCalculatorFactory.newCalculator(np);
    final RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    //double lambda = 10;
    for (double background : testbackground) // Peak 1
    for (double signal1 : testsignal1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double w1 : testw1) {
        double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
        f0.initialise(p);
        f0.forEach(new ValueProcedure() {

            int i = 0;

            public void execute(double value) {
                // Poisson data 
                y[i++] = rdg.nextPoisson(value);
            }
        });
        double[][] alpha = new double[np][np];
        double[] beta = new double[np];
        //double ss = 
        calc.findLinearised(n, y, p, alpha, beta, f0);
        //System.out.printf("SS = %f\n", ss);
        // As per the LVM algorithm
        //for (int i = 0; i < np; i++)
        //	alpha[i][i] *= lambda;
        aList.add(EJMLLinearSolver.toA(alpha));
        bList.add(EJMLLinearSolver.toB(beta));
    }
    DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[aList.size()]);
    DenseMatrix64F[] b = bList.toArray(new DenseMatrix64F[bList.size()]);
    int runs = 100000 / a.length;
    TimingService ts = new TimingService(runs);
    TurboList<SolverTimingTask> tasks = new TurboList<SolverTimingTask>();
    tasks.add(new PseudoInverseSolverTimingTask(a, b));
    tasks.add(new LinearSolverTimingTask(a, b));
    tasks.add(new CholeskySolverTimingTask(a, b));
    tasks.add(new CholeskyLDLTSolverTimingTask(a, b));
    tasks.add(new DirectInversionSolverTimingTask(a, b));
    for (SolverTimingTask task : tasks) if (!task.badSolver)
        ts.execute(task);
    ts.repeat();
    ts.report();
}