diff --git a/Bindings/LevMar/CurryFriendly.hs b/Bindings/LevMar/CurryFriendly.hs new file mode 100644 --- /dev/null +++ b/Bindings/LevMar/CurryFriendly.hs @@ -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 diff --git a/Demo.hs b/Demo.hs deleted file mode 100644 --- a/Demo.hs +++ /dev/null @@ -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: !! --- !! : (==): No overloading for function !! --- !! : interrupted !! --- !! : 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 . - -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 --------------------------------------------------------------------- diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -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. diff --git a/LevMar.hs b/LevMar.hs deleted file mode 100644 --- a/LevMar.hs +++ /dev/null @@ -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: --- --- --------------------------------------------------------------------------------- - -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: - -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 --------------------------------------------------------------------- diff --git a/LevMar/AD.hs b/LevMar/AD.hs deleted file mode 100644 --- a/LevMar/AD.hs +++ /dev/null @@ -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: --- --- --------------------------------------------------------------------------------- - - -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 --------------------------------------------------------------------- diff --git a/LevMar/Fitting.hs b/LevMar/Fitting.hs deleted file mode 100644 --- a/LevMar/Fitting.hs +++ /dev/null @@ -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: --- --- --------------------------------------------------------------------------------- - -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: - -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 --------------------------------------------------------------------- diff --git a/LevMar/Fitting/AD.hs b/LevMar/Fitting/AD.hs deleted file mode 100644 --- a/LevMar/Fitting/AD.hs +++ /dev/null @@ -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: --- --- --------------------------------------------------------------------------------- - -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 --------------------------------------------------------------------- diff --git a/LevMar/Intermediate.hs b/LevMar/Intermediate.hs deleted file mode 100644 --- a/LevMar/Intermediate.hs +++ /dev/null @@ -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: --- --- --------------------------------------------------------------------------------- - -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: - - * 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 --------------------------------------------------------------------- diff --git a/LevMar/Intermediate/AD.hs b/LevMar/Intermediate/AD.hs deleted file mode 100644 --- a/LevMar/Intermediate/AD.hs +++ /dev/null @@ -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: --- --- --------------------------------------------------------------------------------- - -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 --------------------------------------------------------------------- diff --git a/LevMar/Intermediate/Fitting.hs b/LevMar/Intermediate/Fitting.hs deleted file mode 100644 --- a/LevMar/Intermediate/Fitting.hs +++ /dev/null @@ -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: --- --- --------------------------------------------------------------------------------- - -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: - - * 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 --------------------------------------------------------------------- diff --git a/LevMar/Intermediate/Fitting/AD.hs b/LevMar/Intermediate/Fitting/AD.hs deleted file mode 100644 --- a/LevMar/Intermediate/Fitting/AD.hs +++ /dev/null @@ -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: --- --- --------------------------------------------------------------------------------- - -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 --------------------------------------------------------------------- diff --git a/LevMar/Utils.hs b/LevMar/Utils.hs deleted file mode 100644 --- a/LevMar/Utils.hs +++ /dev/null @@ -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 - ) diff --git a/LevMar/Utils/AD.hs b/LevMar/Utils/AD.hs deleted file mode 100644 --- a/LevMar/Utils/AD.hs +++ /dev/null @@ -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 diff --git a/NFunction.hs b/NFunction.hs deleted file mode 100644 --- a/NFunction.hs +++ /dev/null @@ -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 . --} diff --git a/Numeric/LevMar.hs b/Numeric/LevMar.hs new file mode 100644 --- /dev/null +++ b/Numeric/LevMar.hs @@ -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 +-- Bas van Dijk +-- Stability: Experimental +-- +-- For additional documentation see the +-- +-- 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 +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> 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> [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) diff --git a/README.markdown b/README.markdown new file mode 100644 --- /dev/null +++ b/README.markdown @@ -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 diff --git a/SizedList.hs b/SizedList.hs deleted file mode 100644 --- a/SizedList.hs +++ /dev/null @@ -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 --------------------------------------------------------------------- diff --git a/TypeLevelNat.hs b/TypeLevelNat.hs deleted file mode 100644 --- a/TypeLevelNat.hs +++ /dev/null @@ -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: -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) - diff --git a/levmar.cabal b/levmar.cabal --- a/levmar.cabal +++ b/levmar.cabal @@ -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 + Bas van Dijk +maintainer: Roel van Dijk + Bas van Dijk +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: - . - . - 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 + 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