1/*2* Licensed to the Apache Software Foundation (ASF) under one or more3* contributor license agreements. See the NOTICE file distributed with4* this work for additional information regarding copyright ownership.5* The ASF licenses this file to You under the Apache License, Version 2.06* (the "License"); you may not use this file except in compliance with7* the License. You may obtain a copy of the License at8*9* http://www.apache.org/licenses/LICENSE-2.010*11* Unless required by applicable law or agreed to in writing, software12* distributed under the License is distributed on an "AS IS" BASIS,13* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.14* See the License for the specific language governing permissions and15* limitations under the License.16*/17 18packageorg.apache.commons.math3.optimization.direct; 19 20importorg.apache.commons.math3.analysis.MultivariateFunction; 21importorg.apache.commons.math3.exception.DimensionMismatchException; 22importorg.apache.commons.math3.exception.NumberIsTooSmallException; 23importorg.apache.commons.math3.util.FastMath; 24importorg.apache.commons.math3.util.MathUtils; 25 26/**27* <p>Adapter extending bounded {@link MultivariateFunction} to an unbouded28* domain using a penalty function.</p>29*30* <p>31* This adapter can be used to wrap functions subject to simple bounds on32* parameters so they can be used by optimizers that do <em>not</em> directly33* support simple bounds.34* </p>35* <p>36* The principle is that the user function that will be wrapped will see its37* parameters bounded as required, i.e when its {@code value} method is called38* with argument array {@code point}, the elements array will fulfill requirement39* {@code lower[i] <= point[i] <= upper[i]} for all i. Some of the components40* may be unbounded or bounded only on one side if the corresponding bound is41* set to an infinite value. The optimizer will not manage the user function by42* itself, but it will handle this adapter and it is this adapter that will take43* care the bounds are fulfilled. The adapter {@link #value(double[])} method will44* be called by the optimizer with unbound parameters, and the adapter will check45* if the parameters is within range or not. If it is in range, then the underlying46* user function will be called, and if it is not the value of a penalty function47* will be returned instead.48* </p>49* <p>50* This adapter is only a poor man solution to simple bounds optimization constraints51* that can be used with simple optimizers like {@link SimplexOptimizer} with {@link52* NelderMeadSimplex} or {@link MultiDirectionalSimplex}. A better solution is to use53* an optimizer that directly supports simple bounds like {@link CMAESOptimizer} or54* {@link BOBYQAOptimizer}. One caveat of this poor man solution is that if start point55* or start simplex is completely outside of the allowed range, only the penalty function56* is used, and the optimizer may converge without ever entering the range.57* </p>58*59* @see MultivariateFunctionMappingAdapter60*61* @deprecated As of 3.1 (to be removed in 4.0).62* @since 3.063*/64 65 @Deprecated 66publicclassMultivariateFunctionPenaltyAdapterimplementsMultivariateFunction { 67 68/** Underlying bounded function. */69privatefinalMultivariateFunction bounded; 70 71/** Lower bounds. */72privatefinaldouble[] lower; 73 74/** Upper bounds. */75privatefinaldouble[] upper; 76 77/** Penalty offset. */78privatefinaldoubleoffset; 79 80/** Penalty scales. */81privatefinaldouble[] scale; 82 83/** Simple constructor.84* <p>85* When the optimizer provided points are out of range, the value of the86* penalty function will be used instead of the value of the underlying87* function. In order for this penalty to be effective in rejecting this88* point during the optimization process, the penalty function value should89* be defined with care. This value is computed as:90* <pre>91* penalty(point) = offset + ∑<sub>i</sub>[scale[i] * √|point[i]-boundary[i]|]92* </pre>93* where indices i correspond to all the components that violates their boundaries.94* </p>95* <p>96* So when attempting a function minimization, offset should be larger than97* the maximum expected value of the underlying function and scale components98* should all be positive. When attempting a function maximization, offset99* should be lesser than the minimum expected value of the underlying function100* and scale components should all be negative.101* minimization, and lesser than the minimum expected value of the underlying102* function when attempting maximization.103* </p>104* <p>105* These choices for the penalty function have two properties. First, all out106* of range points will return a function value that is worse than the value107* returned by any in range point. Second, the penalty is worse for large108* boundaries violation than for small violations, so the optimizer has an hint109* about the direction in which it should search for acceptable points.110* </p>111* @param bounded bounded function112* @param lower lower bounds for each element of the input parameters array113* (some elements may be set to {@code Double.NEGATIVE_INFINITY} for114* unbounded values)115* @param upper upper bounds for each element of the input parameters array116* (some elements may be set to {@code Double.POSITIVE_INFINITY} for117* unbounded values)118* @param offset base offset of the penalty function119* @param scale scale of the penalty function120* @exception DimensionMismatchException if lower bounds, upper bounds and121* scales are not consistent, either according to dimension or to bounadary122* values123*/124publicMultivariateFunctionPenaltyAdapter(finalMultivariateFunction bounded, 125finaldouble[] lower,finaldouble[] upper, 126finaldoubleoffset,finaldouble[] scale) { 127 128// safety checks129 MathUtils.checkNotNull(lower); 130 MathUtils.checkNotNull(upper); 131 MathUtils.checkNotNull(scale); 132if(lower.length != upper.length) { 133thrownewDimensionMismatchException(lower.length, upper.length); 134 } 135if(lower.length != scale.length) { 136thrownewDimensionMismatchException(lower.length, scale.length); 137 } 138for(inti = 0; i < lower.length; ++i) { 139// note the following test is written in such a way it also fails for NaN140if(!(upper[i] >= lower[i])) { 141thrownewNumberIsTooSmallException(upper[i], lower[i],true); 142 } 143 } 144 145this.bounded = bounded; 146this.lower = lower.clone(); 147this.upper = upper.clone(); 148this.offset = offset; 149this.scale = scale.clone(); 150 151 } 152 153/** Compute the underlying function value from an unbounded point.154* <p>155* This method simply returns the value of the underlying function156* if the unbounded point already fulfills the bounds, and compute157* a replacement value using the offset and scale if bounds are158* violated, without calling the function at all.159* </p>160* @param point unbounded point161* @return either underlying function value or penalty function value162*/163publicdoublevalue(double[] point) { 164 165for(inti = 0; i < scale.length; ++i) { 166if((point[i] < lower[i]) || (point[i] > upper[i])) { 167// bound violation starting at this component168doublesum = 0; 169for(intj = i; j < scale.length; ++j) { 170finaldoubleovershoot; 171if(point[j] < lower[j]) { 172 overshoot = scale[j] * (lower[j] - point[j]); 173 }elseif(point[j] > upper[j]) { 174 overshoot = scale[j] * (point[j] - upper[j]); 175 }else{ 176 overshoot = 0; 177 } 178 sum += FastMath.sqrt(overshoot); 179 } 180returnoffset + sum; 181 } 182 } 183 184// all boundaries are fulfilled, we are in the expected185// domain of the underlying function186returnbounded.value(point); 187 188 } 189 190 }