summaryrefslogtreecommitdiff
path: root/numlib/dhsx.c
diff options
context:
space:
mode:
authorJörg Frings-Fürst <debian@jff-webhosting.net>2014-09-01 13:56:46 +0200
committerJörg Frings-Fürst <debian@jff-webhosting.net>2014-09-01 13:56:46 +0200
commit22f703cab05b7cd368f4de9e03991b7664dc5022 (patch)
tree6f4d50beaa42328e24b1c6b56b6ec059e4ef21a5 /numlib/dhsx.c
Initial import of argyll version 1.5.1-8debian/1.5.1-8
Diffstat (limited to 'numlib/dhsx.c')
-rw-r--r--numlib/dhsx.c302
1 files changed, 302 insertions, 0 deletions
diff --git a/numlib/dhsx.c b/numlib/dhsx.c
new file mode 100644
index 0000000..180d15c
--- /dev/null
+++ b/numlib/dhsx.c
@@ -0,0 +1,302 @@
+
+/* This is better for stocastic optimisation, where the function */
+/* being evaluated may have a random component, or is not smooth. */
+
+/*
+ * Copyright 1999 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ */
+
+/* A general purpose downhill simplex multivariate optimser, */
+/* based on the Nelder and Mead algorithm. */
+/* Code is an original expression of the algorithms decsribed in */
+/* "Numerical Recipes in C", by W.H.Press, B.P.Flannery, */
+/* S.A.Teukolsky & W.T.Vetterling. */
+
+#include "numsup.h"
+
+#undef DEBUG
+
+static void simplexinit(int di, double *cp, double *r,double **p);
+static double trypoint(int di,double *cp, double **p, double *y, int hix, double hpfac,
+ double (*funk)(void *fdata, double *tp), void *fdata, double *tryp);
+
+#ifdef NEVER /* Experimental */
+
+#define ALPHA 0.7 /* Extrapolate hight point through oposite face factor */
+#define GAMMA 1.4 /* Aditional extrapolation if ALPHA is good */
+#define BETA 0.4 /* One dimensional contraction factor */
+#define NONEXP 2 /* non expanding passes */
+
+#else /* Standard tuning values */
+
+#define ALPHA 1.0 /* Extrapolate hight point through oposite face factor */
+#define GAMMA 2.0 /* Aditional extrapolation if ALPHA is good */
+#define BETA 0.5 /* One dimensional contraction factor */
+#define NONEXP 2 /* non expanding passes */
+
+#endif
+
+
+/* Down hill simplex function */
+/* return 0 on sucess, 1 on failure due to excessive itterations */
+/* Result will be in cp */
+/* Arrays start at 0 */
+int dhsx(
+double *rv, /* If not NULL, return the residual error */
+int di, /* Dimentionality */
+double *cp, /* Initial starting point */
+double *s, /* Size of initial search area */
+double ftol, /* Finishing tollerance of error change */
+double atol, /* Absolute return value tollerance */
+int maxit, /* Maximum iterations allowed */
+double (*funk)(void *fdata, double *tp), /* Error function to evaluate */
+void *fdata /* Data needed by function */
+) {
+ int i, j;
+ int lox, hix, nhix; /* Lowest point index, highest point, next highest point */
+ int nsp = di+1; /* Number of simplex verticy points */
+ int nit; /* Number of iterations */
+ double tryy, ysave;
+ double tol;
+ double **p; /* Simplex array */
+ double *y; /* values of func at verticies */
+ double *tryp; /* Temporary used by trypoint() */
+
+ /* Allocate array arrays */
+ y = dvector(0, di); /* Value of function at verticies */
+ tryp = dvector(0, di-1);
+ p = dmatrix(0, di+1, 0, di); /* Vertex array of dimentions */
+
+ /* Init the search simplex */
+ simplexinit(di, cp, s, p);
+
+ /* Compute current center point position */
+ for (j = 0; j < di; j++) { /* For all dimensions */
+ double sum;
+ for (i = 0, sum = 0.0; i < nsp; i++) /* For all verticies */
+ sum += p[i][j];
+ cp[j] = sum;
+ }
+
+ /* Compute initial y (function) values at verticies */
+ for (i = 0; i < nsp; i++) /* For all verticies */
+ y[i] = (*funk)(fdata,p[i]); /* Compute error function */
+
+ /* Untill we find a solution or give up */
+ for (nit = 0; nit < maxit; nit++) {
+ /* Find highest, next highest and lowest vertex */
+
+ lox = nhix = hix = 0;
+ for (i = 0; i < nsp; i++) {
+ if (y[i] < y[lox])
+ lox = i;
+ if (y[i] > y[hix]) {
+ nhix = hix;
+ hix = i;
+ } else if (y[i] > y[nhix]) {
+ nhix = i;
+ }
+ }
+
+ tol = y[hix] - y[lox];
+
+#ifdef DEBUG /* 2D */
+ printf("Current vs = %f,%f %f,%f %f,%f\n",
+ p[0].c[0],p[0].c[1],p[1].c[0],p[1].c[1],p[2].c[0],p[2].c[1]);
+ printf("Current errs = %e %e %e\n",y[0],y[1],y[2]);
+ printf("Current sxs = %d %d %d\n",sy[0]->no,sy[1]->no,sy[2]->no);
+ printf("Current y[hix] = %e\n",y[hix]);
+#endif /* DEBUG */
+
+ if (tol < ftol && y[lox] < atol) { /* Found an adequate solution */
+ /* (allow 10 x range for disambiguator) */
+ /* set cp[] to best point, and return error value of that point */
+ tryy = 1.0/(di+1);
+ for (j = 0; j < di; j++) /* For all dimensions */
+ cp[j] *= tryy; /* Set cp to center point of simplex */
+#ifdef DEBUG
+ printf("C point = %f,%f\n",cp[0],cp[1]);
+#endif
+ tryy = (*funk)(fdata,cp); /* Compute error function */
+
+ if (tryy > y[lox]) { /* Center point is not the best */
+ tryy = y[lox];
+ for (j = 0; j < di; j++)
+ cp[j] = p[lox][j];
+ }
+ free_dmatrix(p, 0, di+1, 0, di);
+ free_dvector(tryp, 0, di-1);
+ free_dvector(y, 0, di);
+printf("~1 itterations = %d\n",nit);
+ if (rv != NULL)
+ *rv = tryy;
+ return 0;
+ }
+
+ if (nit > NONEXP) { /* Only try expanding after a couple of iterations */
+ /* Try moving the high point through the oposite face by ALPHA */
+#ifdef DEBUG
+ printf("dhsx: try moving high point %d through oposite face",hix);
+#endif
+ tryy = trypoint(di, cp, p, y, hix, -ALPHA, funk, fdata, tryp);
+ }
+
+ if (nit > NONEXP && tryy <= y[lox]) {
+#ifdef DEBUG
+ verbose(4,"dhsx: moving high through oposite face worked");
+#endif
+ /* Gave good result, so continue on in that direction */
+ tryy=trypoint(di,cp,p,y,hix,GAMMA,funk,fdata,tryp);
+
+
+ } else if (nit <= NONEXP || tryy >= y[nhix]) {
+
+ /* else if ALPHA move made things worse, do a one dimensional */
+ /* contraction by a factor BETA */
+#ifdef DEBUG
+ verbose(4,"dhsx: else try moving contracting point %d, y[ini] = %f",hix,y[hix]);
+#endif
+ ysave = y[hix];
+ tryy = trypoint(di, cp, p, y, hix, BETA, funk, fdata, tryp);
+
+ if (tryy >= ysave) {
+#ifdef DEBUG
+ verbose(4,"dhsx: contracting didn't work, try contracting other points to low");
+#endif
+ /* That still didn't help us, so move all the */
+ /* other points towards the high point */
+ for (i = 0; i < nsp; i++) { /* For all verts except low */
+ if (i != lox) {
+ for (j = 0; j < di; j++) { /* For all dimensions */
+ cp[j] = 0.5 * (p[i][j] + p[lox][j]);
+ p[i][j] = cp[j];
+ }
+ /* Compute function value for new point */
+ y[i] = (*funk)(fdata,p[i]); /* Compute error function */
+ }
+ }
+ /* Re-compute current center point value */
+ for (j = 0; j < di; j++) {
+ double sum;
+ for (i = 0,sum = 0.0;i<nsp;i++)
+ sum += p[i][j];
+ cp[j] = sum;
+ }
+ } else {
+#ifdef DEBUG
+ verbose(4,"dhsx: contracting point %d worked, tryy = %e, ysave = %e",hix,tryy,ysave);
+#endif
+ }
+ }
+ }
+ free_dmatrix(p, 0, di+1, 0, di);
+ free_dvector(tryp, 0, di-1);
+ free_dvector(y, 0, di);
+ if (rv != NULL)
+ *rv = tryy;
+ return 1; /* Failed */
+}
+
+/* Try moving the high point through the opposite face */
+/* by a factor of fac, and replaces the high point if */
+/* that proves to be better. */
+static double trypoint(
+int di, /* Dimentionality */
+double *cp, /* nsp * center coord/Returned coordinate */
+double **p, /* Starting/Current simplex (modified by dhsx) */
+double *y, /* values of func at verticies */
+int hix, /* Index of high point we are moving */
+double hpfac, /* factor to move high point */
+double (*funk)(void *fdata, double tp[]), /* Error function to evaluate */
+void *fdata, /* Data needed by function */
+double *tryp /* temporary array of size di-1 */
+) {
+ int j;
+ double tt, tryy;
+
+ /* Compute trial high point */
+ tt = (1.0 - hpfac)/di;
+ for (j = 0; j < di; j++)
+ tryp[j] = cp[j] * tt - p[hix][j] * (tt - hpfac);
+
+ /* Evaluate trial point */
+ tryy = (*funk)(fdata, tryp); /* Compute error function */
+
+ /* If new high point pos. is better */
+ if (tryy < y[hix]) {
+#ifdef DEBUG
+ printf("Try gave improved %e from sx %d",tryy,lsxp->last_fwd->no);
+#endif
+ y[hix] = tryy; /* Replace func val of hi with trial */
+
+ for (j = 0; j < di; j++) {
+ cp[j] += tryp[j] - p[hix][j]; /* Recompute cp */
+ p[hix][j] = tryp[j]; /* Replace co-ords of hi with trial */
+ }
+ } else {
+#ifdef DEBUG
+ verbose(4,"Try gave worse %e from sx %d",tryy,
+ lsxp->last_fwd == NULL ? -999 : lsxp->last_fwd->no);
+#endif
+ }
+ return tryy; /* Function value of trial point */
+}
+
+
+/* Make up an initial simplex for dhsx routine */
+static void
+simplexinit(
+int di, /* Dimentionality */
+double *cp, /* Initial solution values */
+double *s, /* initial radius for each dimention */
+double **p /* Simplex to initialize */
+) {
+ double bb;
+ double hh = 0.5; /* Constant */
+ double rr = sqrt(3.0)/2.0; /* Constant */
+ int i,j;
+ for (i = 0; i <= di; i++) { /* For each vertex */
+ /* The bounding points form a equalateral simplex */
+ /* whose vertexes are on a spehere about the data */
+ /* point center. The coordinate sequence is: */
+ /* B = sphere radius */
+ /* H = 0.5 */
+ /* R = sqrt(3)/2 */
+ /* 0 0 0 +B */
+ /* 0 0 0 -B */
+
+ /* 0 0 0 +B */
+ /* 0 0 +RB -HB */
+ /* 0 0 -RB -HB */
+
+ /* 0 0 0 +B */
+ /* 0 0 +RB -HB */
+ /* 0 +RRb -HRB -HB */
+ /* 0 -RRb -HRB -HB */
+
+ /* 0 0 0 +B */
+ /* 0 0 +RB -HB */
+ /* 0 +RRb -HRB -HB */
+ /* +RRRb -HRRb -HRB -HB */
+ /* -RRRb -HRRb -HRB -HB */
+
+ /* etc. */
+ bb = 1.0; /* Initial unscaled radius */
+ for (j = 0; j < di; j++) { /* For each coordinate in vertex */
+ if (j > i)
+ p[i][j] = cp[j] + s[j] * 0.0; /* If beyond last */
+ else if (j == i) /* If last non-zero */
+ p[i][j] = cp[j] + s[j] * bb;
+ else if (i == di && j == (di-1)) /* If last of last */
+ p[i][j] = cp[j] + s[j] * -1.0 * bb;
+ else /* If before last */
+ p[i][j] = cp[j] + s[j] * -hh * bb;
+ bb *= rr;
+ }
+ }
+}
+