levmar 0.2.1 → 1.2.1.8
raw patch · 19 files changed
Files
- Bindings/LevMar/CurryFriendly.hs +148/−0
- Demo.hs +0/−1451
- LICENSE +1/−1
- LevMar.hs +0/−150
- LevMar/AD.hs +0/−142
- LevMar/Fitting.hs +0/−151
- LevMar/Fitting/AD.hs +0/−139
- LevMar/Intermediate.hs +0/−409
- LevMar/Intermediate/AD.hs +0/−106
- LevMar/Intermediate/Fitting.hs +0/−131
- LevMar/Intermediate/Fitting/AD.hs +0/−106
- LevMar/Utils.hs +0/−44
- LevMar/Utils/AD.hs +0/−42
- NFunction.hs +0/−54
- Numeric/LevMar.hs +584/−0
- README.markdown +26/−0
- SizedList.hs +0/−109
- TypeLevelNat.hs +0/−103
- levmar.cabal +43/−90
+ Bindings/LevMar/CurryFriendly.hs view
@@ -0,0 +1,148 @@+{-# LANGUAGE NoImplicitPrelude #-}++module Bindings.LevMar.CurryFriendly+ ( -- * Handy type synonyms used in the curry friendly types.+ BoxConstraints+ , LinearConstraints++ -- * Curry friendly types of the Levenberg-Marquardt algorithms.+ , LevMarDer+ , LevMarDif+ , LevMarBCDer+ , LevMarBCDif+ , LevMarLecDer+ , LevMarLecDif+ , LevMarBLecDer+ , LevMarBLecDif++ -- * Curry friendly variants of the Levenberg-Marquardt+ -- algorithms in 'Bindings.Levmar'.+ , dlevmar_der, slevmar_der+ , dlevmar_dif, slevmar_dif+ , dlevmar_bc_der, slevmar_bc_der+ , dlevmar_bc_dif, slevmar_bc_dif+ , dlevmar_lec_der, slevmar_lec_der+ , dlevmar_lec_dif, slevmar_lec_dif+ , dlevmar_blec_der, slevmar_blec_der+ , dlevmar_blec_dif, slevmar_blec_dif+ ) where++import Prelude ( Double, Float )+import Foreign.Ptr ( FunPtr )++import qualified Bindings.LevMar as BLM+++--------------------------------------------------------------------------------+-- Handy type synonyms used in the curry friendly types.+--------------------------------------------------------------------------------++type BoxConstraints r a = BLM.LowerBounds r+ -> BLM.UpperBounds r+ -> a++type LinearConstraints r a = BLM.ConstraintsMatrix r+ -> BLM.ConstraintsVector r+ -> BLM.NrOfConstraints+ -> a+++--------------------------------------------------------------------------------+-- Curry friendly types of the Levenberg-Marquardt algorithms.+--------------------------------------------------------------------------------++type LevMarDif r = BLM.LevMarDif r+type LevMarDer r = FunPtr (BLM.Jacobian r) -> LevMarDif r+type LevMarBCDif r = BoxConstraints r (LevMarDif r)+type LevMarBCDer r = BoxConstraints r (LevMarDer r)+type LevMarLecDif r = LinearConstraints r (LevMarDif r)+type LevMarLecDer r = LinearConstraints r (LevMarDer r)+type LevMarBLecDif r = BoxConstraints r (LinearConstraints r (BLM.Weights r -> LevMarDif r))+type LevMarBLecDer r = BoxConstraints r (LinearConstraints r (BLM.Weights r -> LevMarDer r))+++--------------------------------------------------------------------------------+-- Reordering arguments to create curry friendly variants.+--------------------------------------------------------------------------------++mk_levmar_der :: BLM.LevMarDer r -> LevMarDer r+mk_levmar_der lma j f+ = lma f j++mk_levmar_bc_dif :: BLM.LevMarBCDif r -> LevMarBCDif r+mk_levmar_bc_dif lma lb ub f p x m n+ = lma f p x m n lb ub++mk_levmar_bc_der :: BLM.LevMarBCDer r -> LevMarBCDer r+mk_levmar_bc_der lma lb ub j f p x m n+ = lma f j p x m n lb ub++mk_levmar_lec_dif :: BLM.LevMarLecDif r -> LevMarLecDif r+mk_levmar_lec_dif lma a b k f p x m n+ = lma f p x m n a b k++mk_levmar_lec_der :: BLM.LevMarLecDer r -> LevMarLecDer r+mk_levmar_lec_der lma a b k j f p x m n+ = lma f j p x m n a b k++mk_levmar_blec_dif :: BLM.LevMarBLecDif r -> LevMarBLecDif r+mk_levmar_blec_dif lma lb ub a b k wghts f p x m n+ = lma f p x m n lb ub a b k wghts++mk_levmar_blec_der :: BLM.LevMarBLecDer r -> LevMarBLecDer r+mk_levmar_blec_der lma lb ub a b k wghts j f p x m n+ = lma f j p x m n lb ub a b k wghts+++--------------------------------------------------------------------------------+-- Curry friendly variants of the Levenberg-Marquardt algorithms in+-- 'Bindings.Levmar'.+--------------------------------------------------------------------------------++slevmar_dif :: LevMarDif Float+slevmar_dif = BLM.c'slevmar_dif++dlevmar_dif :: LevMarDif Double+dlevmar_dif = BLM.c'dlevmar_dif++slevmar_der :: LevMarDer Float+slevmar_der = mk_levmar_der BLM.c'slevmar_der++dlevmar_der :: LevMarDer Double+dlevmar_der = mk_levmar_der BLM.c'dlevmar_der++slevmar_bc_dif :: LevMarBCDif Float+slevmar_bc_dif = mk_levmar_bc_dif BLM.c'slevmar_bc_dif++dlevmar_bc_dif :: LevMarBCDif Double+dlevmar_bc_dif = mk_levmar_bc_dif BLM.c'dlevmar_bc_dif++slevmar_bc_der :: LevMarBCDer Float+slevmar_bc_der = mk_levmar_bc_der BLM.c'slevmar_bc_der++dlevmar_bc_der :: LevMarBCDer Double+dlevmar_bc_der = mk_levmar_bc_der BLM.c'dlevmar_bc_der++slevmar_lec_dif :: LevMarLecDif Float+slevmar_lec_dif = mk_levmar_lec_dif BLM.c'slevmar_lec_dif++dlevmar_lec_dif :: LevMarLecDif Double+dlevmar_lec_dif = mk_levmar_lec_dif BLM.c'dlevmar_lec_dif++slevmar_lec_der :: LevMarLecDer Float+slevmar_lec_der = mk_levmar_lec_der BLM.c'slevmar_lec_der++dlevmar_lec_der :: LevMarLecDer Double+dlevmar_lec_der = mk_levmar_lec_der BLM.c'dlevmar_lec_der++slevmar_blec_dif :: LevMarBLecDif Float+slevmar_blec_dif = mk_levmar_blec_dif BLM.c'slevmar_blec_dif++dlevmar_blec_dif :: LevMarBLecDif Double+dlevmar_blec_dif = mk_levmar_blec_dif BLM.c'dlevmar_blec_dif++slevmar_blec_der :: LevMarBLecDer Float+slevmar_blec_der = mk_levmar_blec_der BLM.c'slevmar_blec_der++dlevmar_blec_der :: LevMarBLecDer Double+dlevmar_blec_der = mk_levmar_blec_der BLM.c'dlevmar_blec_der
− Demo.hs
@@ -1,1451 +0,0 @@--- This module is a Haskell translation of lmdemo.c from the C levmar library.--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 ---------------------------------------------------------------------
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2009 Roel van Dijk, Bas van Dijk+Copyright (c) 2009-2014 Roel van Dijk, Bas van Dijk All rights reserved.
− LevMar.hs
@@ -1,150 +0,0 @@-{-# 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 LevMar.Utils ( LinearConstraints- , noLinearConstraints- , Matrix- , CovarMatrix- , convertLinearConstraints- , convertResult- )--import TypeLevelNat ( Z, S, Nat )-import SizedList ( SizedList(..), toList, unsafeFromList )-import NFunction ( NFunction, ($*) )--import Data.Either-------------------------------------------------------------------------------------- Model & Jacobian.-----------------------------------------------------------------------------------{- | 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 N4 Double-hatfldc p0 p1 p2 p3 = p0 - 1.0- ::: p0 - sqrt p1- ::: p1 - sqrt p2- ::: p3 - 1.0- ::: Nil-@--}-type Model m n r = NFunction m r (SizedList n 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>--For example the jacobian of the above @hatfldc@ model is:--@-type N4 = 'S' ('S' ('S' ('S' 'Z')))--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 m n r = NFunction m r (Matrix n m r)-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm.-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 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)- (toList ys)- itMax- opts- (fmap toList mLowBs)- (fmap toList mUpBs)- (fmap convertLinearConstraints mLinC)- (fmap toList mWghts)- where- 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
@@ -1,142 +0,0 @@-{-# 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
@@ -1,151 +0,0 @@-{-# 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- , SimpleModel- , Jacobian- , SimpleJacobian-- -- * 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.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, ($*) )-------------------------------------------------------------------------------------- Model & Jacobian.-----------------------------------------------------------------------------------{- | 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:--@-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 m r a = NFunction m r (a -> 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 +-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 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 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_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- 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
@@ -1,139 +0,0 @@-{-# 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
@@ -1,409 +0,0 @@-{-# 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.-----------------------------------------------------------------------------------{- | 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]--{- | 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]]-------------------------------------------------------------------------------------- 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) $ 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--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- = 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.- | 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, 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)- , (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/AD.hs
@@ -1,106 +0,0 @@-{-# 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- , LMA_I.Jacobian- , jacobianOf-- -- * 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---- | Compute the 'LMA_I.Jacobian' of the 'LMA_I.Model' using Automatic--- Differentiation.-jacobianOf :: (HasBasis r, Basis r ~ (), VectorSpace (Scalar r))- => 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
@@ -1,131 +0,0 @@------------------------------------------------------------------------------------ |--- 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- , SimpleModel- , Jacobian- , SimpleJacobian-- -- * 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.-----------------------------------------------------------------------------------{- | 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---- | 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.------------------------------------------------------------------------------------- | 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 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
@@ -1,106 +0,0 @@-{-# 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- , LMA_I.Jacobian- , LMA_I.SimpleJacobian- , jacobianOf-- -- * 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)---- | Compute the 'LMA_I.Jacobian' of the 'LMA_I.Model' using Automatic--- Differentiation.-jacobianOf :: (HasBasis r, Basis r ~ (), VectorSpace (Scalar r))- => 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
@@ -1,44 +0,0 @@-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
@@ -1,42 +0,0 @@-{-# 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
@@ -1,54 +0,0 @@-{-# 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 .--}
+ Numeric/LevMar.hs view
@@ -0,0 +1,584 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE UndecidableInstances #-}++--------------------------------------------------------------------------------+-- |+-- Module: Numeric.LevMar+-- Copyright: (c) 2009 - 2014 Roel van Dijk & Bas van Dijk+-- License: BSD-style (see the file LICENSE)+-- Maintainer: Roel van Dijk <vandijk.roel@gmail.com>+-- Bas van Dijk <v.dijk.bas@gmail.com>+-- Stability: Experimental+--+-- For additional documentation see the+-- <http://www.ics.forth.gr/~lourakis/levmar/ documentation of the levmar C>+-- library which this library is based on:+--+--------------------------------------------------------------------------------++module Numeric.LevMar+ ( -- * Model & Jacobian.+ Params, Samples+ , Model, Jacobian++ -- * Levenberg-Marquardt algorithm.+ , LevMarable(levmar)++ -- * Minimization options.+ , Options(..), defaultOpts++ -- * Constraints+ , Constraints(..), LinearConstraints++ -- * Output+ , Info(..), StopReason(..), LevMarError(..)+ ) where+++-------------------------------------------------------------------------------+-- Imports+-------------------------------------------------------------------------------++-- from base:+import Control.Monad ( return, mplus )+import Control.Exception ( Exception )+import Data.Bool ( (&&), (||), otherwise )+import Data.Data ( Data )+import Data.Typeable ( Typeable )+import Data.Either ( Either(Left, Right) )+import Data.Eq ( Eq, (==), (/=) )+import Data.Function ( (.), ($) )+import Data.Functor ( (<$>) )+import Data.Int ( Int )+import Data.List ( lookup, (++) )+import Data.Maybe ( Maybe(Nothing, Just), isJust, fromJust, fromMaybe )+import Data.Monoid ( Monoid, mempty, mappend )+import Data.Ord ( Ord, (<) )+import Foreign.C.Types ( CInt )+import Foreign.Marshal.Array ( allocaArray, withArray, peekArray, copyArray )+import Foreign.Ptr ( Ptr, nullPtr )+import Foreign.ForeignPtr ( ForeignPtr, newForeignPtr_, withForeignPtr )+import Foreign.Storable ( Storable )+import Prelude ( Num, Enum, Fractional, RealFrac, Float, Double+ , fromIntegral, toEnum, (-), (*), error, floor+ )+import System.IO ( IO )+import System.IO.Unsafe ( unsafePerformIO )+import Text.Read ( Read )+import Text.Show ( Show, show )++#if __GLASGOW_HASKELL__ >= 605+import GHC.ForeignPtr ( mallocPlainForeignPtrBytes )+import Prelude ( undefined )+import Foreign.Storable ( sizeOf )+#else+import Foreign.ForeignPtr ( mallocForeignPtrArray )+#endif++#if __GLASGOW_HASKELL__ < 700+import Prelude ( fromInteger, (>>=), (>>), fail )+#endif++-- from hmatrix:+#if MIN_VERSION_hmatrix(0,17,0)+import Numeric.LinearAlgebra.Data ( Matrix, flatten, rows, reshape )+import Numeric.LinearAlgebra ( Container, Element )+#else+import Data.Packed.Matrix ( Matrix, Element, flatten, rows, reshape )+import Numeric.Container ( Container )+import Numeric.LinearAlgebra ( {- Instances for Matrix -} )+#endif++-- from vector:+import Data.Vector.Storable ( Vector )+import qualified Data.Vector.Storable as VS ( unsafeWith, length+ , unsafeFromForeignPtr+ , length+ )++-- from bindings-levmar:+import Bindings.LevMar ( c'LM_INFO_SZ++ , withModel+ , withJacobian++ , c'LM_ERROR+ , c'LM_ERROR_LAPACK_ERROR+ , c'LM_ERROR_FAILED_BOX_CHECK+ , c'LM_ERROR_MEMORY_ALLOCATION_FAILURE+ , c'LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS+ , c'LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK+ , c'LM_ERROR_TOO_FEW_MEASUREMENTS+ , c'LM_ERROR_SINGULAR_MATRIX+ , c'LM_ERROR_SUM_OF_SQUARES_NOT_FINITE++ , c'LM_INIT_MU+ , c'LM_STOP_THRESH+ , c'LM_DIFF_DELTA+ )+import qualified Bindings.LevMar ( Model, Jacobian )++-- from levmar (this package):+import Bindings.LevMar.CurryFriendly ( LevMarDer, LevMarDif+ , LevMarBCDer, LevMarBCDif+ , LevMarLecDer, LevMarLecDif+ , LevMarBLecDer, LevMarBLecDif++ , dlevmar_der, slevmar_der+ , dlevmar_dif, slevmar_dif+ , dlevmar_bc_der, slevmar_bc_der+ , dlevmar_bc_dif, slevmar_bc_dif+ , dlevmar_lec_der, slevmar_lec_der+ , dlevmar_lec_dif, slevmar_lec_dif+ , dlevmar_blec_der, slevmar_blec_der+ , dlevmar_blec_dif, slevmar_blec_dif+ )+++--------------------------------------------------------------------------------+-- Model & Jacobian.+--------------------------------------------------------------------------------++-- | Parameter vector of length @m@.+--+-- Ensure that @m <= n@ where @n@ is the length of the 'Samples' vector.+type Params r = Vector r++-- | Sample vector of length @n@.+--+-- Ensure that @n >= m@ where @m@ is the length of the 'Params' vector.+type Samples r = Vector r++{-| A functional relation describing measurements represented as a function+from a vector of parameters to a vector of expected samples.++ * Ensure that the length @m@ of the parameter vector equals the length of the+ initial parameter vector in 'levmar'.++ * Ensure that the length @n@ of the output sample vector equals the length of+ the sample vector in 'levmar'.++ * Ensure that the length @n@ of the output sample vector vector is bigger than or+ equal to the length @m@ of the parameter vector.+-}+type Model r = Params r -> Samples r++{-| The <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant jacobian>+of the 'Model' function. Expressed as a function from a vector of+parameters to a matrix which for each expected sample describes the+partial derivatives of the parameters.++ * Ensure that the length @m@ of the parameter vector equals the length of the+ initial parameter vector in 'levmar'.++ * Ensure that the output matrix has the dimension @n><m@ where @n@ is the+ number of samples and @m@ is the number of parameters.+-}+type Jacobian r = Params r -> Matrix r++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm is overloaded to work on 'Double' and 'Float'.+class LevMarable r where++ -- | The Levenberg-Marquardt algorithm.+ --+ -- Returns a triple of the found parameters, a structure containing+ -- information about the minimization and the covariance matrix+ -- corresponding to LS solution.+ --+ -- Ensure that @n >= m@.+ levmar :: Model r -- ^ Model+ -> Maybe (Jacobian r) -- ^ Optional jacobian+ -> Params r -- ^ Initial parameters of length @m@+ -> Samples r -- ^ Sample vector of length @n@+ -> Int -- ^ Maximum iterations+ -> Options r -- ^ Minimization options+ -> Constraints r -- ^ Constraints+ -> Either LevMarError (Params r, Info r, Matrix r)++instance LevMarable Float where+ levmar = gen_levmar slevmar_der+ slevmar_dif+ slevmar_bc_der+ slevmar_bc_dif+ slevmar_lec_der+ slevmar_lec_dif+ slevmar_blec_der+ slevmar_blec_dif++instance LevMarable Double where+ levmar = gen_levmar dlevmar_der+ dlevmar_dif+ dlevmar_bc_der+ dlevmar_bc_dif+ dlevmar_lec_der+ dlevmar_lec_dif+ dlevmar_blec_der+ 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 r. (RealFrac r, Element r)+ => LevMarDer r+ -> LevMarDif r+ -> LevMarBCDer r+ -> LevMarBCDif r+ -> LevMarLecDer r+ -> LevMarLecDif r+ -> LevMarBLecDer r+ -> LevMarBLecDif r++ -> Model r -- ^ Model+ -> Maybe (Jacobian r) -- ^ Optional jacobian+ -> Params r -- ^ Initial parameters+ -> Samples r -- ^ Samples+ -> Int -- ^ Maximum iterations+ -> Options r -- ^ Options+ -> Constraints r -- ^ Constraints+ -> Either LevMarError (Params r, Info r, Matrix 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 (Constraints mLowBs mUpBs mWeights mLinC)+ | m == 0 = Left LevMarError -- LAPACK will crash otherwise!+ | otherwise =+ -- All effects are contained, so we can safely perform:+ unsafePerformIO $ do++ -- We need to pass the initial parameters 'ps' to the C function.+ -- However, we can't just pass a pointer to them because the C function+ -- will modify the parameters during execution which will violate+ -- referential transparanency. Instead we allocate new space+ -- and copy the parameters to it.+ --+ -- Note that, in the end, the array is returned from this function.+ -- This means that the only way to guarantee its finalisation+ -- is to allocate it using a ForeignPtr:+ psFP <- fastMallocForeignPtrArray m+ withForeignPtr psFP $ \psPtr -> do+ VS.unsafeWith ps $ \psPtrInp ->+ copyArray psPtr psPtrInp m++ -- Retrieve the (read-only) pointer 'ysPtr' to the samples vector 'ys'+ -- so we can pass it to the C function:+ VS.unsafeWith ys $ \ysPtr ->++ -- Convert the Options 'opts' to a list and then to an array+ -- so we can pass the (read-only) pointer 'optsPtr' to the C function:+ withArray (optsToList opts) $ \optsPtr ->++ -- Allocate space for the info array+ -- so we can pass it to the C function.+ -- Note that, in the end, this array is converted to an Info value+ -- and returned from this function.+ allocaArray c'LM_INFO_SZ $ \infoPtr -> do++ -- Allocate space for the covariance matrix+ -- so we can pass it to the C function.+ -- Like the parameters array the matrix+ -- needs to be returned from this function.+ -- So we also allocate it using a ForeignPtr:+ covarFP <- fastMallocForeignPtrArray mm+ withForeignPtr covarFP $ \covarPtr ->++ -- 'cmodel' is the low-level model function which is converted+ -- to the FunPtr 'modelFunPtr' and passed to the C function.+ -- 'cmodel' will first convert the parameters pointer 'parPtr'+ -- into a Vector after converting it into a ForeignPtr+ -- (without a finalizer).+ -- Then it will apply the high-level 'model' function+ -- to this parameter vector. The resulting vector is then copied+ -- to the output buffer 'hxPtr':+ let cmodel :: Bindings.LevMar.Model r+ cmodel parPtr hxPtr _ _ _ = do+ parFP <- newForeignPtr_ parPtr+ let psV = VS.unsafeFromForeignPtr parFP 0 m+ vector = model psV+ VS.unsafeWith vector $ \p -> copyArray hxPtr p (VS.length vector)+ in withModel cmodel $ \modelFunPtr -> do++ -- All the low-level C functions share a common set of arguments.+ -- 'runDif' applies these arguments to the given C function 'f':+ let runDif :: LevMarDif r -> IO CInt+ runDif f = f modelFunPtr+ psPtr+ ysPtr+ (fromIntegral m)+ (fromIntegral n)+ (fromIntegral itMax)+ optsPtr+ infoPtr+ nullPtr+ covarPtr+ nullPtr++ err <- case mJac of+ 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++ Just jac ->+ let cjacobian :: Bindings.LevMar.Jacobian r+ cjacobian parPtr jPtr _ _ _ = do+ parFP <- newForeignPtr_ parPtr+ let psV = VS.unsafeFromForeignPtr parFP 0 m+ matrix = jac psV+ vector = flatten matrix+ VS.unsafeWith vector $ \p ->+ copyArray jPtr p (VS.length vector)+ in withJacobian cjacobian $ \jacobPtr ->++ let runDer :: LevMarDer r -> 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++ -- Handling errors:+ if err < 0+ -- we don't treat the following two as an error:+ && err /= c'LM_ERROR_SINGULAR_MATRIX+ && err /= c'LM_ERROR_SUM_OF_SQUARES_NOT_FINITE+ then return $ Left $ convertLevMarError err++ else do -- Converting results:+ info <- listToInfo <$> peekArray c'LM_INFO_SZ infoPtr+ let psV = VS.unsafeFromForeignPtr psFP 0 m+ let covarM = reshape m $ VS.unsafeFromForeignPtr covarFP 0 mm++ return $ Right (psV, info, covarM)+ where+ m = VS.length ps+ n = VS.length ys+ mm = m*m++ -- Whether the parameters are constrained by a linear equation.+ linConstrained = isJust mLinC+ (cMat, rhcVec) = fromJust mLinC++ -- Whether the parameters are constrained by a bounding box.+ boxConstrained = isJust mLowBs || isJust mUpBs++ withBoxConstraints f g =+ maybeWithArray mLowBs $ \lBsPtr ->+ maybeWithArray mUpBs $ \uBsPtr ->+ f $ g lBsPtr uBsPtr++ withLinConstraints f g =+ VS.unsafeWith (flatten cMat) $ \cMatPtr ->+ VS.unsafeWith rhcVec $ \rhcVecPtr ->+ f . g cMatPtr rhcVecPtr $ fromIntegral $ rows cMat++ withWeights f g = maybeWithArray mWeights $ f . g++maybeWithArray :: (Storable a) => Maybe (Vector a) -> (Ptr a -> IO β) -> IO β+maybeWithArray Nothing f = f nullPtr+maybeWithArray (Just v) f = VS.unsafeWith v f++#if __GLASGOW_HASKELL__ >= 605+{-# INLINE fastMallocForeignPtrArray #-}+fastMallocForeignPtrArray :: forall a. Storable a => Int -> IO (ForeignPtr a)+fastMallocForeignPtrArray n = mallocPlainForeignPtrBytes+ (n * sizeOf (undefined :: a))+#else+fastMallocForeignPtrArray :: forall a. Storable a => Int -> IO (ForeignPtr a)+fastMallocForeignPtrArray = mallocForeignPtrArray+#endif+++--------------------------------------------------------------------------------+-- 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 (Eq, Ord, Read, Show, Data, Typeable)++-- | Default minimization options+defaultOpts :: Fractional r => Options r+defaultOpts = Opts { optScaleInitMu = c'LM_INIT_MU+ , optStopNormInfJacTe = c'LM_STOP_THRESH+ , optStopNorm2Dp = c'LM_STOP_THRESH+ , optStopNorm2E = c'LM_STOP_THRESH+ , optDelta = c'LM_DIFF_DELTA+ }++optsToList :: Options r -> [r]+optsToList (Opts mu eps1 eps2 eps3 delta) =+ [mu, eps1, eps2, eps3, delta]+++--------------------------------------------------------------------------------+-- Constraints+--------------------------------------------------------------------------------++-- | Ensure that these vectors have the same length as the number of parameters.+data Constraints r = Constraints+ { lowerBounds :: !(Maybe (Params r)) -- ^ Optional lower bounds+ , upperBounds :: !(Maybe (Params r)) -- ^ Optional upper bounds+ , weights :: !(Maybe (Params r)) -- ^ Optional weights+ , linearConstraints :: !(Maybe (LinearConstraints r)) -- ^ Optional linear constraints+ } deriving (Read, Show, Typeable)++deriving instance (Eq r, Container Vector r, Num r) => Eq (Constraints r)++-- | Linear constraints consisting of a constraints matrix, @k><m@ and+-- a right hand constraints vector, of length @k@ where @m@ is the number of+-- parameters and @k@ is the number of constraints.+type LinearConstraints r = (Matrix r, Vector r)++-- | * 'mempty' is defined as a 'Constraints' where all fields are 'Nothing'.+--+-- * 'mappend' merges two 'Constraints' by taking the first non-'Nothing' value+-- for each field.+instance Monoid (Constraints r) where+ mempty = Constraints Nothing Nothing Nothing Nothing+ mappend (Constraints lb1 ub1 w1 l1)+ (Constraints lb2 ub2 w2 l2) = Constraints (lb1 `mplus` lb2)+ (ub1 `mplus` ub2)+ (w1 `mplus` w2)+ (l1 `mplus` l2)+++--------------------------------------------------------------------------------+-- 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 :: !Int -- ^ Number of iterations.+ , infStopReason :: !StopReason -- ^ Reason for terminating.+ , infNumFuncEvals :: !Int -- ^ Number of function evaluations.+ , infNumJacobEvals :: !Int -- ^ Number of jacobian evaluations.+ , infNumLinSysSolved :: !Int -- ^ Number of linear systems solved,+ -- i.e. attempts for reducing error.+ } deriving (Eq, Ord, Read, Show, Data, Typeable)++listToInfo :: (RealFrac r) => [r] -> Info r+listToInfo [a,b,c,d,e,f,g,h,i,j] =+ Info { infNorm2initE = a+ , infNorm2E = b+ , infNormInfJacTe = c+ , infNorm2Dp = d+ , infMuDivMax = 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 (Eq, Ord, Read, Show, Data, Typeable, Enum)+++--------------------------------------------------------------------------------+-- Error+--------------------------------------------------------------------------------++data LevMarError+ = 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.+ | 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 (Eq, Ord, Read, Show, Data, Typeable)++-- Handy in case you want to thow a LevMarError as an exception:+instance Exception LevMarError++levmarCErrorToLevMarError :: [(CInt, LevMarError)]+levmarCErrorToLevMarError =+ [ (c'LM_ERROR, LevMarError)+ , (c'LM_ERROR_LAPACK_ERROR, LapackError)+ --, (c'LM_ERROR_NO_JACOBIAN, can never happen)+ --, (c'LM_ERROR_NO_BOX_CONSTRAINTS, can never happen)+ , (c'LM_ERROR_FAILED_BOX_CHECK, FailedBoxCheck)+ , (c'LM_ERROR_MEMORY_ALLOCATION_FAILURE, MemoryAllocationFailure)+ , (c'LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS, ConstraintMatrixRowsGtCols)+ , (c'LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK, ConstraintMatrixNotFullRowRank)+ , (c'LM_ERROR_TOO_FEW_MEASUREMENTS, TooFewMeasurements)+ --, (c'LM_ERROR_SINGULAR_MATRIX, we don't treat this as an error)+ --, (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: " ++ show err)+ (lookup err levmarCErrorToLevMarError)
+ README.markdown view
@@ -0,0 +1,26 @@+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.++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].++[bindings-levmar]: http://hackage.haskell.org/package/bindings-levmar+[C library]: http://www.ics.forth.gr/~lourakis/levmar
− SizedList.hs
@@ -1,109 +0,0 @@-{-# LANGUAGE GADTs #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE Rank2Types #-}--module SizedList- ( SizedList(..)- , foldr- , foldrN- , toList- , length- , fromList- , unsafeFromList- , replicate- ) where---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--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--------------------------------------------------------------------------------------- | 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- 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---- | 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---- | Convert a @SizedList@ to a normal list.-toList :: SizedList n a -> [a]-toList = foldr (:) []---- | 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 = unFL $ induction (witnessNat :: n) (FL flZ) (FL . flS . unFL)- where- flZ [] = Just Nil- flZ _ = 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-------------------------------------------------------------------------------------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
@@ -1,103 +0,0 @@--- 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(..)- , nat- ) 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)--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-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
@@ -1,102 +1,55 @@ name: levmar-version: 0.2.1+version: 1.2.1.8 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+author: Roel van Dijk <vandijk.roel@gmail.com>+ Bas van Dijk <v.dijk.bas@gmail.com>+maintainer: Roel van Dijk <vandijk.roel@gmail.com>+ Bas van Dijk <v.dijk.bas@gmail.com>+copyright: (c) 2009 - 2014 Roel van Dijk & Bas van Dijk license: BSD3 license-file: LICENSE+homepage: https://github.com/basvandijk/levmar+bug-reports: https://github.com/basvandijk/levmar/issues 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- 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 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.- .- 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+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+ <http://www.ics.forth.gr/~lourakis/levmar/ C library> which is+ bundled with @bindings-levmar@ which this package depends on.+ .+ 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: README.markdown+ source-repository head- Type: darcs- Location: http://code.haskell.org/levmar+ Type: git+ Location: git://github.com/basvandijk/levmar.git library- build-depends: base >= 3 && < 4.2- , bindings-levmar >= 0.1.1.1 && < 0.2- , 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+ build-depends: base >= 3 && < 5+ , bindings-levmar >= 1.1+ , hmatrix >= 0.12+ , vector >= 0.8+ exposed-modules: Numeric.LevMar+ other-modules: Bindings.LevMar.CurryFriendly+ ghc-options: -Wall