The R Project SVN R

Rev

Rev 24068 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 24068 Rev 27922
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  R : A Computer Language for Statistical Data Analysis
 2 
 *  R : A Computer Language for Statistical Data Analysis
3 
 
 3 
 
4 
 *  Copyright (C) 1997-2001   the R Development Core Team
 4 
 *  Copyright (C) 1997-2001  The R Development Core Team
- 
 
 5 
 *  Copyright (C) 2003-2004  The R Foundation
5 
 *
 6 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 7 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
 8 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
 9 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 10 
 *  (at your option) any later version.
Line 19... Line 20...
19 
 */
 20 
 */
20 
 
 21 
 
21 
/* fmin.f -- translated by f2c (version 19990503).
 22 
/* fmin.f -- translated by f2c (version 19990503).
22 
*/
 23 
*/
23 
 
 24 
 
24 
 
 - 
 
25 
/* R's  optimize() :   function	fmin(ax,bx,f,tol)
 25 
/* R's  optimize() :   function	fmin(ax,bx,f,tol)
26 
   =    ==========		~~~~~~~~~~~~~~~~~
 26 
   =    ==========		~~~~~~~~~~~~~~~~~
27 
 
 27 
 
28 
        an approximation  x  to the point where  f  attains a minimum  on
 28 
        an approximation  x  to the point where  f  attains a minimum  on
29 
    the interval  (ax,bx)  is determined.
 29 
    the interval  (ax,bx)  is determined.
30 
 
 30 
 
31 
    input..
 31 
    INPUT..
32 
 
 32 
 
33 
    ax    left endpoint of initial interval
 33 
    ax    left endpoint of initial interval
34 
    bx    right endpoint of initial interval
 34 
    bx    right endpoint of initial interval
35 
    f     function which evaluates  f(x, info)  for any  x
 35 
    f     function which evaluates  f(x, info)  for any  x
36 
          in the interval  (ax,bx)
 36 
          in the interval  (ax,bx)
37 
    tol   desired length of the interval of uncertainty of the final
 37 
    tol   desired length of the interval of uncertainty of the final
38 
          result (.ge.0.)
 38 
          result ( >= 0.)
39 
 
 39 
 
40 
    output..
 40 
    OUTPUT..
41 
 
 41 
 
42 
    fmin  abcissa approximating the point where  f  attains a
 42 
    fmin  abcissa approximating the point where  f  attains a minimum
43 
          minimum
 - 
 
44 
 
 43 
 
45 
        The method used is a combination of  golden  section  search  and
 44 
        The method used is a combination of  golden  section  search  and
46 
    successive parabolic interpolation.  convergence is never much slower
 45 
    successive parabolic interpolation.  convergence is never much slower
47 
    than  that  for  a  Fibonacci search.  If  f  has a continuous second
 46 
    than  that  for  a  Fibonacci search.  If  f  has a continuous second
48 
    derivative which is positive at the minimum (which is not  at  ax  or
 47 
    derivative which is positive at the minimum (which is not  at  ax  or
Line 76... Line 75...
76 
 
 75 
 
77 
    /* Local variables */
 76 
    /* Local variables */
78 
    double a, b, d, e, p, q, r, u, v, w, x;
 77 
    double a, b, d, e, p, q, r, u, v, w, x;
79 
    double t2, fu, fv, fw, fx, xm, eps, tol1, tol3;
 78 
    double t2, fu, fv, fw, fx, xm, eps, tol1, tol3;
80 
 
 79 
 
81 
/*  eps is approximately the square root of the relative machine
 80 
/*  eps is approximately the square root of the relative machine precision. */
82 
    precision. */
 - 
 
83 
    eps = d1mach(4);
 81 
    eps = d1mach(4);
84 
    tol1 = eps + 1.;
 82 
    tol1 = eps + 1.;/* the smallest 1.000... > 1 */
85 
    eps = sqrt(eps);
 83 
    eps = sqrt(eps);
86 
 
 84 
 
87 
    a = ax;
 85 
    a = ax;
88 
    b = bx;
 86 
    b = bx;
89 
    v = a + c * (b - a);
 87 
    v = a + c * (b - a);
90 
    w = v;
 88 
    w = v;
91 
    x = v;
 89 
    x = v;
92 
 
 90 
 
93 
/*  -Wall indicates that d may be used before being assigned */
 - 
 
94 
 
 - 
 
95 
    d = 0.;
 91 
    d = 0.;/* -Wall */
96 
    e = 0.;
 92 
    e = 0.;
97 
    fx = (*f)(x, info);
 93 
    fx = (*f)(x, info);
98 
    fv = fx;
 94 
    fv = fx;
99 
    fw = fx;
 95 
    fw = fx;
100 
    tol3 = tol / 3.;
 96 
    tol3 = tol / 3.;
Line 104... Line 100...
104 
    for(;;) {
 100 
    for(;;) {
105 
	xm = (a + b) * .5;
 101 
	xm = (a + b) * .5;
106 
	tol1 = eps * fabs(x) + tol3;
 102 
	tol1 = eps * fabs(x) + tol3;
107 
	t2 = tol1 * 2.;
 103 
	t2 = tol1 * 2.;
108 
 
 104 
 
109 
/*  check stopping criterion */
 105 
	/* check stopping criterion */
110 
 
 106 
 
111 
	if (fabs(x - xm) <= t2 - (b - a) * .5) break;
 107 
	if (fabs(x - xm) <= t2 - (b - a) * .5) break;
112 
	p = 0.;
 108 
	p = 0.;
113 
	q = 0.;
 109 
	q = 0.;
114 
	r = 0.;
 110 
	r = 0.;
115 
	if (fabs(e) > tol1) {
 111 
	if (fabs(e) > tol1) { /* fit parabola */
116 
 
 - 
 
117 
/*     fit parabola */
 - 
 
118 
 
 112 
 
119 
	    r = (x - w) * (fx - fv);
 113 
	    r = (x - w) * (fx - fv);
120 
	    q = (x - v) * (fx - fw);
 114 
	    q = (x - v) * (fx - fw);
121 
	    p = (x - v) * q - (x - w) * r;
 115 
	    p = (x - v) * q - (x - w) * r;
122 
	    q = (q - r) * 2.;
 116 
	    q = (q - r) * 2.;
123 
	    if (q > 0.) p = -p; else q = -q;
 117 
	    if (q > 0.) p = -p; else q = -q;
124 
	    r = e;
 118 
	    r = e;
125 
	    e = d;
 119 
	    e = d;
126 
	}
 120 
	}
127 
	if (fabs(p) >= fabs(q * .5 * r) ||
 - 
 
128 
	    p <= q * (a - x) || p >= q * (b - x)) {
 - 
 
129 
 
 121 
 
- 
 
 122 
	if (fabs(p) >= fabs(q * .5 * r) ||
130 
/*     a golden-section step */
 123 
	    p <= q * (a - x) || p >= q * (b - x)) { /* a golden-section step */
131 
 
 124 
 
132 
	    if (x < xm) e = b - x; else e = a - x;
 125 
	    if (x < xm) e = b - x; else e = a - x;
133 
	    d = c * e;
 126 
	    d = c * e;
134 
	} else {
 - 
 
135 
 
 127 
	}
136 
/*     a parabolic-interpolation step */
 128 
	else { /* a parabolic-interpolation step */
137 
 
 129 
 
138 
	    d = p / q;
 130 
	    d = p / q;
139 
	    u = x + d;
 131 
	    u = x + d;
140 
 
 132 
 
141 
/*     f must not be evaluated too close to ax or bx */
 133 
	    /* f must not be evaluated too close to ax or bx */
142 
 
 134 
 
143 
	    if (u - a >= t2 && b - u >= t2) goto L90;
 135 
	    if (u - a < t2 || b - u < t2) {
144 
	    d = tol1;
 136 
		d = tol1;
145 
	    if (x >= xm) d = -d;
 137 
		if (x >= xm) d = -d;
- 
 
 138 
	    }
146 
	}
 139 
	}
147 
 
 140 
 
148 
/* f must not be evaluated too close to x */
 141 
	/* f must not be evaluated too close to x */
149 
 
 142 
 
150 
    L90:
 - 
 
151 
	if (fabs(d) >= tol1) {
 143 
	if (fabs(d) >= tol1)
152 
	    u = x + d;
 144 
	    u = x + d;
153 
	} else {
 - 
 
154 
	    if (d > 0.) {
 145 
	else if (d > 0.)
155 
		u = x + tol1;
 146 
	    u = x + tol1;
156 
	    } else {
 147 
	else
157 
		u = x - tol1;
 148 
	    u = x - tol1;
158 
	    }
 - 
 
159 
	}
 149 
 
160 
	fu = (*f)(u, info);
 150 
	fu = (*f)(u, info);
161 
 
 151 
 
162 
/*  update  a, b, v, w, and x */
 152 
	/*  update  a, b, v, w, and x */
163 
 
 153 
 
164 
	if (fx <= fu) {
 154 
	if (fx <= fu) {
165 
	    if (u < x) {
 155 
	    if (u < x)
166 
		a = u;
 156 
		a = u;
167 
	    } else {
 157 
	    else
168 
		b = u;
 158 
		b = u;
169 
	    }
 - 
 
170 
	}
 159 
	}
171 
	if (fu <= fx) {
 160 
	if (fu <= fx) {
172 
	    if (u < x) {
 161 
	    if (u < x)
173 
		b = x;
 162 
		b = x;
174 
	    } else {
 163 
	    else
175 
		a = x;
 164 
		a = x;
176 
	    }
 - 
 
177 
	    v = w;
 - 
 
178 
	    fv = fw;
 - 
 
179 
	    w = x;
 165 
	    v = w;    w = x;   x = u;
180 
	    fw = fx;
 166 
	    fv = fw; fw = fx; fx = fu;
181 
	    x = u;
 167 
	}
182 
	    fx = fu;
 - 
 
183 
	} else {
 168 
	else {
184 
	    if (fu <= fw || w == x) {
 169 
	    if (fu <= fw || w == x) {
185 
		v = w;
 - 
 
186 
		fv = fw;
 170 
		v = w; fv = fw;
187 
		w = u;
 - 
 
188 
		fw = fu;
 171 
		w = u; fw = fu;
189 
	    } else {
 172 
	    }
190 
		if (fu <= fv || v == x || v == w) {
 173 
	    else if (fu <= fv || v == x || v == w) {
191 
		    v = u;
 - 
 
192 
		    fv = fu;
 174 
		v = u; fv = fu;
193 
		}
 - 
 
194 
	    }
 175 
	    }
195 
	}
 176 
	}
196 
    }
 177 
    }
197 
/*  end of main loop */
 178 
    /* end of main loop */
198 
 
 179 
 
199 
    return x;
 180 
    return x;
200 
}
 181 
}
201 
 
 -