List of usage examples for org.apache.commons.math3.linear LUDecomposition LUDecomposition
public LUDecomposition(RealMatrix matrix, double singularityThreshold)
From source file:edu.cmu.tetrad.util.TetradMatrix1.java
public TetradMatrix1 inverse() { if (!isSquare()) throw new IllegalArgumentException("I can only invert square matrices."); // Trying for a speedup by not having to construct the matrix factorization. if (rows() == 0) { return new TetradMatrix1(0, 0); } else if (rows() == 1) { TetradMatrix1 m = new TetradMatrix1(1, 1); m.set(0, 0, 1.0 / apacheData.getEntry(0, 0)); return m; } else if (rows() == 2) { double a = apacheData.getEntry(0, 0); double b = apacheData.getEntry(0, 1); double c = apacheData.getEntry(1, 0); double d = apacheData.getEntry(1, 1); double delta = a * d - b * c; TetradMatrix1 inverse = new TetradMatrix1(2, 2); inverse.set(0, 0, d);/* w w w . j av a 2 s . c o m*/ inverse.set(0, 1, -b); inverse.set(1, 0, -c); inverse.set(1, 1, a); return inverse.scalarMult(1.0 / delta); } else if (rows() == 3) { RealMatrix m = apacheData; double a11 = m.getEntry(0, 0); double a12 = m.getEntry(0, 1); double a13 = m.getEntry(0, 2); double a21 = m.getEntry(1, 0); double a22 = m.getEntry(1, 1); double a23 = m.getEntry(1, 2); double a31 = m.getEntry(2, 0); double a32 = m.getEntry(2, 1); double a33 = m.getEntry(2, 2); final double denom = -a12 * a21 * a33 + a11 * a22 * a33 - a13 * a22 * a31 + a12 * a23 * a31 + a13 * a21 * a32 - a11 * a23 * a32; double[][] inverse = new double[][] { { (a22 * a33 - a23 * a32) / denom, (-a12 * a33 + a13 * a32) / denom, (-a13 * a22 + a12 * a23) / denom }, { (-a21 * a33 + a23 * a31) / denom, (a11 * a33 - a13 * a31) / denom, (a13 * a21 - a11 * a23) / denom }, { (-a22 * a31 + a21 * a32) / denom, (a12 * a31 - a11 * a32) / denom, (-a12 * a21 + a11 * a22) / denom } }; return new TetradMatrix1(inverse); } else if (rows() == 4) { RealMatrix m = apacheData; double a11 = m.getEntry(0, 0); double a12 = m.getEntry(0, 1); double a13 = m.getEntry(0, 2); double a14 = m.getEntry(0, 3); double a21 = m.getEntry(1, 0); double a22 = m.getEntry(1, 1); double a23 = m.getEntry(1, 2); double a24 = m.getEntry(1, 3); double a31 = m.getEntry(2, 0); double a32 = m.getEntry(2, 1); double a33 = m.getEntry(2, 2); double a34 = m.getEntry(2, 3); double a41 = m.getEntry(3, 0); double a42 = m.getEntry(3, 1); double a43 = m.getEntry(3, 2); double a44 = m.getEntry(3, 3); final double denom = a14 * a23 * a32 * a41 - a13 * a24 * a32 * a41 - a14 * a22 * a33 * a41 + a12 * a24 * a33 * a41 + a13 * a22 * a34 * a41 - a12 * a23 * a34 * a41 - a14 * a23 * a31 * a42 + a13 * a24 * a31 * a42 + a14 * a21 * a33 * a42 - a11 * a24 * a33 * a42 - a13 * a21 * a34 * a42 + a11 * a23 * a34 * a42 + a14 * a22 * a31 * a43 - a12 * a24 * a31 * a43 - a14 * a21 * a32 * a43 + a11 * a24 * a32 * a43 + a12 * a21 * a34 * a43 - a11 * a22 * a34 * a43 - a13 * a22 * a31 * a44 + a12 * a23 * a31 * a44 + a13 * a21 * a32 * a44 - a11 * a23 * a32 * a44 - a12 * a21 * a33 * a44 + a11 * a22 * a33 * a44; double[][] inverse = new double[][] { { (-a24 * a33 * a42 + a23 * a34 * a42 + a24 * a32 * a43 - a22 * a34 * a43 - a23 * a32 * a44 + a22 * a33 * a44) / denom, (a14 * a33 * a42 - a13 * a34 * a42 - a14 * a32 * a43 + a12 * a34 * a43 + a13 * a32 * a44 - a12 * a33 * a44) / denom, (-a14 * a23 * a42 + a13 * a24 * a42 + a14 * a22 * a43 - a12 * a24 * a43 - a13 * a22 * a44 + a12 * a23 * a44) / denom, (a14 * a23 * a32 - a13 * a24 * a32 - a14 * a22 * a33 + a12 * a24 * a33 + a13 * a22 * a34 - a12 * a23 * a34) / denom }, { (a24 * a33 * a41 - a23 * a34 * a41 - a24 * a31 * a43 + a21 * a34 * a43 + a23 * a31 * a44 - a21 * a33 * a44) / denom, (-a14 * a33 * a41 + a13 * a34 * a41 + a14 * a31 * a43 - a11 * a34 * a43 - a13 * a31 * a44 + a11 * a33 * a44) / denom, (a14 * a23 * a41 - a13 * a24 * a41 - a14 * a21 * a43 + a11 * a24 * a43 + a13 * a21 * a44 - a11 * a23 * a44) / denom, (-a14 * a23 * a31 + a13 * a24 * a31 + a14 * a21 * a33 - a11 * a24 * a33 - a13 * a21 * a34 + a11 * a23 * a34) / denom }, { (-a24 * a32 * a41 + a22 * a34 * a41 + a24 * a31 * a42 - a21 * a34 * a42 - a22 * a31 * a44 + a21 * a32 * a44) / denom, (a14 * a32 * a41 - a12 * a34 * a41 - a14 * a31 * a42 + a11 * a34 * a42 + a12 * a31 * a44 - a11 * a32 * a44) / denom, (-a14 * a22 * a41 + a12 * a24 * a41 + a14 * a21 * a42 - a11 * a24 * a42 - a12 * a21 * a44 + a11 * a22 * a44) / denom, (a14 * a22 * a31 - a12 * a24 * a31 - a14 * a21 * a32 + a11 * a24 * a32 + a12 * a21 * a34 - a11 * a22 * a34) / denom }, { (a23 * a32 * a41 - a22 * a33 * a41 - a23 * a31 * a42 + a21 * a33 * a42 + a22 * a31 * a43 - a21 * a32 * a43) / denom, (-a13 * a32 * a41 + a12 * a33 * a41 + a13 * a31 * a42 - a11 * a33 * a42 - a12 * a31 * a43 + a11 * a32 * a43) / denom, (a13 * a22 * a41 - a12 * a23 * a41 - a13 * a21 * a42 + a11 * a23 * a42 + a12 * a21 * a43 - a11 * a22 * a43) / denom, (-a13 * a22 * a31 + a12 * a23 * a31 + a13 * a21 * a32 - a11 * a23 * a32 - a12 * a21 * a33 + a11 * a22 * a33) / denom } }; return new TetradMatrix1(inverse); } else { // Using LUDecomposition. // other options: QRDecomposition, CholeskyDecomposition, EigenDecomposition, QRDecomposition, // RRQRDDecomposition, SingularValueDecomposition. Very cool. Also MatrixUtils.blockInverse, // though that can't handle matrices of size 1. Many ways to invert. // Note CholeskyDecomposition only takes inverses of symmetric matrices. // return new TetradMatrix(new CholeskyDecomposition(apacheData).getSolver().getInverse()); // return new TetradMatrix(new EigenDecomposition(apacheData).getSolver().getInverse()); // return new TetradMatrix(new QRDecomposition(apacheData).getSolver().getInverse()); // // return new TetradMatrix(new SingularValueDecomposition(apacheData).getSolver().getInverse()); return new TetradMatrix1(new LUDecomposition(apacheData, 1e-9).getSolver().getInverse()); } }
From source file:nova.core.render.model.MeshModel.java
public Set<Model> flatten(MatrixStack matrixStack) { Set<Model> models = new HashSet<>(); matrixStack.pushMatrix();/*from ww w .j a va2 s . c o m*/ matrixStack.transform(matrix.getMatrix()); //Create a new model with transformation applied. MeshModel transformedModel = clone(); // correct formula for Normal Matrix is transpose(inverse(mat3(model_mat)) // we have to augemnt that to 4x4 RealMatrix normalMatrix3x3 = new LUDecomposition(matrixStack.getMatrix().getSubMatrix(0, 2, 0, 2), 1e-5) .getSolver().getInverse().transpose(); RealMatrix normalMatrix = MatrixUtils.createRealMatrix(4, 4); normalMatrix.setSubMatrix(normalMatrix3x3.getData(), 0, 0); normalMatrix.setEntry(3, 3, 1); transformedModel.faces.stream().forEach(f -> { f.normal = TransformUtil.transform(f.normal, normalMatrix); f.vertices.forEach(v -> v.vec = matrixStack.apply(v.vec)); }); models.add(transformedModel); //Flatten child models models.addAll(children.stream().flatMap(m -> m.flatten(matrixStack).stream()).collect(Collectors.toSet())); matrixStack.popMatrix(); return models; }
From source file:org.drugis.mtc.yadas.MultivariateGaussian.java
public double compute(double[] x, double[] mu, double[][] sigma) { int d = x.length; if (d != mu.length || d != sigma.length) { throw new IllegalArgumentException("All arguments need to be of equal length"); }//from ww w . ja v a 2s. com Array2DRowRealMatrix sigmaM = new Array2DRowRealMatrix(sigma); Array2DRowRealMatrix xM = new Array2DRowRealMatrix(x); Array2DRowRealMatrix muM = new Array2DRowRealMatrix(mu); // Return the log of: // 1/sqrt(2pi^d * det(sigma)) * e^(-.5 (x - mu)' inv(sigma) (x - mu)) // Which is: // -log(sqrt(2pi^d * det(sigma))) + -.5 (x - mu)' inv(sigma) (x - mu) Array2DRowRealMatrix dM = xM.subtract(muM); LUDecomposition sigmaD = new LUDecomposition(sigmaM, Precision.SAFE_MIN); try { RealMatrix sigmaInv = sigmaD.getSolver().getInverse(); return -0.5 * (Math.log(2 * Math.PI) * d + Math.log(sigmaD.getDeterminant()) + dM.transpose().multiply(sigmaInv).multiply(dM).getEntry(0, 0)); } catch (RuntimeException e) { System.out.println(sigmaM); throw e; } /* RealMatrix sigmaInv = sigmaM.inverse(); return -0.5 * ( Math.log(2 * Math.PI) * d + Math.log(sigmaM.getDeterminant()) + dM.transpose().multiply(sigmaInv).multiply(dM).getEntry(0, 0)); */ }
From source file:org.mcisb.util.math.MathUtils.java
/** * /* w w w . j a va 2 s . c o m*/ * @param array * @param singularityThreshold * @return double[] */ public static double[][] inverse(final double[][] array, final double singularityThreshold) { return new LUDecomposition(MatrixUtils.createRealMatrix(array), singularityThreshold).getSolver() .getInverse().getData(); }