List of usage examples for org.apache.commons.math3.linear EigenDecomposition EigenDecomposition
public EigenDecomposition(final double[] main, final double[] secondary)
From source file:ellipsoidFit.FitPoints.java
/** * 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./* w w w. j ava2s . co m*/ * * @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); }
From source file:com.trickl.math.lanczos.TridiagonalMatrix.java
protected void compute() { err.clear();//w w w.j a va2 s . co m eigval_distinct.clear(); multiplicty.clear(); err_noghost.clear(); eigval_distinct_noghost.clear(); multiplicty_noghost.clear(); computed = true; int n = alpha.length; EigenDecomposition eigenDecomposition = new EigenDecomposition(alpha, beta); double[] eval = eigenDecomposition.getRealEigenvalues(); Arrays.sort(eval); // Consistent with IETL // tolerance values: multipleTolerance = Math.max(alpha_max, beta_max) * 2 * epsilon * (1000 + n); threshold = Math.max(eval[0], eval[n - 1]); threshold = Math.max(errorTolerance * threshold, 5 * multipleTolerance); // error estimates of eigen values starts: // the unique eigen values selection, their multiplicities and corresponding errors calculation follows: double temp = eval[0]; eigval_distinct.add(eval[0]); int multiple = 1; for (int i = 1; i < n; i++) { double[] eigenvector = eigenDecomposition.getEigenvector(eval.length - i).toArray(); if (Math.abs(eval[i] - temp) > threshold) { eigval_distinct.add(eval[i]); temp = eval[i]; multiplicty.add(multiple); if (multiple > 1) { err.add(0.); } else { err.add(Math.abs(beta[beta.length - 1] * eigenvector[n - 1])); // *beta.rbegin() = betaMplusOne. } multiple = 1; } else { multiple++; } } // for last eigen value. multiplicty.add(multiple); if (multiple > 1) { err.add(.0); } else { double[] eigenvector = eigenDecomposition.getEigenvector(eval.length - n).toArray(); err.add(Math.abs(beta[beta.length - 1] * eigenvector[n - 1])); // *beta.rbegin() = betaMplusOne. } // the unique eigen values selection, their multiplicities and corresponding errors calculation ends. // ghosts calculations starts: double[] beta_g = Arrays.copyOfRange(beta, 1, beta.length); double[] alpha_g = Arrays.copyOfRange(alpha, 1, alpha.length); eigenDecomposition = new EigenDecomposition(alpha_g, beta_g); double[] eval_g = eigenDecomposition.getRealEigenvalues(); Arrays.sort(eval_g); // Consistent with IETL int i = 0, t2 = 0; for (double eigval : eigval_distinct) { if (multiplicty.get(i) == 1) { // test of spuriousness for the eigenvalues whose multiplicity is one. for (int j = t2; j < n - 1; j++, t2++) { // since size of reduced matrix is n-1 if (eval_g[j] - eigval >= multipleTolerance) { break; } if (Math.abs(eigval - eval_g[j]) < multipleTolerance) { multiplicty.set(i, 0); err.set(i, .0); // if eigen value is a ghost => error calculation not required, 0=> ignore error. t2++; break; } } } i++; } i = 0; for (double eigval : eigval_distinct) { if (multiplicty.get(i) != 0) { eigval_distinct_noghost.add(eigval); multiplicty_noghost.add(multiplicty.get(i)); err_noghost.add(err.get(i)); } i++; } }
From source file:com.joptimizer.functions.SDPLogarithmicBarrier.java
/** * Calculates the initial value for the s parameter in Phase I. * Return s so that F(x)-s.I is negative definite * @see "S.Boyd and L.Vandenberghe, Convex Optimization, 11.6.2" * @see "S.Boyd and L.Vandenberghe, Semidefinite programming, 6.1" */// www . ja v a 2s .co m public double calculatePhase1InitialFeasiblePoint(double[] originalNotFeasiblePoint, double tolerance) { RealMatrix F = this.buildS(originalNotFeasiblePoint).scalarMultiply(-1); RealMatrix S = F.scalarMultiply(-1); try { new CholeskyDecomposition(S); //already feasible return -1; } catch (NonPositiveDefiniteMatrixException ee) { //it does NOT mean that F is negative, it can be not definite EigenDecomposition eFact = new EigenDecomposition(F, Double.NaN); double[] eValues = eFact.getRealEigenvalues(); double minEigenValue = eValues[Utils.getMinIndex(eValues)]; return -Math.min(minEigenValue * Math.pow(tolerance, -0.5), 0.); } }
From source file:com.trickl.math.lanczos.LanczosSolver.java
public Pair<DoubleMatrix2D, Info> getEigenvectorsWithInfo(double[] eigenvalues, RandomGenerator randomGenerator, int maxIterations) { List<DoubleMatrix1D> eigvectors = new LinkedList<>(); // contains ritz vectors. List<List<Double>> Tvectors = new LinkedList<>(); // contains // calculation of eigen vectors of T matrix(consists of alphas & betas): int start = 0; int end = eigenvalues.length; int n1 = alpha.length; double mamax, error, lambda; Pair<Double, Double> a_and_b; int ma = 0, deltam; int nth, maMax = 0, count; List<Double> eigenval_a = new LinkedList<>(); List<Double> residuum = new LinkedList<>(); List<ErrorInfo> status = new LinkedList<>(); findM1M2(eigenvalues);//from w w w .ja v a 2 s. c o m int M1_itr = 0; int M2_itr = 0; for (int pos = start; pos < end; ++pos) { int maxcount = 10; lambda = 0; count = 0; ErrorInfo errInf = ErrorInfo.OK; // calculation of ma starts: if (M1[M1_itr] != 0 && M2[M2_itr] != 0) { ma = (3 * (M1[M1_itr]) + M2[M2_itr]) / 4 + 1; deltam = ((3 * (M1[M1_itr]) + 5 * (M2[M2_itr])) / 8 + 1 - ma) / 10 + 1; } else if (M1[M1_itr] != 0 && M2[M2_itr] == 0) { ma = (5 * (M1[M1_itr])) / 4 + 1; mamax = Math.min((11 * n1) / 8 + 12, (13 * (M1[M1_itr])) / 8 + 1); deltam = (int) ((mamax - ma) / 10) + 1; if (maxIterations > 0) { maxcount = maxIterations / deltam; } } else { errInf = ErrorInfo.NO_EIGENVALUE; deltam = 0; ma = 0; } // calculation of ma ends. eigvectors.add(DoubleFactory1D.dense.make(startvector.size())); // new ritz vector is being added in eigvectors. List<Double> Tvector = new LinkedList<>(); Tvectors.add(Tvector); // new T matrix vector is being added in Tvectors. if (ma == 0) { eigvectors.get(eigvectors.size() - 1).assign(0); } if (ma != 0) { double[] eval; double[] nthEigenvector = null; do { if (ma > alpha.length) { // size of T matrix is to be increased. LanczosIteration iter = new LanczosIterationFixed(ma); generateTMatrix(iter); } count++; // on return, z contains all orthonormal eigen vectors of T matrix. EigenDecomposition eigenDecomposition = new EigenDecomposition(Arrays.copyOfRange(alpha, 0, ma), Arrays.copyOfRange(beta, 0, ma)); eval = eigenDecomposition.getRealEigenvalues(); Arrays.sort(eval); // Consistent with IETL // search for the value of nth starts, where nth is the nth eigen vector in z. for (nth = ma - 1; nth >= 0; nth--) { if (Math.abs(eval[nth] - eigenvalues[pos]) <= threshold) { break; } } // search for the value of ith ends, where ith is the ith eigen vector in z. if (nth == -1) { error = 0; ma = 0; eigvectors.get(eigvectors.size() - 1).assign(0); errInf = ErrorInfo.NO_EIGENVALUE; } else { nthEigenvector = eigenDecomposition.getEigenvector(ma - 1 - nth).toArray(); error = Math.abs(beta[ma - 1] * nthEigenvector[ma - 1]); // beta[ma - 1] = betaMplusOne. if (error > errorTolerance) { ma += deltam; } } // end of else } while (error > errorTolerance && count < maxcount); if (error > errorTolerance) { eigvectors.get(eigvectors.size() - 1).assign(0); errInf = ErrorInfo.NOT_CALCULATED; } else { // if error is small enough. if (ma != 0) { for (int i = 0; i < ma; i++) { (Tvectors.get(Tvectors.size() - 1)).add(nthEigenvector[i]); } if (ma > maMax) { maMax = ma; } lambda = eval[nth]; } } } eigenval_a.add(lambda); // for Info object. Ma.add(ma); // for Info object. status.add(errInf); M1_itr++; M2_itr++; } // end of while(in_eigvals_start != in_eigvals_end) // basis transformation of eigen vectors of T. These vectors are good // approximation of eigen vectors of actual matrix. int eigenvectors_itr; int Tvectors_itr; a_and_b = makeFirstStep(randomGenerator); if (a_and_b.getFirst() != alpha[0] || a_and_b.getSecond() != beta[0]) { throw new RuntimeException("T-matrix problem at first step"); } eigenvectors_itr = 0; Tvectors_itr = 0; while (eigenvectors_itr != end) { if (!Tvectors.get(Tvectors_itr).isEmpty()) { DoubleMatrix1D eigvector = eigvectors.get(eigenvectors_itr); eigvector.assign(0); eigvector.assign(startvector, PlusMult.plusMult(Tvectors.get(Tvectors_itr).get(0))); eigvector.assign(vec2, PlusMult.plusMult(Tvectors.get(Tvectors_itr).get(1))); } eigenvectors_itr++; Tvectors_itr++; } vec2index = 2; for (int j = 2; j < maMax; j++) { a_and_b = makeStep(j - 1); if (a_and_b.getFirst() != alpha[j - 1] || a_and_b.getSecond() != beta[j - 1]) { throw new RuntimeException("T-matrix problem"); } ++vec2index; eigenvectors_itr = 0; Tvectors_itr = 0; while (eigenvectors_itr != end) { if (Tvectors.get(Tvectors_itr).size() > j) { DoubleMatrix1D eigvector = eigvectors.get(eigenvectors_itr); eigvector.assign(vec2, PlusMult.plusMult(Tvectors.get(Tvectors_itr).get(j))); } // vec2 is being added in one vector of eigvectors. eigenvectors_itr++; Tvectors_itr++; } } // end of basis transformation. // end of basis transformation. // copying to the output iterator & residuum calculation starts: int i = 0; DoubleMatrix2D eigenvectors = DoubleFactory2D.dense.make(startvector.size(), end - start); for (eigenvectors_itr = start; eigenvectors_itr != end; eigenvectors_itr++) { DoubleMatrix1D eigvector = eigvectors.get(eigenvectors_itr); eigenvectors.viewColumn(i).assign(eigvector); matrix.zMult(eigvector, vec3); vec3.assign(eigvector, PlusMult.minusMult(eigenval_a.get(i))); // now vec3 is (A*v - eigenval_a*v); *eigenvectors_itr) is being added in vec3. residuum.add(norm2(vec3)); i++; } // copying to the output iterator ends. Info info = new Info(M1, M2, Ma, eigenval_a, residuum, status); return new Pair<>(eigenvectors, info); }
From source file:ffx.autoparm.Superpose.java
/** * <p>// www . jav a 2 s .c o m * quatfit</p> */ public void quatfit() { double xxyx = 0.0; double xxyy = 0.0; double xxyz = 0.0; double xyyx = 0.0; double xyyy = 0.0; double xyyz = 0.0; double xzyx = 0.0; double xzyy = 0.0; double xzyz = 0.0; double c[][] = new double[4][4]; for (int i = 0; i < nfit; i++) { xxyx = xxyx + xyz1[i][0] * xyz2[i][0]; xxyy = xxyy + xyz1[i][1] * xyz2[i][0]; xxyz = xxyz + xyz1[i][2] * xyz2[i][0]; xyyx = xyyx + xyz1[i][0] * xyz2[i][1]; xyyy = xyyy + xyz1[i][1] * xyz2[i][1]; xyyz = xyyz + xyz1[i][2] * xyz2[i][1]; xzyx = xzyx + xyz1[i][0] * xyz2[i][2]; xzyy = xzyy + xyz1[i][1] * xyz2[i][2]; xzyz = xzyz + xyz1[i][2] * xyz2[i][2]; } c[0][0] = xxyx + xyyy + xzyz; c[0][1] = xzyy - xyyz; c[1][1] = xxyx - xyyy - xzyz; c[0][2] = xxyz - xzyx; c[1][2] = xxyy + xyyx; c[2][2] = xyyy - xzyz - xxyx; c[0][3] = xyyx - xxyy; c[1][3] = xzyx + xxyz; c[2][3] = xyyz + xzyy; c[3][3] = xzyz - xxyx - xyyy; RealMatrix a = new Array2DRowRealMatrix( new double[][] { { c[0][0], c[0][1], c[0][2], c[0][3] }, { c[0][1], c[1][1], c[1][2], c[1][3] }, { c[0][2], c[1][2], c[2][2], c[2][3] }, { c[0][3], c[1][3], c[2][3], c[3][3] } }); EigenDecomposition e = new EigenDecomposition(a, 1); a = e.getV(); double[] q = { a.getEntry(0, 0), a.getEntry(1, 0), a.getEntry(2, 0), a.getEntry(3, 0) }; double rot[][] = new double[3][3]; rot[0][0] = q[0] * q[0] + q[1] * q[1] - q[2] * q[2] - q[3] * q[3]; rot[1][0] = 2.0 * (q[1] * q[2] - q[0] * q[3]); rot[2][0] = 2.0 * (q[1] * q[3] + q[0] * q[2]); rot[0][1] = 2.0 * (q[2] * q[1] + q[0] * q[3]); rot[1][1] = q[0] * q[0] - q[1] * q[1] + q[2] * q[2] - q[3] * q[3]; rot[2][1] = 2.0 * (q[2] * q[3] - q[0] * q[1]); rot[0][2] = 2.0 * (q[3] * q[1] - q[0] * q[2]); rot[1][2] = 2.0 * (q[3] * q[2] + q[0] * q[1]); rot[2][2] = q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3] * q[3]; double xrot, yrot, zrot; for (int i = 0; i < n2; i++) { xrot = xyz2[i][0] * rot[0][0] + xyz2[i][1] * rot[0][1] + xyz2[i][2] * rot[0][2]; yrot = xyz2[i][0] * rot[1][0] + xyz2[i][1] * rot[1][1] + xyz2[i][2] * rot[1][2]; zrot = xyz2[i][0] * rot[2][0] + xyz2[i][1] * rot[2][1] + xyz2[i][2] * rot[2][2]; xyz2[i][0] = xrot; xyz2[i][1] = yrot; xyz2[i][2] = zrot; } }
From source file:edu.ucla.stat.SOCR.analyses.gui.PrincipalComponentAnalysis.java
/**This method defines the specific statistical Analysis to be carried our on the user specified data. ANOVA is done in this case. */ public void doAnalysis() { if (dataTable.isEditing()) dataTable.getCellEditor().stopCellEditing(); if (!hasExample) { JOptionPane.showMessageDialog(this, DATA_MISSING_MESSAGE); return;/*from w ww. j a va 2 s . com*/ } Data data = new Data(); String firstMessage = "Would you like to use all the data columns in this Principal Components Analysis?"; String title = "SOCR - Principal Components Analysis"; String secondMessage = "Please enter the column numbers (seperated by comma) that you would like to use."; String columnNumbers = ""; int reply = JOptionPane.showConfirmDialog(null, firstMessage, title, JOptionPane.YES_NO_OPTION); if (reply == JOptionPane.YES_OPTION) { String cellValue = null; int originalRow = 0; for (int k = 0; k < dataTable.getRowCount(); k++) { cellValue = ((String) dataTable.getValueAt(k, 0)); if (cellValue != null && !cellValue.equals("")) { originalRow++; } } cellValue = null; int originalColumn = 0; for (int k = 0; k < dataTable.getColumnCount(); k++) { cellValue = ((String) dataTable.getValueAt(0, k)); if (cellValue != null && !cellValue.equals("")) { originalColumn++; } } dataRow = originalRow; dataColumn = originalColumn; String PCA_Data1[][] = new String[originalRow][originalColumn]; double PCA_Data[][] = new double[originalRow][originalColumn]; for (int k = 0; k < originalRow; k++) for (int j = 0; j < originalColumn; j++) { if (dataTable.getValueAt(k, j) != null && !dataTable.getValueAt(k, j).equals("")) { PCA_Data1[k][j] = (String) dataTable.getValueAt(k, j); PCA_Data[k][j] = Double.parseDouble(PCA_Data1[k][j]); } } double PCA_Adjusted_Data[][] = new double[originalRow][originalColumn]; double column_Total = 0; double column_Mean = 0; for (int j = 0; j < originalColumn; j++) for (int k = 0; k < originalRow; k++) { column_Total += PCA_Data[k][j]; if (k == (originalRow - 1)) { column_Mean = column_Total / originalRow; for (int p = 0; p < originalRow; p++) { PCA_Adjusted_Data[p][j] = PCA_Data[p][j] - column_Mean; } column_Total = 0; column_Mean = 0; } } Covariance cov = new Covariance(PCA_Adjusted_Data); RealMatrix matrix = cov.getCovarianceMatrix(); EigenDecomposition eigenDecomp = new EigenDecomposition(matrix, 0); storedData = eigenDecomp; RealMatrix eigenvectorMatrix = eigenDecomp.getV(); EValueArray = eigenDecomp.getRealEigenvalues(); eigenvectorMatrix = eigenvectorMatrix.transpose(); double eigenvectorArray[][] = eigenvectorMatrix.getData(); /*for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { System.out.println(eigenvectorArray[j][k] + " "); } */ Matrix matrix1 = new Matrix(eigenvectorArray); Matrix matrix2 = new Matrix(PCA_Adjusted_Data); matrix2 = matrix2.transpose(); Matrix finalProduct = matrix1.times(matrix2); finalProduct = finalProduct.transpose(); double finalArray[][] = finalProduct.getArrayCopy(); for (int j = 0; j < originalRow; j++) for (int k = 0; k < originalColumn; k++) { PCATable.setValueAt(finalArray[j][k], j, k); } xData = new double[originalRow]; yData = new double[originalRow]; for (int i = 0; i < originalRow; i++) { xData[i] = finalArray[i][0]; } for (int i = 0; i < originalRow; i++) { yData[i] = finalArray[i][1]; } dependentHeader = "C1"; independentHeader = "C2"; } else { // start here try { columnNumbers = JOptionPane.showInputDialog(secondMessage); } catch (Exception e) { } String columnNumbersFinal = "," + columnNumbers.replaceAll("\\s", "") + ","; Vector<Integer> locationOfComma = new Vector<Integer>(100); for (int i = 0; i < columnNumbersFinal.length(); i++) { char d = columnNumbersFinal.charAt(i); if (d == ',') locationOfComma.add(i); } Vector<Integer> vector = new Vector<Integer>(100); // vector containing column selected numbers for (int i = 0; i < locationOfComma.size() - 1; i++) { String temp = columnNumbersFinal.substring(locationOfComma.get(i) + 1, locationOfComma.get(i + 1)); if (temp == "") continue; vector.add((Integer.parseInt(temp) - 1)); } dependentHeader = "C" + (vector.get(0) + 1); independentHeader = "C" + (vector.get(1) + 1); // System.out.println("Vector size is: " + vector.size() + "\n"); String cellValue = null; int originalRow = 0; for (int k = 0; k < dataTable.getRowCount(); k++) { cellValue = ((String) dataTable.getValueAt(k, 0)); if (cellValue != null && !cellValue.equals("")) { originalRow++; } } int originalColumn = vector.size(); dataRow = originalRow; dataColumn = originalColumn; String PCA_Data1[][] = new String[originalRow][originalColumn]; double PCA_Data[][] = new double[originalRow][originalColumn]; for (int k = 0; k < originalRow; k++) for (int j = 0; j < originalColumn; j++) { if (dataTable.getValueAt(k, vector.get(j)) != null && !dataTable.getValueAt(k, vector.get(j)).equals("")) { PCA_Data1[k][j] = (String) dataTable.getValueAt(k, vector.get(j)); PCA_Data[k][j] = Double.parseDouble(PCA_Data1[k][j]); } } double PCA_Adjusted_Data[][] = new double[originalRow][originalColumn]; double column_Total = 0; double column_Mean = 0; for (int j = 0; j < originalColumn; j++) for (int k = 0; k < originalRow; k++) { column_Total += PCA_Data[k][j]; if (k == (originalRow - 1)) { column_Mean = column_Total / originalRow; for (int p = 0; p < originalRow; p++) { PCA_Adjusted_Data[p][j] = PCA_Data[p][j] - column_Mean; } column_Total = 0; column_Mean = 0; } } Covariance cov = new Covariance(PCA_Adjusted_Data); RealMatrix matrix = cov.getCovarianceMatrix(); EigenDecomposition eigenDecomp = new EigenDecomposition(matrix, 0); storedData = eigenDecomp; RealMatrix eigenvectorMatrix = eigenDecomp.getV(); EValueArray = eigenDecomp.getRealEigenvalues(); // added eigenvectorMatrix = eigenvectorMatrix.transpose(); double eigenvectorArray[][] = eigenvectorMatrix.getData(); /*for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { System.out.println(eigenvectorArray[j][k] + " "); } */ Matrix matrix1 = new Matrix(eigenvectorArray); Matrix matrix2 = new Matrix(PCA_Adjusted_Data); matrix2 = matrix2.transpose(); Matrix finalProduct = matrix1.times(matrix2); finalProduct = finalProduct.transpose(); double finalArray[][] = finalProduct.getArrayCopy(); /* for (int j = 0; j < dataTable.getColumnCount(); j++) for (int k = 0; k < dataTable.getRowCount(); k++) { System.out.println(finalArray[j][k] + " "); }*/ for (int j = 0; j < originalRow; j++) for (int k = 0; k < originalColumn; k++) { PCATable.setValueAt(finalArray[j][k], j, k); } xData = new double[originalRow]; yData = new double[originalRow]; for (int i = 0; i < originalRow; i++) { xData[i] = finalArray[i][0]; } for (int i = 0; i < originalRow; i++) { yData[i] = finalArray[i][1]; } } map = new HashMap<Double, double[]>(); for (int i = 0; i < dataColumn; i++) { map.put(EValueArray[i], storedData.getEigenvector(i).toArray()); } Arrays.sort(EValueArray); xData1 = new double[EValueArray.length]; // for Scree Plot yData1 = new double[EValueArray.length]; for (int i = 0; i < EValueArray.length; i++) { xData1[i] = i + 1; } for (int i = 0; i < EValueArray.length; i++) { yData1[i] = EValueArray[i]; } for (int i = 0; i < yData1.length / 2; i++) { double temp = yData1[i]; yData1[i] = yData1[yData1.length - i - 1]; yData1[yData1.length - i - 1] = temp; } for (int i = 0; i < xData1.length; i++) { System.out.println("xData1 contains: " + xData1[i] + "\n"); } for (int i = 0; i < yData1.length; i++) { System.out.println("yData1 contains: " + yData1[i] + "\n"); } for (int i = 0; i < EValueArray.length / 2; i++) { double temp = EValueArray[i]; EValueArray[i] = EValueArray[EValueArray.length - i - 1]; EValueArray[EValueArray.length - i - 1] = temp; } resultPanelTextArea.append( "Click on \"PCA RESULT\" panel to view the transformed data (Eigenvector Transposed * Adjusted Data Transposed)"); resultPanelTextArea.append("\n\nThe real eigenvalues (in descending order) are: \n\n"); resultPanelTextArea.append("" + round(EValueArray[0], 3)); for (int i = 1; i < EValueArray.length; i++) { resultPanelTextArea.append("\n" + round(EValueArray[i], 3)); } resultPanelTextArea.append("\n\nThe corresponding eigenvectors (in columns) are: \n\n"); double temp[] = new double[100]; for (int j = 0; j < temp.length; j++) for (int i = 0; i < EValueArray.length; i++) { temp = map.get(EValueArray[i]); resultPanelTextArea.append("" + round(temp[j], 3) + "\t"); if (i == EValueArray.length - 1) { resultPanelTextArea.append("\n"); } } doGraph(); }
From source file:com.trickl.math.lanczos.LanczosSolver.java
private void findM1M2(double[] eigenvalues) { int m2counter = 0; int n = 1;// ww w . j ava 2 s .co m int pos = 0; M1 = new int[eigenvalues.length]; Arrays.fill(M1, 0); M2 = new int[eigenvalues.length]; Arrays.fill(M2, 0); while (m2counter < (eigenvalues.length) && (n < alpha.length)) { n++; // n++ == 2, at first time in this loop. EigenDecomposition eigenDecomposition = new EigenDecomposition(alpha, beta); double[] eval = eigenDecomposition.getRealEigenvalues(); Arrays.sort(eval); // Consistent with IETL int M1_itr = 0; int M2_itr = 0; while (pos != eigenvalues.length) { if (M1[M1_itr] == 0 || M2[M2_itr] == 0) { double eigenvalue = eigenvalues[pos]; int ub = Sorting.LowerBound(eval, eigenvalue + threshold); int lb = Sorting.UpperBound(eval, eigenvalue - threshold); if (M1[M1_itr] == 0 && ub - lb >= 1) { M1[M1_itr] = n; } if (M2[M2_itr] == 0 && ub - lb >= 2) { M2[M2_itr] = n; m2counter++; } } pos++; M1_itr++; M2_itr++; } } }
From source file:ffx.autoparm.Energy.java
/** * <p>/*from ww w. ja v a2s. c o m*/ * system_mpoles</p> */ public void system_mpoles() { //Find center of mass. double weigh = 0; double xyzmid[] = { 0, 0, 0 }; double xyzcm[][] = new double[nAtoms][3]; for (int i = 0; i < nAtoms; i++) { weigh = weigh + atoms[i].getMass(); double coords[] = atoms[i].getXYZ(null); for (int j = 0; j < 3; j++) { xyzmid[j] = xyzmid[j] + coords[j] * atoms[i].getMass(); } } if (weigh != 0) { for (int j = 0; j < 3; j++) { xyzmid[j] = xyzmid[j] / weigh; } } for (int i = 0; i < nAtoms; i++) { for (int j = 0; j < 3; j++) { double coords[] = atoms[i].getXYZ(null); xyzcm[i][j] = coords[j] - xyzmid[j]; } } addInducedToGlobal(); double netchg = 0, xdpl = 0, ydpl = 0, zdpl = 0, xxqdp = 0, xyqdp = 0, xzqdp = 0, yxqdp = 0, yyqdp = 0, yzqdp = 0, zxqdp = 0, zyqdp = 0, zzqdp = 0; for (int i = 0; i < nAtoms; i++) { double charge = atoms[i].getMultipoleType().charge; double[] dipole = { pme2.globalMultipole[0][i][1], pme2.globalMultipole[0][i][2], pme2.globalMultipole[0][i][3] }; //double[] dipole = atoms[i].getMultipoleType().dipole; netchg = netchg + charge; xdpl = xdpl + xyzcm[i][0] * charge + dipole[0]; ydpl = ydpl + xyzcm[i][1] * charge + dipole[1]; zdpl = zdpl + xyzcm[i][2] * charge + dipole[2]; xxqdp = xxqdp + xyzcm[i][0] * xyzcm[i][0] * charge + 2 * xyzcm[i][0] * dipole[0]; xyqdp = xyqdp + xyzcm[i][0] * xyzcm[i][1] * charge + xyzcm[i][0] * dipole[1] + xyzcm[i][1] * dipole[0]; xzqdp = xzqdp + xyzcm[i][0] * xyzcm[i][2] * charge + xyzcm[i][0] * dipole[2] + xyzcm[i][2] * dipole[0]; yxqdp = yxqdp + xyzcm[i][0] * xyzcm[i][1] * charge + xyzcm[i][0] * dipole[1] + xyzcm[i][1] * dipole[0]; yyqdp = yyqdp + xyzcm[i][1] * xyzcm[i][1] * charge + 2 * xyzcm[i][1] * dipole[1]; yzqdp = yzqdp + xyzcm[i][1] * xyzcm[i][2] * charge + xyzcm[i][1] * dipole[2] + xyzcm[i][2] * dipole[1]; //zxqdp = zxqdp + xyzcm[i][2] * xyzcm[i][0] * charge + xyzcm[i][2] * dipole[0] + xyzcm[i][0] * dipole[2]; zxqdp = zxqdp + xyzcm[i][0] * xyzcm[i][2] * charge + xyzcm[i][0] * dipole[2] + xyzcm[i][2] * dipole[0]; //zyqdp = zyqdp + xyzcm[i][2] * xyzcm[i][1] * charge + xyzcm[i][2] * dipole[1] + xyzcm[i][1] * dipole[2]; zyqdp = zyqdp + xyzcm[i][1] * xyzcm[i][2] * charge + xyzcm[i][1] * dipole[2] + xyzcm[i][2] * dipole[1]; zzqdp = zzqdp + xyzcm[i][2] * xyzcm[i][2] * charge + 2 * xyzcm[i][2] * dipole[2]; } double qave = (xxqdp + yyqdp + zzqdp) / 3; xxqdp = 1.5 * (xxqdp - qave); xyqdp = 1.5 * xyqdp; xzqdp = 1.5 * xzqdp; yxqdp = 1.5 * yxqdp; yyqdp = 1.5 * (yyqdp - qave); yzqdp = 1.5 * yzqdp; zxqdp = 1.5 * zxqdp; zyqdp = 1.5 * zyqdp; zzqdp = 1.5 * (zzqdp - qave); for (int i = 0; i < nAtoms; i++) { double[][] quadrupole = { { pme2.globalMultipole[0][i][4], pme2.globalMultipole[0][i][7], pme2.globalMultipole[0][i][8] }, { pme2.globalMultipole[0][i][7], pme2.globalMultipole[0][i][5], pme2.globalMultipole[0][i][9] }, { pme2.globalMultipole[0][i][8], pme2.globalMultipole[0][i][9], pme2.globalMultipole[0][i][6] } }; //double[][] quadrupole = atoms[i].getMultipoleType().quadrupole; xxqdp = xxqdp + 3 * quadrupole[0][0]; xyqdp = xyqdp + 3 * quadrupole[0][1]; xzqdp = xzqdp + 3 * quadrupole[0][2]; yxqdp = yxqdp + 3 * quadrupole[1][0]; yyqdp = yyqdp + 3 * quadrupole[1][1]; yzqdp = yzqdp + 3 * quadrupole[1][2]; zxqdp = zxqdp + 3 * quadrupole[2][0]; zyqdp = zyqdp + 3 * quadrupole[2][1]; zzqdp = zzqdp + 3 * quadrupole[2][2]; } xdpl = MultipoleType.DEBYE * xdpl; ydpl = MultipoleType.DEBYE * ydpl; zdpl = MultipoleType.DEBYE * zdpl; xxqdp = MultipoleType.DEBYE * xxqdp; xyqdp = MultipoleType.DEBYE * xyqdp; xzqdp = MultipoleType.DEBYE * xzqdp; yxqdp = MultipoleType.DEBYE * yxqdp; yyqdp = MultipoleType.DEBYE * yyqdp; yzqdp = MultipoleType.DEBYE * yzqdp; zxqdp = MultipoleType.DEBYE * zxqdp; zyqdp = MultipoleType.DEBYE * zyqdp; zzqdp = MultipoleType.DEBYE * zzqdp; double netdpl = Math.sqrt(xdpl * xdpl + ydpl * ydpl + zdpl * zdpl); RealMatrix a = new Array2DRowRealMatrix( new double[][] { { xxqdp, xyqdp, xzqdp }, { yxqdp, yyqdp, yzqdp }, { zxqdp, zyqdp, zzqdp } }); EigenDecomposition e = new EigenDecomposition(a, 1); a = e.getD(); double[] netqdp = { a.getColumn(0)[0], a.getColumn(1)[1], a.getColumn(2)[2] }; DecimalFormat myFormatter = new DecimalFormat(" ##########0.00000;-##########0.00000"); String output; output = String.format(" Total Electric Charge: %13s %s Electrons\n", " ", myFormatter.format(netchg)); System.out.println(output); output = String.format(" Dipole Moment Magnitude: %13s %s Debyes\n", " ", myFormatter.format(netdpl)); System.out.println(output); output = String.format(" Dipole X,Y,Z-Components: %13s %s %s %s\n", " ", myFormatter.format(xdpl), myFormatter.format(ydpl), myFormatter.format(zdpl)); System.out.println(output); output = String.format(" Quadrupole Moment Tensor:%13s %s %s %s", " ", myFormatter.format(xxqdp), myFormatter.format(xyqdp), myFormatter.format(xzqdp)); System.out.println(output); output = String.format(" (Buckinghams) %13s %s %s %s", " ", myFormatter.format(yxqdp), myFormatter.format(yyqdp), myFormatter.format(yzqdp)); System.out.println(output); output = String.format(" %13s %s %s %s\n", " ", myFormatter.format(zxqdp), myFormatter.format(zyqdp), myFormatter.format(zzqdp)); System.out.println(output); output = String.format("Principle Axes Quadrupole:%13s %s %s %s", " ", myFormatter.format(netqdp[2]), myFormatter.format(netqdp[1]), myFormatter.format(netqdp[0])); System.out.println(output); }
From source file:gamlss.smoothing.PB.java
/** * GCV smoothing method.//from w w w . j av a2 s .co m * @param lambda - smoothing parameter * @return fitted values */ private ArrayRealVector functionGCV(Double lambda) { //QR <-qr(sqrt(w)*X) //Rinv <- solve(qr.R(QR)) QRDecomposition qr = new QRDecomposition( MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(w), basisM)); rM = (BlockRealMatrix) qr.getR(); rM = rM.getSubMatrix(0, rM.getColumnDimension() - 1, 0, rM.getColumnDimension() - 1); rMinv = (BlockRealMatrix) new QRDecomposition(rM).getSolver().getInverse(); //wy <- sqrt(w)*y ArrayRealVector wy = MatrixFunctions.sqrtVec(w).ebeMultiply(y); //y.y <- sum(wy^2) double y_y = MatrixFunctions.sumV(wy.ebeMultiply(wy)); //S <- t(D)%*%D sM = penaltyM.transpose().multiply(penaltyM); //UDU <- eigen(t(Rinv)%*%S%*%Rinv) uduM = new BlockRealMatrix( new EigenDecomposition(rMinv.transpose().multiply(sM).multiply(rMinv), Controls.SPLIT_TOLERANCE) .getV().getData()); //yy <- t(UDU$vectors)%*%t(qr.Q(QR))%*%wy BlockRealMatrix qM = (BlockRealMatrix) qr.getQ(); //SingularValueDecomposition svd = new SingularValueDecomposition( //MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(w), basisM)); //BlockRealMatrix qM = new BlockRealMatrix(svd.getV().getData()); //.... to finish !!!!!!! MatrixFunctions.matrixPrint(qM); return null; }
From source file:org.graphstream.algorithm.Spectrum.java
public void compute() { if (graph == null) throw new NotInitializedException(this); int m = graph.getNodeCount(); RealMatrix a = new Array2DRowRealMatrix(m, m); Edge e;/*ww w. j a va 2 s . co m*/ for (int idx1 = 0; idx1 < m; idx1++) for (int idx2 = 0; idx2 < m; idx2++) { e = graph.getNode(idx1).getEdgeToward(idx2); a.setEntry(idx1, idx2, e != null ? 1 : 0); } decomposition = new EigenDecomposition(a, 0); }