levmar 0.1 → 0.2
raw patch · 15 files changed
+2326/−1148 lines, 15 filesdep +MemoTriedep +vector-spacedep ~bindings-levmar
Dependencies added: MemoTrie, vector-space
Dependency ranges changed: bindings-levmar
Files
- Demo.hs +1447/−990
- LevMar.hs +45/−59
- LevMar/AD.hs +142/−0
- LevMar/Fitting.hs +60/−45
- LevMar/Fitting/AD.hs +139/−0
- LevMar/Intermediate.hs +59/−16
- LevMar/Intermediate/AD.hs +101/−0
- LevMar/Intermediate/Fitting.hs +54/−5
- LevMar/Intermediate/Fitting/AD.hs +100/−0
- LevMar/Utils.hs +44/−0
- LevMar/Utils/AD.hs +42/−0
- NFunction.hs +2/−2
- SizedList.hs +60/−27
- TypeLevelNat.hs +4/−0
- levmar.cabal +27/−4
Demo.hs view
@@ -1,994 +1,1451 @@ -- This module is a Haskell translation of lmdemo.c from the C levmar library. -module Demo where--import LevMar ( levmar-- , Model- , Jacobian-- , Options(..), defaultOpts-- , LinearConstraints, noLinearConstraints-- , LevMarError-- , Info, CovarMatrix-- , S, Z- , SizedList(..)- )--import qualified LevMar.Fitting as Fitting--type Result n = Either LevMarError- ( SizedList n Double- , Info Double- , CovarMatrix n Double- )------------------------------------------------------------------------------------- Handy type synonyms for type-level naturals:--type N0 = Z-type N1 = S N0-type N2 = S N1-type N3 = S N2-type N4 = S N3-type N5 = S N4------------------------------------------------------------------------------------- Default options:--opts :: Options Double-opts = defaultOpts { optStopNormInfJacTe = 1e-15- , optStopNorm2Dp = 1e-15- , optStopNorm2E = 1e-20- }------------------------------------------------------------------------------------- Rosenbrock function,--- global minimum at (1, 1)--ros :: Model N2 Double-ros p0 p1 = replicate ros_n ((1.0 - p0)**2 + ros_d*m**2)- where- m = p1 - p0**2--ros_jac :: Jacobian N2 Double-ros_jac p0 p1 = replicate ros_n (p0d ::: p1d ::: Nil)- where- p0d = -2 + 2*p0 - 4*ros_d*m*p0- p1d = 2*ros_d*m- m = p1 - p0**2--ros_d :: Double-ros_d = 105.0--ros_n :: Int-ros_n = 2--ros_params :: SizedList N2 Double-ros_params = -1.2 ::: 1.0 ::: Nil--ros_samples :: [Double]-ros_samples = replicate ros_n 0.0--run_ros :: Result N2-run_ros = levmar ros- (Just ros_jac)- ros_params- ros_samples- 1000- opts- Nothing- Nothing- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Modified Rosenbrock problem,--- global minimum at (1, 1)--modros :: Model N2 Double-modros p0 p1 = [ 10*(p1 - p0**2)- , 1.0 - p0- , modros_lam- ]--modros_jac :: Jacobian N2 Double-modros_jac p0 _ = [ -20*p0 ::: 10.0 ::: Nil- , -1.0 ::: 0.0 ::: Nil- , 0.0 ::: 0.0 ::: Nil- ]--modros_lam :: Double-modros_lam = 1e02--modros_n :: Int-modros_n = 3--modros_params :: SizedList N2 Double-modros_params = -1.2 ::: 1.0 ::: Nil--modros_samples :: [Double]-modros_samples = replicate modros_n 0.0--run_modros :: Result N2-run_modros = levmar modros- (Just modros_jac)- modros_params- modros_samples- 1000- opts- Nothing- Nothing- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Powell's function,--- minimum at (0, 0)--powell :: Model N2 Double-powell p0 p1 = [ p0- , 10.0*p0 / m + 2*p1**2- ]- where- m = p0 + 0.1--powell_jac :: Jacobian N2 Double-powell_jac p0 p1 = [ 1.0 ::: 0.0 ::: Nil- , 1.0 / m**2 ::: 4.0*p1 ::: Nil- ]- where- m = p0 + 0.1--powell_n :: Int-powell_n = 2--powell_params :: SizedList N2 Double-powell_params = -1.2 ::: 1.0 ::: Nil--powell_samples :: [Double]-powell_samples = replicate powell_n 0.0--run_powell :: Result N2-run_powell = levmar powell- (Just powell_jac)- powell_params- powell_samples- 1000- opts- Nothing- Nothing- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Wood's function,--- minimum at (1, 1, 1, 1)--wood :: Model N4 Double-wood p0 p1 p2 p3 = [ 10.0*(p1 - p0**2)- , 1.0 - p0- , sqrt 90.0*(p3 - p2**2)- , 1.0 - p2- , sqrt 10.0*(p1 + p3 - 2.0)- , (p1 - p3) / sqrt 10.0- ]--wood_n :: Int-wood_n = 6--wood_params :: SizedList N4 Double-wood_params = -3.0 ::: -1.0 ::: -3.0 ::: -1.0 ::: Nil--wood_samples :: [Double]-wood_samples = replicate wood_n 0.0--run_wood :: Result N4-run_wood = levmar wood- Nothing- wood_params- wood_samples- 1000- opts- Nothing- Nothing- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Meyer's (reformulated) data fitting problem,--- minimum at (2.48, 6.18, 3.45)--meyer :: Fitting.Model N3 Double Double-meyer p0 p1 p2 x = p0*exp (10.0*p1 / (ui + p2) - 13.0)- where- ui = 0.45 + 0.05*x--meyer_jac :: Fitting.Jacobian N3 Double Double-meyer_jac p0 p1 p2 x = tmp- ::: 10.0*p0*tmp / (ui + p2)- ::: -10.0*p0*p1*tmp / ((ui + p2)*(ui + p2))- ::: Nil- where- tmp = exp (10.0*p1 / (ui + p2) - 13.0)- ui = 0.45 + 0.05*x--meyer_n :: Int-meyer_n = 16--meyer_params :: SizedList N3 Double-meyer_params = 8.85 ::: 4.0 ::: 2.5 ::: Nil--meyer_samples :: [(Double, Double)]-meyer_samples = zip [0..] [ 34.780- , 28.610- , 23.650- , 19.630- , 16.370- , 13.720- , 11.540- , 9.744- , 8.261- , 7.030- , 6.005- , 5.147- , 4.427- , 3.820- , 3.307- , 2.872- ]--run_meyer_jac :: Result N3-run_meyer_jac = Fitting.levmar meyer- (Just meyer_jac)- meyer_params- meyer_samples- 1000- opts- Nothing- Nothing- noLinearConstraints- Nothing--run_meyer :: Result N3-run_meyer = Fitting.levmar meyer- Nothing- meyer_params- meyer_samples- 1000- opts- Nothing- Nothing- noLinearConstraints- Nothing------------------------------------------------------------------------------------- helical valley function,--- minimum at (1.0, 0.0, 0.0)--helval :: Model N3 Double-helval p0 p1 p2 = [ 10.0*(p2 - 10.0*theta)- , 10.0*sqrt tmp - 1.0- , p2- ]- where- m = atan (p1 / p0) / (2.0*pi)-- tmp = p0**2 + p1**2-- theta | p0 < 0.0 = m + 0.5- | 0.0 < p0 = m- | p1 >= 0 = 0.25- | otherwise = -0.25--heval_jac :: Jacobian N3 Double-heval_jac p0 p1 _ = [ 50.0*p1 / (pi*tmp) ::: -50.0*p0 / (pi*tmp) ::: 10.0 ::: Nil- , 10.0*p0 / sqrt tmp ::: 10.0*p1 / sqrt tmp ::: 0.0 ::: Nil- , 0.0 ::: 0.0 ::: 1.0 ::: Nil- ]- where- tmp = p0**2 + p1**2--helval_n :: Int-helval_n = 3--helval_params :: SizedList N3 Double-helval_params = -1.0 ::: 0.0 ::: 0.0 ::: Nil--helval_samples :: [Double]-helval_samples = replicate helval_n 0.0--run_helval :: Result N3-run_helval = levmar helval- (Just heval_jac)- helval_params- helval_samples- 1000- opts- Nothing- Nothing- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Boggs - Tolle problem 3 (linearly constrained),--- minimum at (-0.76744, 0.25581, 0.62791, -0.11628, 0.25581)------ constr1: p0 + 3*p1 = 0--- constr2: p2 + p3 - 2*p4 = 0--- constr3: p1 - p4 = 0--bt3 :: Model N5 Double-bt3 p0 p1 p2 p3 p4 = replicate bt3_n (t1**2 + t2**2 + t3**2 + t4**2)- where- t1 = p0 - p1- t2 = p1 + p2 - 2.0- t3 = p3 - 1.0- t4 = p4 - 1.0--bt3_jac :: Jacobian N5 Double-bt3_jac p0 p1 p2 p3 p4 = replicate bt3_n ( 2.0*t1- ::: 2.0*(t2 - t1)- ::: 2.0*t2- ::: 2.0*t3- ::: 2.0*t4- ::: Nil- )- where- t1 = p0 - p1- t2 = p1 + p2 - 2.0- t3 = p3 - 1.0- t4 = p4 - 1.0--bt3_n :: Int-bt3_n = 5--bt3_params :: SizedList N5 Double-bt3_params = 2.0 ::: 2.0 ::: 2.0 :::2.0 ::: 2.0 ::: Nil--bt3_samples :: [Double]-bt3_samples = replicate bt3_n 0.0--bt3_linear_constraints :: LinearConstraints N3 N5 Double-bt3_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)- ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)- ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)- ::: Nil- , 0.0 ::: 0.0 ::: 0.0 ::: Nil- )--run_bt3 :: Result N5-run_bt3 = levmar bt3- (Just bt3_jac)- bt3_params- bt3_samples- 1000- opts- Nothing- Nothing- (Just bt3_linear_constraints)- Nothing------------------------------------------------------------------------------------- Hock - Schittkowski problem 28 (linearly constrained),--- minimum at (0.5, -0.5, 0.5)------ constr1: p0 + 2*p1 + 3*p2 = 1--hs28 :: Model N3 Double-hs28 p0 p1 p2 = replicate hs28_n (t1**2 + t2**2)- where- t1 = p0 + p1- t2 = p1 + p2--hs28_jac :: Jacobian N3 Double-hs28_jac p0 p1 p2 = replicate hs28_n ( 2.0*t1- ::: 2.0*(t1 + t2)- ::: 2.0*t2- ::: Nil- )- where- t1 = p0 + p1- t2 = p1 + p2--hs28_n :: Int-hs28_n = 3--hs28_params :: SizedList N3 Double-hs28_params = -4.0 ::: 1.0 ::: 1.0 ::: Nil--hs28_samples :: [Double]-hs28_samples = replicate hs28_n 0.0--hs28_linear_constraints :: LinearConstraints N1 N3 Double-hs28_linear_constraints = ( ((1.0 ::: 2.0 ::: 3.0 ::: Nil) ::: Nil)- , 1.0 ::: Nil- )--run_hs28 :: Result N3-run_hs28 = levmar hs28- (Just hs28_jac)- hs28_params- hs28_samples- 1000- opts- Nothing- Nothing- (Just hs28_linear_constraints)- Nothing------------------------------------------------------------------------------------- Hock - Schittkowski problem 48 (linearly constrained),--- minimum at (1.0, 1.0, 1.0, 1.0, 1.0)------ constr1: sum [p0, p1, p2, p3, p4] = 5--- constr2: p2 - 2*(p3 + p4) = -3--hs48 :: Model N5 Double-hs48 p0 p1 p2 p3 p4 = replicate hs48_n (t1**2 + t2**2 + t3**2)- where- t1 = p0 - 1.0- t2 = p1 - p2- t3 = p3 - p4--hs48_jac :: Jacobian N5 Double-hs48_jac p0 p1 p2 p3 p4 = replicate hs48_n ( 2.0*t1- ::: 2.0*t2- ::: 2.0*t2- ::: 2.0*t3- ::: 2.0*t3- ::: Nil- )- where- t1 = p0 - 1.0- t2 = p1 - p2- t3 = p3 - p4--hs48_n :: Int-hs48_n = 3--hs48_params :: SizedList N5 Double-hs48_params = 3.0 ::: 5.0 ::: -3.0 ::: 2.0 ::: -2.0 ::: Nil--hs48_samples :: [Double]-hs48_samples = replicate hs48_n 0.0--hs48_linear_constraints :: LinearConstraints N2 N5 Double-hs48_linear_constraints = ( (1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: Nil)- ::: (0.0 ::: 0.0 ::: 1.0 ::: -2.0 ::: -2.0 ::: Nil)- ::: Nil- , 5.0 ::: -3.0 ::: Nil- )--run_hs48 :: Result N5-run_hs48 = levmar hs48- (Just hs48_jac)- hs48_params- hs48_samples- 1000- opts- Nothing- Nothing- (Just hs48_linear_constraints)- Nothing------------------------------------------------------------------------------------- Hock - Schittkowski problem 51 (linearly constrained),--- minimum at (1.0, 1.0, 1.0, 1.0, 1.0)------ constr1: p0 + 3*p1 = 4--- constr2: p2 + p3 - 2*p4 = 0--- constr3: p1 - p4 = 0--hs51 :: Model N5 Double-hs51 p0 p1 p2 p3 p4 = replicate hs51_n (t1**2 + t2**2 + t3**2 + t4**2)- where- t1 = p0 - p1- t2 = p1 + p2 - 2.0- t3 = p3 - 1.0- t4 = p4 - 1.0--hs51_jac :: Jacobian N5 Double-hs51_jac p0 p1 p2 p3 p4 = replicate hs51_n ( 2.0*t1- ::: 2.0*(t2 - t1)- ::: 2.0*t2- ::: 2.0*t3- ::: 2.0*t4- ::: Nil- )- where- t1 = p0 - p1- t2 = p1 + p2 - 2.0- t3 = p3 - 1.0- t4 = p4 - 1.0--hs51_n :: Int-hs51_n = 5--hs51_params :: SizedList N5 Double-hs51_params = 2.5 ::: 0.5 ::: 2.0 ::: -1.0 ::: 0.5 ::: Nil--hs51_samples :: [Double]-hs51_samples = replicate hs51_n 0.0--hs51_linear_constraints :: LinearConstraints N3 N5 Double-hs51_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)- ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)- ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)- ::: Nil- , 4.0 ::: 0.0 ::: 0.0 ::: Nil- )--run_hs51 :: Result N5-run_hs51 = levmar hs51- (Just hs51_jac)- hs51_params- hs51_samples- 1000- opts- Nothing- Nothing- (Just hs51_linear_constraints)- Nothing------------------------------------------------------------------------------------- Hock - Schittkowski problem 01 (box constrained),--- minimum at (1.0, 1.0)------ constr1: p1 >= -1.5--hs01 :: Model N2 Double-hs01 p0 p1 = [ 10.0*(p1 - p0**2)- , 1.0 - p0- ]--hs01_jac :: Jacobian N2 Double-hs01_jac p0 _ = [ -20.0*p0 ::: 10.0 ::: Nil- , -1.0 ::: 0.0 ::: Nil- ]--hs01_n :: Int-hs01_n = 2--hs01_params :: SizedList N2 Double-hs01_params = -2.0 ::: 1.0 ::: Nil--hs01_samples :: [Double]-hs01_samples = replicate hs01_n 0.0--hs01_lb, hs01_ub :: SizedList N2 Double-hs01_lb = -_DBL_MAX ::: -1.5 ::: Nil-hs01_ub = _DBL_MAX ::: _DBL_MAX ::: Nil--_DBL_MAX :: Double-_DBL_MAX = 1e+37 -- TODO: Get this directly from <float.h>.--run_hs01 :: Result N2-run_hs01 = levmar hs01- (Just hs01_jac)- hs01_params- hs01_samples- 1000- opts- (Just hs01_lb)- (Just hs01_ub)- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Hock - Schittkowski MODIFIED problem 21 (box constrained),--- minimum at (2.0, 0.0)------ constr1: 2 <= p0 <=50--- constr2: -50 <= p1 <=50------ Original HS21 has the additional constraint 10*p0 - p1 >= 10--- which is inactive at the solution, so it is dropped here.--hs21 :: Model N2 Double-hs21 p0 p1 = [p0 / 10.0, p1]--hs21_jac :: Jacobian N2 Double-hs21_jac _ _ = [ 0.1 ::: 0.0 ::: Nil- , 0.0 ::: 1.0 ::: Nil- ]--hs21_n :: Int-hs21_n = 2--hs21_params :: SizedList N2 Double-hs21_params = -1.0 ::: -1.0 ::: Nil--hs21_samples :: [Double]-hs21_samples = replicate hs21_n 0.0--hs21_lb, hs21_ub :: SizedList N2 Double-hs21_lb = 2.0 ::: -50.0 ::: Nil-hs21_ub = 50.0 ::: 50.0 ::: Nil--run_hs21 :: Result N2-run_hs21 = levmar hs21- (Just hs21_jac)- hs21_params- hs21_samples- 1000- opts- (Just hs21_lb)- (Just hs21_ub)- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Problem hatfldb (box constrained),--- minimum at (0.947214, 0.8, 0.64, 0.4096)------ constri: pi >= 0.0 (i=1..4)--- constr5: p1 <= 0.8--hatfldb :: Model N4 Double-hatfldb p0 p1 p2 p3 = [ p0 - 1.0- , p0 - sqrt p1- , p1 - sqrt p2- , p2 - sqrt p3- ]--hatfldb_jac :: Jacobian N4 Double-hatfldb_jac _ p1 p2 p3 = [ 1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil- , 1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil- , 0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil- , 0.0 ::: 0.0 ::: 1.0 ::: -0.5 / sqrt p3 ::: Nil- ]--hatfldb_n :: Int-hatfldb_n = 4--hatfldb_params :: SizedList N4 Double-hatfldb_params = 0.1 ::: 0.1 ::: 0.1 ::: 0.1 ::: Nil--hatfldb_samples :: [Double]-hatfldb_samples = replicate hatfldb_n 0.0--hatfldb_lb, hatfldb_ub :: SizedList N4 Double-hatfldb_lb = 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil-hatfldb_ub = _DBL_MAX ::: 0.8 ::: _DBL_MAX ::: _DBL_MAX ::: Nil--run_hatfldb :: Result N4-run_hatfldb = levmar hatfldb- (Just hatfldb_jac)- hatfldb_params- hatfldb_samples- 1000- opts- (Just hatfldb_lb)- (Just hatfldb_ub)- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Problem hatfldc (box constrained),--- minimum at (1.0, 1.0, 1.0, 1.0)------ constri: pi >= 0.0 (i=1..4)--- constri+4: pi <= 10.0 (i=1..4)--hatfldc :: Model N4 Double-hatfldc p0 p1 p2 p3 = [ p0 - 1.0- , p0 - sqrt p1- , p1 - sqrt p2- , p3 - 1.0- ]--hatfldc_jac :: Jacobian N4 Double-hatfldc_jac _ p1 p2 _ = [ 1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil- , 1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil- , 0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil- , 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil- ]--hatfldc_n :: Int-hatfldc_n = 4--hatfldc_params :: SizedList N4 Double-hatfldc_params = 0.9 ::: 0.9 ::: 0.9 ::: 0.9 ::: Nil--hatfldc_samples :: [Double]-hatfldc_samples = replicate hatfldc_n 0.0--hatfldc_lb, hatfldc_ub :: SizedList N4 Double-hatfldc_lb = 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil-hatfldc_ub = 10.0 ::: 10.0 ::: 10.0 ::: 10.0 ::: Nil--run_hatfldc :: Result N4-run_hatfldc = levmar hatfldc- (Just hatfldc_jac)- hatfldc_params- hatfldc_samples- 1000- opts- (Just hatfldc_lb)- (Just hatfldc_ub)- noLinearConstraints- Nothing------------------------------------------------------------------------------------- Hock - Schittkowski (modified) problem 52 (box/linearly constrained),--- minimum at (-0.09, 0.03, 0.25, -0.19, 0.03)------ constr1: p0 + 3*p1 = 0--- constr2: p2 + p3 - 2*p4 = 0--- constr3: p1 - p4 = 0------ To the above 3 constraints, we add the following 5:--- constr4: -0.09 <= p0--- constr5: 0.0 <= p1 <= 0.3--- constr6: p2 <= 0.25--- constr7: -0.2 <= p3 <= 0.3--- constr8: 0.0 <= p4 <= 0.3--modhs52 :: Model N5 Double-modhs52 p0 p1 p2 p3 p4 = [ 4.0*p0 - p1- , p1 + p2 - 2.0- , p3 - 1.0- , p4 - 1.0- ]--modhs52_jac :: Jacobian N5 Double-modhs52_jac _ _ _ _ _ = [ 4.0 ::: -1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil- , 0.0 ::: 1.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: Nil- , 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: 0.0 ::: Nil- , 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil- ]--modhs52_n :: Int-modhs52_n = 4--modhs52_params :: SizedList N5 Double-modhs52_params = 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: Nil--modhs52_samples :: [Double]-modhs52_samples = replicate modhs52_n 0.0--modhs52_linear_constraints :: LinearConstraints N3 N5 Double-modhs52_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)- ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)- ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)- ::: Nil- , 0.0 ::: 0.0 ::: 0.0 ::: Nil- )--modhs52_weights :: SizedList N5 Double-modhs52_weights = 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: Nil--modhs52_lb, modhs52_ub :: SizedList N5 Double-modhs52_lb = -0.09 ::: 0.0 ::: -_DBL_MAX ::: -0.2 ::: 0.0 ::: Nil-modhs52_ub = _DBL_MAX ::: 0.3 ::: 0.25 ::: 0.3 ::: 0.3 ::: Nil--run_modhs52 :: Result N5-run_modhs52 = levmar modhs52- (Just modhs52_jac)- modhs52_params- modhs52_samples- 1000- opts- (Just modhs52_lb)- (Just modhs52_ub)- (Just modhs52_linear_constraints)- (Just modhs52_weights)------------------------------------------------------------------------------------- Schittkowski (modified) problem 235 (box/linearly constrained),--- minimum at (-1.725, 2.9, 0.725)------ constr1: p0 + p2 = -1.0;------ To the above constraint, we add the following 2:--- constr2: p1 - 4*p2 = 0--- constr3: 0.1 <= p1 <= 2.9--- constr4: 0.7 <= p2--mods235 :: Model N3 Double-mods235 p0 p1 _ = [ 0.1*(p0 - 1.0)- , p1 - p0**2- ]--mods235_jac :: Jacobian N3 Double-mods235_jac p0 _ _ = [ 0.1 ::: 0.0 ::: 0.0 ::: Nil- , -2.0*p0 ::: 1.0 ::: 0.0 ::: Nil- ]--mods235_n :: Int-mods235_n = 2--mods235_params :: SizedList N3 Double-mods235_params = -2.0 ::: 3.0 ::: 1.0 ::: Nil--mods235_samples :: [Double]-mods235_samples = replicate mods235_n 0.0--mods235_linear_constraints :: LinearConstraints N2 N3 Double-mods235_linear_constraints = ( (1.0 ::: 0.0 ::: 1.0 ::: Nil)- ::: (0.0 ::: 1.0 ::: -4.0 ::: Nil)- ::: Nil- , -1.0 ::: 0.0 ::: Nil- )--mods235_lb, mods235_ub :: SizedList N3 Double-mods235_lb = -_DBL_MAX ::: 0.1 ::: 0.7 ::: Nil-mods235_ub = _DBL_MAX ::: 2.9 ::: _DBL_MAX ::: Nil--run_mods235 :: Result N3-run_mods235 = levmar mods235- (Just mods235_jac)- mods235_params- mods235_samples- 1000- opts- (Just mods235_lb)- (Just mods235_ub)- (Just mods235_linear_constraints)- Nothing------------------------------------------------------------------------------------- Boggs and Tolle modified problem 7 (box/linearly constrained),--- minimum at (0.7, 0.49, 0.19, 1.19, -0.2)------ We keep the original objective function & starting point and use the--- following constraints:------ subject to cons1:--- x[1]+x[2] - x[3] = 1.0;--- subject to cons2:--- x[2] - x[4] + x[1] = 0.0;--- subject to cons3:--- x[5] + x[1] = 0.5;--- subject to cons4:--- x[5]>=-0.3;--- subject to cons5:--- x[1]<=0.7;--modbt7 :: Model N5 Double-modbt7 p0 p1 _ _ _ = replicate modbt7_n (100.0*m**2 + n**2)- where- m = p1 - p0**2- n = p0 - 1.0--modbt7_jac :: Jacobian N5 Double-modbt7_jac p0 p1 _ _ _ = replicate modbt7_n- ( -400.0*m*p0 + 2.0*p0 - 2.0- ::: 200.0*m- ::: 0.0- ::: 0.0- ::: 0.0- ::: Nil- )- where- m = p1 - p0**2--modbt7_n :: Int-modbt7_n = 5--modbt7_params :: SizedList N5 Double-modbt7_params = -2.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: Nil--modbt7_samples :: [Double]-modbt7_samples = replicate modbt7_n 0.0--modbt7_linear_constraints :: LinearConstraints N3 N5 Double-modbt7_linear_constraints = ( (1.0 ::: 1.0 ::: -1.0 ::: 0.0 ::: 0.0 ::: Nil)- ::: (1.0 ::: 1.0 ::: 0.0 ::: -1.0 ::: 0.0 ::: Nil)- ::: (1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil)- ::: Nil- , 1.0 ::: 0.0 ::: 0.5 ::: Nil- )--modbt7_lb, modbt7_ub :: SizedList N5 Double-modbt7_lb = -_DBL_MAX ::: -_DBL_MAX ::: -_DBL_MAX ::: -_DBL_MAX ::: -0.3 ::: Nil-modbt7_ub = 0.7 ::: _DBL_MAX ::: _DBL_MAX ::: _DBL_MAX ::: _DBL_MAX ::: Nil--run_modbt7 :: Result N5-run_modbt7 = levmar modbt7- (Just modbt7_jac)- modbt7_params- modbt7_samples- 1000- opts- (Just modbt7_lb)- (Just modbt7_ub)- (Just modbt7_linear_constraints)- Nothing---- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--- !! TODO: This returns with: infStopReason = MaxIterations !!--- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!------------------------------------------------------------------------------------- Equilibrium combustion problem, constrained nonlinear equation from the book--- by Floudas et al.------ Minimum at (0.0034, 31.3265, 0.0684, 0.8595, 0.0370)------ constri: pi>=0.0001 (i=1..5)--- constri+5: pi<=100.0 (i=1..5)--combust :: Model N5 Double-combust p0 p1 p2 p3 p4 =- [ p0*p1 + p0 - 3*p4- , 2*p0*p1 + p0 + 3*r10*p1*p1 + p1*p2*p2 + r7*p1*p2 + r9*p1*p3 + r8*p1 - r*p4- , 2*p1*p2*p2 + r7*p1*p2 + 2*r5*p2*p2 + r6*p2-8*p4- , r9*p1*p3 + 2*p3*p3 - 4*r*p4- , p0*p1 + p0 + r10*p1*p1 + p1*p2*p2 + r7*p1*p2 + r9*p1*p3 + r8*p1 + r5*p2*p2 + r6*p2 + p3*p3 - 1.0- ]--r, r5, r6, r7, r8, r9, r10 :: Double-r = 10-r5 = 0.193-r6 = 4.10622*1e-4-r7 = 5.45177*1e-4-r8 = 4.4975 *1e-7-r9 = 3.40735*1e-5-r10 = 9.615 *1e-7--combust_jac :: Jacobian N5 Double-combust_jac p0 p1 p2 p3 _ =- [ p1 + 1- ::: p0- ::: 0.0- ::: 0.0- ::: -3- ::: Nil- , 2*p1 + 1- ::: 2*p0 + 6*r10*p1 + p2*p2 + r7*p2 + r9*p3 + r8- ::: 2*p1*p2 + r7*p1- ::: r9*p1- ::: -r- ::: Nil- , 0.0- ::: 2*p2*p2 + r7*p2- ::: 4*p1*p2 + r7*p1 + 4*r5*p2 + r6- ::: 0.0- ::: -8- ::: Nil- , 0.0- ::: r9*p3- ::: 0.0- ::: r9*p1 + 4*p3- ::: -4*r- ::: Nil- , p1 + 1- ::: p0 + 2*r10*p1 + p2*p2 + r7*p2 + r9*p3 + r8- ::: 2*p1*p2 + r7*p1 + 2*r5*p2 + r6- ::: r9*p1 + 2*p3- ::: 0.0- ::: Nil- ]--combust_n :: Int-combust_n = 5--combust_params :: SizedList N5 Double-combust_params = 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: Nil--combust_samples :: [Double]-combust_samples = replicate combust_n 0.0--combust_lb, combust_ub :: SizedList N5 Double-combust_lb = 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: Nil-combust_ub = 100.0 ::: 100.0 ::: 100.0 ::: 100.0 ::: 100.0 ::: Nil--run_combust :: Result N5-run_combust = levmar combust- (Just combust_jac)- combust_params- combust_samples- 1000- opts- (Just combust_lb)- (Just combust_ub)- noLinearConstraints- Nothing+module Main where++import LevMar ( levmar++ , Model+ , Jacobian++ , Options(..), defaultOpts++ , LinearConstraints, noLinearConstraints++ , LevMarError++ , Info(..), CovarMatrix++ , S, Z+ , SizedList(..)+ )+++import qualified LevMar.AD as AD+import qualified LevMar.Fitting as Fitting+import qualified LevMar.Fitting.AD as Fitting.AD++import qualified SizedList as SL (replicate)+++--------------------------------------------------------------------------------++type Result n = Either LevMarError+ ( SizedList n Double+ , Info Double+ , CovarMatrix n Double+ )++printInteresting :: Result n -> IO ()+printInteresting (Left err) = putStrLn ("Error: " ++ show err)+printInteresting (Right (ps, inf, covar)) =+ do putStrLn ("infStopReason = " ++ show (infStopReason inf))+ putStrLn ("infNorm2E = " ++ show (infNorm2E inf))+ putStrLn ("infNumIter = " ++ show (infNumIter inf))+ putStrLn ("ps = " ++ show ps)++sqr :: Num a => a -> a+sqr x = x*x++--------------------------------------------------------------------------------+-- Handy type synonyms for type-level naturals:++type N0 = Z+type N1 = S N0+type N2 = S N1+type N3 = S N2+type N4 = S N3+type N5 = S N4+type N6 = S N5++--------------------------------------------------------------------------------+-- Default options:++opts :: Options Double+opts = defaultOpts { optStopNormInfJacTe = 1e-15+ , optStopNorm2Dp = 1e-15+ , optStopNorm2E = 1e-20+ }++--------------------------------------------------------------------------------+-- Rosenbrock function,+-- global minimum at (1, 1)++ros :: Floating r => Model N2 N2 r+ros p0 p1 = SL.replicate (sqr (1.0 - p0) + ros_d*sqr m)+ where+ m = p1 - sqr p0++ros_jac :: Floating r => Jacobian N2 N2 r+ros_jac p0 p1 = SL.replicate ( -2 + 2*p0 - 4*ros_d*m*p0+ ::: 2*ros_d*m+ ::: Nil+ )+ where+ m = p1 - sqr p0++ros_d :: Floating r => r+ros_d = 105.0++ros_params :: Floating r => SizedList N2 r+ros_params = -1.2 ::: 1.0 ::: Nil++ros_samples :: Floating r => SizedList N2 r+ros_samples = SL.replicate 0.0++-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+-- !! TODO: Find out why these return with: infStopReason = MaxIterations !!+-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!++run_ros :: IO ()+run_ros = printInteresting $+ levmar ros+ Nothing+ ros_params+ ros_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++run_ros_jac :: IO ()+run_ros_jac = printInteresting $+ levmar ros+ (Just ros_jac)+ ros_params+ ros_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++run_ros_autojac :: IO ()+run_ros_autojac = printInteresting $+ AD.levmar ros+ ros_params+ ros_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Modified Rosenbrock problem,+-- global minimum at (1, 1)++modros :: Floating r => Model N2 N3 r+modros p0 p1 = 10*(p1 - sqr p0)+ ::: 1.0 - p0+ ::: modros_lam+ ::: Nil++modros_jac :: Floating r => Jacobian N2 N3 r+modros_jac p0 _ = (-20*p0 ::: 10.0 ::: Nil)+ ::: (-1.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: Nil)+ ::: Nil++modros_lam :: Floating r => r+modros_lam = 1e02++modros_params :: Floating r => SizedList N2 r+modros_params = -1.2 ::: 1.0 ::: Nil++modros_samples :: Floating r => SizedList N3 r+modros_samples = SL.replicate 0.0++run_modros :: IO ()+run_modros = printInteresting $+ levmar modros+ Nothing+ modros_params+ modros_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++run_modros_jac :: IO ()+run_modros_jac = printInteresting $+ levmar modros+ (Just modros_jac)+ modros_params+ modros_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++run_modros_autojac :: IO ()+run_modros_autojac = printInteresting $+ AD.levmar modros+ modros_params+ modros_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Powell's function,+-- minimum at (0, 0)++powell :: Floating r => Model N2 N2 r+powell p0 p1 = p0+ ::: 10.0*p0 / m + 2*sqr p1+ ::: Nil+ where+ m = p0 + 0.1++powell_jac :: Floating r => Jacobian N2 N2 r+powell_jac p0 p1 = (1.0 ::: 0.0 ::: Nil)+ ::: (1.0 / sqr m ::: 4.0*p1 ::: Nil)+ ::: Nil+ where+ m = p0 + 0.1++powell_params :: Floating r => SizedList N2 r+powell_params = -1.2 ::: 1.0 ::: Nil++powell_samples :: Floating r => SizedList N2 r+powell_samples = SL.replicate 0.0++run_powell :: IO ()+run_powell = printInteresting $+ levmar powell+ Nothing+ powell_params+ powell_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++run_powell_jac :: IO ()+run_powell_jac = printInteresting $+ levmar powell+ (Just powell_jac)+ powell_params+ powell_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+-- !! TODO: Here the automatic jacobian does not seem right because !!+-- !! infNorm2E is very high compared to the manual jacobian! !!+-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!++run_powell_autojac :: IO ()+run_powell_autojac = printInteresting $+ AD.levmar powell+ powell_params+ powell_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Wood's function,+-- minimum at (1, 1, 1, 1)++wood :: Floating r => Model N4 N6 r+wood p0 p1 p2 p3 = 10.0*(p1 - sqr p0)+ ::: 1.0 - p0+ ::: sqrt 90.0*(p3 - sqr p2)+ ::: 1.0 - p2+ ::: sqrt 10.0*(p1 + p3 - 2.0)+ ::: (p1 - p3) / sqrt 10.0+ ::: Nil++wood_params :: Floating r => SizedList N4 r+wood_params = -3.0 ::: -1.0 ::: -3.0 ::: -1.0 ::: Nil++wood_samples :: Floating r => SizedList N6 r+wood_samples = SL.replicate 0.0++run_wood :: IO ()+run_wood = printInteresting $+ levmar wood+ Nothing+ wood_params+ wood_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++run_wood_autojac :: IO ()+run_wood_autojac = printInteresting $+ AD.levmar wood+ wood_params+ wood_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Meyer's (reformulated) data fitting problem,+-- minimum at (2.48, 6.18, 3.45)++meyer :: Floating r => Fitting.SimpleModel N3 r+meyer p0 p1 p2 x = p0*exp (10.0*p1 / (ui + p2) - 13.0)+ where+ ui = 0.45 + 0.05*x++meyer_jac :: Floating r => Fitting.SimpleJacobian N3 r+meyer_jac p0 p1 p2 x = tmp+ ::: 10.0*p0*tmp / (ui + p2)+ ::: -10.0*p0*p1*tmp / ((ui + p2)*(ui + p2))+ ::: Nil+ where+ tmp = exp (10.0*p1 / (ui + p2) - 13.0)+ ui = 0.45 + 0.05*x++meyer_params :: Floating r => SizedList N3 r+meyer_params = 8.85 ::: 4.0 ::: 2.5 ::: Nil++-- TODO: Unfortunately 'zip [0..] ...' won't work because (:~>)+-- doesn't have an Enum instance:+meyer_samples :: (Num a, Floating r) => [(a, r)]+meyer_samples = [ ( 0, 34.780)+ , ( 1, 28.610)+ , ( 2, 23.650)+ , ( 3, 19.630)+ , ( 4, 16.370)+ , ( 5, 13.720)+ , ( 6, 11.540)+ , ( 7, 9.744)+ , ( 8, 8.261)+ , ( 9, 7.030)+ , (10, 6.005)+ , (11, 5.147)+ , (12, 4.427)+ , (13, 3.820)+ , (14, 3.307)+ , (15, 2.872)+ ]++run_meyer :: IO ()+run_meyer = printInteresting $+ Fitting.levmar meyer+ Nothing+ meyer_params+ meyer_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++run_meyer_jac :: IO ()+run_meyer_jac = printInteresting $+ Fitting.levmar meyer+ (Just meyer_jac)+ meyer_params+ meyer_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+-- !! TODO: Here the automatic jacobian does not seem right because !!+-- !! infNorm2E is very high compared to the manual jacobian! !!+-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!++run_meyer_autojac :: IO ()+run_meyer_autojac = printInteresting $+ Fitting.AD.levmar meyer+ meyer_params+ meyer_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- helical valley function,+-- minimum at (1.0, 0.0, 0.0)++helval :: (Ord r, Floating r) => Model N3 N3 r+helval p0 p1 p2 = 10.0*(p2 - 10.0*theta)+ ::: 10.0*sqrt tmp - 1.0+ ::: p2+ ::: Nil+ where+ m = atan (p1 / p0) / (2.0*pi)++ tmp = sqr p0 + sqr p1++ theta | p0 < 0.0 = m + 0.5+ | 0.0 < p0 = m+ | p1 >= 0 = 0.25+ | otherwise = -0.25++heval_jac :: Floating r => Jacobian N3 N3 r+heval_jac p0 p1 _ = (50.0*p1 / (pi*tmp) ::: -50.0*p0 / (pi*tmp) ::: 10.0 ::: Nil)+ ::: (10.0*p0 / sqrt tmp ::: 10.0*p1 / sqrt tmp ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 1.0 ::: Nil)+ ::: Nil+ where+ tmp = sqr p0 + sqr p1++helval_params :: Floating r => SizedList N3 r+helval_params = -1.0 ::: 0.0 ::: 0.0 ::: Nil++helval_samples :: Floating r => SizedList N3 r+helval_samples = SL.replicate 0.0++run_helval :: IO ()+run_helval = printInteresting $+ levmar helval+ Nothing+ helval_params+ helval_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++run_helval_jac :: IO ()+run_helval_jac = printInteresting $+ levmar helval+ (Just heval_jac)+ helval_params+ helval_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+-- !! TODO: This function exits with the following error: !!+-- !! <interactive>: (==): No overloading for function !!+-- !! <interactive>: interrupted !!+-- !! <interactive>: warning: too many hs_exit()s !!+-- !! !!+-- !! Process haskell exited abnormally with code 252 !!+-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!++run_helval_autojac :: IO ()+run_helval_autojac = printInteresting $+ AD.levmar helval+ helval_params+ helval_samples+ 1000+ opts+ Nothing+ Nothing+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Boggs - Tolle problem 3 (linearly constrained),+-- minimum at (-0.76744, 0.25581, 0.62791, -0.11628, 0.25581)+--+-- constr1: p0 + 3*p1 = 0+-- constr2: p2 + p3 - 2*p4 = 0+-- constr3: p1 - p4 = 0++bt3 :: Floating r => Model N5 N5 r+bt3 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1+ + sqr t2+ + sqr t3+ + sqr t4+ )+ where+ t1 = p0 - p1+ t2 = p1 + p2 - 2.0+ t3 = p3 - 1.0+ t4 = p4 - 1.0++bt3_jac :: Floating r => Jacobian N5 N5 r+bt3_jac p0 p1 p2 p3 p4 = SL.replicate ( 2.0*t1+ ::: 2.0*(t2 - t1)+ ::: 2.0*t2+ ::: 2.0*t3+ ::: 2.0*t4+ ::: Nil+ )+ where+ t1 = p0 - p1+ t2 = p1 + p2 - 2.0+ t3 = p3 - 1.0+ t4 = p4 - 1.0++bt3_params :: Floating r => SizedList N5 r+bt3_params = 2.0 ::: 2.0 ::: 2.0 :::2.0 ::: 2.0 ::: Nil++bt3_samples :: Floating r => SizedList N5 r+bt3_samples = SL.replicate 0.0++bt3_linear_constraints :: Floating r => LinearConstraints N3 N5 r+bt3_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)+ ::: Nil+ , 0.0 ::: 0.0 ::: 0.0 ::: Nil+ )++run_bt3 :: IO ()+run_bt3 = printInteresting $+ levmar bt3+ Nothing+ bt3_params+ bt3_samples+ 1000+ opts+ Nothing+ Nothing+ (Just bt3_linear_constraints)+ Nothing++run_bt3_jac :: IO ()+run_bt3_jac = printInteresting $+ levmar bt3+ (Just bt3_jac)+ bt3_params+ bt3_samples+ 1000+ opts+ Nothing+ Nothing+ (Just bt3_linear_constraints)+ Nothing++run_bt3_autojac :: IO ()+run_bt3_autojac = printInteresting $+ AD.levmar bt3+ bt3_params+ bt3_samples+ 1000+ opts+ Nothing+ Nothing+ (Just bt3_linear_constraints)+ Nothing++--------------------------------------------------------------------------------+-- Hock - Schittkowski problem 28 (linearly constrained),+-- minimum at (0.5, -0.5, 0.5)+--+-- constr1: p0 + 2*p1 + 3*p2 = 1++hs28 :: Floating r => Model N3 N3 r+hs28 p0 p1 p2 = SL.replicate ( sqr t1+ + sqr t2+ )+ where+ t1 = p0 + p1+ t2 = p1 + p2++hs28_jac :: Floating r => Jacobian N3 N3 r+hs28_jac p0 p1 p2 = SL.replicate ( 2.0*t1+ ::: 2.0*(t1 + t2)+ ::: 2.0*t2+ ::: Nil+ )+ where+ t1 = p0 + p1+ t2 = p1 + p2++hs28_params :: Floating r => SizedList N3 r+hs28_params = -4.0 ::: 1.0 ::: 1.0 ::: Nil++hs28_samples :: Floating r => SizedList N3 r+hs28_samples = SL.replicate 0.0++hs28_linear_constraints :: Floating r => LinearConstraints N1 N3 r+hs28_linear_constraints = ( ((1.0 ::: 2.0 ::: 3.0 ::: Nil) ::: Nil)+ , 1.0 ::: Nil+ )++run_hs28 :: IO ()+run_hs28 = printInteresting $+ levmar hs28+ Nothing+ hs28_params+ hs28_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs28_linear_constraints)+ Nothing++run_hs28_jac :: IO ()+run_hs28_jac = printInteresting $+ levmar hs28+ (Just hs28_jac)+ hs28_params+ hs28_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs28_linear_constraints)+ Nothing++run_hs28_autojac :: IO ()+run_hs28_autojac = printInteresting $+ AD.levmar hs28+ hs28_params+ hs28_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs28_linear_constraints)+ Nothing++--------------------------------------------------------------------------------+-- Hock - Schittkowski problem 48 (linearly constrained),+-- minimum at (1.0, 1.0, 1.0, 1.0, 1.0)+--+-- constr1: sum [p0, p1, p2, p3, p4] = 5+-- constr2: p2 - 2*(p3 + p4) = -3++hs48 :: Floating r => Model N5 N5 r+hs48 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1+ + sqr t2+ + sqr t3+ )+ where+ t1 = p0 - 1.0+ t2 = p1 - p2+ t3 = p3 - p4++hs48_jac :: Floating r => Jacobian N5 N5 r+hs48_jac p0 p1 p2 p3 p4 = SL.replicate ( 2.0*t1+ ::: 2.0*t2+ ::: -2.0*t2+ ::: 2.0*t3+ ::: -2.0*t3+ ::: Nil+ )+ where+ t1 = p0 - 1.0+ t2 = p1 - p2+ t3 = p3 - p4++hs48_params :: Floating r => SizedList N5 r+hs48_params = 3.0 ::: 5.0 ::: -3.0 ::: 2.0 ::: -2.0 ::: Nil++hs48_samples :: Floating r => SizedList N5 r+hs48_samples = SL.replicate 0.0++hs48_linear_constraints :: Floating r => LinearConstraints N2 N5 r+hs48_linear_constraints = ( (1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 1.0 ::: -2.0 ::: -2.0 ::: Nil)+ ::: Nil+ , 5.0 ::: -3.0 ::: Nil+ )++run_hs48 :: IO ()+run_hs48 = printInteresting $+ levmar hs48+ Nothing+ hs48_params+ hs48_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs48_linear_constraints)+ Nothing++run_hs48_jac :: IO ()+run_hs48_jac = printInteresting $+ levmar hs48+ (Just hs48_jac)+ hs48_params+ hs48_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs48_linear_constraints)+ Nothing++run_hs48_autojac :: IO ()+run_hs48_autojac = printInteresting $+ AD.levmar hs48+ hs48_params+ hs48_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs48_linear_constraints)+ Nothing++--------------------------------------------------------------------------------+-- Hock - Schittkowski problem 51 (linearly constrained),+-- minimum at (1.0, 1.0, 1.0, 1.0, 1.0)+--+-- constr1: p0 + 3*p1 = 4+-- constr2: p2 + p3 - 2*p4 = 0+-- constr3: p1 - p4 = 0++hs51 :: Floating r => Model N5 N5 r+hs51 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1+ + sqr t2+ + sqr t3+ + sqr t4+ )+ where+ t1 = p0 - p1+ t2 = p1 + p2 - 2.0+ t3 = p3 - 1.0+ t4 = p4 - 1.0++hs51_jac :: Floating r => Jacobian N5 N5 r+hs51_jac p0 p1 p2 p3 p4 = SL.replicate ( 2.0*t1+ ::: 2.0*(t2 - t1)+ ::: 2.0*t2+ ::: 2.0*t3+ ::: 2.0*t4+ ::: Nil+ )+ where+ t1 = p0 - p1+ t2 = p1 + p2 - 2.0+ t3 = p3 - 1.0+ t4 = p4 - 1.0++hs51_params :: Floating r => SizedList N5 r+hs51_params = 2.5 ::: 0.5 ::: 2.0 ::: -1.0 ::: 0.5 ::: Nil++hs51_samples :: Floating r => SizedList N5 r+hs51_samples = SL.replicate 0.0++hs51_linear_constraints :: Floating r => LinearConstraints N3 N5 r+hs51_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)+ ::: Nil+ , 4.0 ::: 0.0 ::: 0.0 ::: Nil+ )++run_hs51 :: IO ()+run_hs51 = printInteresting $+ levmar hs51+ Nothing+ hs51_params+ hs51_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs51_linear_constraints)+ Nothing++run_hs51_jac :: IO ()+run_hs51_jac = printInteresting $+ levmar hs51+ (Just hs51_jac)+ hs51_params+ hs51_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs51_linear_constraints)+ Nothing++run_hs51_autojac :: IO ()+run_hs51_autojac = printInteresting $+ AD.levmar hs51+ hs51_params+ hs51_samples+ 1000+ opts+ Nothing+ Nothing+ (Just hs51_linear_constraints)+ Nothing++--------------------------------------------------------------------------------+-- Hock - Schittkowski problem 01 (box constrained),+-- minimum at (1.0, 1.0)+--+-- constr1: p1 >= -1.5++hs01 :: Floating r => Model N2 N2 r+hs01 p0 p1 = 10.0*(p1 - sqr p0)+ ::: 1.0 - p0+ ::: Nil++hs01_jac :: Floating r => Jacobian N2 N2 r+hs01_jac p0 _ = (-20.0*p0 ::: 10.0 ::: Nil)+ ::: (-1.0 ::: 0.0 ::: Nil)+ ::: Nil++hs01_params :: Floating r => SizedList N2 r+hs01_params = -2.0 ::: 1.0 ::: Nil++hs01_samples :: Floating r => SizedList N2 r+hs01_samples = SL.replicate 0.0++hs01_lb, hs01_ub :: Floating r => SizedList N2 r+hs01_lb = -_DBL_MAX ::: -1.5 ::: Nil+hs01_ub = _DBL_MAX ::: _DBL_MAX ::: Nil++_DBL_MAX :: Floating r => r+_DBL_MAX = 1e+37 -- TODO: Get this directly from <float.h>.++run_hs01 :: IO ()+run_hs01 = printInteresting $+ levmar hs01+ Nothing+ hs01_params+ hs01_samples+ 1000+ opts+ (Just hs01_lb)+ (Just hs01_ub)+ noLinearConstraints+ Nothing++run_hs01_jac :: IO ()+run_hs01_jac = printInteresting $+ levmar hs01+ (Just hs01_jac)+ hs01_params+ hs01_samples+ 1000+ opts+ (Just hs01_lb)+ (Just hs01_ub)+ noLinearConstraints+ Nothing++run_hs01_autojac :: IO ()+run_hs01_autojac = printInteresting $+ AD.levmar hs01+ hs01_params+ hs01_samples+ 1000+ opts+ (Just hs01_lb)+ (Just hs01_ub)+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Hock - Schittkowski MODIFIED problem 21 (box constrained),+-- minimum at (2.0, 0.0)+--+-- constr1: 2 <= p0 <=50+-- constr2: -50 <= p1 <=50+--+-- Original HS21 has the additional constraint 10*p0 - p1 >= 10+-- which is inactive at the solution, so it is dropped here.++hs21 :: Floating r => Model N2 N2 r+hs21 p0 p1 = p0 / 10.0+ ::: p1+ ::: Nil++hs21_jac :: Floating r => Jacobian N2 N2 r+hs21_jac _ _ = (0.1 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: Nil)+ ::: Nil++hs21_params :: Floating r => SizedList N2 r+hs21_params = -1.0 ::: -1.0 ::: Nil++hs21_samples :: Floating r => SizedList N2 r+hs21_samples = SL.replicate 0.0++hs21_lb, hs21_ub :: Floating r => SizedList N2 r+hs21_lb = 2.0 ::: -50.0 ::: Nil+hs21_ub = 50.0 ::: 50.0 ::: Nil++run_hs21 :: IO ()+run_hs21 = printInteresting $+ levmar hs21+ Nothing+ hs21_params+ hs21_samples+ 1000+ opts+ (Just hs21_lb)+ (Just hs21_ub)+ noLinearConstraints+ Nothing++run_hs21_jac :: IO ()+run_hs21_jac = printInteresting $+ levmar hs21+ (Just hs21_jac)+ hs21_params+ hs21_samples+ 1000+ opts+ (Just hs21_lb)+ (Just hs21_ub)+ noLinearConstraints+ Nothing++run_hs21_autojac :: IO ()+run_hs21_autojac = printInteresting $+ AD.levmar hs21+ hs21_params+ hs21_samples+ 1000+ opts+ (Just hs21_lb)+ (Just hs21_ub)+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Problem hatfldb (box constrained),+-- minimum at (0.947214, 0.8, 0.64, 0.4096)+--+-- constri: pi >= 0.0 (i=1..4)+-- constr5: p1 <= 0.8++hatfldb :: Floating r => Model N4 N4 r+hatfldb p0 p1 p2 p3 = p0 - 1.0+ ::: p0 - sqrt p1+ ::: p1 - sqrt p2+ ::: p2 - sqrt p3+ ::: Nil++hatfldb_jac :: Floating r => Jacobian N4 N4 r+hatfldb_jac _ p1 p2 p3 = (1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 1.0 ::: -0.5 / sqrt p3 ::: Nil)+ ::: Nil++hatfldb_params :: Floating r => SizedList N4 r+hatfldb_params = 0.1 ::: 0.1 ::: 0.1 ::: 0.1 ::: Nil++hatfldb_samples :: Floating r => SizedList N4 r+hatfldb_samples = SL.replicate 0.0++hatfldb_lb, hatfldb_ub :: Floating r => SizedList N4 r+hatfldb_lb = 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil+hatfldb_ub = _DBL_MAX ::: 0.8 ::: _DBL_MAX ::: _DBL_MAX ::: Nil++run_hatfldb :: IO ()+run_hatfldb = printInteresting $+ levmar hatfldb+ Nothing+ hatfldb_params+ hatfldb_samples+ 1000+ opts+ (Just hatfldb_lb)+ (Just hatfldb_ub)+ noLinearConstraints+ Nothing++run_hatfldb_jac :: IO ()+run_hatfldb_jac = printInteresting $+ levmar hatfldb+ (Just hatfldb_jac)+ hatfldb_params+ hatfldb_samples+ 1000+ opts+ (Just hatfldb_lb)+ (Just hatfldb_ub)+ noLinearConstraints+ Nothing++run_hatfldb_autojac :: IO ()+run_hatfldb_autojac = printInteresting $+ AD.levmar hatfldb+ hatfldb_params+ hatfldb_samples+ 1000+ opts+ (Just hatfldb_lb)+ (Just hatfldb_ub)+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Problem hatfldc (box constrained),+-- minimum at (1.0, 1.0, 1.0, 1.0)+--+-- constri: pi >= 0.0 (i=1..4)+-- constri+4: pi <= 10.0 (i=1..4)++hatfldc :: Floating r => Model N4 N4 r+hatfldc p0 p1 p2 p3 = p0 - 1.0+ ::: p0 - sqrt p1+ ::: p1 - sqrt p2+ ::: p3 - 1.0+ ::: Nil++hatfldc_jac :: Floating r => Jacobian N4 N4 r+hatfldc_jac _ p1 p2 _ = (1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil)+ ::: Nil++hatfldc_params :: Floating r => SizedList N4 r+hatfldc_params = 0.9 ::: 0.9 ::: 0.9 ::: 0.9 ::: Nil++hatfldc_samples :: Floating r => SizedList N4 r+hatfldc_samples = SL.replicate 0.0++hatfldc_lb, hatfldc_ub :: Floating r => SizedList N4 r+hatfldc_lb = 0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil+hatfldc_ub = 10.0 ::: 10.0 ::: 10.0 ::: 10.0 ::: Nil++run_hatfldc :: IO ()+run_hatfldc = printInteresting $+ levmar hatfldc+ Nothing+ hatfldc_params+ hatfldc_samples+ 1000+ opts+ (Just hatfldc_lb)+ (Just hatfldc_ub)+ noLinearConstraints+ Nothing++run_hatfldc_jac :: IO ()+run_hatfldc_jac = printInteresting $+ levmar hatfldc+ (Just hatfldc_jac)+ hatfldc_params+ hatfldc_samples+ 1000+ opts+ (Just hatfldc_lb)+ (Just hatfldc_ub)+ noLinearConstraints+ Nothing++run_hatfldc_autojac :: IO ()+run_hatfldc_autojac = printInteresting $+ AD.levmar hatfldc+ hatfldc_params+ hatfldc_samples+ 1000+ opts+ (Just hatfldc_lb)+ (Just hatfldc_ub)+ noLinearConstraints+ Nothing++--------------------------------------------------------------------------------+-- Hock - Schittkowski (modified) problem 52 (box/linearly constrained),+-- minimum at (-0.09, 0.03, 0.25, -0.19, 0.03)+--+-- constr1: p0 + 3*p1 = 0+-- constr2: p2 + p3 - 2*p4 = 0+-- constr3: p1 - p4 = 0+--+-- To the above 3 constraints, we add the following 5:+-- constr4: -0.09 <= p0+-- constr5: 0.0 <= p1 <= 0.3+-- constr6: p2 <= 0.25+-- constr7: -0.2 <= p3 <= 0.3+-- constr8: 0.0 <= p4 <= 0.3++modhs52 :: Floating r => Model N5 N4 r+modhs52 p0 p1 p2 p3 p4 = 4.0*p0 - p1+ ::: p1 + p2 - 2.0+ ::: p3 - 1.0+ ::: p4 - 1.0+ ::: Nil++modhs52_jac :: Floating r => Jacobian N5 N4 r+modhs52_jac _ _ _ _ _ = (4.0 ::: -1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil)+ ::: Nil++modhs52_params :: Floating r => SizedList N5 r+modhs52_params = 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: Nil++modhs52_samples :: Floating r => SizedList N4 r+modhs52_samples = SL.replicate 0.0++modhs52_linear_constraints :: Floating r => LinearConstraints N3 N5 r+modhs52_linear_constraints = ( (1.0 ::: 3.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)+ ::: Nil+ , 0.0 ::: 0.0 ::: 0.0 ::: Nil+ )++modhs52_weights :: Floating r => SizedList N5 r+modhs52_weights = 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: Nil++modhs52_lb, modhs52_ub :: Floating r => SizedList N5 r+modhs52_lb = -0.09 ::: 0.0 ::: -_DBL_MAX ::: -0.2 ::: 0.0 ::: Nil+modhs52_ub = _DBL_MAX ::: 0.3 ::: 0.25 ::: 0.3 ::: 0.3 ::: Nil++run_modhs52 :: IO ()+run_modhs52 = printInteresting $+ levmar modhs52+ Nothing+ modhs52_params+ modhs52_samples+ 1000+ opts+ (Just modhs52_lb)+ (Just modhs52_ub)+ (Just modhs52_linear_constraints)+ (Just modhs52_weights)++run_modhs52_jac :: IO ()+run_modhs52_jac = printInteresting $+ levmar modhs52+ (Just modhs52_jac)+ modhs52_params+ modhs52_samples+ 1000+ opts+ (Just modhs52_lb)+ (Just modhs52_ub)+ (Just modhs52_linear_constraints)+ (Just modhs52_weights)++run_modhs52_autojac :: IO ()+run_modhs52_autojac = printInteresting $+ AD.levmar modhs52+ modhs52_params+ modhs52_samples+ 1000+ opts+ (Just modhs52_lb)+ (Just modhs52_ub)+ (Just modhs52_linear_constraints)+ (Just modhs52_weights)++--------------------------------------------------------------------------------+-- Schittkowski (modified) problem 235 (box/linearly constrained),+-- minimum at (-1.725, 2.9, 0.725)+--+-- constr1: p0 + p2 = -1.0;+--+-- To the above constraint, we add the following 2:+-- constr2: p1 - 4*p2 = 0+-- constr3: 0.1 <= p1 <= 2.9+-- constr4: 0.7 <= p2++mods235 :: Floating r => Model N3 N2 r+mods235 p0 p1 _ = 0.1*(p0 - 1.0)+ ::: p1 - sqr p0+ ::: Nil++mods235_jac :: Floating r => Jacobian N3 N2 r+mods235_jac p0 _ _ = (0.1 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (-2.0*p0 ::: 1.0 ::: 0.0 ::: Nil)+ ::: Nil++mods235_params :: Floating r => SizedList N3 r+mods235_params = -2.0 ::: 3.0 ::: 1.0 ::: Nil++mods235_samples :: Floating r => SizedList N2 r+mods235_samples = SL.replicate 0.0++mods235_linear_constraints :: Floating r => LinearConstraints N2 N3 r+mods235_linear_constraints = ( (1.0 ::: 0.0 ::: 1.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: -4.0 ::: Nil)+ ::: Nil+ , -1.0 ::: 0.0 ::: Nil+ )++mods235_lb, mods235_ub :: Floating r => SizedList N3 r+mods235_lb = -_DBL_MAX ::: 0.1 ::: 0.7 ::: Nil+mods235_ub = _DBL_MAX ::: 2.9 ::: _DBL_MAX ::: Nil++run_mods235 :: IO ()+run_mods235 = printInteresting $+ levmar mods235+ Nothing+ mods235_params+ mods235_samples+ 1000+ opts+ (Just mods235_lb)+ (Just mods235_ub)+ (Just mods235_linear_constraints)+ Nothing++run_mods235_jac :: IO ()+run_mods235_jac = printInteresting $+ levmar mods235+ (Just mods235_jac)+ mods235_params+ mods235_samples+ 1000+ opts+ (Just mods235_lb)+ (Just mods235_ub)+ (Just mods235_linear_constraints)+ Nothing+++run_mods235_autojac :: IO ()+run_mods235_autojac = printInteresting $+ AD.levmar mods235+ mods235_params+ mods235_samples+ 1000+ opts+ (Just mods235_lb)+ (Just mods235_ub)+ (Just mods235_linear_constraints)+ Nothing++--------------------------------------------------------------------------------+-- Boggs and Tolle modified problem 7 (box/linearly constrained),+-- minimum at (0.7, 0.49, 0.19, 1.19, -0.2)+--+-- We keep the original objective function & starting point and use the+-- following constraints:+--+-- subject to cons1:+-- x[1]+x[2] - x[3] = 1.0;+-- subject to cons2:+-- x[2] - x[4] + x[1] = 0.0;+-- subject to cons3:+-- x[5] + x[1] = 0.5;+-- subject to cons4:+-- x[5]>=-0.3;+-- subject to cons5:+-- x[1]<=0.7;++modbt7 :: Floating r => Model N5 N5 r+modbt7 p0 p1 _ _ _ = SL.replicate (100.0*sqr m + sqr n)+ where+ m = p1 - sqr p0+ n = p0 - 1.0++modbt7_jac :: Floating r => Jacobian N5 N5 r+modbt7_jac p0 p1 _ _ _ = SL.replicate+ ( -400.0*m*p0 + 2.0*p0 - 2.0+ ::: 200.0*m+ ::: 0.0+ ::: 0.0+ ::: 0.0+ ::: Nil+ )+ where+ m = p1 - sqr p0++modbt7_params :: Floating r => SizedList N5 r+modbt7_params = -2.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: Nil++modbt7_samples :: Floating r => SizedList N5 r+modbt7_samples = SL.replicate 0.0++modbt7_linear_constraints :: Floating r => LinearConstraints N3 N5 r+modbt7_linear_constraints = ( (1.0 ::: 1.0 ::: -1.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (1.0 ::: 1.0 ::: 0.0 ::: -1.0 ::: 0.0 ::: Nil)+ ::: (1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil)+ ::: Nil+ , 1.0 ::: 0.0 ::: 0.5 ::: Nil+ )++modbt7_lb, modbt7_ub :: Floating r => SizedList N5 r+modbt7_lb = -_DBL_MAX ::: -_DBL_MAX ::: -_DBL_MAX ::: -_DBL_MAX ::: -0.3 ::: Nil+modbt7_ub = 0.7 ::: _DBL_MAX ::: _DBL_MAX ::: _DBL_MAX ::: _DBL_MAX ::: Nil++-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+-- !! TODO: Find out why these return with: infStopReason = MaxIterations !!+-- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!++run_modbt7 :: IO ()+run_modbt7 = printInteresting $+ levmar modbt7+ Nothing+ modbt7_params+ modbt7_samples+ 1000+ opts+ (Just modbt7_lb)+ (Just modbt7_ub)+ (Just modbt7_linear_constraints)+ Nothing++run_modbt7_jac :: IO ()+run_modbt7_jac = printInteresting $+ levmar modbt7+ (Just modbt7_jac)+ modbt7_params+ modbt7_samples+ 1000+ opts+ (Just modbt7_lb)+ (Just modbt7_ub)+ (Just modbt7_linear_constraints)+ Nothing++run_modbt7_autojac :: IO ()+run_modbt7_autojac = printInteresting $+ AD.levmar modbt7+ modbt7_params+ modbt7_samples+ 1000+ opts+ (Just modbt7_lb)+ (Just modbt7_ub)+ (Just modbt7_linear_constraints)+ Nothing++--------------------------------------------------------------------------------+-- Equilibrium combustion problem, constrained nonlinear equation from the book+-- by Floudas et al.+--+-- Minimum at (0.0034, 31.3265, 0.0684, 0.8595, 0.0370)+--+-- constri: pi>=0.0001 (i=1..5)+-- constri+5: pi<=100.0 (i=1..5)++combust :: Floating r => Model N5 N5 r+combust p0 p1 p2 p3 p4 =+ p0*p1 + p0 - 3*p4+ ::: 2*p0*p1 + p0 + 3*r10*p1*p1 + p1*p2*p2 + r7*p1*p2 + r9*p1*p3 + r8*p1 - r*p4+ ::: 2*p1*p2*p2 + r7*p1*p2 + 2*r5*p2*p2 + r6*p2-8*p4+ ::: r9*p1*p3 + 2*p3*p3 - 4*r*p4+ ::: p0*p1 + p0 + r10*p1*p1 + p1*p2*p2 + r7*p1*p2 + r9*p1*p3 + r8*p1 + r5*p2*p2 + r6*p2 + p3*p3 - 1.0+ ::: Nil++r, r5, r6, r7, r8, r9, r10 :: Floating r => r+r = 10+r5 = 0.193+r6 = 4.10622*1e-4+r7 = 5.45177*1e-4+r8 = 4.4975 *1e-7+r9 = 3.40735*1e-5+r10 = 9.615 *1e-7++combust_jac :: Floating r => Jacobian N5 N5 r+combust_jac p0 p1 p2 p3 _ =+ ( p1 + 1+ ::: p0+ ::: 0.0+ ::: 0.0+ ::: -3+ ::: Nil+ )+ ::: ( 2*p1 + 1+ ::: 2*p0 + 6*r10*p1 + p2*p2 + r7*p2 + r9*p3 + r8+ ::: 2*p1*p2 + r7*p1+ ::: r9*p1+ ::: -r+ ::: Nil+ )+ ::: ( 0.0+ ::: 2*p2*p2 + r7*p2+ ::: 4*p1*p2 + r7*p1 + 4*r5*p2 + r6+ ::: 0.0+ ::: -8+ ::: Nil+ )+ ::: ( 0.0+ ::: r9*p3+ ::: 0.0+ ::: r9*p1 + 4*p3+ ::: -4*r+ ::: Nil+ )+ ::: ( p1 + 1+ ::: p0 + 2*r10*p1 + p2*p2 + r7*p2 + r9*p3 + r8+ ::: 2*p1*p2 + r7*p1 + 2*r5*p2 + r6+ ::: r9*p1 + 2*p3+ ::: 0.0+ ::: Nil+ )+ ::: Nil++combust_params :: Floating r => SizedList N5 r+combust_params = 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: Nil++combust_samples :: Floating r => SizedList N5 r+combust_samples = SL.replicate 0.0++combust_lb, combust_ub :: Floating r => SizedList N5 r+combust_lb = 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: Nil+combust_ub = 100.0 ::: 100.0 ::: 100.0 ::: 100.0 ::: 100.0 ::: Nil++run_combust :: IO ()+run_combust = printInteresting $+ levmar combust+ Nothing+ combust_params+ combust_samples+ 1000+ opts+ (Just combust_lb)+ (Just combust_ub)+ noLinearConstraints+ Nothing++run_combust_jac :: IO ()+run_combust_jac = printInteresting $+ levmar combust+ (Just combust_jac)+ combust_params+ combust_samples+ 1000+ opts+ (Just combust_lb)+ (Just combust_ub)+ noLinearConstraints+ Nothing++run_combust_autojac :: IO ()+run_combust_autojac = printInteresting $+ AD.levmar combust+ combust_params+ combust_samples+ 1000+ opts+ (Just combust_lb)+ (Just combust_ub)+ noLinearConstraints+ Nothing+ -- The End ---------------------------------------------------------------------
LevMar.hs view
@@ -52,10 +52,18 @@ import qualified LevMar.Intermediate as LMA_I -import TypeLevelNat (Z, S, Nat)-import SizedList (SizedList(..), toList, unsafeFromList)-import NFunction (NFunction, ($*))+import LevMar.Utils ( LinearConstraints+ , noLinearConstraints+ , Matrix+ , CovarMatrix+ , convertLinearConstraints+ , convertResult+ ) +import TypeLevelNat ( Z, S, Nat )+import SizedList ( SizedList(..), toList, unsafeFromList )+import NFunction ( NFunction, ($*) )+ import Data.Either @@ -63,25 +71,28 @@ -- Model & Jacobian. -------------------------------------------------------------------------------- -{- | A function from @n@ parameters of type @r@ to a list of @r@.+{- | A functional relation describing measurements represented as a function+from @m@ parameters to @n@ expected measurements. An example from /Demo.hs/: @ type N4 = 'S' ('S' ('S' ('S' 'Z'))) -hatfldc :: Model N4 Double-hatfldc p0 p1 p2 p3 = [ p0 - 1.0- , p0 - sqrt p1- , p1 - sqrt p2- , p3 - 1.0- ]+hatfldc :: Model N4 N4 Double+hatfldc p0 p1 p2 p3 = p0 - 1.0+ ::: p0 - sqrt p1+ ::: p1 - sqrt p2+ ::: p3 - 1.0+ ::: Nil @ -}-type Model n r = NFunction n r [r]+type Model m n r = NFunction m r (SizedList n r) -{- | The jacobian of the 'Model' function. Expressed as a function from-@n@ parameters of type @r@ to a list of @n@-sized lists of @r@+{- | The jacobian of the 'Model' function. Expressed as a function+from @m@ parameters to a @n@/x/@m@ matrix which for each of the @n@+expected measurement describes the @m@ partial derivatives of the+parameters. See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant> @@ -90,16 +101,16 @@ @ type N4 = 'S' ('S' ('S' ('S' 'Z'))) -hatfldc_jac :: Jacobian N4 Double-hatfldc_jac _ p1 p2 _ = [ 1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil- , 1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil- , 0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil- , 0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil- ]+hatfldc_jac :: Jacobian N4 N4 Double+hatfldc_jac _ p1 p2 _ = (1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (1.0 ::: -0.5 / sqrt p1 ::: 0.0 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 1.0 ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil)+ ::: (0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil)+ ::: Nil @ -} -type Jacobian n r = NFunction n r [SizedList n r]+type Jacobian m n r = NFunction m r (Matrix n m r) --------------------------------------------------------------------------------@@ -107,58 +118,33 @@ -------------------------------------------------------------------------------- -- | The Levenberg-Marquardt algorithm.-levmar :: forall n k r. (Nat n, Nat k, LMA_I.LevMarable r)- => (Model n r) -- ^ Model- -> Maybe (Jacobian n r) -- ^ Optional jacobian- -> SizedList n r -- ^ Initial parameters- -> [r] -- ^ Samples+levmar :: forall m n k r. (Nat m, Nat n, Nat k, LMA_I.LevMarable r)+ => (Model m n r) -- ^ Model+ -> Maybe (Jacobian m n r) -- ^ Optional jacobian+ -> SizedList m r -- ^ Initial parameters+ -> SizedList n r -- ^ Samples -> Integer -- ^ Maximum number of iterations -> LMA_I.Options r -- ^ Minimization options- -> Maybe (SizedList n r) -- ^ Optional lower bounds- -> Maybe (SizedList n r) -- ^ Optional upper bounds- -> Maybe (LinearConstraints k n r) -- ^ Optional linear constraints- -> Maybe (SizedList n r) -- ^ Optional weights- -> Either LMA_I.LevMarError (SizedList n r, LMA_I.Info r, CovarMatrix n r)+ -> Maybe (SizedList m r) -- ^ Optional lower bounds+ -> Maybe (SizedList m r) -- ^ Optional upper bounds+ -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints+ -> Maybe (SizedList m r) -- ^ Optional weights+ -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r) levmar model mJac params ys itMax opts mLowBs mUpBs mLinC mWghts = fmap convertResult $ LMA_I.levmar (convertModel model) (fmap convertJacob mJac) (toList params)- ys+ (toList ys) itMax opts (fmap toList mLowBs) (fmap toList mUpBs)- (fmap convertLinC mLinC)+ (fmap convertLinearConstraints mLinC) (fmap toList mWghts) where- convertModel f = \ps -> f $* (unsafeFromList ps :: SizedList n r)- convertJacob f = \ps -> map toList ((f $* (unsafeFromList ps :: SizedList n r)) :: [SizedList n r])- convertLinC (cMat, rhcVec) = ( map toList $ toList cMat- , toList rhcVec- )- convertResult (psResult, info, covar) = ( unsafeFromList psResult- , info- , unsafeFromList $ map unsafeFromList covar- )---- | Linear constraints consisting of a constraints matrix, /kxn/ and--- a right hand constraints vector, /kx1/ where /n/ is the number of--- parameters and /k/ is the number of constraints.-type LinearConstraints k n r = (Matrix k n r, SizedList k r)---- |Value to denote the absense of any linear constraints over the--- parameters of the model function. Use this instead of 'Nothing'--- because the type parameter which contains the number of constraints--- can't be inferred.-noLinearConstraints :: Nat n => Maybe (LinearConstraints Z n r)-noLinearConstraints = Nothing---- | A /nxm/ matrix is a sized list of /n/ sized lists of length /m/.-type Matrix n m r = SizedList n (SizedList m r)---- | Covariance matrix corresponding to LS solution.-type CovarMatrix n r = Matrix n n r+ convertModel f = \ps -> toList (f $* (unsafeFromList ps :: SizedList m r) :: SizedList n r)+ convertJacob f = \ps -> toList (fmap toList (f $* (unsafeFromList ps :: SizedList m r) :: Matrix n m r)) -- The End ---------------------------------------------------------------------
+ LevMar/AD.hs view
@@ -0,0 +1,142 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}++--------------------------------------------------------------------------------+-- |+-- Module : LevMar.AD+-- Copyright : (c) 2009 Roel van Dijk & Bas van Dijk+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : vandijk.roel@gmail.com, v.dijk.bas@gmail.com+-- Stability : Experimental+--+-- A levmar variant that uses Automatic Differentiation to+-- automatically compute the Jacobian.+--+-- For additional documentation see the documentation of the levmar C+-- library which this library is based on:+-- <http://www.ics.forth.gr/~lourakis/levmar/>+--+--------------------------------------------------------------------------------+++module LevMar.AD+ ( -- * Model+ LMA.Model++ -- * Levenberg-Marquardt algorithm.+ , LMA_I.LevMarable+ , levmar++ , LinearConstraints+ , noLinearConstraints+ , Matrix++ -- * Minimization options.+ , LMA_I.Options(..)+ , LMA_I.defaultOpts++ -- * Output+ , LMA_I.Info(..)+ , LMA_I.StopReason(..)+ , CovarMatrix++ , LMA_I.LevMarError(..)++ -- *Type-level machinery+ , Z, S, Nat+ , SizedList(..)+ , NFunction+ )+ where+++import qualified LevMar as LMA+import qualified LevMar.Intermediate as LMA_I++import LevMar.Utils ( LinearConstraints+ , noLinearConstraints+ , Matrix+ , CovarMatrix+ , convertLinearConstraints+ , convertResult+ )++import TypeLevelNat ( Z, S, Nat )+import SizedList ( SizedList(..), toList, unsafeFromList )+import NFunction ( NFunction, ($*) )++import LevMar.Utils.AD ( value, firstDeriv, constant, idDAt )++-- From vector-space:+import Data.Derivative ( (:~>) )+import Data.VectorSpace ( VectorSpace, Scalar )+import Data.Basis ( HasBasis, Basis )++import Data.List ( transpose )+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm that automatically computes the+-- 'Jacobian' using automatic differentiation of the model function.+--+-- /Warning/: Don't apply 'levmar' to 'LMA.Model's that apply methods of+-- the 'Eq' and 'Ord' classes to the parameters. These methods are+-- undefined for ':~>'!!!+levmar :: forall m n k r.+ ( Nat m+ , Nat n+ , Nat k+ , HasBasis r+ , Basis r ~ ()+ , VectorSpace (Scalar r)+ , LMA_I.LevMarable r+ )+ => (LMA.Model m n (r :~> r)) -- ^ Model. Note that ':~>'+ -- is overloaded for all the+ -- numeric classes.+ -> SizedList m r -- ^ Initial parameters+ -> SizedList n r -- ^ Samples+ -> Integer -- ^ Maximum number of iterations+ -> LMA_I.Options r -- ^ Minimization options+ -> Maybe (SizedList m r) -- ^ Optional lower bounds+ -> Maybe (SizedList m r) -- ^ Optional upper bounds+ -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints+ -> Maybe (SizedList m r) -- ^ Optional weights+ -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r)++levmar model params ys itMax opts mLowBs mUpBs mLinC mWghts =+ fmap convertResult $ LMA_I.levmar (convertModel model)+ (Just $ jacobianOf model)+ (toList params)+ (toList ys)+ itMax+ opts+ (fmap toList mLowBs)+ (fmap toList mUpBs)+ (fmap convertLinearConstraints mLinC)+ (fmap toList mWghts)+ where+ convertModel :: LMA.Model m n (r :~> r) -> LMA_I.Model r+ (convertModel mdl) ps = fmap value $ toList+ (mdl $* pDs :: SizedList n (r :~> r))+ where+ pDs :: SizedList m (r :~> r)+ pDs = unsafeFromList $ fmap constant ps++ jacobianOf :: LMA.Model m n (r :~> r) -> LMA_I.Jacobian r+ (jacobianOf mdl) ps = fmap (\fs -> zipWith (firstDeriv .) fs ps)+ . transpose+ . fmap (\pD -> toList (mdl $* (pD :: SizedList m (r :~> r)) :: SizedList n (r :~> r)))+ $ pDs+ where+ pDs :: [SizedList m (r :~> r)]+ pDs = [unsafeFromList $ idDAt n ps | n <- [0 .. length ps - 1]]+++-- The End ---------------------------------------------------------------------
LevMar/Fitting.hs view
@@ -21,48 +21,56 @@ module LevMar.Fitting ( -- * Model & Jacobian. Model+ , SimpleModel , Jacobian+ , SimpleJacobian -- * Levenberg-Marquardt algorithm.- , LMA.LevMarable+ , LMA_I.LevMarable , levmar - , LMA.LinearConstraints- , LMA.noLinearConstraints- , LMA.Matrix+ , LinearConstraints+ , noLinearConstraints+ , Matrix -- * Minimization options.- , LMA.Options(..)- , LMA.defaultOpts+ , LMA_I.Options(..)+ , LMA_I.defaultOpts -- * Output- , LMA.Info(..)- , LMA.StopReason(..)- , LMA.CovarMatrix+ , LMA_I.Info(..)+ , LMA_I.StopReason(..)+ , CovarMatrix - , LMA.LevMarError(..)+ , LMA_I.LevMarError(..) -- *Type-level machinery , Z, S, Nat , SizedList(..) , NFunction- , ComposeN ) where -import qualified LevMar as LMA+import qualified LevMar.Intermediate.Fitting as LMA_I+import LevMar.Utils ( LinearConstraints+ , noLinearConstraints+ , convertLinearConstraints+ , Matrix+ , CovarMatrix+ , convertResult+ ) -import TypeLevelNat (Z, S, Nat, witnessNat)-import SizedList (SizedList)-import NFunction (NFunction, ComposeN, compose)+import TypeLevelNat ( Z, S, Nat )+import SizedList ( SizedList(..), toList, unsafeFromList )+import NFunction ( NFunction, ($*) ) -------------------------------------------------------------------------------- -- Model & Jacobian. -------------------------------------------------------------------------------- -{- | A function from @n@ parameters of type @r@ and an x-value of type-@a@ to a value of type @r@.+{- | A functional relation describing measurements represented as a function+from @m@ parameters and an x-value to an expected measurement. For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be written as:@@ -74,12 +82,15 @@ quad a b c x = a*x^2 + b*x + c @ -}-type Model n r a = NFunction n r (a -> r)+type Model m r a = NFunction m r (a -> r) -{- | The jacobian of the 'Model' function. Expressed as a function from-@n@ parameters of type @r@ and an x-value of type @a@ to a vector-of @n@ values of type @r@.+-- | This type synonym expresses that usually the @a@ in @'Model' m r a@+-- equals the type of the parameters.+type SimpleModel m r = Model m r r +{- | The jacobian of the 'Model' function. Expressed as a function from @n@+parameters and an x-value to the @m@ partial derivatives of the parameters.+ See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant> For example, the jacobian of the quadratic function @f(x) = a*x^2 +@@ -97,40 +108,44 @@ Notice you don't have to differentiate for @x@. -}-type Jacobian n r a = NFunction n r (a -> SizedList n r)+type Jacobian m r a = NFunction m r (a -> SizedList m r) +-- | This type synonym expresses that usually the @a@ in @'Jacobian' m r a@+-- equals the type of the parameters.+type SimpleJacobian m r = Jacobian m r r + -------------------------------------------------------------------------------- -- Levenberg-Marquardt algorithm. -------------------------------------------------------------------------------- -- | The Levenberg-Marquardt algorithm specialised for curve-fitting.-levmar :: forall n k r a. (Nat n, ComposeN n, Nat k, LMA.LevMarable r)- => (Model n r a) -- ^ Model- -> Maybe (Jacobian n r a) -- ^ Optional jacobian- -> SizedList n r -- ^ Initial parameters+levmar :: forall m k r a. (Nat m, Nat k, LMA_I.LevMarable r)+ => (Model m r a) -- ^ Model+ -> Maybe (Jacobian m r a) -- ^ Optional jacobian+ -> SizedList m r -- ^ Initial parameters -> [(a, r)] -- ^ Samples -> Integer -- ^ Maximum number of iterations- -> LMA.Options r -- ^ Options- -> Maybe (SizedList n r) -- ^ Optional lower bounds- -> Maybe (SizedList n r) -- ^ Optional upper bounds- -> Maybe (LMA.LinearConstraints k n r) -- ^ Optional linear constraints- -> Maybe (SizedList n r) -- ^ Optional weights- -> Either LMA.LevMarError (SizedList n r, LMA.Info r, LMA.CovarMatrix n r)-levmar model mJac params samples = LMA.levmar (convertModel model)- (fmap convertJacob mJac)- params- ys+ -> LMA_I.Options r -- ^ Minimization options+ -> Maybe (SizedList m r) -- ^ Optional lower bounds+ -> Maybe (SizedList m r) -- ^ Optional upper bounds+ -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints+ -> Maybe (SizedList m r) -- ^ Optional weights+ -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r)+levmar model mJac params ys itMax opts mLowBs mUpBs mLinC mWghts =+ fmap convertResult $ LMA_I.levmar (convertModel model)+ (fmap convertJacob mJac)+ (toList params)+ ys+ itMax+ opts+ (fmap toList mLowBs)+ (fmap toList mUpBs)+ (fmap convertLinearConstraints mLinC)+ (fmap toList mWghts) where- (xs, ys) = unzip samples-- convertModel :: Model n r a -> LMA.Model n r- convertModel = compose (witnessNat :: n) (undefined :: r)- (\(f :: a -> r) -> map f xs)-- convertJacob :: Jacobian n r a -> LMA.Jacobian n r- convertJacob = compose (witnessNat :: n) (undefined :: r)- (\(f :: a -> SizedList n r) -> map f xs)+ convertModel mdl = \ps -> mdl $* (unsafeFromList ps :: SizedList m r)+ convertJacob jac = \ps x -> toList ((jac $* (unsafeFromList ps :: SizedList m r)) x :: SizedList m r) -- The End ---------------------------------------------------------------------
+ LevMar/Fitting/AD.hs view
@@ -0,0 +1,139 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}++--------------------------------------------------------------------------------+-- |+-- Module : LevMar.Fitting.AD+-- Copyright : (c) 2009 Roel van Dijk & Bas van Dijk+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : vandijk.roel@gmail.com, v.dijk.bas@gmail.com+-- Stability : Experimental+--+-- This module provides the Levenberg-Marquardt algorithm specialised+-- for curve-fitting that uses Automatic Differentiation to+-- automatically compute the Jacobian.+--+-- For additional documentation see the documentation of the levmar C+-- library which this library is based on:+-- <http://www.ics.forth.gr/~lourakis/levmar/>+--+--------------------------------------------------------------------------------++module LevMar.Fitting.AD+ ( -- * Model.+ LMA.Model+ , LMA.SimpleModel++ -- * Levenberg-Marquardt algorithm.+ , LMA_I.LevMarable+ , levmar++ , LinearConstraints+ , noLinearConstraints+ , Matrix++ -- * Minimization options.+ , LMA_I.Options(..)+ , LMA_I.defaultOpts++ -- * Output+ , LMA_I.Info(..)+ , LMA_I.StopReason(..)+ , CovarMatrix++ , LMA_I.LevMarError(..)++ -- *Type-level machinery+ , Z, S, Nat+ , SizedList(..)+ , NFunction+ ) where+++import qualified LevMar.Fitting as LMA+import qualified LevMar.Intermediate.Fitting as LMA_I++import LevMar.Utils ( LinearConstraints+ , noLinearConstraints+ , convertLinearConstraints+ , Matrix+ , CovarMatrix+ , convertResult+ )++import TypeLevelNat ( Z, S, Nat )+import SizedList ( SizedList(..), toList, unsafeFromList )+import NFunction ( NFunction, ($*) )++import LevMar.Utils.AD ( value, firstDeriv, constant, idDAt )++-- From vector-space:+import Data.Derivative ( (:~>) )+import Data.VectorSpace ( VectorSpace, Scalar )+import Data.Basis ( HasBasis, Basis )+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm specialised for curve-fitting+-- that automatically computes the 'Jacobian' using automatic+-- differentiation of the model function.+--+-- /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of+-- the 'Eq' and 'Ord' classes to the parameters. These methods are+-- undefined for ':~>'!!!+levmar :: forall m k r a.+ ( Nat m+ , Nat k+ , HasBasis r+ , Basis r ~ ()+ , VectorSpace (Scalar r)+ , LMA_I.LevMarable r+ )+ => LMA.Model m (r :~> r) a -- ^ Model. Note that+ -- ':~>' is overloaded+ -- for all the numeric+ -- classes.+ -> SizedList m r -- ^ Initial parameters+ -> [(a, r)] -- ^ Samples+ -> Integer -- ^ Maximum number of iterations+ -> LMA_I.Options r -- ^ Minimization options+ -> Maybe (SizedList m r) -- ^ Optional lower bounds+ -> Maybe (SizedList m r) -- ^ Optional upper bounds+ -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints+ -> Maybe (SizedList m r) -- ^ Optional weights+ -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r)++levmar model params ys itMax opts mLowBs mUpBs mLinC mWghts =+ fmap convertResult $ LMA_I.levmar (convertModel model)+ (Just $ jacobianOf model)+ (toList params)+ ys+ itMax+ opts+ (fmap toList mLowBs)+ (fmap toList mUpBs)+ (fmap convertLinearConstraints mLinC)+ (fmap toList mWghts)+ where+ convertModel :: LMA.Model m (r :~> r) a -> LMA_I.Model r a+ (convertModel f) ps x = value $ (f $* pDs :: a -> r :~> r) x+ where+ pDs :: SizedList m (r :~> r)+ pDs = unsafeFromList $ fmap constant ps++ jacobianOf :: LMA.Model m (r :~> r) a -> LMA_I.Jacobian r a+ (jacobianOf f) ps x = fmap combine $ zip [0..] ps+ where+ combine (ix, p) = firstDeriv $ (f $* pDs :: a -> r :~> r) x p+ where+ pDs :: SizedList m (r :~> r)+ pDs = unsafeFromList $ idDAt ix ps+++-- The End ---------------------------------------------------------------------
LevMar/Intermediate.hs view
@@ -42,12 +42,12 @@ ) where -import Foreign.Marshal.Array (allocaArray, peekArray, pokeArray, withArray)-import Foreign.Ptr (Ptr, nullPtr, plusPtr)-import Foreign.Storable (Storable)-import Foreign.C.Types (CInt)-import System.IO.Unsafe (unsafePerformIO)-import Data.Maybe (fromJust, fromMaybe, isJust)+import Foreign.Marshal.Array ( allocaArray, peekArray, pokeArray, withArray )+import Foreign.Ptr ( Ptr, nullPtr, plusPtr )+import Foreign.Storable ( Storable )+import Foreign.C.Types ( CInt )+import System.IO.Unsafe ( unsafePerformIO )+import Data.Maybe ( fromJust, fromMaybe, isJust ) import Control.Monad.Instances -- for 'instance Functor (Either a)' import qualified Bindings.LevMar.CurryFriendly as LMA_C@@ -57,9 +57,51 @@ -- Model & Jacobian. -------------------------------------------------------------------------------- +{- | A functional relation describing measurements represented as a function+from a list of parameters to a list of expected measurements.++ * Ensure that the length of the parameters list equals the length of the+ initial parameters list in 'levmar'.++ * Ensure that the length of the ouput list equals the length of the samples+ list in 'levmar'.++For example:++@+hatfldc :: Model Double+hatfldc [p0, p1, p2, p3] = [ p0 - 1.0+ , p0 - sqrt p1+ , p1 - sqrt p2+ , p3 - 1.0+ ]+@+-} type Model r = [r] -> [r] --- | See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>+{- | The jacobian of the 'Model' function. Expressed as a function from a list+of parameters to a list of lists which for each expected measurement describes+the partial derivatives of the parameters.++See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>++ * Ensure that the length of the parameter list equals the length of the initial+ parameter list in 'levmar'.++ * Ensure that the output matrix has the dimension @n@/x/@m@ where @n@ is the+ number of samples and @m@ is the number of parameters.++For example the jacobian of the above @hatfldc@ model is:++@+hatfldc_jac :: Jacobian Double+hatfldc_jac _ p1 p2 _ = [ [1.0, 0.0, 0.0, 0.0]+ , [1.0, -0.5 / sqrt p1, 0.0, 0.0]+ , [0.0, 1.0, -0.5 / sqrt p2, 0.0]+ , [0.0, 0.0, 0.0, 1.0]+ ]+@+-} type Jacobian r = [r] -> [[r]] @@ -144,7 +186,7 @@ f_blec_der f_blec_dif model mJac ps ys itMax opts mLowBs mUpBs mLinC mWeights- = unsafePerformIO $+ = unsafePerformIO . withArray (map realToFrac ps) $ \psPtr -> withArray (map realToFrac ys) $ \ysPtr -> withArray (map realToFrac $ optsToList opts) $ \optsPtr ->@@ -188,7 +230,7 @@ if r < 0 && r /= LMA_C._LM_ERROR_SINGULAR_MATRIX -- we don't treat these two as an error && r /= LMA_C._LM_ERROR_SUM_OF_SQUARES_NOT_FINITE- then return $ Left $ convertLevMarError r+ then return . Left $ convertLevMarError r else do result <- peekArray lenPs psPtr info <- peekArray LMA_C._LM_INFO_SZ infoPtr @@ -223,22 +265,21 @@ withLinConstraints f g = withArray (map realToFrac $ concat cMat) $ \cMatPtr -> withArray (map realToFrac rhcVec) $ \rhcVecPtr ->- f $ g cMatPtr rhcVecPtr $ fromIntegral $ length cMat+ f . g cMatPtr rhcVecPtr . fromIntegral $ length cMat - withWeights f g = maybeWithArray ((fmap . fmap) realToFrac mWeights) $ \weightsPtr ->- f $ g weightsPtr+ withWeights f g = maybeWithArray ((fmap . fmap) realToFrac mWeights) $ f . g convertModel :: (Real r, Fractional r, Storable c, Real c, Fractional c) => Model r -> LMA_C.Model c convertModel model = \parPtr hxPtr numPar _ _ -> do params <- peekArray (fromIntegral numPar) parPtr- pokeArray hxPtr $ map realToFrac $ model $ map realToFrac params+ pokeArray hxPtr . map realToFrac . model $ map realToFrac params convertJacobian :: (Real r, Fractional r, Storable c, Real c, Fractional c) => Jacobian r -> LMA_C.Jacobian c convertJacobian jac = \parPtr jPtr numPar _ _ -> do params <- peekArray (fromIntegral numPar) parPtr- pokeArray jPtr $ concatMap (map realToFrac) $ jac $ map realToFrac params+ pokeArray jPtr . concatMap (map realToFrac) . jac $ map realToFrac params maybeWithArray :: Storable a => Maybe [a] -> (Ptr a -> IO b) -> IO b maybeWithArray Nothing f = f nullPtr@@ -333,7 +374,8 @@ -------------------------------------------------------------------------------- data LevMarError- = LapackError -- ^ A call to a lapack subroutine failed in the underlying C levmar library.+ = LevMarError -- ^ Generic error (not one of the others)+ | LapackError -- ^ A call to a lapack subroutine failed in the underlying C levmar library. | FailedBoxCheck -- ^ At least one lower bound exceeds the upper one. | MemoryAllocationFailure -- ^ A call to @malloc@ failed in the underlying C levmar library. | ConstraintMatrixRowsGtCols -- ^ The matrix of constraints cannot have more rows than columns.@@ -346,7 +388,8 @@ levmarCErrorToLevMarError :: [(CInt, LevMarError)] levmarCErrorToLevMarError =- [ (LMA_C._LM_ERROR_LAPACK_ERROR, LapackError)+ [ (LMA_C._LM_ERROR, LevMarError)+ , (LMA_C._LM_ERROR_LAPACK_ERROR, LapackError) --, (LMA_C._LM_ERROR_NO_JACOBIAN, can never happen) --, (LMA_C._LM_ERROR_NO_BOX_CONSTRAINTS, can never happen) , (LMA_C._LM_ERROR_FAILED_BOX_CHECK, FailedBoxCheck)
+ LevMar/Intermediate/AD.hs view
@@ -0,0 +1,101 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}++--------------------------------------------------------------------------------+-- |+-- Module : LevMar.Intermediate.AD+-- Copyright : (c) 2009 Roel van Dijk & Bas van Dijk+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : vandijk.roel@gmail.com, v.dijk.bas@gmail.com+-- Stability : Experimental+--+-- A levmar variant that uses Automatic Differentiation to+-- automatically compute the Jacobian.+--+-- For additional documentation see the documentation of the levmar C+-- library which this library is based on:+-- <http://www.ics.forth.gr/~lourakis/levmar/>+--+--------------------------------------------------------------------------------++module LevMar.Intermediate.AD+ ( -- * Model.+ LMA_I.Model++ -- * Levenberg-Marquardt algorithm.+ , LMA_I.LevMarable+ , levmar++ , LMA_I.LinearConstraints++ -- * Minimization options.+ , LMA_I.Options(..)+ , LMA_I.defaultOpts++ -- * Output+ , LMA_I.Info(..)+ , LMA_I.StopReason(..)+ , LMA_I.CovarMatrix++ , LMA_I.LevMarError(..)+ ) where+++import qualified LevMar.Intermediate as LMA_I++import LevMar.Utils.AD ( value, firstDeriv, constant, idDAt )++-- From vector-space:+import Data.Derivative ( (:~>) )+import Data.VectorSpace ( VectorSpace, Scalar )+import Data.Basis ( HasBasis, Basis )++import Data.List ( transpose )+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm that automatically computes the+-- 'Jacobian' using automatic differentiation of the model function.+--+-- /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of+-- the 'Eq' and 'Ord' classes to the parameters. These methods are+-- undefined for ':~>'!!!+levmar :: forall r.+ ( HasBasis r+ , Basis r ~ ()+ , VectorSpace (Scalar r)+ , LMA_I.LevMarable r+ )+ => LMA_I.Model (r :~> r) -- ^ Model. Note that+ -- ':~>' is overloaded+ -- for all the numeric+ -- classes.+ -> [r] -- ^ Initial parameters+ -> [r] -- ^ Samples+ -> Integer -- ^ Maximum iterations+ -> LMA_I.Options r -- ^ Minimization options+ -> Maybe [r] -- ^ Optional lower bounds+ -> Maybe [r] -- ^ Optional upper bounds+ -> Maybe (LMA_I.LinearConstraints r) -- ^ Optional linear constraints+ -> Maybe [r] -- ^ Optional weights+ -> Either LMA_I.LevMarError ([r], LMA_I.Info r, LMA_I.CovarMatrix r)++levmar model = LMA_I.levmar (convertModel model) . Just $ jacobianOf model+ where+ convertModel :: LMA_I.Model (r :~> r) -> LMA_I.Model r+ convertModel mdl = map value . mdl . map constant++ jacobianOf :: LMA_I.Model (r :~> r) -> LMA_I.Jacobian r+ (jacobianOf mdl) ps = map (\fs -> zipWith (firstDeriv .) fs ps)+ . transpose $ map mdl pDs+ where+ pDs = [idDAt n ps | n <- [0 .. length ps - 1]]+++-- The End ---------------------------------------------------------------------
LevMar/Intermediate/Fitting.hs view
@@ -19,7 +19,9 @@ module LevMar.Intermediate.Fitting ( -- * Model & Jacobian. Model+ , SimpleModel , Jacobian+ , SimpleJacobian -- * Levenberg-Marquardt algorithm. , LMA_I.LevMarable@@ -47,12 +49,56 @@ -- Model & Jacobian. -------------------------------------------------------------------------------- +{- | A functional relation describing measurements represented as a function+from a list of parameters and an x-value to an expected measurement.++ * Ensure that the length of the parameters list equals the lenght of the initial+ parameters list in 'levmar'.++For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be+written as:++@+quad :: 'Num' r => 'Model' r r+quad [a, b, c] x = a*x^2 + b*x + c+@+-} type Model r a = [r] -> a -> r --- | See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>+-- | This type synonym expresses that usually the @a@ in @'Model' r a@+-- equals the type of the parameters.+type SimpleModel r = Model r r++{- | The jacobian of the 'Model' function. Expressed as a function from a list+of parameters and an x-value to the partial derivatives of the parameters.++See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>++ * Ensure that the length of the parameters list equals the lenght of the initial+ parameters list in 'levmar'.++ * Ensure that the length of the output parameter derivatives list equals the+ length of the input parameters list.++For example, the jacobian of the above @quad@ model can be written as:++@+quadJacob :: 'Num' r => 'Jacobian' N3 r r+quadJacob [_, _, _] x = [ x^2 -- with respect to a+ , x -- with respect to b+ , 1 -- with respect to c+ ]+@++Notice you don't have to differentiate for @x@.+-} type Jacobian r a = [r] -> a -> [r] +-- | This type synonym expresses that usually the @a@ in @'Jacobian' r a@+-- equals the type of the parameters.+type SimpleJacobian r = Jacobian r r + -------------------------------------------------------------------------------- -- Levenberg-Marquardt algorithm. --------------------------------------------------------------------------------@@ -70,13 +116,16 @@ -> Maybe (LMA_I.LinearConstraints r) -- ^ Optional linear constraints -> Maybe [r] -- ^ Optional weights -> Either LMA_I.LevMarError ([r], LMA_I.Info r, LMA_I.CovarMatrix r)-levmar model mJac ps samples =- LMA_I.levmar (\ps' -> map (model ps') xs)- (fmap (\jac -> \ps' -> map (jac ps') xs) mJac)- ps+levmar model mJac params samples =+ LMA_I.levmar (convertModel model)+ (fmap convertJacob mJac)+ params ys where (xs, ys) = unzip samples++ convertModel mdl = \ps -> map (mdl ps) xs+ convertJacob jac = \ps -> map (jac ps) xs -- The End ---------------------------------------------------------------------
+ LevMar/Intermediate/Fitting/AD.hs view
@@ -0,0 +1,100 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}++--------------------------------------------------------------------------------+-- |+-- Module : LevMar.Intermediate.Fitting.AD+-- Copyright : (c) 2009 Roel van Dijk & Bas van Dijk+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : vandijk.roel@gmail.com, v.dijk.bas@gmail.com+-- Stability : Experimental+--+-- A levmar variant specialised for curve-fitting that uses Automatic+-- Differentiation to automatically compute the Jacobian.+--+-- For additional documentation see the documentation of the levmar C+-- library which this library is based on:+-- <http://www.ics.forth.gr/~lourakis/levmar/>+--+--------------------------------------------------------------------------------++module LevMar.Intermediate.Fitting.AD+ ( -- * Model.+ LMA_I.Model+ , LMA_I.SimpleModel++ -- * Levenberg-Marquardt algorithm.+ , LMA_I.LevMarable+ , levmar++ , LMA_I.LinearConstraints++ -- * Minimization options.+ , LMA_I.Options(..)+ , LMA_I.defaultOpts++ -- * Output+ , LMA_I.Info(..)+ , LMA_I.StopReason(..)+ , LMA_I.CovarMatrix++ , LMA_I.LevMarError(..)+ ) where+++import qualified LevMar.Intermediate.Fitting as LMA_I++import LevMar.Utils.AD ( value, firstDeriv, constant, idDAt )++-- From vector-space:+import Data.Derivative ( (:~>) )+import Data.VectorSpace ( VectorSpace, Scalar )+import Data.Basis ( HasBasis, Basis )+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm specialised for curve-fitting+-- that automatically computes the 'Jacobian' using automatic+-- differentiation of the model function.+--+-- /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of+-- the 'Eq' and 'Ord' classes to the parameters. These methods are+-- undefined for ':~>'!!!+levmar :: forall r a.+ ( HasBasis r+ , Basis r ~ ()+ , VectorSpace (Scalar r)+ , LMA_I.LevMarable r+ )+ => LMA_I.Model (r :~> r) a -- ^ Model. Note that+ -- ':~>' is overloaded+ -- for all the numeric+ -- classes.+ -> [r] -- ^ Initial parameters+ -> [(a, r)] -- ^ Samples+ -> Integer -- ^ Maximum iterations+ -> LMA_I.Options r -- ^ Minimization options+ -> Maybe [r] -- ^ Optional lower bounds+ -> Maybe [r] -- ^ Optional upper bounds+ -> Maybe (LMA_I.LinearConstraints r) -- ^ Optional linear constraints+ -> Maybe [r] -- ^ Optional weights+ -> Either LMA_I.LevMarError ([r], LMA_I.Info r, LMA_I.CovarMatrix r)++levmar model = LMA_I.levmar (convertModel model) . Just $ jacobianOf model+ where+ convertModel :: LMA_I.Model (r :~> r) a -> LMA_I.Model r a+ convertModel mdl = \ps -> value . mdl (map constant ps)++ jacobianOf :: LMA_I.Model (r :~> r) a -> LMA_I.Jacobian r a+ jacobianOf mdl =+ \ps x -> map (\(ix, p) -> firstDeriv $ mdl (idDAt ix ps) x p)+ $ zip [0..] ps+++-- The End ---------------------------------------------------------------------
+ LevMar/Utils.hs view
@@ -0,0 +1,44 @@+module LevMar.Utils+ ( LinearConstraints+ , noLinearConstraints+ , Matrix+ , CovarMatrix+ , convertLinearConstraints+ , convertResult+ ) where++import qualified LevMar.Intermediate as LMA_I++import TypeLevelNat ( Nat, Z )+import SizedList ( SizedList, toList, unsafeFromList )++-- | Linear constraints consisting of a constraints matrix, /kxn/ and+-- a right hand constraints vector, /kx1/ where /n/ is the number of+-- parameters and /k/ is the number of constraints.+type LinearConstraints k n r = (Matrix k n r, SizedList k r)++-- |Value to denote the absense of any linear constraints over the+-- parameters of the model function. Use this instead of 'Nothing'+-- because the type parameter which contains the number of constraints+-- can't be inferred.+noLinearConstraints :: Nat n => Maybe (LinearConstraints Z n r)+noLinearConstraints = Nothing++-- | A /nxm/ matrix is a sized list of /n/ sized lists of length /m/.+type Matrix n m r = SizedList n (SizedList m r)++-- | Covariance matrix corresponding to LS solution.+type CovarMatrix n r = Matrix n n r++convertLinearConstraints :: (Nat k, Nat n) => LinearConstraints k n r -> LMA_I.LinearConstraints r+convertLinearConstraints (cMat, rhcVec) = ( map toList $ toList cMat+ , toList rhcVec+ )++convertResult :: (Nat n)+ => ([r], LMA_I.Info r, LMA_I.CovarMatrix r)+ -> (SizedList n r, LMA_I.Info r, CovarMatrix n r)+convertResult (psResult, info, covar) = ( unsafeFromList psResult+ , info+ , unsafeFromList $ map unsafeFromList covar+ )
+ LevMar/Utils/AD.hs view
@@ -0,0 +1,42 @@+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}++module LevMar.Utils.AD where++import Data.Derivative ( (:~>), (:>), powVal, idD, pureD, derivAtBasis )+import Data.VectorSpace ( VectorSpace, Scalar, AdditiveGroup )+import Data.Basis ( HasBasis, Basis )+import Data.MemoTrie ( HasTrie )+++value :: a :~> b -> b+value m = powVal $ m undefined++-- | @firstDeriv f@ returns the first derivative of @f@.+firstDeriv :: (HasBasis a, Basis a ~ (), AdditiveGroup b)+ => (a :> b) -> b+firstDeriv f = powVal $ derivAtBasis f ()++-- | A constant infinitely differentiable function.+constant :: (AdditiveGroup b, HasBasis a, HasTrie (Basis a))+ => b -> a:~>b+constant = const . pureD++-- | @idDAt n ps@ maps each parameter in @ps@ to a /constant/+-- infinitely differentiable function (@const . pureD@), except the @n@th+-- parameter is replaced with the differentiable /identity/ function+-- (@idD@).+idDAt :: (HasBasis r, HasTrie (Basis r), VectorSpace (Scalar r))+ => Int -> [r] -> [r :~> r]+idDAt n = replace n idD . map constant++-- | @replace i r xs@ replaces the @i@th element in @xs@ with @r@.+replace :: Int -> a -> [a] -> [a]+replace i r xs+ | i < 0 = xs+ | otherwise = rep i xs+ where rep _ [] = []+ rep j (y:ys)+ | j > 0 = y : rep (j - 1) ys+ | otherwise = r : ys
NFunction.hs view
@@ -8,8 +8,8 @@ , compose ) where -import TypeLevelNat (Z(..), S(..), Nat)-import SizedList (SizedList(..))+import TypeLevelNat ( Z(..), S(..), Nat )+import SizedList ( SizedList(..) ) -- | A @NFunction n a b@ is a function which takes @n@ arguments of -- type @a@ and returns a @b@.
SizedList.hs view
@@ -1,27 +1,41 @@ {-# LANGUAGE GADTs #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-} module SizedList ( SizedList(..)+ , foldr+ , foldrN , toList+ , length , fromList , unsafeFromList- , length , replicate ) where -import Prelude hiding (replicate, length)-import Data.Maybe (fromMaybe)-import TypeLevelNat (Z(..), S(..), Nat, induction, witnessNat, N(..)) +import Prelude hiding ( foldr, replicate, length )+import Data.Maybe ( fromMaybe )+import TypeLevelNat ( Z(..), S(..), Nat, induction, witnessNat, N(..) )+++--------------------------------------------------------------------------------+ -- | A list which is indexed with a type-level natural that denotes the size of -- the list. data SizedList n a where Nil :: SizedList Z a (:::) :: a -> SizedList n a -> SizedList (S n) a +instance Functor (SizedList n) where+ fmap _ Nil = Nil+ fmap f (x ::: xs) = f x ::: fmap f xs+ infixr 5 ::: -- Same precedence and associativity as (:) ++--------------------------------------------------------------------------------+ consPrecedence :: Int consPrecedence = 5 @@ -32,31 +46,46 @@ . showString " ::: " . showsPrec consPrecedence xs -newtype ToList a n = ToList { unToList :: SizedList n a -> [a] } --- | Convert a @SizedList@ to a normal list.-toList :: forall a n. Nat n => SizedList n a -> [a]-toList = unToList $ induction (witnessNat :: n)- (ToList tl0)- (ToList . tlS . unToList)+--------------------------------------------------------------------------------++-- | Fold a binary operator over a @SizedList@.+foldr :: forall a b n. (a -> b -> b) -> b -> SizedList n a -> b+foldr f z = foldr_f_z where- tl0 :: SizedList Z a -> [a]- tl0 Nil = []+ foldr_f_z :: forall k. SizedList k a -> b+ foldr_f_z Nil = z+ foldr_f_z (x ::: xs) = f x $ foldr_f_z xs - tlS :: forall x. Nat x => (SizedList x a -> [a]) -> SizedList (S x) a -> [a]- tlS f (x ::: xs) = x : f xs+-- | Fold a binary operator yielding a value with a natural number+-- indexed type over a @SizedList@.+foldrN :: forall a b n. (forall m. a -> b m -> b (S m)) -> b Z -> SizedList n a -> b n+foldrN f z = foldrN_f_z+ where+ foldrN_f_z :: forall k. SizedList k a -> b k+ foldrN_f_z Nil = z+ foldrN_f_z (x ::: xs) = f x $ foldrN_f_z xs -newtype FromList a n = FromList { unFromList :: [a] -> Maybe (SizedList n a) }+-- | Convert a @SizedList@ to a normal list.+toList :: SizedList n a -> [a]+toList = foldr (:) [] --- | Convert a normal list to a @SizeList@. If the length of the given+-- | Returns the length of the @SizedList@.+length :: SizedList n a -> N n+length = foldrN (const Succ) Zero+++--------------------------------------------------------------------------------++newtype FromList a n = FL { unFL :: [a] -> Maybe (SizedList n a) }++-- | Convert a normal list to a @SizedList@. If the length of the given -- list does not equal @n@, @Nothing@ is returned. fromList :: forall a n. Nat n => [a] -> Maybe (SizedList n a)-fromList = unFromList $ induction (witnessNat :: n)- (FromList fl0)- (FromList . flS . unFromList)+fromList = unFL $ induction (witnessNat :: n) (FL flZ) (FL . flS . unFL) where- fl0 [] = Just Nil- fl0 _ = Nothing+ flZ [] = Just Nil+ flZ _ = Nothing flS _ [] = Nothing flS k (x:xs) = fmap (x :::) $ k xs@@ -67,10 +96,14 @@ unsafeFromList = fromMaybe (error "unsafeFromList xs: xs does not have the right length ") . fromList -replicate :: N n -> a -> SizedList n a-replicate Zero _ = Nil-replicate (Succ n) x = x ::: replicate n x -length :: SizedList n a -> N n-length Nil = Zero-length (_ ::: xs) = Succ $ length xs+--------------------------------------------------------------------------------++newtype Replicate a n = R { unR :: SizedList n a}++-- | @replicate x :: SizedList n a@ returns a @SizedList@ of @n@ @x@s.+replicate :: forall a n. Nat n => a -> SizedList n a+replicate x = unR $ induction (witnessNat :: n) (R Nil) (R . (x :::) . unR)+++-- The End ---------------------------------------------------------------------
TypeLevelNat.hs view
@@ -16,6 +16,7 @@ , witnessNat , N(..)+ , nat ) where @@ -76,6 +77,9 @@ data N n where Zero :: N Z Succ :: N n -> N (S n)++nat :: forall n. Nat n => n -> N n+nat n = induction n Zero Succ {- Template Haskell code to construct a type synonym for an arbitrary
levmar.cabal view
@@ -1,14 +1,15 @@ name: levmar-version: 0.1+version: 0.2 cabal-version: >= 1.6 build-type: Simple stability: experimental+tested-with: GHC ==6.10.4 author: Roel van Dijk & Bas van Dijk maintainer: vandijk.roel@gmail.com, v.dijk.bas@gmail.com copyright: (c) 2009 Roel van Dijk & Bas van Dijk license: BSD3 license-file: LICENSE-category: numerical+category: Numerical, Math synopsis: An implementation of the Levenberg-Marquardt algorithm description: The Levenberg-Marquardt algorithm is an iterative technique that finds a local minimum of a function that@@ -40,12 +41,26 @@ * LevMar: A high-level layer that uses type-level programming to add extra type safety. .- Each layer also has special data-fitting variants:+ Each layer also has special curve-fitting variants: . * LevMar.Intermediate.Fitting . * LevMar.Fitting .+ Each layer also has special variants that automatically compute+ the jacobian using automatic differentiation using Conal+ Elliott's vector-space library:+ .+ * LevMar.Intermediate.AD+ .+ * LevMar.Intermediate.Fitting.AD+ .+ * LevMar.AD+ .+ * LevMar.Fitting.AD+ .+ Note however that this feature is still very experimental!+ . All modules are self-contained; i.e. each module re-exports all the things you need to work with it. .@@ -68,12 +83,20 @@ library build-depends: base >= 3 && < 4.2- , bindings-levmar < 0.2+ , bindings-levmar == 0.1.*+ , vector-space >= 0.5.7 && < 0.6+ , MemoTrie >= 0.4.5 && < 0.5 exposed-modules: LevMar+ , LevMar.AD , LevMar.Fitting+ , LevMar.Fitting.AD , LevMar.Intermediate+ , LevMar.Intermediate.AD , LevMar.Intermediate.Fitting+ , LevMar.Intermediate.Fitting.AD , TypeLevelNat , SizedList , NFunction+ other-modules: LevMar.Utils+ , LevMar.Utils.AD ghc-options: -Wall -O2