List of usage examples for org.apache.commons.math3.linear RealMatrix transpose
RealMatrix transpose();
From source file:com.itemanalysis.psychometrics.factoranalysis.GPArotation.java
/** * Conducts orthogonal rotation of factor loadings. * * @param A matrix of orthogonal factor loadings * @return a matrix of rotated factor loadings. * @throws ConvergenceException//from ww w. ja v a 2 s. c om */ private RotationResults GPForth(RealMatrix A, boolean normalize, int maxIter, double eps) throws ConvergenceException { int ncol = A.getColumnDimension(); if (normalize) { //elementwise division by normalizing weights final RealMatrix W = getNormalizingWeights(A, true); A.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return value / W.getEntry(row, column); } }); } RealMatrix Tmat = new IdentityMatrix(ncol); double alpha = 1; RealMatrix L = A.multiply(Tmat); gpFunction.computeValues(L); double f = gpFunction.getValue(); RealMatrix VgQ = gpFunction.getGradient(); RealMatrix G = A.transpose().multiply(VgQ); double VgQtF = gpFunction.getValue(); RealMatrix VgQt = gpFunction.getGradient(); RealMatrix Tmatt = null; int iter = 0; double s = eps + 0.5; double s2 = 0; int innnerIter = 11; while (iter < maxIter) { RealMatrix M = Tmat.transpose().multiply(G); RealMatrix S = (M.add(M.transpose())); S = S.scalarMultiply(0.5); RealMatrix Gp = G.subtract(Tmat.multiply(S)); s = Math.sqrt((Gp.transpose().multiply(Gp)).getTrace()); s2 = Math.pow(s, 2); if (s < eps) break; alpha *= 2.0; for (int j = 0; j < innnerIter; j++) { Gp = Gp.scalarMultiply(alpha); RealMatrix X = (Tmat.subtract(Gp)); SingularValueDecomposition SVD = new SingularValueDecomposition(X); Tmatt = SVD.getU().multiply(SVD.getV().transpose()); L = A.multiply(Tmatt); gpFunction.computeValues(L); VgQt = gpFunction.getGradient(); VgQtF = gpFunction.getValue(); if (VgQtF < f - 0.5 * s2 * alpha) { break; } alpha /= 2.0; } Tmat = Tmatt; f = VgQtF; G = A.transpose().multiply(VgQt); iter++; } boolean convergence = s < eps; if (!convergence) { throw new ConvergenceException(); } if (normalize) { //elementwise multiplication by normalizing weights final RealMatrix W = getNormalizingWeights(A, true); A.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return value * W.getEntry(row, column); } }); } RealMatrix Phi = Tmat.transpose().multiply(Tmat); RotationResults result = new RotationResults(gpFunction.getValue(), L, Phi, Tmat, rotationMethod); return result; }
From source file:lirmm.inria.fr.math.BigSparseRealMatrixTest.java
/** * test transpose//from w w w . ja v a2 s . c om */ @Test public void testTranspose() { RealMatrix m = new BigSparseRealMatrix(testData); RealMatrix mIT = new LUDecomposition(m).getSolver().getInverse().transpose(); RealMatrix mTI = new LUDecomposition(m.transpose()).getSolver().getInverse(); TestUtils.assertEquals("inverse-transpose", mIT, mTI, normTolerance); m = new BigSparseRealMatrix(testData2); RealMatrix mt = new BigSparseRealMatrix(testData2T); TestUtils.assertEquals("transpose", mt, m.transpose(), normTolerance); }
From source file:eagle.security.userprofile.impl.UserProfileAnomalyEigenEvaluator.java
@Override public List<MLCallbackResult> detect(final String user, final String algorithm, UserActivityAggModel userActivity, UserProfileEigenModel aModel) { RealMatrix inputData = userActivity.matrix(); LOG.warn("EigenBasedAnomalyDetection predictAnomaly called with dimension: " + inputData.getRowDimension() + "x" + inputData.getColumnDimension()); if (aModel == null) { LOG.warn(/*from w w w . j a v a 2s. c o m*/ "nothing to do as the input model does not have required values, returning from evaluating this algorithm.."); return null; } List<MLCallbackResult> mlCallbackResults = new ArrayList<MLCallbackResult>(); RealMatrix normalizedMat = normalizeData(inputData, aModel); UserCommandStatistics[] listStats = aModel.statistics(); int colWithHighVariant = 0; for (int j = 0; j < normalizedMat.getColumnDimension(); j++) { if (listStats[j].isLowVariant() == false) { colWithHighVariant++; } } final Map<String, String> context = new HashMap<String, String>() { { put(UserProfileConstants.USER_TAG, user); put(UserProfileConstants.ALGORITHM_TAG, algorithm); } }; Map<Integer, String> lineNoWithVariantBasedAnomalyDetection = new HashMap<Integer, String>(); for (int i = 0; i < normalizedMat.getRowDimension(); i++) { MLCallbackResult aResult = new MLCallbackResult(); aResult.setAnomaly(false); aResult.setContext(context); for (int j = 0; j < normalizedMat.getColumnDimension(); j++) { //LOG.info("mean for j=" + j + " is:" + listStats[j].getMean()); //LOG.info("stddev for j=" + j + " is:" + listStats[j].getStddev()); if (listStats[j].isLowVariant() == true) { //LOG.info(listOfCmds[j] + " is low variant"); if (normalizedMat.getEntry(i, j) > listStats[j].getMean()) { lineNoWithVariantBasedAnomalyDetection.put(i, "lowVariantAnomaly"); aResult.setAnomaly(true); aResult.setTimestamp(userActivity.timestamp()); aResult.setFeature(listStats[j].getCommandName()); aResult.setAlgorithm(UserProfileConstants.EIGEN_DECOMPOSITION_ALGORITHM); List<String> datapoints = new ArrayList<String>(); double[] rowVals = inputData.getRow(i); for (double rowVal : rowVals) datapoints.add(rowVal + ""); aResult.setDatapoints(datapoints); aResult.setId(user); //mlCallbackResults.add(aResult); } /*else { aResult.setAnomaly(false); aResult.setTimestamp(userActivity.timestamp()); mlCallbackResults.add(aResult); }*/ } } mlCallbackResults.add(i, aResult); //return results; } //LOG.info("results size here: " + results.length); //LOG.info("col with high variant: " + colWithHighVariant); RealMatrix finalMatWithoutLowVariantFeatures = new Array2DRowRealMatrix(normalizedMat.getRowDimension(), colWithHighVariant); //LOG.info("size of final test data: " + finalMatWithoutLowVariantFeatures.getRowDimension() +"x"+ finalMatWithoutLowVariantFeatures.getColumnDimension()); int finalMatrixRow = 0; int finalMatrixCol = 0; for (int i = 0; i < normalizedMat.getRowDimension(); i++) { for (int j = 0; j < normalizedMat.getColumnDimension(); j++) { if (listStats[j].isLowVariant() == false) { finalMatWithoutLowVariantFeatures.setEntry(finalMatrixRow, finalMatrixCol, normalizedMat.getEntry(i, j)); finalMatrixCol++; } } finalMatrixCol = 0; finalMatrixRow++; } RealVector[] pcs = aModel.principalComponents(); //LOG.info("pc size: " + pcs.getRowDimension() +"x" + pcs.getColumnDimension()); RealMatrix finalInputMatTranspose = finalMatWithoutLowVariantFeatures.transpose(); for (int i = 0; i < finalMatWithoutLowVariantFeatures.getRowDimension(); i++) { if (lineNoWithVariantBasedAnomalyDetection.get(i) == null) { //MLCallbackResult result = new MLCallbackResult(); MLCallbackResult result = mlCallbackResults.get(i); //result.setContext(context); for (int sz = 0; sz < pcs.length; sz++) { double[] pc1 = pcs[sz].toArray(); RealMatrix pc1Mat = new Array2DRowRealMatrix(pc1); RealMatrix transposePC1Mat = pc1Mat.transpose(); RealMatrix testData = pc1Mat.multiply(transposePC1Mat) .multiply(finalInputMatTranspose.getColumnMatrix(i)); //LOG.info("testData size: " + testData.getRowDimension() + "x" + testData.getColumnDimension()); RealMatrix testDataTranspose = testData.transpose(); //LOG.info("testData transpose size: " + testDataTranspose.getRowDimension() + "x" + testDataTranspose.getColumnDimension()); RealVector iRowVector = testDataTranspose.getRowVector(0); //RealVector pc1Vector = transposePC1Mat.getRowVector(sz); RealVector pc1Vector = transposePC1Mat.getRowVector(0); double distanceiRowAndPC1 = iRowVector.getDistance(pc1Vector); //LOG.info("distance from pc sz: " + sz + " " + distanceiRowAndPC1 + " " + model.getMaxL2Norm().getEntry(sz)); //LOG.info("model.getMaxL2Norm().getEntry(sz):" + model.getMaxL2Norm().getEntry(sz)); if (distanceiRowAndPC1 > aModel.maximumL2Norm().getEntry(sz)) { //LOG.info("distance from pc sz: " + sz + " " + distanceiRowAndPC1 + " " + model.getMaxL2Norm().getEntry(sz)); result.setAnomaly(true); result.setFeature(aModel.statistics()[sz].getCommandName()); result.setTimestamp(System.currentTimeMillis()); result.setAlgorithm(UserProfileConstants.EIGEN_DECOMPOSITION_ALGORITHM); List<String> datapoints = new ArrayList<String>(); double[] rowVals = inputData.getRow(i); for (double rowVal : rowVals) datapoints.add(rowVal + ""); result.setDatapoints(datapoints); result.setId(user); } } mlCallbackResults.set(i, result); } } return mlCallbackResults; }
From source file:com.google.location.lbs.gnss.gps.pseudorange.UserPositionVelocityWeightedLeastSquare.java
/** * Calculates the position uncertainty in meters and the velocity uncertainty * in meters per second solution in local ENU system. * * <p> Reference: Global Positioning System: Signals, Measurements, and Performance * by Pratap Misra, Per Enge, Page 206 - 209. * * @param velocityWeightMatrix the velocity weight matrix * @param positionWeightMatrix the position weight matrix * @param positionVelocitySolution the position and velocity solution in ECEF * @return an array containing the position and velocity uncertainties in ENU coordinate system. * [0-2] Enu uncertainty of position solution in meters. * [3-5] Enu uncertainty of velocity solution in meters per second. *//* www. j av a 2s. c o m*/ public double[] calculatePositionVelocityUncertaintyEnu(RealMatrix velocityWeightMatrix, RealMatrix positionWeightMatrix, double[] positionVelocitySolution) { if (geometryMatrix == null) { return null; } RealMatrix velocityH = calculateHMatrix(velocityWeightMatrix, geometryMatrix); RealMatrix positionH = calculateHMatrix(positionWeightMatrix, geometryMatrix); // Calculate the rotation Matrix to convert to local ENU system. RealMatrix rotationMatrix = new Array2DRowRealMatrix(4, 4); GeodeticLlaValues llaValues = Ecef2LlaConverter.convertECEFToLLACloseForm(positionVelocitySolution[0], positionVelocitySolution[1], positionVelocitySolution[2]); rotationMatrix.setSubMatrix(Ecef2EnuConverter .getRotationMatrix(llaValues.longitudeRadians, llaValues.latitudeRadians).getData(), 0, 0); rotationMatrix.setEntry(3, 3, 1); // Convert to local ENU by pre-multiply rotation matrix and multiply rotation matrix transposed velocityH = rotationMatrix.multiply(velocityH).multiply(rotationMatrix.transpose()); positionH = rotationMatrix.multiply(positionH).multiply(rotationMatrix.transpose()); // Return the square root of diagonal entries return new double[] { Math.sqrt(positionH.getEntry(0, 0)), Math.sqrt(positionH.getEntry(1, 1)), Math.sqrt(positionH.getEntry(2, 2)), Math.sqrt(velocityH.getEntry(0, 0)), Math.sqrt(velocityH.getEntry(1, 1)), Math.sqrt(velocityH.getEntry(2, 2)) }; }
From source file:ellipsoidFit.FitPoints.java
/** * Solve the polynomial expression Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + * 2Gx + 2Hy + 2Iz from the provided points. * //from ww w. j ava2 s. com * @param points * the points that will be fit to the polynomial expression. * @return the solution vector to the polynomial expression. */ private RealVector solveSystem(ArrayList<ThreeSpacePoint> points) { // determine the number of points int numPoints = points.size(); // the design matrix // size: numPoints x 9 RealMatrix d = new Array2DRowRealMatrix(numPoints, 9); // Fit the ellipsoid in the form of // Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz for (int i = 0; i < d.getRowDimension(); i++) { double xx = Math.pow(points.get(i).x, 2); double yy = Math.pow(points.get(i).y, 2); double zz = Math.pow(points.get(i).z, 2); double xy = 2 * (points.get(i).x * points.get(i).y); double xz = 2 * (points.get(i).x * points.get(i).z); double yz = 2 * (points.get(i).y * points.get(i).z); double x = 2 * points.get(i).x; double y = 2 * points.get(i).y; double z = 2 * points.get(i).z; d.setEntry(i, 0, xx); d.setEntry(i, 1, yy); d.setEntry(i, 2, zz); d.setEntry(i, 3, xy); d.setEntry(i, 4, xz); d.setEntry(i, 5, yz); d.setEntry(i, 6, x); d.setEntry(i, 7, y); d.setEntry(i, 8, z); } // solve the normal system of equations // v = (( d' * d )^-1) * ( d' * ones.mapAddToSelf(1)); // Multiply: d' * d RealMatrix dtd = d.transpose().multiply(d); // Create a vector of ones. RealVector ones = new ArrayRealVector(numPoints); ones.mapAddToSelf(1); // Multiply: d' * ones.mapAddToSelf(1) RealVector dtOnes = d.transpose().operate(ones); // Find ( d' * d )^-1 DecompositionSolver solver = new SingularValueDecomposition(dtd).getSolver(); RealMatrix dtdi = solver.getInverse(); // v = (( d' * d )^-1) * ( d' * ones.mapAddToSelf(1)); RealVector v = dtdi.operate(dtOnes); return v; }
From source file:edu.cudenver.bios.power.GLMMPowerCalculator.java
/** * Perform any preliminary calculations / updates on the input matrices * @param params input parameters// w ww. ja va 2s .c o m */ private void initialize(GLMMPowerParameters params) { debug("entering initialize"); // TODO: why isn't this done in PowerResourceHelper.studyDesignToPowerParameters? // if no power methods are specified, add conditional as a default if (params.getPowerMethodList().size() <= 0) params.addPowerMethod(PowerMethod.CONDITIONAL_POWER); // update sigma error and beta if we have a baseline covariate RealMatrix sigmaG = params.getSigmaGaussianRandom(); debug("sigmaG", sigmaG); int numRandom = sigmaG != null ? sigmaG.getRowDimension() : 0; if (numRandom == 1) { // set the sigma error matrix to [sigmaY - sigmaYG * sigmaG-1 * sigmaGY] RealMatrix sigmaY = params.getSigmaOutcome(); debug("sigmaY", sigmaY); RealMatrix sigmaYG = params.getSigmaOutcomeGaussianRandom(); debug("sigmaYG", sigmaYG); RealMatrix sigmaGY = sigmaYG.transpose(); debug("sigmaYG transpose", sigmaGY); RealMatrix sigmaGInverse = new LUDecomposition(sigmaG).getSolver().getInverse(); debug("sigmaG inverse", sigmaGInverse); RealMatrix sigmaError = forceSymmetric( sigmaY.subtract(sigmaYG.multiply(sigmaGInverse.multiply(sigmaGY)))); debug("sigmaError = sigmaY - sigmaYG * sigmaG inverse * sigmaYG transpose", sigmaError); if (!MatrixUtils.isPositiveSemidefinite(sigmaError)) { throw new IllegalArgumentException(SIGMA_ERROR_NOT_POSITIVE_SEMIDEFINITE_MESSAGE); } params.setSigmaError(sigmaError); // calculate the betaG matrix and fill in the placeholder row for the random predictor FixedRandomMatrix beta = params.getBeta(); beta.updateRandomMatrix(sigmaGInverse.multiply(sigmaGY)); debug("beta with random part updated to sigmaG inverse * sigmaYG transpose", beta.getCombinedMatrix()); } debug("exiting initialize"); }
From source file:edu.cudenver.bios.power.GLMMPowerCalculator.java
/** * Generate a sample of powers for a design with a random covariate * @param params power parameters/*from w w w. j av a 2s .c o m*/ * @param iterations size of the sample * @return */ private SimulatedPower[] getSimulatedPowerSampleValue(GLMMPowerParameters params, Test test, PowerMethod method, double alpha, double sigmaScale, double betaScale, int sampleSize, double quantile, int iterations) { // scale beta and sigma RealMatrix scaledBeta = params.getBeta().scalarMultiply(betaScale, true); RealMatrix scaledSigmaError = params.getSigmaError().scalarMultiply(sigmaScale); // get random errors RandomErrorMatrix randomErrors = new RandomErrorMatrix( MatrixUtils.getTotalSampleSize(params.getDesignEssence(), sampleSize), scaledBeta.getColumnDimension(), scaledSigmaError); randomErrors.setPositivityThreshold(positivityThreshold); randomErrors.setSymmetryThreshold(symmetryThreshold); SimulatedPower[] simPower = new SimulatedPower[iterations]; for (int gInstance = 0; gInstance < iterations; gInstance++) { // force a new realization of the design matrix (i.e. a new covariate column) RealMatrix X = MatrixUtils.getFullDesignMatrix(params.getDesignEssence(), params.getSigmaGaussianRandom(), sampleSize, seed); RealMatrix XtXinverse = new LUDecomposition(X.transpose().multiply(X)).getSolver().getInverse(); int rejectionCount = 0; simPower[gInstance] = new SimulatedPower(scaledBeta.getRowDimension(), scaledBeta.getColumnDimension()); for (int i = 0; i < SIMULATION_ITERATIONS_QUANTILE_UNCONDITIONAL; i++) { SimulationFit fit = simulateAndFitModel(params, test, X, XtXinverse, scaledBeta, randomErrors, sampleSize, alpha); if (fit.Pvalue <= alpha) rejectionCount++; simPower[gInstance].accumulateBeta(fit.beta); } simPower[gInstance] .setPower(((double) rejectionCount) / ((double) SIMULATION_ITERATIONS_QUANTILE_UNCONDITIONAL)); } return simPower; }
From source file:edu.cudenver.bios.power.GLMMPowerCalculator.java
private GLMMPower getSimulatedPowerValue(GLMMPowerParameters params, Test test, PowerMethod method, double alpha, double sigmaScale, double betaScale, int sampleSize, double quantile, int iterations) { // scale beta and sigma RealMatrix scaledBeta = params.getBeta().scalarMultiply(betaScale, true); RealMatrix scaledSigmaError = params.getSigmaError().scalarMultiply(sigmaScale); // get random errors RandomErrorMatrix randomErrors = new RandomErrorMatrix( MatrixUtils.getTotalSampleSize(params.getDesignEssence(), sampleSize), scaledSigmaError.getColumnDimension(), scaledSigmaError); randomErrors.setPositivityThreshold(positivityThreshold); randomErrors.setSymmetryThreshold(symmetryThreshold); double power = Double.NaN; if (params.getSigmaGaussianRandom() != null && method != PowerMethod.CONDITIONAL_POWER) { double[] powerValues = new double[SIMULATION_ITERATIONS_QUANTILE_UNCONDITIONAL]; for (int gInstance = 0; gInstance < SIMULATION_ITERATIONS_QUANTILE_UNCONDITIONAL; gInstance++) { int rejectionCount = 0; // force a new realization of the design matrix (i.e. a new covariate column) RealMatrix X = MatrixUtils.getFullDesignMatrix(params.getDesignEssence(), params.getSigmaGaussianRandom(), sampleSize, seed); RealMatrix XtXinverse = new LUDecomposition(X.transpose().multiply(X)).getSolver().getInverse(); for (int i = 0; i < SIMULATION_ITERATIONS_QUANTILE_UNCONDITIONAL; i++) { SimulationFit fit = simulateAndFitModel(params, test, X, XtXinverse, scaledBeta, randomErrors, sampleSize, alpha); if (fit.Pvalue <= alpha) rejectionCount++;/*ww w . j a v a 2 s .co m*/ } powerValues[gInstance] = (((double) rejectionCount) / ((double) SIMULATION_ITERATIONS_QUANTILE_UNCONDITIONAL)); } switch (method) { case UNCONDITIONAL_POWER: power = StatUtils.mean(powerValues); break; case QUANTILE_POWER: power = StatUtils.percentile(powerValues, quantile); break; default: throw new IllegalArgumentException("Unknown power method"); } } else { // GLMM(F) design, conditional power // run the simulations RealMatrix X = MatrixUtils.getFullDesignMatrix(params.getDesignEssence(), sampleSize); RealMatrix XtXInverse = params.getXtXInverse().scalarMultiply(1 / (double) sampleSize); int rejectionCount = 0; for (int i = 0; i < iterations; i++) { SimulationFit fit = simulateAndFitModel(params, test, X, XtXInverse, scaledBeta, randomErrors, sampleSize, alpha); if (fit.Pvalue <= alpha) rejectionCount++; } power = ((double) rejectionCount) / ((double) iterations); } return new GLMMPower(test, alpha, power, power, MatrixUtils.getTotalSampleSize(params.getDesignEssence(), sampleSize), betaScale, sigmaScale, method, quantile, null); }
From source file:ffx.crystal.Crystal.java
/** * Update all Crystal variables that are a function of unit cell parameters. */// w w w . jav a2 s. co 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:edu.cudenver.bios.power.glmm.NonCentralityDistribution.java
/** * Pre-calculate intermediate matrices, perform setup, etc. */// w w w . j av a 2 s. com private void initialize(Test test, RealMatrix FEssence, RealMatrix FtFinverse, int perGroupN, FixedRandomMatrix CFixedRand, RealMatrix U, RealMatrix thetaNull, RealMatrix beta, RealMatrix sigmaError, RealMatrix sigmaG, boolean exact) throws PowerException { debug("entering initialize"); // reset member variables this.T1 = null; this.FT1 = null; this.S = null; this.mzSq = null; this.H0 = 0; this.sStar = 0; // cache inputs this.test = test; this.FEssence = FEssence; this.FtFinverse = FtFinverse; this.perGroupN = perGroupN; this.CFixedRand = CFixedRand; this.U = U; this.thetaNull = thetaNull; this.beta = beta; this.sigmaError = sigmaError; this.sigmaG = sigmaG; // calculate intermediate matrices // RealMatrix FEssence = params.getDesignEssence().getFullDesignMatrixFixed(); // TODO: do we ever get here with values that can cause integer overflow, // and if so, does it matter? this.N = (double) FEssence.getRowDimension() * perGroupN; this.exact = exact; try { // TODO: need to calculate H0, need to adjust H1 for Unirep // get design matrix for fixed parameters only qF = FEssence.getColumnDimension(); // a = CFixedRand.getCombinedMatrix().getRowDimension(); // get fixed contrasts RealMatrix Cfixed = CFixedRand.getFixedMatrix(); RealMatrix CGaussian = CFixedRand.getRandomMatrix(); // build intermediate terms h1, S if (FtFinverse == null) { FtFinverse = new LUDecomposition(FEssence.transpose().multiply(FEssence)).getSolver().getInverse(); debug("FEssence", FEssence); debug("FtFinverse = (FEssence transpose * FEssence) inverse", FtFinverse); } else { debug("FtFinverse", FtFinverse); } RealMatrix PPt = Cfixed.multiply(FtFinverse.scalarMultiply(1 / (double) perGroupN)) .multiply(Cfixed.transpose()); debug("Cfixed", Cfixed); debug("n = " + perGroupN); debug("PPt = Cfixed * FtF inverse * (1/n) * Cfixed transpose", PPt); T1 = forceSymmetric(new LUDecomposition(PPt).getSolver().getInverse()); debug("T1 = PPt inverse", T1); FT1 = new CholeskyDecomposition(T1).getL(); debug("FT1 = Cholesky decomposition (L) of T1", FT1); // calculate theta difference // RealMatrix thetaNull = params.getTheta(); RealMatrix C = CFixedRand.getCombinedMatrix(); // RealMatrix beta = params.getScaledBeta(); // RealMatrix U = params.getWithinSubjectContrast(); // thetaHat = C * beta * U RealMatrix thetaHat = C.multiply(beta.multiply(U)); debug("C", C); debug("beta", beta); debug("U", U); debug("thetaHat = C * beta * U", thetaHat); // thetaDiff = thetaHat - thetaNull RealMatrix thetaDiff = thetaHat.subtract(thetaNull); debug("thetaNull", thetaNull); debug("thetaDiff = thetaHat - thetaNull", thetaDiff); // TODO: specific to HLT or UNIREP RealMatrix sigmaStarInverse = getSigmaStarInverse(U, sigmaError, test); debug("sigmaStarInverse", sigmaStarInverse); RealMatrix H1matrix = thetaDiff.transpose().multiply(T1).multiply(thetaDiff).multiply(sigmaStarInverse); debug("H1matrix = thetaDiff transpose * T1 * thetaDiff * sigmaStarInverse", H1matrix); H1 = H1matrix.getTrace(); debug("H1 = " + H1); // Matrix which represents the non-centrality parameter as a linear combination of chi-squared r.v.'s. S = FT1.transpose().multiply(thetaDiff).multiply(sigmaStarInverse).multiply(thetaDiff.transpose()) .multiply(FT1).scalarMultiply(1 / H1); debug("S = FT1 transpose * thetaDiff * sigmaStar inverse * thetaDiff transpose * FT1 * (1/H1)", S); // We use the S matrix to generate the F-critical, numerical df's, and denominator df's // for a central F distribution. The resulting F distribution is used as an approximation // for the distribution of the non-centrality parameter. // See formulas 18-21 and A8,A10 from Glueck & Muller (2003) for details. EigenDecomposition sEigenDecomp = new EigenDecomposition(S); sEigenValues = sEigenDecomp.getRealEigenvalues(); // calculate H0 if (sEigenValues.length > 0) H0 = H1 * (1 - sEigenValues[0]); if (H0 <= 0) H0 = 0; // count the # of positive eigen values for (double value : sEigenValues) { if (value > 0) sStar++; } // TODO: throw error if sStar is <= 0 // TODO: NO: throw error if sStar != sEigenValues.length instead??? double stddevG = Math.sqrt(sigmaG.getEntry(0, 0)); RealMatrix svec = sEigenDecomp.getVT(); mzSq = svec.multiply(FT1.transpose()).multiply(CGaussian).scalarMultiply(1 / stddevG); for (int i = 0; i < mzSq.getRowDimension(); i++) { for (int j = 0; j < mzSq.getColumnDimension(); j++) { double entry = mzSq.getEntry(i, j); mzSq.setEntry(i, j, entry * entry); // TODO: is there an apache function to do this? } } debug("exiting initialize normally"); } catch (RuntimeException e) { LOGGER.warn("exiting initialize abnormally", e); throw new PowerException(e.getMessage(), PowerErrorEnum.INVALID_DISTRIBUTION_NONCENTRALITY_PARAMETER); } }