From: Attractive Chaos Date: Thu, 13 Jan 2011 17:42:40 +0000 (-0500) Subject: Added the kmin library X-Git-Tag: ksprintf-final~75 X-Git-Url: http://git.kaiwu.me/postgresql/log/contrib/postgres_fdw/postgres_fdw.c?a=commitdiff_plain;h=4493780d0aac84db60176d0d8f1779127ee92a87;p=klib.git Added the kmin library --- diff --git a/kmin.c b/kmin.c new file mode 100644 index 0000000..8cb0acf --- /dev/null +++ b/kmin.c @@ -0,0 +1,107 @@ +/* The MIT License + + Copyright (c) 2008, by Heng Li + + 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 +#include +#include +#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 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 index 0000000..cb1a4e8 --- /dev/null +++ b/kmin_test.c @@ -0,0 +1,48 @@ +#include +#include +#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; +}