1 /* Minimizers.cpp
2 *
3 * Copyright (C) 2001-2020 David Weenink
4 *
5 * This code is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This code is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this work. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include "melder.h"
20 #include "NUM2.h"
21 #include "Graphics.h"
22 #include "Minimizers.h"
23
24 Thing_implement (Minimizer, Thing, 0);
25
classMinimizer_afterHook(Minimizer me,Thing)26 static void classMinimizer_afterHook (Minimizer me, Thing /* boss */) {
27 if (my success || ! my gmonitor)
28 return;
29 Graphics_beginMovieFrame (my gmonitor, nullptr);
30 Graphics_clearWs (my gmonitor);
31 Minimizer_drawHistory (me, my gmonitor, 0, my maximumNumberOfIterations, 0.0, 1.1 * my history [1], 1);
32 Graphics_textTop (my gmonitor, false, Melder_cat (U"Dimension of search space: ", my numberOfParameters));
33 Graphics_endMovieFrame (my gmonitor, 0.0);
34 Melder_monitor ((double) (my iteration) / my maximumNumberOfIterations, U"Iterations: ", my iteration,
35 U", Function calls: ", my numberOfFunctionCalls, U", Cost: ", my minimum);
36 }
37
Minimizer_init(Minimizer me,integer numberOfParameters,Daata object)38 void Minimizer_init (Minimizer me, integer numberOfParameters, Daata object) {
39 my numberOfParameters = numberOfParameters;
40 my p = zero_VEC (numberOfParameters);
41 my object = object;
42 my minimum = 1e308;
43 my afterHook = classMinimizer_afterHook;
44 autoVEC my_p;
45 Minimizer_reset (me, my_p.get()); // do weights initialization if my_p.size == 0
46 }
47
monitor_off(Minimizer me)48 static void monitor_off (Minimizer me) {
49 Melder_monitor (1.0);
50 my gmonitor = nullptr;
51 }
52
Minimizer_minimize(Minimizer me,integer maximumNumberOfIterations,double tolerance,int monitor)53 void Minimizer_minimize (Minimizer me, integer maximumNumberOfIterations, double tolerance, int monitor) {
54 try {
55 my tolerance = tolerance;
56 if (maximumNumberOfIterations <= 0)
57 return;
58 if (my iteration + maximumNumberOfIterations > my maximumNumberOfIterations) {
59 my maximumNumberOfIterations += maximumNumberOfIterations;
60 my history. resize (my maximumNumberOfIterations);
61 }
62 if (monitor)
63 my gmonitor = (Graphics) Melder_monitor (0.0, U"Starting...");
64 my v_minimize ();
65 if (monitor)
66 monitor_off (me);
67 if (my success)
68 Melder_casual (U"Minimizer_minimize:", U" minimum ", my minimum, U" reached \nafter ", my iteration,
69 U" iterations and ", my numberOfFunctionCalls, U" function calls.");
70 } catch (MelderError) {
71 Melder_clearError(); // memory error in history mechanism is not fatal
72 if (monitor)
73 monitor_off (me); // temporarily until better monitor facilities
74 }
75 }
76
Minimizer_minimizeManyTimes(Minimizer me,integer maxIterationsPerTime,integer numberOfTimes,double tolerance)77 void Minimizer_minimizeManyTimes (Minimizer me, integer maxIterationsPerTime, integer numberOfTimes, double tolerance) {
78 double fopt = my minimum;
79 int monitorSingle = numberOfTimes == 1;
80
81 autoVEC popt = raw_VEC (my numberOfParameters);
82 popt.all() <<= my p.all();
83
84 if (! monitorSingle)
85 Melder_progress (0.0, U"Minimize many times");
86
87 /* on first iteration start with current parameters 27/11/97 */
88 for (integer iter = 1; iter <= numberOfTimes; iter ++) {
89 Minimizer_minimize (me, maxIterationsPerTime, tolerance, monitorSingle);
90 Melder_casual (U"Current ", iter, U": minimum = ", my minimum);
91 if (my minimum < fopt) {
92 my p.all() <<= popt.all();
93 fopt = my minimum;
94 }
95 VEC p;
96 Minimizer_reset (me, p); // do initialization if p.size == 0
97 if (! monitorSingle) {
98 try {
99 Melder_progress ((double) iter / numberOfTimes, iter, U" from ", numberOfTimes);
100 } catch (MelderError) {
101 Melder_clearError (); // interrupted, no error
102 break;
103 }
104 }
105 }
106 if (! monitorSingle)
107 Melder_progress (1.0);
108 Minimizer_reset (me, popt.get());
109 }
110
Minimizer_reset(Minimizer me,constVEC const & guess)111 void Minimizer_reset (Minimizer me, constVEC const& guess) {
112 Melder_assert (guess.size == 0 || guess.size >= my numberOfParameters);
113 if (guess.size > 0)
114 my p.all() <<= guess;
115 else
116 for (integer i = 1; i <= my numberOfParameters; i ++)
117 my p [i] = NUMrandomUniform (-1.0, 1.0);
118
119 my history. resize (0);
120 my maximumNumberOfIterations = my numberOfFunctionCalls = my iteration = 0;
121 my success = false;
122 my minimum = 1.0e38;
123 my v_reset ();
124 }
125
Minimizer_drawHistory(Minimizer me,Graphics g,integer iFrom,integer iTo,double hmin,double hmax,bool garnish)126 void Minimizer_drawHistory (Minimizer me, Graphics g, integer iFrom, integer iTo, double hmin, double hmax, bool garnish) {
127 if (my history.size == 0)
128 return;
129 if (iTo <= iFrom) {
130 iFrom = 1;
131 iTo = my iteration;
132 }
133 integer itmin = iFrom, itmax = iTo;
134 Melder_clipLeft (1_integer, & itmin);
135 Melder_clipRight (& itmax, my iteration);
136 if (hmax <= hmin)
137 NUMextrema (my history.part (itmin, itmax), & hmin, & hmax);
138 if (hmax <= hmin) {
139 hmin -= 0.5 * fabs (hmin);
140 hmax += 0.5 * fabs (hmax);
141 }
142 Graphics_setInner (g);
143 Graphics_setWindow (g, iFrom, iTo, hmin, hmax);
144 Graphics_function (g, my history.asArgumentToFunctionThatExpectsOneBasedArray(), itmin, itmax, itmin, itmax);
145 Graphics_unsetInner (g);
146 if (garnish) {
147 Graphics_drawInnerBox (g);
148 Graphics_textBottom (g, true, U"Number of iterations");
149 Graphics_marksBottom (g, 2, true, true, false);
150 Graphics_marksLeft (g, 2, true, true, false);
151 }
152 }
153
Minimizer_getMinimum(Minimizer me)154 double Minimizer_getMinimum (Minimizer me) {
155 return my minimum;
156 }
157
158 /************** class SteepestDescentMinimizer **********************/
159
160 Thing_implement (SteepestDescentMinimizer, Minimizer, 0);
161
v_minimize()162 void structSteepestDescentMinimizer :: v_minimize () {
163 autoVEC dp = raw_VEC (numberOfParameters);
164 autoVEC dpp = raw_VEC (numberOfParameters);
165 double fret = func (object, p.get());
166 while (iteration < maximumNumberOfIterations) {
167 dfunc (object, p.get(), dp.get());
168 for (integer i = 1; i <= numberOfParameters; i ++) {
169 dpp [i] = - eta * dp [i] + momentum * dpp [i];
170 p [i] += dpp [i];
171 }
172 history [++ iteration] = minimum = func (object, p.get());
173 success = 2.0 * fabs (fret - minimum) < tolerance * (fabs (fret) + fabs (minimum));
174 if (our afterHook) {
175 try {
176 our afterHook (this, our afterBoss);
177 } catch (MelderError) {
178 Melder_casual (U"Interrupted after ", iteration, U" iterations.");
179 Melder_clearError ();
180 break;
181 }
182 }
183 if (success) break;
184 fret = minimum;
185 }
186 }
187
SteepestDescentMinimizer_create(integer numberOfParameters,Daata object,double (* func)(Daata object,VEC const & p),void (* dfunc)(Daata object,VEC const & p,VEC const &))188 autoSteepestDescentMinimizer SteepestDescentMinimizer_create (integer numberOfParameters, Daata object, double (*func) (Daata object, VEC const& p), void (*dfunc) (Daata object, VEC const& p, VEC const&)) {
189 try {
190 autoSteepestDescentMinimizer me = Thing_new (SteepestDescentMinimizer);
191 Minimizer_init (me.get(), numberOfParameters, object);
192 my func = func;
193 my dfunc = dfunc;
194 return me;
195 } catch (MelderError) {
196 Melder_throw (U"SteepestDescentMinimizer not created.");
197 }
198 }
199
200 /***************** class VDSmagtMinimizer ******************************/
201
202 Thing_implement (VDSmagtMinimizer, Minimizer, 0);
203
v_minimize()204 void structVDSmagtMinimizer :: v_minimize () {
205 int decrease_direction_found = 1;
206 int l_iteration = 1; // yes, we can iterate in steps, therefore local and global counter
207 longdouble rtemp, rtemp2;
208 /*
209 df is estimate of function reduction obtainable during line search
210 restart = 2 => line search in direction of steepest descent
211 restart = 1 => line search with Powell-restart.
212 flag = 1 => no decrease in function value during previous line search;
213 flag = 2 => line search did not decrease gradient
214 OK; must restart
215 */
216 if (restart_flag) {
217 minimum = func (object, p.get());
218 dfunc (object, p.get(), dp.get());
219 df = minimum;
220 restart = 2;
221 one_up = flag = 0;
222 gcg0 = gopt_sq = 0.0;
223 }
224 restart_flag = true;
225 while (++ this -> iteration <= maximumNumberOfIterations) {
226 if (flag & 1) {
227 if (one_up) {
228 decrease_direction_found = 0;
229 this -> iteration --;
230 break;
231 } else
232 one_up = 1;
233 } else
234 one_up = 0;
235
236 if (flag & 2)
237 restart = 2; /* flag & 1 ??? */
238 else if (fabs ((double) gcg0) > 0.2 * gopt_sq)
239 restart = 1;
240
241 if (restart == 0) {
242 rtemp = rtemp2 = 0.0;
243 for (integer i = 1; i <= numberOfParameters; i ++) {
244 rtemp += gc [i] * grst [i];
245 rtemp2 += gc [i] * srst [i];
246 }
247 gamma = rtemp / gamma_in;
248 if (fabs (beta * gropt - gamma * rtemp2) > 0.2 * gopt_sq)
249 restart = 1;
250 else
251 for (integer i = 1; i <= numberOfParameters; i ++)
252 s [i] = beta * s [i] + gamma * srst [i] - gc [i];
253 }
254 if (restart == 2) {
255 for (integer i = 1; i <= numberOfParameters; i ++)
256 s [i] = - dp [i];
257 restart = 1;
258 } else if (restart == 1) {
259 gamma_in = gropt - gr0;
260 for (integer i = 1; i <= numberOfParameters; i ++) {
261 srst [i] = s [i];
262 s [i] = beta * s [i] - gc [i];
263 grst [i] = gc [i] - g0 [i];
264 }
265 restart = 0;
266 }
267 /*
268 Begin line search
269 lineSearch_iteration = #iterations during current line search
270 */
271 flag = 0;
272 lineSearch_iteration = 0;
273 rtemp = 0.0;
274 for (integer i = 1; i <= numberOfParameters; i ++) {
275 rtemp += dp [i] * s [i];
276 g0 [i] = dp [i];
277 }
278 gr0 = gropt = rtemp;
279 if (l_iteration == 1)
280 alphamin = fabs (df / gropt);
281 if (gr0 > 0) {
282 flag = 1;
283 restart = 2;
284 continue;
285 }
286 f0 = minimum;
287 /*
288 alpha = length of step along line;
289 dalpha = change in alpha
290 alphamin = position of min along line
291 */
292 alplim = -1;
293 again = -1;
294 rtemp = fabs (df / gropt);
295 dalpha = std::min (alphamin, (double) rtemp);
296 alphamin = 0;
297 do {
298 do {
299 if (lineSearch_iteration) {
300 if (! (fch == 0))
301 gr2s += (temp + temp) / dalpha;
302
303 if (alplim < -0.5)
304 dalpha = 9.0 * alphamin;
305 else
306 dalpha = 0.5 * (alplim - alphamin);
307
308 grs = gropt + dalpha * gr2s;
309 if (gropt * grs < 0)
310 dalpha *= gropt / (gropt - grs);
311 }
312 alpha = alphamin + dalpha;
313 for (integer i = 1; i <= numberOfParameters; i ++)
314 pc [i] = p [i] + dalpha * s [i];
315 fc = func (object, pc.get());
316 dfunc (object, pc.get(), gc.get());
317 l_iteration ++;
318 lineSearch_iteration ++;
319 gsq = grc = 0.0;
320 for (integer i = 1; i <= numberOfParameters; i ++) {
321 gsq += gc [i] * gc [i];
322 grc += gc [i] * s [i];
323 }
324 fch = fc - minimum;
325 gr2s = (grc - gropt) / dalpha;
326 temp = (fch + fch) / dalpha - grc - gropt;
327 if ((fc < minimum) || ((fc == minimum) && (grc / gropt > -1))) {
328 gopt_sq = gsq;
329 history [this -> iteration] = minimum = fc;
330 std::swap (p, pc);
331 std::swap (dp, gc);
332 if (grc * gropt <= 0.0)
333 alplim = alphamin;
334 alphamin = alpha;
335 gropt = grc;
336 dalpha = - dalpha;
337 success = ( gsq < tolerance );
338 if (our afterHook) {
339 try {
340 our afterHook (this, our afterBoss);
341 } catch (MelderError) {
342 Melder_casual (U"Interrupted after ", this -> iteration, U" iterations.");
343 Melder_clearError ();
344 break;
345 }
346 }
347 if (success)
348 return;
349 if (fabs (gropt / gr0) < lineSearchGradient)
350 break;
351 } else {
352 alplim = alpha;
353 }
354 } while (lineSearch_iteration <= lineSearchMaxNumOfIterations);
355
356 fc = history [this -> iteration] = minimum;
357 rtemp = 0.0;
358 for (integer i = 1; i <= numberOfParameters; i ++) {
359 pc [i] = p [i];
360 gc [i] = dp [i];
361 rtemp += gc [i] * g0 [i];
362 }
363 gcg0 = rtemp;
364 if (fabs (gropt - gr0) > tolerance) {
365 beta = (gopt_sq - gcg0) / (gropt - gr0);
366 if (fabs (beta * gropt) < 0.2 * gopt_sq) break;
367 }
368 again ++;
369 if (again > 0) flag += 2;
370 } while (flag < 1);
371
372 if (f0 <= minimum)
373 flag += 1;
374 df = gr0 * alphamin;
375 }
376 if (this -> iteration > maximumNumberOfIterations)
377 this -> iteration = maximumNumberOfIterations;
378 if (decrease_direction_found)
379 restart_flag = false;
380 }
381
v_reset()382 void structVDSmagtMinimizer :: v_reset () {
383 restart_flag = true;
384 }
385
VDSmagtMinimizer_create(integer numberOfParameters,Daata object,double (* func)(Daata object,VEC const & x),void (* dfunc)(Daata object,VEC const & x,VEC const & dx))386 autoVDSmagtMinimizer VDSmagtMinimizer_create (integer numberOfParameters, Daata object, double (*func) (Daata object, VEC const& x), void (*dfunc) (Daata object, VEC const& x, VEC const& dx)) {
387 try {
388 autoVDSmagtMinimizer me = Thing_new (VDSmagtMinimizer);
389 Minimizer_init (me.get(), numberOfParameters, object);
390 my dp = zero_VEC (numberOfParameters);
391 my pc = zero_VEC (numberOfParameters);
392 my gc = zero_VEC (numberOfParameters);
393 my g0 = zero_VEC (numberOfParameters);
394 my s = zero_VEC (numberOfParameters);
395 my srst = zero_VEC (numberOfParameters);
396 my grst = zero_VEC (numberOfParameters);
397 my func = func;
398 my dfunc = dfunc;
399 my lineSearchGradient = 0.9;
400 my lineSearchMaxNumOfIterations = 5;
401 return me;
402 } catch (MelderError) {
403 Melder_throw (U"VDSmagtMinimizer not created.");
404 }
405 }
406
407 /************ class LineMinimizer *******************************/
408
409 Thing_implement (LineMinimizer, Minimizer, 0);
410
LineMinimizer_init(LineMinimizer me,integer numberOfParameters,Daata object,double (* func)(Daata,VEC const & p))411 void LineMinimizer_init (LineMinimizer me, integer numberOfParameters, Daata object, double (*func) (Daata, VEC const& p)) {
412 Minimizer_init (me, numberOfParameters, object);
413 my direction = zero_VEC (numberOfParameters);
414 my ptry = zero_VEC (numberOfParameters);
415 my func = func;
416 my maxLineStep = 100;
417 }
418
419 /* End of file Minimizers.cpp */
420