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