View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.optim.nonlinear.scalar.noderiv;
23  
24  import org.hipparchus.analysis.MultivariateFunction;
25  import org.hipparchus.exception.MathIllegalArgumentException;
26  import org.hipparchus.optim.InitialGuess;
27  import org.hipparchus.optim.MaxEval;
28  import org.hipparchus.optim.PointValuePair;
29  import org.hipparchus.optim.SimpleBounds;
30  import org.hipparchus.optim.nonlinear.scalar.GoalType;
31  import org.hipparchus.optim.nonlinear.scalar.ObjectiveFunction;
32  import org.hipparchus.random.MersenneTwister;
33  import org.hipparchus.util.FastMath;
34  import org.junit.jupiter.api.Test;
35  import org.junitpioneer.jupiter.RetryingTest;
36  
37  import java.util.Arrays;
38  import java.util.Random;
39  
40  import static org.junit.jupiter.api.Assertions.assertEquals;
41  import static org.junit.jupiter.api.Assertions.assertThrows;
42  import static org.junit.jupiter.api.Assertions.assertTrue;
43  
44  /**
45   * Test for {@link CMAESOptimizer}.
46   */
47  public class CMAESOptimizerTest {
48  
49      static final int DIM = 13;
50      static final int LAMBDA = 4 + (int)(3.*FastMath.log(DIM));
51  
52      @Test
53      void testInitOutofbounds1() {
54          assertThrows(MathIllegalArgumentException.class, () -> {
55              double[] startPoint = point(DIM, 3);
56              double[] insigma = point(DIM, 0.3);
57              double[][] boundaries = boundaries(DIM, -1, 2);
58              PointValuePair expected =
59                  new PointValuePair(point(DIM, 1.0), 0.0);
60              doTest(new Rosen(), startPoint, insigma, boundaries,
61                  GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
62                  1e-13, 1e-6, 100000, expected);
63          });
64      }
65  
66      @Test
67      void testInitOutofbounds2() {
68          assertThrows(MathIllegalArgumentException.class, () -> {
69              double[] startPoint = point(DIM, -2);
70              double[] insigma = point(DIM, 0.3);
71              double[][] boundaries = boundaries(DIM, -1, 2);
72              PointValuePair expected =
73                  new PointValuePair(point(DIM, 1.0), 0.0);
74              doTest(new Rosen(), startPoint, insigma, boundaries,
75                  GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
76                  1e-13, 1e-6, 100000, expected);
77          });
78      }
79  
80      @Test
81      void testBoundariesDimensionMismatch() {
82          assertThrows(MathIllegalArgumentException.class, () -> {
83              double[] startPoint = point(DIM, 0.5);
84              double[] insigma = point(DIM, 0.3);
85              double[][] boundaries = boundaries(DIM + 1, -1, 2);
86              PointValuePair expected =
87                  new PointValuePair(point(DIM, 1.0), 0.0);
88              doTest(new Rosen(), startPoint, insigma, boundaries,
89                  GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
90                  1e-13, 1e-6, 100000, expected);
91          });
92      }
93  
94      @Test
95      void testInputSigmaNegative() {
96          assertThrows(MathIllegalArgumentException.class, () -> {
97              double[] startPoint = point(DIM, 0.5);
98              double[] insigma = point(DIM, -0.5);
99              double[][] boundaries = null;
100             PointValuePair expected =
101                 new PointValuePair(point(DIM, 1.0), 0.0);
102             doTest(new Rosen(), startPoint, insigma, boundaries,
103                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
104                 1e-13, 1e-6, 100000, expected);
105         });
106     }
107 
108     @Test
109     void testInputSigmaOutOfRange() {
110         assertThrows(MathIllegalArgumentException.class, () -> {
111             double[] startPoint = point(DIM, 0.5);
112             double[] insigma = point(DIM, 1.1);
113             double[][] boundaries = boundaries(DIM, -0.5, 0.5);
114             PointValuePair expected =
115                 new PointValuePair(point(DIM, 1.0), 0.0);
116             doTest(new Rosen(), startPoint, insigma, boundaries,
117                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
118                 1e-13, 1e-6, 100000, expected);
119         });
120     }
121 
122     @Test
123     void testInputSigmaDimensionMismatch() {
124         assertThrows(MathIllegalArgumentException.class, () -> {
125             double[] startPoint = point(DIM, 0.5);
126             double[] insigma = point(DIM + 1, 0.5);
127             double[][] boundaries = null;
128             PointValuePair expected =
129                 new PointValuePair(point(DIM, 1.0), 0.0);
130             doTest(new Rosen(), startPoint, insigma, boundaries,
131                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
132                 1e-13, 1e-6, 100000, expected);
133         });
134     }
135 
136     @RetryingTest(3)
137     public void testRosen() {
138         double[] startPoint = point(DIM,0.1);
139         double[] insigma = point(DIM,0.1);
140         double[][] boundaries = null;
141         PointValuePair expected =
142             new PointValuePair(point(DIM,1.0),0.0);
143         doTest(new Rosen(), startPoint, insigma, boundaries,
144                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
145                 1e-13, 1e-6, 100000, expected);
146         doTest(new Rosen(), startPoint, insigma, boundaries,
147                 GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
148                 1e-13, 1e-6, 100000, expected);
149     }
150 
151 
152     @RetryingTest(3)
153     void testMaximize() {
154         double[] startPoint = point(DIM,1.0);
155         double[] insigma = point(DIM,0.1);
156         double[][] boundaries = null;
157         PointValuePair expected =
158             new PointValuePair(point(DIM,0.0),1.0);
159         doTest(new MinusElli(), startPoint, insigma, boundaries,
160                 GoalType.MAXIMIZE, LAMBDA, true, 0, 1.0-1e-13,
161                 2e-10, 5e-6, 100000, expected);
162         doTest(new MinusElli(), startPoint, insigma, boundaries,
163                 GoalType.MAXIMIZE, LAMBDA, false, 0, 1.0-1e-13,
164                 2e-10, 5e-6, 100000, expected);
165         boundaries = boundaries(DIM,-0.3,0.3);
166         startPoint = point(DIM,0.1);
167         doTest(new MinusElli(), startPoint, insigma, boundaries,
168                 GoalType.MAXIMIZE, LAMBDA, true, 0, 1.0-1e-13,
169                 2e-10, 5e-6, 100000, expected);
170     }
171 
172     @Test
173     void testEllipse() {
174         double[] startPoint = point(DIM,1.0);
175         double[] insigma = point(DIM,0.1);
176         double[][] boundaries = null;
177         PointValuePair expected =
178             new PointValuePair(point(DIM,0.0),0.0);
179         doTest(new Elli(), startPoint, insigma, boundaries,
180                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
181                 1e-13, 1e-6, 100000, expected);
182         doTest(new Elli(), startPoint, insigma, boundaries,
183                 GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
184                 1e-13, 1e-6, 100000, expected);
185     }
186 
187     @Test
188     void testElliRotated() {
189         double[] startPoint = point(DIM,1.0);
190         double[] insigma = point(DIM,0.1);
191         double[][] boundaries = null;
192         PointValuePair expected =
193             new PointValuePair(point(DIM,0.0),0.0);
194         doTest(new ElliRotated(), startPoint, insigma, boundaries,
195                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
196                 1e-13, 1e-6, 100000, expected);
197         doTest(new ElliRotated(), startPoint, insigma, boundaries,
198                 GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
199                 1e-13, 1e-6, 100000, expected);
200     }
201 
202     @Test
203     void testCigar() {
204         double[] startPoint = point(DIM,1.0);
205         double[] insigma = point(DIM,0.1);
206         double[][] boundaries = null;
207         PointValuePair expected =
208             new PointValuePair(point(DIM,0.0),0.0);
209         doTest(new Cigar(), startPoint, insigma, boundaries,
210                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
211                 1e-13, 1e-6, 200000, expected);
212         doTest(new Cigar(), startPoint, insigma, boundaries,
213                 GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
214                 1e-13, 1e-6, 100000, expected);
215     }
216 
217     @Test
218     void testCigarWithBoundaries() {
219         double[] startPoint = point(DIM,1.0);
220         double[] insigma = point(DIM,0.1);
221         double[][] boundaries = boundaries(DIM, -1e100, Double.POSITIVE_INFINITY);
222         PointValuePair expected =
223             new PointValuePair(point(DIM,0.0),0.0);
224         doTest(new Cigar(), startPoint, insigma, boundaries,
225                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
226                 1e-13, 1e-6, 200000, expected);
227         doTest(new Cigar(), startPoint, insigma, boundaries,
228                 GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
229                 1e-13, 1e-6, 100000, expected);
230     }
231 
232     @Test
233     void testTwoAxes() {
234         double[] startPoint = point(DIM,1.0);
235         double[] insigma = point(DIM,0.1);
236         double[][] boundaries = null;
237         PointValuePair expected =
238             new PointValuePair(point(DIM,0.0),0.0);
239         doTest(new TwoAxes(), startPoint, insigma, boundaries,
240                 GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
241                 1e-13, 1e-6, 200000, expected);
242         doTest(new TwoAxes(), startPoint, insigma, boundaries,
243                 GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
244                 1e-8, 1e-3, 200000, expected);
245     }
246 
247     @Test
248     void testCigTab() {
249         double[] startPoint = point(DIM,1.0);
250         double[] insigma = point(DIM,0.3);
251         double[][] boundaries = null;
252         PointValuePair expected =
253             new PointValuePair(point(DIM,0.0),0.0);
254         doTest(new CigTab(), startPoint, insigma, boundaries,
255                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
256                 1e-13, 5e-5, 100000, expected);
257         doTest(new CigTab(), startPoint, insigma, boundaries,
258                 GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
259                 1e-13, 5e-5, 100000, expected);
260     }
261 
262     @Test
263     void testSphere() {
264         double[] startPoint = point(DIM,1.0);
265         double[] insigma = point(DIM,0.1);
266         double[][] boundaries = null;
267         PointValuePair expected =
268             new PointValuePair(point(DIM,0.0),0.0);
269         doTest(new Sphere(), startPoint, insigma, boundaries,
270                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
271                 1e-13, 1e-6, 100000, expected);
272         doTest(new Sphere(), startPoint, insigma, boundaries,
273                 GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
274                 1e-13, 1e-6, 100000, expected);
275     }
276 
277     @Test
278     void testTablet() {
279         double[] startPoint = point(DIM,1.0);
280         double[] insigma = point(DIM,0.1);
281         double[][] boundaries = null;
282         PointValuePair expected =
283             new PointValuePair(point(DIM,0.0),0.0);
284         doTest(new Tablet(), startPoint, insigma, boundaries,
285                 GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
286                 1e-13, 1e-6, 100000, expected);
287         doTest(new Tablet(), startPoint, insigma, boundaries,
288                 GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
289                 1e-13, 1e-6, 100000, expected);
290     }
291 
292     @Test
293     void testDiffPow() {
294         double[] startPoint = point(DIM,1.0);
295         double[] insigma = point(DIM,0.1);
296         double[][] boundaries = null;
297         PointValuePair expected =
298             new PointValuePair(point(DIM,0.0),0.0);
299         doTest(new DiffPow(), startPoint, insigma, boundaries,
300                 GoalType.MINIMIZE, 10, true, 0, 1e-13,
301                 1e-8, 1e-1, 100000, expected);
302         doTest(new DiffPow(), startPoint, insigma, boundaries,
303                 GoalType.MINIMIZE, 10, false, 0, 1e-13,
304                 1e-8, 2e-1, 100000, expected);
305     }
306 
307     @Test
308     void testSsDiffPow() {
309         double[] startPoint = point(DIM,1.0);
310         double[] insigma = point(DIM,0.1);
311         double[][] boundaries = null;
312         PointValuePair expected =
313             new PointValuePair(point(DIM,0.0),0.0);
314         doTest(new SsDiffPow(), startPoint, insigma, boundaries,
315                 GoalType.MINIMIZE, 10, true, 0, 1e-13,
316                 1e-4, 1e-1, 200000, expected);
317         doTest(new SsDiffPow(), startPoint, insigma, boundaries,
318                 GoalType.MINIMIZE, 10, false, 0, 1e-13,
319                 1e-4, 1e-1, 200000, expected);
320     }
321 
322     @Test
323     void testAckley() {
324         double[] startPoint = point(DIM,1.0);
325         double[] insigma = point(DIM,1.0);
326         double[][] boundaries = null;
327         PointValuePair expected =
328             new PointValuePair(point(DIM,0.0),0.0);
329         doTest(new Ackley(), startPoint, insigma, boundaries,
330                 GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
331                 1e-9, 1e-5, 100000, expected);
332         doTest(new Ackley(), startPoint, insigma, boundaries,
333                 GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
334                 1e-9, 1e-5, 100000, expected);
335     }
336 
337     @Test
338     void testRastrigin() {
339         double[] startPoint = point(DIM,0.1);
340         double[] insigma = point(DIM,0.1);
341         double[][] boundaries = null;
342         PointValuePair expected =
343             new PointValuePair(point(DIM,0.0),0.0);
344         doTest(new Rastrigin(), startPoint, insigma, boundaries,
345                 GoalType.MINIMIZE, (int)(200*FastMath.sqrt(DIM)), true, 0, 1e-13,
346                 1e-13, 1e-6, 200000, expected);
347         doTest(new Rastrigin(), startPoint, insigma, boundaries,
348                 GoalType.MINIMIZE, (int)(200*FastMath.sqrt(DIM)), false, 0, 1e-13,
349                 1e-13, 1e-6, 200000, expected);
350     }
351 
352     @Test
353     void testConstrainedRosen() {
354         double[] startPoint = point(DIM, 0.1);
355         double[] insigma = point(DIM, 0.1);
356         double[][] boundaries = boundaries(DIM, -1, 2);
357         PointValuePair expected =
358             new PointValuePair(point(DIM,1.0),0.0);
359         doTest(new Rosen(), startPoint, insigma, boundaries,
360                 GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
361                 1e-13, 1e-6, 100000, expected);
362         doTest(new Rosen(), startPoint, insigma, boundaries,
363                 GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
364                 1e-13, 1e-6, 100000, expected);
365     }
366 
367     @Test
368     void testDiagonalRosen() {
369         double[] startPoint = point(DIM,0.1);
370         double[] insigma = point(DIM,0.1);
371         double[][] boundaries = null;
372         PointValuePair expected =
373             new PointValuePair(point(DIM,1.0),0.0);
374         doTest(new Rosen(), startPoint, insigma, boundaries,
375                 GoalType.MINIMIZE, LAMBDA, false, 1, 1e-13,
376                 1e-10, 1e-4, 1000000, expected);
377      }
378 
379     @Test
380     void testMath864() {
381         final CMAESOptimizer optimizer
382             = new CMAESOptimizer(30000, 0, true, 10,
383                                  0, new MersenneTwister(), false, null);
384         final MultivariateFunction fitnessFunction = new MultivariateFunction() {
385                 public double value(double[] parameters) {
386                     final double target = 1;
387                     final double error = target - parameters[0];
388                     return error * error;
389                 }
390             };
391 
392         final double[] start = { 0 };
393         final double[] lower = { -1e6 };
394         final double[] upper = { 1.5 };
395         final double[] sigma = { 1e-1 };
396         final double[] result = optimizer.optimize(new MaxEval(10000),
397                                                    new ObjectiveFunction(fitnessFunction),
398                                                    GoalType.MINIMIZE,
399                                                    new CMAESOptimizer.PopulationSize(5),
400                                                    new CMAESOptimizer.Sigma(sigma),
401                                                    new InitialGuess(start),
402                                                    new SimpleBounds(lower, upper)).getPoint();
403         assertTrue(result[0] <= upper[0],
404                           "Out of bounds (" + result[0] + " > " + upper[0] + ")");
405     }
406 
407     /**
408      * Cf. MATH-867
409      */
410     @Test
411     void testFitAccuracyDependsOnBoundary() {
412         final CMAESOptimizer optimizer
413             = new CMAESOptimizer(30000, 0, true, 10,
414                                  0, new MersenneTwister(), false, null);
415         final MultivariateFunction fitnessFunction = new MultivariateFunction() {
416                 public double value(double[] parameters) {
417                     final double target = 11.1;
418                     final double error = target - parameters[0];
419                     return error * error;
420                 }
421             };
422 
423         final double[] start = { 1 };
424 
425         // No bounds.
426         PointValuePair result = optimizer.optimize(new MaxEval(100000),
427                                                    new ObjectiveFunction(fitnessFunction),
428                                                    GoalType.MINIMIZE,
429                                                    SimpleBounds.unbounded(1),
430                                                    new CMAESOptimizer.PopulationSize(5),
431                                                    new CMAESOptimizer.Sigma(new double[] { 1e-1 }),
432                                                    new InitialGuess(start));
433         final double resNoBound = result.getPoint()[0];
434 
435         // Optimum is near the lower bound.
436         final double[] lower = { -20 };
437         final double[] upper = { 5e16 };
438         final double[] sigma = { 10 };
439         result = optimizer.optimize(new MaxEval(100000),
440                                     new ObjectiveFunction(fitnessFunction),
441                                     GoalType.MINIMIZE,
442                                     new CMAESOptimizer.PopulationSize(5),
443                                     new CMAESOptimizer.Sigma(sigma),
444                                     new InitialGuess(start),
445                                     new SimpleBounds(lower, upper));
446         final double resNearLo = result.getPoint()[0];
447 
448         // Optimum is near the upper bound.
449         lower[0] = -5e16;
450         upper[0] = 20;
451         result = optimizer.optimize(new MaxEval(100000),
452                                     new ObjectiveFunction(fitnessFunction),
453                                     GoalType.MINIMIZE,
454                                     new CMAESOptimizer.PopulationSize(5),
455                                     new CMAESOptimizer.Sigma(sigma),
456                                     new InitialGuess(start),
457                                     new SimpleBounds(lower, upper));
458         final double resNearHi = result.getPoint()[0];
459 
460         // System.out.println("resNoBound=" + resNoBound +
461         //                    " resNearLo=" + resNearLo +
462         //                    " resNearHi=" + resNearHi);
463 
464         // The two values currently differ by a substantial amount, indicating that
465         // the bounds definition can prevent reaching the optimum.
466         assertEquals(resNoBound, resNearLo, 1e-3);
467         assertEquals(resNoBound, resNearHi, 1e-3);
468     }
469 
470     /**
471      * @param func Function to optimize.
472      * @param startPoint Starting point.
473      * @param inSigma Individual input sigma.
474      * @param boundaries Upper / lower point limit.
475      * @param goal Minimization or maximization.
476      * @param lambda Population size used for offspring.
477      * @param isActive Covariance update mechanism.
478      * @param diagonalOnly Simplified covariance update.
479      * @param stopValue Termination criteria for optimization.
480      * @param fTol Tolerance relative error on the objective function.
481      * @param pointTol Tolerance for checking that the optimum is correct.
482      * @param maxEvaluations Maximum number of evaluations.
483      * @param expected Expected point / value.
484      */
485     private void doTest(MultivariateFunction func,
486                         double[] startPoint,
487                         double[] inSigma,
488                         double[][] boundaries,
489                         GoalType goal,
490                         int lambda,
491                         boolean isActive,
492                         int diagonalOnly,
493                         double stopValue,
494                         double fTol,
495                         double pointTol,
496                         int maxEvaluations,
497                         PointValuePair expected) {
498         int dim = startPoint.length;
499         // test diagonalOnly = 0 - slow but normally fewer feval#
500         CMAESOptimizer optim = new CMAESOptimizer(30000, stopValue, isActive, diagonalOnly,
501                                                   0, new MersenneTwister(), false, null);
502         PointValuePair result = boundaries == null ?
503             optim.optimize(new MaxEval(maxEvaluations),
504                            new ObjectiveFunction(func),
505                            goal,
506                            new InitialGuess(startPoint),
507                            SimpleBounds.unbounded(dim),
508                            new CMAESOptimizer.Sigma(inSigma),
509                            new CMAESOptimizer.PopulationSize(lambda)) :
510             optim.optimize(new MaxEval(maxEvaluations),
511                            new ObjectiveFunction(func),
512                            goal,
513                            new SimpleBounds(boundaries[0],
514                                             boundaries[1]),
515                            new InitialGuess(startPoint),
516                            new CMAESOptimizer.Sigma(inSigma),
517                            new CMAESOptimizer.PopulationSize(lambda));
518 
519         // System.out.println("sol=" + Arrays.toString(result.getPoint()));
520         assertEquals(expected.getValue(), result.getValue(), fTol);
521         for (int i = 0; i < dim; i++) {
522             assertEquals(expected.getPoint()[i], result.getPoint()[i], pointTol);
523         }
524 
525         assertTrue(optim.getIterations() > 0);
526     }
527 
528     private static double[] point(int n, double value) {
529         double[] ds = new double[n];
530         Arrays.fill(ds, value);
531         return ds;
532     }
533 
534     private static double[][] boundaries(int dim,
535             double lower, double upper) {
536         double[][] boundaries = new double[2][dim];
537         for (int i = 0; i < dim; i++)
538             boundaries[0][i] = lower;
539         for (int i = 0; i < dim; i++)
540             boundaries[1][i] = upper;
541         return boundaries;
542     }
543 
544     private static class Sphere implements MultivariateFunction {
545 
546         public double value(double[] x) {
547             double f = 0;
548             for (int i = 0; i < x.length; ++i)
549                 f += x[i] * x[i];
550             return f;
551         }
552     }
553 
554     private static class Cigar implements MultivariateFunction {
555         private double factor;
556 
557         Cigar() {
558             this(1e3);
559         }
560 
561         Cigar(double axisratio) {
562             factor = axisratio * axisratio;
563         }
564 
565         public double value(double[] x) {
566             double f = x[0] * x[0];
567             for (int i = 1; i < x.length; ++i)
568                 f += factor * x[i] * x[i];
569             return f;
570         }
571     }
572 
573     private static class Tablet implements MultivariateFunction {
574         private double factor;
575 
576         Tablet() {
577             this(1e3);
578         }
579 
580         Tablet(double axisratio) {
581             factor = axisratio * axisratio;
582         }
583 
584         public double value(double[] x) {
585             double f = factor * x[0] * x[0];
586             for (int i = 1; i < x.length; ++i)
587                 f += x[i] * x[i];
588             return f;
589         }
590     }
591 
592     private static class CigTab implements MultivariateFunction {
593         private double factor;
594 
595         CigTab() {
596             this(1e4);
597         }
598 
599         CigTab(double axisratio) {
600             factor = axisratio;
601         }
602 
603         public double value(double[] x) {
604             int end = x.length - 1;
605             double f = x[0] * x[0] / factor + factor * x[end] * x[end];
606             for (int i = 1; i < end; ++i)
607                 f += x[i] * x[i];
608             return f;
609         }
610     }
611 
612     private static class TwoAxes implements MultivariateFunction {
613 
614         private double factor;
615 
616         TwoAxes() {
617             this(1e6);
618         }
619 
620         TwoAxes(double axisratio) {
621             factor = axisratio * axisratio;
622         }
623 
624         public double value(double[] x) {
625             double f = 0;
626             for (int i = 0; i < x.length; ++i)
627                 f += (i < x.length / 2 ? factor : 1) * x[i] * x[i];
628             return f;
629         }
630     }
631 
632     private static class ElliRotated implements MultivariateFunction {
633         private Basis B = new Basis();
634         private double factor;
635 
636         ElliRotated() {
637             this(1e3);
638         }
639 
640         ElliRotated(double axisratio) {
641             factor = axisratio * axisratio;
642         }
643 
644         public double value(double[] x) {
645             double f = 0;
646             x = B.Rotate(x);
647             for (int i = 0; i < x.length; ++i)
648                 f += FastMath.pow(factor, i / (x.length - 1.)) * x[i] * x[i];
649             return f;
650         }
651     }
652 
653     private static class Elli implements MultivariateFunction {
654 
655         private double factor;
656 
657         Elli() {
658             this(1e3);
659         }
660 
661         Elli(double axisratio) {
662             factor = axisratio * axisratio;
663         }
664 
665         public double value(double[] x) {
666             double f = 0;
667             for (int i = 0; i < x.length; ++i)
668                 f += FastMath.pow(factor, i / (x.length - 1.)) * x[i] * x[i];
669             return f;
670         }
671     }
672 
673     private static class MinusElli implements MultivariateFunction {
674 
675         public double value(double[] x) {
676             return 1.0-(new Elli().value(x));
677         }
678     }
679 
680     private static class DiffPow implements MultivariateFunction {
681 
682         public double value(double[] x) {
683             double f = 0;
684             for (int i = 0; i < x.length; ++i)
685                 f += FastMath.pow(FastMath.abs(x[i]), 2. + 10 * (double) i
686                         / (x.length - 1.));
687             return f;
688         }
689     }
690 
691     private static class SsDiffPow implements MultivariateFunction {
692 
693         public double value(double[] x) {
694             double f = FastMath.pow(new DiffPow().value(x), 0.25);
695             return f;
696         }
697     }
698 
699     private static class Rosen implements MultivariateFunction {
700 
701         public double value(double[] x) {
702             double f = 0;
703             for (int i = 0; i < x.length - 1; ++i)
704                 f += 1e2 * (x[i] * x[i] - x[i + 1]) * (x[i] * x[i] - x[i + 1])
705                 + (x[i] - 1.) * (x[i] - 1.);
706             return f;
707         }
708     }
709 
710     private static class Ackley implements MultivariateFunction {
711         private double axisratio;
712 
713         Ackley(double axra) {
714             axisratio = axra;
715         }
716 
717         public Ackley() {
718             this(1);
719         }
720 
721         public double value(double[] x) {
722             double f = 0;
723             double res2 = 0;
724             double fac = 0;
725             for (int i = 0; i < x.length; ++i) {
726                 fac = FastMath.pow(axisratio, (i - 1.) / (x.length - 1.));
727                 f += fac * fac * x[i] * x[i];
728                 res2 += FastMath.cos(2. * FastMath.PI * fac * x[i]);
729             }
730             f = (20. - 20. * FastMath.exp(-0.2 * FastMath.sqrt(f / x.length))
731                     + FastMath.exp(1.) - FastMath.exp(res2 / x.length));
732             return f;
733         }
734     }
735 
736     private static class Rastrigin implements MultivariateFunction {
737 
738         private double axisratio;
739         private double amplitude;
740 
741         Rastrigin() {
742             this(1, 10);
743         }
744 
745         Rastrigin(double axisratio, double amplitude) {
746             this.axisratio = axisratio;
747             this.amplitude = amplitude;
748         }
749 
750         public double value(double[] x) {
751             double f = 0;
752             double fac;
753             for (int i = 0; i < x.length; ++i) {
754                 fac = FastMath.pow(axisratio, (i - 1.) / (x.length - 1.));
755                 if (i == 0 && x[i] < 0)
756                     fac *= 1.;
757                 f += fac * fac * x[i] * x[i] + amplitude
758                 * (1. - FastMath.cos(2. * FastMath.PI * fac * x[i]));
759             }
760             return f;
761         }
762     }
763 
764     private static class Basis {
765         double[][] basis;
766         Random rand = new Random(2); // use not always the same basis
767 
768         double[] Rotate(double[] x) {
769             GenBasis(x.length);
770             double[] y = new double[x.length];
771             for (int i = 0; i < x.length; ++i) {
772                 y[i] = 0;
773                 for (int j = 0; j < x.length; ++j)
774                     y[i] += basis[i][j] * x[j];
775             }
776             return y;
777         }
778 
779         void GenBasis(int DIM) {
780             if (basis != null ? basis.length == DIM : false)
781                 return;
782 
783             double sp;
784             int i, j, k;
785 
786             /* generate orthogonal basis */
787             basis = new double[DIM][DIM];
788             for (i = 0; i < DIM; ++i) {
789                 /* sample components gaussian */
790                 for (j = 0; j < DIM; ++j)
791                     basis[i][j] = rand.nextGaussian();
792                 /* substract projection of previous vectors */
793                 for (j = i - 1; j >= 0; --j) {
794                     for (sp = 0., k = 0; k < DIM; ++k)
795                         sp += basis[i][k] * basis[j][k]; /* scalar product */
796                     for (k = 0; k < DIM; ++k)
797                         basis[i][k] -= sp * basis[j][k]; /* substract */
798                 }
799                 /* normalize */
800                 for (sp = 0., k = 0; k < DIM; ++k)
801                     sp += basis[i][k] * basis[i][k]; /* squared norm */
802                 for (k = 0; k < DIM; ++k)
803                     basis[i][k] /= FastMath.sqrt(sp);
804             }
805         }
806     }
807 }