levmar (empty) → 0.1
raw patch · 11 files changed
+2085/−0 lines, 11 filesdep +basedep +bindings-levmarsetup-changed
Dependencies added: base, bindings-levmar
Files
- Demo.hs +994/−0
- LICENSE +32/−0
- LevMar.hs +164/−0
- LevMar/Fitting.hs +136/−0
- LevMar/Intermediate.hs +366/−0
- LevMar/Intermediate/Fitting.hs +82/−0
- NFunction.hs +54/−0
- Setup.hs +3/−0
- SizedList.hs +76/−0
- TypeLevelNat.hs +99/−0
- levmar.cabal +79/−0
+ Demo.hs view
@@ -0,0 +1,994 @@+-- 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++-- The End ---------------------------------------------------------------------
+ LICENSE view
@@ -0,0 +1,32 @@+Copyright (c) 2009 Roel van Dijk, Bas van Dijk++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are+met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * The name of Roel van Dijk and Bas van Dijk and the names of+ contributors may NOT be used to endorse or promote products+ derived from this software without specific prior written+ permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ LevMar.hs view
@@ -0,0 +1,164 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}++--------------------------------------------------------------------------------+-- |+-- Module : LevMar+-- 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+--+--+--+-- 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+ ( -- * Model & Jacobian.+ Model+ , Jacobian++ -- * 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.Intermediate as LMA_I++import TypeLevelNat (Z, S, Nat)+import SizedList (SizedList(..), toList, unsafeFromList)+import NFunction (NFunction, ($*))++import Data.Either+++--------------------------------------------------------------------------------+-- Model & Jacobian.+--------------------------------------------------------------------------------++{- | A function from @n@ parameters of type @r@ to a list of @r@.++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+ ]+@+-}+type Model n r = NFunction n r [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@++See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>++For example the jacobian of the above @hatfldc@ model is:++@+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+ ]+@+-}++type Jacobian n r = NFunction n r [SizedList n r]+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | 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+ -> 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)++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 convertLinC 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+++-- The End ---------------------------------------------------------------------
+ LevMar/Fitting.hs view
@@ -0,0 +1,136 @@+{-# LANGUAGE ScopedTypeVariables #-}++--------------------------------------------------------------------------------+-- |+-- Module : LevMar.Fitting+-- 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.+--+-- 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+ ( -- * Model & Jacobian.+ Model+ , Jacobian++ -- * Levenberg-Marquardt algorithm.+ , LMA.LevMarable+ , levmar++ , LMA.LinearConstraints+ , LMA.noLinearConstraints+ , LMA.Matrix++ -- * Minimization options.+ , LMA.Options(..)+ , LMA.defaultOpts++ -- * Output+ , LMA.Info(..)+ , LMA.StopReason(..)+ , LMA.CovarMatrix++ , LMA.LevMarError(..)++ -- *Type-level machinery+ , Z, S, Nat+ , SizedList(..)+ , NFunction+ , ComposeN+ ) where+++import qualified LevMar as LMA++import TypeLevelNat (Z, S, Nat, witnessNat)+import SizedList (SizedList)+import NFunction (NFunction, ComposeN, compose)+++--------------------------------------------------------------------------------+-- Model & Jacobian.+--------------------------------------------------------------------------------++{- | A function from @n@ parameters of type @r@ and an x-value of type+@a@ to a value of type @r@.++For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be+written as:++@+type N3 = 'S' ('S' ('S' 'Z'))++quad :: 'Num' r => 'Model' N3 r r+quad a b c x = a*x^2 + b*x + c+@+-}+type Model n r a = NFunction n 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@.++See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>++For example, the jacobian of the quadratic function @f(x) = a*x^2 ++b*x + c@ can be written as:++@+type N3 = 'S' ('S' ('S' 'Z'))++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+ ::: 'Nil'+@++Notice you don't have to differentiate for @x@.+-}+type Jacobian n r a = NFunction n r (a -> SizedList n 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+ -> [(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+ 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)+++-- The End ---------------------------------------------------------------------
+ LevMar/Intermediate.hs view
@@ -0,0 +1,366 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleInstances #-}++--------------------------------------------------------------------------------+-- |+-- Module : LevMar.Intermediate+-- 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+--+--+--+-- 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+ ( -- * Model & Jacobian.+ Model+ , Jacobian++ -- * Levenberg-Marquardt algorithm.+ , LevMarable+ , levmar++ , LinearConstraints++ -- * Minimization options.+ , Options(..)+ , defaultOpts++ -- * Output+ , Info(..)+ , StopReason(..)+ , CovarMatrix++ , LevMarError(..)+ ) 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 Control.Monad.Instances -- for 'instance Functor (Either a)'++import qualified Bindings.LevMar.CurryFriendly as LMA_C+++--------------------------------------------------------------------------------+-- Model & Jacobian.+--------------------------------------------------------------------------------++type Model r = [r] -> [r]++-- | See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>+type Jacobian r = [r] -> [[r]]+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm is overloaded to work on 'Double' and 'Float'.+class LevMarable r where++ -- | The Levenberg-Marquardt algorithm.+ levmar :: Model r -- ^ Model+ -> Maybe (Jacobian r) -- ^ Optional jacobian+ -> [r] -- ^ Initial parameters+ -> [r] -- ^ Samples+ -> Integer -- ^ Maximum iterations+ -> Options r -- ^ Minimization options+ -> Maybe [r] -- ^ Optional lower bounds+ -> Maybe [r] -- ^ Optional upper bounds+ -> Maybe (LinearConstraints r) -- ^ Optional linear constraints+ -> Maybe [r] -- ^ Optional weights+ -> Either LevMarError ([r], Info r, CovarMatrix r)++instance LevMarable Float where+ levmar = gen_levmar LMA_C.slevmar_der+ LMA_C.slevmar_dif+ LMA_C.slevmar_bc_der+ LMA_C.slevmar_bc_dif+ LMA_C.slevmar_lec_der+ LMA_C.slevmar_lec_dif+ LMA_C.slevmar_blec_der+ LMA_C.slevmar_blec_dif++instance LevMarable Double where+ levmar = gen_levmar LMA_C.dlevmar_der+ LMA_C.dlevmar_dif+ LMA_C.dlevmar_bc_der+ LMA_C.dlevmar_bc_dif+ LMA_C.dlevmar_lec_der+ LMA_C.dlevmar_lec_dif+ LMA_C.dlevmar_blec_der+ LMA_C.dlevmar_blec_dif++{- | @gen_levmar@ takes the low-level C functions as arguments and+executes one of them depending on the optional jacobian and constraints.++Preconditions:+ length ys >= length ps++ isJust mLowBs && length (fromJust mLowBs) == length ps+ && isJust mUpBs && length (fromJust mUpBs) == length ps++ boxConstrained && (all $ zipWith (<=) (fromJust mLowBs) (fromJust mUpBs))+-}+gen_levmar :: forall cr r. (Storable cr, RealFrac cr, Real r, Fractional r)+ => LMA_C.LevMarDer cr+ -> LMA_C.LevMarDif cr+ -> LMA_C.LevMarBCDer cr+ -> LMA_C.LevMarBCDif cr+ -> LMA_C.LevMarLecDer cr+ -> LMA_C.LevMarLecDif cr+ -> LMA_C.LevMarBLecDer cr+ -> LMA_C.LevMarBLecDif cr++ -> Model r -- ^ Model+ -> Maybe (Jacobian r) -- ^ Optional jacobian+ -> [r] -- ^ Initial parameters+ -> [r] -- ^ Samples+ -> Integer -- ^ Maximum iterations+ -> Options r -- ^ Options+ -> Maybe [r] -- ^ Optional lower bounds+ -> Maybe [r] -- ^ Optional upper bounds+ -> Maybe (LinearConstraints r) -- ^ Optional linear constraints+ -> Maybe [r] -- ^ Optional weights+ -> Either LevMarError ([r], Info r, CovarMatrix r)+gen_levmar f_der+ f_dif+ f_bc_der+ f_bc_dif+ f_lec_der+ f_lec_dif+ f_blec_der+ f_blec_dif+ model mJac ps ys itMax opts mLowBs mUpBs mLinC mWeights+ = unsafePerformIO $+ withArray (map realToFrac ps) $ \psPtr ->+ withArray (map realToFrac ys) $ \ysPtr ->+ withArray (map realToFrac $ optsToList opts) $ \optsPtr ->+ allocaArray LMA_C._LM_INFO_SZ $ \infoPtr ->+ allocaArray covarLen $ \covarPtr ->+ LMA_C.withModel (convertModel model) $ \modelPtr -> do++ let runDif :: LMA_C.LevMarDif cr -> IO CInt+ runDif f = f modelPtr+ psPtr+ ysPtr+ (fromIntegral lenPs)+ (fromIntegral lenYs)+ (fromIntegral itMax)+ optsPtr+ infoPtr+ nullPtr+ covarPtr+ nullPtr++ r <- case mJac of+ Just jac -> LMA_C.withJacobian (convertJacobian jac) $ \jacobPtr ->+ let runDer :: LMA_C.LevMarDer cr -> IO CInt+ runDer f = runDif $ f jacobPtr+ in if boxConstrained+ then if linConstrained+ then withBoxConstraints (withLinConstraints $ withWeights runDer) f_blec_der+ else withBoxConstraints runDer f_bc_der+ else if linConstrained+ then withLinConstraints runDer f_lec_der+ else runDer f_der++ Nothing -> if boxConstrained+ then if linConstrained+ then withBoxConstraints (withLinConstraints $ withWeights runDif) f_blec_dif+ else withBoxConstraints runDif f_bc_dif+ else if linConstrained+ then withLinConstraints runDif f_lec_dif+ else runDif f_dif++ 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+ else do result <- peekArray lenPs psPtr+ info <- peekArray LMA_C._LM_INFO_SZ infoPtr++ let covarPtrEnd = plusPtr covarPtr covarLen+ let convertCovarMatrix ptr+ | ptr == covarPtrEnd = return []+ | otherwise = do row <- peekArray lenPs ptr+ rows <- convertCovarMatrix $ plusPtr ptr lenPs+ return $ row : rows++ covar <- convertCovarMatrix covarPtr++ return $ Right ( map realToFrac result+ , listToInfo info+ , map (map realToFrac) covar+ )+ where+ lenPs = length ps+ lenYs = length ys+ covarLen = lenPs * lenPs+ (cMat, rhcVec) = fromJust mLinC++ -- Whether the parameters are constrained by a linear equation.+ linConstrained = isJust mLinC++ -- Whether the parameters are constrained by a bounding box.+ boxConstrained = isJust mLowBs || isJust mUpBs++ withBoxConstraints f g = maybeWithArray ((fmap . fmap) realToFrac mLowBs) $ \lBsPtr ->+ maybeWithArray ((fmap . fmap) realToFrac mUpBs) $ \uBsPtr ->+ f $ g lBsPtr uBsPtr++ withLinConstraints f g = withArray (map realToFrac $ concat cMat) $ \cMatPtr ->+ withArray (map realToFrac rhcVec) $ \rhcVecPtr ->+ f $ g cMatPtr rhcVecPtr $ fromIntegral $ length cMat++ withWeights f g = maybeWithArray ((fmap . fmap) realToFrac mWeights) $ \weightsPtr ->+ f $ g weightsPtr++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++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++maybeWithArray :: Storable a => Maybe [a] -> (Ptr a -> IO b) -> IO b+maybeWithArray Nothing f = f nullPtr+maybeWithArray (Just xs) f = withArray xs f+++-- | Linear constraints consisting of a constraints matrix, /kxm/ and+-- a right hand constraints vector, /kx1/ where /m/ is the number of+-- parameters and /k/ is the number of constraints.+type LinearConstraints r = ([[r]], [r])+++--------------------------------------------------------------------------------+-- Minimization options.+--------------------------------------------------------------------------------++-- | Minimization options+data Options r =+ Opts { optScaleInitMu :: r -- ^ Scale factor for initial /mu/.+ , optStopNormInfJacTe :: r -- ^ Stopping thresholds for @||J^T e||_inf@.+ , optStopNorm2Dp :: r -- ^ Stopping thresholds for @||Dp||_2@.+ , optStopNorm2E :: r -- ^ Stopping thresholds for @||e||_2@.+ , optDelta :: r -- ^ Step used in the difference approximation to the Jacobian.+ -- If @optDelta<0@, the Jacobian is approximated+ -- with central differences which are more accurate+ -- (but slower!) compared to the forward differences+ -- employed by default.+ } deriving Show++-- | Default minimization options+defaultOpts :: Fractional r => Options r+defaultOpts = Opts { optScaleInitMu = LMA_C._LM_INIT_MU+ , optStopNormInfJacTe = LMA_C._LM_STOP_THRESH+ , optStopNorm2Dp = LMA_C._LM_STOP_THRESH+ , optStopNorm2E = LMA_C._LM_STOP_THRESH+ , optDelta = LMA_C._LM_DIFF_DELTA+ }++optsToList :: Options r -> [r]+optsToList (Opts mu eps1 eps2 eps3 delta) =+ [mu, eps1, eps2, eps3, delta]+++--------------------------------------------------------------------------------+-- Output+--------------------------------------------------------------------------------++-- | Information regarding the minimization.+data Info r = Info { infNorm2initE :: r -- ^ @||e||_2@ at initial parameters.+ , infNorm2E :: r -- ^ @||e||_2@ at estimated parameters.+ , infNormInfJacTe :: r -- ^ @||J^T e||_inf@ at estimated parameters.+ , infNorm2Dp :: r -- ^ @||Dp||_2@ at estimated parameters.+ , infMuDivMax :: r -- ^ @\mu/max[J^T J]_ii ]@ at estimated parameters.+ , infNumIter :: Integer -- ^ Number of iterations.+ , infStopReason :: StopReason -- ^ Reason for terminating.+ , infNumFuncEvals :: Integer -- ^ Number of function evaluations.+ , infNumJacobEvals :: Integer -- ^ Number of jacobian evaluations.+ , infNumLinSysSolved :: Integer -- ^ Number of linear systems solved, i.e. attempts for reducing error.+ } deriving Show++listToInfo :: (RealFrac cr, Fractional r) => [cr] -> Info r+listToInfo [a,b,c,d,e,f,g,h,i,j] =+ Info { infNorm2initE = realToFrac a+ , infNorm2E = realToFrac b+ , infNormInfJacTe = realToFrac c+ , infNorm2Dp = realToFrac d+ , infMuDivMax = realToFrac e+ , infNumIter = floor f+ , infStopReason = toEnum $ floor g - 1+ , infNumFuncEvals = floor h+ , infNumJacobEvals = floor i+ , infNumLinSysSolved = floor j+ }+listToInfo _ = error "liftToInfo: wrong list length"++-- | Reason for terminating.+data StopReason = SmallGradient -- ^ Stopped because of small gradient @J^T e@.+ | SmallDp -- ^ Stopped because of small Dp.+ | MaxIterations -- ^ Stopped because maximum iterations was reached.+ | SingularMatrix -- ^ Stopped because of singular matrix. Restart from current estimated parameters with increased 'optScaleInitMu'.+ | SmallestError -- ^ Stopped because no further error reduction is possible. Restart with increased 'optScaleInitMu'.+ | SmallNorm2E -- ^ Stopped because of small @||e||_2@.+ | InvalidValues -- ^ Stopped because model function returned invalid values (i.e. NaN or Inf). This is a user error.+ deriving (Show, Enum)++-- | Covariance matrix corresponding to LS solution.+type CovarMatrix r = [[r]]+++--------------------------------------------------------------------------------+-- Error+--------------------------------------------------------------------------------++data LevMarError+ = 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.+ | ConstraintMatrixNotFullRowRank -- ^ Constraints matrix is not of full row rank.+ | TooFewMeasurements -- ^ Cannot solve a problem with fewer measurements than unknowns.+ -- In case linear constraints are provided, this error is also returned+ -- when the number of measurements is smaller than the number of unknowns+ -- minus the number of equality constraints.+ deriving Show++levmarCErrorToLevMarError :: [(CInt, LevMarError)]+levmarCErrorToLevMarError =+ [ (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)+ , (LMA_C._LM_ERROR_MEMORY_ALLOCATION_FAILURE, MemoryAllocationFailure)+ , (LMA_C._LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS, ConstraintMatrixRowsGtCols)+ , (LMA_C._LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK, ConstraintMatrixNotFullRowRank)+ , (LMA_C._LM_ERROR_TOO_FEW_MEASUREMENTS, TooFewMeasurements)+ --, (LMA_C._LM_ERROR_SINGULAR_MATRIX, we don't treat this as an error)+ --, (LMA_C._LM_ERROR_SUM_OF_SQUARES_NOT_FINITE, we don't treat this as an error)+ ]++convertLevMarError :: CInt -> LevMarError+convertLevMarError err = fromMaybe (error "Unknown levmar error") $+ lookup err levmarCErrorToLevMarError+++-- The End ---------------------------------------------------------------------
+ LevMar/Intermediate/Fitting.hs view
@@ -0,0 +1,82 @@+--------------------------------------------------------------------------------+-- |+-- Module : LevMar.Intermediate.Fitting+-- 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.+--+-- 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+ ( -- * Model & Jacobian.+ Model+ , Jacobian++ -- * 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+++--------------------------------------------------------------------------------+-- Model & Jacobian.+--------------------------------------------------------------------------------++type Model r a = [r] -> a -> r++-- | See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>+type Jacobian r a = [r] -> a -> [r]+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm specialised for curve-fitting.+levmar :: LMA_I.LevMarable r+ => Model r a -- ^ Model+ -> Maybe (Jacobian r a) -- ^ Optional jacobian+ -> [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 mJac ps samples =+ LMA_I.levmar (\ps' -> map (model ps') xs)+ (fmap (\jac -> \ps' -> map (jac ps') xs) mJac)+ ps+ ys+ where+ (xs, ys) = unzip samples+++-- The End ---------------------------------------------------------------------
+ NFunction.hs view
@@ -0,0 +1,54 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}++module NFunction+ ( NFunction+ , ($*)+ , ComposeN+ , compose+ ) where++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@.+-- For example: @NFunction (S (S (S Z))) a b ~ (a -> a -> a -> b)@+type family NFunction n a b :: *++type instance NFunction Z a b = b+type instance NFunction (S n) a b = a -> NFunction n a b++-- | @f $* xs@ applies the /n/-arity function @f@ to each of the arguments in+-- the /n/-sized list @xs@.+($*) :: NFunction n a b -> SizedList n a -> b+f $* Nil = f+f $* (x ::: xs) = f x $* xs++infixr 0 $* -- same as $++class Nat n => ComposeN n where+ -- | Composition of NFunctions.+ --+ -- Note that the @n@ and @a@ arguments are used by the type+ -- checker to select the right @ComposeN@ instance. They are+ -- usally given as @(witnessNat :: n)@ and @(undefined :: a)@.+ compose :: forall a b c. n -> a+ -> (b -> c) -> NFunction n a b -> NFunction n a c++instance ComposeN Z where+ compose Z _ = ($)++instance ComposeN n => ComposeN (S n) where+ compose (S n) (_ :: a) f g = compose n (undefined :: a) f . g++{-+TODO: The following does not work as expected.+See: http://www.haskell.org/pipermail/haskell-cafe/2009-August/065850.html++-- | @f .* g@ composes @f@ with the /n/-arity function @g@.+(.*) :: forall n a b c. (ComposeN n) => (b -> c) -> NFunction n a b -> NFunction n a c+(.*) = compose (witnessNat :: n) (undefined :: a)++infixr 9 .* -- same as .+-}
+ Setup.hs view
@@ -0,0 +1,3 @@+import Distribution.Simple++main = defaultMain
+ SizedList.hs view
@@ -0,0 +1,76 @@+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ScopedTypeVariables #-}++module SizedList+ ( SizedList(..)+ , toList+ , fromList+ , unsafeFromList+ , length+ , replicate+ ) where++import Prelude hiding (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++infixr 5 ::: -- Same precedence and associativity as (:)++consPrecedence :: Int+consPrecedence = 5++instance Show a => Show (SizedList n a) where+ showsPrec _ Nil = showString "Nil"+ showsPrec p (x ::: xs) = showParen (p > consPrecedence)+ $ showsPrec (consPrecedence + 1) x+ . 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)+ where+ tl0 :: SizedList Z a -> [a]+ tl0 Nil = []++ tlS :: forall x. Nat x => (SizedList x a -> [a]) -> SizedList (S x) a -> [a]+ tlS f (x ::: xs) = x : f xs++newtype FromList a n = FromList { unFromList :: [a] -> Maybe (SizedList n a) }++-- | Convert a normal list to a @SizeList@. 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)+ where+ fl0 [] = Just Nil+ fl0 _ = Nothing++ flS _ [] = Nothing+ flS k (x:xs) = fmap (x :::) $ k xs++-- | Convert a normal list to a @SizeList@. If the length of the given+-- list does not equal @n@, an error is thrown.+unsafeFromList :: forall a n. Nat n => [a] -> SizedList n a+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
+ TypeLevelNat.hs view
@@ -0,0 +1,99 @@+-- Thanks to Ryan Ingram who wrote most of this module.+-- See: http://www.haskell.org/pipermail/haskell-cafe/2009-August/065674.html++{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}++module TypeLevelNat+ ( Z(..)+ , S(..)+ , Nat+ , caseNat+ , induction+ , witnessNat++ , N(..)+ ) where+++-- | Type-level natural denoting zero+data Z = Z deriving Show++-- | Type-level natural denoting the /S/uccessor of another type-level natural.+newtype S n = S n deriving Show++-- | Class of all type-level naturals.+class Nat n where+ -- | Case analysis on natural numbers.+ caseNat :: forall r.+ n -- ^ The natural number to case analyse.+ -> (n ~ Z => r) -- ^ The result @r@ when @n@ equals zero.+ -> (forall p. (n ~ S p, Nat p) => p -> r) -- ^ Function to apply to the predecessor+ -- of @n@ to yield the result @r@.+ -> r++instance Nat Z where+ caseNat _ z _ = z++instance Nat p => Nat (S p) where+ caseNat (S p) _ s = s p++-- | The axiom of induction on natural numbers.+-- See: <http://en.wikipedia.org/wiki/Mathematical_induction#Axiom_of_induction>+induction :: forall p n. Nat n+ => n+ -> p Z+ -> (forall m. Nat m => p m -> p (S m))+ -> p n+induction n z s = caseNat n isZ isS+ where+ isZ :: n ~ Z => p n+ isZ = z++ isS :: forall m. (n ~ S m, Nat m) => m -> p n+ isS m = s (induction m z s)++newtype Witness x = Witness { unWitness :: x }++-- | The value of @witnessNat :: n@ is the natural number of type @n@.+-- For example:+--+-- @+-- *TypeLevelNat> witnessNat :: S (S (S Z))+-- S (S (S Z))+-- @+witnessNat :: forall n. Nat n => n+witnessNat = theWitness+ where+ theWitness = unWitness $ induction (undefined `asTypeOf` theWitness)+ (Witness Z)+ (Witness . S . unWitness)++-- | A value-level natural indexed with an equivalent type-level natural.+data N n where+ Zero :: N Z+ Succ :: N n -> N (S n)++{-+Template Haskell code to construct a type synonym for an arbitrary+type level natural number.++Instead of++> type N6 = S (S (S (S (S (S Z)))))++you can write++> $(mkNat "N6" 6)+-}++-- import Language.Haskell.TH.Syntax++-- mkNat :: String -> Int -> Q [Dec]+-- mkNat syn = runQ . return . (:[]) . TySynD (mkName syn) [] . go+-- where go 0 = ConT $ mkName "Z"+-- go n = AppT (ConT $ mkName "S") $ go (n - 1)+
+ levmar.cabal view
@@ -0,0 +1,79 @@+name: levmar+version: 0.1+cabal-version: >= 1.6+build-type: Simple+stability: experimental+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+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+ is expressed as the sum of squares of nonlinear+ functions. It has become a standard technique for+ nonlinear least-squares problems and can be thought of+ as a combination of steepest descent and the+ Gauss-Newton method. When the current solution is far+ from the correct one, the algorithm behaves like a+ steepest descent method: slow, but guaranteed to+ converge. When the current solution is close to the+ correct solution, it becomes a Gauss-Newton method.+ .+ Optional box- and linear constraints can be given. Both+ single and double precision floating point types are+ supported.+ .+ The actual algorithm is implemented in a C library+ which is bundled with bindings-levmar which this+ package depends on. See:+ <http://www.ics.forth.gr/~lourakis/levmar/>.+ .+ This library consists of two layers:+ .+ * LevMar.Intermediate: A medium-level layer that wraps+ the low-level functions from bindings-levmar to+ provide a more Haskell friendly interface.+ .+ * LevMar: A high-level layer that uses type-level+ programming to add extra type safety.+ .+ Each layer also has special data-fitting variants:+ .+ * LevMar.Intermediate.Fitting+ .+ * LevMar.Fitting+ .+ All modules are self-contained; i.e. each module+ re-exports all the things you need to work with it.+ .+ For an example how to use this library see Demo.hs+ which is included in this package. Demo.hs is a Haskell+ translation of lmdemo.c from the C levmar library.+ .+ A note regarding the license:+ .+ This library depends on bindings-levmar which is+ bundled together with a C library which falls under the+ GPL. Please be aware of this when distributing programs+ linked with this library. For details see the+ description and license of bindings-levmar.+extra-source-files: Demo.hs++source-repository head+ Type: darcs+ Location: http://code.haskell.org/levmar++library+ build-depends: base >= 3 && < 4.2+ , bindings-levmar < 0.2+ exposed-modules: LevMar+ , LevMar.Fitting+ , LevMar.Intermediate+ , LevMar.Intermediate.Fitting+ , TypeLevelNat+ , SizedList+ , NFunction+ ghc-options: -Wall -O2