tpcnlopt.h File Reference
Header file for library libtpcnlopt. More...
#include "tpcclibConfig.h"#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>#include "tpcextensions.h"#include "tpcrand.h"Go to the source code of this file.
Data Structures | |
| struct | NLOPT |
| struct | NLOPT_DATA |
Functions | |
| void | nloptInit (NLOPT *nlo) |
| void | nloptFree (NLOPT *nlo) |
| int | nloptAllocate (NLOPT *nlo, unsigned int parNr) |
| int | nloptDuplicate (NLOPT *nlo1, NLOPT *nlo2) |
| int | nloptAddP (NLOPT *nlo, double *p, double funval) |
| int | nloptSortP (NLOPT *nlo) |
| int | nloptMeanP (NLOPT *nlo, unsigned int nr, double *meanp, double *sdp) |
| void | nloptPrintP (NLOPT *nlo, unsigned int nr, FILE *fp) |
| void | nloptWrite (NLOPT *d, FILE *fp) |
| unsigned int | nloptLimitFixedNr (NLOPT *d) |
| unsigned int | nloptFixedNr (NLOPT *d) |
| void | nloptRemoveEmpties (NLOPT *d) |
| void | nloptdataInit (NLOPT_DATA *d) |
| int | nloptRandomPoint (double *p, double *low, double *up, unsigned int n, MERTWI *mt) |
| int | nloptGaussianPoint (double *p, double *mean, double *sd, double *low, double *up, unsigned int n, MERTWI *mt) |
| unsigned int | nloptCheckParameters (unsigned int n, double *pLower, double *pUpper, double *p, double *pAccept, double *penalty) |
| Check that model parameters are within given limits. | |
| unsigned int | nloptForceLimits (unsigned int n, double *pLower, double *pUpper, double *p) |
| Enforce the model parameters within given limits. | |
| int | nlopt1D (NLOPT *nlo, unsigned int maxIter, TPCSTATUS *) |
| int | nloptSimplex (NLOPT *nlo, unsigned int maxIter, TPCSTATUS *) |
| int | nloptSimplexARRS (NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status) |
| int | nloptSimplexMS (NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status) |
| int | nloptPowellBrent (NLOPT *nlo, unsigned int ktm, TPCSTATUS *) |
| Powell-Brent (Praxis) non-linear unconstrained optimization. | |
| int | nloptMPSO (NLOPT *nlo, unsigned int maxIter, unsigned int nSwarms, unsigned int nParticles, double wInertia, double wParticle, double wSwarm, double wGlobal, double pDeath, double pImmigration, const int doLocal, TPCSTATUS *) |
| int | nloptITGO1 (NLOPT *nlo, const int doLocal, unsigned int tgoNr, unsigned int sampleNr, unsigned int neighNr, TPCSTATUS *status) |
| int | nloptITGO2 (NLOPT *nlo, const int doLocal, unsigned int tgoNr, unsigned int sampleNr, unsigned int neighNr, TPCSTATUS *status) |
| int | nloptIATGO (NLOPT *nlo, const int doLocal, unsigned int maxIterNr, double neighFract, TPCSTATUS *status) |
Detailed Description
Function Documentation
◆ nlopt1D()
Local one-dimensional minimization by bracketing.
- Precondition
- Initiate the contents of the nlo data structure.
- Returns
- enum tpcerror (TPCERROR_OK when successful).
- See also
- nloptSimplex, nloptPowellBrent
- Parameters
-
nlo Pointer to NLOPT structure.
- Precondition
- Initial guess must be given in x[]. Initial step size must be given in xdelta[]. Constraints xlower[] and xupper[] are required. Parameter tolerances are used as stopping criteria; if set to zero, then only stopping rule is the number of function evaluations.
- Parameters
-
maxIter Maximum number of function evaluations; set to zero to use the default. status Pointer to status data; enter NULL if not needed
Definition at line 25 of file nlopt1d.c.
38 {
40 if(verbose>1) {printf("%s(NLOPT, %d, status)\n", __func__, maxIter); fflush(stdout);}
41 if(nlo==NULL) {
44 }
48 }
49
50 /* Check the number of fixed parameters, based on constraints or delta */
54 }
55 /* Which parameter is the free one? */
57 unsigned int fpi;
61 break;
62 }
66 }
67 if(verbose>3) printf(" index_of_free_parameter := %u\n", fpi);
68
69 /* Check that initial value is inside its limits */
74 if(verbose>2)
79 }
80 /* Check that limits are wider than 2*delta */
81 if(!( (pmax-pmin) > 2.0*delta)) {
84 }
85
86 /* Set the maximum number of iterations, if not given by user */
87 if(maxIter<100) maxIter=200+(pmax-pmin)/delta;
88 if(verbose>1) printf(" maxIter := %u\n", maxIter);
89
90 /* Set 3 bracketing points */
91 double p1=0., p2=0., p3=0., f1=0., f2=0., f3=0.;
92 double minf, minp;
93 /* Set an array function parameters */
96
97 /* Calculate function value with initial parameter value */
98 unsigned int iterNr=1;
101 if(!isfinite(f2)) {
104 }
105 nlo->funCalls++;
109 }
110
111 /* Set the two other bracketing points and compute their function values */
113 nlo->funCalls++;
117 }
118 if(f1<minf) {minf=f1; minp=p1;}
119
121 nlo->funCalls++;
125 }
126 if(f3<minf) {minf=f3; minp=p3;}
127
128
129 /* Now we have 3 points on the function.
130 Start looking for a bracketing set such that f1 > f2 < f3 is the case. */
131 double jump_size=delta;
132 while(!(f1>f2 && f2<f3)) {
133 if(verbose>3) {
134 printf(" jump_size := %g\n", jump_size);
135 printf(" bracketing points: %g %g %g -> %g %g %g\n", p1, p2, p3, f1, f2, f3);
136 printf(" min: %g -> %g\n", minp, minf);
137 }
138 /* check against stopping rules */
139 if(iterNr>=maxIter) break;
140 if((p3-p1)<tol) break;
141 /* if f1 is small then take a step to the left */
142 if(f1<f3) {
143 /* check if the minimum is colliding against the bounds. If so then pick a point between
144 p1 and p2 in the hopes that shrinking the interval will be a good thing to do.
145 Or if p1 and p2 aren't differentiated then try and get them to obtain different values. */
146 if(p1==pmin || (fabs(f1-f2)<macheps && (pmax-pmin)<jump_size )) {
147 p3=p2; f3=f2; p2=0.5*(p1+p2);
149 nlo->funCalls++;
153 }
154 if(f2<minf) {minf=f2; minp=p2;}
155 } else {
156 /* pick a new point to the left of our current bracket */
157 p3=p2; f3=f2; p2=p1; f2=f1;
158 p1-=jump_size; if(p1<pmin) p1=pmin;
160 nlo->funCalls++;
164 }
165 if(f1<minf) {minf=f1; minp=p1;}
166 jump_size*=2.0;
167 }
168 } else { // otherwise f3 is small and we should take a step to the right.
169 /* check if the minimum is colliding against the bounds. If so then pick a point between
170 p2 and p3 in the hopes that shrinking the interval will be a good thing to do.
171 Or if p2 and p3 aren't differentiated then try and get them to obtain different values. */
172 if(p3==pmax || (fabs(f2-f3)<macheps && (pmax-pmin)<jump_size)) {
173 p1=p2; f1=f2; p2=0.5*(p3+p2);
175 nlo->funCalls++;
179 }
180 if(f2<minf) {minf=f2; minp=p2;}
181 } else {
182 /* pick a new point to the right of our current bracket */
183 p1=p2; f1=f2; p2=p3; f2=f3;
184 p3+=jump_size; if(p3>pmax) p3=pmax;
186 nlo->funCalls++;
190 }
191 if(f3<minf) {minf=f3; minp=p3;}
192 jump_size*=2.0;
193 }
194 }
195 }
196 if(verbose>3) {
197 printf(" brackets ready.\n\n");
198 if(verbose>3) {
199 printf(" bracketing points: %g %g %g -> %g %g %g\n", p1, p2, p3, f1, f2, f3);
200 printf(" min: %g -> %g\n", minp, minf);
201 }
202 }
203
204 /* Loop until we have done the max allowable number of iterations or the bracketing window is
205 smaller than required tolerance.
206 Within this loop we maintain the invariant that: f1 > f2 < f3 and p1 < p2 < p3. */
207 const double tau=0.05;
208 while(iterNr<maxIter && (p3-p1)>tol) {
209
210 //p_min = lagrange_poly_min_extrap(p1,p2,p3, f1,f2,f3);
211 double d=f1*(p3*p3-p2*p2) + f2*(p1*p1-p3*p3) + f3*(p2*p2-p1*p1);
212 double e=2.0*(f1*(p3-p2) + f2*(p1-p3) + f3*(p2-p1));
213 if(fabs(e)<macheps || !isfinite(d/=e)) {
214 minp=p2;
215 } else {
216 if(p1<=d && d<=p3) {minp=d;} else {minp=d; if(p1>minp) minp=p1; if(p3<minp) minp=p3;}
217 }
218
219 /* make sure min p isn't too close to the three points we already have */
220 if(minp<p2) {
221 d=(p2-p1)*tau; if(fabs(p1-minp)<d) minp=p1+d; else if(fabs(p2-minp)<d) minp=p2-d;
222 } else {
223 d=(p3-p2)*tau; if(fabs(p2-minp)<d) minp=p2+d; else if(fabs(p3-minp)<d) minp=p3-d;
224 }
225
226 /* make sure one side of the bracket isn't super huge compared to the other side.
227 If it is then contract it. */
228 double bracket_ratio=fabs(p1-p2)/fabs(p2-p3);
229 if(!(bracket_ratio<100.0 && bracket_ratio>0.01)) {
230 /* Force min p to be on a reasonable side. */
231 if(bracket_ratio>1.0 && minp>p2) minp=0.5*(p1+p2); else if(minp<p2) minp=0.5*(p2+p3);
232 }
233
234 /* Compute function value at min p */
236 nlo->funCalls++;
240 }
241
242 /* Remove one of the endpoints of our bracket depending on where the new point falls */
243 if(minp<p2) {
244 if(f1>minf && minf<f2) {p3=p2; f3=f2; p2=minp; f2=minf;} else {p1=minp; f1=minf;}
245 } else {
246 if(f2>minf && minf<f3) {p1=p2; f1=f2; p2=minp; f2=minf;} else {p3=minp; f3=minf;}
247 }
248
249 if(verbose>3) {
250 printf(" bracketing points: %g %g %g -> %g %g %g\n", p1, p2, p3, f1, f2, f3);
251 printf(" min: %g -> %g\n", minp, minf);
252 }
253 }
254
255 /* Copy the estimated parameter into nlo structure */
257 /* Compute function value at the minimum, in case user uses the fitted data in nlo structure */
259
262}
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
Referenced by nloptSimplex().
◆ nloptAddP()
|
extern |
Add parameters and value in NLOPT plist.
Parameters are added starting at index plist[(funCalls-1)*(totalNr+1)], and function value at index plist[(funCalls-1)*(totalNr+1)+totalNr], after plist[] size is increased by (totalNr+1).
- See also
- nloptInit, nloptFree, nloptAllocate, nloptSortP, nloptMeanP
- Returns
- Returns TPCERROR status, TPCERROR_OK (0) when successful.
- Parameters
-
nlo Pointer to initiated NLOPT structure; any old contents are deleted. p Pointer to parameter array of length nlo->totalNr. funval Value of objective function.
Definition at line 149 of file nlopt.c.
156 {
160 /* Allocate memory, if this is the first list item */
162 } else {
163 /* Reallocate memory, if this is not the first list item */
165 }
167 /* Copy the values */
172}
Referenced by nlopt1D(), nloptIATGO(), and nloptSimplex().
◆ nloptAllocate()
|
extern |
Allocate memory for NLOPT data. Any previous contents are deleted.
- See also
- nloptInit, nloptFree, nloptDuplicate
- Returns
- Returns TPCERROR status, TPCERROR_OK (0) when successful.
- Parameters
-
nlo Pointer to initiated NLOPT structure; any old contents are deleted. parNr Nr of parameters.
Definition at line 74 of file nlopt.c.
79 {
81 /* Delete any previous contents */
82 nloptFree(nlo);
83 /* If no memory is requested, then just return */
85
86 /* Allocate memory for arrays */
97 /* Set array length */
98 nlo->totalNr=parNr;
100}
Referenced by nloptDuplicate().
◆ nloptCheckParameters()
|
extern |
Check that model parameters are within given limits.
If one or more parameter(s) are outside of limits, then a penalty factor is optionally computed. Optionally, this function can make a list parameters that are inside the constraints.
- Returns
- Return the number of parameters that are inside or at the limits; 0 in case of error or if all are outside of limits.
- See also
- nloptForceLimits, aicSS
- Parameters
-
n Nr of parameters. pLower Lower limits (array of length par_nr). pUpper Upper limits (array of length par_nr). p Parameters to test (array of length par_nr). pAccept Pointer to corrected parameters (array of length par_nr); NULL if not needed. penalty Pointer to variable in which the possible penalty factor will be written; 1.0 if no penalty was found is necessary, otherwise >1. Set to NULL if not needed.
Definition at line 30 of file constraints.c.
44 {
45 unsigned int pi, nAccept=0;
46 double range;
47
48 if(penalty!=NULL) *penalty=1.0;
49 if(n<1 || p==NULL || pLower==NULL || pUpper==NULL)
50 return(nAccept);
51 for(pi=0; pi<n; pi++) {
52 range=pUpper[pi]-pLower[pi]; if(range<1.0E-30) range=1.0;
53 if(p[pi]<pLower[pi]) {
54 if(pAccept!=NULL) pAccept[pi]=pLower[pi];
55 if(penalty!=NULL) *penalty += (pLower[pi]-p[pi])/range;
56 } else if(p[pi]>pUpper[pi]) {
57 if(pAccept!=NULL) pAccept[pi]=pUpper[pi];
58 if(penalty!=NULL) *penalty += (p[pi]-pUpper[pi])/range;
59 } else {
60 if(pAccept!=NULL) pAccept[pi]=p[pi];
61 nAccept++;
62 }
63 }
64 return nAccept;
65}
◆ nloptdataInit()
|
extern |
Initiate the NLOPT_DATA structure before any use.
- See also
- nloptAllocate, nloptFree, nloptDuplicate
- Parameters
-
d Pointer to NLOPT_DATA.
◆ nloptDuplicate()
Make a duplicate of NLOPT data.
NLOPT plist is not copied but usePList parameter is copied.
- Postcondition
- After the duplicate is no more needed, free the memory using nloptFree().
- See also
- nloptAllocate, nloptFree, nloptInit
- Returns
- Returns TPCERROR status, TPCERROR_OK (0) when successful.
- Parameters
-
nlo1 Pointer to existing NLOPT structure to be duplicated. nlo2 Pointer to the target NLOPT; must be initiated; any old contents are deleted.
Definition at line 112 of file nlopt.c.
117 {
120
123
137}
Referenced by nloptITGO1(), and nloptITGO2().
◆ nloptFixedNr()
|
extern |
Determine the number of fixed parameters, based on constraints or delta.
- See also
- nloptLimitFixedNr, nloptInit, nloptFree, nloptAllocate
- Returns
- Returns the number of fixed parameters.
- Parameters
-
d Pointer to NLOPT
Definition at line 354 of file nlopt.c.
357 {
359 unsigned int i, ln=0, dn=0;
361 /* based on limits */
363 double r;
366 if(!isnan(r) && r<macheps) ln++;
367 }
368 }
369 /* based on deltas */
372//printf(" %g %g\n", d->xdelta[i], macheps);
374 }
375 }
376//printf(" %u %u\n", ln, dn);
377 /* the one that is nonzero is probably used */
378 if(dn==0) return(ln);
379 if(ln==0) return(dn);
380 /* Both seem to be used; thus check both together */
381 unsigned int n=0;
382 double r;
386 if(!isnan(r) && r<macheps) {n++; continue;}
387 }
388 return(n);
389}
Referenced by nlopt1D(), and nloptSimplex().
◆ nloptForceLimits()
|
extern |
Enforce the model parameters within given limits.
- Returns
- Return the number of parameters that are inside or at the limits; 0 in case of error or if all are outside of limits.
- See also
- nloptCheckParameters
- Parameters
-
n Number of parameters. pLower Lower limits (array of length par_nr). pUpper Upper limits (array of length par_nr). p Parameters to enforce inside limits (array of length par_nr).
Definition at line 76 of file constraints.c.
85 {
86 unsigned int pi, nAccept=0;
87
88 if(n<1 || p==NULL || pLower==NULL || pUpper==NULL) return(nAccept);
89 for(pi=0; pi<n; pi++) {
90 if(!(p[pi]>=pLower[pi])) {
91 p[pi]=pLower[pi];
92 } else if(!(p[pi]<=pUpper[pi])) {
93 p[pi]=pUpper[pi];
94 } else {
95 nAccept++;
96 }
97 }
98 return(nAccept);
99}
Referenced by nloptSimplex().
◆ nloptFree()
|
extern |
Free memory allocated for NLOPT data. All contents are destroyed.
- See also
- nloptInit, nloptAllocate, nloptDuplicate
- Precondition
- Before first use initialize the structure with nloptInit().
- Parameters
-
nlo Pointer to initiated NLOPT structure.
Definition at line 52 of file nlopt.c.
55 {
56 if(nlo==NULL) return;
57 free(nlo->xfull);
58 free(nlo->xlower);
59 free(nlo->xupper);
60 free(nlo->xdelta);
61 free(nlo->xtol);
62 free(nlo->plist);
63 // then set everything to zero or NULL again
64 nloptInit(nlo);
65}
Referenced by nloptAddP(), nloptAllocate(), nloptDuplicate(), nloptITGO1(), and nloptITGO2().
◆ nloptGaussianPoint()
|
extern |
Create random parameters with Gaussian distribution.
- Precondition
- Uses rand(), therefore set seed for a new series of pseudo-random numbers; to produce truly random numbers (not just pseudo-random), do srand(time(0)) before calling this function. If no seed is set, then value 1 is used as default seed by rand().
- Returns
- Returns non-zero in case of an error.
- Parameters
-
p Pointer to parameter list to be filled, allocated for length n. mean Pointer to list of mean values for each parameter. sd Pointer to list of standard deviations for each parameter. low Pointer to parameter lower limits; generated pseudo-random numbers falling below the lower limit will be set to the limit; NULL if not needed. up Pointer to parameter upper limits; generated pseudo-random numbers falling above the upper limit will be set to the limit; NULL if not needed. n List length. mt Pointer to initiated and seeded Mersenne Twister MT19937 data structure; enter NULL to use drandGaussian() instead.
- See also
- mertwiInit, mertwiInitWithSeed64
Definition at line 70 of file rndpoint.c.
89 {
90 if(p==NULL || mean==NULL || sd==NULL) return(1);
91 if(n==0) return(0);
92
93 for(unsigned int i=0; i<n; i++) {
96 if(low!=NULL && p[i]<low[i]) p[i]=low[i];
97 if(up!=NULL && p[i]>up[i]) p[i]=up[i];
98 }
99
100 return(0);
101}
double drandGaussian()
Get pseudo-random number with normal (Gaussian) distribution with mean 0 and SD 1.
Definition gaussdev.c:144
double mertwiRandomGaussian(MERTWI *mt)
Generate a 64-bit double precision floating point pseudo-random number with normal (Gaussian) distrib...
Definition mertwi.c:354
Referenced by nloptIATGO(), nloptMPSO(), and nloptSimplexARRS().
◆ nloptIATGO()
|
extern |
Iterative accumulative topographical global optimization algorithm (IATGO).
- Precondition
- Uses drand(), therefore set seed for a new series of pseudo-random numbers with drandSeed(1) before calling this function.
- 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).
- See also
- nloptMPSO, nloptSimplexARRS
- Parameters
-
nlo Pointer to NLOPT data. Counter funCalls is initially set to zero, and then increased here.
- Precondition
- Constraints xlower[] and xupper[] are required. Parameter maxFunCalls is used as one of the stopping criteria. Parameters xtol[] is used for scaling and as one of the stopping criteria. Initial guess can be given in xfull[], but it is not obligatory. Initial step sizes (xdelta[]) are not used.
- Parameters
-
doLocal Local optimization method: 0=no local optimization, 1=Nelder-Mead. maxIterNr Number of TGO iterations; enter 0 to use the default; enter 1 to run TGO just once, ie to use TGO instead of iTGO. neighFract Fraction of samples to define as neighbours (cluster size); enter a large fraction (>0.5) if only few distant local minima are expected. status Pointer to status data; enter NULL if not needed.
Definition at line 681 of file tgo.c.
701 {
702 FILE *fp=stdout;
704 if(verbose>1) {
705 fprintf(fp, "%s(NLOPT, %d, %u, %g, status)\n", __func__, doLocal, maxIterNr, neighFract);
706 fflush(fp);
707 }
708 if(nlo==NULL ) {
711 }
715 }
716
718
719 if(verbose>4) {
720 fprintf(fp, "\nInitial limits and tolerances:\n");
721 for(unsigned int i=0; i<dim; i++)
723 }
724
725 /* Check if any of the parameters are fixed */
727 if(verbose>2 && fixedParNr>0) fprintf(fp, "fixedParNr := %d\n", fixedParNr);
729 if(verbose>2) {fprintf(fp, "fittedParNr := %d\n", fittedParNr); fflush(fp);}
730 if(fittedParNr<1) {
733 }
734
735 /* Check the tolerations */
736 for(unsigned int i=0; i<dim; i++) {
739 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
742 }
743 }
744
745 /* Continue input checking */
746 if(neighFract<0.04) neighFract=0.04; else if(neighFract>0.97) neighFract=0.97;
747 unsigned int sampleNr=50*fittedParNr/neighFract;
748 if(maxIterNr==0) maxIterNr=5+2*fittedParNr;
749 unsigned int neighNr=neighFract*sampleNr;
750 if(verbose>2) {
751 fprintf(fp, "sampleNr := %u\n", sampleNr);
752 fprintf(fp, "neighFract := %g\n", neighFract);
753 fprintf(fp, "neighNr := %u\n", neighNr);
754 fprintf(fp, "maxIterNr := %u\n", maxIterNr);
755 fflush(stdout);
756 }
757
758
759 /*
760 * Initialize the sample list with random points
761 */
762 if(verbose>3) {fprintf(fp, "Making initial random samples\n"); fflush(fp);}
763
764 /* Allocate memory */
766 nlo->funCalls=0;
771 }
772
773 /* Fill the list with random points and compute the object function values */
774 {
775 unsigned int badNr=0;
776 /* Add the random points */
777 for(unsigned int si=0; si<sampleNr; si++) {
781 }
783 nlo->funCalls++;
785 }
786 if(verbose>3 && badNr>0) fprintf(fp, "Number of bad points: %d\n", badNr);
787 /* Not all object function values can be bad */
788 if((sampleNr-badNr)<10) {
789 if(verbose>0) {
790 fprintf(stderr, "Error: invalid values inside parameter range.\n"); fflush(stderr);}
793 }
794 /* Try new points to replace the failed ones */
795 if(badNr>0) {
796 badNr=0;
800 //nlo->funCalls++;
802 }
803 if(verbose>4 && badNr>0) fprintf(fp, "Number of bad points after retry: %d\n", badNr);
804 }
805 }
806
807 /* If initial guess is valid, add it to the list */
808 int initGuessAvailable=1;
809 double initGuessCost=nan("");
810 for(unsigned int i=0; i<dim; i++) {
814 }
815 if(initGuessAvailable) {
817 if(isfinite(initGuessCost)) {
818 nlo->funCalls++;
820 } else {
821 initGuessAvailable=0;
822 }
823 }
824 if(verbose>2) {
825 if(initGuessAvailable) fprintf(fp, "valid initial guess with cost=%g\n", initGuessCost);
826 else fprintf(fp, "no valid initial guess provided.\n");
827 }
828
829
830 /*
831 * Start iterations
832 */
833 if(verbose>2) {fprintf(fp, "\nStarting TGO iterations\n"); fflush(fp);}
834 unsigned int iterNr=0, stopCounter=0;
835 double prevBest=nan("");
836 while(iterNr<maxIterNr && nlo->funCalls<nlo->maxFunCalls) {
837 iterNr++;
839 if(verbose>3) {fprintf(fp, "\nIteration %d with %d samples\n", iterNr, sampleNr); fflush(fp);}
840
841 /* Sort samples based on the evaluated function value */
844
845 /* Print sampled points */
848 else if(verbose>4) {
849 fprintf(fp, "best point so far:");
852 }
853
854 /* If SD of the best points is minimal, and fit is not improving, then stop */
855 if(verbose>4 && !isnan(prevBest))
858 if(verbose>4) fprintf(fp, "Best point did not improve.\n");
860 unsigned int i=0;
862 if(i==dim) {
863 if(verbose>2) fprintf(fp, "Small SD of best points.\n");
864 stopCounter++;
865 } else stopCounter=0;
866 }
867 } else stopCounter=0;
868 prevBest=nlo->plist[dim];
869
870 /* Stop, if stopping rules bumped into several times */
871 if(verbose>4) fprintf(fp, "stopCounter := %d\n", stopCounter);
872 if(stopCounter>2) {
873 if(verbose>2) fprintf(fp, "\n Required tolerance reached.\n");
874 break;
875 }
876
877 /* Find topographic minima (TM) */
878
879 /* Allocate a list specifying whether point is local minimum (1) or not (0) */
880 /* Allocate a list specifying whether point is part of topographic minimum (>0) or not (0) */
881 unsigned int tmlist[sampleNr];
882 for(unsigned int i=0; i<sampleNr; i++) tmlist[i]=0;
883
884 unsigned int topoNr=0;
885 if(verbose>3) {fprintf(fp, "neighNr := %d\n", neighNr); fflush(fp);}
886
887 /* For each point i, find out if it is a "topographic minimum", that is,
888 better than k neighbour points. */
889 for(unsigned int si=0; si<sampleNr-neighNr; si++) if(tmlist[si]==0) {
890
891 /* If function value is invalid, then this cannot be accepted as TM */
893
894 /* Compute the distances to other points, scaled using xtol[] */
895 double sdist[sampleNr];
896 for(unsigned int sj=0; sj<sampleNr; sj++) {
898 sdist[sj]=nan(""); continue;} // valid function value required
899 sdist[sj]=0.0;
900 for(unsigned int k=0; k<dim; k++) {
903 /*if(isfinite(d))*/ sdist[sj]+=d*d;
904 }
905 /* If scaled distance is less than 1, the points are practically equal: leave out */
906 //if(sdist[sj]<1.0) sdist[sj]=nan(""); // would lead to multiple local minima at same point
907 }
908 double sdist2[sampleNr]; // make copy for printing
909 for(unsigned int sj=0; sj<sampleNr; sj++) sdist2[sj]=sdist[sj];
910
911 /* Find its closest neighbours that are not any better */
912 unsigned int nn=0; double prev_mind=nan("");
913 while(1) {
914 /* Find the closest neighbour that was not already found */
915 unsigned int mini=0;
916 double mind=nan("");
917 for(unsigned int sj=0; sj<sampleNr; sj++) {
918 if(!isfinite(mind) || mind>sdist[sj]) {mind=sdist[sj]; mini=sj;}
919 }
920 if(!isfinite(mind)) break;
921 sdist[mini]=nan(""); // this point not used again in search of neighbours
922 if(verbose>100) fprintf(fp, " min_distance[%u] := %g , with fval := %g\n",
923 1+nn, mind, nlo->plist[mini*(dim+1)+dim]);
924 /* If this neighbour is better than the tentative TM, then stop */
926 /* Otherwise, count this as a worse neighbour, belonging to this TM */
927 nn++; tmlist[mini]=1+si;
928 // we do not want to include more than neighNr points into this TM,
929 // but if distance is less than one or not larger than previously then
930 // this is essentially the same point
931 if(isnan(prev_mind)) prev_mind=mind;
932 if(nn>=neighNr && mind>1.0 && mind>prev_mind) break;
933 prev_mind=mind;
934 }
935 /* If the number of worse neighbours is less than required, then stop considering this
936 as a local optimum, and continue with the next sample */
937 if(nn<neighNr) {
938 /* its neighbours then do not belong to TM either */
939 for(unsigned int ni=0; ni<sampleNr; ni++) if(tmlist[ni]==1+si) tmlist[ni]=0;
940 continue;
941 }
942 /* otherwise mark this as topographic minimum (TM) */
943 tmlist[si]=1+si;
944 topoNr++;
945 /* Print TM point */
946 if(verbose>5) {
947 fprintf(fp, "TM point %d\t", 1+si);
950 if(verbose>14) {
951 fprintf(fp, "and its %u neighbours:\n", nn);
952 for(unsigned int sj=0; sj<sampleNr; sj++) if(sj!=si && tmlist[sj]==1+si) {
955 }
956 }
957 fflush(fp);
958 }
959 }
960 if(verbose>3) {fprintf(fp, "topographic minima: %d\n", topoNr); fflush(fp);}
961 if(topoNr==0) {
962 if(verbose>0) {fprintf(stderr, "Error: no topographical minima found\n"); fflush(stderr);}
965 }
966
967 /*
968 * Local optimization from each TM
969 */
970 for(unsigned int si=0; si<sampleNr; si++) if(tmlist[si]==1+si) {
971 if(verbose>5) {fprintf(fp, "LO for TM point %d\n", 1+si); fflush(fp);}
972 if(verbose>6) {
973 fprintf(fp, "TM point %d before LO:\n", 1+si);
976 }
977 /* Compute the mean and SD of points belonging to this TM */
978 double meanp[dim], sdp[dim]; unsigned int nn=0;
979 for(unsigned int i=0; i<dim; i++) {
981 double x[sampleNr];
982 nn=0;
983 for(unsigned int sj=0; sj<sampleNr; sj++)
985 statMeanSD(x, nn, meanp+i, sdp+i, NULL);
986 } else {
987 meanp[i]=nlo->xlower[i]; sdp[i]=0.0;
988 }
989 }
990 if(verbose>7) {
991 fprintf(fp, "TM point (n=%u) mean and SD:\n", nn);
992 for(unsigned int i=0; i<dim; i++) fprintf(fp, "\t%g\t%g\n", meanp[i], sdp[i]);
993 fflush(fp);
994 }
995 /* Make sure that SD is not zero, unless parameter is fixed */
998 }
999 if(doLocal==1) { // downhill simplex from the best point so far
1005 } else {
1006 if(verbose>6) {
1007 fprintf(fp, "Point after LO:");
1010 }
1011 }
1012 } else { // Add random points with Gaussian distribution
1013 unsigned int bi=0; double bval=nan("");
1014 for(unsigned int ni=0; ni<neighNr; ni++) {
1015 double p[dim];
1020 }
1021 if(verbose>6) {
1022 fprintf(fp, "Best TM point %d after LO:\n", 1+si);
1025 }
1026 }
1027 } // next TM
1028
1029 } // end of iterations
1030
1031 /* Check the reason for loop exit */
1032 if(verbose>1) {
1033 if(iterNr>=maxIterNr)
1034 fprintf(fp, "\n Exceeded the maximum number for loops.\n");
1036 fprintf(fp, "\n Exceeded the maximum number for function calls.\n");
1037 fflush(fp);
1038 }
1039
1040 /* Sort samples based on the evaluated function value */
1043 /* Get the best point so far */
1046
1047
1050}
double drandExponential(double mean)
Get pseudo-random number with exponential distribution.
Definition gaussdev.c:178
int statMeanSD(double *data, unsigned int n, double *mean, double *sd, unsigned int *vn)
Definition mean.c:25
int nloptMeanP(NLOPT *nlo, unsigned int nr, double *meanp, double *sdp)
Definition nlopt.c:219
int nloptRandomPoint(double *p, double *low, double *up, unsigned int n, MERTWI *mt)
Definition rndpoint.c:28
int nloptGaussianPoint(double *p, double *mean, double *sd, double *low, double *up, unsigned int n, MERTWI *mt)
Definition rndpoint.c:70
int nloptSimplex(NLOPT *nlo, unsigned int maxIter, TPCSTATUS *status)
Definition simplex.c:32
FILE * fp
File pointer for writing log information during development and testing.
Definition tpcextensions.h:243
◆ nloptInit()
|
extern |
Initiate the NLOPT structure before any use.
- See also
- nloptAllocate, nloptFree, nloptDuplicate
- Parameters
-
nlo Pointer to NLOPT
Definition at line 25 of file nlopt.c.
28 {
29 if(nlo==NULL) return;
30 nlo->totalNr=0;
31 nlo->xfull=NULL;
32 nlo->xlower=NULL;
33 nlo->xupper=NULL;
34 nlo->xdelta=NULL;
35 nlo->xtol=NULL;
36 nlo->_fun=NULL;
37 nlo->fundata=NULL;
38 nlo->maxFunCalls=0;
39 nlo->funCalls=0;
41 nlo->usePList=0;
42 nlo->plist=NULL;
43}
Referenced by nloptFree(), nloptITGO1(), and nloptITGO2().
◆ nloptITGO1()
|
extern |
Iterative Topographical global optimization algorithm 1 (iTGO1).
- Remarks
- Not well tested!
- Precondition
- Uses drand(), therefore set seed for a new series of pseudo-random numbers with drandSeed(1) before calling this function.
- 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).
- Todo
- Finalize the code and tests.
- See also
- nloptIATGO, nloptMPSO, nloptSimplexARRS, nloptITGO2
- Parameters
-
nlo Pointer to NLOPT data. Counter funCalls is initially set to zero, and then increased here. Parameter maxFunCalls is used as one of the stopping criteria. Parameters xtol[] is used for scaling and as one of the stopping criteria. doLocal Local optimization method: 0=no local optimization, 1=Nelder-Mead. tgoNr Number of TGO iterations; enter 0 to use the default; enter 1 to run TGO just once, ie to use TGO instead of iTGO. Iterations are needed if sampleNr (below) would otherwise require too much memory. sampleNr Number of points to sample in one iteration; may need to be large if the number of iterations (above) is small; enter 0 to let function decide it. neighNr Number of neighbours to investigate (cluster size); enter a large number relative to sampleNr (above) if only few local minima are expected; enter 0 to let function decide it. status Pointer to status data; enter NULL if not needed.
Definition at line 59 of file tgo.c.
81 {
83 if(verbose>0) {
84 printf("%s(NLOPT, %d, %u, %u, %u, status)\n", __func__, doLocal, tgoNr, sampleNr, neighNr);
85 fflush(stdout);
86 }
87 if(nlo==NULL ) {
90 }
94 }
96 if(verbose>0) {fprintf(stderr, "Error: too many dimensions to optimize.\n"); fflush(stderr);}
99 }
100
101 /* Check if any of the parameters are fixed */
103 if(verbose>1 && fixedParNr>0) printf("fixedParNr := %d\n", fixedParNr);
105 if(verbose>1) printf("fittedParNr := %d\n", fittedParNr);
106 if(fittedParNr<1) {
109 }
110
111 /* Check the tolerations */
115 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
118 }
119 }
120
121 /* Continue input checking */
122 if(sampleNr<10) sampleNr=50*fittedParNr;
123 if(sampleNr&1) sampleNr++; // If number is odd, then add 1
124 if(neighNr<2) neighNr=sampleNr/2;
125 else if(neighNr>sampleNr-1) neighNr=sampleNr-1; // Make sure "neighNr" isn't too big
126 if(tgoNr==0) tgoNr=1+fittedParNr;
127 if(verbose>2) {
128 printf("sampleNr := %d\n", sampleNr);
129 printf("neighNr := %d\n", neighNr);
130 printf("tgoNr := %d\n", tgoNr);
131 fflush(stdout);
132 }
133
134 /* Check if initial guess is valid */
135 int initGuessAvailable=1;
136 double initGuessCost=nan("");
141 }
142 if(initGuessAvailable) {
144 nlo->funCalls++;
145 if(!isfinite(initGuessCost)) initGuessAvailable=0;
146 else if(verbose>2) {
147 printf("valid initial guess with cost=%g\n", initGuessCost);
148 }
149 }
150
151 /* Allocate memory */
152 TGO_POINT *points=(TGO_POINT*)malloc(sampleNr*sizeof(TGO_POINT));
153 if(points==NULL) {
156 }
157 for(unsigned int i=0; i<sampleNr; i++) {
158 points[i].topomin=0;
159 points[i].fvalue=0.0;
160 }
161
162 /*
163 * Iterative TGO, or non-iterative if tgoNr==1
164 */
165 unsigned int topoNr=0;
166 for(unsigned int tgoi=0; tgoi<tgoNr; tgoi++) {
167
168 if(verbose>3 && tgoNr>1) {printf("TGO iteration # %d: \n", 1+tgoi); fflush(stdout);}
169
170 /* Sample N points in the feasible region and compute the object function values for
171 those points which do not already have it. */
172 unsigned int badNr=0;
173 for(unsigned int si=0; si<sampleNr; si++) if(points[si].topomin==0) {
174 if(tgoi==0 && si==0 && initGuessAvailable) {
175 /* If initial guess was given, then use it as the first point */
177 points[si].fvalue=initGuessCost;
178 continue;
179 }
182 nlo->funCalls++;
183 /* If function return value was not valid, then we'll try twice with new guesses */
184 int tries=0;
185 while(!isfinite(points[si].fvalue) && tries<2) {
186 if(verbose>8) {
187 printf("this point did not give normal return value:\n");
189 printf("\n");
190 }
193 nlo->funCalls++;
194 tries++;
195 }
196 if(!isfinite(points[si].fvalue)) badNr++;
197 }
198 if(verbose>4 && badNr>0) printf("Number of bad points: %d\n", badNr);
199 /* Object functions values must be good for at least NeighNr points */
200 if(tgoi==0 && (sampleNr-badNr)<=neighNr) {
201 if(verbose>0) {
202 fprintf(stderr, "Error: invalid values inside parameter range.\n"); fflush(stderr);}
203 free(points);
206 }
207 /* Print sampled points */
208 if(verbose>10) {
209 printf("Sampled points:\n");
210 for(unsigned int si=0; si<sampleNr; si++) {
211 printf("%d\t", 1+si);
213 printf("=> %e\n", points[si].fvalue);
214 }
215 fflush(stdout);
216 }
217
218 /* For each point i, find out if it is a "topographic minimum", that is,
219 better than k neighbour points. */
220 topoNr=0;
221 for(unsigned int si=0; si<sampleNr; si++) {
222 points[si].topomin=0;
223
224 /* Compute the distances, scaled using xtol[] */
225 double sdist[sampleNr];
226 for(unsigned int sj=0; sj<sampleNr; sj++) {
227 sdist[sj]=0.0;
228 if(si==sj) {sdist[sj]=1.0E+99; continue;}
232 if(isfinite(d)) sdist[sj]+=d*d;
233 }
234 /* Distance should be computed as square root, but it does not affect
235 the order of distances; since sqrt() is slow, it is not done. */
236 // sdist[sj]=sqrt(sdist[sj]);
237 }
238
239 /* Find the closest neighbours for sample i */
240 unsigned int neighs[neighNr]; // collect here the neighbours
241 unsigned int ni=0;
242 for(ni=0; ni<neighNr; ni++) {
243 unsigned int mini=0;
244 double minv=sdist[mini];
245 for(unsigned int sj=0; sj<sampleNr; sj++) {
246 if(sdist[sj]<minv) {minv=sdist[sj]; mini=sj;}
247 }
248 sdist[mini]=1.0E+99; // not used again
249 /* If point i is worse than any of the closest neighbours, then point i
250 is certainly not a topographic minimum; then stop this loop and go to
251 the next point i+1 */
252 if(!(points[si].fvalue<points[mini].fvalue)) break;
253 neighs[ni]=mini;
254 }
255 if(ni<neighNr) continue; // continue with the next i
256 /* otherwise mark this as topographic minimum (TM) */
257 points[si].topomin=1; topoNr++;
258 /* Print TM point */
259 if(verbose>6) {
260 printf("TM point %d:\n", topoNr);
261 printf("%d\t", 1+si);
263 printf("=> %e\n", points[si].fvalue);
264 if(verbose>7) {
265 printf(" points neighbours:");
266 for(ni=0; ni<neighNr; ni++) printf(" %d", 1+neighs[ni]);
267 printf("\n");
268 }
269 fflush(stdout);
270 }
271
272 }
273 if(verbose>4) {printf(" %d topographical minima\n", topoNr); fflush(stdout);}
274 if(topoNr==0) {
275 if(verbose>0) {fprintf(stderr, "Warning: no topographical minima found\n"); fflush(stderr);}
276 continue;
277 }
278
279 } // end of iTGO
280
281
282 /* If requested, do local optimization for the TMs */
283 if(topoNr==0) {
284 if(verbose>0) {fprintf(stderr, "Error: no topographical minima found\n"); fflush(stderr);}
285 free(points);
288 }
289 for(unsigned int si=0; si<sampleNr; si++) if(points[si].topomin>0) {
290 if(verbose>2) printf("LO for TM point %d\n", 1+si);
292 nloptDuplicate(nlo, &lnlo);
294 lnlo.funval=points[si].fvalue;
300 points[si].fvalue=lnlo.funval;
301 if(verbose>6) {
302 printf("TM point %d after LO:\n", 1+si);
303 printf("%d\t", 1+si);
305 printf("=> %e\n", points[si].fvalue); fflush(stdout);
306 }
307 } else if(verbose>7) {printf(" LO failed\n"); fflush(stdout);}
309 nloptFree(&lnlo);
310 }
311
312 /* Find the best minimum and return it */
313 {
314 unsigned int mini=0;
315 double minv=points[mini].fvalue;
316 for(unsigned int si=1; si<sampleNr; si++)
317 if(!(points[si].fvalue>minv)) {minv=points[si].fvalue; mini=si;}
319 nlo->funval=points[mini].fvalue;
320 }
321
322 free(points);
325}
Definition tpcnlopt.h:25
◆ nloptITGO2()
|
extern |
Iterative Topographical global optimization algorithm 2 (iTGO2).
- Remarks
- Not well tested!
- Precondition
- Uses drand(), therefore set seed for a new series of pseudo-random numbers with drandSeed(1) before calling this function.
- 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).
- Todo
- Finalize the code and tests.
- See also
- nloptIATGO, nloptMPSO, nloptSimplexARRS, nloptITGO1
- Parameters
-
nlo Pointer to NLOPT data. Counter funCalls is initially set to zero, and then increased here. Parameter maxFunCalls is used as one of the stopping criteria. Parameters xtol[] is used for scaling and as one of the stopping criteria. doLocal Local optimization method: 0=no local optimization, 1=Nelder-Mead. tgoNr Number of TGO iterations; enter 0 to use the default; enter 1 to run TGO just once, ie to use TGO instead of iTGO. Iterations are needed if sampleNr (below) would otherwise require too much memory. sampleNr Number of points to sample in one iteration; may need to be large if the number of iterations (above) is small; enter 0 to let function decide it. neighNr Number of neighbours to investigate (cluster size); enter a large number relative to sampleNr (above) if only few local minima are expected; enter 0 to let function decide it. status Pointer to status data; enter NULL if not needed.
Definition at line 342 of file tgo.c.
364 {
366 if(verbose>0) {
367 printf("%s(NLOPT, %d, %u, %u, %u, status)\n", __func__, doLocal, tgoNr, sampleNr, neighNr);
368 fflush(stdout);
369 }
370 if(nlo==NULL ) {
373 }
377 }
379 if(verbose>0) {fprintf(stderr, "Error: too many dimensions to optimize.\n"); fflush(stderr);}
382 }
383
384 /* Check if any of the parameters are fixed */
386 if(verbose>1 && fixedParNr>0) printf("fixedParNr := %d\n", fixedParNr);
388 if(verbose>1) printf("fittedParNr := %d\n", fittedParNr);
389 if(fittedParNr<1) {
392 }
393
394 /* Check the tolerations */
398 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
401 }
402 }
403
404 /* Continue input checking */
405 if(sampleNr<10) sampleNr=50*fittedParNr;
406 if(sampleNr&1) sampleNr++; // If number is odd, then add 1
407 if(neighNr<2) neighNr=sampleNr/2;
408 else if(neighNr>sampleNr-1) neighNr=sampleNr-1; // Make sure "neighNr" isn't too big
409 if(tgoNr==0) tgoNr=1+fittedParNr;
410 if(verbose>2) {
411 printf("sampleNr := %d\n", sampleNr);
412 printf("neighNr := %d\n", neighNr);
413 printf("tgoNr := %d\n", tgoNr);
414 fflush(stdout);
415 }
416
417 /* Check if initial guess is valid */
418 int initGuessAvailable=1;
419 double initGuessCost=nan("");
424 }
425 if(initGuessAvailable) {
427 nlo->funCalls++;
428 if(!isfinite(initGuessCost)) initGuessAvailable=0;
429 else if(verbose>2) {
430 printf("valid initial guess with cost=%g\n", initGuessCost);
431 }
432 }
433
434 /* Allocate memory */
435 TGO_POINT *points=(TGO_POINT*)malloc(sampleNr*sizeof(TGO_POINT));
436 if(points==NULL) {
439 }
440 for(unsigned int i=0; i<sampleNr; i++) {
441 points[i].topomin=0;
442 points[i].fvalue=0.0;
443 }
444 TGO_POINT *gpoints=(TGO_POINT*)malloc(sampleNr*sizeof(TGO_POINT));
445 if(gpoints==NULL) {
446 free(points);
449 }
450 for(unsigned int i=0; i<sampleNr; i++) {
451 gpoints[i].topomin=0;
452 gpoints[i].fvalue=0.0;
453 }
454
455 /*
456 * Iterative TGO, or non-iterative if tgoNr==1
457 */
458 unsigned int topoNr=0;
459 if(initGuessAvailable) {
460 /* If initial guess was given, then use it as the first G point */
462 gpoints[0].fvalue=initGuessCost;
463 topoNr=1;
464 }
465 for(unsigned int tgoi=0; tgoi<tgoNr; tgoi++) {
466
467 if(verbose>2 && tgoNr>1) {printf("TGO iteration # %d: \n", 1+tgoi); fflush(stdout);}
468
469 while(topoNr<sampleNr) {
470 if(verbose>3) {printf(" topoNr := %d\n", topoNr); fflush(stdout);}
471
472 /* Sample N points and compute the object function values for all points. */
473 unsigned int badNr=0;
474 for(unsigned int si=0; si<sampleNr; si++) {
475 points[si].topomin=0;
478 nlo->funCalls++;
479 /* If function return value was not valid, then we'll try twice with new guesses */
480 int tries=0;
481 while(!isfinite(points[si].fvalue) && tries<2) {
482 if(verbose>8) {
483 printf("this point did not give normal return value:\n");
485 printf("\n");
486 }
489 nlo->funCalls++;
490 tries++;
491 }
492 if(!isfinite(points[si].fvalue)) badNr++;
493 }
494 if(verbose>4 && badNr>0) printf("Number of bad points: %d\n", badNr);
495 /* Object functions values must be good for at least NeighNr points */
496 if(tgoi==0 && (sampleNr-badNr)<=neighNr) {
497 if(verbose>0) {
498 fprintf(stderr, "Error: invalid values inside parameter range.\n"); fflush(stderr);}
499 free(points); free(gpoints);
502 }
503
504 /* For each point i, find out if it is a "topographic minimum", that is,
505 better than k neighbour points. */
506 for(unsigned int si=0; si<sampleNr; si++) {
507 points[si].topomin=0;
508
509 /* Compute the distances, scaled using xtol[] */
510 double sdist[sampleNr];
511 for(unsigned int sj=0; sj<sampleNr; sj++) {
512 sdist[sj]=0.0;
513 if(si==sj) {sdist[sj]=1.0E+99; continue;}
517 if(isfinite(d)) sdist[sj]+=d*d;
518 }
519 /* Distance should be computed as square root, but it does not affect
520 the order of distances; since sqrt() is slow, it is not done. */
521 // sdist[sj]=sqrt(sdist[sj]);
522 }
523
524 /* Find the closest neighbours for sample i */
525 unsigned int ni=0;
526 for(ni=0; ni<neighNr; ni++) {
527 unsigned int mini=0;
528 double minv=sdist[mini];
529 for(unsigned int sj=0; sj<sampleNr; sj++) {
530 if(sdist[sj]<minv) {minv=sdist[sj]; mini=sj;}
531 }
532 sdist[mini]=1.0E+99; // not used again
533 /* If point i is worse than any of the closest neighbours, then point i is certainly not
534 a topographic minimum; then stop this loop and go to the next point i+1 */
535 if(!(points[si].fvalue<points[mini].fvalue)) break;
536 }
537 if(ni<neighNr) continue; // continue with the next sample i
538 /* otherwise mark this as topographic minimum (TM) */
539 points[si].topomin=1;
540 /* and copy it to G point list (if there is space left) */
541 if(topoNr<sampleNr) {
543 gpoints[topoNr++].fvalue=points[si].fvalue;
544 }
545 /* Print TM point */
546 if(verbose>6) {
547 printf("TM point %d\t", 1+si);
549 printf("=> %e\n", points[si].fvalue);
550 fflush(stdout);
551 }
552 }
553 } // inner loop
554
555
556 if(verbose>3) {printf(" process gpoints.\n"); fflush(stdout);}
557 /* For each point i, find out if it is a "topographic minimum", that is,
558 better than k/2 neighbour points. */
559 for(unsigned int si=0; si<sampleNr; si++) {
560 gpoints[si].topomin=0;
561
562 /* Compute the distances, scaled using xtol[] */
563 double sdist[sampleNr];
564 for(unsigned int sj=0; sj<sampleNr; sj++) {
565 sdist[sj]=0.0;
566 if(si==sj) {sdist[sj]=1.0E+99; continue;}
570 if(isfinite(d)) sdist[sj]+=d*d;
571 }
572 }
573
574 /* Find the closest neighbours for sample i */
575 unsigned int ni=0;
576 for(ni=0; ni<neighNr; ni++) {
577 unsigned int mini=0;
578 double minv=sdist[mini];
579 for(unsigned int sj=0; sj<sampleNr; sj++) {
580 if(sdist[sj]<minv) {minv=sdist[sj]; mini=sj;}
581 }
582 sdist[mini]=1.0E+99; // not used again
583 /* If point i is worse than any of the closest neighbours, then point i is certainly not
584 a topographic minimum; then stop this loop and go to the next point i+1 */
585 if(!(gpoints[si].fvalue<gpoints[mini].fvalue)) break;
586 }
587 if(ni<neighNr) continue; // continue with the next sample i
588 /* otherwise mark this as topographic minimum (TM) */
589 gpoints[si].topomin=1;
590 /* Print TM point */
591 if(verbose>6) {
592 printf("GTM point %d\t", 1+si);
594 printf("=> %e\n", gpoints[si].fvalue);
595 fflush(stdout);
596 }
597 }
598 /* Copy TMs into the start of gpoint list and set topoNr accordingly */
599 topoNr=0;
600 for(unsigned int si=0; si<sampleNr; si++) if(gpoints[si].topomin>0) {
601 if(si>topoNr) {
603 gpoints[topoNr].fvalue=gpoints[si].fvalue;
604 }
605 topoNr++;
606 }
607 if(verbose>3) {printf(" %d topographical minima\n", topoNr); fflush(stdout);}
608 if(topoNr==0) {
609 if(verbose>0) {fprintf(stderr, "Warning: no topographical minima found\n"); fflush(stderr);}
610 continue;
611 }
612
613
614 } // end of iTGO
615
616
617 /* If requested, do local optimization for the TMs */
618 if(topoNr==0) {
619 if(verbose>0) {fprintf(stderr, "Error: no topographical minima found\n"); fflush(stderr);}
620 free(points); free(gpoints);
623 }
624 for(unsigned int si=0; si<sampleNr; si++) if(gpoints[si].topomin>0) {
625 if(verbose>2) printf("LO for TM point %d\n", 1+si);
626 if(verbose>3) {
627 printf("TM point %d before LO:\n", 1+si);
628 printf("%d\t", 1+si);
630 printf("=> %e\n", gpoints[si].fvalue); fflush(stdout);
631 }
633 nloptDuplicate(nlo, &lnlo);
635 lnlo.funval=points[si].fvalue;
641 gpoints[si].fvalue=lnlo.funval;
642 if(verbose>3) {
643 printf("TM point %d after LO:\n", 1+si);
644 printf("%d\t", 1+si);
646 printf("=> %e\n", gpoints[si].fvalue); fflush(stdout);
647 }
648 } else if(verbose>7) {printf(" LO failed\n"); fflush(stdout);}
650 nloptFree(&lnlo);
651 }
652
653 /* Find the best minimum and return it */
654 {
655 unsigned int mini=0;
656 double minv=gpoints[mini].fvalue;
657 for(unsigned int si=1; si<sampleNr; si++)
658 if(!(gpoints[si].fvalue>minv)) {minv=gpoints[si].fvalue; mini=si;}
660 nlo->funval=gpoints[mini].fvalue;
661 }
662
663 free(points); free(gpoints);
666}
◆ nloptLimitFixedNr()
|
extern |
Determine the number of fixed parameters, based on constraints.
- See also
- nloptInit, nloptFree, nloptAllocate
- Returns
- Returns the number of fixed parameters.
- Parameters
-
d Pointer to NLOPT.
Definition at line 333 of file nlopt.c.
336 {
339 unsigned int n=0;
343 if(!isnan(r) && r<macheps) n++;
344 }
345 return(n);
346}
Referenced by nloptIATGO(), nloptITGO1(), nloptITGO2(), nloptMPSO(), nloptSimplexARRS(), and nloptSimplexMS().
◆ nloptMeanP()
|
extern |
Calculate mean point from NLOPT plist.
Missing values are omitted from the mean.
Notice that if plist contains fixed parameters, their mean and SD may not be exactly correct because of limited accuracy of calculation with floating points; you may therefore need to set the mean to the limit and SD to zero afterwards.
- See also
- nloptInit, nloptFree, nloptAllocate, nloptAddP, nloptSortP
- Returns
- Returns TPCERROR status, TPCERROR_OK (0) when successful.
- Parameters
-
nlo Pointer to NLOPT structure. nr Number of points to include; enter 0 to use all. meanp Pointer to array of size totalNr where mean point is saved. sdp Pointer to array of size totalNr where point SD is saved; enter NULL, if not needed.
Definition at line 219 of file nlopt.c.
228 {
233
235 for(unsigned int j=0; j<dim; j++) {
236 unsigned int n=0;
237 meanp[j]=0.0;
238 for(unsigned int i=0; i<nr; i++) {
240 n++;
241 }
243 meanp[j]/=(double)n;
244 }
245
246 /* Calculate SDs if required */
248 for(unsigned int j=0; j<dim; j++) {
249 sdp[j]=0.0;
250 unsigned int n=0;
251 double sqrsum=0.0, sumsqr=0.0;
252 for(unsigned int i=0; i<nr; i++) {
254 sqrsum+=nlo->plist[i*(dim+1)+j];
256 n++;
257 }
258 }
259 if(n<2) continue; // SD=0 if n==1
260 sqrsum*=sqrsum;
261 double ff=sumsqr - sqrsum/(double)n;
262 if(!(ff>0.0)) sdp[j]=0.0;
263 else sdp[j]=sqrt( ff / (double)(n-1) );
264 }
266}
Referenced by nloptIATGO(), and nloptSimplexARRS().
◆ nloptMPSO()
|
extern |
Multi particle swarm optimization (MPSO).
- 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).
- See also
- nloptSimplex, nloptIATGO
- Todo
- Better stopping criteria.
- Parameters
-
nlo Pointer to NLOPT data. Counter funCalls is initially set to zero, and then increased here. Parameter maxFunCalls is used as one of the stopping criteria. Parameters xtol[] is used as one of the stopping criteria. maxIter Maximum number of iterations; set to zero to use the default; 10 is the minimum. nSwarms Number of swarms; set to zero to use the default; 2 is the minimum. nParticles Number of particles per swarm; set to zero to use the default; 5 is the minimum. wInertia Inertia; weight (0-1) for keeping particles own velocity and direction; enter NaN to use the default. wParticle Particle independence; weight how much particle is drawn towards the best point in its own memory; enter NaN to use the default. wSwarm Gravitation to the swarm; weight how much particle is drawn towards the best point in its own swarm memory; enter NaN to use the default. wGlobal Gravitation to the global optimum; weight how much particle is drawn towards the best point the memory of all the swarms. Set to zero for keeping swarms independent on each other. pDeath Probability of particle death and rebirth at random position; must be less than 0.1. pImmigration Probability of two particles switching swarms; must be less than 0.1. doLocal Use local optimization (Nelder-Mead downhill simplex) intermittently; using it (1) does not much increase the speed than not using it (0), but may increase the success rate. status Pointer to status data; enter NULL if not needed.
Definition at line 129 of file mpso.c.
163 {
164 FILE *fp=stdout;
166 if(verbose>0) {
167 fprintf(fp, "\n%s(NLOPT, %d, status)\n", __func__, doLocal);
168 fflush(fp);
169 }
170
171 if(nlo==NULL) {
174 }
178 }
180 if(verbose>0) {fprintf(stderr, "Error: too low limit for function calls.\n"); fflush(stderr);}
183 }
184
185 /* Check if any of the parameters are fixed */
188 if(verbose>1 && fixedParNr>0) fprintf(fp, "fixedParNr := %d\n", fixedParNr);
190 if(verbose>2) fprintf(fp, "fittedParNr := %d\n", fittedParNr);
191 if(fittedParNr<1) {
194 }
195
196 /* Check the tolerations */
200 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
203 }
204 }
205
206
207
208 if(verbose>2) fprintf(fp, "\nInitializing MPSO\n");
209
210 if(maxIter<10) maxIter=1000; // max nr of loops (150)
211 if(nSwarms<2) nSwarms=2+fittedParNr; // nr of swarms (3)
212 if(nParticles<5) nParticles=5+2*fittedParNr; // nr of particles (25)
213 if(!(wInertia>0.0)) wInertia=0.333; // inertia (0.333)
214 if(!(wParticle>0.0)) wParticle=1.00; // particle / cognitive (1.40)
215 if(!(wSwarm>0.0)) wSwarm=1.60; // swarm / social (1.40)
216 if(!isfinite(wGlobal)) wGlobal=0.40; // multi-swarm / global (0.40)
217 if(!isfinite(pDeath) || pDeath>=0.1) // probability of particle death-rebirth (0.005)
218 pDeath=0.005; else if(pDeath<0.0) pDeath=0.0;
219 if(!isfinite(pImmigration) || pImmigration>=0.1) // probability of particle immigration (0.025)
220 pImmigration=0.025; else if(pImmigration<0.0) pImmigration=0.0;
221 /* Initiate the number of function calls */
223 nlo->usePList=0;
225 nlo->plist=NULL;
226 }
227 nlo->funCalls=0;
228
229 unsigned int pFrustration=10; // Nr of non-successful movements when particle becomes frustrated
230 unsigned int nFrustration=15; // Nr of non-successful movements when swarm becomes frustrated
231
232 if(verbose>3) {
233 fprintf(fp, "maxIter := %u\n", maxIter);
235 fprintf(fp, "nSwarms := %u\n", nSwarms);
236 fprintf(fp, "nParticles := %u\n", nParticles);
237 fprintf(fp, "dim := %u\n", dim);
238 fprintf(fp, "wInertia := %g\n", wInertia);
239 fprintf(fp, "wParticle := %g\n", wParticle);
240 fprintf(fp, "swarm/social (wSwarm) := %g\n", wSwarm);
241 fprintf(fp, "multi-swarm/global (wGlobal) := %g\n", wGlobal);
242 fprintf(fp, "probability of death := %g\n", pDeath);
243 fprintf(fp, "probability of immigration := %g\n", pImmigration);
244 fprintf(fp, "pFrustration := %u\n", pFrustration);
245 fprintf(fp, "nFrustration := %u\n", nFrustration);
246 fflush(fp);
247 }
248
249 /* Allocate memory */
250 MPSO_MULTISWARM multi;
251 multi.dim=dim; multi.sn=nSwarms; multi.pn=nParticles;
252 multi.bestPosition=calloc(dim, sizeof(double));
253 multi.swarm=calloc(nSwarms, sizeof(MPSO_SWARM));
254 for(unsigned int si=0; si<nSwarms; si++) {
255 multi.swarm[si].particle=calloc(nParticles, sizeof(MPSO_PARTICLE));
256 multi.swarm[si].bestPosition=calloc(dim, sizeof(double));
257 for(unsigned int pi=0; pi<nParticles; pi++) {
258 multi.swarm[si].particle[pi].position=calloc(dim, sizeof(double));
259 multi.swarm[si].particle[pi].velocity=calloc(dim, sizeof(double));
260 multi.swarm[si].particle[pi].bestPosition=calloc(dim, sizeof(double));
261 }
262 }
263 /* Initialize Mersenne Twister MT19937 */
265
266
267
268 /* Check if initial guess is valid */
269 int initGuessAvailable=1;
270 double initGuessCost=nan("");
271 for(unsigned int i=0; i<dim; i++) {
275 }
276 if(initGuessAvailable) {
278 if(!isfinite(initGuessCost)) initGuessAvailable=0;
279 else if(verbose>3) {
280 fprintf(fp, "valid initial guess with cost=%g\n", initGuessCost); fflush(fp);
281 }
282 }
283
284 int carpetFill=0;
285
286 if(carpetFill) {
287 /* Fill the swarms with random parameters and their costs */
288 multi.bestCost=nan("");
289 for(unsigned int si=0; si<nSwarms; si++) {
290 for(unsigned int pi=0; pi<nParticles; pi++) {
291 nloptRandomPoint(multi.swarm[si].particle[pi].position, nlo->xlower, nlo->xupper, dim, &multi.mt);
292 nloptRandomPoint(multi.swarm[si].particle[pi].velocity, nlo->xlower, nlo->xupper, dim, &multi.mt);
293 for(unsigned int i=0; i<dim; i++)
294 multi.swarm[si].particle[pi].velocity[i]-=multi.swarm[si].particle[pi].position[i];
295 multi.swarm[si].particle[pi].cost=
297 nlo->funCalls++;
298 /* Currently, this is the best position of the particle */
299 multi.swarm[si].particle[pi].bestCost=multi.swarm[si].particle[pi].cost;
300 doubleCopy(multi.swarm[si].particle[pi].bestPosition,
301 multi.swarm[si].particle[pi].position, dim);
302 multi.swarm[si].particle[pi].since=0;
303 /* Set the best cost in this swarm */
304 if(pi==0) {
305 if(si==0 && initGuessAvailable) { // if initial guess available, use it here
306 multi.swarm[si].bestCost=initGuessCost;
308 } else {
309 multi.swarm[si].bestCost=multi.swarm[si].particle[pi].cost;
310 doubleCopy(multi.swarm[si].bestPosition, multi.swarm[si].particle[pi].position, dim);
311 }
312 } else {
313 if(!isfinite(multi.swarm[si].bestCost) ||
314 multi.swarm[si].bestCost>multi.swarm[si].particle[pi].cost)
315 {
316 multi.swarm[si].bestCost=multi.swarm[si].particle[pi].cost;
317 doubleCopy(multi.swarm[si].bestPosition, multi.swarm[si].particle[pi].position, dim);
318 }
319 }
320 } // next particle in this swarm
321 /* Set the best cost of all swarms */
322 if(si==0) {
323 multi.bestCost=multi.swarm[si].bestCost;
324 doubleCopy(multi.bestPosition, multi.swarm[si].bestPosition, dim);
325 } else {
326 if(!isfinite(multi.bestCost) || multi.bestCost>multi.swarm[si].bestCost) {
327 multi.bestCost=multi.swarm[si].bestCost;
328 doubleCopy(multi.bestPosition, multi.swarm[si].bestPosition, dim);
329 }
330 }
331 multi.swarm[si].since=0;
332 } // next swarm
333 multi.since=0;
334 } else {
335 /* For each swarm, make one random particle inside the limits */
336 unsigned int si=0;
337 if(initGuessAvailable) // if initial guess is available, use it as particle for first swarm
339 for(; si<nSwarms; si++)
340 nloptRandomPoint(multi.swarm[si].particle[0].position, nlo->xlower, nlo->xupper, dim, &multi.mt);
341 /* Then add other particles around it with Gaussian distribution */
342 double sd[dim];
344 for(unsigned int si=0; si<nSwarms; si++) {
345 for(unsigned int pi=1; pi<nParticles; pi++) {
346 nloptGaussianPoint(multi.swarm[si].particle[pi].position,
348 &multi.mt);
349 }
350 }
351 /* Add random velocity and direction for each particle */
353 double mvel[dim]; for(unsigned int i=0; i<dim; i++) mvel[i]=0.0;
354 for(unsigned int si=0; si<nSwarms; si++) {
355 for(unsigned int pi=0; pi<nParticles; pi++) {
356 nloptGaussianPoint(multi.swarm[si].particle[pi].velocity,
357 mvel, sd, NULL, NULL, dim, &multi.mt);
358 }
359 }
360 /* Calculate the cost of each particle, and find the best positions */
361 multi.bestCost=nan("");
362 for(unsigned int si=0; si<nSwarms; si++) {
363 multi.swarm[si].bestCost=nan("");
364 for(unsigned int pi=0; pi<nParticles; pi++) {
365 multi.swarm[si].particle[pi].cost=
367 nlo->funCalls++;
368 /* Currently, this is the best position of the particle */
369 multi.swarm[si].particle[pi].bestCost=multi.swarm[si].particle[pi].cost;
370 doubleCopy(multi.swarm[si].particle[pi].bestPosition,
371 multi.swarm[si].particle[pi].position, dim);
372 multi.swarm[si].particle[pi].since=0;
373 /* How about best in the swarm? */
374 if(!isfinite(multi.swarm[si].bestCost) ||
375 multi.swarm[si].bestCost>multi.swarm[si].particle[pi].cost)
376 {
377 multi.swarm[si].bestCost=multi.swarm[si].particle[pi].cost;
378 doubleCopy(multi.swarm[si].bestPosition,
379 multi.swarm[si].particle[pi].position, dim);
380 }
381 }
382 /* Does this swarm has the best particle of all? */
383 if(!isfinite(multi.bestCost) || multi.bestCost>multi.swarm[si].bestCost) {
384 multi.bestCost=multi.swarm[si].bestCost;
385 doubleCopy(multi.bestPosition, multi.swarm[si].bestPosition, dim);
386 }
387 multi.swarm[si].since=0;
388 }
389 multi.since=0;
390 }
391
392
393 /* Inside the loop, save the previous best position to check if fit is improving */
394 double lastBestPosition[dim];
395 unsigned int movementCounter=0; // counter for stopped movement
396 unsigned int convergenceCounter=0; // counter for converged swarms
397
398 unsigned int stopCostLimit=20, stopCost=0;
399 double lastBestCost=multi.bestCost;
400
401 unsigned int iterNr=0; // loop index
402 while(iterNr<maxIter && nlo->funCalls<nlo->maxFunCalls) {
403 iterNr++;
404 if(verbose>4) {
405 fprintf(fp, "-----------------------------\nloop %d\n", iterNr);
406 fprintf(fp, "bestCost=%g since %u iterations\n", multi.bestCost, multi.since);
408 fflush(fp);
409 }
410
411
412#if(0)
413 /* Check if each swarm has concentrated inside tolerance */
414 unsigned int swarmsConcentrated=0;
415 for(unsigned int si=0; si<nSwarms; si++) {
416 double /*pmean[dim],*/ psd[dim], vmean[dim], vsd[dim];
417 nloptMPSOpMean(&multi, si, pmean, psd);
418 nloptMPSOabsvMean(&multi, si, vmean, vsd);
419 if(verbose>5) {
420 fprintf(fp, "swarm %d bestCost=%g since %u\n", 1+si,
421 multi.swarm[si].bestCost, multi.swarm[si].since);
422 if(verbose>6) {
423 fprintf(fp, " bestPosition: %g", multi.swarm[si].bestPosition[0]);
424 for(unsigned int i=1; i<dim; i++) fprintf(fp, ", %g", multi.swarm[si].bestPosition[i]);
425 fprintf(fp, "\n");
426 }
427 if(verbose>9 || iterNr==1) {
428 for(unsigned int i=0; i<dim; i++)
429 fprintf(fp, " parameter %d: meanPosition %g +- %g\n", 1+i, pmean[i], psd[i]);
430 for(unsigned int i=0; i<dim; i++)
431 fprintf(fp, " parameter %d: meanAbsVelocity %g +- %g\n", 1+i, vmean[i], vsd[i]);
432 }
433 fflush(fp);
434 }
435 /* Swarm particle mean and best particle ever inside tolerance? */
436 /* Parameter SD below tolerance? */
437 unsigned int i=0;
438 for(i=0; i<dim; i++) {
441 }
442 if(i==dim) {
443 if(verbose>5) fprintf(fp, " swarm %d particles inside tolerance.\n", 1+si);
444 swarmsConcentrated++;
445 }
446 }
447 if(swarmsConcentrated==nSwarms) {
448 if(verbose>2) fprintf(fp, " all swarm particles shrunk inside swarm tolerance.\n");
449 break;
450 }
451#endif
452
453
454 /* Go through all swarms, and their particles */
455 unsigned int swarmsConcentrated=0;
456 multi.since++;
457 for(unsigned int si=0; si<nSwarms; si++) {
458 multi.swarm[si].since++;
459
460 /* Calculate statistics of the swarm */
461 double pmean[dim], psd[dim], vmean[dim], vsd[dim];
462 nloptMPSOpMean(&multi, si, pmean, psd);
463 nloptMPSOabsvMean(&multi, si, vmean, vsd);
464 if(verbose>5) {
465 fprintf(fp, "swarm %d bestCost=%g since %u\n", 1+si,
466 multi.swarm[si].bestCost, multi.swarm[si].since);
467 if(verbose>6) {
468 fprintf(fp, " bestPosition: %g", multi.swarm[si].bestPosition[0]);
469 for(unsigned int i=1; i<dim; i++) fprintf(fp, ", %g", multi.swarm[si].bestPosition[i]);
470 fprintf(fp, "\n");
471 }
472 if(verbose>7 || iterNr==1) {
473 for(unsigned int i=0; i<dim; i++)
474 fprintf(fp, " parameter %d: meanPosition %g +- %g\n", 1+i, pmean[i], psd[i]);
475 for(unsigned int i=0; i<dim; i++)
476 fprintf(fp, " parameter %d: meanAbsVelocity %g +- %g\n", 1+i, vmean[i], vsd[i]);
477 }
478 fflush(fp);
479 }
480 /* Swarm particle mean and best-ever particle inside tolerance? */
481 /* Parameter SD below tolerance? */
482 {
483 unsigned int i=0;
484 for(i=0; i<dim; i++) {
487 }
488 if(i==dim) {
489 if(verbose>5) fprintf(fp, " swarm %d particles inside tolerance.\n", 1+si);
490 swarmsConcentrated++;
491 }
492 }
493
494
495
496 /* Swarm can get frustrated if no advancement for a long time */
497 int swarmFrustrated=0;
498 if(multi.swarm[si].since>=nFrustration) {
499 swarmFrustrated=1; multi.swarm[si].since=0;
500 if(verbose>5) {fprintf(fp, "swarm is frustrated\n"); fflush(fp);}
501 }
502
503 /* After every 10 loops, try Nelder-Mead, if requested */
504 if(doLocal && !swarmFrustrated && (iterNr%10)==0) {
505 if(verbose>5) {fprintf(fp, "local optimization:\n"); fflush(fp);}
506 /* Best-ever position as starting point */
508 /* Deltas based on swarm particle SD, but making sure that above zero, unless parameter is fixed */
509 for(unsigned int i=0; i<dim; i++) {
513 }
514 if(verbose>7) {
515 fprintf(fp, "initial guesses and deltas, with cost %g:\n", multi.swarm[si].bestCost);
517 fflush(fp);
518 }
521 if(cost<multi.swarm[si].bestCost) {
522 multi.swarm[si].bestCost=cost;
524 multi.swarm[si].since=0;
525 if(verbose>7) {
526 fprintf(fp, "results: %e", multi.swarm[si].bestPosition[0]);
527 for(unsigned int i=1; i<dim; i++) fprintf(fp, ", %e", multi.swarm[si].bestPosition[i]);
528 fprintf(fp, ", with cost: %g\n", cost); fflush(fp);
529 }
530 } else {
531 if(verbose>10) fprintf(fp, "local optimization did not provide lower cost (%e)\n", cost);
532 }
533 } else {
534 if(verbose>10) fprintf(fp, "local optimization failed.\n");
535 }
536 }
537
538 /* Go through the particles of this swarm */
539 for(unsigned int pi=0; pi<nParticles; pi++) {
540
541 if(verbose>14) {
542 fprintf(fp, " particle %d bestCost=%g since %u\n", 1+pi,
543 multi.swarm[si].particle[pi].bestCost, multi.swarm[si].particle[pi].since);
544 fflush(fp);
545 }
546
547 /* Cast lots for whether it is time to kill particle and create a new;
548 but if cost is NaN, then most definitely kill it. */
549 int killed=0;
551 || isnan(multi.swarm[si].particle[pi].cost))
552 {
553 killed=1; multi.swarm[si].particle[pi].since=0;
554 /* Replace current position and velocity with random position and velocity */
555 if(verbose>12) {fprintf(fp, "death of particle %u\n", 1+pi); fflush(fp);}
556
557 nloptRandomPoint(multi.swarm[si].particle[pi].position, nlo->xlower, nlo->xupper, dim, &multi.mt);
558 nloptRandomPoint(multi.swarm[si].particle[pi].velocity, nlo->xlower, nlo->xupper, dim, &multi.mt);
559 for(unsigned int i=0; i<dim; i++)
560 multi.swarm[si].particle[pi].velocity[i]-=multi.swarm[si].particle[pi].position[i];
561 /* cost */
562 multi.swarm[si].particle[pi].cost=
564 nlo->funCalls++;
565 /* Check if this is the best position ever for this particle */
566 if(multi.swarm[si].particle[pi].bestCost>multi.swarm[si].particle[pi].cost) {
567 if(verbose>10) printf("Miracle: resurrection led to a new best!\n");
568 multi.swarm[si].particle[pi].bestCost=multi.swarm[si].particle[pi].cost;
569 doubleCopy(multi.swarm[si].particle[pi].bestPosition,
570 multi.swarm[si].particle[pi].position, dim);
571 }
572 }
573
574 /* Cast lots for whether it is time to immigrate */
575 int immigrated=0;
577 immigrated=1; multi.swarm[si].particle[pi].since=0;
578 /* Swap particle with a random particle in different swarm */
579 if(verbose>12) {fprintf(fp, " immigration\n"); fflush(fp);}
580 unsigned int sj;
581 if(nSwarms<4) {
582 sj=si+1; if(sj>=nSwarms) sj=si-1;
583 } else {
585 }
586 if(verbose>12) fprintf(fp, " swapping particle %d between swarms %d and %d\n", 1+pi, 1+si, 1+sj);
587 for(unsigned int i=0; i<dim; i++) {
588 /* position */
589 double d=multi.swarm[si].particle[pi].position[i];
590 multi.swarm[si].particle[pi].position[i]=multi.swarm[sj].particle[pi].position[i];
591 multi.swarm[sj].particle[pi].position[i]=d;
592 /* velocity */
593 d=multi.swarm[si].particle[pi].velocity[i];
594 multi.swarm[si].particle[pi].velocity[i]=multi.swarm[sj].particle[pi].velocity[i];
595 multi.swarm[sj].particle[pi].velocity[i]=d;
596 }
597 /* cost */
598 double d=multi.swarm[si].particle[pi].cost;
599 multi.swarm[si].particle[pi].cost=multi.swarm[sj].particle[pi].cost;
600 multi.swarm[sj].particle[pi].cost=d;
601 /* also the particles memory of its best ever cost and position follows it to another swarm */
602 d=multi.swarm[si].particle[pi].bestCost;
603 multi.swarm[si].particle[pi].bestCost=multi.swarm[sj].particle[pi].bestCost;
604 multi.swarm[sj].particle[pi].bestCost=d;
605 for(unsigned int i=0; i<dim; i++) {
606 double d=multi.swarm[si].particle[pi].bestPosition[i];
607 multi.swarm[si].particle[pi].bestPosition[i]=multi.swarm[sj].particle[pi].bestPosition[i];
608 multi.swarm[sj].particle[pi].bestPosition[i]=d;
609 }
610 /* Is the immigrated particle the best in its new swarm? */
611 if(multi.swarm[si].bestCost>multi.swarm[si].particle[pi].bestCost) {
612 multi.swarm[si].bestCost=multi.swarm[si].particle[pi].bestCost;
613 doubleCopy(multi.swarm[si].bestPosition, multi.swarm[si].particle[pi].bestPosition, dim);
614 }
615 if(multi.swarm[sj].bestCost>multi.swarm[sj].particle[pi].bestCost) {
616 multi.swarm[sj].bestCost=multi.swarm[sj].particle[pi].bestCost;
617 doubleCopy(multi.swarm[sj].bestPosition, multi.swarm[sj].particle[pi].bestPosition, dim);
618 }
619 }
620
621 /* If particle does not advance, it gets frustrated */
622 int frustrated=0;
623 if(multi.swarm[si].particle[pi].since>=nFrustration && !killed && !immigrated) {
624 frustrated=1;
625 }
626 if(!frustrated && !killed && !immigrated && multi.swarm[si].particle[pi].since>2) {
627 /* If particle has the worst position in the swarm, it will get frustrated, too */
628 frustrated=1;
629 for(unsigned int pj=0; pj<nParticles; pj++) if(pi!=pj)
630 if(multi.swarm[si].particle[pj].cost>multi.swarm[si].particle[pi].cost) {
631 frustrated=0; break;}
632 }
633 if(frustrated) {
634 multi.swarm[si].particle[pi].since=0;
635 /* Particle gets frustrated and will turn towards the swarm mean */
636 if(verbose>12) {
637 fprintf(fp, " particle %u became frustrated\n", 1+pi);
638 fprintf(fp, " current position: %g", multi.swarm[si].particle[pi].position[0]);
639 for(unsigned int i=1; i<dim; i++) fprintf(fp, ", %g", multi.swarm[si].particle[pi].position[i]);
640 fprintf(fp, " -> %g\n", multi.swarm[si].particle[pi].cost);
641 fflush(fp);
642 }
643 }
644
645
646 /* Update velocity in every dimension */
647 if(!frustrated) {
649 multi.swarm[si].particle[pi].velocity[i] =
650 wInertia*multi.swarm[si].particle[pi].velocity[i];
651 multi.swarm[si].particle[pi].velocity[i] += wParticle*mertwiRandomDouble1(&multi.mt)*
652 (multi.swarm[si].particle[pi].bestPosition[i] - multi.swarm[si].particle[pi].position[i]);
653 multi.swarm[si].particle[pi].velocity[i] += wSwarm*mertwiRandomDouble1(&multi.mt)*
654 (multi.swarm[si].bestPosition[i] - multi.swarm[si].particle[pi].position[i]);
655 if(wGlobal>0.0)
656 multi.swarm[si].particle[pi].velocity[i] += wGlobal*mertwiRandomDouble1(&multi.mt)*
657 (multi.bestPosition[i] - multi.swarm[si].particle[pi].position[i]);
658 }
659 if(swarmFrustrated) {
660 /* Additional random panic movement, if swarm is frustrated */
663 if(multi.swarm[si].particle[pi].velocity[i]>0.0)
664 multi.swarm[si].particle[pi].velocity[i]+=v;
665 else multi.swarm[si].particle[pi].velocity[i]-=v;
666 }
667 }
668 } else {
669 /* Frustrated particle turns towards the swarm mean */
671 multi.swarm[si].particle[pi].velocity[i] =
672 0.11*wInertia*multi.swarm[si].particle[pi].velocity[i] +
673 0.89*(pmean[i] - multi.swarm[si].particle[pi].position[i]);
674 }
675 }
676 if(wGlobal<0.0 && nSwarms>1) { // Swarms reject each other; not recommended
677 for(unsigned int sj=0; sj<nSwarms; sj++) if(si!=sj) {
679 multi.swarm[si].particle[pi].velocity[i] += wGlobal*mertwiRandomDouble1(&multi.mt)*
680 (multi.swarm[sj].bestPosition[i] - multi.swarm[si].particle[pi].position[i]);
681 }
682 }
683 }
684 if(verbose>15 || (verbose>12 && frustrated)) {
685 fprintf(fp, " updated velocities: %g", multi.swarm[si].particle[pi].velocity[0]);
686 for(unsigned int i=1; i<dim; i++)
687 fprintf(fp, ", %g", multi.swarm[si].particle[pi].velocity[i]);
688 fprintf(fp, "\n");
689 }
690
691 /* Update position */
692 for(unsigned int i=0; i<dim; i++)
693 multi.swarm[si].particle[pi].position[i]+=multi.swarm[si].particle[pi].velocity[i];
694 /* check that new position is inside limits */
695 for(unsigned int i=0; i<dim; i++)
697 multi.swarm[si].particle[pi].position[i]=nlo->xlower[i];
698 /* reduce velocity and change direction */
699 multi.swarm[si].particle[pi].velocity[i]*=-0.1;
701 multi.swarm[si].particle[pi].position[i]=nlo->xupper[i];
702 /* reduce velocity and change direction */
703 multi.swarm[si].particle[pi].velocity[i]*=-0.1;
704 }
705 if(verbose>15 || (verbose>12 && frustrated)) {
706 fprintf(fp, " updated positions: %g", multi.swarm[si].particle[pi].position[0]);
707 for(unsigned int i=1; i<dim; i++) fprintf(fp, ", %g", multi.swarm[si].particle[pi].position[i]);
708 fprintf(fp, "\n");
709 }
710
711
712 /* Update cost */
713 multi.swarm[si].particle[pi].cost=
715 nlo->funCalls++;
716 if(verbose>15 || (verbose>12 && frustrated))
717 fprintf(fp, " updated_cost=%g\n", multi.swarm[si].particle[pi].cost);
718
719
720 /* Check if this is the best position ever for this particle */
721 if(multi.swarm[si].particle[pi].bestCost>multi.swarm[si].particle[pi].cost) {
722 multi.swarm[si].particle[pi].bestCost=multi.swarm[si].particle[pi].cost;
723 doubleCopy(multi.swarm[si].particle[pi].bestPosition,
724 multi.swarm[si].particle[pi].position, dim);
725 multi.swarm[si].particle[pi].since=0;
726 } else {
727 multi.swarm[si].particle[pi].since++;
728 }
729
730 /* Check if the best position of this particle is the best ever in this swarm?
731 Note: check this separately, because of possible immigrations. */
732 if(multi.swarm[si].bestCost>multi.swarm[si].particle[pi].bestCost) {
733 multi.swarm[si].bestCost=multi.swarm[si].particle[pi].bestCost;
734 doubleCopy(multi.swarm[si].bestPosition, multi.swarm[si].particle[pi].bestPosition, dim);
735 multi.swarm[si].since=0;
736 }
737
738 } // next particle
739
740 /* Is the best particle of this swarm the best particle in all swarms? */
741 if(multi.bestCost>multi.swarm[si].bestCost) {
742 multi.bestCost=multi.swarm[si].bestCost;
743 doubleCopy(multi.bestPosition, multi.swarm[si].bestPosition, dim);
744 multi.since=0;
745 }
746
747 } // next swarm
748
749 /* Check if each swarm has concentrated inside tolerance */
750 if(swarmsConcentrated==nSwarms) {
751 if(verbose>2) fprintf(fp, " all swarm particles shrunk inside swarm tolerance.\n");
752 break;
753 }
754
755
756 /* Stop if the best particle in each swarm is inside tolerances of parameters;
757 all swarms have then converged to the same position. */
758 if(nSwarms>1 && iterNr>10) {
759 int swarmsConverged=1;
760 for(unsigned int i=0; i<dim; i++) {
761 for(unsigned int si=0; si<nSwarms; si++) {
762 double d=fabs(multi.bestPosition[i]-multi.swarm[si].bestPosition[i]);
764 }
765 if(!swarmsConverged) break;
766 }
767 if(swarmsConverged) {
768 convergenceCounter++;
769 if(convergenceCounter>5) {
770 if(verbose>3) fprintf(fp, "all swarms converged into one minimum.\n");
771// break;
772 }
773 }
774 }
775
776
777 /* Stop if the best overall position stops moving */
778 if(iterNr>=20) {
779 int moved=0;
780 for(unsigned int i=0; i<dim; i++) {
781 double d=fabs(lastBestPosition[i]-multi.bestPosition[i]);
783 }
784 if(moved==0) {
785 movementCounter++;
786 if(movementCounter>10) {
787 if(verbose>3)
788 fprintf(fp, "best position did not improve in %d last loops.\n", movementCounter);
789#if(0)
790 if(verbose>80 && movementCounter>10) {
791 fprintf(fp, "Why does the movement not continue?\n");
792 for(unsigned int si=0; si<nSwarms; si++) {
793 fprintf(fp, "Swarm %d : bestCost=%g\n", 1+si, multi.swarm[si].bestCost);
794 /* Search the currently best particle to get its also its velocity */
795 unsigned int pBest=0;
796 double pBestCost=multi.swarm[si].particle[0].cost;
797 for(unsigned int pi=1; pi<nParticles; pi++)
798 if(multi.swarm[si].particle[pi].cost<pBestCost) {
799 pBestCost=multi.swarm[si].particle[pi].cost; pBest=pi;
800 }
801 fprintf(fp, " currentBestCost=%g\n", pBestCost);
802 fprintf(fp, "\tDim\toPos\tcPos\tcVel\n");
803 for(unsigned int i=0; i<dim; i++) {
804 fprintf(fp, "\t%d\t%g\t%g\t%g\n", 1+i, multi.swarm[si].bestPosition[i],
805 multi.swarm[si].particle[pBest].position[i],
806 multi.swarm[si].particle[pBest].velocity[i]);
807 }
808 }
809 }
810#endif
811 }
812 }
813 if(movementCounter>20+dim) {
814 fprintf(fp, "best position did not move markedly in %d last loops.\n", movementCounter);
815 break;
816 }
817 }
818 doubleCopy(lastBestPosition, multi.bestPosition, dim);
819
820
821 /* Check that cost is improving */
822 if(multi.bestCost<lastBestCost) {
823 lastBestCost=multi.bestCost;
824 stopCost=0;
825 } else {
826 stopCost++;
827 }
828 if(stopCost>=stopCostLimit) {
829 if(verbose>1) fprintf(fp, "Cost did not improve in %d last iterations.\n", stopCost);
830 break;
831 }
832
833
834 } // next loop
835
836 /* Check the reason for loop exit */
837 if(iterNr>=maxIter) {
838 if(verbose>1) fprintf(fp, "exceeded the maximum number for loops.\n");
839 }
841 if(verbose>1) fprintf(fp, "exceeded the maximum number for function calls.\n");
842 }
843
844 /* Copy optimized parameters over initial values */
846
847 /* Analyze the swarms */
848 if(verbose>2) {
849 fprintf(fp, "\n---- MPSO end analysis ----\n");
850 fprintf(fp, "loops: %d\n", iterNr);
852 for(unsigned int si=0; si<nSwarms; si++) {
853 fprintf(fp, "\nSwarm %d : bestCost=%g\n", 1+si, multi.swarm[si].bestCost);
854 for(unsigned int i=0; i<dim; i++) {
855 fprintf(fp, "parameter %d: bestPosition %g\n", 1+i, multi.swarm[si].bestPosition[i]);
856 }
857 double pmean[dim], psd[dim];
858 nloptMPSOpMean(&multi, si, pmean, psd);
859 for(unsigned int i=0; i<dim; i++) {
860 fprintf(fp, "parameter %d: meanPosition %g +- %g\n", 1+i, pmean[i], psd[i]);
861 }
862 nloptMPSOabsvMean(&multi, si, pmean, psd);
863 for(unsigned int i=0; i<dim; i++) {
864 fprintf(fp, "parameter %d: meanAbsVelocity %g +- %g\n", 1+i, pmean[i], psd[i]);
865 }
866
867 }
868 fprintf(fp, "\nOverall bestCost=%g\n", multi.bestCost);
869 fprintf(fp, "at position: %g", multi.bestPosition[0]);
870 for(unsigned int i=1; i<dim; i++) fprintf(fp, ", %g", multi.bestPosition[i]);
871 fprintf(fp, "\n"); fflush(fp);
872 }
873
874 /* Free allocated memory */
875 for(unsigned int si=0; si<nSwarms; si++) {
876 for(unsigned int pi=0; pi<nParticles; pi++) {
877 free(multi.swarm[si].particle[pi].position);
878 free(multi.swarm[si].particle[pi].velocity);
879 free(multi.swarm[si].particle[pi].bestPosition);
880 }
881 free(multi.swarm[si].particle);
882 free(multi.swarm[si].bestPosition);
883 }
884 free(multi.bestPosition);
885 free(multi.swarm);
886
887
890}
void mertwiInitWithSeed64(MERTWI *mt, uint64_t seed)
Initialize the state vector mt[] inside data structure for Mersenne Twister MT19937 pseudo-random num...
Definition mertwi.c:94
double mertwiRandomDouble1(MERTWI *mt)
Generate a 64-bit double precision floating point pseudo-random number in the range of [0,...
Definition mertwi.c:216
double mertwiRandomExponential(MERTWI *mt, double mean)
Generate pseudo-random number with exponential distribution and specified mean.
Definition mertwi.c:322
uint64_t mertwiSeed64(void)
Make uint64_t seed for pseudo-random number generators.
Definition mertwi.c:76
int64_t mertwiRandomInt63(MERTWI *mt)
Generate a random number on [0, 2^63-1]-interval using Mersenne Twister MT19937.
Definition mertwi.c:197
void nloptMPSOabsvMean(MPSO_MULTISWARM *ms, unsigned int si, double *mean, double *sd)
Definition mpso.c:97
void nloptMPSOpMean(MPSO_MULTISWARM *ms, unsigned int si, double *mean, double *sd)
Definition mpso.c:76
◆ nloptPowellBrent()
Powell-Brent (Praxis) non-linear unconstrained optimization.
Praxis is a general purpose routine for the minimization of a function in several variables. The algorithm used is a modification of conjugate gradient search method by Powell, with changes by Brent, who gave an algol-w program, which served as a basis for this function.
References:
- Powell, MJD (1964) An efficient method for finding the minimum of a function in several variables without calculating derivatives. Computer Journal 7:155-162.
- Brent, RP (1973) Algorithms for minimization without derivatives. Prentice Hall, Englewood Cliffs. Dover 2013 republication.
- Returns
- enum tpcerror (TPCERROR_OK when successful).
- Parameters
-
nlo Pointer to NLOPT struct. Initial guess must be given in x[]. Initial step sizes must be given in xdelta[]. Note that constraints are not applied here. Parameter tolerances are used as stopping criteria. ktm Parameter tolerances (xtol[]) are given in the previous struct. Praxis stops if the criterion 2 * ||x[k]-x[k-1]|| <= sqrt(macheps) * ||x[k]|| + xtol is fulfilled more than ktm times for every parameter x. You should first try with ktm=1, and usually not higher than 4. status Pointer to status data; enter NULL if not needed.
Definition at line 59 of file praxis.c.
74 {
76 if(verbose>0) printf("%s()\n", __func__);
77 if(nlo==NULL) {
80 }
84 }
85
86
87 /* Check deltas */
88 unsigned int fixedNr=0;
93 }
94
95
96 /* Set a safe machine epsilon and related values */
98 if(verbose>1) printf("macheps := %g\n", macheps);
99 double m2=sqrt(macheps);
100 double m4=sqrt(m2);
101 if(verbose>1) {
102 printf("m2 := %e\n", m2);
103 printf("m4 := %e\n", m4);
104 }
105 double small=macheps*macheps;
106 double vsmall=small*small;
107 double large=1.0/small;
108 double vlarge=1.0/vsmall;
109 if(verbose>1) {
110 printf("small := %e\n", small);
111 printf("vsmall := %e\n", vsmall);
112 printf("large := %e\n", large);
113 printf("vlarge := %e\n", vlarge);
114 }
115
116 /* Set common initial step size and scaling parameter */
117 double h=0.0, scbd=1.0, t=0.0, t2;
118 {
119 unsigned int i;
120 double min=nan(""), max=nan("");
122 h+=nlo->xdelta[i];
123 t+=nlo->xtol[i];
126 }
127 h/=(double)(nlo->totalNr-fixedNr);
128 t/=(double)(nlo->totalNr-fixedNr);
129 if(verbose>2) {printf("xdelta_min := %g\nxdelta_max := %g\n", min, max); fflush(stdout);}
130 scbd+=log10(max/min); // >1 if param scales are different
131 t2=small+fabs(t); t=t2;
132 if(h<100.0*t) h=100.0*t;
133 }
134 if(verbose>1) {
135 printf("step := %g\n", h);
136 printf("scbd := %g\n", scbd);
137 printf("t(2) := %g\n", t);
138 fflush(stdout);
139 }
140 if(isnan(scbd)) { // checking that user provided xdeltas
143 }
144
145
146 /* Is the problem ill-conditioned (1) or not?
147 This variable is automatically set, when Praxis finds the problem to be
148 ill-conditioned during iterations. */
149 int illc=0;
150
151 /* Initialize */
152 if(verbose>2) {printf("initializing\n"); fflush(stdout);}
154 double ldfac; if(illc) ldfac=0.1; else ldfac=0.01;
155 unsigned int nl=0, kt=0;
156 unsigned int nf=0; // nr of objective function calls
157 double dmin=small;
158 double ldt=h;
159 double d[dim], y[dim], z[dim];
160 // double v[dim][dim];
161 double **v=malloc(dim*sizeof(double *)); if(v==NULL) {
164 }
165 for(i=0; i<dim; i++) {
166 v[i]=malloc(dim*sizeof(double));
167 if(v[i]==NULL) {
170 }
171 }
172 /* Make a copy of parameter list */
173 if(verbose>2) {printf("copying parameters\n"); fflush(stdout);}
174 double p[dim];
176
177 for(i=0; i<dim; i++) {
178 d[i]=y[i]=z[i]=0.0;
179 for(j=0; j<dim; j++) v[i][j]=(i==j ? 1.0 : 0.0);
180 }
181 double qd0=0.0, q0[dim], q1[dim];
183
184
185 /* Calculate function value with the initial guess */
186 if(verbose>2) {printf("evaluating initial guess\n"); fflush(stdout);}
187 double fx, qf1;
189 nf++;
190 if(!isfinite(fx)) {
191 for(i=0; i<dim; i++) free(v[i]);
192 free(v);
195 }
196 if(verbose>1) {printf("initial_objfunc := %g\n", fx); fflush(stdout);}
197
198
199 /* The main loop */
200 double sf, s, dn, qd1=0.0;
201 int done=0, nLoop=0, ret=0;
202 while(!done) {
203 nLoop++;
204 if(verbose>2) {
205 printf("\n-------------------------------\nloop %d, with %d fcalls\n", nLoop, nf);
206 fflush(stdout);
207 }
208 sf=d[0]; s=d[0]=0.0;
209
210 /* Minimize along the first direction */
211 ret=praxisMin(0, 2, &d[0], &s, fx, 0, &fx, v, dim, h, t,
212 macheps, m2, m4, small, dmin, ldt, &nl,
213 &nf, qd0, qd1, q0, q1,
215 if(ret) {
216 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
217 for(i=0; i<dim; i++) free(v[i]);
218 free(v);
221 }
222
223 if((sf<=(0.9*d[0])) || ((0.9*sf)>=d[0])) for(i=1; i<dim; i++) d[i]=0.0;
224
225 /* Other directions */
226 unsigned int k, kl=0, k2;
227 double sl, df, f1, lds;
228 for(k=1; k<dim; k++) {
229 if(verbose>5) printf("direction %d\n", k);
230 for(i=0; i<dim; i++) y[i]=p[i];
231 sf=fx; if(kt>0) illc=1;
232 do { // Usually only once, but twice if noticed here that ill-conditioned
233 kl=k; df=0.0;
234 /* If ill-conditioned, take random step to get off resolution valley */
235 if(illc) {
236 for(i=0; i<dim; i++) {
238 s=z[i]; for(j=0; j<dim; j++) p[j]+=s*v[j][i];
239 }
240 /* Calculate function value */
242 if(!isfinite(fx)) {
245 }
246 nf++;
247 }
248 /* Minimize along non-conjugate directions */
249 for(k2=k; k2<dim; k2++) {
250 if(verbose>6) printf(" non-conjugate direction %d\n", k2);
251 sl=fx; s=0.0;
252 ret=praxisMin(k2, 2, &d[k2], &s, fx, 0, &fx, v, dim, h, t,
253 macheps, m2, m4, small, dmin, ldt, &nl,
254 &nf, qd0, qd1, q0, q1,
256 if(ret) {
257 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
258 for(i=0; i<dim; i++) free(v[i]);
259 free(v);
262 }
263 if(illc) s=d[k2]*(s+z[k2])*(s+z[k2]); else s=sl-fx;
264 if(df<s) {df=s; kl=k2;}
265 }
266 if(illc) break;
267 if(df>=fabs(100.0*macheps*fx)) break;
268 illc=1;
269 } while(1);
270
271 /* Minimize along conjugate directions */
272 for(k2=0; k2<k; k2++) { // WHY NOT UNTIL k2==k ??????
273// for(k2=0; k2<=k; k2++) {
274 if(verbose>6) printf(" conjugate direction %d\n", k2);
275 s=0.0;
276 ret=praxisMin(k2, 2, &d[k2], &s, fx, 0, &fx, v, dim, h, t,
277 macheps, m2, m4, small, dmin, ldt, &nl,
278 &nf, qd0, qd1, q0, q1,
280 if(ret) {
281 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
282 for(i=0; i<dim; i++) free(v[i]);
283 free(v);
286 }
287 }
288 f1=fx; fx=sf; lds=0.0;
289 for(i=0; i<dim; i++) {
290 sl=p[i]; p[i]=y[i]; y[i]=sl-y[i]; sl=y[i]; lds+=sl*sl;
291 }
292 lds=sqrt(lds);
293 if(lds>small) {
294 for(i=kl-1; i>=k; i--) {
295 for(j=0; j<dim; j++) v[j][i+1]=v[j][i];
296 d[i+1]=d[i];
297 }
298 d[k]=0.0; for(i=0; i<dim; i++) v[i][k]=y[i]/lds;
299 ret=praxisMin(k, 4, &d[k], &lds, f1, 1, &fx, v, dim, h, t,
300 macheps, m2, m4, small, dmin, ldt, &nl,
301 &nf, qd0, qd1, q0, q1,
303 if(ret) {
304 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
305 for(i=0; i<dim; i++) free(v[i]);
306 free(v);
309 }
310 if(lds<=0.0) {lds=-lds; for(i=0; i<dim; i++) v[i][k]=-v[i][k];}
311 }
312 ldt*=ldfac; if(ldt<lds) ldt=lds;
313 for(i=0, t2=0.0; i<dim; i++) t2+=p[i]*p[i];
314 t2=sqrt(t2); t2*=m2; t2+=t;
315
316 if(verbose>1) printf(" objfunc := %g\n", fx);
317
318 /* Stopping criterion */
319 if(ldt>(0.5*t2)) kt=0; else kt++;
320 if(verbose>3 && kt>0) printf(" kt := %d\n", kt);
321 if(kt>ktm) {done=1; break;}
322
323 } // next direction k
324 if(done) break;
325
326 /*
327 * We are probably in a curved valley, try quadratic extrapolation
328 */
329 /* but only if we have more than one parameter to fit */
330 if(dim<2) {
331 if(nf>1000) done=1;
332 continue;
333 }
334
335 /* Praxis quad:
336 looks for the minimum along a curve defined by q0, q1, and x (q2). */
337 {
338 double qa, qb, qc, qs, ql;
339 qs=fx; fx=qf1; qf1=qs;
340 for(i=0, qd1=0.0; i<dim; i++) {
341 qs=p[i]; ql=q1[i]; p[i]=ql; q1[i]=qs;
342 qd1+=(qs-ql)*(qs-ql);
343 }
344 qd1=sqrt(qd1); ql=qd1; qs=0.0;
345 if(qd0>0.0 && qd1>0.0 && nl>=3*dim*dim) {
346 ret=praxisMin(-1, 2, &qs, &ql, qf1, 1, &fx, v, dim, h, t,
347 macheps, m2, m4, small, dmin, ldt, &nl,
348 &nf, qd0, qd1, q0, q1,
350 if(ret) {
351 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
352 for(i=0; i<dim; i++) free(v[i]);
353 free(v);
356 }
357 qa = ql*(ql-qd1)/(qd0*(qd0+qd1));
358 qb = (ql+qd0)*(qd1-ql)/(qd0*qd1);
359 qc = ql*(ql+qd0)/(qd1*(qd0+qd1));
360 } else {
361 fx=qf1; qa=qb=0.0; qc=1.0;
362 }
363 if(!isfinite(qa)) qa=0.0;
364 if(!isfinite(qb)) qb=0.0;
365 if(!isfinite(qc)) qc=1.0;
366 qd0=qd1;
367 for(i=0; i<dim; i++) {
368 qs=q0[i]; q0[i]=p[i];
369 p[i]=qa*s + qb*p[i] + qc*q1[i];
370 }
371 }
372 dn=0.0;
373 for(i=0; i<dim; i++) {d[i]=1.0/sqrt(d[i]); if(dn<d[i]) dn=d[i];}
374 for(j=0; j<dim; j++) {
375 s=d[j]/dn;
376 if(isfinite(s)) for(i=0; i<dim; i++) v[i][j]*=s;
377 }
378 /* Scale axis to reduce condition number */
379 if(scbd>1.0) {
380 s=vlarge;
381 for(i=0; i<dim; i++) {
382 sl=0.0;
383 for(j=0; j<dim; j++) sl+=v[i][j]*v[i][j];
384 z[i]=sqrt(sl); if(!isfinite(z[i]) || z[i]<m4) z[i]=m4;
385 if(s>z[i]) s=z[i];
386 }
387 for(i=0; i<dim; i++) {
388 sl=s/z[i]; z[i]=1.0/sl;
389 if(z[i]>scbd) {sl=1.0/scbd; z[i]=scbd;}
390 for(j=0; j<dim; j++) v[i][j]*=sl;
391 }
392 }
393 /* Calculate new set of orthogonal directions before repeating
394 the main loop; first, transpose v[][] for the Praxis minfit */
395 for(i=1; i<dim; i++) for(j=0; j<i; j++) {
396 s=v[i][j]; v[i][j]=v[j][i]; v[j][i]=s;
397 }
398 /* Call Praxis minfit to find the SVD of v[][].
399 This gives the principal values and principal directions of approximating
400 quadratic form without squaring the condition number */
401 ret=praxisMinfit((int)dim, macheps, vsmall, v, d);
402 if(ret) {
403 if(verbose>0) printf("Error %d in praxisMinfit()\n", ret);
404 for(i=0; i<dim; i++) free(v[i]);
405 free(v);
408 }
409 if(scbd>1.0) {
410 for(i=0; i<dim; i++) {s=z[i]; for(j=0; j<dim; j++) v[i][j]*=s;}
411 for(i=0; i<dim; i++) {
412 s=0.0; for(j=0; j<dim; j++) s+=v[j][i]*v[j][i];
413 s=sqrt(s); d[i]*=s; s=1.0/s; for(j=0; j<dim; j++) v[j][i]*=s;
414 }
415 }
416 for(i=0; i<dim; i++) {
417 if((dn*d[i]) > large) d[i]=vsmall; // dn was calculated above
418 else if((dn*d[i]) < small) d[i]=vlarge;
419 else d[i]=pow(dn*d[i], -2.0);
420 }
421 /* the new eigenvalues and eigenvectors */
422 praxisSort(d, v, dim);
423 dmin=d[dim-1]; if(dmin<small) dmin=small;
424 illc=(m2*d[0]) > dmin;
425
426
427
428
429 if(nf>1000) done=1;
430 }
431 if(verbose>1) printf("func_evals := %d\n", nf);
432
433
434 /* Free allocated memory */
435 for(i=0; i<dim; i++) free(v[i]);
436 free(v);
437
438 /* If ok, then copy parameters to the provided struct */
440
441
442 // Just to prevent warning for now
443 if(sf>0.0) sf=0.0;
444 if(kt>0) kt=0;
445 if(nl>0) nl=0;
446 if(ldfac>0.0) ldfac=0.0;
447 if(ktm>0) ktm=0;
448
449
451}
◆ nloptPrintP()
|
extern |
Print the contents of NLOPT plist.
- See also
- nloptInit, nloptFree, nloptAllocate, nloptAddP, nloptSortP
- Parameters
-
nlo Pointer to NLOPT structure. nr Number of points to print; enter 0 or a large number to print all. fp File pointer for the output.
Definition at line 274 of file nlopt.c.
281 {
282 if(nlo==NULL || fp==NULL) return;
284 fprintf(fp, "\nSampled points:\n");
287
289 for(unsigned int si=0; si<nr; si++) {
290 fprintf(fp, "%d\t", 1+si);
293 }
294 fflush(fp);
295}
Referenced by nloptIATGO(), and nloptSimplexARRS().
◆ nloptRandomPoint()
|
extern |
Create random parameters between specified limits.
- Precondition
- Uses rand(), therefore set seed for a new series of pseudo-random numbers; to produce truly random numbers (not just pseudo-random), do srand(time(0)) before calling this function. If no seed is set, then value 1 is used as default seed by rand().
- See also
- drand, drandRange, drandGaussian, drandExponential, mertwiRandomBetween
- Returns
- Returns non-zero in case of an error.
- Parameters
-
p Pointer to parameter list to be filled. low Pointer to parameter lower limits. up Pointer to parameter upper limits. n List length. mt Pointer to initiated and seeded Mersenne Twister MT19937 data structure; enter NULL to use drand() instead.
- See also
- mertwiInit, mertwiInitWithSeed64
Definition at line 28 of file rndpoint.c.
41 {
42 if(p==NULL || low==NULL || up==NULL) return(1);
43 if(n==0) return(0);
44
45 if(mt==NULL) {
46 for(unsigned int i=0; i<n; i++) {
47 double dif=up[i]-low[i];
49 }
50 } else {
51 for(unsigned int i=0; i<n; i++) {
52 double dif=up[i]-low[i];
54 }
55 }
56
57 return(0);
58}
Referenced by nloptIATGO(), nloptITGO1(), nloptITGO2(), nloptMPSO(), and nloptSimplexARRS().
◆ nloptRemoveEmpties()
|
extern |
Remove any empty parameters from NLOPT structure.
- See also
- nloptInit, nloptFree, nloptAllocate
Parameters which have NaN in place of xfull[] or xdelta[] are removed from the list; if found, totalNr is reduced accordingly.
Definition at line 398 of file nlopt.c.
400 {
402 unsigned int i=0, j;
403 while(i<d->totalNr) {
411 }
412 d->totalNr--;
413 }
414 return;
415}
◆ nloptSimplex()
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).
- See also
- nlopt1D, nloptPowellBrent, nloptSimplexARRS
- Parameters
-
nlo Pointer 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
-
maxIter Maximum number of iterations; set to zero to use the default. status Pointer to status data; enter NULL if not needed.
Definition at line 32 of file simplex.c.
47 {
48 FILE *fp=stdout;
50 if(verbose>1) fprintf(fp, "%s(NLOPT, %d, status)\n", __func__, maxIter);
51 if(nlo==NULL) {
54 }
58 }
59
60 /* Check that initial values are inside limits */
64 {
65 if(verbose>2) {
66 fprintf(fp, "parameter %d failed: %e <= %e <= %e\n",
68 fflush(fp);
69 }
72 }
73
74 /* If only one parameter free for fitting, then use 1-dimensional method instead */
76 if(verbose>2) fprintf(fp, "going for one-dimensional method.\n");
78 }
79
80 /* Initiate the simplex */
81 if(verbose>10) fprintf(fp, "Simplex initialization\n");
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];
103 simplexP[j][i] += nlo->xdelta[i];
105 /* parameter would exceed limits, try the other direction */
106 simplexP[j][i] -= 2.0*nlo->xdelta[i];
107 /* check again */
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 }
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 }
129 nlo->funCalls++;
130 if(verbose>12) fprintf(fp, " R[%d]=%g\n", j, simplexR[j]);
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++) {
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]);
215 /* and compute the response */
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 }
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]);
236 nlo->funCalls++;
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]);
259 nlo->funCalls++;
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]);
282 nlo->funCalls++;
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 */
339#endif
340
341 /* Check the objective function value */
342 if(!isfinite(simplexR[best])) {
345 }
346
347 /* Copy optimized parameters over initial values */
348 if(simplexR[best]<=initialR) {
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");
356 }
357
360}
unsigned int nloptForceLimits(unsigned int n, double *pLower, double *pUpper, double *p)
Enforce the model parameters within given limits.
Definition constraints.c:76
Referenced by nloptIATGO(), nloptITGO1(), nloptITGO2(), nloptMPSO(), nloptSimplexARRS(), and nloptSimplexMS().
◆ nloptSimplexARRS()
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).
- See also
- drandSeed(), nloptSimplex, nlopt1D, nloptPowellBrent
- Parameters
-
nlo Pointer 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
-
maxIter Maximum number of iterations; set to zero to use the default. status Pointer to status data; enter NULL if not needed.
Definition at line 373 of file simplex.c.
387 {
388 FILE *fp=stdout;
390 //verbose=10;
391
392 if(verbose>1) fprintf(fp, "%s(NLOPT, %d, status)\n", __func__, maxIter);
393 if(nlo==NULL) {
397
398 /* Check if any of the parameters are fixed */
400 if(verbose>2 && fixedParNr>0) fprintf(fp, "fixedParNr := %d\n", fixedParNr);
402 if(verbose>2) fprintf(fp, "fittedParNr := %d\n", fittedParNr);
403 if(fittedParNr<1) {
406 }
407
408 /* Check the tolerations */
412 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
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 */
425 double tol[dim];
427
428 /*
429 * Set initial xdelta based on parameter limits and original tolerances
430 */
431 for(unsigned int i=0; i<dim; i++) {
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);}
448 }
449
450 /* If initial guess was given by user, create random points with Gaussian distribution around it */
451 int initGuessAvailable=1;
456 }
459 if(verbose>2) {
460 fprintf(fp, "Initial guess: ");
463 }
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++) {
475 nlo->funCalls++;
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++) {
490 nlo->funCalls++;
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 */
504
505 /* Downhill simplex from the best point so far */
507 if(verbose>4) {
508 fprintf(fp, "LO: initial point with deltas:\n");
510 }
512 if(verbose>1) {fprintf(fp, " LO failed\n"); fflush(fp);}
514 } else {
515 if(verbose>4) {
516 fprintf(fp, "Point after LO:");
519 }
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 */
535 if(verbose>5) {
536 fprintf(fp, "best point so far:");
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) {
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 }
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++) {
564 nlo->xdelta[i]=sd[i];
565 /* Make sure that individual is not too close to zero */
567 } else {
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");
577 }
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:");
586 }
587 }
588
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 }
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 */
604
605 /* Get the best point so far */
608
609 /* Put back the original tolerances */
611
614}
◆ nloptSimplexMS()
Multi-Start-point Nelder-Mead (downhill simplex) optimization.
- Returns
- enum tpcerror (TPCERROR_OK when successful).
- See also
- nloptSimplex, nlopt1D, nloptPowellBrent
- Parameters
-
nlo Pointer 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
-
spNr Number of start points; set to zero to use the default. status Pointer to status data; enter NULL if not needed.
Definition at line 624 of file simplex.c.
638 {
639 FILE *fp=stdout;
641 //verbose=10;
642
643 if(verbose>1) fprintf(fp, "%s(NLOPT, %d, status)\n", __func__, spNr);
644 if(nlo==NULL) {
648
649 /* Check if any of the parameters are fixed */
651 if(verbose>2 && fixedParNr>0) fprintf(fp, "fixedParNr := %d\n", fixedParNr);
653 if(verbose>2) fprintf(fp, "fittedParNr := %d\n", fittedParNr);
654 if(fittedParNr<1) {
657 }
658
659 /* Check the tolerations */
663 if(verbose>0) {fprintf(stderr, "Error: invalid xtol[].\n"); fflush(stderr);}
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 */
675 double ss[spNr], sp[spNr*xNr], he[fittedParNr];
676 for(unsigned int si=0; si<spNr; si++) {
677 ss[si]=nan("");
680 unsigned int hi=0;
681 for(unsigned int xi=0; xi<xNr; 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++) {
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];
708 if(verbose>1) {fprintf(fp, " LO failed\n"); fflush(fp);}
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 }
731 nlo->funval=ss[bi];
732
733
736}
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
◆ nloptSortP()
|
extern |
Sort NLOPT plist into the order of increasing function values.
The first 'funCalls' items in plist are sorted.
- See also
- nloptInit, nloptFree, nloptAllocate, nloptAddP, nloptMeanP, nloptPrintP
- Returns
- Returns TPCERROR status, TPCERROR_OK (0) when successful.
- Parameters
-
nlo Pointer to NLOPT structure.
Definition at line 182 of file nlopt.c.
185 {
189
192 for(unsigned int i=0; i<nr-1; i++)
193 for(unsigned int j=i+1; j<nr; j++) {
196 {
197 for(unsigned int k=0; k<=dim; k++) {
200 }
201 }
202 }
204}
Referenced by nloptIATGO(), and nloptSimplexARRS().
◆ nloptWrite()
|
extern |
Write the contents of NLOPT structure to the specified file pointer.
- See also
- nloptInit, nloptFree, nloptAllocate
- Parameters
-
d Pointer to NLOPT fp Output file pointer
Definition at line 302 of file nlopt.c.
307 {
308 if(d==NULL || fp==NULL) return;
310 fprintf(fp, "param xfull xlower xupper xdelta xtol\n");
312 fprintf(fp, "%5d", 1+i);
318 }
323 fflush(fp);
324 return;
325}
Generated on for TPCCLIB by
1.15.0