001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.optimization.direct;
019    
020    import java.util.Comparator;
021    
022    import org.apache.commons.math.FunctionEvaluationException;
023    import org.apache.commons.math.optimization.OptimizationException;
024    import org.apache.commons.math.optimization.RealConvergenceChecker;
025    import org.apache.commons.math.optimization.RealPointValuePair;
026    
027    /**
028     * This class implements the multi-directional direct search method.
029     *
030     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
031     * @see NelderMead
032     * @since 1.2
033     */
034    public class MultiDirectional extends DirectSearchOptimizer {
035    
036        /** Expansion coefficient. */
037        private final double khi;
038    
039        /** Contraction coefficient. */
040        private final double gamma;
041    
042        /** Build a multi-directional optimizer with default coefficients.
043         * <p>The default values are 2.0 for khi and 0.5 for gamma.</p>
044         */
045        public MultiDirectional() {
046            this.khi   = 2.0;
047            this.gamma = 0.5;
048        }
049    
050        /** Build a multi-directional optimizer with specified coefficients.
051         * @param khi expansion coefficient
052         * @param gamma contraction coefficient
053         */
054        public MultiDirectional(final double khi, final double gamma) {
055            this.khi   = khi;
056            this.gamma = gamma;
057        }
058    
059        /** {@inheritDoc} */
060        @Override
061        protected void iterateSimplex(final Comparator<RealPointValuePair> comparator)
062            throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
063    
064            final RealConvergenceChecker checker = getConvergenceChecker();
065            while (true) {
066    
067                incrementIterationsCounter();
068    
069                // save the original vertex
070                final RealPointValuePair[] original = simplex;
071                final RealPointValuePair best = original[0];
072    
073                // perform a reflection step
074                final RealPointValuePair reflected = evaluateNewSimplex(original, 1.0, comparator);
075                if (comparator.compare(reflected, best) < 0) {
076    
077                    // compute the expanded simplex
078                    final RealPointValuePair[] reflectedSimplex = simplex;
079                    final RealPointValuePair expanded = evaluateNewSimplex(original, khi, comparator);
080                    if (comparator.compare(reflected, expanded) <= 0) {
081                        // accept the reflected simplex
082                        simplex = reflectedSimplex;
083                    }
084    
085                    return;
086    
087                }
088    
089                // compute the contracted simplex
090                final RealPointValuePair contracted = evaluateNewSimplex(original, gamma, comparator);
091                if (comparator.compare(contracted, best) < 0) {
092                    // accept the contracted simplex
093                    return;
094                }
095    
096                // check convergence
097                final int iter = getIterations();
098                boolean converged = true;
099                for (int i = 0; i < simplex.length; ++i) {
100                    converged &= checker.converged(iter, original[i], simplex[i]);
101                }
102                if (converged) {
103                    return;
104                }
105    
106            }
107    
108        }
109    
110        /** Compute and evaluate a new simplex.
111         * @param original original simplex (to be preserved)
112         * @param coeff linear coefficient
113         * @param comparator comparator to use to sort simplex vertices from best to poorest
114         * @return best point in the transformed simplex
115         * @exception FunctionEvaluationException if the function cannot be evaluated at some point
116         * @exception OptimizationException if the maximal number of evaluations is exceeded
117         */
118        private RealPointValuePair evaluateNewSimplex(final RealPointValuePair[] original,
119                                                  final double coeff,
120                                                  final Comparator<RealPointValuePair> comparator)
121            throws FunctionEvaluationException, OptimizationException {
122    
123            final double[] xSmallest = original[0].getPointRef();
124            final int n = xSmallest.length;
125    
126            // create the linearly transformed simplex
127            simplex = new RealPointValuePair[n + 1];
128            simplex[0] = original[0];
129            for (int i = 1; i <= n; ++i) {
130                final double[] xOriginal    = original[i].getPointRef();
131                final double[] xTransformed = new double[n];
132                for (int j = 0; j < n; ++j) {
133                    xTransformed[j] = xSmallest[j] + coeff * (xSmallest[j] - xOriginal[j]);
134                }
135                simplex[i] = new RealPointValuePair(xTransformed, Double.NaN, false);
136            }
137    
138            // evaluate it
139            evaluateSimplex(comparator);
140            return simplex[0];
141    
142        }
143    
144    }