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
This commit is contained in:
Luc Maisonobe 2011-11-10 20:08:55 +00:00
parent dadf9a70a0
commit bb201d926e
7 changed files with 926 additions and 5 deletions

View File

@ -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;
/**
* <p>Adapter for mapping bounded {@link MultivariateRealFunction} to unbounded ones.</p>
*
* <p>
* This adapter can be used to wrap functions subject to simple bounds on
* parameters so they can be used by optimizers that do <em>not</em> directly
* support simple bounds.
* </p>
* <p>
* 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.
* </p>
* <p>
* 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!</p>
* <p>
* 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.
* </p>
*
* @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.
* <p>
* This method simply bounds the unbounded point using the mappings
* set up at construction and calls the underlying function using
* the bounded point.
* </p>
* @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);
}
}
}

View File

@ -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;
/**
* <p>Adapter extending bounded {@link MultivariateRealFunction} to an unbouded
* domain using a penalty function.</p>
*
* <p>
* This adapter can be used to wrap functions subject to simple bounds on
* parameters so they can be used by optimizers that do <em>not</em> directly
* support simple bounds.
* </p>
* <p>
* 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.
* </p>
* <p>
* 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.
* </p>
*
* @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.
* <p>
* 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:
* <pre>
* penalty(point) = offset + &sum;<sub>i</sub>[scale[i] * &radic;|point[i]-boundary[i]|]
* </pre>
* where indices i correspond to all the components that violates their boundaries.
* </p>
* <p>
* 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.
* </p>
* <p>
* 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.
* </p>
* @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.
* <p>
* This method simply bounds the unbounded point using the mappings
* set up at construction and calls the underlying function using
* the bounded point.
* </p>
* @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);
}
}

View File

@ -66,8 +66,20 @@ import org.apache.commons.math.optimization.MultivariateRealOptimizer;
* previous and current simplex to the convergence checker, not the best
* ones.
* </p>
* <p>
* 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}.
* </p>
*
* @see AbstractSimplex
* @see MultivariateRealFunctionMappingAdapter
* @see MultivariateRealFunctionPenaltyAdapter
* @see CMAESOptimizer
* @see BOBYQAOptimizer
* @version $Id$
* @since 3.0
*/

View File

@ -52,6 +52,11 @@ 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="luc" type="fix" issue="MATH-196" >
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.
</action>
<action dev="luc" type="fix" issue="MATH-702" >
CMA-ES optimizer input sigma is now consistent with boundaries range units.
</action>

View File

@ -154,11 +154,45 @@
already provided by the <code>minimizes</code> method).
</p>
<p>
The <code>direct</code> package provides two solvers. The first one is the classical
<a href="../apidocs/org/apache/commons/math/optimization/direct/NelderMead.html">
Nelder-Mead</a> method. The second one is Virginia Torczon's
<a href="../apidocs/org/apache/commons/math/optimization/direct/MultiDirectional.html">
multi-directional</a> method.
The <code>direct</code> package provides four solvers:
<ul>
<li>the classical <a
href="../apidocs/org/apache/commons/math/optimization/direct/NelderMeadSimplex.html">
Nelder-Mead</a> method,</li>
<li>Virginia Torczon's <a
href="../apidocs/org/apache/commons/math/optimization/direct/MultiDirectionalSimplex.html">
multi-directional</a> method,</li>
<li>Nikolaus Hansen's <a
href="../apidocs/org/apache/commons/math/optimization/direct/CMAESOptimizer.html">
</a>Covariance Matrix Adaptation Evolution Strategy (CMA-ES),</li>
<li>Mike Powell's <a
href="../apidocs/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.html">
BOBYQA</a> method.
</li>
</ul>
</p>
<p>
The first two simplex-based methods do not handle simple bounds constraints by themselves.
However there are two adapters (<a
href="../apidocs/org/apache/commons/math/optimization/direct/MultivariateRealFunctionMappingAdapter.html">
MultivariateRealFunctionMappingAdapter</a> and <a
href="../apidocs/org/apache/commons/math/optimization/direct/MultivariateRealFunctionPenaltyAdapter.html">
MultivariateRealFunctionPenaltyAdapter</a>) 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.
</p>
<p>
The last methods do handle simple bounds constraints directly, so the adapters are not needed
with them.
</p>
</subsection>
<subsection name="12.5 General Case" href="general">

View File

@ -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);
}
}
}

View File

@ -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);
}
}
}