TPCCLIB
Loading...
Searching...
No Matches
Functions
simplex.c File Reference

Downhill simplex nonlinear optimization. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpcextensions.h"
#include "tpcnlopt.h"

Go to the source code of this file.

Functions

int nloptSimplex (NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
int nloptSimplexARRS (NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
int nloptSimplexMS (NLOPT *nlo, unsigned int spNr, TPCSTATUS *status)

Detailed Description

Downhill simplex nonlinear optimization.

Definition in file simplex.c.

Function Documentation

◆ nloptSimplex()

int nloptSimplex ( NLOPT * nlo,
unsigned int maxIter,
TPCSTATUS * status )

Nelder-Mead (downhill simplex) optimization.

If the number of free parameters (not fixed) is one, then bracketing method is used instead of downhill simplex.

Precondition
Initiate the contents of the nlo data structure.
Postcondition
The last objective function call is usually not done with the best parameter estimates; if objective function simulates data that you need, you must call the function with the final parameters.
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
See also
nlopt1D, nloptPowellBrent, nloptSimplexARRS
Parameters
nloPointer to NLOPT structure.
Precondition
Initial guess must be given in x[]. Initial step sizes must be given in xdelta[]. Constraints xlower[] and xupper[] are required, but since the implementation is very simplified, limits should be as wide as possible. Parameter tolerances are used as stopping criteria: Simplex stops if centroid differs from best point less than xtol[], for each parameter; if set to zero, then only stopping rule is the iteration number.
Parameters
maxIterMaximum number of iterations; set to zero to use the default.
statusPointer to status data; enter NULL if not needed.

Definition at line 32 of file simplex.c.

47 {
48 FILE *fp=stdout;
49 int verbose=0; if(status!=NULL) {verbose=status->verbose; fp=status->fp;}
50 if(verbose>1) fprintf(fp, "%s(NLOPT, %d, status)\n", __func__, maxIter);
51 if(nlo==NULL) {
52 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
53 return TPCERROR_FAIL;
54 }
55 if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
56 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
57 return TPCERROR_NO_DATA;
58 }
59
60 /* Check that initial values are inside limits */
61 for(unsigned int i=0; i<nlo->totalNr; i++)
62 if(nlo->xfull[i]<nlo->xlower[i] || nlo->xfull[i]>nlo->xupper[i] ||
63 !isfinite(nlo->xfull[i]) || !isfinite(nlo->xlower[i]) || !isfinite(nlo->xupper[i]))
64 {
65 if(verbose>2) {
66 fprintf(fp, "parameter %d failed: %e <= %e <= %e\n",
67 1+i, nlo->xlower[i], nlo->xupper[i], nlo->xfull[i]);
68 fflush(fp);
69 }
70 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
71 return TPCERROR_INVALID_VALUE;
72 }
73
74 /* If only one parameter free for fitting, then use 1-dimensional method instead */
75 if(nlo->totalNr<2 || nlo->totalNr-nloptFixedNr(nlo)<2) {
76 if(verbose>2) fprintf(fp, "going for one-dimensional method.\n");
77 return(nlopt1D(nlo, maxIter, NULL));
78 }
79
80 /* Initiate the simplex */
81 if(verbose>10) fprintf(fp, "Simplex initialization\n");
82 unsigned int parNr=nlo->totalNr;
83 unsigned int sn=parNr+1; // sn>=3 because parNr >=2 verified before
84 double simplexP[sn+2][parNr]; // including room for first and second new point
85 double simplexC[parNr];
86 double simplexR[sn+2]; // including room for first and second new point
87 for(unsigned int i=0; i<(sn+2); i++) simplexR[i]=nan("");
88 if(verbose>16) fprintf(fp, "simplexR[0..%d]\n", sn+2);
89 unsigned int newPnt=sn;
90 unsigned int newPnt2=newPnt+1;
91 /* Check the maximum iteration number */
92 if(maxIter<2) maxIter=500;
93
94 /* Fill with initial guesses */
95 for(unsigned int j=0; j<sn; j++)
96 for(unsigned int i=0; i<parNr; i++)
97 simplexP[j][i]=nlo->xfull[i];
98
99 for(unsigned int j=1; j<sn; j++) { // initial guess is kept as such in initial simplex
100 for(unsigned int i=0; i<parNr; i++) {
101 simplexP[j][i] = simplexP[j-1][i];
102 if((j-1)==i && nlo->xupper[i]>nlo->xlower[i]) {
103 simplexP[j][i] += nlo->xdelta[i];
104 if(simplexP[j][i]<nlo->xlower[i] || simplexP[j][i]>nlo->xupper[i]) {
105 /* parameter would exceed limits, try the other direction */
106 simplexP[j][i] -= 2.0*nlo->xdelta[i];
107 /* check again */
108 if(simplexP[j][i]<nlo->xlower[i] || simplexP[j][i]>nlo->xupper[i]) {
109 /* stupid limits or delta */
110 if(verbose>2) {
111 fprintf(fp, "failed limits %e,%e or delta: %e\n", nlo->xlower[i], nlo->xupper[i], nlo->xdelta[i]);
112 fflush(fp);
113 }
114 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
115 return TPCERROR_INVALID_VALUE;
116 }
117 }
118 }
119 }
120 }
121 /* Calculate initial simplex function values */
122 for(unsigned int j=0; j<sn; j++) {
123 if(verbose>13) {
124 fprintf(fp, "x[%d][]=", j);
125 for(unsigned int i=0; i<parNr; i++) fprintf(fp, "%g ", simplexP[j][i]);
126 fprintf(fp, "\n");
127 }
128 simplexR[j] = (*nlo->_fun)(parNr, simplexP[j], nlo->fundata);
129 nlo->funCalls++;
130 if(verbose>12) fprintf(fp, " R[%d]=%g\n", j, simplexR[j]);
131 if(nlo->usePList) {
132 int ret=nloptAddP(nlo, simplexP[j], simplexR[j]);
133 if(ret!=TPCERROR_OK) {
134 statusSet(status, __func__, __FILE__, __LINE__, ret);
135 return(ret);
136 }
137 }
138 }
139 double initialR=simplexR[0];
140
141 /* Find the max and min response measured */
142 unsigned int best=0, nextBest=0, worst=0;
143 {
144 double max, min, min2;
145 max=min=simplexR[0]; best=worst=0; // initial guess
146 for(unsigned int j=1; j<sn; j++) {
147 if(!isfinite(max) || simplexR[j] > max) {max=simplexR[j]; worst=j;}
148 if(!isfinite(min) || simplexR[j] < min) {min=simplexR[j]; best=j; }
149 }
150 /* Find the 2nd best, too */
151 min2=max; nextBest=worst;
152 for(unsigned int j=0; j<sn; j++)
153 if(j!=best && (simplexR[j] < min2)) {min2=simplexR[j]; nextBest=j;}
154 if(verbose>14) {
155 fprintf(fp, " best := %g (%d)\n", min, best);
156 fprintf(fp, " next best := %g (%d)\n", min2, nextBest);
157 fprintf(fp, " worst := %g (%d)\n", max, worst);
158 }
159 }
160
161 /* Simplex iterations */
162 if(verbose>10) fprintf(fp, "Simplex iterations\n");
163 double prevBestR=simplexR[best];
164 unsigned int iterNr=0, stopTolerance=0, stopR=0;
165 do {
166 iterNr++;
167 if(verbose>13) {fprintf(fp, "\niteration := %d\n", iterNr); fflush(fp);}
168 if(verbose>16) {
169 for(unsigned int j=0; j<sn; j++) {
170 fprintf(fp, "simplexP[%d][]=", j);
171 for(unsigned int i=0; i<parNr; i++) fprintf(fp, "%g ", simplexP[j][i]);
172 fprintf(fp, " --> %g ", simplexR[j]);
173 if(j==best) fprintf(fp, " (best)\n");
174 else if(j==nextBest) fprintf(fp, " (next best)\n");
175 else if(j==worst) fprintf(fp, " (worst)\n");
176 else fprintf(fp, "\n");
177 fflush(fp);
178 }
179 }
180 /* Calculate centroid of all measurements except the worst */
181 for(unsigned int i=0; i<parNr; i++) {
182 simplexC[i]=0.0; unsigned int n=0;
183 for(unsigned int j=0; j<sn; j++)
184 if(j!=worst && isfinite(simplexP[j][i])) {simplexC[i]+=simplexP[j][i]; n++;}
185 simplexC[i]/=(double)(n);
186 //Centroid is allowed to step over limits, because it is only used to move points
187 }
188 if(verbose>15) {
189 fprintf(fp, "centroid x[]=");
190 for(unsigned int i=0; i<parNr; i++) fprintf(fp, "%g ", simplexC[i]);
191 fprintf(fp, "\n"); fflush(fp);
192 }
193
194 /* Stopping rule: If simplex points do not differ more than tolerance, then that can mean
195 the end is near */
196 {
197 unsigned int k;
198 for(k=0; k<parNr; k++) {
199 if(fabs(simplexC[k]-simplexP[worst][k])>0.5*nlo->xtol[k]) break;
200 if(fabs(simplexC[k]-simplexP[best][k])>0.5*nlo->xtol[k]) break;
201 if(fabs(simplexC[k]-simplexP[nextBest][k])>0.5*nlo->xtol[k]) break;
202 }
203 if(k==parNr) stopTolerance++; else stopTolerance=0;
204 }
205 /* If two stopping rules are filled at the same time, then stop */
206 if(stopTolerance>=parNr && stopR>=parNr) { // break the do - while loop
207 if(verbose>2) {fprintf(fp, " stopping because simplex is not improving.\n"); fflush(fp);}
208 break;
209 }
210 /* Create new point as centroid reflected away from worst */
211 double f=1.0;
212 for(unsigned int i=0; i<parNr; i++)
213 simplexP[newPnt][i] = simplexC[i] + f*(simplexC[i]-simplexP[worst][i]);
214 nloptForceLimits(parNr, nlo->xlower, nlo->xupper, simplexP[newPnt]);
215 /* and compute the response */
216 simplexR[newPnt] = (*nlo->_fun)(parNr, simplexP[newPnt], nlo->fundata);
217 nlo->funCalls++;
218 if(verbose>15) {
219 fprintf(fp, "new point x[%d][]=", newPnt);
220 for(unsigned int i=0; i<parNr; i++) fprintf(fp, "%g ", simplexP[newPnt][i]);
221 fprintf(fp, "\n -> R[%d]=%g\n", newPnt, simplexR[newPnt]);
222 }
223 if(nlo->usePList) {
224 int ret=nloptAddP(nlo, simplexP[newPnt], simplexR[newPnt]);
225 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
226 }
227 /* Depending on the response, decide what to do next */
228 if(simplexR[newPnt] < simplexR[best]) {
229 if(verbose>14) fprintf(fp, " new is better than best\n");
230 /* If this new point is better than previous best, then expand in this direction */
231 f=2.0;
232 for(unsigned int i=0; i<parNr; i++)
233 simplexP[newPnt2][i] = simplexC[i] + f*(simplexC[i]-simplexP[worst][i]);
234 nloptForceLimits(parNr, nlo->xlower, nlo->xupper, simplexP[newPnt2]);
235 simplexR[newPnt2] = (*nlo->_fun)(parNr, simplexP[newPnt2], nlo->fundata);
236 nlo->funCalls++;
237 if(nlo->usePList) {
238 int ret=nloptAddP(nlo, simplexP[newPnt2], simplexR[newPnt2]);
239 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
240 }
241 if(simplexR[newPnt2] < simplexR[newPnt]) {
242 if(verbose>14) fprintf(fp, " expanded (%g) is better than new\n", simplexR[newPnt2]);
243 for(unsigned int i=0; i<parNr; i++) simplexP[worst][i]=simplexP[newPnt2][i];
244 simplexR[worst]=simplexR[newPnt2];
245 } else {
246 if(verbose>14) fprintf(fp, " expanded (%g) is no better than new\n", simplexR[newPnt2]);
247 for(unsigned int i=0; i<parNr; i++) simplexP[worst][i]=simplexP[newPnt][i];
248 simplexR[worst]=simplexR[newPnt];
249 }
250 } else if(!isfinite(simplexR[newPnt]) || simplexR[newPnt] > simplexR[worst]) {
251 if(verbose>14) fprintf(fp, " new is worse than worst\n");
252 /* If new point is worse than previous worst, measure point halfway between
253 the worst and centroid */
254 f=-0.5;
255 for(unsigned int i=0; i<parNr; i++)
256 simplexP[newPnt2][i] = simplexC[i] + f*(simplexC[i]-simplexP[worst][i]);
257 nloptForceLimits(parNr, nlo->xlower, nlo->xupper, simplexP[newPnt2]);
258 simplexR[newPnt2] = (*nlo->_fun)(parNr, simplexP[newPnt2], nlo->fundata);
259 nlo->funCalls++;
260 if(nlo->usePList) {
261 int ret=nloptAddP(nlo, simplexP[newPnt2], simplexR[newPnt2]);
262 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
263 }
264 if(simplexR[newPnt2] < simplexR[newPnt] && isfinite(simplexR[newPnt2])) {
265 if(verbose>14) fprintf(fp, " halfway (%g) is better than new\n", simplexR[newPnt2]);
266 for(unsigned int i=0; i<parNr; i++) simplexP[worst][i]=simplexP[newPnt2][i];
267 simplexR[worst]=simplexR[newPnt2];
268 } else if(isfinite(simplexR[newPnt])) {
269 if(verbose>14) fprintf(fp, " halfway (%g) is no better than new\n", simplexR[newPnt2]);
270 for(unsigned int i=0; i<parNr; i++) simplexP[worst][i]=simplexP[newPnt][i];
271 simplexR[worst]=simplexR[newPnt];
272 }
273 } else if(simplexR[newPnt] > simplexR[nextBest]) {
274 if(verbose>14) fprintf(fp, " new is worse than next best\n");
275 /* If newest response is worse than next best point but better than worst,
276 measure response halfway between centroid and the newest point */
277 f=0.5;
278 for(unsigned int i=0; i<parNr; i++)
279 simplexP[newPnt2][i] = simplexC[i] + f*(simplexP[newPnt][i]-simplexC[i]);
280 nloptForceLimits(parNr, nlo->xlower, nlo->xupper, simplexP[newPnt2]);
281 simplexR[newPnt2] = (*nlo->_fun)(parNr, simplexP[newPnt2], nlo->fundata);
282 nlo->funCalls++;
283 if(nlo->usePList) {
284 int ret=nloptAddP(nlo, simplexP[newPnt2], simplexR[newPnt2]);
285 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
286 }
287 if(simplexR[newPnt2] < simplexR[newPnt]) {
288 if(verbose>14) fprintf(fp, " halfway (%g) is better than new\n", simplexR[newPnt2]);
289 for(unsigned int i=0; i<parNr; i++) simplexP[worst][i]=simplexP[newPnt2][i];
290 simplexR[worst]=simplexR[newPnt2];
291 } else if(isfinite(simplexR[newPnt])) {
292 if(verbose>14) fprintf(fp, " halfway (%g) is worse than new\n", simplexR[newPnt2]);
293 for(unsigned int i=0; i<parNr; i++) simplexP[worst][i]=simplexP[newPnt][i];
294 simplexR[worst]=simplexR[newPnt];
295 }
296 } else if(isfinite(simplexR[newPnt])) {
297 if(verbose>14) fprintf(fp, " new is worse than the best but better than the next best\n");
298 /* If none of the above is true, then replace the second best point with this */
299 for(unsigned int i=0; i<parNr; i++) simplexP[worst][i]=simplexP[newPnt][i];
300 simplexR[worst]=simplexR[newPnt];
301 }
302
303 /* Find the max and min response measured */
304 {
305 double max, min, min2;
306 max=min=simplexR[0]; best=worst=0;
307 for(unsigned int j=1; j<sn; j++) {
308 if(!isfinite(max) || simplexR[j] > max) {max=simplexR[j]; worst=j;}
309 if(!isfinite(min) || simplexR[j] < min) {min=simplexR[j]; best=j; }
310 }
311 /* Find the 2nd best, too */
312 min2=max; nextBest=worst;
313 for(unsigned int j=0; j<sn; j++)
314 if(j!=best && (simplexR[j] < min2)) {min2=simplexR[j]; nextBest=j;}
315 if(verbose>14) {
316 fprintf(fp, " best := %g (%d)\n", min, best);
317 fprintf(fp, " next best := %g (%d)\n", min2, nextBest);
318 fprintf(fp, " worst := %g (%d)\n", max, worst);
319 }
320 }
321
322 /* Stopping rule: count how many consequential times the R has not improved */
323 if(simplexR[best]<prevBestR) stopR=0; else stopR++;
324 prevBestR=simplexR[best];
325
326 } while(iterNr<maxIter);
327
328 if(verbose>2) {
329 if(iterNr>=maxIter) fprintf(fp, " stopped because max iterations reached.\n");
330 fprintf(fp, " %d iterations\n", iterNr);
331 fprintf(fp, " R[%d]=%g\n", best, simplexR[best]);
332 fflush(fp);
333 }
334
335#if(0)
336 /* One more function call with the best parameters, to make sure that
337 if objective function simulates data, it is simulated with the best fit */
338 simplexR[best] = (*nlo->_fun)(parNr, simplexP[best], nlo->fundata);
339#endif
340
341 /* Check the objective function value */
342 if(!isfinite(simplexR[best])) {
343 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
344 return(TPCERROR_BAD_FIT);
345 }
346
347 /* Copy optimized parameters over initial values */
348 if(simplexR[best]<=initialR) {
349 for(unsigned int i=0; i<parNr; i++) nlo->xfull[i]=simplexP[best][i];
350 nlo->funval=simplexR[best];
351 } else {
352 /* This happens only in case of a bug */
353 if(1 || verbose>0) fprintf(fp, "nloptSimplex() gives worse R than initially!\n");
354 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
355 return(TPCERROR_BAD_FIT);
356 }
357
358 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
359 return TPCERROR_OK;
360}
nloptForceLimits
unsigned int nloptForceLimits(unsigned int n, double *pLower, double *pUpper, double *p)
Enforce the model parameters within given limits.
Definition constraints.c:76
nlopt1D
int nlopt1D(NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
Definition nlopt1d.c:25
nloptFixedNr
unsigned int nloptFixedNr(NLOPT *d)
Definition nlopt.c:354
nloptAddP
int nloptAddP(NLOPT *nlo, double *p, double funval)
Definition nlopt.c:149
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
NLOPT::_fun
double(* _fun)(int, double *, void *)
Definition tpcnlopt.h:42
NLOPT::funval
double funval
Definition tpcnlopt.h:50
NLOPT::xupper
double * xupper
Definition tpcnlopt.h:33
NLOPT::xlower
double * xlower
Definition tpcnlopt.h:31
NLOPT::fundata
void * fundata
Definition tpcnlopt.h:44
NLOPT::xfull
double * xfull
Definition tpcnlopt.h:29
NLOPT::funCalls
unsigned int funCalls
Definition tpcnlopt.h:48
NLOPT::xdelta
double * xdelta
Definition tpcnlopt.h:36
NLOPT::usePList
unsigned int usePList
Definition tpcnlopt.h:52
NLOPT::xtol
double * xtol
Definition tpcnlopt.h:39
NLOPT::totalNr
unsigned int totalNr
Definition tpcnlopt.h:27
TPCSTATUS::verbose
int verbose
Verbose level, used by statusPrint() etc.
Definition tpcextensions.h:242
TPCSTATUS::fp
FILE * fp
File pointer for writing log information during development and testing.
Definition tpcextensions.h:243
TPCERROR_INVALID_VALUE
@ TPCERROR_INVALID_VALUE
Invalid value.
Definition tpcextensions.h:220
TPCERROR_BAD_FIT
@ TPCERROR_BAD_FIT
Fitting not successful.
Definition tpcextensions.h:228
TPCERROR_FAIL
@ TPCERROR_FAIL
General error.
Definition tpcextensions.h:192
TPCERROR_OK
@ TPCERROR_OK
No error.
Definition tpcextensions.h:191
TPCERROR_NO_DATA
@ TPCERROR_NO_DATA
File contains no data.
Definition tpcextensions.h:203

Referenced by nloptIATGO(), nloptITGO1(), nloptITGO2(), nloptMPSO(), nloptSimplexARRS(), and nloptSimplexMS().

◆ nloptSimplexARRS()

int nloptSimplexARRS ( NLOPT * nlo,
unsigned int maxIter,
TPCSTATUS * status )

Accumulative Restarting Random start-point Nelder-Mead (downhill simplex) optimization.

Precondition
Uses rand(), therefore set seed for a new series of pseudo-random numbers; to produce different set of pseudo-randoms, call drandSeed(1) before calling this function.
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
See also
drandSeed(), nloptSimplex, nlopt1D, nloptPowellBrent
Parameters
nloPointer to NLOPT structure.
Precondition
Constraints xlower[] and xupper[] are required. Parameter tolerances xtol[] are required, and used as stopping criteria. Other stopping criteria are the iteration number and the NLOPT maxFunCalls. Initial guess can be given in x[], but it is not obligatory. Initial step sizes (xdelta[]) are not used.
See also
nloptInit, nloptFree, nloptAllocate
Parameters
maxIterMaximum number of iterations; set to zero to use the default.
statusPointer to status data; enter NULL if not needed.

Definition at line 373 of file simplex.c.

387 {
388 FILE *fp=stdout;
389 int verbose=0; if(status!=NULL) {verbose=status->verbose; fp=status->fp;}
390 //verbose=10;
391
392 if(verbose>1) fprintf(fp, "%s(NLOPT, %d, status)\n", __func__, maxIter);
393 if(nlo==NULL) {
394 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL); return TPCERROR_FAIL;}
395 if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
396 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA); return TPCERROR_NO_DATA;}
397
398 /* Check if any of the parameters are fixed */
399 unsigned int fixedParNr=nloptLimitFixedNr(nlo);
400 if(verbose>2 && fixedParNr>0) fprintf(fp, "fixedParNr := %d\n", fixedParNr);
401 unsigned int fittedParNr=nlo->totalNr-fixedParNr;
402 if(verbose>2) fprintf(fp, "fittedParNr := %d\n", fittedParNr);
403 if(fittedParNr<1) {
404 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
405 return TPCERROR_INVALID_VALUE;
406 }
407
408 /* Check the tolerations */
409 for(unsigned int i=0; i<nlo->totalNr; i++) {
410 if(nlo->xlower[i]>=nlo->xupper[i]) {nlo->xtol[i]=0.0; continue;}
411 if(!(nlo->xtol[i]>0.0)) {
412 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
413 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
414 return TPCERROR_INVALID_VALUE;
415 }
416 }
417
418
419 /* Set maxIter, if necessary */
420 if(maxIter==0) maxIter=5+fittedParNr;
421 if(verbose>2) {fprintf(fp, "maxIter := %u\n", maxIter); fflush(fp);}
422
423 /* Store the original tolerances, and use smaller tolerance with downhill simplex */
424 unsigned int dim=nlo->totalNr;
425 double tol[dim];
426 for(unsigned int i=0; i<dim; i++) {tol[i]=nlo->xtol[i]; nlo->xtol[i]*=0.1;}
427
428 /*
429 * Set initial xdelta based on parameter limits and original tolerances
430 */
431 for(unsigned int i=0; i<dim; i++) {
432 double d=nlo->xupper[i]-nlo->xlower[i];
433 if(d>0.0) nlo->xdelta[i]=2.0*tol[i]+0.15*d; else nlo->xdelta[i]=0.0;
434 }
435
436 /*
437 * Initialize the sample list with random points
438 */
439
440 /* Allocate memory */
441 unsigned int sampleNr=30*fittedParNr;
442 if(verbose>2) {fprintf(fp, "sampleNr := %u\n", sampleNr); fflush(fp);}
443 nlo->usePList=1; nlo->funCalls=0;
444 nlo->plist=(double*)malloc(sampleNr*(dim+1)*sizeof(double));
445 if(nlo->plist==NULL) { // will be freed with nloptFree()
446 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
447 return TPCERROR_OUT_OF_MEMORY;
448 }
449
450 /* If initial guess was given by user, create random points with Gaussian distribution around it */
451 int initGuessAvailable=1;
452 for(unsigned int i=0; i<nlo->totalNr; i++) {
453 if(!isfinite(nlo->xfull[i])) {initGuessAvailable=0; break;}
454 if(nlo->xfull[i]<nlo->xlower[i]) {initGuessAvailable=0; break;}
455 if(nlo->xfull[i]>nlo->xupper[i]) {initGuessAvailable=0; break;}
456 }
457 nlo->funval=(*nlo->_fun)(dim, nlo->xfull, nlo->fundata);
458 if(isfinite(nlo->funval)) {
459 if(verbose>2) {
460 fprintf(fp, "Initial guess: ");
461 for(unsigned int i=0; i<dim; i++) fprintf(fp, "%g ", nlo->xfull[i]);
462 fprintf(fp, "=> %e\n", nlo->funval); fflush(fp);
463 }
464 doubleCopy(nlo->plist, nlo->xfull, dim);
465 nlo->plist[dim]=nlo->funval; nlo->funCalls++;
466 } else {
467 if(verbose>2) {fprintf(fp, "invalid initial guess\n"); fflush(fp);}
468 initGuessAvailable=0;
469 }
470 if(initGuessAvailable) {
471 for(unsigned int pi=1; pi<sampleNr; pi++) {
472 nloptGaussianPoint(&nlo->plist[pi*(dim+1)], nlo->plist, nlo->xdelta,
473 nlo->xlower, nlo->xupper, dim, NULL);
474 nlo->plist[pi*(dim+1)+dim]=(*nlo->_fun)(dim, &nlo->plist[pi*(dim+1)], nlo->fundata);
475 nlo->funCalls++;
476 while(!isfinite(nlo->plist[pi*(dim+1)+dim])) { // new guess until valid result
477 nloptGaussianPoint(&nlo->plist[pi*(dim+1)], nlo->plist, nlo->xdelta,
478 nlo->xlower, nlo->xupper, dim, NULL);
479 nlo->plist[pi*(dim+1)+dim]=(*nlo->_fun)(dim, &nlo->plist[pi*(dim+1)], nlo->fundata);
480 // Do NOT add funCalls because that is also the length of plist
481 }
482 }
483 }
484
485 /* If valid initial point was not given, fill the list with random points */
486 if(!initGuessAvailable) {
487 for(unsigned int pi=0; pi<sampleNr; pi++) {
488 nloptRandomPoint(&nlo->plist[pi*(dim+1)], nlo->xlower, nlo->xupper, dim, NULL);
489 nlo->plist[pi*(dim+1)+dim]=(*nlo->_fun)(dim, &nlo->plist[pi*(dim+1)], nlo->fundata);
490 nlo->funCalls++;
491 while(!isfinite(nlo->plist[pi*(dim+1)+dim])) { // new guess until valid result
492 nloptGaussianPoint(&nlo->plist[pi*(dim+1)], nlo->plist, nlo->xdelta,
493 nlo->xlower, nlo->xupper, dim, NULL);
494 nlo->plist[pi*(dim+1)+dim]=(*nlo->_fun)(dim, &nlo->plist[pi*(dim+1)], nlo->fundata);
495 // Do NOT add funCalls because that is also the length of plist
496 }
497 }
498 }
499
500 /* Sort samples based on the evaluated function value */
501 if(nloptSortP(nlo)!=TPCERROR_OK) {
502 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL); return TPCERROR_FAIL;}
503 if(verbose>10) nloptPrintP(nlo, 0, fp);
504
505 /* Downhill simplex from the best point so far */
506 doubleCopy(nlo->xfull, nlo->plist, dim);
507 if(verbose>4) {
508 fprintf(fp, "LO: initial point with deltas:\n");
509 for(unsigned int i=0; i<dim; i++) fprintf(fp, "\t%e\t%e\n", nlo->xfull[i], nlo->xdelta[i]);
510 }
511 if(nloptSimplex(nlo, 100*dim, status)!=TPCERROR_OK) {
512 if(verbose>1) {fprintf(fp, " LO failed\n"); fflush(fp);}
513 return(TPCERROR_BAD_FIT);
514 } else {
515 if(verbose>4) {
516 fprintf(fp, "Point after LO:");
517 for(unsigned int i=0; i<dim; i++) fprintf(fp, " %e ", nlo->xfull[i]);
518 fprintf(fp, " => %e\n", nlo->funval); fflush(fp);
519 }
520 if(verbose>5) fprintf(fp, "funCalls=%d\n", nlo->funCalls);
521 }
522
523 /*
524 * Start iterations
525 */
526 unsigned int iterNr=0; // loop index
527 while(iterNr<maxIter && nlo->funCalls<nlo->maxFunCalls) {
528 iterNr++;
529 if(verbose>4) {fprintf(fp, "-----------------------------\nIteration %d\n", iterNr); fflush(fp);}
530
531 /* Sort samples based on the evaluated function value */
532 if(nloptSortP(nlo)!=TPCERROR_OK) {
533 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL); return TPCERROR_FAIL;}
534 if(verbose>12) nloptPrintP(nlo, 20, fp);
535 if(verbose>5) {
536 fprintf(fp, "best point so far:");
537 for(unsigned int i=0; i<dim; i++) fprintf(fp, " %e", nlo->plist[i]);
538 fprintf(fp, " => %e\n", nlo->plist[dim]); fflush(fp);
539 }
540
541 /* Calculate the parameter means and SDs from the best part of points so far */
542 double mean[dim], sd[dim];
543// if(nloptMeanP(nlo, nlo->funCalls/(2*(1+iterNr)), mean, sd)!=TPCERROR_OK) {
544 if(nloptMeanP(nlo, nlo->funCalls/5, mean, sd)!=TPCERROR_OK) {
545 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT); return TPCERROR_BAD_FIT;}
546 if(verbose>10) {
547 fprintf(fp, "Mean and SD of best points so far:\n");
548 for(unsigned int i=0; i<dim; i++) fprintf(fp, "\t%e\t%e\n", mean[i], sd[i]);
549 }
550 if(verbose>20) nloptPrintP(nlo, nlo->funCalls/5, fp);
551
552 /* If SDs were smaller than tolerances, then stop */
553 unsigned int i; for(i=0; i<dim; i++) if(fabs(sd[i])>tol[i]) break;
554 if(i==dim) {
555 if(verbose>1) fprintf(fp, "\n Required tolerance reached.\n");
556 if(iterNr>2) break;
557 }
558
559 /* New start point based on the mean point and the best point so far,
560 and xdeltas based on the SD */
561 for(unsigned int i=0; i<dim; i++) {
562 if(nlo->xupper[i]>nlo->xlower[i]) {
563 nlo->xfull[i]=0.5*mean[i]+0.5*nlo->plist[i];
564 nlo->xdelta[i]=sd[i];
565 /* Make sure that individual is not too close to zero */
566 if(nlo->xdelta[i]<tol[i]) nlo->xdelta[i]=10.0*tol[i];
567 } else {
568 nlo->xfull[i]=nlo->xlower[i];
569 nlo->xdelta[i]=0.0;
570 }
571 }
572
573 /* Downhill simplex */
574 if(verbose>6) {
575 fprintf(fp, "LO: initial point with deltas:\n");
576 for(unsigned int i=0; i<dim; i++) fprintf(fp, "\t%e\t%e\n", nlo->xfull[i], nlo->xdelta[i]);
577 }
578 if(nloptSimplex(nlo, 100*dim, status)!=TPCERROR_OK) {
579 if(verbose>1) {fprintf(fp, " LO failed\n"); fflush(fp);}
580 break;
581 } else {
582 if(verbose>6) {
583 fprintf(fp, "Point after LO:");
584 for(unsigned int i=0; i<dim; i++) fprintf(fp, " %e ", nlo->xfull[i]);
585 fprintf(fp, " => %e\n", nlo->funval); fflush(fp);
586 }
587 }
588
589 if(verbose>5) fprintf(fp, "funCalls=%d\n", nlo->funCalls);
590
591 } // next iteration
592
593 /* Check the reason for loop exit */
594 if(iterNr>=maxIter) {
595 if(verbose>1) fprintf(fp, "\n Exceeded the maximum number for loops.\n");
596 }
597 if(nlo->funCalls>=nlo->maxFunCalls) {
598 if(verbose>1) fprintf(fp, "\n Exceeded the maximum number for function calls.\n");
599 }
600
601 /* Sort samples based on the evaluated function value */
602 if(nloptSortP(nlo)!=TPCERROR_OK) {
603 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL); return TPCERROR_FAIL;}
604
605 /* Get the best point so far */
606 for(unsigned int i=0; i<dim; i++) nlo->xfull[i]=nlo->plist[i];
607 nlo->funval=nlo->plist[dim];
608
609 /* Put back the original tolerances */
610 for(unsigned int i=0; i<dim; i++) nlo->xtol[i]=tol[i];
611
612 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
613 return TPCERROR_OK;
614}
doubleCopy
void doubleCopy(double *t, double *s, const unsigned int n)
Definition doubleutil.c:117
nloptLimitFixedNr
unsigned int nloptLimitFixedNr(NLOPT *d)
Definition nlopt.c:333
nloptMeanP
int nloptMeanP(NLOPT *nlo, unsigned int nr, double *meanp, double *sdp)
Definition nlopt.c:219
nloptPrintP
void nloptPrintP(NLOPT *nlo, unsigned int nr, FILE *fp)
Definition nlopt.c:274
nloptSortP
int nloptSortP(NLOPT *nlo)
Definition nlopt.c:182
nloptRandomPoint
int nloptRandomPoint(double *p, double *low, double *up, unsigned int n, MERTWI *mt)
Definition rndpoint.c:28
nloptGaussianPoint
int nloptGaussianPoint(double *p, double *mean, double *sd, double *low, double *up, unsigned int n, MERTWI *mt)
Definition rndpoint.c:70
nloptSimplex
int nloptSimplex(NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
Definition simplex.c:32
NLOPT::maxFunCalls
unsigned int maxFunCalls
Definition tpcnlopt.h:46
NLOPT::plist
double * plist
Definition tpcnlopt.h:56
TPCERROR_OUT_OF_MEMORY
@ TPCERROR_OUT_OF_MEMORY
Cannot allocate memory.
Definition tpcextensions.h:193

◆ nloptSimplexMS()

int nloptSimplexMS ( NLOPT * nlo,
unsigned int spNr,
TPCSTATUS * status )

Multi-Start-point Nelder-Mead (downhill simplex) optimization.

Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
See also
nloptSimplex, nlopt1D, nloptPowellBrent
Parameters
nloPointer to NLOPT structure.
Precondition
Constraints xlower[] and xupper[] are required. Parameter tolerances xtol[] are required, and used as stopping criteria. Other stopping criteria are the iteration number and the NLOPT maxFunCalls. Initial guess can be given in x[], but it is not obligatory. Initial step sizes (xdelta[]) are not used.
See also
nloptInit, nloptFree, nloptAllocate
Parameters
spNrNumber of start points; set to zero to use the default.
statusPointer to status data; enter NULL if not needed.

Definition at line 624 of file simplex.c.

638 {
639 FILE *fp=stdout;
640 int verbose=0; if(status!=NULL) {verbose=status->verbose; fp=status->fp;}
641 //verbose=10;
642
643 if(verbose>1) fprintf(fp, "%s(NLOPT, %d, status)\n", __func__, spNr);
644 if(nlo==NULL) {
645 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL); return TPCERROR_FAIL;}
646 if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
647 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA); return TPCERROR_NO_DATA;}
648
649 /* Check if any of the parameters are fixed */
650 unsigned int fixedParNr=nloptLimitFixedNr(nlo);
651 if(verbose>2 && fixedParNr>0) fprintf(fp, "fixedParNr := %d\n", fixedParNr);
652 unsigned int fittedParNr=nlo->totalNr-fixedParNr;
653 if(verbose>2) fprintf(fp, "fittedParNr := %d\n", fittedParNr);
654 if(fittedParNr<1) {
655 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
656 return TPCERROR_INVALID_VALUE;
657 }
658
659 /* Check the tolerations */
660 for(unsigned int i=0; i<nlo->totalNr; i++) {
661 if(nlo->xlower[i]>=nlo->xupper[i]) {nlo->xtol[i]=0.0; continue;}
662 if(!(nlo->xtol[i]>0.0)) {
663 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
664 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
665 return TPCERROR_INVALID_VALUE;
666 }
667 }
668
669 /* Set the number of start points, if necessary */
670 if(spNr==0) spNr=100*fittedParNr;
671 if(verbose>2) {fprintf(fp, "spNr := %u\n", spNr); fflush(fp);}
672
673 /* Set start points */
674 unsigned int xNr=nlo->totalNr;
675 double ss[spNr], sp[spNr*xNr], he[fittedParNr];
676 for(unsigned int si=0; si<spNr; si++) {
677 ss[si]=nan("");
678 if(haltonElement(si, fittedParNr, he)) {
679 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL); return TPCERROR_FAIL;}
680 unsigned int hi=0;
681 for(unsigned int xi=0; xi<xNr; xi++) {
682 if(nlo->xtol[xi]>0.0) {
683 sp[si*xNr+xi] = nlo->xlower[xi] + (nlo->xupper[xi]-nlo->xlower[xi])*he[hi++];
684 } else sp[si*xNr+xi] = nlo->xfull[xi];
685 }
686 }
687 if(verbose>3) {
688 printf("Halton-based start points:\n");
689 for(unsigned int si=0; si<spNr; si++) {
690 for(unsigned int xi=0; xi<xNr; xi++) printf(" %g", sp[si*xNr+xi]);
691 printf(" -> %g\n", ss[si]);
692 }
693 }
694
695 /* Perform Simplex optimization from each starting point */
696 /* Set xdeltas based on parameter limits and tolerances
697 */
698 for(unsigned int xi=0; xi<xNr; xi++) {
699 double d=nlo->xupper[xi]-nlo->xlower[xi];
700 if(d>0.0) nlo->xdelta[xi]=2.0*nlo->xtol[xi]+0.01*d; else nlo->xdelta[xi]=0.0;
701 }
702 //status->verbose=3;
703 for(unsigned int si=0; si<spNr; si++) {
704 ss[si]=nan("");
705 for(unsigned int xi=0; xi<xNr; xi++)
706 nlo->xfull[xi]=sp[si*xNr+xi];
707 if(nloptSimplex(nlo, 100*xNr, status)!=TPCERROR_OK) {
708 if(verbose>1) {fprintf(fp, " LO failed\n"); fflush(fp);}
709 return(TPCERROR_BAD_FIT);
710 }
711 for(unsigned int xi=0; xi<xNr; xi++)
712 sp[si*xNr+xi]=nlo->xfull[xi];
713 ss[si]=nlo->funval;
714 }
715 if(verbose>3) {
716 printf("Simplex results from Halton-based start points:\n");
717 for(unsigned int si=0; si<spNr; si++) {
718 for(unsigned int xi=0; xi<xNr; xi++) printf(" %g", sp[si*xNr+xi]);
719 printf(" -> %g\n", ss[si]);
720 }
721 }
722
723 /* Find the best result */
724 unsigned int bi=0;
725 double bv=nan("");
726 for(unsigned int si=0; si<spNr; si++) {
727 if(!isfinite(ss[si])) continue;
728 if(!isfinite(bv) || ss[si]<bv) {bv=ss[si]; bi=si;}
729 }
730 for(unsigned int xi=0; xi<xNr; xi++) nlo->xfull[xi]=sp[bi*xNr+xi];
731 nlo->funval=ss[bi];
732
733
734 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
735 return TPCERROR_OK;
736}
haltonElement
int haltonElement(const int i, const int dim, double *r)
Calculation of an element of quasi-random low-discrepancy Halton sequence.
Definition halton.c:38
Generated on for TPCCLIB by doxygen 1.15.0