← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/optimization/direct/DirectSearchOptimizer.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/optimization/direct/DirectSearchOptimizer.java
 
1
/*
 
2
 * Licensed to the Apache Software Foundation (ASF) under one or more
 
3
 * contributor license agreements.  See the NOTICE file distributed with
 
4
 * this work for additional information regarding copyright ownership.
 
5
 * The ASF licenses this file to You under the Apache License, Version 2.0
 
6
 * (the "License"); you may not use this file except in compliance with
 
7
 * the License.  You may obtain a copy of the License at
 
8
 *
 
9
 *      http://www.apache.org/licenses/LICENSE-2.0
 
10
 *
 
11
 * Unless required by applicable law or agreed to in writing, software
 
12
 * 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 and
 
15
 * limitations under the License.
 
16
 */
 
17
 
 
18
package org.apache.commons.math.optimization.direct;
 
19
 
 
20
import java.util.Arrays;
 
21
import java.util.Comparator;
 
22
 
 
23
import org.apache.commons.math.FunctionEvaluationException;
 
24
import org.apache.commons.math.MathRuntimeException;
 
25
import org.apache.commons.math.MaxEvaluationsExceededException;
 
26
import org.apache.commons.math.MaxIterationsExceededException;
 
27
import org.apache.commons.math.analysis.MultivariateRealFunction;
 
28
import org.apache.commons.math.optimization.GoalType;
 
29
import org.apache.commons.math.optimization.MultivariateRealOptimizer;
 
30
import org.apache.commons.math.optimization.OptimizationException;
 
31
import org.apache.commons.math.optimization.RealConvergenceChecker;
 
32
import org.apache.commons.math.optimization.RealPointValuePair;
 
33
import org.apache.commons.math.optimization.SimpleScalarValueChecker;
 
34
 
 
35
/** 
 
36
 * This class implements simplex-based direct search optimization
 
37
 * algorithms.
 
38
 *
 
39
 * <p>Direct search methods only use objective function values, they don't
 
40
 * need derivatives and don't either try to compute approximation of
 
41
 * the derivatives. According to a 1996 paper by Margaret H. Wright
 
42
 * (<a href="http://cm.bell-labs.com/cm/cs/doc/96/4-02.ps.gz">Direct
 
43
 * Search Methods: Once Scorned, Now Respectable</a>), they are used
 
44
 * when either the computation of the derivative is impossible (noisy
 
45
 * functions, unpredictable discontinuities) or difficult (complexity,
 
46
 * computation cost). In the first cases, rather than an optimum, a
 
47
 * <em>not too bad</em> point is desired. In the latter cases, an
 
48
 * optimum is desired but cannot be reasonably found. In all cases
 
49
 * direct search methods can be useful.</p>
 
50
 *
 
51
 * <p>Simplex-based direct search methods are based on comparison of
 
52
 * the objective function values at the vertices of a simplex (which is a
 
53
 * set of n+1 points in dimension n) that is updated by the algorithms
 
54
 * steps.<p>
 
55
 *
 
56
 * <p>The initial configuration of the simplex can be set using either
 
57
 * {@link #setStartConfiguration(double[])} or {@link
 
58
 * #setStartConfiguration(double[][])}. If neither method has been called
 
59
 * before optimization is attempted, an explicit call to the first method
 
60
 * with all steps set to +1 is triggered, thus building a default
 
61
 * configuration from a unit hypercube. Each call to {@link
 
62
 * #optimize(MultivariateRealFunction, GoalType, double[]) optimize} will reuse
 
63
 * the current start configuration and move it such that its first vertex
 
64
 * is at the provided start point of the optimization. If the same optimizer
 
65
 * is used to solve different problems and the number of parameters change,
 
66
 * the start configuration <em>must</em> be reset or a dimension mismatch
 
67
 * will occur.</p>
 
68
 *
 
69
 * <p>If {@link #setConvergenceChecker(RealConvergenceChecker)} is not called,
 
70
 * a default {@link SimpleScalarValueChecker} is used.</p>
 
71
 *
 
72
 * <p>Convergence is checked by providing the <em>worst</em> points of
 
73
 * previous and current simplex to the convergence checker, not the best ones.</p>
 
74
 *
 
75
 * <p>This class is the base class performing the boilerplate simplex
 
76
 * initialization and handling. The simplex update by itself is
 
77
 * performed by the derived classes according to the implemented
 
78
 * algorithms.</p>
 
79
 *
 
80
 * implements MultivariateRealOptimizer since 2.0
 
81
 * 
 
82
 * @see MultivariateRealFunction
 
83
 * @see NelderMead
 
84
 * @see MultiDirectional
 
85
 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
 
86
 * @since 1.2
 
87
 */
 
88
public abstract class DirectSearchOptimizer implements MultivariateRealOptimizer {
 
89
 
 
90
    /** Simplex. */
 
91
    protected RealPointValuePair[] simplex;
 
92
 
 
93
    /** Objective function. */
 
94
    private MultivariateRealFunction f;
 
95
 
 
96
    /** Convergence checker. */
 
97
    private RealConvergenceChecker checker;
 
98
 
 
99
    /** Maximal number of iterations allowed. */
 
100
    private int maxIterations;
 
101
 
 
102
    /** Number of iterations already performed. */
 
103
    private int iterations;
 
104
 
 
105
    /** Maximal number of evaluations allowed. */
 
106
    private int maxEvaluations;
 
107
 
 
108
    /** Number of evaluations already performed. */
 
109
    private int evaluations;
 
110
 
 
111
    /** Start simplex configuration. */
 
112
    private double[][] startConfiguration;
 
113
 
 
114
    /** Simple constructor.
 
115
     */
 
116
    protected DirectSearchOptimizer() {
 
117
        setConvergenceChecker(new SimpleScalarValueChecker());
 
118
        setMaxIterations(Integer.MAX_VALUE);
 
119
        setMaxEvaluations(Integer.MAX_VALUE);
 
120
    }
 
121
 
 
122
    /** Set start configuration for simplex.
 
123
     * <p>The start configuration for simplex is built from a box parallel to
 
124
     * the canonical axes of the space. The simplex is the subset of vertices
 
125
     * of a box parallel to the canonical axes. It is built as the path followed
 
126
     * while traveling from one vertex of the box to the diagonally opposite
 
127
     * vertex moving only along the box edges. The first vertex of the box will
 
128
     * be located at the start point of the optimization.</p>
 
129
     * <p>As an example, in dimension 3 a simplex has 4 vertices. Setting the
 
130
     * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
 
131
     * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
 
132
     * The first vertex would be set to the start point at (1, 1, 1) and the
 
133
     * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).</p>
 
134
     * @param steps steps along the canonical axes representing box edges,
 
135
     * they may be negative but not null
 
136
     * @exception IllegalArgumentException if one step is null
 
137
     */
 
138
    public void setStartConfiguration(final double[] steps)
 
139
        throws IllegalArgumentException {
 
140
        // only the relative position of the n final vertices with respect
 
141
        // to the first one are stored
 
142
        final int n = steps.length;
 
143
        startConfiguration = new double[n][n];
 
144
        for (int i = 0; i < n; ++i) {
 
145
            final double[] vertexI = startConfiguration[i];
 
146
            for (int j = 0; j < i + 1; ++j) {
 
147
                if (steps[j] == 0.0) {
 
148
                    throw MathRuntimeException.createIllegalArgumentException(
 
149
                            "equals vertices {0} and {1} in simplex configuration",
 
150
                            j, j + 1);
 
151
                }
 
152
                System.arraycopy(steps, 0, vertexI, 0, j + 1);
 
153
            }
 
154
        }
 
155
    }
 
156
 
 
157
    /** Set start configuration for simplex.
 
158
     * <p>The real initial simplex will be set up by moving the reference
 
159
     * simplex such that its first point is located at the start point of the
 
160
     * optimization.</p>
 
161
     * @param referenceSimplex reference simplex
 
162
     * @exception IllegalArgumentException if the reference simplex does not
 
163
     * contain at least one point, or if there is a dimension mismatch
 
164
     * in the reference simplex or if one of its vertices is duplicated
 
165
     */
 
166
    public void setStartConfiguration(final double[][] referenceSimplex)
 
167
        throws IllegalArgumentException {
 
168
 
 
169
        // only the relative position of the n final vertices with respect
 
170
        // to the first one are stored
 
171
        final int n = referenceSimplex.length - 1;
 
172
        if (n < 0) {
 
173
            throw MathRuntimeException.createIllegalArgumentException(
 
174
                    "simplex must contain at least one point");
 
175
        }
 
176
        startConfiguration = new double[n][n];
 
177
        final double[] ref0 = referenceSimplex[0];
 
178
 
 
179
        // vertices loop
 
180
        for (int i = 0; i < n + 1; ++i) {
 
181
 
 
182
            final double[] refI = referenceSimplex[i];
 
183
 
 
184
            // safety checks
 
185
            if (refI.length != n) {
 
186
                throw MathRuntimeException.createIllegalArgumentException(
 
187
                        "dimension mismatch {0} != {1}",
 
188
                        refI.length, n);
 
189
            }
 
190
            for (int j = 0; j < i; ++j) {
 
191
                final double[] refJ = referenceSimplex[j];
 
192
                boolean allEquals = true;
 
193
                for (int k = 0; k < n; ++k) {
 
194
                    if (refI[k] != refJ[k]) {
 
195
                        allEquals = false;
 
196
                        break;
 
197
                    }
 
198
                }
 
199
                if (allEquals) {
 
200
                    throw MathRuntimeException.createIllegalArgumentException(
 
201
                            "equals vertices {0} and {1} in simplex configuration",
 
202
                            i, j);
 
203
                }
 
204
            }
 
205
 
 
206
            // store vertex i position relative to vertex 0 position
 
207
            if (i > 0) {
 
208
                final double[] confI = startConfiguration[i - 1];
 
209
                for (int k = 0; k < n; ++k) {
 
210
                    confI[k] = refI[k] - ref0[k];
 
211
                }
 
212
            }
 
213
 
 
214
        }
 
215
 
 
216
    }
 
217
 
 
218
    /** {@inheritDoc} */
 
219
    public void setMaxIterations(int maxIterations) {
 
220
        this.maxIterations = maxIterations;
 
221
    }
 
222
 
 
223
    /** {@inheritDoc} */
 
224
    public int getMaxIterations() {
 
225
        return maxIterations;
 
226
    }
 
227
 
 
228
    /** {@inheritDoc} */
 
229
    public void setMaxEvaluations(int maxEvaluations) {
 
230
        this.maxEvaluations = maxEvaluations;
 
231
    }
 
232
 
 
233
    /** {@inheritDoc} */
 
234
    public int getMaxEvaluations() {
 
235
        return maxEvaluations;
 
236
    }
 
237
 
 
238
    /** {@inheritDoc} */
 
239
    public int getIterations() {
 
240
        return iterations;
 
241
    }
 
242
 
 
243
    /** {@inheritDoc} */
 
244
    public int getEvaluations() {
 
245
        return evaluations;
 
246
    }
 
247
 
 
248
    /** {@inheritDoc} */
 
249
    public void setConvergenceChecker(RealConvergenceChecker checker) {
 
250
        this.checker = checker;
 
251
    }
 
252
 
 
253
    /** {@inheritDoc} */
 
254
    public RealConvergenceChecker getConvergenceChecker() {
 
255
        return checker;
 
256
    }
 
257
 
 
258
    /** {@inheritDoc} */
 
259
    public RealPointValuePair optimize(final MultivariateRealFunction f,
 
260
                                         final GoalType goalType,
 
261
                                         final double[] startPoint)
 
262
        throws FunctionEvaluationException, OptimizationException,
 
263
        IllegalArgumentException {
 
264
 
 
265
        if (startConfiguration == null) {
 
266
            // no initial configuration has been set up for simplex
 
267
            // build a default one from a unit hypercube
 
268
            final double[] unit = new double[startPoint.length];
 
269
            Arrays.fill(unit, 1.0);
 
270
            setStartConfiguration(unit);
 
271
        }
 
272
 
 
273
        this.f = f;
 
274
        final Comparator<RealPointValuePair> comparator =
 
275
            new Comparator<RealPointValuePair>() {
 
276
                public int compare(final RealPointValuePair o1,
 
277
                                   final RealPointValuePair o2) {
 
278
                    final double v1 = o1.getValue();
 
279
                    final double v2 = o2.getValue();
 
280
                    return (goalType == GoalType.MINIMIZE) ?
 
281
                            Double.compare(v1, v2) : Double.compare(v2, v1);
 
282
                }
 
283
            };
 
284
 
 
285
        // initialize search
 
286
        iterations  = 0;
 
287
        evaluations = 0;
 
288
        buildSimplex(startPoint);
 
289
        evaluateSimplex(comparator);
 
290
 
 
291
        RealPointValuePair[] previous = new RealPointValuePair[simplex.length];
 
292
        while (true) {
 
293
 
 
294
            if (iterations > 0) {
 
295
                boolean converged = true;
 
296
                for (int i = 0; i < simplex.length; ++i) {
 
297
                    converged &= checker.converged(iterations, previous[i], simplex[i]);
 
298
                }
 
299
                if (converged) {
 
300
                    // we have found an optimum
 
301
                    return simplex[0];
 
302
                }
 
303
            }
 
304
 
 
305
            // we still need to search
 
306
            System.arraycopy(simplex, 0, previous, 0, simplex.length);
 
307
            iterateSimplex(comparator);
 
308
 
 
309
        }
 
310
 
 
311
    }
 
312
 
 
313
    /** Increment the iterations counter by 1.
 
314
     * @exception OptimizationException if the maximal number
 
315
     * of iterations is exceeded
 
316
     */
 
317
    protected void incrementIterationsCounter()
 
318
        throws OptimizationException {
 
319
        if (++iterations > maxIterations) {
 
320
            throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
 
321
        }
 
322
    }
 
323
 
 
324
    /** Compute the next simplex of the algorithm.
 
325
     * @param comparator comparator to use to sort simplex vertices from best to worst
 
326
     * @exception FunctionEvaluationException if the function cannot be evaluated at
 
327
     * some point
 
328
     * @exception OptimizationException if the algorithm fails to converge
 
329
     * @exception IllegalArgumentException if the start point dimension is wrong
 
330
     */
 
331
    protected abstract void iterateSimplex(final Comparator<RealPointValuePair> comparator)
 
332
        throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
 
333
 
 
334
    /** Evaluate the objective function on one point.
 
335
     * <p>A side effect of this method is to count the number of
 
336
     * function evaluations</p>
 
337
     * @param x point on which the objective function should be evaluated
 
338
     * @return objective function value at the given point
 
339
     * @exception FunctionEvaluationException if no value can be computed for the
 
340
     * parameters or if the maximal number of evaluations is exceeded
 
341
     * @exception IllegalArgumentException if the start point dimension is wrong
 
342
     */
 
343
    protected double evaluate(final double[] x)
 
344
        throws FunctionEvaluationException, IllegalArgumentException {
 
345
        if (++evaluations > maxEvaluations) {
 
346
            throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations),
 
347
                                                  x);
 
348
        }
 
349
        return f.value(x);
 
350
    }
 
351
 
 
352
    /** Build an initial simplex.
 
353
     * @param startPoint the start point for optimization
 
354
     * @exception IllegalArgumentException if the start point does not match
 
355
     * simplex dimension
 
356
     */
 
357
    private void buildSimplex(final double[] startPoint)
 
358
        throws IllegalArgumentException {
 
359
 
 
360
        final int n = startPoint.length;
 
361
        if (n != startConfiguration.length) {
 
362
            throw MathRuntimeException.createIllegalArgumentException(
 
363
                    "dimension mismatch {0} != {1}",
 
364
                    n, startConfiguration.length);
 
365
        }
 
366
 
 
367
        // set first vertex
 
368
        simplex = new RealPointValuePair[n + 1];
 
369
        simplex[0] = new RealPointValuePair(startPoint, Double.NaN);
 
370
 
 
371
        // set remaining vertices
 
372
        for (int i = 0; i < n; ++i) {
 
373
            final double[] confI   = startConfiguration[i];
 
374
            final double[] vertexI = new double[n];
 
375
            for (int k = 0; k < n; ++k) {
 
376
                vertexI[k] = startPoint[k] + confI[k];
 
377
            }
 
378
            simplex[i + 1] = new RealPointValuePair(vertexI, Double.NaN);
 
379
        }
 
380
 
 
381
    }
 
382
 
 
383
    /** Evaluate all the non-evaluated points of the simplex.
 
384
     * @param comparator comparator to use to sort simplex vertices from best to worst
 
385
     * @exception FunctionEvaluationException if no value can be computed for the parameters
 
386
     * @exception OptimizationException if the maximal number of evaluations is exceeded
 
387
     */
 
388
    protected void evaluateSimplex(final Comparator<RealPointValuePair> comparator)
 
389
        throws FunctionEvaluationException, OptimizationException {
 
390
 
 
391
        // evaluate the objective function at all non-evaluated simplex points
 
392
        for (int i = 0; i < simplex.length; ++i) {
 
393
            final RealPointValuePair vertex = simplex[i];
 
394
            final double[] point = vertex.getPointRef();
 
395
            if (Double.isNaN(vertex.getValue())) {
 
396
                simplex[i] = new RealPointValuePair(point, evaluate(point), false);
 
397
            }
 
398
        }
 
399
 
 
400
        // sort the simplex from best to worst
 
401
        Arrays.sort(simplex, comparator);
 
402
 
 
403
    }
 
404
 
 
405
    /** Replace the worst point of the simplex by a new point.
 
406
     * @param pointValuePair point to insert
 
407
     * @param comparator comparator to use to sort simplex vertices from best to worst
 
408
     */
 
409
    protected void replaceWorstPoint(RealPointValuePair pointValuePair,
 
410
                                     final Comparator<RealPointValuePair> comparator) {
 
411
        int n = simplex.length - 1;
 
412
        for (int i = 0; i < n; ++i) {
 
413
            if (comparator.compare(simplex[i], pointValuePair) > 0) {
 
414
                RealPointValuePair tmp = simplex[i];
 
415
                simplex[i]         = pointValuePair;
 
416
                pointValuePair     = tmp;
 
417
            }
 
418
        }
 
419
        simplex[n] = pointValuePair;
 
420
    }
 
421
 
 
422
}