OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
amoeba.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdint.h>
3 #include <math.h>
4 #include "amoeba.h"
5 
6 #define GET_PSUM \
7 for (j=0;j<ndim;j++) \
8 { \
9  for (sum=0.0,i=0;i<ndim+1;i++) \
10  sum += *(pnts + i*ndim +j); \
11  psum[j] = sum;\
12 }
13 
14 short amoeba(double *pnts, FITSTRUCT * auxdata,
15  double (*func)(FITSTRUCT *, double []), float tol) {
16  short i, j, k, isml = 0, inxt, ibig, ndim;
17  double y[16], ptry[16], psum[16], ytry, ysave, sum, rtol;
18  double reflect(double *, double [], double [], short, FITSTRUCT *,
19  double (*func)(FITSTRUCT *, double []), short, float);
20  int32_t itermax = auxdata->niter;
21 
22 
23  ndim = auxdata->nfunc;
24 
25 
26  GET_PSUM
27 
28  /* Determine initial functional values */
29  for (i = 0; i < ndim + 1; i++) {
30  for (j = 0; j < ndim; j++) ptry[j] = pnts[i * ndim + j];
31  y[i] = (*func)(auxdata, ptry);
32  }
33 
34  for (k = 0; k < itermax; k++) {
35 
36  auxdata->niter = k + 1;
37 
38  /* Order by smallest, second-largest, largest */
39  if (y[0] > y[1]) ibig = 0;
40  else ibig = 1;
41  isml = 1 - ibig;
42  inxt = isml;
43 
44  for (i = 2; i < ndim + 1; i++) {
45  if (y[i] > y[ibig]) {
46  inxt = ibig;
47  ibig = i;
48  } else if (y[i] < y[isml]) {
49  isml = i;
50  } else if (y[i] > y[inxt] && y[i] < y[ibig]) inxt = i;
51  }
52 
53  rtol = 2.0 * fabs(y[ibig] - y[isml]) /
54  (fabs(y[ibig]) + fabs(y[isml]));
55 
56  /* printf("%f %f\n", rtol,tol);*/
57 
58  if (rtol < tol) break;
59 
60  ytry = reflect(pnts, y, psum, ndim, auxdata, func, ibig, -1.0);
61 
62  if (ytry <= y[isml]) {
63  ytry = reflect(pnts, y, psum, ndim, auxdata, func, ibig, 2.0);
64  }
65  else if (ytry >= y[inxt]) {
66  ysave = y[ibig];
67  ytry = reflect(pnts, y, psum, ndim, auxdata, func, ibig, 0.5);
68  if (ytry >= ysave) {
69  for (i = 0; i < ndim + 1; i++) {
70  if (i != isml) {
71  for (j = 0; j < ndim; j++)
72  pnts[i * ndim + j] = psum[j] =
73  0.5 * (pnts[i * ndim + j] + pnts[isml * ndim + j]);
74  }
75  }
76  GET_PSUM
77  }
78  }
79 
80  }
81 
82  auxdata->niter = k + 1;
83 
84  return (isml);
85 }
86 
87 double reflect(double *pnts, double y[], double psum[], short ndim,
88  FITSTRUCT *auxdata, double (*func)(FITSTRUCT *, double []),
89  short ibig, float fac) {
90 
91  short j;
92  float fac1, fac2;
93  double ytry, ptry[16];
94 
95  fac1 = (1.0 - fac) / ndim;
96  fac2 = fac1 - fac;
97 
98 
99  for (j = 0; j < ndim; j++) ptry[j] = psum[j] * fac1 - pnts[ibig * ndim + j] * fac2;
100 
101  ytry = (*func)(auxdata, ptry);
102 
103  if (ytry < y[ibig]) {
104  y[ibig] = ytry;
105  for (j = 0; j < ndim; j++) {
106  psum[j] += ptry[j] - pnts[ibig * ndim + j];
107  pnts[ibig * ndim + j] = ptry[j];
108  }
109  }
110 
111  return (ytry);
112 
113 }
114 
115 
116 
117 
118 
119 
120 
121 
122 
int j
Definition: decode_rs.h:73
int niter
Definition: amoeba.h:9
double reflect(double *pnts, double y[], double psum[], short ndim, FITSTRUCT *auxdata, double(*func)(FITSTRUCT *, double[]), short ibig, float fac)
Definition: amoeba.c:87
short amoeba(double *pnts, FITSTRUCT *auxdata, double(*func)(FITSTRUCT *, double[]), float tol)
Definition: amoeba.c:14
#define fac
#define GET_PSUM
Definition: amoeba.c:6
subroutine func(x, conec, n, bconecno, bn, units, u, inno, i, outno, o, Input, Targ, p, sqerr)
Definition: ffnet.f:287
float tol
short int nfunc
Definition: amoeba.h:2
#define fabs(a)
Definition: misc.h:93
int i
Definition: decode_rs.h:71
int k
Definition: decode_rs.h:73