List of usage examples for org.apache.commons.math3.linear RealMatrix getColumnMatrix
RealMatrix getColumnMatrix(int column) throws OutOfRangeException;
From source file:edu.cudenver.bios.matrix.OrthogonalPolynomials.java
/** * Create a within or between subject contrast (C) for polynomial trends * for an arbitrary list of factors. The returned collection includes * the grand mean contrast, all 1-factor main effect contrasts, and * all possible interaction contrasts//from w w w . j a v a 2s.c om * * @param factorList list of factors, including name and value information * @return polynomial contrast collection. * @throws IllegalArgumentException */ public static OrthogonalPolynomialContrastCollection buildContrastCollection(List<Factor> factorList, boolean between) throws IllegalArgumentException { if (factorList == null || factorList.size() <= 0) throw new IllegalArgumentException("no factors specified"); ArrayList<RealMatrix> zeroTrendList = new ArrayList<RealMatrix>(factorList.size()); ArrayList<RealMatrix> factorTrendList = new ArrayList<RealMatrix>(factorList.size()); for (Factor factor : factorList) { double[] centered = centerAndScale(factor.getValues()); RealMatrix poly = OrthogonalPolynomials.orthogonalPolynomialCoefficients(centered, centered.length - 1); zeroTrendList.add(poly.getColumnMatrix(0)); factorTrendList.add(poly.getSubMatrix(0, poly.getRowDimension() - 1, 1, poly.getColumnDimension() - 1)); } OrthogonalPolynomialContrastCollection results = new OrthogonalPolynomialContrastCollection(); /* * We need to create contrasts for every possible combination of * zero and factor trends. The possible combinations can be * enumerated by the binary representation of the numbers * 0 - 2^(#factors). In this case, 0 indicates that the zero-trend contrast should * be included in the Kronecker product for a given factor, and 1 indicates that the * factor trend contrast should be included in the Kronecker product. */ ArrayList<RealMatrix> kroneckerList = new ArrayList<RealMatrix>(factorList.size()); ArrayList<Factor> activeFactorList = new ArrayList<Factor>(factorList.size()); // build the grand mean for (RealMatrix zeroTrend : zeroTrendList) kroneckerList.add(zeroTrend); if (between) results.addContrast( new OrthogonalPolynomialContrast(MatrixUtils.getKroneckerProduct(kroneckerList).transpose())); else results.addContrast(new OrthogonalPolynomialContrast(MatrixUtils.getKroneckerProduct(kroneckerList))); // loop over the remaining contrasts int totalContrasts = (int) Math.pow(2.0, (double) factorList.size()); for (int i = 1; i < totalContrasts; i++) { kroneckerList.clear(); activeFactorList.clear(); int mask = 1; for (int factorIdx = 0; factorIdx < factorList.size(); factorIdx++, mask *= 2) { if ((i & mask) != 0) { kroneckerList.add(factorTrendList.get(factorIdx)); activeFactorList.add(factorList.get(factorIdx)); } else { kroneckerList.add(zeroTrendList.get(factorIdx)); } } // add the appropriate contrast type // note that if "i" is a power of 2 then we have a main effect contrast, else interaction RealMatrix contrast = null; if (between) contrast = MatrixUtils.getKroneckerProduct(kroneckerList).transpose(); else contrast = MatrixUtils.getKroneckerProduct(kroneckerList); results.addContrast(new OrthogonalPolynomialContrast( ((i & (i - 1)) == 0 ? ContrastType.MAIN_EFFECT : ContrastType.INTERACTION), activeFactorList, contrast)); } return results; }
From source file:com.yahoo.egads.utilities.SpectralMethods.java
public static RealMatrix mFilter(RealMatrix data, int windowSize, FilteringMethod method, double methodParameter) { int n = data.getRowDimension(); int m = data.getColumnDimension(); int k = n - windowSize + 1; int i = 0, ind = 0; double[] temp; double sum = 0; RealMatrix hankelMat = SpectralMethods.createHankelMatrix(data, windowSize); SingularValueDecomposition svd = new SingularValueDecomposition(hankelMat); double[] singularValues = svd.getSingularValues(); switch (method) { case VARIANCE: temp = new double[singularValues.length - 1]; for (i = 1; i < singularValues.length; ++i) { sum += (singularValues[i] * singularValues[i]); }// w w w. j av a 2s .co m for (i = 0; i < temp.length; ++i) { temp[i] = (singularValues[i + 1] * singularValues[i + 1]) / sum; } sum = 0; for (i = temp.length - 1; i >= 0; --i) { sum += temp[i]; if (sum >= 1 - methodParameter) { ind = i; break; } } break; case EXPLICIT: ind = (int) Math.max(Math.min(methodParameter - 1, singularValues.length - 1), 0); break; case K_GAP: final double[] eigenGaps = new double[singularValues.length - 1]; Integer[] index = new Integer[singularValues.length - 1]; for (i = 0; i < eigenGaps.length; ++i) { eigenGaps[i] = singularValues[i] - singularValues[i + 1]; index[i] = i; } Arrays.sort(index, new Comparator<Integer>() { @Override public int compare(Integer o1, Integer o2) { return Double.compare(eigenGaps[o1], eigenGaps[o2]); } }); int maxIndex = 0; for (i = index.length - (int) methodParameter; i < index.length; ++i) { if (index[i] > maxIndex) { maxIndex = index[i]; } } ind = Math.min(maxIndex, singularValues.length / 3); break; case SMOOTHNESS: double[] variances = new double[singularValues.length]; for (i = 1; i < singularValues.length; ++i) { variances[i] = (singularValues[i] * singularValues[i]); } variances[0] = variances[1]; double smoothness = SpectralMethods .computeSmoothness(Arrays.copyOfRange(variances, 1, variances.length)); if (methodParameter - smoothness < 0.01) { methodParameter += 0.01; } double invalidS = smoothness; int validIndex = 1, invalidIndex = singularValues.length; while (true) { if (invalidS >= methodParameter) { ind = invalidIndex - 1; break; } else if (invalidIndex - validIndex <= 1) { ind = validIndex - 1; break; } int ii = (validIndex + invalidIndex) / 2; double[] tempVariances = Arrays.copyOf(Arrays.copyOfRange(variances, 0, ii + 1), singularValues.length); double s = SpectralMethods.computeSmoothness(tempVariances); if (s >= methodParameter) { validIndex = ii; } else { invalidIndex = ii; invalidS = s; } } break; case EIGEN_RATIO: int startIndex = 0, endIndex = singularValues.length - 1; if (singularValues[endIndex] / singularValues[0] >= methodParameter) { ind = endIndex; } else { while (true) { int midIndex = (startIndex + endIndex) / 2; if (singularValues[midIndex] / singularValues[0] >= methodParameter) { if (singularValues[midIndex + 1] / singularValues[0] < methodParameter) { ind = midIndex; break; } else { startIndex = midIndex; } } else { endIndex = midIndex; } } } break; case GAP_RATIO: double[] gaps = new double[singularValues.length - 1]; for (i = 0; i < gaps.length; ++i) { gaps[i] = singularValues[i] - singularValues[i + 1]; } ind = 0; for (i = gaps.length - 1; i >= 0; --i) { if (gaps[i] / singularValues[0] >= methodParameter) { ind = i; break; } } break; default: ind = singularValues.length - 1; break; } ind = Math.max(0, Math.min(ind, singularValues.length - 1)); RealMatrix truncatedHankelMatrix = MatrixUtils.createRealMatrix(k, m * windowSize); RealMatrix mU = svd.getU(); RealMatrix mVT = svd.getVT(); for (i = 0; i <= ind; ++i) { truncatedHankelMatrix = truncatedHankelMatrix .add(mU.getColumnMatrix(i).multiply(mVT.getRowMatrix(i)).scalarMultiply(singularValues[i])); } return SpectralMethods.averageHankelMatrix(truncatedHankelMatrix, windowSize); }
From source file:edu.ucdenver.bios.powersvc.resource.ContrastHelper.java
/** * Create a trend contrast of the specified type * @param spacing list of integer positions representing spacing of measurements * @param trendType type of trend contrast * @param transpose if true, return the transpose of the contrast * @return trend contrast/* w w w . j a v a 2 s. c o m*/ */ private static RealMatrix getTrendContrast(double[] spacing, HypothesisTrendTypeEnum trendType, boolean transpose) { int levels = spacing.length; // get all possible polynomial trends RealMatrix allTrendContrast = OrthogonalPolynomials.orthogonalPolynomialCoefficients(spacing, levels - 1); // select a subset of the polynomial trends based on the hypothesis of interest RealMatrix trendContrast = null; switch (trendType) { case NONE: RealMatrix negIdentity = org.apache.commons.math3.linear.MatrixUtils .createRealIdentityMatrix(levels - 1).scalarMultiply(-1); RealMatrix column1s = MatrixUtils.getRealMatrixWithFilledValue(levels - 1, 1, 1); trendContrast = MatrixUtils.getHorizontalAppend(column1s, negIdentity).transpose(); break; case CHANGE_FROM_BASELINE: trendContrast = MatrixUtils.getRealMatrixWithFilledValue(levels, 1, 0); trendContrast.setEntry(0, 0, 1); trendContrast.setEntry(levels - 1, 0, -1); break; case ALL_POLYNOMIAL: trendContrast = allTrendContrast; break; case LINEAR: if (allTrendContrast.getRowDimension() > 1) { trendContrast = allTrendContrast.getColumnMatrix(1); } break; case QUADRATIC: if (allTrendContrast.getRowDimension() > 2) { trendContrast = allTrendContrast.getColumnMatrix(2); } break; case CUBIC: if (allTrendContrast.getRowDimension() > 3) { trendContrast = allTrendContrast.getColumnMatrix(3); } break; } if (transpose) { return trendContrast.transpose(); } else { return trendContrast; } }
From source file:edu.cudenver.bios.matrix.GramSchmidtOrthonormalization.java
/** * Perform Gram Schmidt Orthonormalization on the specified * matrix. The matrix A (mxn) is decomposed into two matrices * Q (mxn), R (nxn) such that/*w w w . j ava 2 s . c o m*/ * <ul> * <li>A = QR * <li>Q'Q = Identity * <li>R is upper triangular * </ul> * The resulting Q, R matrices can be retrieved with the getQ() * and getR() functions. * * @param matrix */ public GramSchmidtOrthonormalization(RealMatrix matrix) { if (matrix == null) throw new IllegalArgumentException("Null matrix"); // create the Q, R matrices int m = matrix.getRowDimension(); int n = matrix.getColumnDimension(); Q = MatrixUtils.createRealMatrix(m, n); R = MatrixUtils.createRealMatrix(n, n); // perform Gram Schmidt process using the following algorithm // let w<n> be the resulting orthonormal column vectors // let v<n> be the columns of the incoming matrix // w1 = (1/norm(v1))*v1 // ... // wj = 1/norm(vj - projectionVj-1Vj)*[vj - projectionVj-1Vj] // where projectionVj-1Vj = (w1 * vj) * w1 + (w2 * vj) * w2 + ... + (wj-1 * vj) * wj-1 // for (int i = 0; i < n; i++) { RealMatrix v = matrix.getColumnMatrix(i); for (int j = 0; j < i; j++) { RealMatrix Qj = Q.getColumnMatrix(j); double value = Qj.transpose().multiply(v).getEntry(0, 0); R.setEntry(j, i, value); v = v.subtract(Qj.scalarMultiply(value)); } double norm = v.getFrobeniusNorm(); R.setEntry(i, i, norm); Q.setColumnMatrix(i, v.scalarMultiply(1 / norm)); } }
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(/* w w w . j a v a 2 s . co 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:lirmm.inria.fr.math.BigSparseRealMatrixTest.java
@Test public void testSetColumnMatrix() { RealMatrix m = new BigSparseRealMatrix(subTestData); RealMatrix mColumn3 = new BigSparseRealMatrix(subColumn3); Assert.assertNotSame(mColumn3, m.getColumnMatrix(1)); m.setColumnMatrix(1, mColumn3);// w w w. ja va2s .c o m Assert.assertEquals(mColumn3, m.getColumnMatrix(1)); try { m.setColumnMatrix(-1, mColumn3); Assert.fail("Expecting OutOfRangeException"); } catch (OutOfRangeException ex) { // expected } try { m.setColumnMatrix(0, m); Assert.fail("Expecting MatrixDimensionMismatchException"); } catch (MatrixDimensionMismatchException ex) { // expected } }
From source file:lirmm.inria.fr.math.BigSparseRealMatrixTest.java
@Test public void testGetColumnMatrix() { RealMatrix m = new BigSparseRealMatrix(subTestData); RealMatrix mColumn1 = new BigSparseRealMatrix(subColumn1); RealMatrix mColumn3 = new BigSparseRealMatrix(subColumn3); Assert.assertEquals("Column1", mColumn1, m.getColumnMatrix(1)); Assert.assertEquals("Column3", mColumn3, m.getColumnMatrix(3)); try {/*from w w w . java 2s .c om*/ m.getColumnMatrix(-1); Assert.fail("Expecting OutOfRangeException"); } catch (OutOfRangeException ex) { // expected } try { m.getColumnMatrix(4); Assert.fail("Expecting OutOfRangeException"); } catch (OutOfRangeException ex) { // expected } }
From source file:org.akvo.caddisfly.sensor.colorimetry.strip.calibration.CalibrationCard.java
@NonNull private static Mat doIlluminationCorrection(@NonNull Mat imgLab, @NonNull CalibrationData calData) { // create HLS image for homogeneous illumination calibration int pHeight = imgLab.rows(); int pWidth = imgLab.cols(); RealMatrix points = createWhitePointMatrix(imgLab, calData); // create coefficient matrix for all three variables L,A,B // the model for all three is y = ax + bx^2 + cy + dy^2 + exy + f // 6th row is the constant 1 RealMatrix coefficient = new Array2DRowRealMatrix(points.getRowDimension(), 6); coefficient.setColumnMatrix(0, points.getColumnMatrix(0)); coefficient.setColumnMatrix(2, points.getColumnMatrix(1)); //create constant, x^2, y^2 and xy terms for (int i = 0; i < points.getRowDimension(); i++) { coefficient.setEntry(i, 1, Math.pow(coefficient.getEntry(i, 0), 2)); // x^2 coefficient.setEntry(i, 3, Math.pow(coefficient.getEntry(i, 2), 2)); // y^2 coefficient.setEntry(i, 4, coefficient.getEntry(i, 0) * coefficient.getEntry(i, 2)); // xy coefficient.setEntry(i, 5, 1d); // constant = 1 }// w w w.ja v a 2 s . c om // create vectors RealVector L = points.getColumnVector(2); RealVector A = points.getColumnVector(3); RealVector B = points.getColumnVector(4); // solve the least squares problem for all three variables DecompositionSolver solver = new SingularValueDecomposition(coefficient).getSolver(); RealVector solutionL = solver.solve(L); RealVector solutionA = solver.solve(A); RealVector solutionB = solver.solve(B); // get individual coefficients float La = (float) solutionL.getEntry(0); float Lb = (float) solutionL.getEntry(1); float Lc = (float) solutionL.getEntry(2); float Ld = (float) solutionL.getEntry(3); float Le = (float) solutionL.getEntry(4); float Lf = (float) solutionL.getEntry(5); float Aa = (float) solutionA.getEntry(0); float Ab = (float) solutionA.getEntry(1); float Ac = (float) solutionA.getEntry(2); float Ad = (float) solutionA.getEntry(3); float Ae = (float) solutionA.getEntry(4); float Af = (float) solutionA.getEntry(5); float Ba = (float) solutionB.getEntry(0); float Bb = (float) solutionB.getEntry(1); float Bc = (float) solutionB.getEntry(2); float Bd = (float) solutionB.getEntry(3); float Be = (float) solutionB.getEntry(4); float Bf = (float) solutionB.getEntry(5); // compute mean (the luminosity value of the plane in the middle of the image) float L_mean = (float) (0.5 * La * pWidth + 0.5 * Lc * pHeight + Lb * pWidth * pWidth / 3.0 + Ld * pHeight * pHeight / 3.0 + Le * 0.25 * pHeight * pWidth + Lf); float A_mean = (float) (0.5 * Aa * pWidth + 0.5 * Ac * pHeight + Ab * pWidth * pWidth / 3.0 + Ad * pHeight * pHeight / 3.0 + Ae * 0.25 * pHeight * pWidth + Af); float B_mean = (float) (0.5 * Ba * pWidth + 0.5 * Bc * pHeight + Bb * pWidth * pWidth / 3.0 + Bd * pHeight * pHeight / 3.0 + Be * 0.25 * pHeight * pWidth + Bf); // Correct image // we do this per row. We tried to do it in one block, but there is no speed difference. byte[] temp = new byte[imgLab.cols() * imgLab.channels()]; int valL, valA, valB; int ii, ii3; float iiSq, iSq; int imgCols = imgLab.cols(); int imgRows = imgLab.rows(); // use lookup tables to speed up computation // create lookup tables float[] L_aii = new float[imgCols]; float[] L_biiSq = new float[imgCols]; float[] A_aii = new float[imgCols]; float[] A_biiSq = new float[imgCols]; float[] B_aii = new float[imgCols]; float[] B_biiSq = new float[imgCols]; float[] Lci = new float[imgRows]; float[] LdiSq = new float[imgRows]; float[] Aci = new float[imgRows]; float[] AdiSq = new float[imgRows]; float[] Bci = new float[imgRows]; float[] BdiSq = new float[imgRows]; for (ii = 0; ii < imgCols; ii++) { iiSq = ii * ii; L_aii[ii] = La * ii; L_biiSq[ii] = Lb * iiSq; A_aii[ii] = Aa * ii; A_biiSq[ii] = Ab * iiSq; B_aii[ii] = Ba * ii; B_biiSq[ii] = Bb * iiSq; } for (int i = 0; i < imgRows; i++) { iSq = i * i; Lci[i] = Lc * i; LdiSq[i] = Ld * iSq; Aci[i] = Ac * i; AdiSq[i] = Ad * iSq; Bci[i] = Bc * i; BdiSq[i] = Bd * iSq; } // We can also improve the performance of the i,ii term, if we want, but it won't make much difference. for (int i = 0; i < imgRows; i++) { // y imgLab.get(i, 0, temp); ii3 = 0; for (ii = 0; ii < imgCols; ii++) { //x valL = capValue( Math.round((temp[ii3] & 0xFF) - (L_aii[ii] + L_biiSq[ii] + Lci[i] + LdiSq[i] + Le * i * ii + Lf) + L_mean), 0, 255); valA = capValue( Math.round((temp[ii3 + 1] & 0xFF) - (A_aii[ii] + A_biiSq[ii] + Aci[i] + AdiSq[i] + Ae * i * ii + Af) + A_mean), 0, 255); valB = capValue( Math.round((temp[ii3 + 2] & 0xFF) - (B_aii[ii] + B_biiSq[ii] + Bci[i] + BdiSq[i] + Be * i * ii + Bf) + B_mean), 0, 255); temp[ii3] = (byte) valL; temp[ii3 + 1] = (byte) valA; temp[ii3 + 2] = (byte) valB; ii3 += 3; } imgLab.put(i, 0, temp); } return imgLab; }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * {@inheritDoc}/* ww w . j av a 2s .c om*/ */ @Override protected PointValuePair doOptimize() { // -------------------- Initialization -------------------------------- isMinimize = getGoalType().equals(GoalType.MINIMIZE); final double[] guess = getStartPoint(); // number of objective variables/problem dimension dimension = guess.length; initializeCMA(guess); iterations = 0; double bestValue = (isMinimize ? Double.MAX_VALUE : Double.MIN_VALUE); push(fitnessHistory, bestValue); PointValuePair optimum = new PointValuePair(getStartPoint(), isMinimize ? bestValue : -bestValue); PointValuePair lastResult = null; // -------------------- Generation Loop -------------------------------- EvaluatedPopulation<double[]> evaluatedPopulation = null; Stopwatch stopwatch = Stopwatch.createUnstarted(); generationLoop: for (iterations = 1; iterations <= maxIterations; iterations++) { stopwatch.reset(); stopwatch.start(); incrementIterationCount(); // Generate and evaluate lambda offspring final RealMatrix arz = randn1(dimension, lambda); final RealMatrix arx = zeros(dimension, lambda); final double[] fitness = new double[lambda]; // generate random offspring for (int k = 0; k < lambda; k++) { RealMatrix arxk = null; for (int i = 0; i < checkFeasableCount + 1; i++) { if (diagonalOnly <= 0) { arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k)).scalarMultiply(sigma)); // m + sig * Normal(0,C) } else { arxk = xmean.add(times(diagD, arz.getColumnMatrix(k)).scalarMultiply(sigma)); } //if (i >= checkFeasableCount || // fitfun.isFeasible(arxk.getColumn(0))) { // break; //} // regenerate random arguments for row arz.setColumn(k, randn(dimension)); } copyColumn(arxk, 0, arx, k); //try { // valuePenaltyPairs[k] = fitfun.value(arx.getColumn(k)); // compute fitness //} catch (TooManyEvaluationsException e) { // break generationLoop; //} } double newPopTime = stopwatch.elapsed(TimeUnit.MILLISECONDS) / 1000.0; stopwatch.reset(); stopwatch.start(); ArrayList<double[]> population = new ArrayList<>(lambda); // This is mine. I ignore constraints. for (int k = 0; k < lambda; ++k) { population.add(arx.getColumn(k)); } evaluatedPopulation = populationEvaluator.evaluate(population, iterations - 1, random); final ValuePenaltyPair[] valuePenaltyPairs = new ValuePenaltyPair[lambda]; for (int k = 0; k < lambda; ++k) { valuePenaltyPairs[k] = new ValuePenaltyPair(evaluatedPopulation.getPopulation().get(k).getFitness(), 0.0); } // Compute fitnesses by adding value and penalty after scaling by value range. double valueRange = valueRange(valuePenaltyPairs); for (int iValue = 0; iValue < valuePenaltyPairs.length; iValue++) { fitness[iValue] = valuePenaltyPairs[iValue].value + valuePenaltyPairs[iValue].penalty * valueRange; if (!isMinimize) fitness[iValue] = -fitness[iValue]; } double evalTime = stopwatch.elapsed(TimeUnit.MILLISECONDS) / 1000.0; stopwatch.reset(); stopwatch.start(); // Sort by fitness and compute weighted mean into xmean final int[] arindex = sortedIndices(fitness); // Calculate new xmean, this is selection and recombination final RealMatrix xold = xmean; // for speed up of Eq. (2) and (3) final RealMatrix bestArx = selectColumns(arx, MathArrays.copyOf(arindex, mu)); xmean = bestArx.multiply(weights); final RealMatrix bestArz = selectColumns(arz, MathArrays.copyOf(arindex, mu)); final RealMatrix zmean = bestArz.multiply(weights); final boolean hsig = updateEvolutionPaths(zmean, xold); if (diagonalOnly <= 0) { updateCovariance(hsig, bestArx, arz, arindex, xold); } else { updateCovarianceDiagonalOnly(hsig, bestArz); } // Adapt step size sigma - Eq. (5) sigma *= FastMath.exp(FastMath.min(1, (normps / chiN - 1) * cs / damps)); final double bestFitness = fitness[arindex[0]]; final double worstFitness = fitness[arindex[arindex.length - 1]]; if (bestValue > bestFitness) { bestValue = bestFitness; lastResult = optimum; optimum = new PointValuePair(bestArx.getColumn(0), isMinimize ? bestFitness : -bestFitness); if (getConvergenceChecker() != null && lastResult != null && getConvergenceChecker().converged(iterations, optimum, lastResult)) { break generationLoop; } } // handle termination criteria // Break, if fitness is good enough if (stopFitness != 0 && bestFitness < (isMinimize ? stopFitness : -stopFitness)) { break generationLoop; } final double[] sqrtDiagC = sqrt(diagC).getColumn(0); final double[] pcCol = pc.getColumn(0); for (int i = 0; i < dimension; i++) { if (sigma * FastMath.max(FastMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) { break; } if (i >= dimension - 1) { break generationLoop; } } for (int i = 0; i < dimension; i++) { if (sigma * sqrtDiagC[i] > stopTolUpX) { break generationLoop; } } final double historyBest = min(fitnessHistory); final double historyWorst = max(fitnessHistory); if (iterations > 2 && FastMath.max(historyWorst, worstFitness) - FastMath.min(historyBest, bestFitness) < stopTolFun) { break generationLoop; } if (iterations > fitnessHistory.length && historyWorst - historyBest < stopTolHistFun) { break generationLoop; } // condition number of the covariance matrix exceeds 1e14 if (max(diagD) / min(diagD) > 1e7) { break generationLoop; } // user defined termination if (getConvergenceChecker() != null) { final PointValuePair current = new PointValuePair(bestArx.getColumn(0), isMinimize ? bestFitness : -bestFitness); if (lastResult != null && getConvergenceChecker().converged(iterations, current, lastResult)) { break generationLoop; } lastResult = current; } // Adjust step size in case of equal function values (flat fitness) if (bestValue == fitness[arindex[(int) (0.1 + lambda / 4.)]]) { sigma *= FastMath.exp(0.2 + cs / damps); } if (iterations > 2 && FastMath.max(historyWorst, bestFitness) - FastMath.min(historyBest, bestFitness) == 0) { sigma *= FastMath.exp(0.2 + cs / damps); } // store best in history push(fitnessHistory, bestFitness); if (generateStatistics) { statisticsSigmaHistory.add(sigma); statisticsFitnessHistory.add(bestFitness); statisticsMeanHistory.add(xmean.transpose()); statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5)); } double cmaesTime = stopwatch.elapsed(TimeUnit.MILLISECONDS) / 1000.0; stopwatch.reset(); stopwatch.start(); listener.onNextIteraction(evaluatedPopulation); double listernerTime = stopwatch.elapsed(TimeUnit.MILLISECONDS) / 1000.0; logger.info(String.format("NewPop: %.2f, Eval: %.2f, CMAES: %.2f, Listerner: %.2f", newPopTime, evalTime, cmaesTime, listernerTime)); } listener.onLastIteraction(evaluatedPopulation); return optimum; }