levmar 0.2.1 → 0.3

raw patch · 18 files changed

+863/−3225 lines, 18 filesdep +base-unicode-symbolsdep −MemoTriedep −vector-spacedep ~basedep ~bindings-levmarPVP ok

version bump matches the API change (PVP)

Dependencies added: base-unicode-symbols

Dependencies removed: MemoTrie, vector-space

Dependency ranges changed: base, bindings-levmar

API changes (from Hackage documentation)

Files

• Bindings/LevMar/CurryFriendly.hs +152/−0
• Demo.hs +0/−1451
• LevMar.hs +0/−150
• LevMar/AD.hs +0/−142
• LevMar/Fitting.hs +0/−151
• LevMar/Fitting/AD.hs +0/−139
• LevMar/Intermediate.hs +0/−409
• LevMar/Intermediate/AD.hs +0/−106
• LevMar/Intermediate/Fitting.hs +0/−131
• LevMar/Intermediate/Fitting/AD.hs +0/−106
• LevMar/Utils.hs +0/−44
• LevMar/Utils/AD.hs +0/−42
• NFunction.hs +0/−54
• Numeric/LevMar.hs +527/−0
• Numeric/LevMar/Fitting.hs +142/−0
• SizedList.hs +0/−109
• TypeLevelNat.hs +0/−103
• levmar.cabal +42/−88

@@ -0,0 +1,152 @@+{-# LANGUAGE NoImplicitPrelude, UnicodeSyntax #-}++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 Foreign.C.Types ( CFloat, CDouble )+import Foreign.Ptr     ( FunPtr )++import qualified Bindings.LevMar as BLM+++--------------------------------------------------------------------------------+-- Handy type synonyms used in the curry friendly types.+--------------------------------------------------------------------------------++type BoxConstraints cr a = BLM.LowerBounds cr+                         → BLM.UpperBounds cr+                         → a++type LinearConstraints cr a = BLM.ConstraintsMatrix cr+                            → BLM.ConstraintsVector cr+                            → BLM.NrOfConstraints+                            → a+++--------------------------------------------------------------------------------+-- Curry friendly types of the Levenberg-Marquardt algorithms.+--------------------------------------------------------------------------------++type LevMarDif     cr = BLM.LevMarDif cr+type LevMarDer     cr = FunPtr (BLM.Jacobian cr) → LevMarDif cr+type LevMarBCDif   cr = BoxConstraints cr (LevMarDif cr)+type LevMarBCDer   cr = BoxConstraints cr (LevMarDer cr)+type LevMarLecDif  cr = LinearConstraints cr (LevMarDif cr)+type LevMarLecDer  cr = LinearConstraints cr (LevMarDer cr)+type LevMarBLecDif cr = BoxConstraints cr (LinearConstraints cr (BLM.Weights cr → LevMarDif cr))+type LevMarBLecDer cr = BoxConstraints cr (LinearConstraints cr (BLM.Weights cr → LevMarDer cr))+++--------------------------------------------------------------------------------+-- Reordering arguments to create curry friendly variants.+--------------------------------------------------------------------------------++mk_levmar_der ∷ BLM.LevMarDer cr → LevMarDer cr+mk_levmar_der lma j f+            = lma f j++mk_levmar_bc_dif ∷ BLM.LevMarBCDif cr → LevMarBCDif cr+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 cr → LevMarBCDer cr+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 cr → LevMarLecDif cr+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 cr → LevMarLecDer cr+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 cr → LevMarBLecDif cr+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 cr → LevMarBLecDer cr+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 CFloat+slevmar_dif = BLM.c'slevmar_dif++dlevmar_dif ∷ LevMarDif CDouble+dlevmar_dif = BLM.c'dlevmar_dif++slevmar_der ∷ LevMarDer CFloat+slevmar_der = mk_levmar_der BLM.c'slevmar_der++dlevmar_der ∷ LevMarDer CDouble+dlevmar_der = mk_levmar_der BLM.c'dlevmar_der++slevmar_bc_dif ∷ LevMarBCDif CFloat+slevmar_bc_dif = mk_levmar_bc_dif BLM.c'slevmar_bc_dif++dlevmar_bc_dif ∷ LevMarBCDif CDouble+dlevmar_bc_dif = mk_levmar_bc_dif BLM.c'dlevmar_bc_dif++slevmar_bc_der ∷ LevMarBCDer CFloat+slevmar_bc_der = mk_levmar_bc_der BLM.c'slevmar_bc_der++dlevmar_bc_der ∷ LevMarBCDer CDouble+dlevmar_bc_der = mk_levmar_bc_der BLM.c'dlevmar_bc_der++slevmar_lec_dif ∷ LevMarLecDif CFloat+slevmar_lec_dif = mk_levmar_lec_dif BLM.c'slevmar_lec_dif++dlevmar_lec_dif ∷ LevMarLecDif CDouble+dlevmar_lec_dif = mk_levmar_lec_dif BLM.c'dlevmar_lec_dif++slevmar_lec_der ∷ LevMarLecDer CFloat+slevmar_lec_der = mk_levmar_lec_der BLM.c'slevmar_lec_der++dlevmar_lec_der ∷ LevMarLecDer CDouble+dlevmar_lec_der = mk_levmar_lec_der BLM.c'dlevmar_lec_der++slevmar_blec_dif ∷ LevMarBLecDif CFloat+slevmar_blec_dif = mk_levmar_blec_dif BLM.c'slevmar_blec_dif++dlevmar_blec_dif ∷ LevMarBLecDif CDouble+dlevmar_blec_dif = mk_levmar_blec_dif BLM.c'dlevmar_blec_dif++slevmar_blec_der ∷ LevMarBLecDer CFloat+slevmar_blec_der = mk_levmar_blec_der BLM.c'slevmar_blec_der++dlevmar_blec_der ∷ LevMarBLecDer CDouble+dlevmar_blec_der = mk_levmar_blec_der BLM.c'dlevmar_blec_der+++-- The End ---------------------------------------------------------------------

@@ -1,1451 +0,0 @@--- This module is a Haskell translation of lmdemo.c from the C levmar library.--module Main where--import LevMar ( levmar--              , Model-              , Jacobian--              , Options(..), defaultOpts--              , LinearConstraints, noLinearConstraints--              , LevMarError--              , Info(..), CovarMatrix--              , S, Z-              , SizedList(..)-              )---import qualified LevMar.AD         as AD-import qualified LevMar.Fitting    as Fitting-import qualified LevMar.Fitting.AD as Fitting.AD--import qualified SizedList as SL (replicate)-------------------------------------------------------------------------------------type Result n = Either LevMarError-                       ( SizedList n Double-                       , Info Double-                       , CovarMatrix n Double-                       )--printInteresting :: Result n -> IO ()-printInteresting (Left err) = putStrLn ("Error: " ++ show err)-printInteresting (Right (ps, inf, covar)) =-    do putStrLn ("infStopReason = " ++ show (infStopReason inf))-       putStrLn ("infNorm2E     = " ++ show (infNorm2E     inf))-       putStrLn ("infNumIter    = " ++ show (infNumIter    inf))-       putStrLn ("ps            = " ++ show ps)--sqr :: Num a => a -> a-sqr x = x*x------------------------------------------------------------------------------------- Handy type synonyms for type-level naturals:--type N0 = Z-type N1 = S N0-type N2 = S N1-type N3 = S N2-type N4 = S N3-type N5 = S N4-type N6 = S N5------------------------------------------------------------------------------------- Default options:--opts :: Options Double-opts = defaultOpts { optStopNormInfJacTe = 1e-15-                   , optStopNorm2Dp      = 1e-15-                   , optStopNorm2E       = 1e-20-                   }------------------------------------------------------------------------------------- Rosenbrock function,--- global minimum at (1, 1)--ros :: Floating r => Model N2 N2 r-ros p0 p1 = SL.replicate (sqr (1.0 - p0) + ros_d*sqr m)-    where-      m = p1 - sqr p0--ros_jac :: Floating r => Jacobian N2 N2 r-ros_jac p0 p1 = SL.replicate (  -2 + 2*p0 - 4*ros_d*m*p0-                             ::: 2*ros_d*m-                             ::: Nil-                             )-    where-      m = p1 - sqr p0--ros_d :: Floating r => r-ros_d = 105.0--ros_params :: Floating r => SizedList N2 r-ros_params = -1.2 ::: 1.0 ::: Nil--ros_samples :: Floating r => SizedList N2 r-ros_samples = SL.replicate 0.0---- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--- !! TODO: Find out why these return with: infStopReason = MaxIterations !!--- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--run_ros :: IO ()-run_ros = printInteresting $-          levmar ros-                 Nothing-                 ros_params-                 ros_samples-                 1000-                 opts-                 Nothing-                 Nothing-                 noLinearConstraints-                 Nothing--run_ros_jac :: IO ()-run_ros_jac = printInteresting $-              levmar ros-                     (Just ros_jac)-                     ros_params-                     ros_samples-                     1000-                     opts-                     Nothing-                     Nothing-                     noLinearConstraints-                     Nothing--run_ros_autojac :: IO ()-run_ros_autojac = printInteresting $-                  AD.levmar ros-                            ros_params-                            ros_samples-                            1000-                            opts-                            Nothing-                            Nothing-                            noLinearConstraints-                            Nothing------------------------------------------------------------------------------------- Modified Rosenbrock problem,--- global minimum at (1, 1)--modros :: Floating r => Model N2 N3 r-modros p0 p1 =     10*(p1 - sqr p0)-               ::: 1.0 - p0-               ::: modros_lam-               ::: Nil--modros_jac :: Floating r => Jacobian N2 N3 r-modros_jac p0 _ =     (-20*p0 ::: 10.0 ::: Nil)-                  ::: (-1.0   ::: 0.0  ::: Nil)-                  ::: (0.0    ::: 0.0  ::: Nil)-                  ::: Nil--modros_lam :: Floating r => r-modros_lam = 1e02--modros_params :: Floating r => SizedList N2 r-modros_params = -1.2 ::: 1.0 ::: Nil--modros_samples :: Floating r => SizedList N3 r-modros_samples = SL.replicate 0.0--run_modros :: IO ()-run_modros = printInteresting $-             levmar modros-                    Nothing-                    modros_params-                    modros_samples-                    1000-                    opts-                    Nothing-                    Nothing-                    noLinearConstraints-                    Nothing--run_modros_jac :: IO ()-run_modros_jac = printInteresting $-                 levmar modros-                        (Just modros_jac)-                        modros_params-                        modros_samples-                        1000-                        opts-                        Nothing-                        Nothing-                        noLinearConstraints-                        Nothing--run_modros_autojac :: IO ()-run_modros_autojac = printInteresting $-                     AD.levmar modros-                               modros_params-                               modros_samples-                               1000-                               opts-                               Nothing-                               Nothing-                               noLinearConstraints-                               Nothing------------------------------------------------------------------------------------- Powell's function,--- minimum at (0, 0)--powell :: Floating r => Model N2 N2 r-powell p0 p1 =     p0-               ::: 10.0*p0 / m + 2*sqr p1-               ::: Nil-    where-      m = p0 + 0.1--powell_jac :: Floating r => Jacobian N2 N2 r-powell_jac p0 p1 =     (1.0         ::: 0.0    ::: Nil)-                   ::: (1.0 / sqr m ::: 4.0*p1 ::: Nil)-                   ::: Nil-    where-      m = p0 + 0.1--powell_params :: Floating r => SizedList N2 r-powell_params = -1.2 ::: 1.0 ::: Nil--powell_samples :: Floating r => SizedList N2 r-powell_samples = SL.replicate 0.0--run_powell :: IO ()-run_powell = printInteresting $-             levmar powell-                    Nothing-                    powell_params-                    powell_samples-                    1000-                    opts-                    Nothing-                    Nothing-                    noLinearConstraints-                    Nothing--run_powell_jac :: IO ()-run_powell_jac = printInteresting $-                 levmar powell-                        (Just powell_jac)-                        powell_params-                        powell_samples-                        1000-                        opts-                        Nothing-                        Nothing-                        noLinearConstraints-                        Nothing---- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--- !! TODO: Here the automatic jacobian does not seem right because !!--- !! infNorm2E is very high compared to the manual jacobian!       !!--- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--run_powell_autojac :: IO ()-run_powell_autojac = printInteresting $-                     AD.levmar powell-                               powell_params-                               powell_samples-                               1000-                               opts-                               Nothing-                               Nothing-                               noLinearConstraints-                               Nothing------------------------------------------------------------------------------------- Wood's function,--- minimum at (1, 1, 1, 1)--wood :: Floating r => Model N4 N6 r-wood p0 p1 p2 p3 =     10.0*(p1 - sqr p0)-                   ::: 1.0 - p0-                   ::: sqrt 90.0*(p3 - sqr p2)-                   ::: 1.0 - p2-                   ::: sqrt 10.0*(p1 + p3 - 2.0)-                   ::: (p1 - p3) / sqrt 10.0-                   ::: Nil--wood_params :: Floating r => SizedList N4 r-wood_params =  -3.0 ::: -1.0 ::: -3.0 ::: -1.0 ::: Nil--wood_samples :: Floating r => SizedList N6 r-wood_samples = SL.replicate 0.0--run_wood :: IO ()-run_wood = printInteresting $-           levmar wood-                  Nothing-                  wood_params-                  wood_samples-                  1000-                  opts-                  Nothing-                  Nothing-                  noLinearConstraints-                  Nothing--run_wood_autojac :: IO ()-run_wood_autojac = printInteresting $-                   AD.levmar wood-                             wood_params-                             wood_samples-                             1000-                             opts-                             Nothing-                             Nothing-                             noLinearConstraints-                             Nothing------------------------------------------------------------------------------------- Meyer's (reformulated) data fitting problem,--- minimum at (2.48, 6.18, 3.45)--meyer :: Floating r => Fitting.SimpleModel N3 r-meyer p0 p1 p2 x = p0*exp (10.0*p1 / (ui + p2) - 13.0)-    where-      ui = 0.45 + 0.05*x--meyer_jac :: Floating r => Fitting.SimpleJacobian N3 r-meyer_jac p0 p1 p2 x =     tmp-                       ::: 10.0*p0*tmp / (ui + p2)-                       ::: -10.0*p0*p1*tmp / ((ui + p2)*(ui + p2))-                       ::: Nil-    where-      tmp = exp (10.0*p1 / (ui + p2) - 13.0)-      ui = 0.45 + 0.05*x--meyer_params :: Floating r => SizedList N3 r-meyer_params = 8.85 ::: 4.0 ::: 2.5 ::: Nil---- TODO: Unfortunately 'zip [0..] ...' won't work because (:~>)--- doesn't have an Enum instance:-meyer_samples :: (Num a, Floating r) => [(a, r)]-meyer_samples = [ ( 0, 34.780)-                , ( 1, 28.610)-                , ( 2, 23.650)-                , ( 3, 19.630)-                , ( 4, 16.370)-                , ( 5, 13.720)-                , ( 6, 11.540)-                , ( 7,  9.744)-                , ( 8,  8.261)-                , ( 9,  7.030)-                , (10,  6.005)-                , (11,  5.147)-                , (12,  4.427)-                , (13,  3.820)-                , (14,  3.307)-                , (15,  2.872)-                ]--run_meyer :: IO ()-run_meyer = printInteresting $-            Fitting.levmar meyer-                           Nothing-                           meyer_params-                           meyer_samples-                           1000-                           opts-                           Nothing-                           Nothing-                           noLinearConstraints-                           Nothing--run_meyer_jac :: IO ()-run_meyer_jac = printInteresting $-                Fitting.levmar meyer-                               (Just meyer_jac)-                               meyer_params-                               meyer_samples-                               1000-                               opts-                               Nothing-                               Nothing-                               noLinearConstraints-                               Nothing---- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--- !! TODO: Here the automatic jacobian does not seem right because !!--- !! infNorm2E is very high compared to the manual jacobian!       !!--- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--run_meyer_autojac :: IO ()-run_meyer_autojac = printInteresting $-                    Fitting.AD.levmar meyer-                                      meyer_params-                                      meyer_samples-                                      1000-                                      opts-                                      Nothing-                                      Nothing-                                      noLinearConstraints-                                      Nothing------------------------------------------------------------------------------------- helical valley function,--- minimum at (1.0, 0.0, 0.0)--helval :: (Ord r, Floating r) => Model N3 N3 r-helval p0 p1 p2 =     10.0*(p2 - 10.0*theta)-                  ::: 10.0*sqrt tmp - 1.0-                  ::: p2-                  ::: Nil-    where-      m = atan (p1 / p0) / (2.0*pi)--      tmp = sqr p0 + sqr p1--      theta | p0 < 0.0  = m + 0.5-            | 0.0 < p0  = m-            | p1 >= 0   = 0.25-            | otherwise = -0.25--heval_jac :: Floating r => Jacobian N3 N3 r-heval_jac p0 p1 _ =     (50.0*p1 / (pi*tmp) ::: -50.0*p0 / (pi*tmp) ::: 10.0 ::: Nil)-                    ::: (10.0*p0 / sqrt tmp :::  10.0*p1 / sqrt tmp ::: 0.0  ::: Nil)-                    ::: (0.0                :::  0.0                ::: 1.0  ::: Nil)-                    ::: Nil-    where-      tmp = sqr p0 + sqr p1--helval_params :: Floating r => SizedList N3 r-helval_params = -1.0 ::: 0.0 ::: 0.0 ::: Nil--helval_samples :: Floating r => SizedList N3 r-helval_samples = SL.replicate 0.0--run_helval :: IO ()-run_helval = printInteresting $-             levmar helval-                    Nothing-                    helval_params-                    helval_samples-                    1000-                    opts-                    Nothing-                    Nothing-                    noLinearConstraints-                    Nothing--run_helval_jac :: IO ()-run_helval_jac = printInteresting $-                 levmar helval-                        (Just heval_jac)-                        helval_params-                        helval_samples-                        1000-                        opts-                        Nothing-                        Nothing-                        noLinearConstraints-                        Nothing---- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--- !! TODO: This function exits with the following error: !!--- !! <interactive>: (==): No overloading for function    !!--- !! <interactive>: interrupted                          !!--- !! <interactive>: warning: too many hs_exit()s         !!--- !!                                                     !!--- !! Process haskell exited abnormally with code 252     !!--- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--run_helval_autojac :: IO ()-run_helval_autojac = printInteresting $-                     AD.levmar helval-                               helval_params-                               helval_samples-                               1000-                               opts-                               Nothing-                               Nothing-                               noLinearConstraints-                               Nothing------------------------------------------------------------------------------------- Boggs - Tolle problem 3 (linearly constrained),--- minimum at (-0.76744, 0.25581, 0.62791, -0.11628, 0.25581)------ constr1: p0 + 3*p1      = 0--- constr2: p2 + p3 - 2*p4 = 0--- constr3: p1 - p4          = 0--bt3 :: Floating r => Model N5 N5 r-bt3 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1-                                  + sqr t2-                                  + sqr t3-                                  + sqr t4-                                  )-    where-      t1 = p0 - p1-      t2 = p1 + p2 - 2.0-      t3 = p3 - 1.0-      t4 = p4 - 1.0--bt3_jac :: Floating r => Jacobian N5 N5 r-bt3_jac p0 p1 p2 p3 p4 = SL.replicate (   2.0*t1-                                      ::: 2.0*(t2 - t1)-                                      ::: 2.0*t2-                                      ::: 2.0*t3-                                      ::: 2.0*t4-                                      ::: Nil-                                      )-    where-      t1 = p0 - p1-      t2 = p1 + p2 - 2.0-      t3 = p3 - 1.0-      t4 = p4 - 1.0--bt3_params :: Floating r => SizedList N5 r-bt3_params = 2.0 ::: 2.0 ::: 2.0 :::2.0 ::: 2.0 ::: Nil--bt3_samples :: Floating r => SizedList N5 r-bt3_samples = SL.replicate 0.0--bt3_linear_constraints :: Floating r => LinearConstraints N3 N5 r-bt3_linear_constraints = (     (1.0 ::: 3.0 ::: 0.0 ::: 0.0 :::  0.0 ::: Nil)-                           ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)-                           ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)-                           ::: Nil-                         , 0.0 ::: 0.0 ::: 0.0 ::: Nil-                         )--run_bt3 :: IO ()-run_bt3 = printInteresting $-          levmar bt3-                 Nothing-                 bt3_params-                 bt3_samples-                 1000-                 opts-                 Nothing-                 Nothing-                 (Just bt3_linear_constraints)-                 Nothing--run_bt3_jac :: IO ()-run_bt3_jac = printInteresting $-              levmar bt3-                     (Just bt3_jac)-                     bt3_params-                     bt3_samples-                     1000-                     opts-                     Nothing-                     Nothing-                     (Just bt3_linear_constraints)-                     Nothing--run_bt3_autojac :: IO ()-run_bt3_autojac = printInteresting $-                  AD.levmar bt3-                            bt3_params-                            bt3_samples-                            1000-                            opts-                            Nothing-                            Nothing-                            (Just bt3_linear_constraints)-                            Nothing------------------------------------------------------------------------------------- Hock - Schittkowski problem 28 (linearly constrained),--- minimum at (0.5, -0.5, 0.5)------ constr1: p0 + 2*p1 + 3*p2 = 1--hs28 :: Floating r => Model N3 N3 r-hs28 p0 p1 p2 = SL.replicate ( sqr t1-                             + sqr t2-                             )-    where-      t1 = p0 + p1-      t2 = p1 + p2--hs28_jac :: Floating r => Jacobian N3 N3 r-hs28_jac p0 p1 p2 = SL.replicate (   2.0*t1-                                 ::: 2.0*(t1 + t2)-                                 ::: 2.0*t2-                                 ::: Nil-                                 )-    where-      t1 = p0 + p1-      t2 = p1 + p2--hs28_params :: Floating r => SizedList N3 r-hs28_params = -4.0 ::: 1.0 ::: 1.0 ::: Nil--hs28_samples :: Floating r => SizedList N3 r-hs28_samples = SL.replicate 0.0--hs28_linear_constraints :: Floating r => LinearConstraints N1 N3 r-hs28_linear_constraints = ( ((1.0 ::: 2.0 ::: 3.0 ::: Nil) ::: Nil)-                          , 1.0 ::: Nil-                          )--run_hs28 :: IO ()-run_hs28 = printInteresting $-           levmar hs28-                  Nothing-                  hs28_params-                  hs28_samples-                  1000-                  opts-                  Nothing-                  Nothing-                  (Just hs28_linear_constraints)-                  Nothing--run_hs28_jac :: IO ()-run_hs28_jac = printInteresting $-               levmar hs28-                      (Just hs28_jac)-                      hs28_params-                      hs28_samples-                      1000-                      opts-                      Nothing-                      Nothing-                      (Just hs28_linear_constraints)-                      Nothing--run_hs28_autojac :: IO ()-run_hs28_autojac = printInteresting $-                   AD.levmar hs28-                             hs28_params-                             hs28_samples-                             1000-                             opts-                             Nothing-                             Nothing-                             (Just hs28_linear_constraints)-                             Nothing------------------------------------------------------------------------------------- Hock - Schittkowski problem 48 (linearly constrained),--- minimum at (1.0, 1.0, 1.0, 1.0, 1.0)------ constr1: sum [p0, p1, p2, p3, p4] = 5--- constr2: p2 - 2*(p3 + p4)       = -3--hs48 :: Floating r => Model N5 N5 r-hs48 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1-                                   + sqr t2-                                   + sqr t3-                                   )-    where-      t1 = p0 - 1.0-      t2 = p1 - p2-      t3 = p3 - p4--hs48_jac :: Floating r => Jacobian N5 N5 r-hs48_jac p0 p1 p2 p3 p4 = SL.replicate (    2.0*t1-                                       :::  2.0*t2-                                       ::: -2.0*t2-                                       :::  2.0*t3-                                       ::: -2.0*t3-                                       ::: Nil-                                       )-    where-      t1 = p0 - 1.0-      t2 = p1 - p2-      t3 = p3 - p4--hs48_params :: Floating r => SizedList N5 r-hs48_params = 3.0 ::: 5.0 ::: -3.0 ::: 2.0 ::: -2.0 ::: Nil--hs48_samples :: Floating r => SizedList N5 r-hs48_samples = SL.replicate 0.0--hs48_linear_constraints :: Floating r => LinearConstraints N2 N5 r-hs48_linear_constraints = (     (1.0 ::: 1.0 ::: 1.0 :::  1.0 :::  1.0 ::: Nil)-                            ::: (0.0 ::: 0.0 ::: 1.0 ::: -2.0 ::: -2.0 ::: Nil)-                            ::: Nil-                          , 5.0 ::: -3.0 ::: Nil-                          )--run_hs48 :: IO ()-run_hs48 = printInteresting $-           levmar hs48-                  Nothing-                  hs48_params-                  hs48_samples-                  1000-                  opts-                  Nothing-                  Nothing-                  (Just hs48_linear_constraints)-                  Nothing--run_hs48_jac :: IO ()-run_hs48_jac = printInteresting $-               levmar hs48-                      (Just hs48_jac)-                      hs48_params-                      hs48_samples-                      1000-                      opts-                      Nothing-                      Nothing-                      (Just hs48_linear_constraints)-                      Nothing--run_hs48_autojac :: IO ()-run_hs48_autojac = printInteresting $-                   AD.levmar hs48-                             hs48_params-                             hs48_samples-                             1000-                             opts-                             Nothing-                             Nothing-                             (Just hs48_linear_constraints)-                             Nothing------------------------------------------------------------------------------------- Hock - Schittkowski problem 51 (linearly constrained),--- minimum at (1.0, 1.0, 1.0, 1.0, 1.0)------ constr1: p0 + 3*p1      = 4--- constr2: p2 + p3 - 2*p4 = 0--- constr3: p1 - p4          = 0--hs51 :: Floating r => Model N5 N5 r-hs51 p0 p1 p2 p3 p4 = SL.replicate ( sqr t1-                                   + sqr t2-                                   + sqr t3-                                   + sqr t4-                                   )-    where-      t1 = p0 - p1-      t2 = p1 + p2 - 2.0-      t3 = p3 - 1.0-      t4 = p4 - 1.0--hs51_jac :: Floating r => Jacobian N5 N5 r-hs51_jac p0 p1 p2 p3 p4 = SL.replicate (   2.0*t1-                                       ::: 2.0*(t2 - t1)-                                       ::: 2.0*t2-                                       ::: 2.0*t3-                                       ::: 2.0*t4-                                       ::: Nil-                                       )-    where-      t1 = p0 - p1-      t2 = p1 + p2 - 2.0-      t3 = p3 - 1.0-      t4 = p4 - 1.0--hs51_params :: Floating r => SizedList N5 r-hs51_params = 2.5 ::: 0.5 ::: 2.0 ::: -1.0 ::: 0.5 ::: Nil--hs51_samples :: Floating r => SizedList N5 r-hs51_samples = SL.replicate 0.0--hs51_linear_constraints :: Floating r => LinearConstraints N3 N5 r-hs51_linear_constraints = (     (1.0 ::: 3.0 ::: 0.0 ::: 0.0 :::  0.0 ::: Nil)-                            ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)-                            ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)-                            ::: Nil-                          , 4.0 ::: 0.0 ::: 0.0 ::: Nil-                          )--run_hs51 :: IO ()-run_hs51 = printInteresting $-           levmar hs51-                  Nothing-                  hs51_params-                  hs51_samples-                  1000-                  opts-                  Nothing-                  Nothing-                  (Just hs51_linear_constraints)-                  Nothing--run_hs51_jac :: IO ()-run_hs51_jac = printInteresting $-               levmar hs51-                      (Just hs51_jac)-                      hs51_params-                      hs51_samples-                      1000-                      opts-                      Nothing-                      Nothing-                      (Just hs51_linear_constraints)-                      Nothing--run_hs51_autojac :: IO ()-run_hs51_autojac = printInteresting $-                   AD.levmar hs51-                             hs51_params-                             hs51_samples-                             1000-                             opts-                             Nothing-                             Nothing-                             (Just hs51_linear_constraints)-                             Nothing------------------------------------------------------------------------------------- Hock - Schittkowski problem 01 (box constrained),--- minimum at (1.0, 1.0)------ constr1: p1 >= -1.5--hs01 :: Floating r => Model N2 N2 r-hs01 p0 p1 =     10.0*(p1 - sqr p0)-             ::: 1.0 - p0-             ::: Nil--hs01_jac :: Floating r => Jacobian N2 N2 r-hs01_jac p0 _ =     (-20.0*p0 ::: 10.0 ::: Nil)-                ::: (-1.0     ::: 0.0  ::: Nil)-                ::: Nil--hs01_params :: Floating r => SizedList N2 r-hs01_params = -2.0 ::: 1.0 ::: Nil--hs01_samples :: Floating r => SizedList N2 r-hs01_samples = SL.replicate 0.0--hs01_lb, hs01_ub :: Floating r => SizedList N2 r-hs01_lb = -_DBL_MAX ::: -1.5     ::: Nil-hs01_ub =  _DBL_MAX ::: _DBL_MAX ::: Nil--_DBL_MAX :: Floating r => r-_DBL_MAX = 1e+37 -- TODO: Get this directly from <float.h>.--run_hs01 :: IO ()-run_hs01 = printInteresting $-           levmar hs01-                  Nothing-                  hs01_params-                  hs01_samples-                  1000-                  opts-                  (Just hs01_lb)-                  (Just hs01_ub)-                  noLinearConstraints-                  Nothing--run_hs01_jac :: IO ()-run_hs01_jac = printInteresting $-               levmar hs01-                      (Just hs01_jac)-                      hs01_params-                      hs01_samples-                      1000-                      opts-                      (Just hs01_lb)-                      (Just hs01_ub)-                      noLinearConstraints-                      Nothing--run_hs01_autojac :: IO ()-run_hs01_autojac = printInteresting $-                   AD.levmar hs01-                             hs01_params-                             hs01_samples-                             1000-                             opts-                             (Just hs01_lb)-                             (Just hs01_ub)-                             noLinearConstraints-                             Nothing------------------------------------------------------------------------------------- Hock - Schittkowski MODIFIED problem 21 (box constrained),--- minimum at (2.0, 0.0)------ constr1: 2 <= p0 <=50--- constr2: -50 <= p1 <=50------ Original HS21 has the additional constraint 10*p0 - p1 >= 10--- which is inactive at the solution, so it is dropped here.--hs21 :: Floating r => Model N2 N2 r-hs21 p0 p1 =     p0 / 10.0-             ::: p1-             ::: Nil--hs21_jac :: Floating r => Jacobian N2 N2 r-hs21_jac _ _ =     (0.1 ::: 0.0 ::: Nil)-               ::: (0.0 ::: 1.0 ::: Nil)-               ::: Nil--hs21_params :: Floating r => SizedList N2 r-hs21_params = -1.0 ::: -1.0 ::: Nil--hs21_samples :: Floating r => SizedList N2 r-hs21_samples = SL.replicate 0.0--hs21_lb, hs21_ub :: Floating r => SizedList N2 r-hs21_lb = 2.0  ::: -50.0 ::: Nil-hs21_ub = 50.0 :::  50.0 ::: Nil--run_hs21 :: IO ()-run_hs21 = printInteresting $-           levmar hs21-                  Nothing-                  hs21_params-                  hs21_samples-                  1000-                  opts-                  (Just hs21_lb)-                  (Just hs21_ub)-                  noLinearConstraints-                  Nothing--run_hs21_jac :: IO ()-run_hs21_jac = printInteresting $-               levmar hs21-                      (Just hs21_jac)-                      hs21_params-                      hs21_samples-                      1000-                      opts-                      (Just hs21_lb)-                      (Just hs21_ub)-                      noLinearConstraints-                      Nothing--run_hs21_autojac :: IO ()-run_hs21_autojac = printInteresting $-                   AD.levmar hs21-                             hs21_params-                             hs21_samples-                             1000-                             opts-                             (Just hs21_lb)-                             (Just hs21_ub)-                             noLinearConstraints-                             Nothing------------------------------------------------------------------------------------- Problem hatfldb (box constrained),--- minimum at (0.947214, 0.8, 0.64, 0.4096)------ constri: pi >= 0.0 (i=1..4)--- constr5: p1 <= 0.8--hatfldb :: Floating r => Model N4 N4 r-hatfldb p0 p1 p2 p3 =     p0 - 1.0-                      ::: p0 - sqrt p1-                      ::: p1 - sqrt p2-                      ::: p2 - sqrt p3-                      ::: Nil--hatfldb_jac :: Floating r => Jacobian N4 N4 r-hatfldb_jac _ p1 p2 p3 =     (1.0 ::: 0.0            ::: 0.0            ::: 0.0            ::: Nil)-                         ::: (1.0 ::: -0.5 / sqrt p1 ::: 0.0            ::: 0.0            ::: Nil)-                         ::: (0.0 ::: 1.0            ::: -0.5 / sqrt p2 ::: 0.0            ::: Nil)-                         ::: (0.0 ::: 0.0            ::: 1.0            ::: -0.5 / sqrt p3 ::: Nil)-                         ::: Nil--hatfldb_params :: Floating r => SizedList N4 r-hatfldb_params = 0.1 ::: 0.1 ::: 0.1 ::: 0.1 ::: Nil--hatfldb_samples :: Floating r => SizedList N4 r-hatfldb_samples = SL.replicate 0.0--hatfldb_lb, hatfldb_ub :: Floating r => SizedList N4 r-hatfldb_lb = 0.0      ::: 0.0 ::: 0.0      ::: 0.0      ::: Nil-hatfldb_ub = _DBL_MAX ::: 0.8 ::: _DBL_MAX ::: _DBL_MAX ::: Nil--run_hatfldb :: IO ()-run_hatfldb = printInteresting $-              levmar hatfldb-                     Nothing-                     hatfldb_params-                     hatfldb_samples-                     1000-                     opts-                     (Just hatfldb_lb)-                     (Just hatfldb_ub)-                     noLinearConstraints-                     Nothing--run_hatfldb_jac :: IO ()-run_hatfldb_jac = printInteresting $-                  levmar hatfldb-                         (Just hatfldb_jac)-                         hatfldb_params-                         hatfldb_samples-                         1000-                         opts-                         (Just hatfldb_lb)-                         (Just hatfldb_ub)-                         noLinearConstraints-                         Nothing--run_hatfldb_autojac :: IO ()-run_hatfldb_autojac = printInteresting $-                      AD.levmar hatfldb-                                hatfldb_params-                                hatfldb_samples-                                1000-                                opts-                                (Just hatfldb_lb)-                                (Just hatfldb_ub)-                                noLinearConstraints-                                Nothing------------------------------------------------------------------------------------- Problem hatfldc (box constrained),--- minimum at (1.0, 1.0, 1.0, 1.0)------ constri:   pi >= 0.0  (i=1..4)--- constri+4: pi <= 10.0 (i=1..4)--hatfldc :: Floating r => Model N4 N4 r-hatfldc p0 p1 p2 p3 =     p0 - 1.0-                      ::: p0 - sqrt p1-                      ::: p1 - sqrt p2-                      ::: p3 - 1.0-                      ::: Nil--hatfldc_jac :: Floating r => Jacobian N4 N4 r-hatfldc_jac _ p1 p2 _ =     (1.0 ::: 0.0            ::: 0.0            ::: 0.0 ::: Nil)-                        ::: (1.0 ::: -0.5 / sqrt p1 ::: 0.0            ::: 0.0 ::: Nil)-                        ::: (0.0 ::: 1.0            ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil)-                        ::: (0.0 ::: 0.0            ::: 0.0            ::: 1.0 ::: Nil)-                        ::: Nil--hatfldc_params :: Floating r => SizedList N4 r-hatfldc_params = 0.9 ::: 0.9 ::: 0.9 ::: 0.9 ::: Nil--hatfldc_samples :: Floating r => SizedList N4 r-hatfldc_samples = SL.replicate 0.0--hatfldc_lb, hatfldc_ub :: Floating r => SizedList N4 r-hatfldc_lb =  0.0 :::  0.0 :::  0.0 :::  0.0 ::: Nil-hatfldc_ub = 10.0 ::: 10.0 ::: 10.0 ::: 10.0 ::: Nil--run_hatfldc :: IO ()-run_hatfldc = printInteresting $-              levmar hatfldc-                     Nothing-                     hatfldc_params-                     hatfldc_samples-                     1000-                     opts-                     (Just hatfldc_lb)-                     (Just hatfldc_ub)-                     noLinearConstraints-                     Nothing--run_hatfldc_jac :: IO ()-run_hatfldc_jac = printInteresting $-                  levmar hatfldc-                         (Just hatfldc_jac)-                         hatfldc_params-                         hatfldc_samples-                         1000-                         opts-                         (Just hatfldc_lb)-                         (Just hatfldc_ub)-                         noLinearConstraints-                         Nothing--run_hatfldc_autojac :: IO ()-run_hatfldc_autojac = printInteresting $-                      AD.levmar hatfldc-                                hatfldc_params-                                hatfldc_samples-                                1000-                                opts-                                (Just hatfldc_lb)-                                (Just hatfldc_ub)-                                noLinearConstraints-                                Nothing------------------------------------------------------------------------------------- Hock - Schittkowski (modified) problem 52 (box/linearly constrained),--- minimum at (-0.09, 0.03, 0.25, -0.19, 0.03)------ constr1: p0 + 3*p1 = 0--- constr2: p2 +   p3 - 2*p4 = 0--- constr3: p1 -   p4 = 0------ To the above 3 constraints, we add the following 5:--- constr4: -0.09 <= p0--- constr5:   0.0 <= p1 <= 0.3--- constr6:          p2 <= 0.25--- constr7:  -0.2 <= p3 <= 0.3--- constr8:   0.0 <= p4 <= 0.3--modhs52 :: Floating r => Model N5 N4 r-modhs52 p0 p1 p2 p3 p4 =     4.0*p0 - p1-                         ::: p1 + p2 - 2.0-                         ::: p3 - 1.0-                         ::: p4 - 1.0-                         ::: Nil--modhs52_jac :: Floating r => Jacobian N5 N4 r-modhs52_jac _ _ _ _ _ =     (4.0 ::: -1.0 ::: 0.0 ::: 0.0 ::: 0.0 ::: Nil)-                        ::: (0.0 :::  1.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: Nil)-                        ::: (0.0 :::  0.0 ::: 0.0 ::: 1.0 ::: 0.0 ::: Nil)-                        ::: (0.0 :::  0.0 ::: 0.0 ::: 0.0 ::: 1.0 ::: Nil)-                        ::: Nil--modhs52_params :: Floating r => SizedList N5 r-modhs52_params = 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: 2.0 ::: Nil--modhs52_samples :: Floating r => SizedList N4 r-modhs52_samples = SL.replicate 0.0--modhs52_linear_constraints :: Floating r => LinearConstraints N3 N5 r-modhs52_linear_constraints = (     (1.0 ::: 3.0 ::: 0.0 ::: 0.0 :::  0.0 ::: Nil)-                               ::: (0.0 ::: 0.0 ::: 1.0 ::: 1.0 ::: -2.0 ::: Nil)-                               ::: (0.0 ::: 1.0 ::: 0.0 ::: 0.0 ::: -1.0 ::: Nil)-                               ::: Nil-                             , 0.0 ::: 0.0 ::: 0.0 ::: Nil-                             )--modhs52_weights :: Floating r => SizedList N5 r-modhs52_weights = 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: 2000.0 ::: Nil--modhs52_lb, modhs52_ub :: Floating r => SizedList N5 r-modhs52_lb = -0.09    ::: 0.0 ::: -_DBL_MAX ::: -0.2 ::: 0.0 ::: Nil-modhs52_ub = _DBL_MAX ::: 0.3 ::: 0.25      :::  0.3 ::: 0.3 ::: Nil--run_modhs52 :: IO ()-run_modhs52 = printInteresting $-              levmar modhs52-                     Nothing-                     modhs52_params-                     modhs52_samples-                     1000-                     opts-                     (Just modhs52_lb)-                     (Just modhs52_ub)-                     (Just modhs52_linear_constraints)-                     (Just modhs52_weights)--run_modhs52_jac :: IO ()-run_modhs52_jac = printInteresting $-                  levmar modhs52-                         (Just modhs52_jac)-                         modhs52_params-                         modhs52_samples-                         1000-                         opts-                         (Just modhs52_lb)-                         (Just modhs52_ub)-                         (Just modhs52_linear_constraints)-                         (Just modhs52_weights)--run_modhs52_autojac :: IO ()-run_modhs52_autojac = printInteresting $-                      AD.levmar modhs52-                                modhs52_params-                                modhs52_samples-                                1000-                                opts-                                (Just modhs52_lb)-                                (Just modhs52_ub)-                                (Just modhs52_linear_constraints)-                                (Just modhs52_weights)------------------------------------------------------------------------------------- Schittkowski (modified) problem 235 (box/linearly constrained),--- minimum at (-1.725, 2.9, 0.725)------ constr1: p0 + p2 = -1.0;------ To the above constraint, we add the following 2:--- constr2: p1 - 4*p2 = 0--- constr3: 0.1 <= p1 <= 2.9--- constr4: 0.7 <= p2--mods235 :: Floating r => Model N3 N2 r-mods235 p0 p1 _ =     0.1*(p0 - 1.0)-                  ::: p1 - sqr p0-                  ::: Nil--mods235_jac :: Floating r => Jacobian N3 N2 r-mods235_jac p0 _ _ =     (0.1     ::: 0.0 ::: 0.0 ::: Nil)-                     ::: (-2.0*p0 ::: 1.0 ::: 0.0 ::: Nil)-                     ::: Nil--mods235_params :: Floating r => SizedList N3 r-mods235_params = -2.0 ::: 3.0 ::: 1.0 ::: Nil--mods235_samples :: Floating r => SizedList N2 r-mods235_samples = SL.replicate 0.0--mods235_linear_constraints :: Floating r => LinearConstraints N2 N3 r-mods235_linear_constraints = (     (1.0 ::: 0.0 :::  1.0 ::: Nil)-                               ::: (0.0 ::: 1.0 ::: -4.0 ::: Nil)-                               ::: Nil-                             , -1.0 ::: 0.0 ::: Nil-                             )--mods235_lb, mods235_ub :: Floating r => SizedList N3 r-mods235_lb = -_DBL_MAX ::: 0.1 ::: 0.7      ::: Nil-mods235_ub =  _DBL_MAX ::: 2.9 ::: _DBL_MAX ::: Nil--run_mods235 :: IO ()-run_mods235 = printInteresting $-              levmar mods235-                     Nothing-                     mods235_params-                     mods235_samples-                     1000-                     opts-                     (Just mods235_lb)-                     (Just mods235_ub)-                     (Just mods235_linear_constraints)-                     Nothing--run_mods235_jac :: IO ()-run_mods235_jac = printInteresting $-                  levmar mods235-                         (Just mods235_jac)-                         mods235_params-                         mods235_samples-                         1000-                         opts-                         (Just mods235_lb)-                         (Just mods235_ub)-                         (Just mods235_linear_constraints)-                         Nothing---run_mods235_autojac :: IO ()-run_mods235_autojac = printInteresting $-                      AD.levmar mods235-                             mods235_params-                             mods235_samples-                             1000-                             opts-                             (Just mods235_lb)-                             (Just mods235_ub)-                             (Just mods235_linear_constraints)-                             Nothing------------------------------------------------------------------------------------- Boggs and Tolle modified problem 7 (box/linearly constrained),--- minimum at (0.7, 0.49, 0.19, 1.19, -0.2)------ We keep the original objective function & starting point and use the--- following constraints:------ subject to cons1:---  x[1]+x[2] - x[3] = 1.0;--- subject to cons2:---   x[2] - x[4] + x[1] = 0.0;--- subject to cons3:---   x[5] + x[1] = 0.5;--- subject to cons4:---   x[5]>=-0.3;--- subject to cons5:---    x[1]<=0.7;--modbt7 :: Floating r => Model N5 N5 r-modbt7 p0 p1 _ _ _ = SL.replicate (100.0*sqr m + sqr n)-    where-      m = p1 - sqr p0-      n = p0 - 1.0--modbt7_jac :: Floating r => Jacobian N5 N5 r-modbt7_jac p0 p1 _ _ _ = SL.replicate-                         (    -400.0*m*p0 + 2.0*p0 - 2.0-                           ::: 200.0*m-                           ::: 0.0-                           ::: 0.0-                           ::: 0.0-                           ::: Nil-                         )-    where-      m = p1 - sqr p0--modbt7_params :: Floating r => SizedList N5 r-modbt7_params = -2.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: 1.0 ::: Nil--modbt7_samples :: Floating r => SizedList N5 r-modbt7_samples = SL.replicate 0.0--modbt7_linear_constraints :: Floating r => LinearConstraints N3 N5 r-modbt7_linear_constraints = (     (1.0 ::: 1.0 ::: -1.0 :::  0.0 ::: 0.0 ::: Nil)-                              ::: (1.0 ::: 1.0 :::  0.0 ::: -1.0 ::: 0.0 ::: Nil)-                              ::: (1.0 ::: 0.0 :::  0.0 :::  0.0 ::: 1.0 ::: Nil)-                              ::: Nil-                            , 1.0 ::: 0.0 ::: 0.5 ::: Nil-                            )--modbt7_lb, modbt7_ub :: Floating r => SizedList N5 r-modbt7_lb = -_DBL_MAX ::: -_DBL_MAX ::: -_DBL_MAX ::: -_DBL_MAX ::: -0.3     ::: Nil-modbt7_ub = 0.7       ::: _DBL_MAX  ::: _DBL_MAX  ::: _DBL_MAX  ::: _DBL_MAX ::: Nil---- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--- !! TODO: Find out why these return with: infStopReason = MaxIterations !!--- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!--run_modbt7 :: IO ()-run_modbt7 = printInteresting $-             levmar modbt7-                     Nothing-                     modbt7_params-                     modbt7_samples-                     1000-                     opts-                     (Just modbt7_lb)-                     (Just modbt7_ub)-                     (Just modbt7_linear_constraints)-                     Nothing--run_modbt7_jac :: IO ()-run_modbt7_jac = printInteresting $-                 levmar modbt7-                        (Just modbt7_jac)-                        modbt7_params-                        modbt7_samples-                        1000-                        opts-                        (Just modbt7_lb)-                        (Just modbt7_ub)-                        (Just modbt7_linear_constraints)-                        Nothing--run_modbt7_autojac :: IO ()-run_modbt7_autojac = printInteresting $-                     AD.levmar modbt7-                               modbt7_params-                               modbt7_samples-                               1000-                               opts-                               (Just modbt7_lb)-                               (Just modbt7_ub)-                               (Just modbt7_linear_constraints)-                               Nothing------------------------------------------------------------------------------------- Equilibrium combustion problem, constrained nonlinear equation from the book--- by Floudas et al.------ Minimum at (0.0034, 31.3265, 0.0684, 0.8595, 0.0370)------ constri:   pi>=0.0001 (i=1..5)--- constri+5: pi<=100.0  (i=1..5)--combust :: Floating r => Model N5 N5 r-combust p0 p1 p2 p3 p4 =-        p0*p1 + p0 - 3*p4-    ::: 2*p0*p1 + p0 + 3*r10*p1*p1 + p1*p2*p2 + r7*p1*p2 + r9*p1*p3 + r8*p1 - r*p4-    ::: 2*p1*p2*p2 + r7*p1*p2 + 2*r5*p2*p2 + r6*p2-8*p4-    ::: r9*p1*p3 + 2*p3*p3 - 4*r*p4-    ::: p0*p1 + p0 + r10*p1*p1 + p1*p2*p2 + r7*p1*p2 + r9*p1*p3 + r8*p1 + r5*p2*p2 + r6*p2 + p3*p3 - 1.0-    ::: Nil--r, r5, r6, r7, r8, r9, r10 :: Floating r => r-r   = 10-r5  = 0.193-r6  = 4.10622*1e-4-r7  = 5.45177*1e-4-r8  = 4.4975 *1e-7-r9  = 3.40735*1e-5-r10 = 9.615  *1e-7--combust_jac :: Floating r => Jacobian N5 N5 r-combust_jac p0 p1 p2 p3 _ =-        (   p1 + 1-        ::: p0-        ::: 0.0-        ::: 0.0-        ::: -3-        ::: Nil-        )-    ::: (   2*p1 + 1-        ::: 2*p0 + 6*r10*p1 + p2*p2 + r7*p2 + r9*p3 + r8-        ::: 2*p1*p2 + r7*p1-        ::: r9*p1-        ::: -r-        ::: Nil-        )-    ::: (   0.0-        ::: 2*p2*p2 + r7*p2-        ::: 4*p1*p2 + r7*p1 + 4*r5*p2 + r6-        ::: 0.0-        ::: -8-        ::: Nil-        )-    ::: (   0.0-        ::: r9*p3-        ::: 0.0-        ::: r9*p1 + 4*p3-        ::: -4*r-        ::: Nil-        )-    ::: (   p1 + 1-        ::: p0 + 2*r10*p1 + p2*p2 + r7*p2 + r9*p3 + r8-        ::: 2*p1*p2 + r7*p1 + 2*r5*p2 + r6-        ::: r9*p1 + 2*p3-        ::: 0.0-        ::: Nil-        )-    ::: Nil--combust_params :: Floating r => SizedList N5 r-combust_params = 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: 0.0001 ::: Nil--combust_samples :: Floating r => SizedList N5 r-combust_samples = SL.replicate 0.0--combust_lb, combust_ub :: Floating r => SizedList N5 r-combust_lb =   0.0001 :::   0.0001 :::   0.0001 :::   0.0001 :::   0.0001 ::: Nil-combust_ub = 100.0    ::: 100.0    ::: 100.0    ::: 100.0    ::: 100.0    ::: Nil--run_combust :: IO ()-run_combust = printInteresting $-              levmar combust-                     Nothing-                     combust_params-                     combust_samples-                     1000-                     opts-                     (Just combust_lb)-                     (Just combust_ub)-                     noLinearConstraints-                     Nothing--run_combust_jac :: IO ()-run_combust_jac = printInteresting $-                  levmar combust-                         (Just combust_jac)-                         combust_params-                         combust_samples-                         1000-                         opts-                         (Just combust_lb)-                         (Just combust_ub)-                         noLinearConstraints-                         Nothing--run_combust_autojac :: IO ()-run_combust_autojac = printInteresting $-                      AD.levmar combust-                                combust_params-                                combust_samples-                                1000-                                opts-                                (Just combust_lb)-                                (Just combust_ub)-                                noLinearConstraints-                                Nothing----- The End ---------------------------------------------------------------------

@@ -1,150 +0,0 @@-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE ScopedTypeVariables #-}------------------------------------------------------------------------------------- |--- Module      :  LevMar--- Copyright   :  (c) 2009 Roel van Dijk & Bas van Dijk--- License     :  BSD-style (see the file LICENSE)------ Maintainer  :  vandijk.roel@gmail.com, v.dijk.bas@gmail.com--- Stability   :  Experimental------------ For additional documentation see the documentation of the levmar C--- library which this library is based on:--- <http://www.ics.forth.gr/~lourakis/levmar/>--------------------------------------------------------------------------------------module LevMar-    ( -- * Model & Jacobian.-      Model-    , Jacobian--      -- * Levenberg-Marquardt algorithm.-    , LMA_I.LevMarable-    , levmar--    , LinearConstraints-    , noLinearConstraints-    , Matrix--      -- * Minimization options.-    , LMA_I.Options(..)-    , LMA_I.defaultOpts--      -- * Output-    , LMA_I.Info(..)-    , LMA_I.StopReason(..)-    , CovarMatrix--    , LMA_I.LevMarError(..)--      -- *Type-level machinery-    , Z, S, Nat-    , SizedList(..)-    , NFunction-    )-    where---import qualified LevMar.Intermediate as LMA_I--import LevMar.Utils ( LinearConstraints-                    , noLinearConstraints-                    , Matrix-                    , CovarMatrix-                    , convertLinearConstraints-                    , convertResult-                    )--import TypeLevelNat ( Z, S, Nat )-import SizedList    ( SizedList(..), toList, unsafeFromList )-import NFunction    ( NFunction, ($*) )--import Data.Either-------------------------------------------------------------------------------------- Model & Jacobian.-----------------------------------------------------------------------------------{- | A functional relation describing measurements represented as a function-from @m@ parameters to @n@ expected measurements.--An example from /Demo.hs/:--@-type N4 = 'S' ('S' ('S' ('S' 'Z')))--hatfldc :: Model N4 N4 Double-hatfldc p0 p1 p2 p3 =     p0 - 1.0-                      ::: p0 - sqrt p1-                      ::: p1 - sqrt p2-                      ::: p3 - 1.0-                      ::: Nil-@--}-type Model m n r = NFunction m r (SizedList n r)--{- | The jacobian of the 'Model' function. Expressed as a function-from @m@ parameters to a @n@/x/@m@ matrix which for each of the @n@-expected measurement describes the @m@ partial derivatives of the-parameters.--See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>--For example the jacobian of the above @hatfldc@ model is:--@-type N4 = 'S' ('S' ('S' ('S' 'Z')))--hatfldc_jac :: Jacobian N4 N4 Double-hatfldc_jac _ p1 p2 _ =     (1.0 ::: 0.0            ::: 0.0            ::: 0.0 ::: Nil)-                        ::: (1.0 ::: -0.5 / sqrt p1 ::: 0.0            ::: 0.0 ::: Nil)-                        ::: (0.0 ::: 1.0            ::: -0.5 / sqrt p2 ::: 0.0 ::: Nil)-                        ::: (0.0 ::: 0.0            ::: 0.0            ::: 1.0 ::: Nil)-                        ::: Nil-@--}--type Jacobian m n r = NFunction m r (Matrix n m r)-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm.-levmar :: forall m n k r. (Nat m, Nat n, Nat k, LMA_I.LevMarable r)-       => (Model m n r)                   -- ^ Model-       -> Maybe (Jacobian m n r)          -- ^ Optional jacobian-       -> SizedList m r                   -- ^ Initial parameters-       -> SizedList n r                   -- ^ Samples-       -> Integer                         -- ^ Maximum number of iterations-       -> LMA_I.Options r                 -- ^ Minimization options-       -> Maybe (SizedList m r)           -- ^ Optional lower bounds-       -> Maybe (SizedList m r)           -- ^ Optional upper bounds-       -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints-       -> Maybe (SizedList m r)           -- ^ Optional weights-       -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r)--levmar model mJac params ys itMax opts mLowBs mUpBs mLinC mWghts =-    fmap convertResult $ LMA_I.levmar (convertModel model)-                                      (fmap convertJacob mJac)-                                      (toList params)-                                      (toList ys)-                                      itMax-                                      opts-                                      (fmap toList mLowBs)-                                      (fmap toList mUpBs)-                                      (fmap convertLinearConstraints mLinC)-                                      (fmap toList mWghts)-    where-      convertModel f = \ps -> toList (f $* (unsafeFromList ps :: SizedList m r) :: SizedList n r)-      convertJacob f = \ps -> toList (fmap toList (f $* (unsafeFromList ps :: SizedList m r) :: Matrix n m r))----- The End ---------------------------------------------------------------------

@@ -1,142 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE FlexibleContexts #-}------------------------------------------------------------------------------------- |--- Module      :  LevMar.AD--- Copyright   :  (c) 2009 Roel van Dijk & Bas van Dijk--- License     :  BSD-style (see the file LICENSE)------ Maintainer  :  vandijk.roel@gmail.com, v.dijk.bas@gmail.com--- Stability   :  Experimental------ A levmar variant that uses Automatic Differentiation to--- automatically compute the Jacobian.------ For additional documentation see the documentation of the levmar C--- library which this library is based on:--- <http://www.ics.forth.gr/~lourakis/levmar/>---------------------------------------------------------------------------------------module LevMar.AD-    ( -- * Model-      LMA.Model--      -- * Levenberg-Marquardt algorithm.-    , LMA_I.LevMarable-    , levmar--    , LinearConstraints-    , noLinearConstraints-    , Matrix--      -- * Minimization options.-    , LMA_I.Options(..)-    , LMA_I.defaultOpts--      -- * Output-    , LMA_I.Info(..)-    , LMA_I.StopReason(..)-    , CovarMatrix--    , LMA_I.LevMarError(..)--      -- *Type-level machinery-    , Z, S, Nat-    , SizedList(..)-    , NFunction-    )-    where---import qualified LevMar              as LMA-import qualified LevMar.Intermediate as LMA_I--import LevMar.Utils ( LinearConstraints-                    , noLinearConstraints-                    , Matrix-                    , CovarMatrix-                    , convertLinearConstraints-                    , convertResult-                    )--import TypeLevelNat ( Z, S, Nat )-import SizedList    ( SizedList(..), toList, unsafeFromList )-import NFunction    ( NFunction, ($*) )--import LevMar.Utils.AD  ( value, firstDeriv, constant, idDAt )---- From vector-space:-import Data.Derivative  ( (:~>) )-import Data.VectorSpace ( VectorSpace, Scalar )-import Data.Basis       ( HasBasis, Basis )--import Data.List        ( transpose )-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm that automatically computes the--- 'Jacobian' using automatic differentiation of the model function.------ /Warning/: Don't apply 'levmar' to 'LMA.Model's that apply methods of--- the 'Eq' and 'Ord' classes to the parameters. These methods are--- undefined for ':~>'!!!-levmar :: forall m n k r.-          ( Nat m-          , Nat n-          , Nat k-          , HasBasis r-          , Basis r ~ ()-          , VectorSpace (Scalar r)-          , LMA_I.LevMarable r-          )-       => (LMA.Model m n (r :~> r))       -- ^ Model. Note that ':~>'-                                          --   is overloaded for all the-                                          --   numeric classes.-       -> SizedList m r                   -- ^ Initial parameters-       -> SizedList n r                   -- ^ Samples-       -> Integer                         -- ^ Maximum number of iterations-       -> LMA_I.Options r                 -- ^ Minimization options-       -> Maybe (SizedList m r)           -- ^ Optional lower bounds-       -> Maybe (SizedList m r)           -- ^ Optional upper bounds-       -> Maybe (LinearConstraints k m r) -- ^ Optional linear constraints-       -> Maybe (SizedList m r)           -- ^ Optional weights-       -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r)--levmar model params ys itMax opts mLowBs mUpBs mLinC mWghts =-    fmap convertResult $ LMA_I.levmar (convertModel model)-                                      (Just $ jacobianOf model)-                                      (toList params)-                                      (toList ys)-                                      itMax-                                      opts-                                      (fmap toList mLowBs)-                                      (fmap toList mUpBs)-                                      (fmap convertLinearConstraints mLinC)-                                      (fmap toList mWghts)-    where-      convertModel :: LMA.Model m n (r :~> r) -> LMA_I.Model r-      (convertModel mdl) ps = fmap value $ toList-                              (mdl $* pDs :: SizedList n (r :~> r))-          where-            pDs :: SizedList m (r :~> r)-            pDs = unsafeFromList $ fmap constant ps--      jacobianOf :: LMA.Model m n (r :~> r) -> LMA_I.Jacobian r-      (jacobianOf mdl) ps = fmap (\fs -> zipWith (firstDeriv .) fs ps)-                          . transpose-                          . fmap (\pD -> toList (mdl $* (pD :: SizedList m (r :~> r)) :: SizedList n (r :~> r)))-                          $ pDs-          where-            pDs :: [SizedList m (r :~> r)]-            pDs = [unsafeFromList $ idDAt n ps | n <- [0 .. length ps - 1]]----- The End ---------------------------------------------------------------------

@@ -1,151 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}------------------------------------------------------------------------------------- |--- Module      :  LevMar.Fitting--- Copyright   :  (c) 2009 Roel van Dijk & Bas van Dijk--- License     :  BSD-style (see the file LICENSE)------ Maintainer  :  vandijk.roel@gmail.com, v.dijk.bas@gmail.com--- Stability   :  Experimental------ This module provides the Levenberg-Marquardt algorithm specialised--- for curve-fitting.------ For additional documentation see the documentation of the levmar C--- library which this library is based on:--- <http://www.ics.forth.gr/~lourakis/levmar/>--------------------------------------------------------------------------------------module LevMar.Fitting-    ( -- * Model & Jacobian.-      Model-    , SimpleModel-    , Jacobian-    , SimpleJacobian--      -- * Levenberg-Marquardt algorithm.-    , LMA_I.LevMarable-    , levmar--    , LinearConstraints-    , noLinearConstraints-    , Matrix--    -- * Minimization options.-    , LMA_I.Options(..)-    , LMA_I.defaultOpts--      -- * Output-    , LMA_I.Info(..)-    , LMA_I.StopReason(..)-    , CovarMatrix--    , LMA_I.LevMarError(..)--      -- *Type-level machinery-    , Z, S, Nat-    , SizedList(..)-    , NFunction-    ) where---import qualified LevMar.Intermediate.Fitting as LMA_I-import LevMar.Utils ( LinearConstraints-                    , noLinearConstraints-                    , convertLinearConstraints-                    , Matrix-                    , CovarMatrix-                    , convertResult-                    )--import TypeLevelNat ( Z, S, Nat )-import SizedList    ( SizedList(..), toList, unsafeFromList )-import NFunction    ( NFunction, ($*) )-------------------------------------------------------------------------------------- Model & Jacobian.-----------------------------------------------------------------------------------{- | A functional relation describing measurements represented as a function-from @m@ parameters and an x-value to an expected measurement.--For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be-written as:--@-type N3 = 'S' ('S' ('S' 'Z'))--quad :: 'Num' r => 'Model' N3 r r-quad a b c x = a*x^2 + b*x + c-@--}-type Model m r a = NFunction m r (a -> r)---- | This type synonym expresses that usually the @a@ in @'Model' m r a@--- equals the type of the parameters.-type SimpleModel m r = Model m r r--{- | The jacobian of the 'Model' function. Expressed as a function from @n@-parameters and an x-value to the @m@ partial derivatives of the parameters.--See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>--For example, the jacobian of the quadratic function @f(x) = a*x^2 +-b*x + c@ can be written as:--@-type N3 = 'S' ('S' ('S' 'Z'))--quadJacob :: 'Num' r => 'Jacobian' N3 r r-quadJacob _ _ _ x =   x^2   -- with respect to a-                  ::: x     -- with respect to b-                  ::: 1     -- with respect to c-                  ::: 'Nil'-@--Notice you don't have to differentiate for @x@.--}-type Jacobian m r a = NFunction m r (a -> SizedList m r)---- | This type synonym expresses that usually the @a@ in @'Jacobian' m r a@--- equals the type of the parameters.-type SimpleJacobian m r = Jacobian m r r-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm specialised for curve-fitting.-levmar :: forall m k r a. (Nat m, Nat k, LMA_I.LevMarable r)-       => (Model m r a)                          -- ^ Model-       -> Maybe (Jacobian m r a)                 -- ^ Optional jacobian-       -> SizedList m r                          -- ^ Initial parameters-       -> [(a, r)]                               -- ^ Samples-       -> Integer                                -- ^ Maximum number of iterations-       -> LMA_I.Options r                        -- ^ Minimization options-       -> Maybe (SizedList m r)                  -- ^ Optional lower bounds-       -> Maybe (SizedList m r)                  -- ^ Optional upper bounds-       -> Maybe (LinearConstraints k m r)        -- ^ Optional linear constraints-       -> Maybe (SizedList m r)                  -- ^ Optional weights-       -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r)-levmar model mJac params ys itMax opts mLowBs mUpBs mLinC mWghts =-    fmap convertResult $ LMA_I.levmar (convertModel model)-                                      (fmap convertJacob mJac)-                                      (toList params)-                                      ys-                                      itMax-                                      opts-                                      (fmap toList mLowBs)-                                      (fmap toList mUpBs)-                                      (fmap convertLinearConstraints mLinC)-                                      (fmap toList mWghts)-    where-      convertModel mdl = \ps   ->          mdl $* (unsafeFromList ps :: SizedList m r)-      convertJacob jac = \ps x -> toList ((jac $* (unsafeFromList ps :: SizedList m r)) x :: SizedList m r)----- The End ---------------------------------------------------------------------

@@ -1,139 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE FlexibleContexts #-}------------------------------------------------------------------------------------- |--- Module      :  LevMar.Fitting.AD--- Copyright   :  (c) 2009 Roel van Dijk & Bas van Dijk--- License     :  BSD-style (see the file LICENSE)------ Maintainer  :  vandijk.roel@gmail.com, v.dijk.bas@gmail.com--- Stability   :  Experimental------ This module provides the Levenberg-Marquardt algorithm specialised--- for curve-fitting that uses Automatic Differentiation to--- automatically compute the Jacobian.------ For additional documentation see the documentation of the levmar C--- library which this library is based on:--- <http://www.ics.forth.gr/~lourakis/levmar/>--------------------------------------------------------------------------------------module LevMar.Fitting.AD-    ( -- * Model.-      LMA.Model-    , LMA.SimpleModel--      -- * Levenberg-Marquardt algorithm.-    , LMA_I.LevMarable-    , levmar--    , LinearConstraints-    , noLinearConstraints-    , Matrix--    -- * Minimization options.-    , LMA_I.Options(..)-    , LMA_I.defaultOpts--      -- * Output-    , LMA_I.Info(..)-    , LMA_I.StopReason(..)-    , CovarMatrix--    , LMA_I.LevMarError(..)--      -- *Type-level machinery-    , Z, S, Nat-    , SizedList(..)-    , NFunction-    ) where---import qualified LevMar.Fitting              as LMA-import qualified LevMar.Intermediate.Fitting as LMA_I--import LevMar.Utils ( LinearConstraints-                    , noLinearConstraints-                    , convertLinearConstraints-                    , Matrix-                    , CovarMatrix-                    , convertResult-                    )--import TypeLevelNat ( Z, S, Nat )-import SizedList    ( SizedList(..), toList, unsafeFromList )-import NFunction    ( NFunction, ($*) )--import LevMar.Utils.AD  ( value, firstDeriv, constant, idDAt )---- From vector-space:-import Data.Derivative  ( (:~>) )-import Data.VectorSpace ( VectorSpace, Scalar )-import Data.Basis       ( HasBasis, Basis )-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm specialised for curve-fitting--- that automatically computes the 'Jacobian' using automatic--- differentiation of the model function.------ /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of--- the 'Eq' and 'Ord' classes to the parameters. These methods are--- undefined for ':~>'!!!-levmar :: forall m k r a.-          ( Nat m-          , Nat k-          , HasBasis r-          , Basis r ~ ()-          , VectorSpace (Scalar r)-          , LMA_I.LevMarable r-          )-       => LMA.Model m (r :~> r) a             -- ^ Model. Note that-                                              --   ':~>' is overloaded-                                              --   for all the numeric-                                              --   classes.-       -> SizedList m r                       -- ^ Initial parameters-       -> [(a, r)]                            -- ^ Samples-       -> Integer                             -- ^ Maximum number of iterations-       -> LMA_I.Options r                       -- ^ Minimization options-       -> Maybe (SizedList m r)               -- ^ Optional lower bounds-       -> Maybe (SizedList m r)               -- ^ Optional upper bounds-       -> Maybe (LinearConstraints k m r)     -- ^ Optional linear constraints-       -> Maybe (SizedList m r)               -- ^ Optional weights-       -> Either LMA_I.LevMarError (SizedList m r, LMA_I.Info r, CovarMatrix m r)--levmar model params ys itMax opts mLowBs mUpBs mLinC mWghts =-    fmap convertResult $ LMA_I.levmar (convertModel model)-                                      (Just $ jacobianOf model)-                                      (toList params)-                                      ys-                                      itMax-                                      opts-                                      (fmap toList mLowBs)-                                      (fmap toList mUpBs)-                                      (fmap convertLinearConstraints mLinC)-                                      (fmap toList mWghts)-    where-      convertModel :: LMA.Model m (r :~> r) a -> LMA_I.Model r a-      (convertModel f) ps x = value $ (f $* pDs :: a -> r :~> r) x-          where-            pDs :: SizedList m (r :~> r)-            pDs = unsafeFromList $ fmap constant ps--      jacobianOf :: LMA.Model m (r :~> r) a -> LMA_I.Jacobian r a-      (jacobianOf f) ps x = fmap combine $ zip [0..] ps-          where-            combine (ix, p) = firstDeriv $ (f $* pDs :: a -> r :~> r) x p-                where-                  pDs :: SizedList m (r :~> r)-                  pDs = unsafeFromList $ idDAt ix ps----- The End ---------------------------------------------------------------------

@@ -1,409 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE FlexibleInstances #-}------------------------------------------------------------------------------------- |--- Module      :  LevMar.Intermediate--- Copyright   :  (c) 2009 Roel van Dijk & Bas van Dijk--- License     :  BSD-style (see the file LICENSE)------ Maintainer  :  vandijk.roel@gmail.com, v.dijk.bas@gmail.com--- Stability   :  Experimental------------ For additional documentation see the documentation of the levmar C--- library which this library is based on:--- <http://www.ics.forth.gr/~lourakis/levmar/>--------------------------------------------------------------------------------------module LevMar.Intermediate-    ( -- * Model & Jacobian.-       Model-    , Jacobian--      -- * Levenberg-Marquardt algorithm.-    , LevMarable-    , levmar--    , LinearConstraints--      -- * Minimization options.-    , Options(..)-    , defaultOpts--      -- * Output-    , Info(..)-    , StopReason(..)-    , CovarMatrix--    , LevMarError(..)-    ) where---import Foreign.Marshal.Array ( allocaArray, peekArray, pokeArray, withArray )-import Foreign.Ptr           ( Ptr, nullPtr, plusPtr )-import Foreign.Storable      ( Storable )-import Foreign.C.Types       ( CInt )-import System.IO.Unsafe      ( unsafePerformIO )-import Data.Maybe            ( fromJust, fromMaybe, isJust )-import Control.Monad.Instances -- for 'instance Functor (Either a)'--import qualified Bindings.LevMar.CurryFriendly as LMA_C-------------------------------------------------------------------------------------- Model & Jacobian.-----------------------------------------------------------------------------------{- | A functional relation describing measurements represented as a function-from a list of parameters to a list of expected measurements.-- * Ensure that the length of the parameters list equals the length of the-   initial parameters list in 'levmar'.-- * Ensure that the length of the ouput list equals the length of the samples-   list in 'levmar'.--For example:--@-hatfldc :: Model Double-hatfldc [p0, p1, p2, p3] = [ p0 - 1.0-                           , p0 - sqrt p1-                           , p1 - sqrt p2-                           , p3 - 1.0-                           ]-@--}-type Model r = [r] -> [r]--{- | The jacobian of the 'Model' function. Expressed as a function from a list-of parameters to a list of lists which for each expected measurement describes-the partial derivatives of the parameters.--See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>-- * Ensure that the length of the parameter list equals the length of the initial-   parameter list in 'levmar'.-- * Ensure that the output matrix has the dimension @n@/x/@m@ where @n@ is the-   number of samples and @m@ is the number of parameters.--For example the jacobian of the above @hatfldc@ model is:--@-hatfldc_jac :: Jacobian Double-hatfldc_jac _ p1 p2 _ = [ [1.0,  0.0,           0.0,           0.0]-                        , [1.0, -0.5 / sqrt p1, 0.0,           0.0]-                        , [0.0,  1.0,          -0.5 / sqrt p2, 0.0]-                        , [0.0,  0.0,           0.0,           1.0]-                        ]-@--}-type Jacobian r = [r] -> [[r]]-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm is overloaded to work on 'Double' and 'Float'.-class LevMarable r where--    -- | The Levenberg-Marquardt algorithm.-    levmar :: Model r                     -- ^ Model-           -> Maybe (Jacobian r)          -- ^ Optional jacobian-           -> [r]                         -- ^ Initial parameters-           -> [r]                         -- ^ Samples-           -> Integer                     -- ^ Maximum iterations-           -> Options r                   -- ^ Minimization options-           -> Maybe [r]                   -- ^ Optional lower bounds-           -> Maybe [r]                   -- ^ Optional upper bounds-           -> Maybe (LinearConstraints r) -- ^ Optional linear constraints-           -> Maybe [r]                   -- ^ Optional weights-           -> Either LevMarError ([r], Info r, CovarMatrix r)--instance LevMarable Float where-    levmar = gen_levmar LMA_C.slevmar_der-                        LMA_C.slevmar_dif-                        LMA_C.slevmar_bc_der-                        LMA_C.slevmar_bc_dif-                        LMA_C.slevmar_lec_der-                        LMA_C.slevmar_lec_dif-                        LMA_C.slevmar_blec_der-                        LMA_C.slevmar_blec_dif--instance LevMarable Double where-    levmar = gen_levmar LMA_C.dlevmar_der-                        LMA_C.dlevmar_dif-                        LMA_C.dlevmar_bc_der-                        LMA_C.dlevmar_bc_dif-                        LMA_C.dlevmar_lec_der-                        LMA_C.dlevmar_lec_dif-                        LMA_C.dlevmar_blec_der-                        LMA_C.dlevmar_blec_dif--{- | @gen_levmar@ takes the low-level C functions as arguments and-executes one of them depending on the optional jacobian and constraints.--Preconditions:-  length ys >= length ps--     isJust mLowBs && length (fromJust mLowBs) == length ps-  && isJust mUpBs  && length (fromJust mUpBs)  == length ps--  boxConstrained && (all $ zipWith (<=) (fromJust mLowBs) (fromJust mUpBs))--}-gen_levmar :: forall cr r. (Storable cr, RealFrac cr, Real r, Fractional r)-           => LMA_C.LevMarDer cr-           -> LMA_C.LevMarDif cr-           -> LMA_C.LevMarBCDer cr-           -> LMA_C.LevMarBCDif cr-           -> LMA_C.LevMarLecDer cr-           -> LMA_C.LevMarLecDif cr-           -> LMA_C.LevMarBLecDer cr-           -> LMA_C.LevMarBLecDif cr--           -> Model r                     -- ^ Model-           -> Maybe (Jacobian r)          -- ^ Optional jacobian-           -> [r]                         -- ^ Initial parameters-           -> [r]                         -- ^ Samples-           -> Integer                     -- ^ Maximum iterations-           -> Options r                   -- ^ Options-           -> Maybe [r]                   -- ^ Optional lower bounds-           -> Maybe [r]                   -- ^ Optional upper bounds-           -> Maybe (LinearConstraints r) -- ^ Optional linear constraints-           -> Maybe [r]                   -- ^ Optional weights-           -> Either LevMarError ([r], Info r, CovarMatrix r)-gen_levmar f_der-           f_dif-           f_bc_der-           f_bc_dif-           f_lec_der-           f_lec_dif-           f_blec_der-           f_blec_dif-           model mJac ps ys itMax opts mLowBs mUpBs mLinC mWeights-    = unsafePerformIO .-        withArray (map realToFrac ps) $ \psPtr ->-        withArray (map realToFrac ys) $ \ysPtr ->-        withArray (map realToFrac $ optsToList opts) $ \optsPtr ->-        allocaArray LMA_C._LM_INFO_SZ $ \infoPtr ->-        allocaArray covarLen $ \covarPtr ->-        LMA_C.withModel (convertModel model) $ \modelPtr -> do--          let runDif :: LMA_C.LevMarDif cr -> IO CInt-              runDif f = f modelPtr-                           psPtr-                           ysPtr-                           (fromIntegral lenPs)-                           (fromIntegral lenYs)-                           (fromIntegral itMax)-                           optsPtr-                           infoPtr-                           nullPtr-                           covarPtr-                           nullPtr--          r <- case mJac of-                 Just jac -> LMA_C.withJacobian (convertJacobian jac) $ \jacobPtr ->-                               let runDer :: LMA_C.LevMarDer cr -> IO CInt-                                   runDer f = runDif $ f jacobPtr-                               in if boxConstrained-                                  then if linConstrained-                                       then withBoxConstraints (withLinConstraints $ withWeights runDer) f_blec_der-                                       else withBoxConstraints runDer f_bc_der-                                  else if linConstrained-                                       then withLinConstraints runDer f_lec_der-                                       else runDer f_der--                 Nothing -> if boxConstrained-                            then if linConstrained-                                 then withBoxConstraints (withLinConstraints $ withWeights runDif) f_blec_dif-                                 else withBoxConstraints runDif f_bc_dif-                            else if linConstrained-                                 then withLinConstraints runDif f_lec_dif-                                 else runDif f_dif--          if    r < 0-             && r /= LMA_C._LM_ERROR_SINGULAR_MATRIX -- we don't treat these two as an error-             && r /= LMA_C._LM_ERROR_SUM_OF_SQUARES_NOT_FINITE-            then return . Left $ convertLevMarError r-            else do result <- peekArray lenPs psPtr-                    info   <- peekArray LMA_C._LM_INFO_SZ infoPtr--                    let covarPtrEnd = plusPtr covarPtr covarLen-                    let convertCovarMatrix ptr-                            | ptr == covarPtrEnd = return []-                            | otherwise = do row <- peekArray lenPs ptr-                                             rows <- convertCovarMatrix $ plusPtr ptr lenPs-                                             return $ row : rows--                    covar  <- convertCovarMatrix covarPtr--                    return $ Right ( map realToFrac result-                                   , listToInfo info-                                   , map (map realToFrac) covar-                                   )-    where-      lenPs          = length ps-      lenYs          = length ys-      covarLen       = lenPs * lenPs-      (cMat, rhcVec) = fromJust mLinC--      -- Whether the parameters are constrained by a linear equation.-      linConstrained = isJust mLinC--      -- Whether the parameters are constrained by a bounding box.-      boxConstrained = isJust mLowBs || isJust mUpBs--      withBoxConstraints f g = maybeWithArray ((fmap . fmap) realToFrac mLowBs) $ \lBsPtr ->-                                 maybeWithArray ((fmap . fmap) realToFrac mUpBs) $ \uBsPtr ->-                                   f $ g lBsPtr uBsPtr--      withLinConstraints f g = withArray (map realToFrac $ concat cMat) $ \cMatPtr ->-                                 withArray (map realToFrac rhcVec) $ \rhcVecPtr ->-                                   f . g cMatPtr rhcVecPtr . fromIntegral $ length cMat--      withWeights f g = maybeWithArray ((fmap . fmap) realToFrac mWeights) $ f . g--convertModel :: (Real r, Fractional r, Storable c, Real c, Fractional c)-             =>  Model r -> LMA_C.Model c-convertModel model = \parPtr hxPtr numPar _ _ -> do-                       params <- peekArray (fromIntegral numPar) parPtr-                       pokeArray hxPtr . map realToFrac . model $ map realToFrac params--convertJacobian :: (Real r, Fractional r, Storable c, Real c, Fractional c)-                => Jacobian r -> LMA_C.Jacobian c-convertJacobian jac = \parPtr jPtr numPar _ _ -> do-                        params <- peekArray (fromIntegral numPar) parPtr-                        pokeArray jPtr . concatMap (map realToFrac) . jac $ map realToFrac params--maybeWithArray :: Storable a => Maybe [a] -> (Ptr a -> IO b) -> IO b-maybeWithArray Nothing   f = f nullPtr-maybeWithArray (Just xs) f = withArray xs f----- | Linear constraints consisting of a constraints matrix, /kxm/ and---   a right hand constraints vector, /kx1/ where /m/ is the number of---   parameters and /k/ is the number of constraints.-type LinearConstraints r = ([[r]], [r])-------------------------------------------------------------------------------------- Minimization options.------------------------------------------------------------------------------------- | Minimization options-data Options r =-    Opts { optScaleInitMu      :: r -- ^ Scale factor for initial /mu/.-         , optStopNormInfJacTe :: r -- ^ Stopping thresholds for @||J^T e||_inf@.-         , optStopNorm2Dp      :: r -- ^ Stopping thresholds for @||Dp||_2@.-         , optStopNorm2E       :: r -- ^ Stopping thresholds for @||e||_2@.-         , optDelta            :: r -- ^ Step used in the difference approximation to the Jacobian.-                                    --   If @optDelta<0@, the Jacobian is approximated-                                    --   with central differences which are more accurate-                                    --   (but slower!) compared to the forward differences-                                    --   employed by default.-         } deriving Show---- | Default minimization options-defaultOpts :: Fractional r => Options r-defaultOpts = Opts { optScaleInitMu      = LMA_C._LM_INIT_MU-                   , optStopNormInfJacTe = LMA_C._LM_STOP_THRESH-                   , optStopNorm2Dp      = LMA_C._LM_STOP_THRESH-                   , optStopNorm2E       = LMA_C._LM_STOP_THRESH-                   , optDelta            = LMA_C._LM_DIFF_DELTA-                   }--optsToList :: Options r -> [r]-optsToList (Opts mu  eps1  eps2  eps3  delta) =-                [mu, eps1, eps2, eps3, delta]-------------------------------------------------------------------------------------- Output------------------------------------------------------------------------------------- | Information regarding the minimization.-data Info r = Info { infNorm2initE      :: r          -- ^ @||e||_2@             at initial   parameters.-                   , infNorm2E          :: r          -- ^ @||e||_2@             at estimated parameters.-                   , infNormInfJacTe    :: r          -- ^ @||J^T e||_inf@       at estimated parameters.-                   , infNorm2Dp         :: r          -- ^ @||Dp||_2@            at estimated parameters.-                   , infMuDivMax        :: r          -- ^ @\mu/max[J^T J]_ii ]@ at estimated parameters.-                   , infNumIter         :: Integer    -- ^ Number of iterations.-                   , infStopReason      :: StopReason -- ^ Reason for terminating.-                   , infNumFuncEvals    :: Integer    -- ^ Number of function evaluations.-                   , infNumJacobEvals   :: Integer    -- ^ Number of jacobian evaluations.-                   , infNumLinSysSolved :: Integer    -- ^ Number of linear systems solved, i.e. attempts for reducing error.-                   } deriving Show--listToInfo :: (RealFrac cr, Fractional r) => [cr] -> Info r-listToInfo [a,b,c,d,e,f,g,h,i,j] =-    Info { infNorm2initE      = realToFrac a-         , infNorm2E          = realToFrac b-         , infNormInfJacTe    = realToFrac c-         , infNorm2Dp         = realToFrac d-         , infMuDivMax        = realToFrac e-         , infNumIter         = floor f-         , infStopReason      = toEnum $ floor g - 1-         , infNumFuncEvals    = floor h-         , infNumJacobEvals   = floor i-         , infNumLinSysSolved = floor j-         }-listToInfo _ = error "liftToInfo: wrong list length"---- | Reason for terminating.-data StopReason = SmallGradient  -- ^ Stopped because of small gradient @J^T e@.-                | SmallDp        -- ^ Stopped because of small Dp.-                | MaxIterations  -- ^ Stopped because maximum iterations was reached.-                | SingularMatrix -- ^ Stopped because of singular matrix. Restart from current estimated parameters with increased 'optScaleInitMu'.-                | SmallestError  -- ^ Stopped because no further error reduction is possible. Restart with increased 'optScaleInitMu'.-                | SmallNorm2E    -- ^ Stopped because of small @||e||_2@.-                | InvalidValues  -- ^ Stopped because model function returned invalid values (i.e. NaN or Inf). This is a user error.-                  deriving (Show, Enum)---- | Covariance matrix corresponding to LS solution.-type CovarMatrix r = [[r]]-------------------------------------------------------------------------------------- Error-----------------------------------------------------------------------------------data LevMarError-    = LevMarError                    -- ^ Generic error (not one of the others)-    | LapackError                    -- ^ A call to a lapack subroutine failed in the underlying C levmar library.-    | FailedBoxCheck                 -- ^ At least one lower bound exceeds the upper one.-    | MemoryAllocationFailure        -- ^ A call to @malloc@ failed in the underlying C levmar library.-    | ConstraintMatrixRowsGtCols     -- ^ The matrix of constraints cannot have more rows than columns.-    | ConstraintMatrixNotFullRowRank -- ^ Constraints matrix is not of full row rank.-    | TooFewMeasurements             -- ^ Cannot solve a problem with fewer measurements than unknowns.-                                     --   In case linear constraints are provided, this error is also returned-                                     --   when the number of measurements is smaller than the number of unknowns-                                     --   minus the number of equality constraints.-      deriving Show--levmarCErrorToLevMarError :: [(CInt, LevMarError)]-levmarCErrorToLevMarError =-    [ (LMA_C._LM_ERROR,                                     LevMarError)-    , (LMA_C._LM_ERROR_LAPACK_ERROR,                        LapackError)-  --, (LMA_C._LM_ERROR_NO_JACOBIAN,                         can never happen)-  --, (LMA_C._LM_ERROR_NO_BOX_CONSTRAINTS,                  can never happen)-    , (LMA_C._LM_ERROR_FAILED_BOX_CHECK,                    FailedBoxCheck)-    , (LMA_C._LM_ERROR_MEMORY_ALLOCATION_FAILURE,           MemoryAllocationFailure)-    , (LMA_C._LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS,      ConstraintMatrixRowsGtCols)-    , (LMA_C._LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK, ConstraintMatrixNotFullRowRank)-    , (LMA_C._LM_ERROR_TOO_FEW_MEASUREMENTS,                TooFewMeasurements)-  --, (LMA_C._LM_ERROR_SINGULAR_MATRIX,                     we don't treat this as an error)-  --, (LMA_C._LM_ERROR_SUM_OF_SQUARES_NOT_FINITE,           we don't treat this as an error)-    ]--convertLevMarError :: CInt -> LevMarError-convertLevMarError err = fromMaybe (error "Unknown levmar error") $-                         lookup err levmarCErrorToLevMarError----- The End ---------------------------------------------------------------------

@@ -1,106 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE FlexibleContexts #-}------------------------------------------------------------------------------------- |--- Module      :  LevMar.Intermediate.AD--- Copyright   :  (c) 2009 Roel van Dijk & Bas van Dijk--- License     :  BSD-style (see the file LICENSE)------ Maintainer  :  vandijk.roel@gmail.com, v.dijk.bas@gmail.com--- Stability   :  Experimental------ A levmar variant that uses Automatic Differentiation to--- automatically compute the Jacobian.------ For additional documentation see the documentation of the levmar C--- library which this library is based on:--- <http://www.ics.forth.gr/~lourakis/levmar/>--------------------------------------------------------------------------------------module LevMar.Intermediate.AD-    ( -- * Model.-      LMA_I.Model-    , LMA_I.Jacobian-    , jacobianOf--      -- * Levenberg-Marquardt algorithm.-    , LMA_I.LevMarable-    , levmar--    , LMA_I.LinearConstraints--      -- * Minimization options.-    , LMA_I.Options(..)-    , LMA_I.defaultOpts--      -- * Output-    , LMA_I.Info(..)-    , LMA_I.StopReason(..)-    , LMA_I.CovarMatrix--    , LMA_I.LevMarError(..)-    ) where---import qualified LevMar.Intermediate as LMA_I--import LevMar.Utils.AD  ( value, firstDeriv, constant, idDAt )---- From vector-space:-import Data.Derivative  ( (:~>) )-import Data.VectorSpace ( VectorSpace, Scalar )-import Data.Basis       ( HasBasis, Basis )--import Data.List        ( transpose )-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm that automatically computes the--- 'Jacobian' using automatic differentiation of the model function.------ /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of--- the 'Eq' and 'Ord' classes to the parameters. These methods are--- undefined for ':~>'!!!-levmar :: forall r.-          ( HasBasis r-          , Basis r ~ ()-          , VectorSpace (Scalar r)-          , LMA_I.LevMarable r-          )-       => LMA_I.Model (r :~> r)             -- ^ Model. Note that-                                            --   ':~>' is overloaded-                                            --   for all the numeric-                                            --   classes.-       -> [r]                               -- ^ Initial parameters-       -> [r]                               -- ^ Samples-       -> Integer                           -- ^ Maximum iterations-       -> LMA_I.Options r                   -- ^ Minimization options-       -> Maybe [r]                         -- ^ Optional lower bounds-       -> Maybe [r]                         -- ^ Optional upper bounds-       -> Maybe (LMA_I.LinearConstraints r) -- ^ Optional linear constraints-       -> Maybe [r]                         -- ^ Optional weights-       -> Either LMA_I.LevMarError ([r], LMA_I.Info r, LMA_I.CovarMatrix r)--levmar model = LMA_I.levmar (convertModel model) . Just $ jacobianOf model-    where-      convertModel :: LMA_I.Model (r :~> r) -> LMA_I.Model r-      convertModel mdl = map value . mdl . map constant---- | Compute the 'LMA_I.Jacobian' of the 'LMA_I.Model' using Automatic--- Differentiation.-jacobianOf :: (HasBasis r, Basis r ~ (), VectorSpace (Scalar r))-           => LMA_I.Model (r :~> r) -> LMA_I.Jacobian r-(jacobianOf mdl) ps = map (\fs -> zipWith (firstDeriv .) fs ps)-                    . transpose $ map mdl pDs-    where-      pDs = [idDAt n ps | n <- [0 .. length ps - 1]]----- The End ---------------------------------------------------------------------

@@ -1,131 +0,0 @@------------------------------------------------------------------------------------ |--- Module      :  LevMar.Intermediate.Fitting--- Copyright   :  (c) 2009 Roel van Dijk & Bas van Dijk--- License     :  BSD-style (see the file LICENSE)------ Maintainer  :  vandijk.roel@gmail.com, v.dijk.bas@gmail.com--- Stability   :  Experimental------ This module provides the Levenberg-Marquardt algorithm specialised--- for curve-fitting.------ For additional documentation see the documentation of the levmar C--- library which this library is based on:--- <http://www.ics.forth.gr/~lourakis/levmar/>--------------------------------------------------------------------------------------module LevMar.Intermediate.Fitting-    ( -- * Model & Jacobian.-      Model-    , SimpleModel-    , Jacobian-    , SimpleJacobian--      -- * Levenberg-Marquardt algorithm.-    , LMA_I.LevMarable-    , levmar--    , LMA_I.LinearConstraints--      -- * Minimization options.-    , LMA_I.Options(..)-    , LMA_I.defaultOpts--      -- * Output-    , LMA_I.Info(..)-    , LMA_I.StopReason(..)-    , LMA_I.CovarMatrix--    , LMA_I.LevMarError(..)-    ) where---import qualified LevMar.Intermediate as LMA_I-------------------------------------------------------------------------------------- Model & Jacobian.-----------------------------------------------------------------------------------{- | A functional relation describing measurements represented as a function-from a list of parameters and an x-value to an expected measurement.-- * Ensure that the length of the parameters list equals the lenght of the initial-   parameters list in 'levmar'.--For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be-written as:--@-quad :: 'Num' r => 'Model' r r-quad [a, b, c] x = a*x^2 + b*x + c-@--}-type Model r a = [r] -> a -> r---- | This type synonym expresses that usually the @a@ in @'Model' r a@--- equals the type of the parameters.-type SimpleModel r = Model r r--{- | The jacobian of the 'Model' function. Expressed as a function from a list-of parameters and an x-value to the partial derivatives of the parameters.--See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>-- * Ensure that the length of the parameters list equals the lenght of the initial-   parameters list in 'levmar'.-- * Ensure that the length of the output parameter derivatives list equals the-   length of the input parameters list.--For example, the jacobian of the above @quad@ model can be written as:--@-quadJacob :: 'Num' r => 'Jacobian' N3 r r-quadJacob [_, _, _] x = [ x^2   -- with respect to a-                        , x     -- with respect to b-                        , 1     -- with respect to c-                        ]-@--Notice you don't have to differentiate for @x@.--}-type Jacobian r a = [r] -> a -> [r]---- | This type synonym expresses that usually the @a@ in @'Jacobian' r a@--- equals the type of the parameters.-type SimpleJacobian r = Jacobian r r-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm specialised for curve-fitting.-levmar :: LMA_I.LevMarable r-       => Model r a                         -- ^ Model-       -> Maybe (Jacobian r a)              -- ^ Optional jacobian-       -> [r]                               -- ^ Initial parameters-       -> [(a, r)]                          -- ^ Samples-       -> Integer                           -- ^ Maximum iterations-       -> LMA_I.Options r                   -- ^ Minimization options-       -> Maybe [r]                         -- ^ Optional lower bounds-       -> Maybe [r]                         -- ^ Optional upper bounds-       -> Maybe (LMA_I.LinearConstraints r) -- ^ Optional linear constraints-       -> Maybe [r]                         -- ^ Optional weights-       -> Either LMA_I.LevMarError ([r], LMA_I.Info r, LMA_I.CovarMatrix r)-levmar model mJac params samples =-    LMA_I.levmar (convertModel model)-                 (fmap convertJacob mJac)-                 params-                 ys-        where-          (xs, ys) = unzip samples--          convertModel mdl = \ps -> map (mdl ps) xs-          convertJacob jac = \ps -> map (jac ps) xs----- The End ---------------------------------------------------------------------

@@ -1,106 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE FlexibleContexts #-}------------------------------------------------------------------------------------- |--- Module      :  LevMar.Intermediate.Fitting.AD--- Copyright   :  (c) 2009 Roel van Dijk & Bas van Dijk--- License     :  BSD-style (see the file LICENSE)------ Maintainer  :  vandijk.roel@gmail.com, v.dijk.bas@gmail.com--- Stability   :  Experimental------ A levmar variant specialised for curve-fitting that uses Automatic--- Differentiation to automatically compute the Jacobian.------ For additional documentation see the documentation of the levmar C--- library which this library is based on:--- <http://www.ics.forth.gr/~lourakis/levmar/>--------------------------------------------------------------------------------------module LevMar.Intermediate.Fitting.AD-    ( -- * Model.-      LMA_I.Model-    , LMA_I.SimpleModel-    , LMA_I.Jacobian-    , LMA_I.SimpleJacobian-    , jacobianOf--      -- * Levenberg-Marquardt algorithm.-    , LMA_I.LevMarable-    , levmar--    , LMA_I.LinearConstraints--      -- * Minimization options.-    , LMA_I.Options(..)-    , LMA_I.defaultOpts--      -- * Output-    , LMA_I.Info(..)-    , LMA_I.StopReason(..)-    , LMA_I.CovarMatrix--    , LMA_I.LevMarError(..)-    ) where---import qualified LevMar.Intermediate.Fitting as LMA_I--import LevMar.Utils.AD  ( value, firstDeriv, constant, idDAt )---- From vector-space:-import Data.Derivative  ( (:~>) )-import Data.VectorSpace ( VectorSpace, Scalar )-import Data.Basis       ( HasBasis, Basis )-------------------------------------------------------------------------------------- Levenberg-Marquardt algorithm.------------------------------------------------------------------------------------- | The Levenberg-Marquardt algorithm specialised for curve-fitting--- that automatically computes the 'Jacobian' using automatic--- differentiation of the model function.------ /Warning/: Don't apply 'levmar' to 'LMA_I.Model's that apply methods of--- the 'Eq' and 'Ord' classes to the parameters. These methods are--- undefined for ':~>'!!!-levmar :: forall r a.-          ( HasBasis r-          , Basis r ~ ()-          , VectorSpace (Scalar r)-          , LMA_I.LevMarable r-          )-       => LMA_I.Model (r :~> r) a           -- ^ Model. Note that-                                            --   ':~>' is overloaded-                                            --   for all the numeric-                                            --   classes.-       -> [r]                               -- ^ Initial parameters-       -> [(a, r)]                          -- ^ Samples-       -> Integer                           -- ^ Maximum iterations-       -> LMA_I.Options r                   -- ^ Minimization options-       -> Maybe [r]                         -- ^ Optional lower bounds-       -> Maybe [r]                         -- ^ Optional upper bounds-       -> Maybe (LMA_I.LinearConstraints r) -- ^ Optional linear constraints-       -> Maybe [r]                         -- ^ Optional weights-       -> Either LMA_I.LevMarError ([r], LMA_I.Info r, LMA_I.CovarMatrix r)--levmar model = LMA_I.levmar (convertModel model) . Just $ jacobianOf model-    where-      convertModel :: LMA_I.Model (r :~> r) a -> LMA_I.Model r a-      convertModel mdl = \ps -> value . mdl (map constant ps)---- | Compute the 'LMA_I.Jacobian' of the 'LMA_I.Model' using Automatic--- Differentiation.-jacobianOf :: (HasBasis r, Basis r ~ (), VectorSpace (Scalar r))-           => LMA_I.Model (r :~> r) a -> LMA_I.Jacobian r a-jacobianOf mdl =-    \ps x -> map (\(ix, p) -> firstDeriv $ mdl (idDAt ix ps) x p) $-                 zip [0..] ps----- The End ---------------------------------------------------------------------

@@ -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-                                        )

@@ -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

@@ -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 .--}

@@ -0,0 +1,527 @@+{-# LANGUAGE CPP+           , NoImplicitPrelude+           , UnicodeSyntax+           , ScopedTypeVariables+           , DeriveDataTypeable+  #-}++--------------------------------------------------------------------------------+-- |+-- Module:     Numeric.LevMar+-- Copyright:  (c) 2009 - 2010 Roel van Dijk & Bas van Dijk+-- License:    BSD-style (see the file LICENSE)+-- Maintainer: Roel van Dijk <vandijk.roel@gmail.com>+--             Bas van Dijk <v.dijk.bas@gmail.com>+-- Stability:  Experimental+--+-- For additional documentation see the documentation of the levmar C+-- library which this library is based on:+-- <http://www.ics.forth.gr/~lourakis/levmar/>+--+--------------------------------------------------------------------------------++module Numeric.LevMar+    ( -- * Model & Jacobian.+      Model+    , Jacobian++      -- * Levenberg-Marquardt algorithm.+    , LevMarable(levmar)++    , LinearConstraints++      -- * Minimization options.+    , Options(..)+    , defaultOpts++      -- * Constraints+    , Constraints(..)+    , noConstraints++      -- * Output+    , Info(..)+    , StopReason(..)+    , CovarMatrix++    , LevMarError(..)+    ) where+++-------------------------------------------------------------------------------+-- Imports+-------------------------------------------------------------------------------++-- from base:+import Control.Monad.Instances -- for 'instance Functor (Either a)'+import Control.Exception     ( Exception )+import Data.Typeable         ( Typeable )+import Data.Bool             ( otherwise )+import Data.Either           ( Either(Left, Right) )+import Data.Function         ( ($) )+import Data.List             ( lookup, map, concat, concatMap, length )+import Data.Maybe            ( Maybe(Nothing, Just)+                             , isJust, fromJust, fromMaybe+                             )+import Data.Ord              ( (<) )+import Foreign.Marshal.Array ( allocaArray, peekArray, pokeArray, withArray )+import Foreign.Ptr           ( Ptr, nullPtr, plusPtr )+import Foreign.Storable      ( Storable )+import Foreign.C.Types       ( CInt )+import Prelude               ( Enum, Fractional, Real, RealFrac+                             , Integer, Float, Double+                             , fromIntegral, realToFrac, toEnum+                             , (-), error, floor+                             )+import System.IO             ( IO )+import System.IO.Unsafe      ( unsafePerformIO )+import Text.Read             ( Read )+import Text.Show             ( Show )++#if __GLASGOW_HASKELL__ < 700+import Prelude               ( fromInteger )+#endif++-- from base-unicode-symbols:+import Data.Bool.Unicode     ( (∧), (∨) )+import Data.Eq.Unicode       ( (≡), (≢) )+import Data.Function.Unicode ( (∘) )+import Prelude.Unicode       ( (⋅) )++-- 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.+--------------------------------------------------------------------------------++{-| A functional relation describing measurements represented as a function+from a list of parameters to a list of expected measurements.++ * Ensure that the length of the parameters list equals the length of the+   initial parameters list in 'levmar'.++ * Ensure that the length of the ouput list equals the length of the samples+   list in 'levmar'.++For example:++@+hatfldc :: Model Double+hatfldc [p0, p1, p2, p3] = [ p0 - 1.0+                           , p0 - sqrt p1+                           , p1 - sqrt p2+                           , p3 - 1.0+                           ]+@+-}+type Model r = [r] → [r]++{-| The jacobian of the 'Model' function. Expressed as a function from a list+of parameters to a list of lists which for each expected measurement describes+the partial derivatives of the parameters.++See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>++ * Ensure that the length of the parameter list equals the length of the initial+   parameter list in 'levmar'.++ * Ensure that the output matrix has the dimension @n@/x/@m@ where @n@ is the+   number of samples and @m@ is the number of parameters.++For example the jacobian of the above @hatfldc@ model is:++@+hatfldc_jac :: Jacobian Double+hatfldc_jac _ p1 p2 _ = [ [1.0,  0.0,           0.0,           0.0]+                        , [1.0, -0.5 / sqrt p1, 0.0,           0.0]+                        , [0.0,  1.0,          -0.5 / sqrt p2, 0.0]+                        , [0.0,  0.0,           0.0,           1.0]+                        ]+@+-}+type Jacobian r = [r] → [[r]]+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm is overloaded to work on 'Double' and 'Float'.+class LevMarable r where++    -- | The Levenberg-Marquardt algorithm.+    levmar ∷ Model r            -- ^ Model+           → Maybe (Jacobian r) -- ^ Optional jacobian+           → [r]                -- ^ Initial parameters+           → [r]                -- ^ Samples+           → Integer            -- ^ Maximum iterations+           → Options r          -- ^ Minimization options+           → Constraints r      -- ^ Constraints+           → Either LevMarError ([r], Info r, CovarMatrix 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 ∷ ∀ cr r. (Storable cr, RealFrac cr, Real r, Fractional r)+           ⇒ LevMarDer cr+           → LevMarDif cr+           → LevMarBCDer cr+           → LevMarBCDif cr+           → LevMarLecDer cr+           → LevMarLecDif cr+           → LevMarBLecDer cr+           → LevMarBLecDif cr++           → Model r            -- ^ Model+           → Maybe (Jacobian r) -- ^ Optional jacobian+           → [r]                -- ^ Initial parameters+           → [r]                -- ^ Samples+           → Integer            -- ^ Maximum iterations+           → Options r          -- ^ Options+           → Constraints r      -- ^ Constraints+           → 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 (Constraints mLowBs mUpBs mWeights mLinC)+    = unsafePerformIO ∘++       -- Allocation:+       withArray (map realToFrac ps) $ \psPtr →+        withArray (map realToFrac ys) $ \ysPtr →+         withArray (map realToFrac $ optsToList opts) $ \optsPtr →+          allocaArray c'LM_INFO_SZ $ \infoPtr →+           allocaArray covarLen $ \covarPtr →+            withModel (convertModel model) $ \modelPtr → do++              -- Calling the correct low-level levmar function:+              let runDif ∷ 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 → withJacobian (convertJacobian jac) $ \jacobPtr →+                                   let runDer ∷ 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++              -- Handling errors:+              if r < 0+                 ∧ r ≢ c'LM_ERROR_SINGULAR_MATRIX -- we don't treat these two as an error+                 ∧ r ≢ c'LM_ERROR_SUM_OF_SQUARES_NOT_FINITE+                then return $ Left $ convertLevMarError r+                else -- Converting results:+                     do result ← peekArray lenPs psPtr+                        info   ← peekArray c'LM_INFO_SZ infoPtr++                        let convertCovarMatrix ptr+                                | ptr ≡ covarPtr `plusPtr` covarLen = return []+                                | otherwise = do row ← peekArray lenPs ptr+                                                 rows ← convertCovarMatrix $ ptr `plusPtr` lenPs+                                                 return $ map realToFrac row : rows++                        covar ← convertCovarMatrix covarPtr+                        return $ Right ( map realToFrac result+                                       , listToInfo info+                                       , 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 mLowBs $ \lBsPtr →+            maybeWithArray 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 mWeights $ f ∘ g++convertModel ∷ (Real r, Fractional r, Storable c, Real c, Fractional c)+             ⇒ Model r → Bindings.LevMar.Model c+convertModel model =+    \parPtr hxPtr numPar _ _ →+      peekArray (fromIntegral numPar) parPtr >>=+        pokeArray hxPtr ∘ map realToFrac ∘ model ∘ map realToFrac++convertJacobian ∷ (Real r, Fractional r, Storable c, Real c, Fractional c)+                ⇒ Jacobian r → Bindings.LevMar.Jacobian c+convertJacobian jac =+    \parPtr jPtr numPar _ _ →+      peekArray (fromIntegral numPar) parPtr >>=+        pokeArray jPtr ∘ concatMap (map realToFrac) ∘ jac ∘ map realToFrac++-- | 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 (Read, Show)++-- | 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+--------------------------------------------------------------------------------++data Constraints r = Constraints+    { lowerBounds       ∷ Maybe [r]                   -- ^ Optional lower bounds+    , upperBounds       ∷ Maybe [r]                   -- ^ Optional upper bounds+    , weights           ∷ Maybe [r]                   -- ^ Optional weights+    , linearConstraints ∷ Maybe (LinearConstraints r) -- ^ Optional linear constraints+    }++-- | Constraints where all fields are 'Nothing'.+noConstraints ∷ Constraints r+noConstraints = Constraints Nothing Nothing Nothing Nothing++maybeWithArray ∷ (Real α, Fractional r, Storable r)+               ⇒ Maybe [α] → (Ptr r → IO β) → IO β+maybeWithArray Nothing   f = f nullPtr+maybeWithArray (Just xs) f = withArray (map realToFrac xs) f+++--------------------------------------------------------------------------------+-- 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 (Read, 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 (Read, 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, 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") $+                         lookup err levmarCErrorToLevMarError+++-- The End ---------------------------------------------------------------------

@@ -0,0 +1,142 @@+{-# LANGUAGE NoImplicitPrelude, UnicodeSyntax #-}++--------------------------------------------------------------------------------+-- |+-- Module:     Numeric.LevMar.Fitting+-- Copyright:  (c) 2009 - 2010 Roel van Dijk & Bas van Dijk+-- License:    BSD-style (see the file LICENSE)+-- Maintainer: Roel van Dijk <vandijk.roel@gmail.com>+--             Bas van Dijk <v.dijk.bas@gmail.com>+-- Stability:  Experimental+--+-- This module provides the Levenberg-Marquardt algorithm specialised+-- for curve-fitting.+--+-- For additional documentation see the documentation of the levmar C+-- library which this library is based on:+-- <http://www.ics.forth.gr/~lourakis/levmar/>+--+--------------------------------------------------------------------------------++module Numeric.LevMar.Fitting+    ( -- * Model & Jacobian.+      Model+    , SimpleModel+    , Jacobian+    , SimpleJacobian++      -- * Levenberg-Marquardt algorithm.+    , LevMar.LevMarable+    , levmar++    , LevMar.LinearConstraints++      -- * Minimization options.+    , LevMar.Options(..)+    , LevMar.defaultOpts++      -- * Output+    , LevMar.Info(..)+    , LevMar.StopReason(..)+    , LevMar.CovarMatrix++    , LevMar.LevMarError(..)+    ) where+++--------------------------------------------------------------------------------+-- Imports+--------------------------------------------------------------------------------++-- from base:+import Data.Functor  ( fmap )+import Data.Either   ( Either )+import Data.List     ( map, unzip )+import Data.Maybe    ( Maybe )+import Prelude       ( Integer )++-- from levmar:+import qualified Numeric.LevMar as LevMar+++--------------------------------------------------------------------------------+-- Model & Jacobian.+--------------------------------------------------------------------------------++{-| A functional relation describing measurements represented as a function+from a list of parameters and an x-value to an expected measurement.++ * Ensure that the length of the parameters list equals the lenght of the initial+   parameters list in 'levmar'.++For example, the quadratic function @f(x) = a*x^2 + b*x + c@ can be+written as:++@+quad :: 'Num' r => 'Model' r r+quad [a, b, c] x = a*x^2 + b*x + c+@+-}+type Model r a = [r] → (a → r)++-- | This type synonym expresses that usually the @a@ in @'Model' r a@+-- equals the type of the parameters.+type SimpleModel r = Model r r++{-| The jacobian of the 'Model' function. Expressed as a function from a list+of parameters and an x-value to the partial derivatives of the parameters.++See: <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant>++ * Ensure that the length of the parameters list equals the lenght of the initial+   parameters list in 'levmar'.++ * Ensure that the length of the output parameter derivatives list equals the+   length of the input parameters list.++For example, the jacobian of the above @quad@ model can be written as:++@+quadJacob :: 'Num' r => 'Jacobian' N3 r r+quadJacob [_, _, _] x = [ x^2   -- with respect to a+                        , x     -- with respect to b+                        , 1     -- with respect to c+                        ]+@++(Notice you don't have to differentiate for @x@.)+-}+type Jacobian r a = [r] → (a → [r])++-- | This type synonym expresses that usually the @a@ in @'Jacobian' r a@+-- equals the type of the parameters.+type SimpleJacobian r = Jacobian r r+++--------------------------------------------------------------------------------+-- Levenberg-Marquardt algorithm.+--------------------------------------------------------------------------------++-- | The Levenberg-Marquardt algorithm specialised for curve-fitting.+levmar ∷ LevMar.LevMarable r+       ⇒ Model r a                          -- ^ Model+       → Maybe (Jacobian r a)               -- ^ Optional jacobian+       → [r]                                -- ^ Initial parameters+       → [(a, r)]                           -- ^ Samples+       → Integer                            -- ^ Maximum iterations+       → LevMar.Options r                   -- ^ Minimization options+       → LevMar.Constraints r               -- ^ Constraints+       → Either LevMar.LevMarError ([r], LevMar.Info r, LevMar.CovarMatrix r)+levmar model mJac params samples =+    LevMar.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 ---------------------------------------------------------------------

@@ -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 ---------------------------------------------------------------------

@@ -1,103 +0,0 @@--- Thanks to Ryan Ingram who wrote most of this module.--- See: http://www.haskell.org/pipermail/haskell-cafe/2009-August/065674.html--{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE RankNTypes #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeFamilies #-}--module TypeLevelNat-    ( Z(..)-    , S(..)-    , Nat-    , caseNat-    , induction-    , witnessNat--    , N(..)-    , nat-    ) where----- | Type-level natural denoting zero-data Z = Z deriving Show---- | Type-level natural denoting the /S/uccessor of another type-level natural.-newtype S n = S n deriving Show---- | Class of all type-level naturals.-class Nat n where-   -- | Case analysis on natural numbers.-   caseNat :: forall r.-              n                                      -- ^ The natural number to case analyse.-           -> (n ~ Z => r)                           -- ^ The result @r@ when @n@ equals zero.-           -> (forall p. (n ~ S p, Nat p) => p -> r) -- ^ Function to apply to the predecessor-                                                     --   of @n@ to yield the result @r@.-           -> r--instance Nat Z where-   caseNat _ z _ = z--instance Nat p => Nat (S p) where-   caseNat (S p) _ s = s p---- | The axiom of induction on natural numbers.--- See: <http://en.wikipedia.org/wiki/Mathematical_induction#Axiom_of_induction>-induction :: forall p n. Nat n-          => n-          -> p Z-          -> (forall m. Nat m => p m -> p (S m))-          -> p n-induction n z s = caseNat n isZ isS-    where-      isZ :: n ~ Z => p n-      isZ = z--      isS :: forall m. (n ~ S m, Nat m) => m -> p n-      isS m = s (induction m z s)--newtype Witness x = Witness { unWitness :: x }---- | The value of @witnessNat :: n@ is the natural number of type @n@.--- For example:------ @--- *TypeLevelNat> witnessNat :: S (S (S Z))--- S (S (S Z))--- @-witnessNat :: forall n. Nat n => n-witnessNat = theWitness-    where-      theWitness = unWitness $ induction (undefined `asTypeOf` theWitness)-                                         (Witness Z)-                                         (Witness . S . unWitness)---- | A value-level natural indexed with an equivalent type-level natural.-data N n where-    Zero :: N Z-    Succ :: N n -> N (S n)--nat :: forall n. Nat n => n -> N n-nat n = induction n Zero Succ--{--Template Haskell code to construct a type synonym for an arbitrary-type level natural number.--Instead of--> type N6 = S (S (S (S (S (S Z)))))--you can write--> $(mkNat "N6" 6)--}---- import Language.Haskell.TH.Syntax---- mkNat :: String -> Int -> Q [Dec]--- mkNat syn = runQ . return . (:[]) . TySynD (mkName syn) [] . go---     where go 0 = ConT $ mkName "Z"---           go n = AppT (ConT $ mkName "S") $ go (n - 1)-

@@ -1,102 +1,56 @@ name:          levmar-version:       0.2.1+version:       0.3 cabal-version: >= 1.6 build-type:    Simple stability:     experimental-tested-with:   GHC ==6.10.4-author:        Roel van Dijk & Bas van Dijk-maintainer:    vandijk.roel@gmail.com, v.dijk.bas@gmail.com-copyright:     (c) 2009 Roel van Dijk & Bas van Dijk+author:        Roel van Dijk <vandijk.roel@gmail.com>+               Bas van Dijk <v.dijk.bas@gmail.com>+maintainer:    Roel van Dijk <vandijk.roel@gmail.com>+               Bas van Dijk <v.dijk.bas@gmail.com>+copyright:     (c) 2009 - 2010 Roel van Dijk & Bas van Dijk license:       BSD3 license-file:  LICENSE category:      Numerical, Math synopsis:      An implementation of the Levenberg-Marquardt algorithm-description:   The Levenberg-Marquardt algorithm is an iterative-               technique that finds a local minimum of a function that-               is expressed as the sum of squares of nonlinear-               functions. It has become a standard technique for-               nonlinear least-squares problems and can be thought of-               as a combination of steepest descent and the-               Gauss-Newton method. When the current solution is far-               from the correct one, the algorithm behaves like a-               steepest descent method: slow, but guaranteed to-               converge. When the current solution is close to the-               correct solution, it becomes a Gauss-Newton method.-               .-               Optional box- and linear constraints can be given. Both-               single and double precision floating point types are-               supported.-               .-               The actual algorithm is implemented in a C library-               which is bundled with bindings-levmar which this-               package depends on. See:-               <http://www.ics.forth.gr/~lourakis/levmar/>.-               .-               This library consists of two layers:-               .-               * LevMar.Intermediate: A medium-level layer that wraps-                 the low-level functions from bindings-levmar to-                 provide a more Haskell friendly interface.-               .-	       * LevMar: A high-level layer that uses type-level-                 programming to add extra type safety.-               .-               Each layer also has special curve-fitting variants:-               .-	       * LevMar.Intermediate.Fitting-               .-               * LevMar.Fitting-               .-               Each layer also has special variants that automatically compute-               the jacobian using automatic differentiation using Conal-               Elliott's vector-space library:-               .-               * LevMar.Intermediate.AD-               .-	       * LevMar.Intermediate.Fitting.AD-               .-               * LevMar.AD-               .-	       * LevMar.Fitting.AD-	       .-               Note however that this feature is still very experimental!-               .-	       All modules are self-contained; i.e. each module-	       re-exports all the things you need to work with it.-	       .-	       For an example how to use this library see Demo.hs-	       which is included in this package. Demo.hs is a Haskell-	       translation of lmdemo.c from the C levmar library.-	       .-               A note regarding the license:-               .-               This library depends on bindings-levmar which is-               bundled together with a C library which falls under the-               GPL. Please be aware of this when distributing programs-               linked with this library. For details see the-               description and license of bindings-levmar.-extra-source-files: Demo.hs+description:+  The Levenberg-Marquardt algorithm is an iterative technique that+  finds a local minimum of a function that is expressed as the sum of+  squares of nonlinear functions. It has become a standard technique+  for nonlinear least-squares problems and can be thought of as a+  combination of steepest descent and the Gauss-Newton method. When+  the current solution is far from the correct one, the algorithm+  behaves like a steepest descent method: slow, but guaranteed to+  converge. When the current solution is close to the correct+  solution, it becomes a Gauss-Newton method.+  .+  Optional box- and linear constraints can be given. Both single and+  double precision floating point types are supported.+  .+  The actual algorithm is implemented in a C library which is bundled+  with bindings-levmar which this package depends on. See:+  <http://www.ics.forth.gr/~lourakis/levmar/>.+  .+  All modules are self-contained; i.e. each module re-exports all the+  things you need to work with it.+  .+  Also see the @levmar-safe@ package which adds extra type-safety on+  top of this package.+  .+  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.  source-repository head   Type: darcs   Location: http://code.haskell.org/levmar  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     && < 4.4+               , base-unicode-symbols >= 0.1.1 && < 0.3+               , bindings-levmar      >= 0.2   && < 0.3+  exposed-modules: Numeric.LevMar, Numeric.LevMar.Fitting+  other-modules: Bindings.LevMar.CurryFriendly+  ghc-options: -Wall