List of usage examples for org.apache.commons.math3.util MathArrays copyOf
public static double[] copyOf(double[] source, int len)
From source file:org.pmad.gmm.MyMixMNDEM.java
/** * Creates an object to fit a multivariate normal mixture model to data. * * @param data Data to use in fitting procedure * @throws NotStrictlyPositiveException if data has no rows * @throws DimensionMismatchException if rows of data have different numbers * of columns/*w w w . ja v a 2s.c o m*/ * @throws NumberIsTooSmallException if the number of columns in the data is * less than 2 */ public MyMixMNDEM(double[][] data) throws NotStrictlyPositiveException, DimensionMismatchException, NumberIsTooSmallException { if (data.length < 1) { throw new NotStrictlyPositiveException(data.length); } this.data = new double[data.length][data[0].length]; for (int i = 0; i < data.length; i++) { if (data[i].length != data[0].length) { // Jagged arrays not allowed throw new DimensionMismatchException(data[i].length, data[0].length); } if (data[i].length < 2) { throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_TOO_SMALL, data[i].length, 2, true); } this.data[i] = MathArrays.copyOf(data[i], data[i].length); } }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * {@inheritDoc}/*from ww w .j a va 2 s . c o m*/ */ @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; }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * Update of the covariance matrix C.// w w w.ja v a 2s . c om * * @param hsig Flag indicating a small correction. * @param bestArx Fitness-sorted matrix of the argument vectors producing the current offspring. * @param arz Unsorted matrix containing the gaussian random values of the current offspring. * @param arindex Indices indicating the fitness-order of the current offspring. * @param xold xmean matrix of the previous generation. */ private void updateCovariance(boolean hsig, final RealMatrix bestArx, final RealMatrix arz, final int[] arindex, final RealMatrix xold) { double negccov = 0; if (ccov1 + ccovmu > 0) { final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu)).scalarMultiply(1 / sigma); // mu difference vectors final RealMatrix roneu = pc.multiply(pc.transpose()).scalarMultiply(ccov1); // rank one update // minor correction if hsig==false double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc); oldFac += 1 - ccov1 - ccovmu; if (isActiveCMA) { // Adapt covariance matrix C active CMA negccov = (1 - ccovmu) * 0.25 * mueff / (FastMath.pow(dimension + 2, 1.5) + 2 * mueff); // keep at least 0.66 in all directions, small popsize are most // critical final double negminresidualvariance = 0.66; // where to make up for the variance loss final double negalphaold = 0.5; // prepare vectors, compute negative updating matrix Cneg final int[] arReverseIndex = reverse(arindex); RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu)); RealMatrix arnorms = sqrt(sumRows(square(arzneg))); final int[] idxnorms = sortedIndices(arnorms.getRow(0)); final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms); final int[] idxReverse = reverse(idxnorms); final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse); arnorms = divide(arnormsReverse, arnormsSorted); final int[] idxInv = inverse(idxnorms); final RealMatrix arnormsInv = selectColumns(arnorms, idxInv); // check and set learning rate negccov final double negcovMax = (1 - negminresidualvariance) / square(arnormsInv).multiply(weights).getEntry(0, 0); if (negccov > negcovMax) { negccov = negcovMax; } arzneg = times(arzneg, repmat(arnormsInv, dimension, 1)); final RealMatrix artmp = BD.multiply(arzneg); final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose()); oldFac += negalphaold * negccov; C = C.scalarMultiply(oldFac).add(roneu) // regard old matrix .add(arpos.scalarMultiply( // plus rank one update ccovmu + (1 - negalphaold) * negccov) // plus rank mu update .multiply(times(repmat(weights, 1, dimension), arpos.transpose()))) .subtract(Cneg.scalarMultiply(negccov)); } else { // Adapt covariance matrix C - nonactive C = C.scalarMultiply(oldFac) // regard old matrix .add(roneu) // plus rank one update .add(arpos.scalarMultiply(ccovmu) // plus rank mu update .multiply(times(repmat(weights, 1, dimension), arpos.transpose()))); } } updateBDFast(negccov); }