List of usage examples for org.apache.commons.math3.linear RealMatrix preMultiply
RealVector preMultiply(RealVector v) throws DimensionMismatchException;
From source file:ffx.crystal.Crystal.java
/** * Update all Crystal variables that are a function of unit cell parameters. */// ww w. j a va 2s .c o m private void updateCrystal() { switch (crystalSystem) { case CUBIC: case ORTHORHOMBIC: case TETRAGONAL: cos_alpha = 0.0; sin_beta = 1.0; cos_beta = 0.0; sin_gamma = 1.0; cos_gamma = 0.0; beta_term = 0.0; gamma_term = 1.0; volume = a * b * c; dVdA = b * c; dVdB = a * c; dVdC = a * b; dVdAlpha = 0.0; dVdBeta = 0.0; dVdGamma = 0.0; break; case MONOCLINIC: cos_alpha = 0.0; sin_beta = sin(toRadians(beta)); cos_beta = cos(toRadians(beta)); sin_gamma = 1.0; cos_gamma = 0.0; beta_term = 0.0; gamma_term = sin_beta; volume = sin_beta * a * b * c; dVdA = sin_beta * b * c; dVdB = sin_beta * a * c; dVdC = sin_beta * a * b; dVdAlpha = 0.0; dVdBeta = cos_beta * a * b * c; dVdGamma = 0.0; break; case HEXAGONAL: case TRICLINIC: case TRIGONAL: default: double sin_alpha = sin(toRadians(alpha)); cos_alpha = cos(toRadians(alpha)); sin_beta = sin(toRadians(beta)); cos_beta = cos(toRadians(beta)); sin_gamma = sin(toRadians(gamma)); cos_gamma = cos(toRadians(gamma)); beta_term = (cos_alpha - cos_beta * cos_gamma) / sin_gamma; gamma_term = sqrt(sin_beta * sin_beta - beta_term * beta_term); volume = sin_gamma * gamma_term * a * b * c; dVdA = sin_gamma * gamma_term * b * c; dVdB = sin_gamma * gamma_term * a * c; dVdC = sin_gamma * gamma_term * a * b; double dbeta = 2.0 * sin_beta * cos_beta - (2.0 * (cos_alpha - cos_beta * cos_gamma) * sin_beta * cos_gamma) / (sin_gamma * sin_gamma); double dgamma1 = -2.0 * (cos_alpha - cos_beta * cos_gamma) * cos_beta / sin_gamma; double dgamma2 = cos_alpha - cos_beta * cos_gamma; dgamma2 *= dgamma2 * 2.0 * cos_gamma / (sin_gamma * sin_gamma * sin_gamma); dVdAlpha = (cos_alpha - cos_beta * cos_gamma) * sin_alpha / (sin_gamma * gamma_term) * a * b * c; dVdBeta = 0.5 * sin_gamma * dbeta / gamma_term * a * b * c; dVdGamma = (cos_gamma * gamma_term + 0.5 * sin_gamma * (dgamma1 + dgamma2) / gamma_term) * a * b * c; break; } // a is the first row of A^(-1). Ai[0][0] = a; Ai[0][1] = 0.0; Ai[0][2] = 0.0; // b is the second row of A^(-1). Ai[1][0] = b * cos_gamma; Ai[1][1] = b * sin_gamma; Ai[1][2] = 0.0; // c is the third row of A^(-1). Ai[2][0] = c * cos_beta; Ai[2][1] = c * beta_term; Ai[2][2] = c * gamma_term; Ai00 = Ai[0][0]; Ai01 = Ai[0][1]; Ai02 = Ai[0][2]; Ai10 = Ai[1][0]; Ai11 = Ai[1][1]; Ai12 = Ai[1][2]; Ai20 = Ai[2][0]; Ai21 = Ai[2][1]; Ai22 = Ai[2][2]; // Invert A^-1 to get A RealMatrix m = new Array2DRowRealMatrix(Ai, true); m = new LUDecomposition(m).getSolver().getInverse(); A = m.getData(); // The columns of A are the reciprocal basis vectors A00 = A[0][0]; A10 = A[1][0]; A20 = A[2][0]; A01 = A[0][1]; A11 = A[1][1]; A21 = A[2][1]; A02 = A[0][2]; A12 = A[1][2]; A22 = A[2][2]; // Reciprocal basis vector lengths aStar = 1.0 / sqrt(A00 * A00 + A10 * A10 + A20 * A20); bStar = 1.0 / sqrt(A01 * A01 + A11 * A11 + A21 * A21); cStar = 1.0 / sqrt(A02 * A02 + A12 * A12 + A22 * A22); if (logger.isLoggable(Level.FINEST)) { logger.finest(String.format(" Reciprocal Lattice Lengths: (%8.3f, %8.3f, %8.3f)", aStar, bStar, cStar)); } // Interfacial diameters from the dot product of the real and reciprocal vectors interfacialRadiusA = (Ai00 * A00 + Ai01 * A10 + Ai02 * A20) * aStar; interfacialRadiusB = (Ai10 * A01 + Ai11 * A11 + Ai12 * A21) * bStar; interfacialRadiusC = (Ai20 * A02 + Ai21 * A12 + Ai22 * A22) * cStar; // Divide by 2 to get radii. interfacialRadiusA /= 2.0; interfacialRadiusB /= 2.0; interfacialRadiusC /= 2.0; if (logger.isLoggable(Level.FINEST)) { logger.finest(String.format(" Interfacial radii: (%8.3f, %8.3f, %8.3f)", interfacialRadiusA, interfacialRadiusB, interfacialRadiusC)); } G[0][0] = a * a; G[0][1] = a * b * cos_gamma; G[0][2] = a * c * cos_beta; G[1][0] = G[0][1]; G[1][1] = b * b; G[1][2] = b * c * cos_alpha; G[2][0] = G[0][2]; G[2][1] = G[1][2]; G[2][2] = c * c; // invert G to yield Gstar m = new Array2DRowRealMatrix(G, true); m = new LUDecomposition(m).getSolver().getInverse(); Gstar = m.getData(); List<SymOp> symOps = spaceGroup.symOps; int nSymm = symOps.size(); if (symOpsCartesian == null) { symOpsCartesian = new ArrayList<>(nSymm); } else { symOpsCartesian.clear(); } RealMatrix toFrac = new Array2DRowRealMatrix(A); RealMatrix toCart = new Array2DRowRealMatrix(Ai); for (int i = 0; i < nSymm; i++) { SymOp symOp = symOps.get(i); m = new Array2DRowRealMatrix(symOp.rot); // rot_c = A^(-1).rot_f.A RealMatrix rotMat = m.preMultiply(toCart.transpose()).multiply(toFrac.transpose()); // tr_c = tr_f.A^(-1) double tr[] = toCart.preMultiply(symOp.tr); symOpsCartesian.add(new SymOp(rotMat.getData(), tr)); } }
From source file:uk.ac.cam.cl.dtg.teaching.programmingtest.java.FittedCurve.java
private double fit(List<Point> mapped, double proportion) { Collections.shuffle(mapped);//from w ww .j a v a 2s . co m int total = (int) (mapped.size() * proportion); double[][] x = new double[total][]; double[] y = new double[total]; Iterator<Point> step = mapped.iterator(); for (int i = 0; i < total; ++i) { Point next = step.next(); x[i] = next.xcoeff; y[i] = next.yvalue; } OLSMultipleLinearRegression ols = new OLSMultipleLinearRegression(); ols.setNoIntercept(true); ols.newSampleData(y, x); RealMatrix coef = MatrixUtils.createColumnRealMatrix(ols.estimateRegressionParameters()); double rsqTotal = 0.0; for (Point p : mapped) { double yhat = coef.preMultiply(p.xcoeff)[0]; rsqTotal += unmapY(Math.pow(p.yvalue - yhat, 2.0)); } if (DEBUG) { List<Point> newMaps = mapValues(this.xs, this.ys); LinkedList<Double> xs = new LinkedList<>(); LinkedList<Double> ys = new LinkedList<>(); int i = 0; for (Point p : newMaps) { double yhat = unmapY(coef.preMultiply(p.xcoeff)[0]); ys.addLast(yhat); xs.addLast(this.xs[i++]); } System.out.println("x=[" + Joiner.on(",").join(xs) + "]"); System.out.println("y=[" + Joiner.on(",").join(ys) + "]"); } return rsqTotal; }