List of usage examples for org.apache.commons.math3.linear EigenDecomposition getD
public RealMatrix getD()
From source file:ffx.autoparm.Energy.java
/** * <p>/*from w w w . j ava 2 s .c om*/ * 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:msi.gama.util.matrix.GamaFloatMatrix.java
@Override public IList<Double> getEigen(final IScope scope) throws GamaRuntimeException { final RealMatrix rm = toApacheMatrix(scope); final EigenDecomposition ed = new EigenDecomposition(rm); return fromApacheMatrixtoDiagList(scope, ed.getD()); }
From source file:org.briljantframework.array.netlib.NetlibLinearAlgebraRoutinesTest.java
@Test public void testCi() throws Exception { final int n = 3; // DoubleArray a = // bj.newDoubleVector(-1.01, 3.98, 3.30, 4.43, 7.31, 0.86, 0.53, 8.26, 4.96, -6.43, -4.60, // -7.04, -3.89, -7.66, -6.16, 3.31, 5.29, 8.20, -7.33, 2.47, -4.81, 3.55, -1.51, 6.18, // 5.58).reshape(n, n); // System.out.println(a); DoubleArray a = bj.newDoubleMatrix(new double[][] { { 0, 2, 3 }, { 2, 0, 2 }, { 3, 2, 0 } }); DoubleArray wr = bj.newDoubleArray(n); DoubleArray wi = bj.newDoubleArray(n); DoubleArray vl = bj.newDoubleArray(n, n); DoubleArray vr = bj.newDoubleArray(n, n); DoubleArray copy = a.copy();//from ww w . jav a 2s. c o m // linalg.geev('v', 'v', copy, wr, wi, vl, vr); linalg.syev('v', 'u', copy, wr); wr.sort((i, j) -> Double.compare(j, i)); System.out.println(copy); System.out.println(wr); ComplexArray v = bj.newComplexArray(n, n); for (int i = 0; i < n; i++) { if (Precision.equals(wi.get(i), 0, 1e-6)) { v.setColumn(i, vr.getColumn(i).asComplex()); } else { DoubleArray real = vr.getColumn(i); DoubleArray imag = vr.getColumn(i + 1); v.setColumn(i, toComplex(real, imag)); v.setColumn(i + 1, toComplex(real, imag.negate())); i++; } } // System.out.println(v); RealMatrix matrix = Matrices.asRealMatrix(a); EigenDecomposition d = new EigenDecomposition(matrix); System.out.println(Matrices.toArray(d.getD())); System.out.println(Matrices.toArray(d.getV())); // System.out.println(d.getEigenvector(0)); System.out.println(Arrays.toString(d.getRealEigenvalues())); // System.out.println(Arrays.toString(d.getImagEigenvalues())); }
From source file:org.ojalgo.benchmark.lab.library.ACM.java
@Override public DecompositionOperation<RealMatrix, RealMatrix> getOperationEvD(final int dim) { final RealMatrix[] ret = this.makeArray(3); return (matrix) -> { final EigenDecomposition svd = new EigenDecomposition(matrix); ret[0] = svd.getV();/* ww w.ja v a 2 s . c om*/ ret[1] = svd.getD(); ret[2] = svd.getVT(); return ret; }; }
From source file:org.ujmp.commonsmath.AbstractCommonsMathDenseDoubleMatrix2D.java
public Matrix[] eig() { EigenDecomposition evd = new EigenDecomposition(matrix); Matrix v = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(evd.getV()); Matrix d = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(evd.getD()); return new Matrix[] { v, d }; }
From source file:playground.sergioo.facilitiesGenerator2012.WorkFacilitiesGeneration.java
private static Set<PointPerson> getPCATransformation(Collection<PointPerson> points) { RealMatrix pointsM = new Array2DRowRealMatrix(points.iterator().next().getDimension(), points.size()); int k = 0;/* w ww.j ava 2 s . co m*/ for (PointND<Double> point : points) { for (int f = 0; f < point.getDimension(); f++) pointsM.setEntry(f, k, point.getElement(f)); k++; } RealMatrix means = new Array2DRowRealMatrix(pointsM.getRowDimension(), 1); for (int r = 0; r < means.getRowDimension(); r++) { double mean = 0; for (int c = 0; c < pointsM.getColumnDimension(); c++) mean += pointsM.getEntry(r, c) / pointsM.getColumnDimension(); means.setEntry(r, 0, mean); } RealMatrix deviations = new Array2DRowRealMatrix(pointsM.getRowDimension(), pointsM.getColumnDimension()); for (int r = 0; r < deviations.getRowDimension(); r++) for (int c = 0; c < deviations.getColumnDimension(); c++) deviations.setEntry(r, c, pointsM.getEntry(r, c) - means.getEntry(r, 0)); RealMatrix covariance = deviations.multiply(deviations.transpose()) .scalarMultiply(1 / (double) pointsM.getColumnDimension()); EigenDecomposition eigenDecomposition = new EigenDecomposition(covariance, 0); RealMatrix eigenVectorsT = eigenDecomposition.getVT(); RealVector eigenValues = new ArrayRealVector(eigenDecomposition.getD().getRowDimension()); for (int r = 0; r < eigenDecomposition.getD().getRowDimension(); r++) eigenValues.setEntry(r, eigenDecomposition.getD().getEntry(r, r)); for (int i = 0; i < eigenValues.getDimension(); i++) { for (int j = i + 1; j < eigenValues.getDimension(); j++) if (eigenValues.getEntry(i) < eigenValues.getEntry(j)) { double tempValue = eigenValues.getEntry(i); eigenValues.setEntry(i, eigenValues.getEntry(j)); eigenValues.setEntry(j, tempValue); RealVector tempVector = eigenVectorsT.getRowVector(i); eigenVectorsT.setRowVector(i, eigenVectorsT.getRowVector(j)); eigenVectorsT.setRowVector(j, tempVector); } eigenVectorsT.setRowVector(i, eigenVectorsT.getRowVector(i).mapMultiply(Math.sqrt(1 / eigenValues.getEntry(i)))); } RealVector standardDeviations = new ArrayRealVector(pointsM.getRowDimension()); for (int r = 0; r < covariance.getRowDimension(); r++) standardDeviations.setEntry(r, Math.sqrt(covariance.getEntry(r, r))); double zValue = standardDeviations.dotProduct(new ArrayRealVector(pointsM.getRowDimension(), 1)); RealMatrix zScore = deviations.scalarMultiply(1 / zValue); pointsM = eigenVectorsT.multiply(zScore); Set<PointPerson> pointsC = new HashSet<PointPerson>(); k = 0; for (PointPerson point : points) { PointPerson pointC = new PointPerson(point.getId(), point.getOccupation(), new Double[] { pointsM.getEntry(0, k), pointsM.getEntry(1, k) }, point.getPlaceType()); pointC.setWeight(point.getWeight()); pointsC.add(pointC); k++; } return pointsC; }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * Update B and D from C./*from w w w . j av a 2 s . c o m*/ * * @param negccov Negative covariance factor. */ private void updateBD(double negccov) { if (ccov1 + ccovmu + negccov > 0 && (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1) { // to achieve O(N^2) C = triu(C, 0).add(triu(C, 1).transpose()); // enforce symmetry to prevent complex numbers final EigenDecomposition eig = new EigenDecomposition(C); B = eig.getV(); // eigen decomposition, B==normalized eigenvectors D = eig.getD(); diagD = diag(D); if (min(diagD) <= 0) { for (int i = 0; i < dimension; i++) { if (diagD.getEntry(i, 0) < 0) { diagD.setEntry(i, 0, 0); } } final double tfac = max(diagD) / 1e14; C = C.add(eye(dimension, dimension).scalarMultiply(tfac)); diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac)); } if (max(diagD) > 1e14 * min(diagD)) { final double tfac = max(diagD) / 1e14 - min(diagD); C = C.add(eye(dimension, dimension).scalarMultiply(tfac)); diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac)); } diagC = diag(C); diagD = sqrt(diagD); // D contains standard deviations now BD = times(B, repmat(diagD.transpose(), dimension, 1)); // O(n^2) } }