List of usage examples for org.apache.commons.math3.linear DecompositionSolver solve
RealMatrix solve(final RealMatrix b);
From source file:org.opentestsystem.airose.regression.orderedprobit.LinearEquationSolver.java
public static double[][] SolveLinear(double[][] XpX, double[][] XpY) { RealMatrix coefficients = new Array2DRowRealMatrix(XpX, false); LUDecomposition luDecomposition = new LUDecomposition(coefficients); DecompositionSolver solver = luDecomposition.getSolver(); double[] yRowVector = new double[XpY.length]; for (int counter1 = 0; counter1 < XpY.length; ++counter1) { yRowVector[counter1] = XpY[counter1][0]; }//ww w. j a v a2s .co m RealVector constants = new ArrayRealVector(yRowVector, false); RealVector solution = solver.solve(constants); double[][] result = new double[solution.getDimension()][1]; for (int counter1 = 0; counter1 < solution.getDimension(); ++counter1) { result[counter1][0] = solution.getEntry(counter1); } return result; }
From source file:org.orekit.propagation.numerical.JacobiansMapper.java
/** Set the Jacobian with respect to state into a one-dimensional additional state array. * <p>/* w ww.j a va 2s .c om*/ * This method converts the Jacobians to cartesian parameters and put the converted data * in the one-dimensional {@code p} array. * </p> * @param state spacecraft state * @param dY1dY0 Jacobian of current state at time t? * with respect to state at some previous time t * @param dY1dP Jacobian of current state at time t? * with respect to parameters (may be null if there are no parameters) * @param p placeholder where to put the one-dimensional additional state * @see #getStateJacobian(SpacecraftState, double[][]) */ void setInitialJacobians(final SpacecraftState state, final double[][] dY1dY0, final double[][] dY1dP, final double[] p) { // set up a converter between state parameters and cartesian parameters final RealMatrix dY1dC1 = new Array2DRowRealMatrix(getdYdC(state), false); final DecompositionSolver solver = new QRDecomposition(dY1dC1).getSolver(); // convert the provided state Jacobian to cartesian parameters final RealMatrix dC1dY0 = solver.solve(new Array2DRowRealMatrix(dY1dY0, false)); // map the converted state Jacobian to one-dimensional array int index = 0; for (int i = 0; i < stateDimension; ++i) { for (int j = 0; j < stateDimension; ++j) { p[index++] = dC1dY0.getEntry(i, j); } } if (parameters > 0) { // convert the provided state Jacobian to cartesian parameters final RealMatrix dC1dP = solver.solve(new Array2DRowRealMatrix(dY1dP, false)); // map the converted parameters Jacobian to one-dimensional array for (int i = 0; i < stateDimension; ++i) { for (int j = 0; j < parameters; ++j) { p[index++] = dC1dP.getEntry(i, j); } } } }
From source file:org.orekit.utils.AngularCoordinates.java
/** Find a vector from two known cross products. * <p>/* w w w . ja v a 2 s . co m*/ * We want to find such that: v? = c? and v = c * </p> * <p> * The first equation ( v? = c?) will always be fulfilled exactly, * and the second one will be fulfilled if possible. * </p> * @param v1 vector forming the first known cross product * @param c1 know vector for cross product v? * @param v2 vector forming the second known cross product * @param c2 know vector for cross product v * @param tolerance relative tolerance factor used to check singularities * @return vector such that: v? = c? and v = c * @exception MathIllegalArgumentException if vectors are inconsistent and * no solution can be found */ private static Vector3D inverseCrossProducts(final Vector3D v1, final Vector3D c1, final Vector3D v2, final Vector3D c2, final double tolerance) throws MathIllegalArgumentException { final double v12 = v1.getNormSq(); final double v1n = FastMath.sqrt(v12); final double v22 = v2.getNormSq(); final double v2n = FastMath.sqrt(v22); final double threshold = tolerance * FastMath.max(v1n, v2n); Vector3D omega; try { // create the over-determined linear system representing the two cross products final RealMatrix m = MatrixUtils.createRealMatrix(6, 3); m.setEntry(0, 1, v1.getZ()); m.setEntry(0, 2, -v1.getY()); m.setEntry(1, 0, -v1.getZ()); m.setEntry(1, 2, v1.getX()); m.setEntry(2, 0, v1.getY()); m.setEntry(2, 1, -v1.getX()); m.setEntry(3, 1, v2.getZ()); m.setEntry(3, 2, -v2.getY()); m.setEntry(4, 0, -v2.getZ()); m.setEntry(4, 2, v2.getX()); m.setEntry(5, 0, v2.getY()); m.setEntry(5, 1, -v2.getX()); final RealVector rhs = MatrixUtils.createRealVector( new double[] { c1.getX(), c1.getY(), c1.getZ(), c2.getX(), c2.getY(), c2.getZ() }); // find the best solution we can final DecompositionSolver solver = new QRDecomposition(m, threshold).getSolver(); final RealVector v = solver.solve(rhs); omega = new Vector3D(v.getEntry(0), v.getEntry(1), v.getEntry(2)); } catch (SingularMatrixException sme) { // handle some special cases for which we can compute a solution final double c12 = c1.getNormSq(); final double c1n = FastMath.sqrt(c12); final double c22 = c2.getNormSq(); final double c2n = FastMath.sqrt(c22); if (c1n <= threshold && c2n <= threshold) { // simple special case, velocities are cancelled return Vector3D.ZERO; } else if (v1n <= threshold && c1n >= threshold) { // this is inconsistent, if v? is zero, c? must be 0 too throw new NumberIsTooLargeException(c1n, 0, true); } else if (v2n <= threshold && c2n >= threshold) { // this is inconsistent, if v is zero, c must be 0 too throw new NumberIsTooLargeException(c2n, 0, true); } else if (Vector3D.crossProduct(v1, v2).getNorm() <= threshold && v12 > threshold) { // simple special case, v is redundant with v?, we just ignore it // use the simplest : orthogonal to both v? and c? omega = new Vector3D(1.0 / v12, Vector3D.crossProduct(v1, c1)); } else { throw sme; } } // check results final double d1 = Vector3D.distance(Vector3D.crossProduct(omega, v1), c1); if (d1 > threshold) { throw new NumberIsTooLargeException(d1, 0, true); } final double d2 = Vector3D.distance(Vector3D.crossProduct(omega, v2), c2); if (d2 > threshold) { throw new NumberIsTooLargeException(d2, 0, true); } return omega; }
From source file:pl.matrix.core.MatrixCalculator.java
public Matrix solveLinearEquation(Matrix a, Matrix b) { RealMatrix aRealMatrix = a.toRealMatrix(); RealMatrix bRealMatrix = b.toRealMatrix(); DecompositionSolver decompositionSolver = new LUDecomposition(aRealMatrix).getSolver(); Matrix result = new Matrix(decompositionSolver.solve(bRealMatrix)); return result; }
From source file:plugins.ImageRectificationPanel.java
public void calculateEquations(double[] imageX, double[] imageY, double[] mapX, double[] mapY) { try {/*from w w w.j av a 2s. c o m*/ int m, i, j, k; int n = mapX.length; // How many coefficients are there? numCoefficients = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { numCoefficients++; } } for (i = 0; i < n; i++) { imageX[i] -= imageXMin; imageY[i] -= imageYMin; mapX[i] -= mapXMin; mapY[i] -= mapYMin; } // Solve the forward transformation equations double[][] forwardCoefficientMatrix = new double[n][numCoefficients]; for (i = 0; i < n; i++) { m = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { forwardCoefficientMatrix[i][m] = Math.pow(imageX[i], j) * Math.pow(imageY[i], k); m++; } } } RealMatrix coefficients = new Array2DRowRealMatrix(forwardCoefficientMatrix, false); //DecompositionSolver solver = new SingularValueDecomposition(coefficients).getSolver(); DecompositionSolver solver = new QRDecomposition(coefficients).getSolver(); // do the x-coordinate first RealVector constants = new ArrayRealVector(mapX, false); RealVector solution = solver.solve(constants); forwardRegressCoeffX = new double[n]; for (int a = 0; a < numCoefficients; a++) { forwardRegressCoeffX[a] = solution.getEntry(a); } double[] residualsX = new double[n]; double SSresidX = 0; for (i = 0; i < n; i++) { double yHat = 0.0; for (j = 0; j < numCoefficients; j++) { yHat += forwardCoefficientMatrix[i][j] * forwardRegressCoeffX[j]; } residualsX[i] = mapX[i] - yHat; SSresidX += residualsX[i] * residualsX[i]; } double sumX = 0; double SSx = 0; for (i = 0; i < n; i++) { SSx += mapX[i] * mapX[i]; sumX += mapX[i]; } double varianceX = (SSx - (sumX * sumX) / n) / n; double SStotalX = (n - 1) * varianceX; double rsqX = 1 - SSresidX / SStotalX; //System.out.println("x-coordinate r-square: " + rsqX); // now the y-coordinate constants = new ArrayRealVector(mapY, false); solution = solver.solve(constants); forwardRegressCoeffY = new double[numCoefficients]; for (int a = 0; a < numCoefficients; a++) { forwardRegressCoeffY[a] = solution.getEntry(a); } double[] residualsY = new double[n]; residualsXY = new double[n]; double SSresidY = 0; for (i = 0; i < n; i++) { double yHat = 0.0; for (j = 0; j < numCoefficients; j++) { yHat += forwardCoefficientMatrix[i][j] * forwardRegressCoeffY[j]; } residualsY[i] = mapY[i] - yHat; SSresidY += residualsY[i] * residualsY[i]; residualsXY[i] = Math.sqrt(residualsX[i] * residualsX[i] + residualsY[i] * residualsY[i]); } double sumY = 0; double sumR = 0; double SSy = 0; double SSr = 0; for (i = 0; i < n; i++) { SSy += mapY[i] * mapY[i]; SSr += residualsXY[i] * residualsXY[i]; sumY += mapY[i]; sumR += residualsXY[i]; } double varianceY = (SSy - (sumY * sumY) / n) / n; double varianceResiduals = (SSr - (sumR * sumR) / n) / n; double SStotalY = (n - 1) * varianceY; double rsqY = 1 - SSresidY / SStotalY; overallRMSE = Math.sqrt(varianceResiduals); //System.out.println("y-coordinate r-square: " + rsqY); // // Print the residuals. // System.out.println("\nResiduals:"); // for (i = 0; i < n; i++) { // System.out.println("Point " + (i + 1) + "\t" + residualsX[i] // + "\t" + residualsY[i] + "\t" + residualsXY[i]); // } // Solve the backward transformation equations double[][] backCoefficientMatrix = new double[n][numCoefficients]; for (i = 0; i < n; i++) { m = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { backCoefficientMatrix[i][m] = Math.pow(mapX[i], j) * Math.pow(mapY[i], k); m++; } } } coefficients = new Array2DRowRealMatrix(backCoefficientMatrix, false); //DecompositionSolver solver = new SingularValueDecomposition(coefficients).getSolver(); solver = new QRDecomposition(coefficients).getSolver(); // do the x-coordinate first constants = new ArrayRealVector(imageX, false); solution = solver.solve(constants); backRegressCoeffX = new double[numCoefficients]; for (int a = 0; a < numCoefficients; a++) { backRegressCoeffX[a] = solution.getEntry(a); } // now the y-coordinate constants = new ArrayRealVector(imageY, false); solution = solver.solve(constants); backRegressCoeffY = new double[n]; for (int a = 0; a < numCoefficients; a++) { backRegressCoeffY[a] = solution.getEntry(a); } } catch (OutOfMemoryError oe) { myHost.showFeedback("An out-of-memory error has occurred during operation."); } catch (Exception e) { myHost.showFeedback("An error has occurred during operation. See log file for details."); myHost.logException("Error in ImageRectification", e); } }
From source file:plugins.TrendSurface.java
public double calculateEquation(double[] X, double[] Y, double[] Z) { try {/*from w w w .j ava 2 s.c om*/ int m, i, j, k; int n = Z.length; // How many coefficients are there? numCoefficients = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { numCoefficients++; } } // Solve the forward transformation equations double[][] forwardCoefficientMatrix = new double[n][numCoefficients]; for (i = 0; i < n; i++) { m = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { forwardCoefficientMatrix[i][m] = Math.pow(X[i], j) * Math.pow(Y[i], k); m++; } } } RealMatrix coefficients = new Array2DRowRealMatrix(forwardCoefficientMatrix, false); //DecompositionSolver solver = new SingularValueDecomposition(coefficients).getSolver(); DecompositionSolver solver = new QRDecomposition(coefficients).getSolver(); // do the x-coordinate first RealVector constants = new ArrayRealVector(Z, false); RealVector solution = solver.solve(constants); regressCoefficents = new double[numCoefficients]; for (int a = 0; a < numCoefficients; a++) { regressCoefficents[a] = solution.getEntry(a); } double[] residuals = new double[n]; double SSresid = 0; for (i = 0; i < n; i++) { double yHat = 0.0; for (j = 0; j < numCoefficients; j++) { yHat += forwardCoefficientMatrix[i][j] * regressCoefficents[j]; } residuals[i] = Z[i] - yHat; SSresid += residuals[i] * residuals[i]; } double sum = 0; double SS = 0; for (i = 0; i < n; i++) { SS += Z[i] * Z[i]; sum += Z[i]; } double variance = (SS - (sum * sum) / n) / n; double SStotal = (n - 1) * variance; double rsq = 1 - SSresid / SStotal; return rsq; } catch (DimensionMismatchException | NoDataException | NullArgumentException | OutOfRangeException e) { showFeedback("Error in TrendSurface.calculateEquation: " + e.toString()); return -1; } }
From source file:plugins.TrendSurfaceVectorPoints.java
public double calculateEquation(double[] X, double[] Y, double[] Z) { try {//from w w w . ja va2 s . co m int m, i, j, k; int n = Z.length; // How many coefficients are there? numCoefficients = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { numCoefficients++; } } // Solve the forward transformation equations double[][] forwardCoefficientMatrix = new double[n][numCoefficients]; for (i = 0; i < n; i++) { m = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { forwardCoefficientMatrix[i][m] = Math.pow(X[i], j) * Math.pow(Y[i], k); m++; } } } RealMatrix coefficients = new Array2DRowRealMatrix(forwardCoefficientMatrix, false); //DecompositionSolver solver = new SingularValueDecomposition(coefficients).getSolver(); DecompositionSolver solver = new QRDecomposition(coefficients).getSolver(); // do the z-coordinate RealVector constants = new ArrayRealVector(Z, false); RealVector solution = solver.solve(constants); regressCoefficents = new double[numCoefficients]; for (int a = 0; a < numCoefficients; a++) { regressCoefficents[a] = solution.getEntry(a); } double[] residuals = new double[n]; double SSresid = 0; for (i = 0; i < n; i++) { double yHat = 0.0; for (j = 0; j < numCoefficients; j++) { yHat += forwardCoefficientMatrix[i][j] * regressCoefficents[j]; } residuals[i] = Z[i] - yHat; SSresid += residuals[i] * residuals[i]; } double sum = 0; double SS = 0; for (i = 0; i < n; i++) { SS += Z[i] * Z[i]; sum += Z[i]; } double variance = (SS - (sum * sum) / n) / n; double SStotal = (n - 1) * variance; double rsq = 1 - SSresid / SStotal; return rsq; } catch (DimensionMismatchException | NoDataException | NullArgumentException | OutOfRangeException e) { showFeedback("Error in TrendSurface.calculateEquation: " + e.toString()); return -1; } }
From source file:rcdemo.math.ClosedHermiteSpline.java
public static RealMatrix generate(RealVector p) { int n = p.getDimension(); RealMatrix A = MatrixUtils.createRealMatrix(n, n); RealMatrix B = MatrixUtils.createRealMatrix(n, n); for (int i = 0; i < n; i++) { int ip = (i + 1) % n; int im = (i - 1 + n) % n; A.setEntry(i, im, 2);//from www. j av a 2 s.c o m A.setEntry(i, i, 8); A.setEntry(i, ip, 2); B.setEntry(i, im, -6); B.setEntry(i, ip, 6); } RealVector Bp = B.operate(p); DecompositionSolver solver = new CholeskyDecomposition(A).getSolver(); RealVector m = solver.solve(Bp); RealMatrix a = MatrixUtils.createRealMatrix(4, n); for (int i = 0; i < n; i++) { int ip = (i + 1) % n; a.setEntry(0, i, p.getEntry(i)); a.setEntry(1, i, m.getEntry(i)); a.setEntry(2, i, p.getEntry(ip)); a.setEntry(3, i, m.getEntry(ip)); } return a; }
From source file:resolution.DifferencesFinies.java
public double[] solveCommonMaths(IFonction fonction, double u0, double u1, double[] mesh) { int n = mesh.length; if (n == 0) { return new double[0]; }/* w w w . ja va 2s. c om*/ double b[]; try { b = fonction.evaluer(mesh); } catch (Exception e) { return null; } double h = mesh[0]; for (int i = 0; i < n - 1; i++) { h = Math.max(h, mesh[i + 1] - mesh[i]); } h = Math.max(h, 1 - mesh[n - 1]); double hCarre = h * h; double a[][] = new double[n][n]; double temp1 = -1.0 / hCarre; double temp2 = 2.0 / hCarre; a[0][0] = temp2; if (n > 1) { for (int i = 1; i < n; i++) { a[i][i] = temp2; a[i][i - 1] = temp1; a[i - 1][i] = temp1; } } //Matrix aMatrix = new Matrix(a); // Cration de la matrice A RealMatrix matrice = new Array2DRowRealMatrix(a, false); DecompositionSolver solver = new LUDecomposition(matrice).getSolver(); //Construction du vecteur b b[0] = b[0] + u0 / hCarre; b[n - 1] = b[n - 1] + u1 / hCarre; RealVector bVect = new ArrayRealVector(b, false); //Matrix bMatrix = new Matrix(b, n); //Rsolution de l'quation //RealVector res = solver.solve(bVect); // double[] res1 = res.toArray(); return solver.solve(bVect).toArray(); }
From source file:resolution.VolumesFinis.java
public double[] solveCommonMaths(IFonction fonction, double u0, double u1, double[] mesh) { double[] result = null; int n = mesh.length; if (n == 0) { return result; }/*from w w w. j a v a2s .c o m*/ // Calcul des h(i) et h(i+1/2) double[] hi2 = new double[n + 1]; double[] hi = new double[n]; hi2[0] = mesh[0]; for (int i = 1; i < n; i++) { hi2[i] = mesh[i] - mesh[i - 1]; hi[i - 1] = 0.5 * (hi2[i] + hi2[i - 1]); } hi2[n] = 1 - mesh[n - 1]; hi[n - 1] = 0.5 * (hi2[n] + hi2[n - 1]); // Construction de la matrice A double temp = 1.0 / hi2[0] + 1.0 / hi2[1]; // Cration de la matrice A RealMatrix matrice = new Array2DRowRealMatrix(n, n); matrice.addToEntry(0, 0, temp); if (n > 1) { //matrice.addToEntry(0, 1, -1.0/hi2[0]); for (int i = 1; i < n; i++) { matrice.addToEntry(i, i, (1.0 / hi2[i] + 1.0 / hi2[i + 1])); matrice.addToEntry(i, i - 1, -1.0 / hi2[i]); matrice.addToEntry(i - 1, i, -1.0 / hi2[i + 1]); } } DecompositionSolver solver = new LUDecomposition(matrice).getSolver(); //Construction du vecteur b double b[]; try { b = fonction.evaluer(mesh); for (int i = 0; i < n; i++) { b[i] = b[i] * hi[i]; } } catch (Exception e) { return null; } b[0] = b[0] + u0 / hi2[0]; b[n - 1] = b[n - 1] + u1 / hi2[n - 1]; /* RealVector bVect = new ArrayRealVector(n); for (int i=0; i<n; i++) bVect.setEntry(i, b[i]); //*/ ArrayRealVector bVect = new ArrayRealVector(b); //Matrix bMatrix = new Matrix(b, n); //Rsolution de l'quation //RealVector res = solver.solve(bVect); // double[] res1 = res.toArray(); return solver.solve(bVect).toArray(); }