/** * @param p * @return the p-th root matrix of this matrix */ public Matrix root (double p) { if (!isSquare()) throw new IllegalArgumentException("No sqaure matrix."); if (p==1) return this; if (p<1) throw new IllegalArgumentException(String.format("Illegal fractional root: %f", p)); EigenDecomposition ev = new EigenDecomposition(this); RealMatrix V = ev.getV(); RealMatrix D = ev.getD(); for (int i=0; i<D.getRowDimension(); i++) { if (D.getEntry(i,i)<0) throw new IllegalStateException("Root of the matrix is not defined."); D.setEntry(i,i,Math.pow(D.getEntry(i,i),1./p)); } RealMatrix B = V.multiply(D).multiply(new LUDecomposition(V).getSolver().getInverse()); return new Matrix(B.getData()); }
@Override public Optional<Ellipsoid> calculate(final Matrix4d quadricSolution) { Vector3d center; try { center = findCenter(quadricSolution); } catch (SingularMatrixException sme) { return Optional.empty(); } final Matrix4d translated = translateToCenter(quadricSolution, center); final EigenDecomposition decomposition = solveEigenDecomposition( translated); if (!isEllipsoid(decomposition.getRealEigenvalues())) { return Optional.empty(); } final double[] radii = Arrays.stream(decomposition.getRealEigenvalues()) .map(ev -> Math.sqrt(1.0 / ev)).toArray(); final Ellipsoid ellipsoid = new Ellipsoid(radii[0], radii[1], radii[2]); ellipsoid.setCentroid(center); final Matrix3d orientation = toOrientationMatrix(decomposition); ellipsoid.setOrientation(orientation); return Optional.of(ellipsoid); }
private Matrix3d toOrientationMatrix(final EigenDecomposition decomposition) { final RealVector e1 = decomposition.getEigenvector(0); final RealVector e2 = decomposition.getEigenvector(1); final RealVector e3 = decomposition.getEigenvector(2); final Vector3d x = new Vector3d(e1.getEntry(0), e1.getEntry(1), e1.getEntry( 2)); final Vector3d y = new Vector3d(e2.getEntry(0), e2.getEntry(1), e2.getEntry( 2)); final Vector3d z = new Vector3d(e3.getEntry(0), e3.getEntry(1), e3.getEntry( 2)); final Matrix3d orientation = new Matrix3d(); orientation.setColumn(0, x); orientation.setColumn(1, y); orientation.setColumn(2, z); return orientation; }
@Before public void setup() { assertEquals(weights.length, showerPixelIds.length); assertEquals(valuesNotWeighted.length, weights.length); Weighted2dStatistics stats2d = Weighted2dStatistics.ofArrays(pixelX, pixelY, weights); cog = stats2d.mean; covarianceMatrix = stats2d.covarianceMatrix; eig = new EigenDecomposition(covarianceMatrix); //get statistics for (double value : valuesWeighted) { s.addValue(value); } }
/** * Find feature space tensor [ ]. * * @param log the log * @param featureVectors the feature vectors * @param components the components * @return the tensor [ ] */ protected Tensor[] findFeatureSpace(final NotebookOutput log, final Supplier<Stream<Tensor[]>> featureVectors, final int components) { return log.code(() -> { final int column = 1; final Tensor[] prototype = featureVectors.get().findAny().get(); final int[] dimensions = prototype[column].getDimensions(); final EigenDecomposition decomposition = new EigenDecomposition(FindPCAFeatures.getCovariance(() -> featureVectors.get().map(x -> x[column].getData()))); final int[] orderedVectors = IntStream.range(0, components).mapToObj(x -> x) .sorted(Comparator.comparing(x -> -decomposition.getRealEigenvalue(x))).mapToInt(x -> x).toArray(); return IntStream.range(0, orderedVectors.length) .mapToObj(i -> { final Tensor src = new Tensor(decomposition.getEigenvector(orderedVectors[i]).toArray(), dimensions).copy(); return src .scale(1.0 / src.rms()) //.scale((decomposition.getRealEigenvalue(orderedVectors[i]) / decomposition.getRealEigenvalue(orderedVectors[orderedVectors.length - 1]))) .scale(Math.sqrt(6. / (components + prototype[column].dim() + 1))) ; } ).toArray(i -> new Tensor[i]); }); }
@Override public void compute(final Mesh input, final DoubleType output) { final RealMatrix it = inertiaTensor.calculate(input); final EigenDecomposition ed = new EigenDecomposition(it); final double l1 = ed.getRealEigenvalue(0) - ed.getRealEigenvalue(2) + ed.getRealEigenvalue(1); final double l2 = ed.getRealEigenvalue(0) - ed.getRealEigenvalue(1) + ed.getRealEigenvalue(2); final double l3 = ed.getRealEigenvalue(2) - ed.getRealEigenvalue(0) + ed.getRealEigenvalue(1); final double g = 1 / (8 * Math.PI / 15); final double a = Math.pow(g * l1 * l1 / Math.sqrt(l2 * l3), 1 / 5d); final double b = Math.pow(g * l2 * l2 / Math.sqrt(l1 * l3), 1 / 5d); final double c = Math.pow(g * l3 * l3 / Math.sqrt(l1 * l2), 1 / 5d); double volumeEllipsoid = (4 / 3d * Math.PI * a * b * c); output.set(volume.calculate(input).get() / volumeEllipsoid); }
@Override public Double[] getRoots(RealPolynomialFunction1D function) { ArgChecker.notNull(function, "function"); double[] coeffs = function.getCoefficients(); int l = coeffs.length - 1; double[][] hessianDeref = new double[l][l]; for (int i = 0; i < l; i++) { hessianDeref[0][i] = -coeffs[l - i - 1] / coeffs[l]; for (int j = 1; j < l; j++) { hessianDeref[j][i] = 0; if (i != l - 1) { hessianDeref[i + 1][i] = 1; } } } RealMatrix hessian = new Array2DRowRealMatrix(hessianDeref); double[] d = new EigenDecomposition(hessian).getRealEigenvalues(); Double[] result = new Double[d.length]; for (int i = 0; i < d.length; i++) { result[i] = d[i]; } return result; }
private void computeCriticalDelay() { EigenDecomposition e = new EigenDecomposition(effortMatrix_T); double[] realParts = e.getRealEigenvalues(); double[] imagParts = e.getImagEigenvalues(); double largestAbsoluteEigenvalue = 0; for (int i = 0; i < realParts.length; i++) { if (imagParts[i] < 0.0000000001 & realParts[i] > largestAbsoluteEigenvalue) { largestAbsoluteEigenvalue = realParts[i]; } } this.criticalDelay = largestAbsoluteEigenvalue; }
/** * Conditions the supplied covariance matrix by enforcing * positive eigenvalues along its main diagonal. * @param S original covariance matrix * @return modified covariance matrix */ private RealMatrix conditionCovarianceMatrix(RealMatrix S) { EigenDecomposition ed = new EigenDecomposition(S); // S -> V . D . V^T RealMatrix V = ed.getV(); RealMatrix D = ed.getD(); // diagonal matrix of eigenvalues RealMatrix VT = ed.getVT(); for (int i = 0; i < D.getRowDimension(); i++) { D.setEntry(i, i, Math.max(D.getEntry(i, i), 10E-6)); // setting eigenvalues to zero is not enough! } return V.multiply(D).multiply(VT); }
/** * Constructor * @param matrix the input matrix */ XEVD(RealMatrix matrix) { final EigenDecomposition evd = new EigenDecomposition(matrix); this.d = toDataFrame(evd.getD()); this.v = toDataFrame(evd.getV()); this.eigenValues = Array.of(evd.getRealEigenvalues()); }
/** * @param a * @return dataset of eigenvalues (can be double or complex double) */ public static Dataset calcEigenvalues(Dataset a) { EigenDecomposition evd = new EigenDecomposition(createRealMatrix(a)); double[] rev = evd.getRealEigenvalues(); if (evd.hasComplexEigenvalues()) { double[] iev = evd.getImagEigenvalues(); return DatasetFactory.createComplexDataset(ComplexDoubleDataset.class, rev, iev); } return DatasetFactory.createFromObject(rev); }
/** * Calculate eigen-decomposition A = V D V^T * @param a * @return array of D eigenvalues (can be double or complex double) and V eigenvectors */ public static Dataset[] calcEigenDecomposition(Dataset a) { EigenDecomposition evd = new EigenDecomposition(createRealMatrix(a)); Dataset[] results = new Dataset[2]; double[] rev = evd.getRealEigenvalues(); if (evd.hasComplexEigenvalues()) { double[] iev = evd.getImagEigenvalues(); results[0] = DatasetFactory.createComplexDataset(ComplexDoubleDataset.class, rev, iev); } else { results[0] = DatasetFactory.createFromObject(rev); } results[1] = createDataset(evd.getV()); return results; }
static double[][] matrixLog (double[][] matrix) { if (matrix.length==1) return new double[][]{{Math.log(matrix[0][0])}}; if (matrix.length!=matrix[0].length) throw new IllegalArgumentException("No sqaure matrix."); RealMatrix A = new Array2DRowRealMatrix(matrix); EigenDecomposition ev = new EigenDecomposition(A); RealMatrix V = ev.getV(); RealMatrix Vinv = (new LUDecomposition(V).getSolver().getInverse()); RealMatrix logA = Vinv.multiply(A).multiply(V); for (int i=0; i<matrix.length; i++) logA.setEntry(i, i,Math.log(logA.getEntry(i, i))); return V.multiply(logA).multiply(Vinv).getData(); }
private static EigenDecomposition solveEigenDecomposition( final Matrix4d quadric) { // @formatter:off final RealMatrix input = new Array2DRowRealMatrix(new double[][]{ {quadric.m00, quadric.m01, quadric.m02}, {quadric.m10, quadric.m11, quadric.m12}, {quadric.m20, quadric.m21, quadric.m22}, }).scalarMultiply(-1.0 / quadric.m33); // @formatter:on return new EigenDecomposition(input); }
public double calculateDelta(EigenDecomposition eig) { // calculate the angle between the eigenvector and the camera axis. // So basicly the angle between the major-axis of the ellipse and the // camrera axis. // this will be written in radians. double longitudinalComponent = eig.getEigenvector(0).getEntry(0); double transverseComponent = eig.getEigenvector(0).getEntry(1); return Math.atan(transverseComponent / longitudinalComponent); }
/** * Computes the square-root of the weight matrix. * * @param m Symmetric, positive-definite (weight) matrix. * @return the square-root of the weight matrix. */ private static RealMatrix squareRoot(final RealMatrix m) { if (m instanceof DiagonalMatrix) { final int dim = m.getRowDimension(); final RealMatrix sqrtM = new DiagonalMatrix(dim); for (int i = 0; i < dim; i++) { sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i))); } return sqrtM; } else { final EigenDecomposition dec = new EigenDecomposition(m); return dec.getSquareRoot(); } }
/** * Computes the square-root of the weight matrix. * * @param m Symmetric, positive-definite (weight) matrix. * @return the square-root of the weight matrix. */ private RealMatrix squareRoot(RealMatrix m) { if (m instanceof DiagonalMatrix) { final int dim = m.getRowDimension(); final RealMatrix sqrtM = new DiagonalMatrix(dim); for (int i = 0; i < dim; i++) { sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i))); } return sqrtM; } else { final EigenDecomposition dec = new EigenDecomposition(m); return dec.getSquareRoot(); } }
/** * Fit points to the polynomial expression Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz * + 2Fyz + 2Gx + 2Hy + 2Iz = 1 and determine the center and radii of the * fit ellipsoid. * * @param points * the points to be fit to the ellipsoid. */ public void fitEllipsoid(ArrayList<ThreeSpacePoint> points) { // Fit the points to Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz // + 2Fyz + 2Gx + 2Hy + 2Iz = 1 and solve the system. // v = (( d' * d )^-1) * ( d' * ones.mapAddToSelf(1)); RealVector v = solveSystem(points); // Form the algebraic form of the ellipsoid. RealMatrix a = formAlgebraicMatrix(v); // Find the center of the ellipsoid. center = findCenter(a); // Translate the algebraic form of the ellipsoid to the center. RealMatrix r = translateToCenter(center, a); // Generate a submatrix of r. RealMatrix subr = r.getSubMatrix(0, 2, 0, 2); // subr[i][j] = subr[i][j] / -r[3][3]). double divr = -r.getEntry(3, 3); for (int i = 0; i < subr.getRowDimension(); i++) { for (int j = 0; j < subr.getRowDimension(); j++) { subr.setEntry(i, j, subr.getEntry(i, j) / divr); } } // Get the eigenvalues and eigenvectors. EigenDecomposition ed = new EigenDecomposition(subr, 0); evals = ed.getRealEigenvalues(); evecs = ed.getEigenvector(0); evecs1 = ed.getEigenvector(1); evecs2 = ed.getEigenvector(2); // Find the radii of the ellipsoid. radii = findRadii(evals); }
/** * Calculates the eigen decomposition of a real matrix. The eigen * decomposition of matrix A is a set of two matrices: V and D such that A = * V × D × VT. A, V and D are all m × m matrices. * * @param a Given matrix. * @return Result W/V arrays. */ public static Array[] eigen_bak(Array a) { int m = a.getShape()[0]; Array Wa; Array Va = Array.factory(DataType.DOUBLE, new int[]{m, m}); double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray(a); RealMatrix matrix = new Array2DRowRealMatrix(aa, false); EigenDecomposition decomposition = new EigenDecomposition(matrix); if (decomposition.hasComplexEigenvalues()) { Wa = Array.factory(DataType.OBJECT, new int[]{m}); double[] rev = decomposition.getRealEigenvalues(); double[] iev = decomposition.getImagEigenvalues(); for (int i = 0; i < m; i++) { Wa.setObject(i, new Complex(rev[i], iev[i])); RealVector v = decomposition.getEigenvector(i); for (int j = 0; j < v.getDimension(); j++) { Va.setDouble(j * m + i, v.getEntry(j)); } } } else { RealMatrix V = decomposition.getV(); RealMatrix D = decomposition.getD(); Wa = Array.factory(DataType.DOUBLE, new int[]{m}); for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { Va.setDouble(i * m + (m - j - 1), V.getEntry(i, j)); if (i == j) { Wa.setDouble(m - i - 1, D.getEntry(i, j)); } } } } return new Array[]{Wa, Va}; }
public void componize(double [][] covariate){ // replace the initial data set the same values minus the means from each variable-column for (int j=0; j < covariate[0].length; j++ ) { DescriptiveStatistics stats = new DescriptiveStatistics(); for (int k=0 ; k < covariate.length ; k++){ stats.addValue(covariate[k][j]);} double mean=stats.getMean(); for (int k=0 ; k < covariate.length ; k++){ covariate[k][j]= covariate[k][j]-mean;} // here ends the iteration that replaces all values with the value minus the mean of the respective column } // We get the Covariance matrix for the "adjusted by mean" data set RealMatrix covariance_matrix= new PearsonsCorrelation(covariate).getCorrelationMatrix(); // we get the Eigen values and the Eigen Vectors Via Eigen Value Decomposition. EigenDecomposition Eig = new EigenDecomposition(covariance_matrix, 1); // get the Eigen Values double Eigenvaluess [] =Eig.getRealEigenvalues(); // Get the Eigen vectors RealMatrix Eigenvec = Eig.getV(); double [][] EigenVecors=Eigenvec.getData(); //Sort everything to get the Vectors with the highest significance SortTableDouble sort = new SortTableDouble(); sort.sorting (Eigenvaluess,EigenVecors ); EigenVecor =sort.value_matrix(); Eigenvalues=sort.scanned_values(); Perchentages= new double[Eigenvalues.length]; for (int k=0 ; k < Eigenvalues.length ; k++){ Perchentages[k]= Eigenvalues[k]/Eigenvalues.length ;} }
public double calculateNthMoment(int n){ Array2DRowRealMatrix gyr = RadiusGyrationTensor2D.getRadiusOfGyrationTensor(t); EigenDecomposition eigdec = new EigenDecomposition(gyr); Vector2d eigv = new Vector2d(eigdec.getEigenvector(0).getEntry(0),eigdec.getEigenvector(0).getEntry(1)); double[] projected = new double[t.size()]; for(int i = 0; i < t.size(); i++){ Vector2d pos = new Vector2d(t.get(i).x,t.get(i).y); double v = eigv.dot(pos); projected[i] = v; } Mean m = new Mean(); StandardDeviation s = new StandardDeviation(); double mean = m.evaluate(projected); double sd = s.evaluate(projected); double sumPowN=0; for(int i = 0; i < projected.length; i++){ sumPowN += Math.pow( (projected[i]-mean)/sd, n); } double nThMoment = sumPowN/projected.length; return nThMoment; }
/** * @return Returns an double array with the following elements [0]=asymmetry */ @Override public double[] evaluate() { Array2DRowRealMatrix gyr = RadiusGyrationTensor2D.getRadiusOfGyrationTensor(t); EigenDecomposition eigdec = new EigenDecomposition(gyr); double e1 = eigdec.getRealEigenvalue(0); double e2 = eigdec.getRealEigenvalue(1); double asym = e2/e1; //-1*Math.log(1-Math.pow(e1-e2, 2)/(2*Math.pow(e1+e2, 2))); result = new double[]{asym}; return result; }
/** * @return Returns an double array with the following elements [0]=Asymmetry */ @Override public double[] evaluate() { Array2DRowRealMatrix gyr = RadiusGyrationTensor2D.getRadiusOfGyrationTensor(t); EigenDecomposition eigdec = new EigenDecomposition(gyr); double e1 = eigdec.getRealEigenvalue(0); double e2 = eigdec.getRealEigenvalue(1); double asym = -1*Math.log(1-Math.pow(e1-e2,2)/(2*Math.pow(e1+e2, 2))); result = new double[]{asym}; return result; }
/** * @return Returns an double array with the following elements [0]=Asymmetry */ @Override public double[] evaluate() { Array2DRowRealMatrix gyr = RadiusGyrationTensor2D.getRadiusOfGyrationTensor(t); EigenDecomposition eigdec = new EigenDecomposition(gyr); double e1 = eigdec.getRealEigenvalue(0); double e2 = eigdec.getRealEigenvalue(1); double asym = Math.pow(e1*e1-e2*e2, 2)/(Math.pow(e1*e1+e2*e2, 2)); //-1*Math.log(1-Math.pow(e1-e2, 2)/(2*Math.pow(e1+e2, 2))); result = new double[]{asym}; return result; }
public Wishart(RealMatrix omega, double nu) { this.omega = omega; this.nu = nu; this.D = omega.getColumnDimension(); if (omega.getRowDimension() == 2) { this.logDetOmega = Math.log(omega.getEntry(0, 0) * omega.getEntry(1, 1) - omega.getEntry(1, 0) * omega.getEntry(0, 1)); } else { this.logDetOmega = Math.log((new EigenDecomposition(omega)).getDeterminant()); } }
private void calculateBandwidthAncillaries() { RealMatrix inverseBandwidth; if (bandwidth.getColumnDimension() > 1) { inverseBandwidth = MatrixUtils.blockInverse(bandwidth, (bandwidth.getColumnDimension() - 1) / 2); } else { // Manually invert size 1 x 1 matrix, because block Inverse requires dimensions > 1 inverseBandwidth = bandwidth.copy(); inverseBandwidth.setEntry(0, 0, 1.0 / inverseBandwidth.getEntry(0, 0)); } this.bandwidthToNegativeHalf = (new EigenDecomposition(inverseBandwidth)).getSquareRoot(); this.bandwidthDeterminantSqrt = Math.sqrt((new EigenDecomposition(bandwidth)).getDeterminant()); }
