Examples with DecompositionSolver org.apache.commons.math3.linear.DecompositionSolver used on opensource projects

Example 16 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project imagingbook-common by imagingbook.

the class AffineFit method fit.

@Override
public void fit(List<double[]> X, List<double[]> Y) {
    // fits n-dimensional data sets with affine model
    if (X.size() != Y.size())
        throw new IllegalArgumentException("point sequences X, Y must have same length");
    this.m = X.size();
    this.n = X.get(0).length;
    RealMatrix M = MatrixUtils.createRealMatrix(2 * m, 2 * (n + 1));
    RealVector b = new ArrayRealVector(2 * m);
    // mount matrix M:
    int row = 0;
    for (double[] x : X) {
        for (int j = 0; j < n; j++) {
            M.setEntry(row, j, x[j]);
            M.setEntry(row, n, 1);
            row++;
        }
        for (int j = 0; j < n; j++) {
            M.setEntry(row, j + n + 1, x[j]);
            M.setEntry(row, 2 * n + 1, 1);
            row++;
        }
    }
    // mount vector b
    row = 0;
    for (double[] y : Y) {
        for (int j = 0; j < n; j++) {
            b.setEntry(row, y[j]);
            row++;
        }
    }
    SingularValueDecomposition svd = new SingularValueDecomposition(M);
    DecompositionSolver solver = svd.getSolver();
    RealVector a = solver.solve(b);
    A = makeTransformationMatrix(a);
}

Example 17 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project imagingbook-common by imagingbook.

the class Matrix method solve.

// use general method, i.e. double[][] inverse(double[][] A)
// @Deprecated
// public static double[][] inverse3x3(final double[][] A) {
// double[][] B = duplicate(A);
// final double det = determinant3x3(B);
// if (Math.abs(det) < Arithmetic.EPSILON_DOUBLE)
// return null;
// else {
// final double a00 = B[0][0];
// final double a01 = B[0][1];
// final double a02 = B[0][2];
// final double a10 = B[1][0];
// final double a11 = B[1][1];
// final double a12 = B[1][2];
// final double a20 = B[2][0];
// final double a21 = B[2][1];
// final double a22 = B[2][2];
// B[0][0] = (a11 * a22 - a12 * a21) / det;
// B[0][1] = (a02 * a21 - a01 * a22) / det;
// B[0][2] = (a01 * a12 - a02 * a11) / det;
// 
// B[1][0] = (a12 * a20 - a10 * a22) / det;
// B[1][1] = (a00 * a22 - a02 * a20) / det;
// B[1][2] = (a02 * a10 - a00 * a12) / det;
// 
// B[2][0] = (a10 * a21 - a11 * a20) / det;
// B[2][1] = (a01 * a20 - a00 * a21) / det;
// B[2][2] = (a00 * a11 - a01 * a10) / det;
// return B;
// }
// }
// numerically stable? should be replaced by standard inversion
// @Deprecated
// public static float[][] inverse3x3(final float[][] A) {
// final float[][] B = duplicate(A);
// final double det = determinant3x3(B);
// // IJ.log("   determinant = " + det);
// if (Math.abs(det) < Arithmetic.EPSILON_DOUBLE)
// return null;
// else {
// final double a00 = B[0][0];
// final double a01 = B[0][1];
// final double a02 = B[0][2];
// final double a10 = B[1][0];
// final double a11 = B[1][1];
// final double a12 = B[1][2];
// final double a20 = B[2][0];
// final double a21 = B[2][1];
// final double a22 = B[2][2];
// B[0][0] = (float) ((a11 * a22 - a12 * a21) / det);
// B[0][1] = (float) ((a02 * a21 - a01 * a22) / det);
// B[0][2] = (float) ((a01 * a12 - a02 * a11) / det);
// 
// B[1][0] = (float) ((a12 * a20 - a10 * a22) / det);
// B[1][1] = (float) ((a00 * a22 - a02 * a20) / det);
// B[1][2] = (float) ((a02 * a10 - a00 * a12) / det);
// 
// B[2][0] = (float) ((a10 * a21 - a11 * a20) / det);
// B[2][1] = (float) ((a01 * a20 - a00 * a21) / det);
// B[2][2] = (float) ((a00 * a11 - a01 * a10) / det);
// return B;
// }
// }
// ------------------------------------------------------------------------
// Finds the EXACT solution x for A.x = b
public static double[] solve(final double[][] A, double[] b) {
    RealMatrix AA = MatrixUtils.createRealMatrix(A);
    RealVector bb = MatrixUtils.createRealVector(b);
    DecompositionSolver solver = new LUDecomposition(AA).getSolver();
    double[] x = null;
    try {
        x = solver.solve(bb).toArray();
    } catch (SingularMatrixException e) {
    }
    return x;
}

Example 18 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project systemml by apache.

the class LibCommonsMath method computeSolve.

/**
 * Function to solve a given system of equations.
 *
 * @param in1 matrix object 1
 * @param in2 matrix object 2
 * @return matrix block
 */
private static MatrixBlock computeSolve(MatrixObject in1, MatrixObject in2) {
    Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in1);
    Array2DRowRealMatrix vectorInput = DataConverter.convertToArray2DRowRealMatrix(in2);
    /*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput);
		DecompositionSolver lusolver = ludecompose.getSolver();
		RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/
    // Setup a solver based on QR Decomposition
    QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
    DecompositionSolver solver = qrdecompose.getSolver();
    // Invoke solve
    RealMatrix solutionMatrix = solver.solve(vectorInput);
    return DataConverter.convertToMatrixBlock(solutionMatrix.getData());
}

Example 19 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project FSensor by KalebKE.

the class FitPoints method solveSystem.

/**
 * Solve the polynomial expression Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz +
 * 2Gx + 2Hy + 2Iz from the provided points.
 *
 * @param points the points that will be fit to the polynomial expression.
 * @return the solution vector to the polynomial expression.
 */
private RealVector solveSystem(ArrayList<ThreeSpacePoint> points) {
    // determine the number of points
    int numPoints = points.size();
    // the design matrix
    // size: numPoints x 9
    RealMatrix d = new Array2DRowRealMatrix(numPoints, 9);
    // Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz
    for (int i = 0; i < d.getRowDimension(); i++) {
        double xx = Math.pow(points.get(i).x, 2);
        double yy = Math.pow(points.get(i).y, 2);
        double zz = Math.pow(points.get(i).z, 2);
        double xy = 2 * (points.get(i).x * points.get(i).y);
        double xz = 2 * (points.get(i).x * points.get(i).z);
        double yz = 2 * (points.get(i).y * points.get(i).z);
        double x = 2 * points.get(i).x;
        double y = 2 * points.get(i).y;
        double z = 2 * points.get(i).z;
        d.setEntry(i, 0, xx);
        d.setEntry(i, 1, yy);
        d.setEntry(i, 2, zz);
        d.setEntry(i, 3, xy);
        d.setEntry(i, 4, xz);
        d.setEntry(i, 5, yz);
        d.setEntry(i, 6, x);
        d.setEntry(i, 7, y);
        d.setEntry(i, 8, z);
    }
    // solve the normal system of equations
    // v = (( d' * d )^-1) * ( d' * ones.mapAddToSelf(1));
    // Multiply: d' * d
    RealMatrix dtd = d.transpose().multiply(d);
    // Create a vector of ones.
    RealVector ones = new ArrayRealVector(numPoints);
    ones.mapAddToSelf(1);
    // Multiply: d' * ones.mapAddToSelf(1)
    RealVector dtOnes = d.transpose().operate(ones);
    // Find ( d' * d )^-1
    DecompositionSolver solver = new SingularValueDecomposition(dtd).getSolver();
    RealMatrix dtdi = solver.getInverse();
    // v = (( d' * d )^-1) * ( d' * ones.mapAddToSelf(1));
    return dtdi.operate(dtOnes);
}

Example 20 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project imagingbook-common by imagingbook.

the class EllipseFit5Points method fit.

private double[] fit(Pnt2d[] points) {
    if (points.length < 5) {
        throw new IllegalArgumentException("5 points are needed for " + this.getClass().getSimpleName());
    }
    final int n = points.length;
    // create design matrix X:
    RealMatrix X = new Array2DRowRealMatrix(n, 5);
    RealVector b = new ArrayRealVector(n);
    for (int i = 0; i < n; i++) {
        double x = points[i].getX();
        double y = points[i].getY();
        X.setEntry(i, 0, sqr(x) - sqr(y));
        X.setEntry(i, 1, x * y);
        X.setEntry(i, 2, x);
        X.setEntry(i, 3, y);
        X.setEntry(i, 4, 1);
        b.setEntry(i, -sqr(y));
    }
    // System.out.println("X = \n" + Matrix.toString(X));
    // System.out.println("b = " + Matrix.toString(b));
    LUDecomposition decomp = new LUDecomposition(X);
    // see https://commons.apache.org/proper/commons-math/userguide/linear.html
    DecompositionSolver solver = decomp.getSolver();
    RealVector x = null;
    try {
        x = solver.solve(b);
    } catch (SingularMatrixException e) {
    }
    if (x == null) {
        return null;
    }
    double A = x.getEntry(0);
    double B = x.getEntry(1);
    double C = 1 - A;
    double D = x.getEntry(2);
    double E = x.getEntry(3);
    double F = x.getEntry(4);
    boolean isEllipse = (4 * A * C - sqr(B)) > 0;
    if (isEllipse) {
        double[] p = { A, B, C, D, E, F };
        return p;
    } else {
        return null;
    }
/*
		double[] p = {A,B,C,D,E,F};
		System.out.println("p = " + Matrix.toString(p));
		
		
		
		if (isEllipse) {	
			AlgebraicEllipse ea = new AlgebraicEllipse(p);
			Ellipse eg = Ellipse.from(ea);
			System.out.println("P1: ellipse = " + eg);
		}
		else {
			System.out.println("P1: NO ellipse!");
		}
		
		EllipseFitAlgebraic fit2 = new EllipseFitFitzgibbonStable(points);
		double[] p2 = fit2.getParameters();
		if (p2 == null) {
			System.out.println("P2: NULL ellipse!");
		}
		else {
			boolean isEllipse2 = (4 * p2[0] * p2[2] - sqr(p2[1])) > 0;
			System.out.println("P2: is ellipse = " +isEllipse2);
			if (isEllipse2) {
				Ellipse ell2 = Ellipse.from(fit2.getEllipse());
				System.out.println("P2: error = " + ell2.getError(points));
			}
			
		}
		System.out.println();
		return p;
		*/
}