From bb201d926ebf225439d0007e7015365a9fbc3578 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Thu, 10 Nov 2011 20:08:55 +0000 Subject: [PATCH] Added adapters for simple bounds constraints optimization. The adapters are useful only for optimizers that do not support simple bounds constraints by themselves (i.e. Nelder-Mead and Torczon's multidirectional). Two adapters are available, one performs a mapping between the whole real range and the bounded range (bounds being set component wise), and one uses a penalty function. JIRA: MATH-196 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1200516 13f79535-47bb-0310-9956-ffa450edef68 --- ...ultivariateRealFunctionMappingAdapter.java | 300 ++++++++++++++++++ ...ultivariateRealFunctionPenaltyAdapter.java | 188 +++++++++++ .../optimization/direct/SimplexOptimizer.java | 12 + src/site/xdoc/changes.xml | 5 + src/site/xdoc/userguide/optimization.xml | 44 ++- ...variateRealFunctionMappingAdapterTest.java | 190 +++++++++++ ...variateRealFunctionPenaltyAdapterTest.java | 192 +++++++++++ 7 files changed, 926 insertions(+), 5 deletions(-) create mode 100644 src/main/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapter.java create mode 100644 src/main/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapter.java create mode 100644 src/test/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapterTest.java create mode 100644 src/test/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapterTest.java diff --git a/src/main/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapter.java b/src/main/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapter.java new file mode 100644 index 000000000..911da318f --- /dev/null +++ b/src/main/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapter.java @@ -0,0 +1,300 @@ +/* + * 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 org.apache.commons.math.analysis.MultivariateRealFunction; +import org.apache.commons.math.analysis.UnivariateRealFunction; +import org.apache.commons.math.analysis.function.Logit; +import org.apache.commons.math.analysis.function.Sigmoid; +import org.apache.commons.math.exception.DimensionMismatchException; +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.NumberIsTooSmallException; +import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.util.MathUtils; + +/** + *

Adapter for mapping bounded {@link MultivariateRealFunction} to unbounded ones.

+ * + *

+ * This adapter can be used to wrap functions subject to simple bounds on + * parameters so they can be used by optimizers that do not directly + * support simple bounds. + *

+ *

+ * The principle is that the user function that will be wrapped will see its + * parameters bounded as required, i.e when its {@code value} method is called + * with argument array {@code point}, the elements array will fulfill requirement + * {@code lower[i] <= point[i] <= upper[i]} for all i. Some of the components + * may be unbounded or bounded only on one side if the corresponding bound is + * set to an infinite value. The optimizer will not manage the user function by + * itself, but it will handle this adapter and it is this adapter that will take + * care the bounds are fulfilled. The adapter {@link #value(double[])} method will + * be called by the optimizer with unbound parameters, and the adapter will map + * the unbounded value to the bounded range using appropriate functions like + * {@link Sigmoid} for double bounded elements for example. + *

+ *

+ * As the optimizer sees only unbounded parameters, it should be noted that the + * start point or simplex expected by the optimizer should be unbounded, so the + * user is responsible for converting his bounded point to unbounded by calling + * {@link #boundedToUnbounded(double[])} before providing them to the optimizer. + * For the same reason, the point returned by the {@link + * org.apache.commons.math.optimization.BaseMultivariateRealOptimizer#optimize(int, + * MultivariateRealFunction, org.apache.commons.math.optimization.GoalType, double[])} + * method is unbounded. So to convert this point to bounded, users must call + * {@link #unboundedToBounded(double[])} by themselves!

+ *

+ * This adapter is only a poor man solution to simple bounds optimization constraints + * that can be used with simple optimizers like {@link SimplexOptimizer} with {@link + * NelderMeadSimplex} or {@link MultiDirectionalSimplex}. A better solution is to use + * an optimizer that directly supports simple bounds like {@link CMAESOptimizer} or + * {@link BOBYQAOptimizer}. One caveat of this poor man solution is that behavior near + * the bounds may be numerically unstable as bounds are mapped from infinite values. + * Another caveat is that convergence values are evaluated by the optimizer with respect + * to unbounded variables, so there will be scales differences when converted to bounded + * variables. + *

+ * + * @see MultivariateRealFunctionPenaltyAdapter + * + * @version $Id$ + * @since 3.0 + */ + +public class MultivariateRealFunctionMappingAdapter implements MultivariateRealFunction { + + /** Underlying bounded function. */ + private final MultivariateRealFunction bounded; + + /** Mapping functions. */ + private final Mapper[] mappers; + + /** Simple constructor. + * @param bounded bounded function + * @param lower lower bounds for each element of the input parameters array + * (some elements may be set to {@code Double.NEGATIVE_INFINITY} for + * unbounded values) + * @param upper upper bounds for each element of the input parameters array + * (some elements may be set to {@code Double.POSITIVE_INFINITY} for + * unbounded values) + * @exception MathIllegalArgumentException if lower and upper bounds are not + * consistent, either according to dimension or to values + */ + public MultivariateRealFunctionMappingAdapter(final MultivariateRealFunction bounded, + final double[] lower, final double[] upper) { + + // safety checks + MathUtils.checkNotNull(lower); + MathUtils.checkNotNull(upper); + if (lower.length != upper.length) { + throw new DimensionMismatchException(lower.length, upper.length); + } + for (int i = 0; i < lower.length; ++i) { + // note the following test is written in such a way it also fails for NaN + if (!(upper[i] >= lower[i])) { + throw new NumberIsTooSmallException(upper[i], lower[i], true); + } + } + + this.bounded = bounded; + this.mappers = new Mapper[lower.length]; + for (int i = 0; i < mappers.length; ++i) { + if (Double.isInfinite(lower[i])) { + if (Double.isInfinite(upper[i])) { + // element is unbounded, no transformation is needed + mappers[i] = new NoBoundsMapper(); + } else { + // element is simple-bounded on the upper side + mappers[i] = new UpperBoundMapper(upper[i]); + } + } else { + if (Double.isInfinite(upper[i])) { + // element is simple-bounded on the lower side + mappers[i] = new LowerBoundMapper(lower[i]); + } else { + // element is double-bounded + mappers[i] = new LowerUpperBoundMapper(lower[i], upper[i]); + } + } + } + + } + + /** Map an array from unbounded to bounded. + * @param x unbounded value + * @return bounded value + */ + public double[] unboundedToBounded(double[] point) { + + // map unbounded input point to bounded point + final double[] mapped = new double[mappers.length]; + for (int i = 0; i < mappers.length; ++i) { + mapped[i] = mappers[i].unboundedToBounded(point[i]); + } + + return mapped; + + } + + /** Map an array from bounded to unbounded. + * @param y bounded value + * @return unbounded value + */ + public double[] boundedToUnbounded(double[] point) { + + // map bounded input point to unbounded point + final double[] mapped = new double[mappers.length]; + for (int i = 0; i < mappers.length; ++i) { + mapped[i] = mappers[i].boundedToUnbounded(point[i]); + } + + // call underlying function + return mapped; + + } + + /** Compute the underlying function value from an unbounded point. + *

+ * This method simply bounds the unbounded point using the mappings + * set up at construction and calls the underlying function using + * the bounded point. + *

+ * @see #unboundedToBounded(double[]) + */ + public double value(double[] point) { + return bounded.value(unboundedToBounded(point)); + } + + /** Mapping interface. */ + private static interface Mapper { + + /** Map a value from unbounded to bounded. + * @param y unbounded value + * @return bounded value + */ + public double unboundedToBounded(double y); + + /** Map a value from bounded to unbounded. + * @param x bounded value + * @return unbounded value + */ + public double boundedToUnbounded(double x); + + } + + /** Local class for no bounds mapping. */ + private static class NoBoundsMapper implements Mapper { + + /** Simple constructor. + */ + public NoBoundsMapper() { + } + + /** {@inheritDoc} */ + public double unboundedToBounded(final double y) { + return y; + } + + /** {@inheritDoc} */ + public double boundedToUnbounded(final double x) { + return x; + } + + } + + /** Local class for lower bounds mapping. */ + private static class LowerBoundMapper implements Mapper { + + /** Low bound. */ + private final double lower; + + /** Simple constructor. + * @param lower lower bound + */ + public LowerBoundMapper(final double lower) { + this.lower = lower; + } + + /** {@inheritDoc} */ + public double unboundedToBounded(final double y) { + return lower + FastMath.exp(y); + } + + /** {@inheritDoc} */ + public double boundedToUnbounded(final double x) { + return FastMath.log(x - lower); + } + + } + + /** Local class for upper bounds mapping. */ + private static class UpperBoundMapper implements Mapper { + + /** Upper bound. */ + private final double upper; + + /** Simple constructor. + * @param upper upper bound + */ + public UpperBoundMapper(final double upper) { + this.upper = upper; + } + + /** {@inheritDoc} */ + public double unboundedToBounded(final double y) { + return upper - FastMath.exp(-y); + } + + /** {@inheritDoc} */ + public double boundedToUnbounded(final double x) { + return -FastMath.log(upper - x); + } + + } + + /** Local class for lower and bounds mapping. */ + private static class LowerUpperBoundMapper implements Mapper { + + /** Function from unbounded to bounded. */ + private final UnivariateRealFunction boundingFunction; + + /** Function from bounded to unbounded. */ + private final UnivariateRealFunction unboundingFunction; + + /** Simple constructor. + * @param lower lower bound + * @param upper upper bound + */ + public LowerUpperBoundMapper(final double lower, final double upper) { + boundingFunction = new Sigmoid(lower, upper); + unboundingFunction = new Logit(lower, upper); + } + + /** {@inheritDoc} */ + public double unboundedToBounded(final double y) { + return boundingFunction.value(y); + } + + /** {@inheritDoc} */ + public double boundedToUnbounded(final double x) { + return unboundingFunction.value(x); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapter.java b/src/main/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapter.java new file mode 100644 index 000000000..56429b494 --- /dev/null +++ b/src/main/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapter.java @@ -0,0 +1,188 @@ +/* + * 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 org.apache.commons.math.analysis.MultivariateRealFunction; +import org.apache.commons.math.exception.DimensionMismatchException; +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.NumberIsTooSmallException; +import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.util.MathUtils; + +/** + *

Adapter extending bounded {@link MultivariateRealFunction} to an unbouded + * domain using a penalty function.

+ * + *

+ * This adapter can be used to wrap functions subject to simple bounds on + * parameters so they can be used by optimizers that do not directly + * support simple bounds. + *

+ *

+ * The principle is that the user function that will be wrapped will see its + * parameters bounded as required, i.e when its {@code value} method is called + * with argument array {@code point}, the elements array will fulfill requirement + * {@code lower[i] <= point[i] <= upper[i]} for all i. Some of the components + * may be unbounded or bounded only on one side if the corresponding bound is + * set to an infinite value. The optimizer will not manage the user function by + * itself, but it will handle this adapter and it is this adapter that will take + * care the bounds are fulfilled. The adapter {@link #value(double[])} method will + * be called by the optimizer with unbound parameters, and the adapter will check + * if the parameters is within range or not. If it is in range, then the underlying + * user function will be called, and if it is not the value of a penalty function + * will be returned instead. + *

+ *

+ * This adapter is only a poor man solution to simple bounds optimization constraints + * that can be used with simple optimizers like {@link SimplexOptimizer} with {@link + * NelderMeadSimplex} or {@link MultiDirectionalSimplex}. A better solution is to use + * an optimizer that directly supports simple bounds like {@link CMAESOptimizer} or + * {@link BOBYQAOptimizer}. One caveat of this poor man solution is that if start point + * or start simplex is completely outside of the allowed range, only the penalty function + * is used, and the optimizer may converge without ever entering the range. + *

+ * + * @see MultivariateRealFunctionMappingAdapter + * + * @version $Id$ + * @since 3.0 + */ + +public class MultivariateRealFunctionPenaltyAdapter implements MultivariateRealFunction { + + /** Underlying bounded function. */ + private final MultivariateRealFunction bounded; + + /** Lower bounds. */ + private final double[] lower; + + /** Upper bounds. */ + private final double[] upper; + + /** Penalty offset. */ + private final double offset; + + /** Penalty scales. */ + private final double[] scale; + + /** Simple constructor. + *

+ * When the optimizer provided points are out of range, the value of the + * penalty function will be used instead of the value of the underlying + * function. In order for this penalty to be effective in rejecting this + * point during the optimization process, the penalty function value should + * be defined with care. This value is computed as: + * 

+     *   penalty(point) = offset + ∑i[scale[i] * √|point[i]-boundary[i]|]
+     *
+ * where indices i correspond to all the components that violates their boundaries. + *

+ *

+ * So when attempting a function minimization, offset should be larger than + * the maximum expected value of the underlying function and scale components + * should all be positive. When attempting a function maximization, offset + * should be lesser than the minimum expected value of the underlying function + * and scale components should all be negative. + * minimization, and lesser than the minimum expected value of the underlying + * function when attempting maximization. + *

+ *

+ * These choices for the penalty function have two properties. First, all out + * of range points will return a function value that is worse than the value + * returned by any in range point. Second, the penalty is worse for large + * boundaries violation than for small violations, so the optimizer has an hint + * about the direction in which it should search for acceptable points. + *

+ * @param bounded bounded function + * @param lower lower bounds for each element of the input parameters array + * (some elements may be set to {@code Double.NEGATIVE_INFINITY} for + * unbounded values) + * @param upper upper bounds for each element of the input parameters array + * (some elements may be set to {@code Double.POSITIVE_INFINITY} for + * unbounded values) + * @param offset base offset of the penalty function + * @param scale scale of the penalty function + * @exception MathIllegalArgumentException if lower bounds, upper bounds and + * scales are not consistent, either according to dimension or to bounadary + * values + */ + public MultivariateRealFunctionPenaltyAdapter(final MultivariateRealFunction bounded, + final double[] lower, final double[] upper, + final double offset, final double[] scale) { + + // safety checks + MathUtils.checkNotNull(lower); + MathUtils.checkNotNull(upper); + MathUtils.checkNotNull(scale); + if (lower.length != upper.length) { + throw new DimensionMismatchException(lower.length, upper.length); + } + if (lower.length != scale.length) { + throw new DimensionMismatchException(lower.length, scale.length); + } + for (int i = 0; i < lower.length; ++i) { + // note the following test is written in such a way it also fails for NaN + if (!(upper[i] >= lower[i])) { + throw new NumberIsTooSmallException(upper[i], lower[i], true); + } + } + + this.bounded = bounded; + this.lower = lower.clone(); + this.upper = upper.clone(); + this.offset = offset; + this.scale = scale.clone(); + + } + + /** Compute the underlying function value from an unbounded point. + *

+ * This method simply bounds the unbounded point using the mappings + * set up at construction and calls the underlying function using + * the bounded point. + *

+ * @see #unboundedToBounded(double[]) + */ + public double value(double[] point) { + + for (int i = 0; i < scale.length; ++i) { + if ((point[i] < lower[i]) || (point[i] > upper[i])) { + // bound violation starting at this component + double sum = 0; + for (int j = i; j < scale.length; ++j) { + final double overshoot; + if (point[j] < lower[j]) { + overshoot = lower[j] - point[j]; + } else if (point[j] > upper[j]) { + overshoot = point[j] - upper[j]; + } else { + overshoot = 0; + } + sum += FastMath.sqrt(overshoot); + } + return offset + sum; + } + } + + // all boundaries are fulfilled, we are in the expected + // domain of the underlying function + return bounded.value(point); + + } + +} diff --git a/src/main/java/org/apache/commons/math/optimization/direct/SimplexOptimizer.java b/src/main/java/org/apache/commons/math/optimization/direct/SimplexOptimizer.java index 1143b1978..d8c37eea7 100644 --- a/src/main/java/org/apache/commons/math/optimization/direct/SimplexOptimizer.java +++ b/src/main/java/org/apache/commons/math/optimization/direct/SimplexOptimizer.java @@ -66,8 +66,20 @@ import org.apache.commons.math.optimization.MultivariateRealOptimizer; * previous and current simplex to the convergence checker, not the best * ones. *

+ *

+ * This simplex optimizer implementation does not directly support constrained + * optimization with simple bounds, so for such optimizations, either a more + * dedicated method must be used like {@link CMAESOptimizer} or {@link + * BOBYQAOptimizer}, or the optimized method must be wrapped in an adapter like + * {@link MultivariateRealFunctionMappingAdapter} or {@link + * MultivariateRealFunctionPenaltyAdapter}. + *

* * @see AbstractSimplex + * @see MultivariateRealFunctionMappingAdapter + * @see MultivariateRealFunctionPenaltyAdapter + * @see CMAESOptimizer + * @see BOBYQAOptimizer * @version $Id$ * @since 3.0 */ diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 27e4e48e6..61fb88757 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,11 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Added adapters for simple bound constraints optimization that can be + used for all direct optimization methods, including the ones that do not + support constraints by themselves. + CMA-ES optimizer input sigma is now consistent with boundaries range units. diff --git a/src/site/xdoc/userguide/optimization.xml b/src/site/xdoc/userguide/optimization.xml index b0f3ef390..83189db20 100644 --- a/src/site/xdoc/userguide/optimization.xml +++ b/src/site/xdoc/userguide/optimization.xml @@ -154,11 +154,45 @@ already provided by the minimizes method).

- The direct package provides two solvers. The first one is the classical - - Nelder-Mead method. The second one is Virginia Torczon's - - multi-directional method. + The direct package provides four solvers: +

+

+

+ The first two simplex-based methods do not handle simple bounds constraints by themselves. + However there are two adapters ( + MultivariateRealFunctionMappingAdapter and + MultivariateRealFunctionPenaltyAdapter) that can be used to wrap the user function in + such a way the wrapped function is unbounded and can be used with these optimizers, despite + the fact the underlying function is still bounded and will be called only with feasible + points that fulfill the constraints. Note however that using these adapters is only a + apter is only a poor man solution to simple bounds optimization constraints. A better solution + is to use an optimizer that directly supports simple bounds. Two caveats of the mapping adapter + solution is that behavior near the bounds may be numerically unstable as bounds are mapped from + infinite values, and that convergence values are evaluated by the optimizer with respect + to unbounded variables, so there will be scales differences when converted to bounded variables. + One caveat of penalty adapter is that if start point or start simplex is outside of the allowed + range, only the penalty function is used, and the optimizer may converge without ever entering + the alowed range. +

+

+ The last methods do handle simple bounds constraints directly, so the adapters are not needed + with them.

diff --git a/src/test/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapterTest.java b/src/test/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapterTest.java new file mode 100644 index 000000000..1b008978a --- /dev/null +++ b/src/test/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapterTest.java @@ -0,0 +1,190 @@ +/* + * 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 org.apache.commons.math.analysis.MultivariateRealFunction; +import org.apache.commons.math.optimization.GoalType; +import org.apache.commons.math.optimization.RealPointValuePair; +import org.junit.Assert; +import org.junit.Test; + +public class MultivariateRealFunctionMappingAdapterTest { + + @Test + public void testStartSimplexInsideRange() { + + final BiQuadratic biQuadratic = new BiQuadratic(2.0, 2.5, 1.0, 3.0, 2.0, 3.0); + final MultivariateRealFunctionMappingAdapter wrapped = + new MultivariateRealFunctionMappingAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper()); + + SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30); + optimizer.setSimplex(new NelderMeadSimplex(new double[][] { + wrapped.boundedToUnbounded(new double[] { 1.5, 2.75 }), + wrapped.boundedToUnbounded(new double[] { 1.5, 2.95 }), + wrapped.boundedToUnbounded(new double[] { 1.7, 2.90 }) + })); + + final RealPointValuePair optimum + = optimizer.optimize(300, wrapped, GoalType.MINIMIZE, + wrapped.boundedToUnbounded(new double[] { 1.5, 2.25 })); + final double[] bounded = wrapped.unboundedToBounded(optimum.getPoint()); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), bounded[0], 2e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), bounded[1], 2e-7); + + } + + @Test + public void testOptimumOutsideRange() { + + final BiQuadratic biQuadratic = new BiQuadratic(4.0, 0.0, 1.0, 3.0, 2.0, 3.0); + final MultivariateRealFunctionMappingAdapter wrapped = + new MultivariateRealFunctionMappingAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper()); + + SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30); + optimizer.setSimplex(new NelderMeadSimplex(new double[][] { + wrapped.boundedToUnbounded(new double[] { 1.5, 2.75 }), + wrapped.boundedToUnbounded(new double[] { 1.5, 2.95 }), + wrapped.boundedToUnbounded(new double[] { 1.7, 2.90 }) + })); + + final RealPointValuePair optimum + = optimizer.optimize(100, wrapped, GoalType.MINIMIZE, + wrapped.boundedToUnbounded(new double[] { 1.5, 2.25 })); + final double[] bounded = wrapped.unboundedToBounded(optimum.getPoint()); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), bounded[0], 2e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), bounded[1], 2e-7); + + } + + @Test + public void testUnbounded() { + + final BiQuadratic biQuadratic = new BiQuadratic(4.0, 0.0, + Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, + Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY); + final MultivariateRealFunctionMappingAdapter wrapped = + new MultivariateRealFunctionMappingAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper()); + + SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30); + optimizer.setSimplex(new NelderMeadSimplex(new double[][] { + wrapped.boundedToUnbounded(new double[] { 1.5, 2.75 }), + wrapped.boundedToUnbounded(new double[] { 1.5, 2.95 }), + wrapped.boundedToUnbounded(new double[] { 1.7, 2.90 }) + })); + + final RealPointValuePair optimum + = optimizer.optimize(300, wrapped, GoalType.MINIMIZE, + wrapped.boundedToUnbounded(new double[] { 1.5, 2.25 })); + final double[] bounded = wrapped.unboundedToBounded(optimum.getPoint()); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), bounded[0], 2e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), bounded[1], 2e-7); + + } + + @Test + public void testHalfBounded() { + + final BiQuadratic biQuadratic = new BiQuadratic(4.0, 4.0, + 1.0, Double.POSITIVE_INFINITY, + Double.NEGATIVE_INFINITY, 3.0); + final MultivariateRealFunctionMappingAdapter wrapped = + new MultivariateRealFunctionMappingAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper()); + + SimplexOptimizer optimizer = new SimplexOptimizer(1e-13, 1e-30); + optimizer.setSimplex(new NelderMeadSimplex(new double[][] { + wrapped.boundedToUnbounded(new double[] { 1.5, 2.75 }), + wrapped.boundedToUnbounded(new double[] { 1.5, 2.95 }), + wrapped.boundedToUnbounded(new double[] { 1.7, 2.90 }) + })); + + final RealPointValuePair optimum + = optimizer.optimize(200, wrapped, GoalType.MINIMIZE, + wrapped.boundedToUnbounded(new double[] { 1.5, 2.25 })); + final double[] bounded = wrapped.unboundedToBounded(optimum.getPoint()); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), bounded[0], 1e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), bounded[1], 1e-7); + + } + + private static class BiQuadratic implements MultivariateRealFunction { + + private final double xOptimum; + private final double yOptimum; + + private final double xMin; + private final double xMax; + private final double yMin; + private final double yMax; + + public BiQuadratic(final double xOptimum, final double yOptimum, + final double xMin, final double xMax, + final double yMin, final double yMax) { + this.xOptimum = xOptimum; + this.yOptimum = yOptimum; + this.xMin = xMin; + this.xMax = xMax; + this.yMin = yMin; + this.yMax = yMax; + } + + public double value(double[] point) { + + // the function should never be called with out of range points + Assert.assertTrue(point[0] >= xMin); + Assert.assertTrue(point[0] <= xMax); + Assert.assertTrue(point[1] >= yMin); + Assert.assertTrue(point[1] <= yMax); + + final double dx = point[0] - xOptimum; + final double dy = point[1] - yOptimum; + return dx * dx + dy * dy; + + } + + public double[] getLower() { + return new double[] { xMin, yMin }; + } + + public double[] getUpper() { + return new double[] { xMax, yMax }; + } + + public double getBoundedXOptimum() { + return (xOptimum < xMin) ? xMin : ((xOptimum > xMax) ? xMax : xOptimum); + } + + public double getBoundedYOptimum() { + return (yOptimum < yMin) ? yMin : ((yOptimum > yMax) ? yMax : yOptimum); + } + + } + +} diff --git a/src/test/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapterTest.java b/src/test/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapterTest.java new file mode 100644 index 000000000..c1a7a1458 --- /dev/null +++ b/src/test/java/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapterTest.java @@ -0,0 +1,192 @@ +/* + * 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 org.apache.commons.math.analysis.MultivariateRealFunction; +import org.apache.commons.math.optimization.GoalType; +import org.apache.commons.math.optimization.RealPointValuePair; +import org.apache.commons.math.optimization.SimpleRealPointChecker; +import org.junit.Assert; +import org.junit.Test; + +public class MultivariateRealFunctionPenaltyAdapterTest { + + @Test + public void testStartSimplexInsideRange() { + + final BiQuadratic biQuadratic = new BiQuadratic(2.0, 2.5, 1.0, 3.0, 2.0, 3.0); + final MultivariateRealFunctionPenaltyAdapter wrapped = + new MultivariateRealFunctionPenaltyAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper(), + 1000.0, new double[] { 100.0, 100.0 }); + + SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30); + optimizer.setSimplex(new NelderMeadSimplex(new double[] { 1.0, 0.5 })); + + final RealPointValuePair optimum + = optimizer.optimize(300, wrapped, GoalType.MINIMIZE, new double[] { 1.5, 2.25 }); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), optimum.getPoint()[0], 2e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), optimum.getPoint()[1], 2e-7); + + } + + @Test + public void testStartSimplexOutsideRange() { + + final BiQuadratic biQuadratic = new BiQuadratic(2.0, 2.5, 1.0, 3.0, 2.0, 3.0); + final MultivariateRealFunctionPenaltyAdapter wrapped = + new MultivariateRealFunctionPenaltyAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper(), + 1000.0, new double[] { 100.0, 100.0 }); + + SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30); + optimizer.setSimplex(new NelderMeadSimplex(new double[] { 1.0, 0.5 })); + + final RealPointValuePair optimum + = optimizer.optimize(300, wrapped, GoalType.MINIMIZE, new double[] { -1.5, 4.0 }); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), optimum.getPoint()[0], 2e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), optimum.getPoint()[1], 2e-7); + + } + + @Test + public void testOptimumOutsideRange() { + + final BiQuadratic biQuadratic = new BiQuadratic(4.0, 0.0, 1.0, 3.0, 2.0, 3.0); + final MultivariateRealFunctionPenaltyAdapter wrapped = + new MultivariateRealFunctionPenaltyAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper(), + 1000.0, new double[] { 100.0, 100.0 }); + + SimplexOptimizer optimizer = new SimplexOptimizer(new SimpleRealPointChecker(1.0e-11, 1.0e-20)); + optimizer.setSimplex(new NelderMeadSimplex(new double[] { 1.0, 0.5 })); + + final RealPointValuePair optimum + = optimizer.optimize(600, wrapped, GoalType.MINIMIZE, new double[] { -1.5, 4.0 }); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), optimum.getPoint()[0], 2e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), optimum.getPoint()[1], 2e-7); + + } + + @Test + public void testUnbounded() { + + final BiQuadratic biQuadratic = new BiQuadratic(4.0, 0.0, + Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, + Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY); + final MultivariateRealFunctionPenaltyAdapter wrapped = + new MultivariateRealFunctionPenaltyAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper(), + 1000.0, new double[] { 100.0, 100.0 }); + + SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30); + optimizer.setSimplex(new NelderMeadSimplex(new double[] { 1.0, 0.5 })); + + final RealPointValuePair optimum + = optimizer.optimize(300, wrapped, GoalType.MINIMIZE, new double[] { -1.5, 4.0 }); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), optimum.getPoint()[0], 2e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), optimum.getPoint()[1], 2e-7); + + } + + @Test + public void testHalfBounded() { + + final BiQuadratic biQuadratic = new BiQuadratic(4.0, 4.0, + 1.0, Double.POSITIVE_INFINITY, + Double.NEGATIVE_INFINITY, 3.0); + final MultivariateRealFunctionPenaltyAdapter wrapped = + new MultivariateRealFunctionPenaltyAdapter(biQuadratic, + biQuadratic.getLower(), + biQuadratic.getUpper(), + 1000.0, new double[] { 100.0, 100.0 }); + + SimplexOptimizer optimizer = new SimplexOptimizer(new SimpleRealPointChecker(1.0e-10, 1.0e-20)); + optimizer.setSimplex(new NelderMeadSimplex(new double[] { 1.0, 0.5 })); + + final RealPointValuePair optimum + = optimizer.optimize(400, wrapped, GoalType.MINIMIZE, new double[] { -1.5, 4.0 }); + + Assert.assertEquals(biQuadratic.getBoundedXOptimum(), optimum.getPoint()[0], 2e-7); + Assert.assertEquals(biQuadratic.getBoundedYOptimum(), optimum.getPoint()[1], 2e-7); + + } + + private static class BiQuadratic implements MultivariateRealFunction { + + private final double xOptimum; + private final double yOptimum; + + private final double xMin; + private final double xMax; + private final double yMin; + private final double yMax; + + public BiQuadratic(final double xOptimum, final double yOptimum, + final double xMin, final double xMax, + final double yMin, final double yMax) { + this.xOptimum = xOptimum; + this.yOptimum = yOptimum; + this.xMin = xMin; + this.xMax = xMax; + this.yMin = yMin; + this.yMax = yMax; + } + + public double value(double[] point) { + + // the function should never be called with out of range points + Assert.assertTrue(point[0] >= xMin); + Assert.assertTrue(point[0] <= xMax); + Assert.assertTrue(point[1] >= yMin); + Assert.assertTrue(point[1] <= yMax); + + final double dx = point[0] - xOptimum; + final double dy = point[1] - yOptimum; + return dx * dx + dy * dy; + + } + + public double[] getLower() { + return new double[] { xMin, yMin }; + } + + public double[] getUpper() { + return new double[] { xMax, yMax }; + } + + public double getBoundedXOptimum() { + return (xOptimum < xMin) ? xMin : ((xOptimum > xMax) ? xMax : xOptimum); + } + + public double getBoundedYOptimum() { + return (yOptimum < yMin) ? yMin : ((yOptimum > yMax) ? yMax : yOptimum); + } + + } + +}