bcftools/kmin.c

   1 /* The MIT License
   2
   3    Copyright (c) 2008, 2010 by Attractive Chaos <attractor@live.co.uk>
   4
   5    Permission is hereby granted, free of charge, to any person obtaining
   6    a copy of this software and associated documentation files (the
   7    "Software"), to deal in the Software without restriction, including
   8    without limitation the rights to use, copy, modify, merge, publish,
   9    distribute, sublicense, and/or sell copies of the Software, and to
  10    permit persons to whom the Software is furnished to do so, subject to
  11    the following conditions:
  12
  13    The above copyright notice and this permission notice shall be
  14    included in all copies or substantial portions of the Software.
  15
  16    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  18    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  20    BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  21    ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  22    CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  23    SOFTWARE.
  24 */
  25
  26 /* Hooke-Jeeves algorithm for nonlinear minimization
  27
  28    Based on the pseudocodes by Bell and Pike (CACM 9(9):684-685), and
  29    the revision by Tomlin and Smith (CACM 12(11):637-638). Both of the
  30    papers are comments on Kaupe's Algorithm 178 "Direct Search" (ACM
  31    6(6):313-314). The original algorithm was designed by Hooke and
  32    Jeeves (ACM 8:212-229). This program is further revised according to
  33    Johnson's implementation at Netlib (opt/hooke.c).
  34
  35    Hooke-Jeeves algorithm is very simple and it works quite well on a
  36    few examples. However, it might fail to converge due to its heuristic
  37    nature. A possible improvement, as is suggested by Johnson, may be to
  38    choose a small r at the beginning to quickly approach to the minimum
  39    and a large r at later step to hit the minimum.
  40  */
  41
  42 #include <stdlib.h>
  43 #include <string.h>
  44 #include <math.h>
  45 #include "kmin.h"
  46
  47 static double __kmin_hj_aux(kmin_f func, int n, double *x1, void *data, double fx1, double *dx, int *n_calls)
  48 {
  49         int k, j = *n_calls;
  50         double ftmp;
  51         for (k = 0; k != n; ++k) {
  52                 x1[k] += dx[k];
  53                 ftmp = func(n, x1, data); ++j;
  54                 if (ftmp < fx1) fx1 = ftmp;
  55                 else { /* search the opposite direction */
  56                         dx[k] = 0.0 - dx[k];
  57                         x1[k] += dx[k] + dx[k];
  58                         ftmp = func(n, x1, data); ++j;
  59                         if (ftmp < fx1) fx1 = ftmp;
  60                         else x1[k] -= dx[k]; /* back to the original x[k] */
  61                 }
  62         }
  63         *n_calls = j;
  64         return fx1; /* here: fx1=f(n,x1) */
  65 }
  66
  67 double kmin_hj(kmin_f func, int n, double *x, void *data, double r, double eps, int max_calls)
  68 {
  69         double fx, fx1, *x1, *dx, radius;
  70         int k, n_calls = 0;
  71         x1 = (double*)calloc(n, sizeof(double));
  72         dx = (double*)calloc(n, sizeof(double));
  73         for (k = 0; k != n; ++k) { /* initial directions, based on MGJ */
  74                 dx[k] = fabs(x[k]) * r;
  75                 if (dx[k] == 0) dx[k] = r;
  76         }
  77         radius = r;
  78         fx1 = fx = func(n, x, data); ++n_calls;
  79         for (;;) {
  80                 memcpy(x1, x, n * sizeof(double)); /* x1 = x */
  81                 fx1 = __kmin_hj_aux(func, n, x1, data, fx, dx, &n_calls);
  82                 while (fx1 < fx) {
  83                         for (k = 0; k != n; ++k) {
  84                                 double t = x[k];
  85                                 dx[k] = x1[k] > x[k]? fabs(dx[k]) : 0.0 - fabs(dx[k]);
  86                                 x[k] = x1[k];
  87                                 x1[k] = x1[k] + x1[k] - t;
  88                         }
  89                         fx = fx1;
  90                         if (n_calls >= max_calls) break;
  91                         fx1 = func(n, x1, data); ++n_calls;
  92                         fx1 = __kmin_hj_aux(func, n, x1, data, fx1, dx, &n_calls);
  93                         if (fx1 >= fx) break;
  94                         for (k = 0; k != n; ++k)
  95                                 if (fabs(x1[k] - x[k]) > .5 * fabs(dx[k])) break;
  96                         if (k == n) break;
  97                 }
  98                 if (radius >= eps) {
  99                         if (n_calls >= max_calls) break;
 100                         radius *= r;
 101                         for (k = 0; k != n; ++k) dx[k] *= r;
 102                 } else break; /* converge */
 103         }
 104         free(x1); free(dx);
 105         return fx1;
 106 }
 107
 108 // I copied this function somewhere several years ago with some of my modifications, but I forgot the source.
 109 double kmin_brent(kmin1_f func, double a, double b, void *data, double tol, double *xmin)
 110 {
 111         double bound, u, r, q, fu, tmp, fa, fb, fc, c;
 112         const double gold1 = 1.6180339887;
 113         const double gold2 = 0.3819660113;
 114         const double tiny = 1e-20;
 115         const int max_iter = 100;
 116
 117         double e, d, w, v, mid, tol1, tol2, p, eold, fv, fw;
 118         int iter;
 119
 120         fa = func(a, data); fb = func(b, data);
 121         if (fb > fa) { // swap, such that f(a) > f(b)
 122                 tmp = a; a = b; b = tmp;
 123                 tmp = fa; fa = fb; fb = tmp;
 124         }
 125         c = b + gold1 * (b - a), fc = func(c, data); // golden section extrapolation
 126         while (fb > fc) {
 127                 bound = b + 100.0 * (c - b); // the farthest point where we want to go
 128                 r = (b - a) * (fb - fc);
 129                 q = (b - c) * (fb - fa);
 130                 if (fabs(q - r) < tiny) { // avoid 0 denominator
 131                         tmp = q > r? tiny : 0.0 - tiny;
 132                 } else tmp = q - r;
 133                 u = b - ((b - c) * q - (b - a) * r) / (2.0 * tmp); // u is the parabolic extrapolation point
 134                 if ((b > u && u > c) || (b < u && u < c)) { // u lies between b and c
 135                         fu = func(u, data);
 136                         if (fu < fc) { // (b,u,c) bracket the minimum
 137                                 a = b; b = u; fa = fb; fb = fu;
 138                                 break;
 139                         } else if (fu > fb) { // (a,b,u) bracket the minimum
 140                                 c = u; fc = fu;
 141                                 break;
 142                         }
 143                         u = c + gold1 * (c - b); fu = func(u, data); // golden section extrapolation
 144                 } else if ((c > u && u > bound) || (c < u && u < bound)) { // u lies between c and bound
 145                         fu = func(u, data);
 146                         if (fu < fc) { // fb > fc > fu
 147                                 b = c; c = u; u = c + gold1 * (c - b);
 148                                 fb = fc; fc = fu; fu = func(u, data);
 149                         } else { // (b,c,u) bracket the minimum
 150                                 a = b; b = c; c = u;
 151                                 fa = fb; fb = fc; fc = fu;
 152                                 break;
 153                         }
 154                 } else if ((u > bound && bound > c) || (u < bound && bound < c)) { // u goes beyond the bound
 155                         u = bound; fu = func(u, data);
 156                 } else { // u goes the other way around, use golden section extrapolation
 157                         u = c + gold1 * (c - b); fu = func(u, data);
 158                 }
 159                 a = b; b = c; c = u;
 160                 fa = fb; fb = fc; fc = fu;
 161         }
 162         if (a > c) u = a, a = c, c = u; // swap
 163
 164         // now, a<b<c, fa>fb and fb<fc, move on to Brent's algorithm
 165         e = d = 0.0;
 166         w = v = b; fv = fw = fb;
 167         for (iter = 0; iter != max_iter; ++iter) {
 168                 mid = 0.5 * (a + c);
 169                 tol2 = 2.0 * (tol1 = tol * fabs(b) + tiny);
 170                 if (fabs(b - mid) <= (tol2 - 0.5 * (c - a))) {
 171                         *xmin = b; return fb; // found
 172                 }
 173                 if (fabs(e) > tol1) {
 174                         // related to parabolic interpolation
 175                         r = (b - w) * (fb - fv);
 176                         q = (b - v) * (fb - fw);
 177                         p = (b - v) * q - (b - w) * r;
 178                         q = 2.0 * (q - r);
 179                         if (q > 0.0) p = 0.0 - p;
 180                         else q = 0.0 - q;
 181                         eold = e; e = d;
 182                         if (fabs(p) >= fabs(0.5 * q * eold) || p <= q * (a - b) || p >= q * (c - b)) {
 183                                 d = gold2 * (e = (b >= mid ? a - b : c - b));
 184                         } else {
 185                                 d = p / q; u = b + d; // actual parabolic interpolation happens here
 186                                 if (u - a < tol2 || c - u < tol2)
 187                                         d = (mid > b)? tol1 : 0.0 - tol1;
 188                         }
 189                 } else d = gold2 * (e = (b >= mid ? a - b : c - b)); // golden section interpolation
 190                 u = fabs(d) >= tol1 ? b + d : b + (d > 0.0? tol1 : -tol1);
 191                 fu = func(u, data);
 192                 if (fu <= fb) { // u is the minimum point so far
 193                         if (u >= b) a = b;
 194                         else c = b;
 195                         v = w; w = b; b = u; fv = fw; fw = fb; fb = fu;
 196                 } else { // adjust (a,c) and (u,v,w)
 197                         if (u < b) a = u;
 198                         else c = u;
 199                         if (fu <= fw || w == b) {
 200                                 v = w; w = u;
 201                                 fv = fw; fw = fu;
 202                         } else if (fu <= fv || v == b || v == w) {
 203                                 v = u; fv = fu;
 204                         }
 205                 }
 206         }
 207         *xmin = b;
 208         return fb;
 209 }