-
Notifications
You must be signed in to change notification settings - Fork 138
/
Univariate.hs
803 lines (708 loc) · 28.4 KB
/
Univariate.hs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
-- | This module provides numerical routines for minimizing/maximizing one dimensional functions.
-- These methods are guaranteed to converge to the correct solution if your function has exactly one minimum.
--
-- The module is divided into simple and advanced interface.
-- The simple interface is reminiscent of most optimization packages.
-- The advanced interface uses the "History" monad to provide complex debugging information.
module HLearn.Optimization.Univariate
(
-- * Simple interface
-- | These simple functions provide an interface similar to Matlab/Ocatave's optimization toolkit.
-- (See for example the [fminunc documentation](http://www.mathworks.com/help/optim/ug/fminunc.html).)
-- They are recommended for quick prototypes and working from within GHCI.
-- Here's an example:
--
-- >>> let f x = (x-5)**2 :: Double
-- >>> fminunc f
-- 5.000000000000001
--
-- In many applications, we want to nest these optimization calls.
-- Debugging these nested optimiztions is notoriously difficult.
-- Consider the following example:
--
-- >>> let g x y = x*x + (y-5)**2
-- >>> fminunc (\x -> fminunc (\y -> g y x))
-- -2.821859578658409
--
-- It's clear that the value of @x@ that minimizes @g x y@ should be zero;
-- but we get a ridiculous answer instead.
-- The advanced interface gives us much more debugging information and control over the optimization process.
-- We can use it to solve problems that can't be solved with the simple interface.
fminunc
, fmaxunc
-- * Advanced interface
-- | The advanced interface is more complex, but has two advantages.
-- First, it lets us specify a detailed level of debugging information using the "History" monad.
-- Second, it lets us specify in more detail the algorithms used in the optimization.
--
-- Let's repeat the first example from the simple interface.
-- This time we will use Brent's method on "Rational" numbers and ask for higher precision than we could get with floating point numbers.
--
-- >>> let f x = (x-5)*(x-5) :: Rational
-- >>> :{
-- >>> traceHistory $ do
-- >>> lb <- lineBracket f 0 1
-- >>> fminuncM_bracketed_brent_ (return . f) lb (maxIterations 10 || stop_brent 1e-10)
-- >>> }:
-- Iterator_brent; x=5%1; f(x)=0%1
--
-- In this case, we get an exact answer.
-- The "traceHistory" command also gives us a record of how our optimization proceeded.
--
-- FIXME: Add an example of nested optimization
, UnivariateMinimizer
, UnivariateMinimizerLB
-- ** Brent's Method
-- | Brent's method uses parabolic interpolation to find the minimum.
-- This is a good black box method.
--
-- See <http://en.wikipedia.org/wiki/Iterator_brent%27s_method wikipedia>
-- and Section 10.2 of "Numerical Recipes in C" for more details.
, fminuncM_brent
, fminuncM_bracketed_brent
, stop_brent
, Iterator_brent
-- ** Brent's method with derivatives
-- | We can make Brent's method a bit faster sometimes by using the function's derivative.
-- This method can be slower than Brent's method for some functions.
--
-- See section 10.3 of "Numerical Recipes in C" for more details.
, fminuncM_bracketed_dbrent
, stop_dbrent
, Iterator_dbrent
-- ** Golden section search
-- | Golden section search is used to find the minimum of "poorly behaved" functions.
-- Other methods almost always converge faster.
--
-- See <http://en.wikipedia.org/wiki/Golden_section_search wikipedia>
-- and Section 10.1 of "Numerical Recipes in C" for more details.
, fminuncM_gss
, fminuncM_bracketed_gss
, stop_gss
, Iterator_gss
-- ** Line bracketing
-- | The univariate optimization methods above require an interval within which a minimum is guaranteed to exist.
--
-- See Section 10.1 of "Numerical Recipes in C" for more details.
--
-- FIXME:
-- Should we switch "LineBracket" to the "Interval" type provided by subhask?
-- This gives us a more consistent interface at the expensive of (sometimes) an extra function call.
, lineBracket
, lineBracketM
, LineBracket
-- * Test Functions
, slopes
, slopesWithDiscontinuity
-- ** Lipshitz Univariate
-- | These problems come from a paper by Famularo, Sergeyev, and Pugliese (hence fsp).
-- See <wwwinfo.deis.unical.it/~yaro/tests.pdf test problems for lipshitz univariate global optimization> for more details.
, fsp_problem1
)
where
import SubHask
import SubHask.Category.Trans.Derivative
import HLearn.History
-------------------------------------------------------------------------------
-- simple interface
-- | Finds the minimum of a function
{-# INLINABLE fminunc #-}
fminunc ::
( OrdField a
) => a -- ^ initial guess of the input with maximum value
-> (a -> a) -- ^ function to minimize
-> a -- ^ input with (approximately) the minimum value
fminunc x0 f = x1 $ evalHistory $ fminuncM_brent (maxIterations 50 || noProgress) x0 (return . f)
-- | Finds the maximum of a function
{-# INLINABLE fmaxunc #-}
fmaxunc ::
( OrdField a
) => a -- ^ initial guess of the input with maximum value
-> (a -> a) -- ^ function to maximize
-> a -- ^ input with (approximately) the maximum value
fmaxunc x0 f = fminunc x0 (-f)
-------------------------------------------------------------------------------
-- helpers
-- | The type of univariate minimization functions that require a starting point
type UnivariateMinimizer cxt opt a
= StopCondition (opt a) -- ^ controls when to stop
-> a -- ^ initial guess of the input with minimum value
-> (a -> History cxt a) -- ^ monadic function we're minimizing
-> History cxt (opt a) -- ^ result of the optimization
-- | The type of univariate minimization functions that require a "LineBracket" bounding the solution
type UnivariateMinimizerLB cxt opt a
= StopCondition (opt a) -- ^ controls when to stop
-> LineBracket a -- ^ search for the solution within this range
-> (a -> History cxt a) -- ^ monadic function we're minimizing
-> History cxt (opt a) -- ^ result of the optimization
-- | Converts an optimization function that expects a "LineBracket" input into one that only requires an initial starting point.
unbracket ::
( cxt (opt a)
, cxt (LineBracket a)
, cxt (Scalar a)
, cxt a
, OrdField a
) => UnivariateMinimizerLB cxt opt a
-> UnivariateMinimizer cxt opt a
unbracket fminuncM_bracketed stop x0 f = do
lb <- lineBracketM f (x0/2) (x0*2)
fminuncM_bracketed stop lb f
-------------------------------------------------------------------------------
-- line bracketing
-- | Variable names correspond to the algorithm in "Numerical Recipes in C"
data LineBracket a = LineBracket
{ _ax :: !a
, _bx :: !a
, _cx :: !a
, _fa :: !a
, _fb :: !a
, _fc :: !a
}
deriving (Typeable)
type instance Scalar (LineBracket a) = Scalar a
instance Show a => Show (LineBracket a) where
show lb = "LineBracket; a="++show (_ax lb)++"; b="++show (_bx lb)++"; c="++show (_cx lb)
-- | finds two points ax and cx between which a minimum is guaranteed to exist
-- this is a transliteration of the 'lineBracket' function from the \"Numerical
-- Recipes\" series
{-# INLINEABLE lineBracket #-}
lineBracket ::
( OrdField a
) => (a -> a) -- ^ the function we're bracketing
-> a -- ^ an initial guess for the lower bound
-> a -- ^ an initial guess for the upper bound
-> History cxt (LineBracket a)
lineBracket f = lineBracketM (return . f)
{-# INLINEABLE lineBracketM #-}
lineBracketM ::
( OrdField a
) => (a -> History cxt a) -- ^ the function we're bracketing
-> a -- ^ an initial guess for the lower bound
-> a -- ^ an initial guess for the upper bound
-> History cxt (LineBracket a)
lineBracketM !f !pt1 !pt2 = beginFunction "lineBracketM" $ do
let gold = 1.618034
fpt1 <- f pt1
fpt2 <- f pt2
let (ax,fax,bx,fbx) = if fpt1 > fpt2
then (pt1,fpt1,pt2,fpt2)
else (pt2,fpt2,pt1,fpt1)
let cx = bx+gold*(bx-ax)
fcx <- f cx
let lb0 = LineBracket
{ _ax = ax
, _bx = bx
, _cx = cx
, _fa = fax
, _fb = fbx
, _fc = fcx
}
iterate (step_LineBracket f) stop_LineBracket lb0
stop_LineBracket :: OrdField a => StopCondition (LineBracket a)
stop_LineBracket _ lb = if _fb lb /= _fb lb
then error "NaN in linebracket"
else return $ _fb lb <= _fc lb
step_LineBracket :: OrdField a
=> (a -> History cxt a)
-> LineBracket a
-> History cxt (LineBracket a)
step_LineBracket !f lb@(LineBracket ax bx cx fa fb fc) = do
let sign a b = if b>0 then abs a else -(abs a)
tiny = 1e-20
glimit = 100
gold = 1.618034
r = (bx-ax)*(fb-fc)
q = (bx-cx)*(fb-fa)
u = bx-((bx-cx)*q-(bx-ax)*r)/2*(sign (max (abs $ q-r) tiny) (q-r))
u' = cx+gold*(cx-bx)
ulim = bx+glimit*(cx-bx)
-- due to laziness, we will only evaluate the function if we absolutely have to
fu <- f u
fu' <- f u'
fulim <- f ulim
return $ if (bx-u)*(u-cx) > 0
then if fu < fc -- Got a minimum between a and c
then lb
{ _ax = bx
, _bx = u
, _fa = fb
, _fb = fu
}
else if fu > fb -- Got a minimum between a and u
then lb
{ _cx = u
, _fc = fu
}
else lb -- Parabolic fit was no use. Use default magnification
{ _ax = bx
, _bx = cx
, _cx = u'
, _fa = fb
, _fb = fc
, _fc = fu'
}
else if (cx-u)*(u-ulim) > 0 -- parabolic fit is between c and its allowed limit
then if fu < fc
then lb
{ _ax = cx
, _bx = u
, _cx = u'
, _fa = fc
, _fb = fu
, _fc = fu'
}
else lb
{ _ax = bx
, _bx = cx
, _cx = u
, _fa = fb
, _fb = fc
, _fc = fu
}
else if (u-ulim)*(ulim-cx) >= 0
then lb -- limit parabolic u to maximum allowed value
{ _ax = bx
, _bx = cx
, _cx = ulim
, _fa = fb
, _fb = fc
, _fc = fulim
}
else lb -- reject parabolic u, use default magnification
{ _ax = bx
, _bx = cx
, _cx = u'
, _fa = fb
, _fb = fc
, _fc = fu'
}
-------------------------------------------------------------------------------
-- golden section search
-- | Variable names correspond to the algorithm in "Numerical Recipes in C"
data Iterator_gss a = Iterator_gss
{ _gss_fxb :: !a
, _gss_fxc :: !a
, _gss_xa :: !a
, _gss_xb :: !a
, _gss_xc :: !a
, _gss_xd :: !a
}
deriving (Typeable)
type instance Scalar (Iterator_gss a) = Scalar a
instance (IsScalar a, POrd a) => Has_x1 Iterator_gss a where
x1 gss = if (_gss_fxb < _gss_fxc) gss then _gss_fxb gss else _gss_fxc gss
instance (IsScalar a, POrd a) => Has_fx1 Iterator_gss a where
fx1 = inf _gss_fxb _gss_fxc
instance (ClassicalLogic a, Show a, OrdField a) => Show (Iterator_gss a) where
show gss = "Iterator_gss; "++"x="++show x++"; f(x)="++show fx
where
x = gss_x gss
fx = gss_fx gss
gss_x gss = if _gss_fxb gss < _gss_fxc gss
then _gss_xb gss
else _gss_xc gss
gss_fx gss = min (_gss_fxb gss) (_gss_fxc gss)
---------------------------------------
-- | Golden section search initialized by a starting point
{-# INLINEABLE fminuncM_gss #-}
fminuncM_gss ::
( cxt (LineBracket a)
, cxt a
, OrdField a
) => UnivariateMinimizer cxt Iterator_gss a
fminuncM_gss = unbracket fminuncM_bracketed_gss
-- | Golden section search that only searches within a specified "LineBracket"
{-# INLINEABLE fminuncM_bracketed_gss #-}
fminuncM_bracketed_gss ::
( OrdField a
, cxt a
) => UnivariateMinimizerLB cxt Iterator_gss a
fminuncM_bracketed_gss stop (LineBracket ax bx cx fa fb fc) f = beginFunction "goldenSectionSearch_" $ do
let r = 0.61803399
c = 1-r
let xb = if abs (cx-bx) > abs (bx-ax) then bx else bx-c*(bx-ax)
xc = if abs (cx-bx) > abs (bx-ax) then bx+c*(cx-bx) else bx
fxb <- f xb
fxc <- f xc
let gss0 = Iterator_gss
{ _gss_fxb = fxb
, _gss_fxc = fxc
, _gss_xa = ax
, _gss_xb = xb
, _gss_xc = xc
, _gss_xd = cx
}
iterate (step_gss f) stop gss0
where
step_gss :: OrdField a => (a -> History cxt a) -> Iterator_gss a -> History cxt (Iterator_gss a)
step_gss f (Iterator_gss f1 f2 x0 x1 x2 x3) = if f2 < f1
then do
let x' = r*x2+c*x3
fx' <- f x'
return $ Iterator_gss f2 fx' x1 x2 x' x3
else do
let x' = r*x1+c*x0
fx' <- f x'
return $ Iterator_gss fx' f1 x0 x' x1 x2
where
r = 0.61803399
c = 1-r
-- | Does not stop until the independent variable is accurate to within the tolerance passed in.
{-# INLINEABLE stop_gss #-}
stop_gss :: OrdField a => a -> StopCondition (Iterator_gss a)
stop_gss tol _ (Iterator_gss _ _ !x0 !x1 !x2 !x3 )
= return $ abs (x3-x0) <= tol*(abs x1+abs x2)
-------------------------------------------------------------------------------
-- Brent's method
-- | Variable names correspond to the algorithm in "Numerical Recipes in C"
data Iterator_brent a = Iterator_brent
{ _brent_a :: !a
, _brent_b :: !a
, _brent_d :: !a
, _brent_e :: !a
, _brent_fv :: !a
, _brent_fw :: !a
, _brent_fx :: !a
, _brent_v :: !a
, _brent_w :: !a
, _brent_x :: !a
}
deriving Typeable
type instance Scalar (Iterator_brent a) = Scalar a
type instance Logic (Iterator_brent a) = Bool
instance Eq a => Eq_ (Iterator_brent a) where
b1==b2
= (_brent_a b1 == _brent_a b2)
&& (_brent_b b1 == _brent_b b2)
&& (_brent_d b1 == _brent_d b2)
&& (_brent_e b1 == _brent_e b2)
&& (_brent_fv b1 == _brent_fv b2)
&& (_brent_fw b1 == _brent_fw b2)
&& (_brent_fx b1 == _brent_fx b2)
&& (_brent_v b1 == _brent_v b2)
&& (_brent_w b1 == _brent_w b2)
&& (_brent_x b1 == _brent_x b2)
instance (Field a, IsScalar a, Show a) => Show (Iterator_brent a) where
show b = "Iterator_brent; x="++show (_brent_x b)++"; f(x)="++show (fx1 b)
instance IsScalar a => Has_x1 Iterator_brent a where
x1 = _brent_x
instance IsScalar a => Has_fx1 Iterator_brent a where
fx1 = _brent_fx
-- | Brent's method initialized by a starting point
{-# INLINEABLE fminuncM_brent #-}
fminuncM_brent ::
( cxt (LineBracket a)
, cxt a
, OrdField a
) => UnivariateMinimizer cxt Iterator_brent a
fminuncM_brent = unbracket fminuncM_bracketed_brent
-- | Brent's method that only searches within a specified "LineBracket"
{-# INLINEABLE fminuncM_bracketed_brent #-}
fminuncM_bracketed_brent ::
( OrdField a
, cxt a
) => UnivariateMinimizerLB cxt Iterator_brent a
fminuncM_bracketed_brent stop (LineBracket ax bx cx fa fb fc) f = beginFunction "brent" $ do
fbx <- f bx
iterate (step_brent_noTolerance f) stop $ Iterator_brent
{ _brent_a = min ax cx
, _brent_b = max ax cx
, _brent_d = zero
, _brent_e = zero
, _brent_v = bx
, _brent_w = bx
, _brent_x = bx
, _brent_fv = fbx
, _brent_fw = fbx
, _brent_fx = fbx
}
-- | This version of brent's method is slightly simpler than the one presented in "Numerical Recipes".
-- In particular, the "Numerical Recipes" version has a minor optimization where it does not reevaluate a function if the input changes below some tolerance.
-- Removing this check actually *improves* numeric stability.
-- We can use types the are more accurate than "Double"s (e.g. "Rational"s) where this check doesn't make sense.
step_brent_noTolerance ::
( OrdField a
) => (a -> History cxt a) -> Iterator_brent a -> History cxt (Iterator_brent a)
step_brent_noTolerance f brent@(Iterator_brent a b d e fv fw fx v w x) = do
let cgold = 0.3819660
xm = 0.5*(a+b)
r' = (x-w)*(fx-fv)
q' = (x-v)*(fx-fw)
p' = (x-v)*q'-(x-w)*r'
q'' = 2*(q'-r')
p'' = if q''>0 then -p' else p'
q''' = abs q''
etemp' = e
(d',e') = if abs p'' >= abs (0.5*q'''*etemp') || p'' <= q'''*(a-x) || p'' >= q'''*(b-x)
then let e'' = if x>=xm then a-x else b-x in (cgold*e'',e'')
else let d''=p''/q'''; u''=x+d'' in (d'',d)
u' = x+d'
fu' <- f u'
return $ if fu' <= fx
then brent
{ _brent_e = e'
, _brent_d = d'
, _brent_a = if u' >= x then x else a
, _brent_b = if u' >= x then b else x
, _brent_v = w
, _brent_w = x
, _brent_x = u'
, _brent_fv = fw
, _brent_fw = fx
, _brent_fx = fu'
}
else brent
{ _brent_e = e'
, _brent_d = d'
, _brent_a = if u' < x then u' else a
, _brent_b = if u' < x then b else u'
, _brent_v = if fu' <= fw || w==x then w else if fu' <= fv || v==x || v==w then u' else v
, _brent_fv = if fu' <= fw || w==x then fw else if fu' <= fv || v==x || v==w then fu' else fv
, _brent_w = if fu' <= fw || w==x then u' else w
, _brent_fw = if fu' <= fw || w==x then fu' else fw
}
-- | This version of Brent's method is a direct translation from "Numerical Recipes".
-- It is very slightly faster than "step_brent_noTolerance" but has worse accuracy.
step_brent_orig ::
( OrdField a
) => (a -> History cxt a) -> Iterator_brent a -> History cxt (Iterator_brent a)
step_brent_orig f brent@(Iterator_brent a b d e fv fw fx v w x) = do
let cgold = 0.3819660
xm = 0.5*(a+b)
sign a b = if b>0 then abs a else -(abs a)
zeps = 1e-10
tol = 1e-12
tol1' = tol*(abs x)+zeps
tol2' = 2*tol1'
(d',e') = if abs e > tol1'
then let
r' = (x-w)*(fx-fv)
q' = (x-v)*(fx-fw)
p' = (x-v)*q'-(x-w)*r'
q'' = 2*(q'-r')
p'' = if q''>0 then -p' else p'
q''' = abs q''
etemp' = e
in if abs p'' >= abs (0.5*q'''*etemp') || p'' <= q'''*(a-x) || p'' >= q'''*(b-x)
then let e'' = if x>=xm then a-x else b-x in (cgold*e'',e'')
else let d''=p''/q'''; u''=x+d'' in
if u''-a < tol2' || b-u'' < tol2'
then (sign tol1' (xm-x),d)
else (d'',d)
else let e'' = if x>=xm then a-x else b-x in (cgold*e'',e'')
u' = if abs d' >= tol1'
then x+d'
else x+sign tol1' d'
fu' <- f u'
return $ if fu' <= fx
then brent
{ _brent_e = e'
, _brent_d = d'
, _brent_a = if u' >= x then x else a
, _brent_b = if u' >= x then b else x
, _brent_v = w
, _brent_w = x
, _brent_x = u'
, _brent_fv = fw
, _brent_fw = fx
, _brent_fx = fu'
}
else brent
{ _brent_e = e'
, _brent_d = d'
, _brent_a = if u' < x then u' else a
, _brent_b = if u' < x then b else u'
, _brent_v = if fu' <= fw || w==x then w else if fu' <= fv || v==x || v==w then u' else v
, _brent_fv = if fu' <= fw || w==x then fw else if fu' <= fv || v==x || v==w then fu' else fv
, _brent_w = if fu' <= fw || w==x then u' else w
, _brent_fw = if fu' <= fw || w==x then fu' else fw
}
-- | Does not stop until the independent variable is accurate to within the tolerance passed in.
{-# INLINEABLE stop_brent #-}
stop_brent :: OrdField a => a -> StopCondition (Iterator_brent a)
stop_brent tol _ opt = return $ abs (x-xm) < tol2'-0.5*(b-a)
where
(Iterator_brent a b d e fv fw fx v w x) = opt
xm = 0.5*(a+b)
tol1' = tol*(abs x)+zeps
tol2' = 2*tol1'
zeps = 0 -- 1e-10
-------------------------------------------------------------------------------
-- Brent's method with derivatives
-- | Variable names correspond to the algorithm in "Numerical Recipes in C"
data Iterator_dbrent a = Iterator_dbrent
{ _dbrent_a :: !a
, _dbrent_b :: !a
, _dbrent_d :: !a
, _dbrent_e :: !a
, _dbrent_fv :: !a
, _dbrent_fw :: !a
, _dbrent_fx :: !a
, _dbrent_dv :: !a
, _dbrent_dw :: !a
, _dbrent_dx :: !a
, _dbrent_v :: !a
, _dbrent_w :: !a
, _dbrent_x :: !a
, _dbrent_break :: Bool
}
deriving Typeable
type instance Scalar (Iterator_dbrent a) = Scalar a
instance Show a => Show (Iterator_dbrent a) where
show b = "Iterator_brent; x="++show (_dbrent_x b)++"; f(x)="++show (_dbrent_fx b)
-- | Brent's method with derivatives that only searches within a specified "LineBracket"
--
-- FIXME:
-- Update this function so it uses the derivative category.
-- This will require some careful thought about:
-- 1) how derivatives interact with monads,
-- 2) restructuring the "UnivariateMinimization" types to handle arbitrary categories.
--
{-# INLINEABLE fminuncM_bracketed_dbrent #-}
fminuncM_bracketed_dbrent ::
( OrdField a
, cxt a
) => StopCondition (Iterator_dbrent a) -- ^ controls when to stop
-> LineBracket a -- ^ search for the solution within this range
-> (a -> History cxt a) -- ^ monadic function we're minimizing
-> (a -> History cxt a) -- ^ derivative of the function we're minimizing
-> History cxt (Iterator_dbrent a) -- ^ result of the optimization
fminuncM_bracketed_dbrent stop (LineBracket ax bx cx fa fb fc) f df = beginFunction "dbrent" $ do
fbx <- f bx
dfx <- df bx
iterate (step_dbrent f df) stop $ Iterator_dbrent
{ _dbrent_a = min ax cx
, _dbrent_b = max ax cx
, _dbrent_d = zero
, _dbrent_e = zero
, _dbrent_v = bx
, _dbrent_w = bx
, _dbrent_x = bx
, _dbrent_fv = fbx
, _dbrent_fw = fbx
, _dbrent_fx = fbx
, _dbrent_dv = dfx
, _dbrent_dw = dfx
, _dbrent_dx = dfx
, _dbrent_break = False
}
where
step_dbrent f df dbrent@(Iterator_dbrent a b d e fv fw fx dv dw dx v w x _) = do
let zeps = 1e-10
sign a b = if b>0 then abs a else -(abs a)
tol = 1e-6
xm = 0.5*(a+b)
tol1' = tol*(abs x)+zeps
tol2' = 2*tol1'
(d',e') = if abs e > tol1'
then let
d1 = if dw /= dx then (w-x)*dx/(dx-dw) else 2*(b-a)
d2 = if dv /= dx then (v-x)*dx/(dx-dv) else 2*(b-a)
u1 = x+d1
u2 = x+d2
ok1 = (a-u1)*(u1-b) > 0 && dx*d1 <= 0
ok2 = (a-u2)*(u2-b) > 0 && dx*d2 <= 0
in if ok1 || ok2
then let
d'' = if ok1 && ok2
then if abs d1 < abs d2 then d1 else d2
else if ok1 then d1 else d2
in if abs d'' <= abs (0.5 * e)
then let u' = x + d''
in if u'-a < tol2' || b-u' < tol2'
then (sign tol1' xm-x, d)
else (d'', d)
else
let e'' = if dx>=0 then a-x else b-x
in (0.5*e'',e'')
else
let e'' = if dx>=0 then a-x else b-x
in (0.5*e'',e'')
else
let e'' = if dx>=0 then a-x else b-x
in (0.5*e'',e'')
u' = if abs d' >= tol1' then x+d' else x+sign tol1' d'
fu' <- f u'
du' <- df u'
return $ if abs d' < tol1' && fu' > fx
then dbrent
{ _dbrent_x = x
, _dbrent_fx = fx
, _dbrent_break = True
}
else
if fu' <= fx
then dbrent
{ _dbrent_e = e'
, _dbrent_d = d'
, _dbrent_a = if u' >= x then x else a
, _dbrent_b = if u' >= x then b else x
, _dbrent_v = w
, _dbrent_w = x
, _dbrent_x = u'
, _dbrent_fv = fw
, _dbrent_fw = fx
, _dbrent_fx = fu'
, _dbrent_dv = dw
, _dbrent_dw = dx
, _dbrent_dx = du'
}
else dbrent
{ _dbrent_e = e'
, _dbrent_d = d'
, _dbrent_a = if u' < x then u' else a
, _dbrent_b = if u' < x then b else u'
, _dbrent_v = if fu' <= fw || w==x then w else if fu' <= fv || v==x || v==w then u' else v
, _dbrent_fv = if fu' <= fw || w==x then fw else if fu' <= fv || v==x || v==w then fu' else fv
, _dbrent_dv = if fu' <= fw || w==x then dw else if fu' <= fv || v==x || v==w then du' else dv
, _dbrent_w = if fu' <= fw || w==x then u' else w
, _dbrent_fw = if fu' <= fw || w==x then fu' else fw
, _dbrent_dw = if fu' <= fw || w==x then du' else dw
}
-- | Does not stop until the independent variable is accurate to within the tolerance passed in.
{-# INLINEABLE stop_dbrent #-}
stop_dbrent :: OrdField a => a -> StopCondition (Iterator_dbrent a)
stop_dbrent tol _ opt = return $ should_break || abs (x-xm) < tol2'-0.5*(b-a)
where
(Iterator_dbrent a b d e fv fw fx dv dw dx v w x should_break) = opt
xm = 0.5*(a+b)
tol1' = tol*(abs x)+zeps
tol2' = 2*tol1'
zeps = 1e-10
--------------------------------------------------------------------------------
-- test functions
-- | A simple piecewise linear function.
{-# INLINEABLE slopes #-}
slopes :: (ClassicalLogic r, OrdField r) => r -> r
slopes x
| x < 0 = -2*x
| x < 30 = -x
| x >= 30 = x-60
-- | This problem is interesting because it has no minimum when the domain is the real numbers.
-- Notice the discontinuity at x==30.
{-# INLINEABLE slopesWithDiscontinuity #-}
slopesWithDiscontinuity :: (ClassicalLogic r, OrdField r) => r -> r
slopesWithDiscontinuity x
| x < 0 = -2*x
| x < 30 = -x
| x >= 30 = x-60+10
{-# INLINEABLE fsp_problem1 #-}
fsp_problem1 :: Real r => r -> r
fsp_problem1 x = -13/6*x + sin (13/4*(2*x+5)) - 53/12
{-
tester :: History Double
tester = do
lb <- lineBracketM f 0 1
gss <- fminuncM_gss_ f lb (maxIterations 5)
return $ gss_x gss
where
f x = do
let f a = abs x + (a+5-x)*(a+5)
lb <- lineBracket f 0 1
b <- fminuncM_brent_
(return . f)
lb
(maxIterations 10)
return $ _brent_fx b
-------------------------------------------------------------------------------
-}