1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
19
20 import java.util.ArrayList;
21 import java.util.Arrays;
22 import java.util.List;
23
24 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
25 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
26 import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
27 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
28 import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
29 import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
30 import org.apache.commons.math4.legacy.linear.EigenDecomposition;
31 import org.apache.commons.math4.legacy.linear.MatrixUtils;
32 import org.apache.commons.math4.legacy.linear.RealMatrix;
33 import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
34 import org.apache.commons.math4.legacy.optim.OptimizationData;
35 import org.apache.commons.math4.legacy.optim.PointValuePair;
36 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
37 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.PopulationSize;
38 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.Sigma;
39 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimizer;
40 import org.apache.commons.rng.UniformRandomProvider;
41 import org.apache.commons.statistics.distribution.ContinuousDistribution;
42 import org.apache.commons.statistics.distribution.NormalDistribution;
43 import org.apache.commons.math4.core.jdkmath.JdkMath;
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84 public class CMAESOptimizer
85 extends MultivariateOptimizer {
86
87
88
89
90
91
92
93
94 private int lambda;
95
96
97
98
99
100
101
102 private final boolean isActiveCMA;
103
104
105
106
107 private final int checkFeasableCount;
108
109
110
111 private double[] inputSigma;
112
113 private int dimension;
114
115
116
117
118
119
120
121
122 private int diagonalOnly;
123
124 private boolean isMinimize = true;
125
126 private final boolean generateStatistics;
127
128
129
130 private final int maxIterations;
131
132 private final double stopFitness;
133
134 private double stopTolUpX;
135
136 private double stopTolX;
137
138 private double stopTolFun;
139
140 private double stopTolHistFun;
141
142
143
144 private int mu;
145
146 private double logMu2;
147
148 private RealMatrix weights;
149
150 private double mueff;
151
152
153
154 private double sigma;
155
156 private double cc;
157
158 private double cs;
159
160 private double damps;
161
162 private double ccov1;
163
164 private double ccovmu;
165
166 private double chiN;
167
168 private double ccov1Sep;
169
170 private double ccovmuSep;
171
172
173
174 private RealMatrix xmean;
175
176 private RealMatrix pc;
177
178 private RealMatrix ps;
179
180 private double normps;
181
182 private RealMatrix B;
183
184 private RealMatrix D;
185
186 private RealMatrix BD;
187
188 private RealMatrix diagD;
189
190 private RealMatrix C;
191
192 private RealMatrix diagC;
193
194 private int iterations;
195
196
197 private double[] fitnessHistory;
198
199 private int historySize;
200
201
202 private final ContinuousDistribution.Sampler random;
203
204
205 private final List<Double> statisticsSigmaHistory = new ArrayList<>();
206
207 private final List<RealMatrix> statisticsMeanHistory = new ArrayList<>();
208
209 private final List<Double> statisticsFitnessHistory = new ArrayList<>();
210
211 private final List<RealMatrix> statisticsDHistory = new ArrayList<>();
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228 public CMAESOptimizer(int maxIterations,
229 double stopFitness,
230 boolean isActiveCMA,
231 int diagonalOnly,
232 int checkFeasableCount,
233 UniformRandomProvider rng,
234 boolean generateStatistics,
235 ConvergenceChecker<PointValuePair> checker) {
236 super(checker);
237 this.maxIterations = maxIterations;
238 this.stopFitness = stopFitness;
239 this.isActiveCMA = isActiveCMA;
240 this.diagonalOnly = diagonalOnly;
241 this.checkFeasableCount = Math.max(0, checkFeasableCount);
242 this.random = NormalDistribution.of(0, 1).createSampler(rng);
243 this.generateStatistics = generateStatistics;
244 }
245
246
247
248
249 public List<Double> getStatisticsSigmaHistory() {
250 return statisticsSigmaHistory;
251 }
252
253
254
255
256 public List<RealMatrix> getStatisticsMeanHistory() {
257 return statisticsMeanHistory;
258 }
259
260
261
262
263 public List<Double> getStatisticsFitnessHistory() {
264 return statisticsFitnessHistory;
265 }
266
267
268
269
270 public List<RealMatrix> getStatisticsDHistory() {
271 return statisticsDHistory;
272 }
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308 @Override
309 public PointValuePair optimize(OptimizationData... optData) {
310
311 return super.optimize(optData);
312 }
313
314
315 @Override
316 protected PointValuePair doOptimize() {
317
318 isMinimize = getGoalType().equals(GoalType.MINIMIZE);
319 final FitnessFunction fitfun = new FitnessFunction();
320 final double[] guess = getStartPoint();
321
322 dimension = guess.length;
323 initializeCMA(guess);
324 iterations = 0;
325 ValuePenaltyPair valuePenalty = fitfun.value(guess);
326 double bestValue = valuePenalty.value+valuePenalty.penalty;
327 push(fitnessHistory, bestValue);
328 PointValuePair optimum
329 = new PointValuePair(getStartPoint(),
330 isMinimize ? bestValue : -bestValue);
331 PointValuePair lastResult = null;
332
333
334 generationLoop:
335 for (iterations = 1; iterations <= maxIterations; iterations++) {
336 incrementIterationCount();
337
338
339 final RealMatrix arz = randn1(dimension, lambda);
340 final RealMatrix arx = zeros(dimension, lambda);
341 final double[] fitness = new double[lambda];
342 final ValuePenaltyPair[] valuePenaltyPairs = new ValuePenaltyPair[lambda];
343
344 for (int k = 0; k < lambda; k++) {
345 RealMatrix arxk = null;
346 for (int i = 0; i <= checkFeasableCount; i++) {
347 if (diagonalOnly <= 0) {
348 arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k))
349 .scalarMultiply(sigma));
350 } else {
351 arxk = xmean.add(times(diagD,arz.getColumnMatrix(k))
352 .scalarMultiply(sigma));
353 }
354 if (i >= checkFeasableCount ||
355 fitfun.isFeasible(arxk.getColumn(0))) {
356 break;
357 }
358
359 arz.setColumn(k, randn(dimension));
360 }
361 copyColumn(arxk, 0, arx, k);
362 try {
363 valuePenaltyPairs[k] = fitfun.value(arx.getColumn(k));
364 } catch (TooManyEvaluationsException e) {
365 break generationLoop;
366 }
367 }
368
369 double valueRange = valueRange(valuePenaltyPairs);
370 for (int iValue=0;iValue<valuePenaltyPairs.length;iValue++) {
371 fitness[iValue] = valuePenaltyPairs[iValue].value + valuePenaltyPairs[iValue].penalty*valueRange;
372 }
373
374 final int[] arindex = sortedIndices(fitness);
375
376 final RealMatrix xold = xmean;
377 final RealMatrix bestArx = selectColumns(arx, Arrays.copyOf(arindex, mu));
378 xmean = bestArx.multiply(weights);
379 final RealMatrix bestArz = selectColumns(arz, Arrays.copyOf(arindex, mu));
380 final RealMatrix zmean = bestArz.multiply(weights);
381 final boolean hsig = updateEvolutionPaths(zmean, xold);
382 if (diagonalOnly <= 0) {
383 updateCovariance(hsig, bestArx, arz, arindex, xold);
384 } else {
385 updateCovarianceDiagonalOnly(hsig, bestArz);
386 }
387
388 sigma *= JdkMath.exp(JdkMath.min(1, (normps/chiN - 1) * cs / damps));
389 final double bestFitness = fitness[arindex[0]];
390 final double worstFitness = fitness[arindex[arindex.length - 1]];
391 if (bestValue > bestFitness) {
392 bestValue = bestFitness;
393 lastResult = optimum;
394 optimum = new PointValuePair(fitfun.repair(bestArx.getColumn(0)),
395 isMinimize ? bestFitness : -bestFitness);
396 if (getConvergenceChecker() != null && lastResult != null &&
397 getConvergenceChecker().converged(iterations, optimum, lastResult)) {
398 break generationLoop;
399 }
400 }
401
402
403 if (stopFitness != 0 && bestFitness < (isMinimize ? stopFitness : -stopFitness)) {
404 break generationLoop;
405 }
406 final double[] sqrtDiagC = sqrt(diagC).getColumn(0);
407 final double[] pcCol = pc.getColumn(0);
408 for (int i = 0; i < dimension; i++) {
409 if (sigma * JdkMath.max(JdkMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) {
410 break;
411 }
412 if (i >= dimension - 1) {
413 break generationLoop;
414 }
415 }
416 for (int i = 0; i < dimension; i++) {
417 if (sigma * sqrtDiagC[i] > stopTolUpX) {
418 break generationLoop;
419 }
420 }
421 final double historyBest = min(fitnessHistory);
422 final double historyWorst = max(fitnessHistory);
423 if (iterations > 2 &&
424 JdkMath.max(historyWorst, worstFitness) -
425 JdkMath.min(historyBest, bestFitness) < stopTolFun) {
426 break generationLoop;
427 }
428 if (iterations > fitnessHistory.length &&
429 historyWorst - historyBest < stopTolHistFun) {
430 break generationLoop;
431 }
432
433 if (max(diagD) / min(diagD) > 1e7) {
434 break generationLoop;
435 }
436
437 if (getConvergenceChecker() != null) {
438 final PointValuePair current
439 = new PointValuePair(bestArx.getColumn(0),
440 isMinimize ? bestFitness : -bestFitness);
441 if (lastResult != null &&
442 getConvergenceChecker().converged(iterations, current, lastResult)) {
443 break generationLoop;
444 }
445 lastResult = current;
446 }
447
448 if (bestValue == fitness[arindex[(int)(0.1+lambda/4.)]]) {
449 sigma *= JdkMath.exp(0.2 + cs / damps);
450 }
451 if (iterations > 2 && JdkMath.max(historyWorst, bestFitness) -
452 JdkMath.min(historyBest, bestFitness) == 0) {
453 sigma *= JdkMath.exp(0.2 + cs / damps);
454 }
455
456 push(fitnessHistory,bestFitness);
457 if (generateStatistics) {
458 statisticsSigmaHistory.add(sigma);
459 statisticsFitnessHistory.add(bestFitness);
460 statisticsMeanHistory.add(xmean.transpose());
461 statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5));
462 }
463 }
464 return optimum;
465 }
466
467
468
469
470
471
472
473
474
475
476
477 @Override
478 protected void parseOptimizationData(OptimizationData... optData) {
479
480 super.parseOptimizationData(optData);
481
482
483
484 for (OptimizationData data : optData) {
485 if (data instanceof Sigma) {
486 inputSigma = ((Sigma) data).getSigma();
487 continue;
488 }
489 if (data instanceof PopulationSize) {
490 lambda = ((PopulationSize) data).getPopulationSize();
491 continue;
492 }
493 }
494
495 checkParameters();
496 }
497
498
499
500
501 private void checkParameters() {
502 if (inputSigma != null) {
503 final double[] init = getStartPoint();
504
505 if (inputSigma.length != init.length) {
506 throw new DimensionMismatchException(inputSigma.length, init.length);
507 }
508
509 final double[] lB = getLowerBound();
510 final double[] uB = getUpperBound();
511
512 for (int i = 0; i < init.length; i++) {
513 if (inputSigma[i] > uB[i] - lB[i]) {
514 throw new OutOfRangeException(inputSigma[i], 0, uB[i] - lB[i]);
515 }
516 }
517 }
518 }
519
520
521
522
523
524
525 private void initializeCMA(double[] guess) {
526 if (lambda <= 0) {
527 throw new NotStrictlyPositiveException(lambda);
528 }
529
530 final double[][] sigmaArray = new double[guess.length][1];
531 for (int i = 0; i < guess.length; i++) {
532 sigmaArray[i][0] = inputSigma[i];
533 }
534 final RealMatrix insigma = new Array2DRowRealMatrix(sigmaArray, false);
535 sigma = max(insigma);
536
537
538 stopTolUpX = 1e3 * max(insigma);
539 stopTolX = 1e-11 * max(insigma);
540 stopTolFun = 1e-12;
541 stopTolHistFun = 1e-13;
542
543
544 mu = lambda / 2;
545 logMu2 = JdkMath.log(mu + 0.5);
546 weights = log(sequence(1, mu, 1)).scalarMultiply(-1).scalarAdd(logMu2);
547 double sumw = 0;
548 double sumwq = 0;
549 for (int i = 0; i < mu; i++) {
550 double w = weights.getEntry(i, 0);
551 sumw += w;
552 sumwq += w * w;
553 }
554 weights = weights.scalarMultiply(1 / sumw);
555 mueff = sumw * sumw / sumwq;
556
557
558 cc = (4 + mueff / dimension) /
559 (dimension + 4 + 2 * mueff / dimension);
560 cs = (mueff + 2) / (dimension + mueff + 3.);
561 damps = (1 + 2 * JdkMath.max(0, JdkMath.sqrt((mueff - 1) /
562 (dimension + 1)) - 1)) *
563 JdkMath.max(0.3,
564 1 - dimension / (1e-6 + maxIterations)) + cs;
565 ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
566 ccovmu = JdkMath.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) /
567 ((dimension + 2) * (dimension + 2) + mueff));
568 ccov1Sep = JdkMath.min(1, ccov1 * (dimension + 1.5) / 3);
569 ccovmuSep = JdkMath.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
570 chiN = JdkMath.sqrt(dimension) *
571 (1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension));
572
573 xmean = MatrixUtils.createColumnRealMatrix(guess);
574 diagD = insigma.scalarMultiply(1 / sigma);
575 diagC = square(diagD);
576 pc = zeros(dimension, 1);
577 ps = zeros(dimension, 1);
578 normps = ps.getFrobeniusNorm();
579
580 B = eye(dimension, dimension);
581 D = ones(dimension, 1);
582 BD = times(B, repmat(diagD.transpose(), dimension, 1));
583 C = B.multiply(diag(square(D)).multiply(B.transpose()));
584 historySize = 10 + (int) (3 * 10 * dimension / (double) lambda);
585 fitnessHistory = new double[historySize];
586 for (int i = 0; i < historySize; i++) {
587 fitnessHistory[i] = Double.MAX_VALUE;
588 }
589 }
590
591
592
593
594
595
596
597
598
599 private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) {
600 ps = ps.scalarMultiply(1 - cs).add(
601 B.multiply(zmean).scalarMultiply(
602 JdkMath.sqrt(cs * (2 - cs) * mueff)));
603 normps = ps.getFrobeniusNorm();
604 final boolean hsig = normps /
605 JdkMath.sqrt(1 - JdkMath.pow(1 - cs, 2 * iterations)) /
606 chiN < 1.4 + 2 / ((double) dimension + 1);
607 pc = pc.scalarMultiply(1 - cc);
608 if (hsig) {
609 pc = pc.add(xmean.subtract(xold).scalarMultiply(JdkMath.sqrt(cc * (2 - cc) * mueff) / sigma));
610 }
611 return hsig;
612 }
613
614
615
616
617
618
619
620
621 private void updateCovarianceDiagonalOnly(boolean hsig,
622 final RealMatrix bestArz) {
623
624 double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc);
625 oldFac += 1 - ccov1Sep - ccovmuSep;
626 diagC = diagC.scalarMultiply(oldFac)
627 .add(square(pc).scalarMultiply(ccov1Sep))
628 .add(times(diagC, square(bestArz).multiply(weights))
629 .scalarMultiply(ccovmuSep));
630 diagD = sqrt(diagC);
631 if (diagonalOnly > 1 &&
632 iterations > diagonalOnly) {
633
634 diagonalOnly = 0;
635 B = eye(dimension, dimension);
636 BD = diag(diagD);
637 C = diag(diagC);
638 }
639 }
640
641
642
643
644
645
646
647
648
649
650
651
652 private void updateCovariance(boolean hsig, final RealMatrix bestArx,
653 final RealMatrix arz, final int[] arindex,
654 final RealMatrix xold) {
655 double negccov = 0;
656 if (ccov1 + ccovmu > 0) {
657 final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu))
658 .scalarMultiply(1 / sigma);
659 final RealMatrix roneu = pc.multiply(pc.transpose())
660 .scalarMultiply(ccov1);
661
662 double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc);
663 oldFac += 1 - ccov1 - ccovmu;
664 if (isActiveCMA) {
665
666 negccov = (1 - ccovmu) * 0.25 * mueff /
667 (JdkMath.pow(dimension + 2.0, 1.5) + 2 * mueff);
668
669
670 final double negminresidualvariance = 0.66;
671
672 final double negalphaold = 0.5;
673
674 final int[] arReverseIndex = reverse(arindex);
675 RealMatrix arzneg = selectColumns(arz, Arrays.copyOf(arReverseIndex, mu));
676 RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
677 final int[] idxnorms = sortedIndices(arnorms.getRow(0));
678 final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);
679 final int[] idxReverse = reverse(idxnorms);
680 final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse);
681 arnorms = divide(arnormsReverse, arnormsSorted);
682 final int[] idxInv = inverse(idxnorms);
683 final RealMatrix arnormsInv = selectColumns(arnorms, idxInv);
684
685 final double negcovMax = (1 - negminresidualvariance) /
686 square(arnormsInv).multiply(weights).getEntry(0, 0);
687 if (negccov > negcovMax) {
688 negccov = negcovMax;
689 }
690 arzneg = times(arzneg, repmat(arnormsInv, dimension, 1));
691 final RealMatrix artmp = BD.multiply(arzneg);
692 final RealMatrix cNeg = artmp.multiply(diag(weights)).multiply(artmp.transpose());
693 oldFac += negalphaold * negccov;
694 C = C.scalarMultiply(oldFac)
695 .add(roneu)
696 .add(arpos.scalarMultiply(
697 ccovmu + (1 - negalphaold) * negccov)
698 .multiply(times(repmat(weights, 1, dimension),
699 arpos.transpose())))
700 .subtract(cNeg.scalarMultiply(negccov));
701 } else {
702
703 C = C.scalarMultiply(oldFac)
704 .add(roneu)
705 .add(arpos.scalarMultiply(ccovmu)
706 .multiply(times(repmat(weights, 1, dimension),
707 arpos.transpose())));
708 }
709 }
710 updateBD(negccov);
711 }
712
713
714
715
716
717
718 private void updateBD(double negccov) {
719 if (ccov1 + ccovmu + negccov > 0 &&
720 (iterations % 1. / (ccov1 + ccovmu + negccov) / dimension / 10.) < 1) {
721
722 C = triu(C, 0).add(triu(C, 1).transpose());
723
724 final EigenDecomposition eig = new EigenDecomposition(C);
725 B = eig.getV();
726 D = eig.getD();
727 diagD = diag(D);
728 if (min(diagD) <= 0) {
729 for (int i = 0; i < dimension; i++) {
730 if (diagD.getEntry(i, 0) < 0) {
731 diagD.setEntry(i, 0, 0);
732 }
733 }
734 final double tfac = max(diagD) / 1e14;
735 C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
736 diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
737 }
738 if (max(diagD) > 1e14 * min(diagD)) {
739 final double tfac = max(diagD) / 1e14 - min(diagD);
740 C = C.add(eye(dimension, dimension).scalarMultiply(tfac));
741 diagD = diagD.add(ones(dimension, 1).scalarMultiply(tfac));
742 }
743 diagC = diag(C);
744 diagD = sqrt(diagD);
745 BD = times(B, repmat(diagD.transpose(), dimension, 1));
746 }
747 }
748
749
750
751
752
753
754
755 private static void push(double[] vals, double val) {
756 if (vals.length - 1 > 0) {
757 System.arraycopy(vals, 0, vals, 1, vals.length - 1);
758 }
759 vals[0] = val;
760 }
761
762
763
764
765
766
767
768 private int[] sortedIndices(final double[] doubles) {
769 final DoubleIndex[] dis = new DoubleIndex[doubles.length];
770 for (int i = 0; i < doubles.length; i++) {
771 dis[i] = new DoubleIndex(doubles[i], i);
772 }
773 Arrays.sort(dis);
774 final int[] indices = new int[doubles.length];
775 for (int i = 0; i < doubles.length; i++) {
776 indices[i] = dis[i].index;
777 }
778 return indices;
779 }
780
781
782
783
784
785
786 private double valueRange(final ValuePenaltyPair[] vpPairs) {
787 double max = Double.NEGATIVE_INFINITY;
788 double min = Double.MAX_VALUE;
789 for (ValuePenaltyPair vpPair:vpPairs) {
790 if (vpPair.value > max) {
791 max = vpPair.value;
792 }
793 if (vpPair.value < min) {
794 min = vpPair.value;
795 }
796 }
797 return max-min;
798 }
799
800
801
802
803
804 private static final class DoubleIndex implements Comparable<DoubleIndex> {
805
806 private final double value;
807
808 private final int index;
809
810
811
812
813
814 DoubleIndex(double value, int index) {
815 this.value = value;
816 this.index = index;
817 }
818
819
820 @Override
821 public int compareTo(DoubleIndex o) {
822 return Double.compare(value, o.value);
823 }
824
825
826 @Override
827 public boolean equals(Object other) {
828
829 if (this == other) {
830 return true;
831 }
832
833 if (other instanceof DoubleIndex) {
834 return Double.compare(value, ((DoubleIndex) other).value) == 0;
835 }
836
837 return false;
838 }
839
840
841 @Override
842 public int hashCode() {
843 long bits = Double.doubleToLongBits(value);
844 return (int) (1438542 ^ (bits >>> 32) ^ bits);
845 }
846 }
847
848
849
850 private static final class ValuePenaltyPair {
851
852 private double value;
853
854 private double penalty;
855
856
857
858
859
860 ValuePenaltyPair(final double value, final double penalty) {
861 this.value = value;
862 this.penalty = penalty;
863 }
864 }
865
866
867
868
869
870
871 private final class FitnessFunction {
872
873
874
875
876 private final boolean isRepairMode;
877
878
879
880 FitnessFunction() {
881 isRepairMode = true;
882 }
883
884
885
886
887
888 public ValuePenaltyPair value(final double[] point) {
889 final MultivariateFunction func = CMAESOptimizer.this.getObjectiveFunction();
890 double value;
891 double penalty = 0;
892 if (isRepairMode) {
893 double[] repaired = repair(point);
894 value = func.value(repaired);
895 penalty = penalty(point, repaired);
896 } else {
897 value = func.value(point);
898 }
899 value = isMinimize ? value : -value;
900 penalty = isMinimize ? penalty : -penalty;
901 return new ValuePenaltyPair(value,penalty);
902 }
903
904
905
906
907
908 public boolean isFeasible(final double[] x) {
909 final double[] lB = CMAESOptimizer.this.getLowerBound();
910 final double[] uB = CMAESOptimizer.this.getUpperBound();
911
912 for (int i = 0; i < x.length; i++) {
913 if (x[i] < lB[i]) {
914 return false;
915 }
916 if (x[i] > uB[i]) {
917 return false;
918 }
919 }
920 return true;
921 }
922
923
924
925
926
927 private double[] repair(final double[] x) {
928 final double[] lB = CMAESOptimizer.this.getLowerBound();
929 final double[] uB = CMAESOptimizer.this.getUpperBound();
930
931 final double[] repaired = new double[x.length];
932 for (int i = 0; i < x.length; i++) {
933 if (x[i] < lB[i]) {
934 repaired[i] = lB[i];
935 } else if (x[i] > uB[i]) {
936 repaired[i] = uB[i];
937 } else {
938 repaired[i] = x[i];
939 }
940 }
941 return repaired;
942 }
943
944
945
946
947
948
949 private double penalty(final double[] x, final double[] repaired) {
950 double penalty = 0;
951 for (int i = 0; i < x.length; i++) {
952 double diff = JdkMath.abs(x[i] - repaired[i]);
953 penalty += diff;
954 }
955 return isMinimize ? penalty : -penalty;
956 }
957 }
958
959
960
961
962
963
964
965 private static RealMatrix log(final RealMatrix m) {
966 final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
967 for (int r = 0; r < m.getRowDimension(); r++) {
968 for (int c = 0; c < m.getColumnDimension(); c++) {
969 d[r][c] = JdkMath.log(m.getEntry(r, c));
970 }
971 }
972 return new Array2DRowRealMatrix(d, false);
973 }
974
975
976
977
978
979 private static RealMatrix sqrt(final RealMatrix m) {
980 final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
981 for (int r = 0; r < m.getRowDimension(); r++) {
982 for (int c = 0; c < m.getColumnDimension(); c++) {
983 d[r][c] = JdkMath.sqrt(m.getEntry(r, c));
984 }
985 }
986 return new Array2DRowRealMatrix(d, false);
987 }
988
989
990
991
992
993 private static RealMatrix square(final RealMatrix m) {
994 final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
995 for (int r = 0; r < m.getRowDimension(); r++) {
996 for (int c = 0; c < m.getColumnDimension(); c++) {
997 double e = m.getEntry(r, c);
998 d[r][c] = e * e;
999 }
1000 }
1001 return new Array2DRowRealMatrix(d, false);
1002 }
1003
1004
1005
1006
1007
1008
1009 private static RealMatrix times(final RealMatrix m, final RealMatrix n) {
1010 final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1011 for (int r = 0; r < m.getRowDimension(); r++) {
1012 for (int c = 0; c < m.getColumnDimension(); c++) {
1013 d[r][c] = m.getEntry(r, c) * n.getEntry(r, c);
1014 }
1015 }
1016 return new Array2DRowRealMatrix(d, false);
1017 }
1018
1019
1020
1021
1022
1023
1024 private static RealMatrix divide(final RealMatrix m, final RealMatrix n) {
1025 final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1026 for (int r = 0; r < m.getRowDimension(); r++) {
1027 for (int c = 0; c < m.getColumnDimension(); c++) {
1028 d[r][c] = m.getEntry(r, c) / n.getEntry(r, c);
1029 }
1030 }
1031 return new Array2DRowRealMatrix(d, false);
1032 }
1033
1034
1035
1036
1037
1038
1039 private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) {
1040 final double[][] d = new double[m.getRowDimension()][cols.length];
1041 for (int r = 0; r < m.getRowDimension(); r++) {
1042 for (int c = 0; c < cols.length; c++) {
1043 d[r][c] = m.getEntry(r, cols[c]);
1044 }
1045 }
1046 return new Array2DRowRealMatrix(d, false);
1047 }
1048
1049
1050
1051
1052
1053
1054 private static RealMatrix triu(final RealMatrix m, int k) {
1055 final double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
1056 for (int r = 0; r < m.getRowDimension(); r++) {
1057 for (int c = 0; c < m.getColumnDimension(); c++) {
1058 d[r][c] = r <= c - k ? m.getEntry(r, c) : 0;
1059 }
1060 }
1061 return new Array2DRowRealMatrix(d, false);
1062 }
1063
1064
1065
1066
1067
1068 private static RealMatrix sumRows(final RealMatrix m) {
1069 final double[][] d = new double[1][m.getColumnDimension()];
1070 for (int c = 0; c < m.getColumnDimension(); c++) {
1071 double sum = 0;
1072 for (int r = 0; r < m.getRowDimension(); r++) {
1073 sum += m.getEntry(r, c);
1074 }
1075 d[0][c] = sum;
1076 }
1077 return new Array2DRowRealMatrix(d, false);
1078 }
1079
1080
1081
1082
1083
1084
1085 private static RealMatrix diag(final RealMatrix m) {
1086 if (m.getColumnDimension() == 1) {
1087 final double[][] d = new double[m.getRowDimension()][m.getRowDimension()];
1088 for (int i = 0; i < m.getRowDimension(); i++) {
1089 d[i][i] = m.getEntry(i, 0);
1090 }
1091 return new Array2DRowRealMatrix(d, false);
1092 } else {
1093 final double[][] d = new double[m.getRowDimension()][1];
1094 for (int i = 0; i < m.getColumnDimension(); i++) {
1095 d[i][0] = m.getEntry(i, i);
1096 }
1097 return new Array2DRowRealMatrix(d, false);
1098 }
1099 }
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109 private static void copyColumn(final RealMatrix m1, int col1,
1110 RealMatrix m2, int col2) {
1111 for (int i = 0; i < m1.getRowDimension(); i++) {
1112 m2.setEntry(i, col2, m1.getEntry(i, col1));
1113 }
1114 }
1115
1116
1117
1118
1119
1120
1121 private static RealMatrix ones(int n, int m) {
1122 final double[][] d = new double[n][m];
1123 for (int r = 0; r < n; r++) {
1124 Arrays.fill(d[r], 1);
1125 }
1126 return new Array2DRowRealMatrix(d, false);
1127 }
1128
1129
1130
1131
1132
1133
1134
1135 private static RealMatrix eye(int n, int m) {
1136 final double[][] d = new double[n][m];
1137 for (int r = 0; r < n; r++) {
1138 if (r < m) {
1139 d[r][r] = 1;
1140 }
1141 }
1142 return new Array2DRowRealMatrix(d, false);
1143 }
1144
1145
1146
1147
1148
1149
1150 private static RealMatrix zeros(int n, int m) {
1151 return new Array2DRowRealMatrix(n, m);
1152 }
1153
1154
1155
1156
1157
1158
1159
1160 private static RealMatrix repmat(final RealMatrix mat, int n, int m) {
1161 final int rd = mat.getRowDimension();
1162 final int cd = mat.getColumnDimension();
1163 final double[][] d = new double[n * rd][m * cd];
1164 for (int r = 0; r < n * rd; r++) {
1165 for (int c = 0; c < m * cd; c++) {
1166 d[r][c] = mat.getEntry(r % rd, c % cd);
1167 }
1168 }
1169 return new Array2DRowRealMatrix(d, false);
1170 }
1171
1172
1173
1174
1175
1176
1177
1178 private static RealMatrix sequence(double start, double end, double step) {
1179 final int size = (int) ((end - start) / step + 1);
1180 final double[][] d = new double[size][1];
1181 double value = start;
1182 for (int r = 0; r < size; r++) {
1183 d[r][0] = value;
1184 value += step;
1185 }
1186 return new Array2DRowRealMatrix(d, false);
1187 }
1188
1189
1190
1191
1192
1193 private static double max(final RealMatrix m) {
1194 double max = -Double.MAX_VALUE;
1195 for (int r = 0; r < m.getRowDimension(); r++) {
1196 for (int c = 0; c < m.getColumnDimension(); c++) {
1197 double e = m.getEntry(r, c);
1198 if (max < e) {
1199 max = e;
1200 }
1201 }
1202 }
1203 return max;
1204 }
1205
1206
1207
1208
1209
1210 private static double min(final RealMatrix m) {
1211 double min = Double.MAX_VALUE;
1212 for (int r = 0; r < m.getRowDimension(); r++) {
1213 for (int c = 0; c < m.getColumnDimension(); c++) {
1214 double e = m.getEntry(r, c);
1215 if (min > e) {
1216 min = e;
1217 }
1218 }
1219 }
1220 return min;
1221 }
1222
1223
1224
1225
1226
1227 private static double max(final double[] m) {
1228 double max = -Double.MAX_VALUE;
1229 for (int r = 0; r < m.length; r++) {
1230 if (max < m[r]) {
1231 max = m[r];
1232 }
1233 }
1234 return max;
1235 }
1236
1237
1238
1239
1240
1241 private static double min(final double[] m) {
1242 double min = Double.MAX_VALUE;
1243 for (int r = 0; r < m.length; r++) {
1244 if (min > m[r]) {
1245 min = m[r];
1246 }
1247 }
1248 return min;
1249 }
1250
1251
1252
1253
1254
1255 private static int[] inverse(final int[] indices) {
1256 final int[] inverse = new int[indices.length];
1257 for (int i = 0; i < indices.length; i++) {
1258 inverse[indices[i]] = i;
1259 }
1260 return inverse;
1261 }
1262
1263
1264
1265
1266
1267 private static int[] reverse(final int[] indices) {
1268 final int[] reverse = new int[indices.length];
1269 for (int i = 0; i < indices.length; i++) {
1270 reverse[i] = indices[indices.length - i - 1];
1271 }
1272 return reverse;
1273 }
1274
1275
1276
1277
1278
1279 private double[] randn(int size) {
1280 final double[] randn = new double[size];
1281 for (int i = 0; i < size; i++) {
1282 randn[i] = random.sample();
1283 }
1284 return randn;
1285 }
1286
1287
1288
1289
1290
1291
1292 private RealMatrix randn1(int size, int popSize) {
1293 final double[][] d = new double[size][popSize];
1294 for (int r = 0; r < size; r++) {
1295 for (int c = 0; c < popSize; c++) {
1296 d[r][c] = random.sample();
1297 }
1298 }
1299 return new Array2DRowRealMatrix(d, false);
1300 }
1301 }