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

Nonlinear one-dimensional optimization. More...

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

Go to the source code of this file.

Functions

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

Detailed Description

Nonlinear one-dimensional optimization.

Author
Vesa Oikonen

Definition in file nlopt1d.c.

Function Documentation

◆ nlopt1D()

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

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
nloPointer 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
maxIterMaximum number of function evaluations; set to zero to use the default.
statusPointer to status data; enter NULL if not needed

Definition at line 25 of file nlopt1d.c.

38 {
39 int verbose=0; if(status!=NULL) verbose=status->verbose;
40 if(verbose>1) {printf("%s(NLOPT, %d, status)\n", __func__, maxIter); fflush(stdout);}
41 if(nlo==NULL) {
42 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
43 return TPCERROR_FAIL;
44 }
45 if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
46 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
47 return TPCERROR_NO_DATA;
48 }
49
50 /* Check the number of fixed parameters, based on constraints or delta */
51 if(nlo->totalNr-nloptFixedNr(nlo)!=1) {
52 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_PARNR);
53 return TPCERROR_INVALID_PARNR;
54 }
55 /* Which parameter is the free one? */
56 double macheps=doubleMachEps();
57 unsigned int fpi;
58 for(fpi=0; fpi<nlo->totalNr; fpi++) {
59 double r=nlo->xupper[fpi]-nlo->xlower[fpi];
60 if(isfinite(nlo->xdelta[fpi]) && fabs(nlo->xdelta[fpi])>macheps && isfinite(r) && r>macheps)
61 break;
62 }
63 if(fpi==nlo->totalNr) {
64 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_PARNR);
65 return TPCERROR_INVALID_PARNR;
66 }
67 if(verbose>3) printf(" index_of_free_parameter := %u\n", fpi);
68
69 /* Check that initial value is inside its limits */
70 double pmin=nlo->xlower[fpi];
71 double pmax=nlo->xupper[fpi];
72 double delta=nlo->xdelta[fpi];
73 double tol=nlo->xtol[fpi];
74 if(verbose>2)
75 printf(" initial_parameter := %g [%g, %g, %g]\n", nlo->xfull[fpi], pmin, pmax, delta);
76 if(!(nlo->xfull[fpi]>=pmin) || !(nlo->xfull[fpi]<=pmax)) {
77 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
78 return TPCERROR_INVALID_VALUE;
79 }
80 /* Check that limits are wider than 2*delta */
81 if(!( (pmax-pmin) > 2.0*delta)) {
82 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
83 return TPCERROR_INVALID_VALUE;
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 */
94 double p[nlo->totalNr];
95 for(unsigned int i=0; i<nlo->totalNr; i++) p[i]=nlo->xfull[i];
96
97 /* Calculate function value with initial parameter value */
98 unsigned int iterNr=1;
99 p2=minp=nlo->xfull[fpi]; if(p2<pmin+delta) p2=pmin+delta; else if(p2>pmax-delta) p2=pmax-delta;
100 p[fpi]=p2; minf=f2=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata);
101 if(!isfinite(f2)) {
102 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
103 return TPCERROR_INVALID_VALUE;
104 }
105 nlo->funCalls++;
106 if(nlo->usePList) {
107 int ret=nloptAddP(nlo, p, f2);
108 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
109 }
110
111 /* Set the two other bracketing points and compute their function values */
112 p1=p2-delta; p[fpi]=p1; f1=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
113 nlo->funCalls++;
114 if(nlo->usePList) {
115 int ret=nloptAddP(nlo, p, f1);
116 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
117 }
118 if(f1<minf) {minf=f1; minp=p1;}
119
120 p3=p2+delta; p[fpi]=p3; f3=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
121 nlo->funCalls++;
122 if(nlo->usePList) {
123 int ret=nloptAddP(nlo, p, f3);
124 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
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);
148 p[fpi]=p2; f2=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
149 nlo->funCalls++;
150 if(nlo->usePList) {
151 int ret=nloptAddP(nlo, p, f2);
152 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
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;
159 p[fpi]=p1; f1=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
160 nlo->funCalls++;
161 if(nlo->usePList) {
162 int ret=nloptAddP(nlo, p, f1);
163 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
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);
174 p[fpi]=p2; f2=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
175 nlo->funCalls++;
176 if(nlo->usePList) {
177 int ret=nloptAddP(nlo, p, f2);
178 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
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;
185 p[fpi]=p3; f3=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
186 nlo->funCalls++;
187 if(nlo->usePList) {
188 int ret=nloptAddP(nlo, p, f3);
189 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
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 */
235 p[fpi]=minp; minf=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
236 nlo->funCalls++;
237 if(nlo->usePList) {
238 int ret=nloptAddP(nlo, p, minf);
239 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
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 */
256 nlo->xfull[fpi]=p2; // Not minp!
257 /* Compute function value at the minimum, in case user uses the fitted data in nlo structure */
258 nlo->funval=(*nlo->_fun)(nlo->totalNr, nlo->xfull, nlo->fundata);
259
260 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
261 return TPCERROR_OK;
262}
doubleMachEps
double doubleMachEps()
Definition doubleutil.c:105
nloptFixedNr
unsigned int nloptFixedNr(NLOPT *d)
Definition nlopt.c:354
nloptAddP
int nloptAddP(NLOPT *nlo, double *p, double funval)
Definition nlopt.c:149
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
NLOPT::_fun
double(* _fun)(int, double *, void *)
Definition tpcnlopt.h:42
NLOPT::funval
double funval
Definition tpcnlopt.h:50
NLOPT::xupper
double * xupper
Definition tpcnlopt.h:33
NLOPT::xlower
double * xlower
Definition tpcnlopt.h:31
NLOPT::fundata
void * fundata
Definition tpcnlopt.h:44
NLOPT::xfull
double * xfull
Definition tpcnlopt.h:29
NLOPT::funCalls
unsigned int funCalls
Definition tpcnlopt.h:48
NLOPT::xdelta
double * xdelta
Definition tpcnlopt.h:36
NLOPT::usePList
unsigned int usePList
Definition tpcnlopt.h:52
NLOPT::xtol
double * xtol
Definition tpcnlopt.h:39
NLOPT::totalNr
unsigned int totalNr
Definition tpcnlopt.h:27
TPCSTATUS::verbose
int verbose
Verbose level, used by statusPrint() etc.
Definition tpcextensions.h:242
TPCERROR_INVALID_VALUE
@ TPCERROR_INVALID_VALUE
Invalid value.
Definition tpcextensions.h:220
TPCERROR_INVALID_PARNR
@ TPCERROR_INVALID_PARNR
Invalid number of parameters.
Definition tpcextensions.h:222
TPCERROR_FAIL
@ TPCERROR_FAIL
General error.
Definition tpcextensions.h:192
TPCERROR_OK
@ TPCERROR_OK
No error.
Definition tpcextensions.h:191
TPCERROR_NO_DATA
@ TPCERROR_NO_DATA
File contains no data.
Definition tpcextensions.h:203

Referenced by nloptSimplex().

Generated on for TPCCLIB by doxygen 1.15.0