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