Added the kmin library

diff --git a/kmin.c b/kmin.c
new file mode 100644 (file)
index 0000000..8cb0acf
--- /dev/null
+++ b/kmin.c
@@ -0,0 +1,107 @@
+/* The MIT License
+
+   Copyright (c) 2008, by Heng Li <lh3@live.co.uk>
+
+   Permission is hereby granted, free of charge, to any person obtaining
+   a copy of this software and associated documentation files (the
+   "Software"), to deal in the Software without restriction, including
+   without limitation the rights to use, copy, modify, merge, publish,
+   distribute, sublicense, and/or sell copies of the Software, and to
+   permit persons to whom the Software is furnished to do so, subject to
+   the following conditions:
+
+   The above copyright notice and this permission notice shall be
+   included in all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+   SOFTWARE.
+*/
+
+/* Hooke-Jeeves algorithm for nonlinear minimization
+   Heng Li, Februay 3, 2008
+ 
+   Based on the pseudocodes by Bell and Pike (CACM 9(9):684-685), and
+   the revision by Tomlin and Smith (CACM 12(11):637-638). Both of the
+   papers are comments on Kaupe's Algorithm 178 "Direct Search" (ACM
+   6(6):313-314). The original algorithm was designed by Hooke and
+   Jeeves (ACM 8:212-229). This program is further revised according to
+   Johnson's implementation at Netlib (opt/hooke.c).
+ 
+   Hooke-Jeeves algorithm is very simple and it works quite well on a
+   few examples. However, it might fail to converge due to its heuristic
+   nature. A possible improvement, as is suggested by Johnson, may be to
+   choose a small r at the beginning to quickly approach to the minimum
+   and a large r at later step to hit the minimum.
+ */
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "kmin.h"
+
+static double __kmin_hj_aux(kmin_f func, int n, double *x1, void *data, double fx1, double *dx, int *n_calls)
+{
+       int k, j = *n_calls;
+       double ftmp;
+       for (k = 0; k != n; ++k) {
+               x1[k] += dx[k];
+               ftmp = func(n, x1, data); ++j;
+               if (ftmp < fx1) fx1 = ftmp;
+               else { /* search the opposite direction */
+                       dx[k] = 0.0 - dx[k];
+                       x1[k] += dx[k] + dx[k];
+                       ftmp = func(n, x1, data); ++j;
+                       if (ftmp < fx1) fx1 = ftmp;
+                       else x1[k] -= dx[k]; /* back to the original x[k] */
+               }
+       }
+       *n_calls = j;
+       return fx1; /* here: fx1=f(n,x1) */
+}
+
+double kmin_hj(kmin_f func, int n, double *x, void *data, double r, double eps, int max_calls)
+{
+       double fx, fx1, *x1, *dx, radius;
+       int k, n_calls = 0;
+       x1 = (double*)calloc(n, sizeof(double));
+       dx = (double*)calloc(n, sizeof(double));
+       for (k = 0; k != n; ++k) { /* initial directions, based on MGJ */
+               dx[k] = fabs(x[k]) * r;
+               if (dx[k] == 0) dx[k] = r;
+       }
+       radius = r;
+       fx1 = fx = func(n, x, data); ++n_calls;
+       for (;;) {
+               memcpy(x1, x, n * sizeof(double)); /* x1 = x */
+               fx1 = __kmin_hj_aux(func, n, x1, data, fx, dx, &n_calls);
+               while (fx1 < fx) {
+                       for (k = 0; k != n; ++k) {
+                               double t = x[k];
+                               dx[k] = x1[k] > x[k]? fabs(dx[k]) : 0.0 - fabs(dx[k]);
+                               x[k] = x1[k];
+                               x1[k] = x1[k] + x1[k] - t;
+                       }
+                       fx = fx1;
+                       if (n_calls >= max_calls) break;
+                       fx1 = func(n, x1, data); ++n_calls;
+                       fx1 = __kmin_hj_aux(func, n, x1, data, fx1, dx, &n_calls);
+                       if (fx1 >= fx) break;
+                       for (k = 0; k != n; ++k)
+                               if (fabs(x1[k] - x[k]) > .5 * fabs(dx[k])) break;
+                       if (k == n) break;
+               }
+               if (radius >= eps) {
+                       if (n_calls >= max_calls) break;
+                       radius *= r;
+                       for (k = 0; k != n; ++k) dx[k] *= r;
+               } else break; /* converge */
+       }
+       free(x1); free(dx);
+       return fx1;
+}

diff --git a/kmin.h b/kmin.h
new file mode 100644 (file)
index 0000000..4298b3f
--- /dev/null
+++ b/kmin.h
@@ -0,0 +1,20 @@
+#ifndef KMIN_H
+#define KMIN_H
+
+#define KMIN_RADIUS  0.5
+#define KMIN_EPS     1e-7
+#define KMIN_MAXCALL 50000
+
+typedef double (*kmin_f)(int, double*, void*);
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+       double kmin_hj(kmin_f func, int n, double *x, void *data, double r, double eps, int max_calls);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

diff --git a/kmin_test.c b/kmin_test.c
new file mode 100644 (file)
index 0000000..cb1a4e8
--- /dev/null
+++ b/kmin_test.c
@@ -0,0 +1,48 @@
+#include <stdio.h>
+#include <math.h>
+#include "kmin.h"
+
+static int n_evals;
+
+double f_Chebyquad(int n, double *x, void *data)
+{
+    int i, j;
+    double y[20][20], f;
+    int np, iw;
+    double sum;
+    for (j = 0; j != n; ++j) {
+               y[0][j] = 1.;
+               y[1][j] = 2. * x[j] - 1.;
+    }
+    for (i = 1; i != n; ++i)
+               for (j = 0; j != n; ++j)
+                       y[i+1][j] = 2. * y[1][j] * y[i][j] - y[i-1][j];
+    f = 0.;
+    np = n + 1;
+    iw = 1;
+    for (i = 0; i != np; ++i) {
+               sum = 0.;
+               for (j = 0; j != n; ++j) sum += y[i][j];
+               sum /= n;
+               if (iw > 0) sum += 1. / ((i - 1) * (i + 1));
+               iw = -iw;
+               f += sum * sum;
+    }
+       ++n_evals;
+    return f;
+}
+
+int main()
+{
+       double x[20], y;
+       int n, i;
+       printf("\nMinimizer: Hooke-Jeeves\n");
+       for (n = 2; n <= 8; n += 2) {
+               for (i = 0; i != n; ++i) x[i] = (double)(i + 1) / n;
+               n_evals = 0;
+               y = kmin_hj(f_Chebyquad, n, x, 0, KMIN_RADIUS, KMIN_EPS, KMIN_MAXCALL);
+               printf("n=%d,min=%.8lg,n_evals=%d\n", n, y, n_evals);
+       }
+       printf("\n");
+       return 0;
+}