TPCCLIB
Loading...
Searching...
No Matches
simplex.c
Go to the documentation of this file.
1
4/*****************************************************************************/
5#include "tpcclibConfig.h"
6/*****************************************************************************/
7#include <stdio.h>
8#include <stdlib.h>
9#include <math.h>
10#include <time.h>
11#include <string.h>
12/*****************************************************************************/
13#include "tpcextensions.h"
14/*****************************************************************************/
15#include "tpcnlopt.h"
16/*****************************************************************************/
17
18/*****************************************************************************/
32int nloptSimplex(
42 NLOPT *nlo,
44 unsigned int maxIter,
46 TPCSTATUS *status
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}
361/*****************************************************************************/
362
363/*****************************************************************************/
373int nloptSimplexARRS(
382 NLOPT *nlo,
384 unsigned int maxIter,
386 TPCSTATUS *status
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}
615/*****************************************************************************/
616
617/*****************************************************************************/
624int nloptSimplexMS(
633 NLOPT *nlo,
635 unsigned int spNr,
637 TPCSTATUS *status
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}
737/*****************************************************************************/
738
739/*****************************************************************************/
nloptForceLimits
unsigned int nloptForceLimits(unsigned int n, double *pLower, double *pUpper, double *p)
Enforce the model parameters within given limits.
Definition constraints.c:76
doubleCopy
void doubleCopy(double *t, double *s, const unsigned int n)
Definition doubleutil.c:117
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
nlopt1D
int nlopt1D(NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
Definition nlopt1d.c:25
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
nloptFixedNr
unsigned int nloptFixedNr(NLOPT *d)
Definition nlopt.c:354
nloptPrintP
void nloptPrintP(NLOPT *nlo, unsigned int nr, FILE *fp)
Definition nlopt.c:274
nloptAddP
int nloptAddP(NLOPT *nlo, double *p, double funval)
Definition nlopt.c:149
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
nloptSimplexMS
int nloptSimplexMS(NLOPT *nlo, unsigned int spNr, TPCSTATUS *status)
Definition simplex.c:624
nloptSimplex
int nloptSimplex(NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
Definition simplex.c:32
nloptSimplexARRS
int nloptSimplexARRS(NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
Definition simplex.c:373
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::maxFunCalls
unsigned int maxFunCalls
Definition tpcnlopt.h:46
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::plist
double * plist
Definition tpcnlopt.h:56
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
tpcextensions.h
Header file for library libtpcextensions.
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_OUT_OF_MEMORY
@ TPCERROR_OUT_OF_MEMORY
Cannot allocate memory.
Definition tpcextensions.h:193
TPCERROR_OK
@ TPCERROR_OK
No error.
Definition tpcextensions.h:191
TPCERROR_NO_DATA
@ TPCERROR_NO_DATA
File contains no data.
Definition tpcextensions.h:203
tpcnlopt.h
Header file for library libtpcnlopt.
Generated on Sat Nov 29 2025 23:04:55 for TPCCLIB by doxygen 1.12.0