← 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/univariate/BrentOptimizer.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/univariate/BrentOptimizer.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
package org.apache.commons.math.optimization.univariate;
 
18
 
 
19
import org.apache.commons.math.FunctionEvaluationException;
 
20
import org.apache.commons.math.MaxIterationsExceededException;
 
21
import org.apache.commons.math.analysis.UnivariateRealFunction;
 
22
import org.apache.commons.math.optimization.GoalType;
 
23
 
 
24
/**
 
25
 * Implements Richard Brent's algorithm (from his book "Algorithms for
 
26
 * Minimization without Derivatives", p. 79) for finding minima of real
 
27
 * univariate functions.
 
28
 *  
 
29
 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
 
30
 * @since 2.0
 
31
 */
 
32
public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
 
33
    
 
34
    /**
 
35
     * Golden section.
 
36
     */
 
37
    private static final double c = 0.5 * (3 - Math.sqrt(5));
 
38
 
 
39
    /**
 
40
     * Construct a solver.
 
41
     */
 
42
    public BrentOptimizer() {
 
43
        super(100, 1E-10);
 
44
    }
 
45
 
 
46
    /** {@inheritDoc} */
 
47
    public double optimize(final UnivariateRealFunction f, final GoalType goalType,
 
48
                           final double min, final double max, final double startValue)
 
49
        throws MaxIterationsExceededException, FunctionEvaluationException {
 
50
        return optimize(f, goalType, min, max);
 
51
    }
 
52
    
 
53
    /** {@inheritDoc} */
 
54
    public double optimize(final UnivariateRealFunction f, final GoalType goalType,
 
55
                           final double min, final double max)
 
56
        throws MaxIterationsExceededException, FunctionEvaluationException {
 
57
        clearResult();
 
58
        return localMin(f, goalType, min, max, relativeAccuracy, absoluteAccuracy);
 
59
    }
 
60
    
 
61
    /**
 
62
     * Find the minimum of the function {@code f} within the interval {@code (a, b)}.
 
63
     *
 
64
     * If the function {@code f} is defined on the interval {@code (a, b)}, then
 
65
     * this method finds an approximation {@code x} to the point at which {@code f}
 
66
     * attains its minimum.<br/>
 
67
     * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} and
 
68
     * {@code f} is never evaluated at two points closer together than {@code tol}.
 
69
     * {@code eps} should be no smaller than <em>2 macheps</em> and preferable not
 
70
     * much less than <em>sqrt(macheps)</em>, where <em>macheps</em> is the relative
 
71
     * machine precision. {@code t} should be positive.
 
72
     * @param f the function to solve
 
73
     * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
 
74
     * or {@link GoalType#MINIMIZE}
 
75
     * @param a Lower bound of the interval
 
76
     * @param b Higher bound of the interval
 
77
     * @param eps Relative accuracy
 
78
     * @param t Absolute accuracy
 
79
     * @return the point at which the function is minimal.
 
80
     * @throws MaxIterationsExceededException if the maximum iteration count
 
81
     * is exceeded.
 
82
     * @throws FunctionEvaluationException if an error occurs evaluating
 
83
     * the function. 
 
84
     */
 
85
    private double localMin(final UnivariateRealFunction f, final GoalType goalType,
 
86
                            double a, double b, final double eps, final double t)
 
87
        throws MaxIterationsExceededException, FunctionEvaluationException {
 
88
        double x = a + c * (b - a);
 
89
        double v = x;
 
90
        double w = x;
 
91
        double e = 0;
 
92
        double fx = computeObjectiveValue(f, x);
 
93
        if (goalType == GoalType.MAXIMIZE) {
 
94
            fx = -fx;
 
95
        }
 
96
        double fv = fx;
 
97
        double fw = fx;
 
98
 
 
99
        int count = 0;
 
100
        while (count < maximalIterationCount) {
 
101
            double m = 0.5 * (a + b);
 
102
            double tol = eps * Math.abs(x) + t;
 
103
            double t2 = 2 * tol;
 
104
 
 
105
            // Check stopping criterion.
 
106
            if (Math.abs(x - m) > t2 - 0.5 * (b - a)) {
 
107
                double p = 0;
 
108
                double q = 0;
 
109
                double r = 0;
 
110
                double d = 0;
 
111
                double u = 0;
 
112
 
 
113
                if (Math.abs(e) > tol) { // Fit parabola.
 
114
                    r = (x - w) * (fx - fv);
 
115
                    q = (x - v) * (fx - fw);
 
116
                    p = (x - v) * q - (x - w) * r;
 
117
                    q = 2 * (q - r);
 
118
 
 
119
                    if (q > 0) {
 
120
                        p = -p;
 
121
                    } else {
 
122
                        q = -q;
 
123
                    }
 
124
 
 
125
                    r = e;
 
126
                    e = d;
 
127
                }
 
128
 
 
129
                if (Math.abs(p) < Math.abs(0.5 * q * r) &&
 
130
                    (p < q * (a - x)) && (p < q * (b - x))) { // Parabolic interpolation step.
 
131
                    d = p / q;
 
132
                    u = x + d;
 
133
 
 
134
                    // f must not be evaluated too close to a or b.
 
135
                    if (((u - a) < t2) || ((b - u) < t2)) {
 
136
                        d = (x < m) ? tol : -tol;
 
137
                    }
 
138
                } else { // Golden section step.
 
139
                    e = ((x < m) ? b : a) - x;
 
140
                    d = c * e;
 
141
                }
 
142
 
 
143
                // f must not be evaluated too close to a or b.
 
144
                u = x + ((Math.abs(d) > tol) ? d : ((d > 0) ? tol : -tol));
 
145
                double fu = computeObjectiveValue(f, u);
 
146
                if (goalType == GoalType.MAXIMIZE) {
 
147
                    fu = -fu;
 
148
                }
 
149
 
 
150
                // Update a, b, v, w and x.
 
151
                if (fu <= fx) {
 
152
                    if (u < x) {
 
153
                        b = x;
 
154
                    } else {
 
155
                        a = x;
 
156
                    }
 
157
                    v = w;
 
158
                    fv = fw;
 
159
                    w = x;
 
160
                    fw = fx;
 
161
                    x = u;
 
162
                    fx = fu;
 
163
                } else {
 
164
                    if (u < x) {
 
165
                        a = u;
 
166
                    } else {
 
167
                        b = u;
 
168
                    }
 
169
                    if ((fu <= fw) || (w == x)) {
 
170
                        v = w;
 
171
                        fv = fw;
 
172
                        w = u;
 
173
                        fw = fu;
 
174
                    } else if ((fu <= fv) || (v == x) || (v == w)) {
 
175
                        v = u;
 
176
                        fv = fu;
 
177
                    }
 
178
                }
 
179
            } else { // termination
 
180
                setResult(x, (goalType == GoalType.MAXIMIZE) ? -fx : fx, count);
 
181
                return x;
 
182
            }
 
183
 
 
184
            ++count;
 
185
        }
 
186
 
 
187
        throw new MaxIterationsExceededException(maximalIterationCount);
 
188
 
 
189
    }
 
190
 
 
191
}