MATH-442

Implementation of the CMA-ES optimization algorithm provided by Dietmar Wolz
and Nikolaus Hansen.


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1071600 13f79535-47bb-0310-9956-ffa450edef68

src/site/xdoc/changes.xml

@@ -52,6 +52,9 @@ The <action> type attribute can be add,update,fix,remove.
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="erans" type="fix" issue="MATH-442" due-to="Dietmar Wolz">
+        Implementation of the CMA-ES optimization algorithm.
+      </action>
       <action dev="erans" type="fix" issue="MATH-513">
         The interface "ParametricRealFunction" (in package "optimization.fitting") has
         been renamed to "ParametricUnivariateRealFunction" and moved to package "analysis".

src/test/java/org/apache/commons/math/optimization/direct/CMAESOptimizerTest.java

@@ -0,0 +1,667 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.optimization.direct;
+
+import java.util.Arrays;
+import java.util.Random;
+
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.analysis.MultivariateRealFunction;
+import org.apache.commons.math.exception.MathUserException;
+import org.apache.commons.math.exception.MultiDimensionMismatchException;
+import org.apache.commons.math.exception.NoDataException;
+import org.apache.commons.math.exception.NotPositiveException;
+import org.apache.commons.math.exception.OutOfRangeException;
+import org.apache.commons.math.optimization.GoalType;
+import org.apache.commons.math.optimization.MultivariateRealOptimizer;
+import org.apache.commons.math.optimization.RealPointValuePair;
+import org.apache.commons.math.random.MersenneTwister;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test for {@link CMAESOptimizer}.
+ */
+public class CMAESOptimizerTest {
+
+    static final int DIM = 13;
+    static final int LAMBDA = 4 + (int)(3.*Math.log(DIM));
+   
+    @Test(expected = OutOfRangeException.class)
+    public void testInitOutofbounds() throws MathUserException, MathException {
+        double[] startPoint = point(DIM,3);
+        double[] insigma = null;
+        double[][] boundaries = boundaries(DIM,-1,2);
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+    
+    @Test(expected = MultiDimensionMismatchException.class)
+    public void testBoundariesDimensionMismatch() throws MathUserException, MathException {
+        double[] startPoint = point(DIM,0.5);
+        double[] insigma = null;
+        double[][] boundaries = boundaries(DIM+1,-1,2);
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test(expected = NoDataException.class)
+    public void testBoundariesNoData() throws MathUserException, MathException {
+        double[] startPoint = point(DIM,0.5);
+        double[] insigma = null;
+        double[][] boundaries = boundaries(DIM,-1,2);
+        boundaries[1] = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test(expected = NotPositiveException.class)
+    public void testInputSigmaNegative() throws MathUserException, MathException {
+        double[] startPoint = point(DIM,0.5);
+        double[] insigma = point(DIM,-0.5);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test(expected = OutOfRangeException.class)
+    public void testInputSigmaOutOfRange() throws MathUserException, MathException {
+        double[] startPoint = point(DIM,0.5);
+        double[] insigma = point(DIM, 1.1);
+        double[][] boundaries = boundaries(DIM,-1,2);
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test(expected = MultiDimensionMismatchException.class)
+    public void testInputSigmaDimensionMismatch() throws MathUserException, MathException {
+        double[] startPoint = point(DIM,0.5);
+        double[] insigma = point(DIM+1,-0.5);
+        double[][] boundaries = null;;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+    
+    @Test
+    public void testRosen() throws MathException {
+        double[] startPoint = point(DIM,0.1);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test
+    public void testMaximize() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),1.0);
+        doTest(new MinusElli(), startPoint, insigma, boundaries,
+                GoalType.MAXIMIZE, LAMBDA, true, 0, 1.0-1e-13,
+                2e-10, 5e-6, 100000, expected);
+        doTest(new MinusElli(), startPoint, insigma, boundaries,
+                GoalType.MAXIMIZE, LAMBDA, false, 0, 1.0-1e-13,
+                2e-10, 5e-6, 100000, expected);
+        boundaries = boundaries(DIM,-0.3,0.3); 
+        startPoint = point(DIM,0.1);
+        doTest(new MinusElli(), startPoint, insigma, boundaries,
+                GoalType.MAXIMIZE, LAMBDA, true, 0, 1.0-1e-13,
+                2e-10, 5e-6, 100000, expected);
+    }
+
+    @Test
+    public void testEllipse() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new Elli(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+        doTest(new Elli(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test
+    public void testElliRotated() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new ElliRotated(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+        doTest(new ElliRotated(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test
+    public void testCigar() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new Cigar(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 200000, expected);
+        doTest(new Cigar(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test
+    public void testTwoAxes() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new TwoAxes(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 200000, expected);
+        doTest(new TwoAxes(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
+                1e-8, 1e-3, 200000, expected);
+    }
+
+    @Test
+    public void testCigTab() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.3);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new CigTab(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 5e-5, 100000, expected);
+        doTest(new CigTab(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+                1e-13, 5e-5, 100000, expected);
+    }
+
+    @Test
+    public void testSphere() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new Sphere(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+        doTest(new Sphere(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test
+    public void testTablet() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new Tablet(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+        doTest(new Tablet(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, false, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test
+    public void testDiffPow() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new DiffPow(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 10, true, 0, 1e-13,
+                1e-8, 1e-1, 100000, expected);
+        doTest(new DiffPow(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 10, false, 0, 1e-13,
+                1e-8, 2e-1, 100000, expected);
+    }
+
+    @Test
+    public void testSsDiffPow() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new SsDiffPow(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 10, true, 0, 1e-13,
+                1e-4, 1e-1, 200000, expected);
+        doTest(new SsDiffPow(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 10, false, 0, 1e-13,
+                1e-4, 1e-1, 200000, expected);
+    }
+
+    @Test
+    public void testAckley() throws MathException {
+        double[] startPoint = point(DIM,1.0);
+        double[] insigma = point(DIM,1.0);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new Ackley(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
+                1e-9, 1e-5, 100000, expected);
+        doTest(new Ackley(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
+                1e-9, 1e-5, 100000, expected);
+    }
+
+    @Test
+    public void testRastrigin() throws MathException {
+        double[] startPoint = point(DIM,0.1);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,0.0),0.0);
+        doTest(new Rastrigin(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, (int)(200*Math.sqrt(DIM)), true, 0, 1e-13,
+                1e-13, 1e-6, 200000, expected);
+        doTest(new Rastrigin(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, (int)(200*Math.sqrt(DIM)), false, 0, 1e-13,
+                1e-13, 1e-6, 200000, expected);
+    }
+
+    @Test
+    public void testConstrainedRosen() throws MathException {
+        double[] startPoint = point(DIM,0.1);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = boundaries(DIM,-1,2);
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 2*LAMBDA, true, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, 2*LAMBDA, false, 0, 1e-13,
+                1e-13, 1e-6, 100000, expected);
+    }
+
+    @Test
+    public void testDiagonalRosen() throws MathException {
+        double[] startPoint = point(DIM,0.1);
+        double[] insigma = point(DIM,0.1);
+        double[][] boundaries = null;
+        RealPointValuePair expected =
+            new RealPointValuePair(point(DIM,1.0),0.0);
+        doTest(new Rosen(), startPoint, insigma, boundaries,
+                GoalType.MINIMIZE, LAMBDA, false, 1, 1e-13,
+                1e-10, 1e-4, 1000000, expected);
+     }
+
+    /**
+     * @param func Function to optimize.
+     * @param startPoint Starting point.
+     * @param inSigma Individual input sigma.
+     * @param boundaries Upper / lower point limit.
+     * @param goal Minimization or maximization.
+     * @param lambda Population size used for offspring.
+     * @param isActive Covariance update mechanism.
+     * @param diagonalOnly Simplified covariance update.
+     * @param stopValue Termination criteria for optimization.
+     * @param fTol Tolerance relative error on the objective function.
+     * @param pointTol Tolerance for checking that the optimum is correct.
+     * @param maxEvaluations Maximum number of evaluations.
+     * @param expected Expected point / value.
+     */
+    private void doTest(MultivariateRealFunction func,
+            double[] startPoint,
+            double[] inSigma,
+            double[][] boundaries,
+            GoalType goal,
+            int lambda,
+            boolean isActive,
+            int diagonalOnly, 
+            double stopValue,
+            double fTol,
+            double pointTol,
+            int maxEvaluations,
+            RealPointValuePair expected)
+    throws MathException {
+        int dim = startPoint.length;
+        // test diagonalOnly = 0 - slow but normally fewer feval#
+        MultivariateRealOptimizer optim =
+            new CMAESOptimizer(
+                    lambda, inSigma, boundaries, 30000,
+                    stopValue, isActive, diagonalOnly, 0, new MersenneTwister(),false);
+        RealPointValuePair result = optim.optimize(maxEvaluations, func, goal, startPoint);
+        Assert.assertEquals(expected.getValue(),
+                result.getValue(), fTol);
+        for (int i = 0; i < dim; i++) {
+            Assert.assertEquals(expected.getPoint()[i],
+                    result.getPoint()[i], pointTol);
+        }
+    }
+
+    private static double[] point(int n, double value) {
+        double[] ds = new double[n];
+        Arrays.fill(ds, value);
+        return ds;
+    }
+
+    private static double[][] boundaries(int dim,
+            double lower, double upper) {
+        double[][] boundaries = new double[2][dim];
+        for (int i = 0; i < dim; i++)
+            boundaries[0][i] = lower;
+        for (int i = 0; i < dim; i++)
+            boundaries[1][i] = upper;
+        return boundaries;
+    }
+
+    private static class Sphere implements MultivariateRealFunction {
+
+        public double value(double[] x) {
+            double f = 0;
+            for (int i = 0; i < x.length; ++i)
+                f += x[i] * x[i];
+            return f;
+        }
+    }
+
+    private static class Cigar implements MultivariateRealFunction {
+        private double factor;
+
+        Cigar() {
+            this(1e3);
+        }
+
+        Cigar(double axisratio) {
+            factor = axisratio * axisratio;
+        }
+
+        public double value(double[] x) {
+            double f = x[0] * x[0];
+            for (int i = 1; i < x.length; ++i)
+                f += factor * x[i] * x[i];
+            return f;
+        }
+    }
+
+    private static class Tablet implements MultivariateRealFunction {
+        private double factor;
+
+        Tablet() {
+            this(1e3);
+        }
+
+        Tablet(double axisratio) {
+            factor = axisratio * axisratio;
+        }
+
+        public double value(double[] x) {
+            double f = factor * x[0] * x[0];
+            for (int i = 1; i < x.length; ++i)
+                f += x[i] * x[i];
+            return f;
+        }
+    }
+
+    private static class CigTab implements MultivariateRealFunction {
+        private double factor;
+
+        CigTab() {
+            this(1e4);
+        }
+
+        CigTab(double axisratio) {
+            factor = axisratio;
+        }
+
+        public double value(double[] x) {
+            int end = x.length - 1;
+            double f = x[0] * x[0] / factor + factor * x[end] * x[end];
+            for (int i = 1; i < end; ++i)
+                f += x[i] * x[i];
+            return f;
+        }
+    }
+
+    private static class TwoAxes implements MultivariateRealFunction {
+
+        private double factor;
+
+        TwoAxes() {
+            this(1e6);
+        }
+
+        TwoAxes(double axisratio) {
+            factor = axisratio * axisratio;
+        }
+
+        public double value(double[] x) {
+            double f = 0;
+            for (int i = 0; i < x.length; ++i)
+                f += (i < x.length / 2 ? factor : 1) * x[i] * x[i];
+            return f;
+        }
+    }
+
+    private static class ElliRotated implements MultivariateRealFunction {
+        private Basis B = new Basis();
+        private double factor;
+
+        ElliRotated() {
+            this(1e3);
+        }
+
+        ElliRotated(double axisratio) {
+            factor = axisratio * axisratio;
+        }
+
+        public double value(double[] x) {
+            double f = 0;
+            x = B.Rotate(x);
+            for (int i = 0; i < x.length; ++i)
+                f += Math.pow(factor, i / (x.length - 1.)) * x[i] * x[i];
+            return f;
+        }
+    }
+
+    private static class Elli implements MultivariateRealFunction {
+
+        private double factor;
+
+        Elli() {
+            this(1e3);
+        }
+
+        Elli(double axisratio) {
+            factor = axisratio * axisratio;
+        }
+
+        public double value(double[] x) {
+            double f = 0;
+            for (int i = 0; i < x.length; ++i)
+                f += Math.pow(factor, i / (x.length - 1.)) * x[i] * x[i];
+            return f;
+        }
+    }
+
+    private static class MinusElli implements MultivariateRealFunction {
+
+        public double value(double[] x) {
+            return 1.0-(new Elli().value(x));
+        }
+    }
+
+    private static class DiffPow implements MultivariateRealFunction {
+
+        public double value(double[] x) {
+            double f = 0;
+            for (int i = 0; i < x.length; ++i)
+                f += Math.pow(Math.abs(x[i]), 2. + 10 * (double) i
+                        / (x.length - 1.));
+            return f;
+        }
+    }
+
+    private static class SsDiffPow implements MultivariateRealFunction {
+
+        public double value(double[] x) {
+            double f = Math.pow(new DiffPow().value(x), 0.25);
+            return f;
+        }
+    }
+
+    private static class Rosen implements MultivariateRealFunction {
+
+        public double value(double[] x) {
+            double f = 0;
+            for (int i = 0; i < x.length - 1; ++i)
+                f += 1e2 * (x[i] * x[i] - x[i + 1]) * (x[i] * x[i] - x[i + 1])
+                + (x[i] - 1.) * (x[i] - 1.);
+            return f;
+        }
+    }
+
+    private static class Ackley implements MultivariateRealFunction {
+        private double axisratio;
+
+        Ackley(double axra) {
+            axisratio = axra;
+        }
+
+        public Ackley() {
+            this(1);
+        }
+
+        public double value(double[] x) {
+            double f = 0;
+            double res2 = 0;
+            double fac = 0;
+            for (int i = 0; i < x.length; ++i) {
+                fac = Math.pow(axisratio, (i - 1.) / (x.length - 1.));
+                f += fac * fac * x[i] * x[i];
+                res2 += Math.cos(2. * Math.PI * fac * x[i]);
+            }
+            f = (20. - 20. * Math.exp(-0.2 * Math.sqrt(f / x.length))
+                    + Math.exp(1.) - Math.exp(res2 / x.length));
+            return f;
+        }
+    }
+
+    private static class Rastrigin implements MultivariateRealFunction {
+
+        private double axisratio;
+        private double amplitude;
+
+        Rastrigin() {
+            this(1, 10);
+        }
+
+        Rastrigin(double axisratio, double amplitude) {
+            this.axisratio = axisratio;
+            this.amplitude = amplitude;
+        }
+
+        public double value(double[] x) {
+            double f = 0;
+            double fac;
+            for (int i = 0; i < x.length; ++i) {
+                fac = Math.pow(axisratio, (i - 1.) / (x.length - 1.));
+                if (i == 0 && x[i] < 0)
+                    fac *= 1.;
+                f += fac * fac * x[i] * x[i] + amplitude
+                * (1. - Math.cos(2. * Math.PI * fac * x[i]));
+            }
+            return f;
+        }
+    }
+
+    private static class Basis {
+        double[][] basis;
+        Random rand = new Random(2); // use not always the same basis
+
+        double[] Rotate(double[] x) {
+            GenBasis(x.length);
+            double[] y = new double[x.length];
+            for (int i = 0; i < x.length; ++i) {
+                y[i] = 0;
+                for (int j = 0; j < x.length; ++j)
+                    y[i] += basis[i][j] * x[j];
+            }
+            return y;
+        }
+
+        void GenBasis(int DIM) {
+            if (basis != null ? basis.length == DIM : false)
+                return;
+
+            double sp;
+            int i, j, k;
+
+            /* generate orthogonal basis */
+            basis = new double[DIM][DIM];
+            for (i = 0; i < DIM; ++i) {
+                /* sample components gaussian */
+                for (j = 0; j < DIM; ++j)
+                    basis[i][j] = rand.nextGaussian();
+                /* substract projection of previous vectors */
+                for (j = i - 1; j >= 0; --j) {
+                    for (sp = 0., k = 0; k < DIM; ++k)
+                        sp += basis[i][k] * basis[j][k]; /* scalar product */
+                    for (k = 0; k < DIM; ++k)
+                        basis[i][k] -= sp * basis[j][k]; /* substract */
+                }
+                /* normalize */
+                for (sp = 0., k = 0; k < DIM; ++k)
+                    sp += basis[i][k] * basis[i][k]; /* squared norm */
+                for (k = 0; k < DIM; ++k)
+                    basis[i][k] /= Math.sqrt(sp);
+            }
+        }
+    }
+}