levmar 0.3 → 1.0
raw patch · 6 files changed
+340/−413 lines, 6 filesdep +hmatrixdep +vectordep ~bindings-levmar
Dependencies added: hmatrix, vector
Dependency ranges changed: bindings-levmar
Files
- Bindings/LevMar/CurryFriendly.hs +40/−41
- LICENSE +1/−1
- Numeric/LevMar.hs +258/−214
- Numeric/LevMar/Fitting.hs +0/−142
- README.markdown +26/−0
- levmar.cabal +15/−15
Bindings/LevMar/CurryFriendly.hs view
@@ -27,9 +27,8 @@ , dlevmar_blec_dif, slevmar_blec_dif ) where --import Foreign.C.Types ( CFloat, CDouble )-import Foreign.Ptr ( FunPtr )+import Prelude ( Double, Float )+import Foreign.Ptr ( FunPtr ) import qualified Bindings.LevMar as BLM @@ -38,59 +37,59 @@ -- Handy type synonyms used in the curry friendly types. -------------------------------------------------------------------------------- -type BoxConstraints cr a = BLM.LowerBounds cr- → BLM.UpperBounds cr- → a+type BoxConstraints r α = BLM.LowerBounds r+ → BLM.UpperBounds r+ → α -type LinearConstraints cr a = BLM.ConstraintsMatrix cr- → BLM.ConstraintsVector cr- → BLM.NrOfConstraints- → a+type LinearConstraints r α = BLM.ConstraintsMatrix r+ → BLM.ConstraintsVector r+ → BLM.NrOfConstraints+ → α -------------------------------------------------------------------------------- -- 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))+type LevMarDif r = BLM.LevMarDif r+type LevMarDer r = FunPtr (BLM.Jacobian r) → LevMarDif r+type LevMarBCDif r = BoxConstraints r (LevMarDif r)+type LevMarBCDer r = BoxConstraints r (LevMarDer r)+type LevMarLecDif r = LinearConstraints r (LevMarDif r)+type LevMarLecDer r = LinearConstraints r (LevMarDer r)+type LevMarBLecDif r = BoxConstraints r (LinearConstraints r (BLM.Weights r → LevMarDif r))+type LevMarBLecDer r = BoxConstraints r (LinearConstraints r (BLM.Weights r → LevMarDer r)) -------------------------------------------------------------------------------- -- Reordering arguments to create curry friendly variants. -------------------------------------------------------------------------------- -mk_levmar_der ∷ BLM.LevMarDer cr → LevMarDer cr+mk_levmar_der ∷ BLM.LevMarDer r → LevMarDer r mk_levmar_der lma j f = lma f j -mk_levmar_bc_dif ∷ BLM.LevMarBCDif cr → LevMarBCDif cr+mk_levmar_bc_dif ∷ BLM.LevMarBCDif r → LevMarBCDif r mk_levmar_bc_dif lma lb ub f p x m n = lma f p x m n lb ub -mk_levmar_bc_der ∷ BLM.LevMarBCDer cr → LevMarBCDer cr+mk_levmar_bc_der ∷ BLM.LevMarBCDer r → LevMarBCDer r mk_levmar_bc_der lma lb ub j f p x m n = lma f j p x m n lb ub -mk_levmar_lec_dif ∷ BLM.LevMarLecDif cr → LevMarLecDif cr+mk_levmar_lec_dif ∷ BLM.LevMarLecDif r → LevMarLecDif r mk_levmar_lec_dif lma a b k f p x m n = lma f p x m n a b k -mk_levmar_lec_der ∷ BLM.LevMarLecDer cr → LevMarLecDer cr+mk_levmar_lec_der ∷ BLM.LevMarLecDer r → LevMarLecDer r mk_levmar_lec_der lma a b k j f p x m n = lma f j p x m n a b k -mk_levmar_blec_dif ∷ BLM.LevMarBLecDif cr → LevMarBLecDif cr+mk_levmar_blec_dif ∷ BLM.LevMarBLecDif r → LevMarBLecDif r mk_levmar_blec_dif lma lb ub a b k wghts f p x m n = lma f p x m n lb ub a b k wghts -mk_levmar_blec_der ∷ BLM.LevMarBLecDer cr → LevMarBLecDer cr+mk_levmar_blec_der ∷ BLM.LevMarBLecDer r → LevMarBLecDer r mk_levmar_blec_der lma lb ub a b k wghts j f p x m n = lma f j p x m n lb ub a b k wghts @@ -100,52 +99,52 @@ -- 'Bindings.Levmar'. -------------------------------------------------------------------------------- -slevmar_dif ∷ LevMarDif CFloat+slevmar_dif ∷ LevMarDif Float slevmar_dif = BLM.c'slevmar_dif -dlevmar_dif ∷ LevMarDif CDouble+dlevmar_dif ∷ LevMarDif Double dlevmar_dif = BLM.c'dlevmar_dif -slevmar_der ∷ LevMarDer CFloat+slevmar_der ∷ LevMarDer Float slevmar_der = mk_levmar_der BLM.c'slevmar_der -dlevmar_der ∷ LevMarDer CDouble+dlevmar_der ∷ LevMarDer Double dlevmar_der = mk_levmar_der BLM.c'dlevmar_der -slevmar_bc_dif ∷ LevMarBCDif CFloat+slevmar_bc_dif ∷ LevMarBCDif Float slevmar_bc_dif = mk_levmar_bc_dif BLM.c'slevmar_bc_dif -dlevmar_bc_dif ∷ LevMarBCDif CDouble+dlevmar_bc_dif ∷ LevMarBCDif Double dlevmar_bc_dif = mk_levmar_bc_dif BLM.c'dlevmar_bc_dif -slevmar_bc_der ∷ LevMarBCDer CFloat+slevmar_bc_der ∷ LevMarBCDer Float slevmar_bc_der = mk_levmar_bc_der BLM.c'slevmar_bc_der -dlevmar_bc_der ∷ LevMarBCDer CDouble+dlevmar_bc_der ∷ LevMarBCDer Double dlevmar_bc_der = mk_levmar_bc_der BLM.c'dlevmar_bc_der -slevmar_lec_dif ∷ LevMarLecDif CFloat+slevmar_lec_dif ∷ LevMarLecDif Float slevmar_lec_dif = mk_levmar_lec_dif BLM.c'slevmar_lec_dif -dlevmar_lec_dif ∷ LevMarLecDif CDouble+dlevmar_lec_dif ∷ LevMarLecDif Double dlevmar_lec_dif = mk_levmar_lec_dif BLM.c'dlevmar_lec_dif -slevmar_lec_der ∷ LevMarLecDer CFloat+slevmar_lec_der ∷ LevMarLecDer Float slevmar_lec_der = mk_levmar_lec_der BLM.c'slevmar_lec_der -dlevmar_lec_der ∷ LevMarLecDer CDouble+dlevmar_lec_der ∷ LevMarLecDer Double dlevmar_lec_der = mk_levmar_lec_der BLM.c'dlevmar_lec_der -slevmar_blec_dif ∷ LevMarBLecDif CFloat+slevmar_blec_dif ∷ LevMarBLecDif Float slevmar_blec_dif = mk_levmar_blec_dif BLM.c'slevmar_blec_dif -dlevmar_blec_dif ∷ LevMarBLecDif CDouble+dlevmar_blec_dif ∷ LevMarBLecDif Double dlevmar_blec_dif = mk_levmar_blec_dif BLM.c'dlevmar_blec_dif -slevmar_blec_der ∷ LevMarBLecDer CFloat+slevmar_blec_der ∷ LevMarBLecDer Float slevmar_blec_der = mk_levmar_blec_der BLM.c'slevmar_blec_der -dlevmar_blec_der ∷ LevMarBLecDer CDouble+dlevmar_blec_der ∷ LevMarBLecDer Double dlevmar_blec_der = mk_levmar_blec_der BLM.c'dlevmar_blec_der
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2009 Roel van Dijk, Bas van Dijk+Copyright (c) 2009-2011 Roel van Dijk, Bas van Dijk All rights reserved.
Numeric/LevMar.hs view
@@ -8,7 +8,7 @@ -------------------------------------------------------------------------------- -- | -- Module: Numeric.LevMar--- Copyright: (c) 2009 - 2010 Roel van Dijk & Bas van Dijk+-- Copyright: (c) 2009 - 2011 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>@@ -28,21 +28,17 @@ -- * Levenberg-Marquardt algorithm. , LevMarable(levmar) - , LinearConstraints- -- * Minimization options. , Options(..) , defaultOpts -- * Constraints , Constraints(..)- , noConstraints+ , LinearConstraints -- * Output , Info(..) , StopReason(..)- , CovarMatrix- , LevMarError(..) ) where @@ -52,41 +48,56 @@ ------------------------------------------------------------------------------- -- from base:-import Control.Monad.Instances -- for 'instance Functor (Either a)'+import Control.Monad ( return, mplus ) 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.Functor ( (<$>) )+import Data.Int ( Int )+import Data.List ( lookup, (++) )+import Data.Maybe ( Maybe(Nothing, Just), isJust, fromJust, fromMaybe )+import Data.Monoid ( Monoid(mempty, mappend) ) import Data.Ord ( (<) )-import Foreign.Marshal.Array ( allocaArray, peekArray, pokeArray, withArray )-import Foreign.Ptr ( Ptr, nullPtr, plusPtr )+import Foreign.Marshal.Array ( allocaArray, withArray, peekArray, copyArray )+import Foreign.Ptr ( Ptr, nullPtr )+import Foreign.ForeignPtr ( ForeignPtr, newForeignPtr_, withForeignPtr ) import Foreign.Storable ( Storable )-import Foreign.C.Types ( CInt )-import Prelude ( Enum, Fractional, Real, RealFrac- , Integer, Float, Double- , fromIntegral, realToFrac, toEnum- , (-), error, floor+import Prelude ( Enum, Fractional, RealFrac, Float, Double+ , toEnum, (-), (*), error, floor ) import System.IO ( IO ) import System.IO.Unsafe ( unsafePerformIO ) import Text.Read ( Read )-import Text.Show ( Show )+import Text.Show ( Show, show ) +#if __GLASGOW_HASKELL__ >= 605+import GHC.ForeignPtr ( mallocPlainForeignPtrBytes )+import Prelude ( undefined )+import Foreign.Storable ( sizeOf )+#else+import Foreign.ForeignPtr ( mallocForeignPtrArray )+#endif+ #if __GLASGOW_HASKELL__ < 700 import Prelude ( fromInteger ) #endif -- from base-unicode-symbols: import Data.Bool.Unicode ( (∧), (∨) )-import Data.Eq.Unicode ( (≡), (≢) )+import Data.Eq.Unicode ( (≢) ) import Data.Function.Unicode ( (∘) )-import Prelude.Unicode ( (⋅) ) +-- from hmatrix:+import Data.Packed.Vector ( Vector )+import Data.Packed.Matrix ( Matrix, Element, flatten, rows, reshape )++-- from vector:+import qualified Data.Vector.Storable as VS ( unsafeWith, length+ , unsafeFromForeignPtr+ , length+ )+ -- from bindings-levmar: import Bindings.LevMar ( c'LM_INFO_SZ @@ -134,51 +145,29 @@ -------------------------------------------------------------------------------- {-| 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'.+from a vector of parameters to a vector of expected measurements. -For example:+ * Ensure that the length of the parameters vector equals the length of the+ initial parameters vector in 'levmar'. -@-hatfldc :: Model Double-hatfldc [p0, p1, p2, p3] = [ p0 - 1.0- , p0 - sqrt p1- , p1 - sqrt p2- , p3 - 1.0- ]-@+ * Ensure that the length of the ouput vector equals the length of the samples+ vector in 'levmar'. -}-type Model r = [r] → [r]+type Model r = Vector r → Vector 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 jacobian of the 'Model' function. Expressed as a function from a vector+of parameters to a matrix 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 length of the parameter vector equals the length of the initial+ parameter vector in 'levmar'. - * Ensure that the output matrix has the dimension @n@/x/@m@ where @n@ is the+ * Ensure that the output matrix has the dimension @n><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]]+type Jacobian r = Vector r → Matrix r --------------------------------------------------------------------------------@@ -189,14 +178,18 @@ class LevMarable r where -- | The Levenberg-Marquardt algorithm.+ --+ -- Returns a tuple of the found parameters, a structure containing+ -- information about the minimization and the covariance matrix+ -- corresponding to LS solution. levmar ∷ Model r -- ^ Model → Maybe (Jacobian r) -- ^ Optional jacobian- → [r] -- ^ Initial parameters- → [r] -- ^ Samples- → Integer -- ^ Maximum iterations+ → Vector r -- ^ Initial parameters+ → Vector r -- ^ Samples+ → Int -- ^ Maximum iterations → Options r -- ^ Minimization options → Constraints r -- ^ Constraints- → Either LevMarError ([r], Info r, CovarMatrix r)+ → Either LevMarError (Vector r, Info r, Matrix r) instance LevMarable Float where levmar = gen_levmar slevmar_der@@ -232,24 +225,24 @@ 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+gen_levmar ∷ ∀ r. (Storable r, RealFrac r, Element r)+ ⇒ LevMarDer r+ → LevMarDif r+ → LevMarBCDer r+ → LevMarBCDif r+ → LevMarLecDer r+ → LevMarLecDif r+ → LevMarBLecDer r+ → LevMarBLecDif r → Model r -- ^ Model → Maybe (Jacobian r) -- ^ Optional jacobian- → [r] -- ^ Initial parameters- → [r] -- ^ Samples- → Integer -- ^ Maximum iterations+ → Vector r -- ^ Initial parameters+ → Vector r -- ^ Samples+ → Int -- ^ Maximum iterations → Options r -- ^ Options → Constraints r -- ^ Constraints- → Either LevMarError ([r], Info r, CovarMatrix r)+ → Either LevMarError (Vector r, Info r, Matrix r) gen_levmar f_der f_dif f_bc_der@@ -258,135 +251,180 @@ f_lec_dif f_blec_der f_blec_dif- model mJac ps ys itMax opts (Constraints mLowBs mUpBs mWeights mLinC)- = unsafePerformIO ∘+ model mJac ps ys itMax opts (Constraints mLowBs mUpBs mWeights mLinC) =+ -- All effects are contained, so we can safely perform:+ unsafePerformIO $ do - -- 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+ -- We need to pass the initial parameters 'ps' to the C function.+ -- However, we can't just pass a pointer to them because the C function+ -- will modify the parameters during execution which will violate+ -- referential transparanency. Instead we allocate new space+ -- and copy the parameters to it.+ --+ -- Note that, in the end, the array is returned from this function.+ -- This means that the only way to guarantee its finalisation+ -- is to allocate it using a ForeignPtr:+ psFP ← fastMallocForeignPtrArray lenPs+ withForeignPtr psFP $ \psPtr → do+ VS.unsafeWith ps $ \psPtrInp →+ copyArray psPtr psPtrInp lenPs - -- 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+ -- Retrieve the (read-only) pointer 'ysPtr' to the samples vector 'ys'+ -- so we can pass it to the C function:+ VS.unsafeWith ys $ \ysPtr → - 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+ -- Convert the Options 'opts' to a list and then to an array+ -- so we can pass the (read-only) pointer 'optsPtr' to the C function:+ withArray (optsToList opts) $ \optsPtr → - 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+ -- Allocate space for the info array+ -- so we can pass it to the C function.+ -- Note that, in the end, this array is converted to an Info value+ -- and returned from this function.+ allocaArray c'LM_INFO_SZ $ \infoPtr → do - -- 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+ -- Allocate space for the covariance matrix+ -- so we can pass it to the C function.+ -- Like the parameters array the matrix+ -- needs to be returned from this function.+ -- So we also allocate it using a ForeignPtr:+ covarFP ← fastMallocForeignPtrArray covarLen+ withForeignPtr covarFP $ \covarPtr → - let convertCovarMatrix ptr- | ptr ≡ covarPtr `plusPtr` covarLen = return []- | otherwise = do row ← peekArray lenPs ptr- rows ← convertCovarMatrix $ ptr `plusPtr` lenPs- return $ map realToFrac row : rows+ -- 'cmodel' is the low-level model function which is converted+ -- to the FunPtr 'modelFunPtr' and passed to the C function.+ -- 'cmodel' will first convert the parameters pointer 'parPtr'+ -- into a Vector after converting it into a ForeignPtr+ -- (without a finalizer).+ -- Then it will apply the high-level 'model' function+ -- to this parameter vector. The resulting vector is then copied+ -- to the output buffer 'hxPtr':+ let cmodel ∷ Bindings.LevMar.Model r+ cmodel parPtr hxPtr _ _ _ = do+ parFP ← newForeignPtr_ parPtr+ let psV = VS.unsafeFromForeignPtr parFP 0 lenPs+ vector = model psV+ VS.unsafeWith vector $ \p → copyArray hxPtr p (VS.length vector)+ in withModel cmodel $ \modelFunPtr → do - 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+ -- All the low-level C functions share a common set of arguments.+ -- 'runDif' applies these arguments to the given C function 'f':+ let runDif ∷ LevMarDif r → IO Int+ runDif f = f modelFunPtr+ psPtr+ ysPtr+ lenPs+ lenYs+ itMax+ optsPtr+ infoPtr+ nullPtr+ covarPtr+ nullPtr - -- Whether the parameters are constrained by a linear equation.- linConstrained = isJust mLinC+ err ← case mJac of+ Nothing → if boxConstrained+ then if linConstrained+ then withBoxConstraints+ (withLinConstraints $ withWeights runDif)+ f_blec_dif+ else withBoxConstraints runDif f_bc_dif+ else if linConstrained+ then withLinConstraints runDif f_lec_dif+ else runDif f_dif - -- Whether the parameters are constrained by a bounding box.- boxConstrained = isJust mLowBs ∨ isJust mUpBs+ Just jac →+ let cjacobian ∷ Bindings.LevMar.Jacobian r+ cjacobian parPtr jPtr _ _ _ = do+ parFP ← newForeignPtr_ parPtr+ let psV = VS.unsafeFromForeignPtr parFP 0 lenPs+ matrix = jac psV+ vector = flatten matrix+ VS.unsafeWith vector $ \p →+ copyArray jPtr p (VS.length vector)+ in withJacobian cjacobian $ \jacobPtr → - withBoxConstraints f g =- maybeWithArray mLowBs $ \lBsPtr →- maybeWithArray mUpBs $ \uBsPtr →- f $ g lBsPtr uBsPtr+ let runDer ∷ LevMarDer r → IO Int+ 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 - withLinConstraints f g =- withArray (map realToFrac $ concat cMat) $ \cMatPtr →- withArray (map realToFrac rhcVec) $ \rhcVecPtr →- f ∘ g cMatPtr rhcVecPtr ∘ fromIntegral $ length cMat+ -- Handling errors:+ if err < 0+ -- we don't treat the following two as an error:+ ∧ err ≢ c'LM_ERROR_SINGULAR_MATRIX+ ∧ err ≢ c'LM_ERROR_SUM_OF_SQUARES_NOT_FINITE+ then return $ Left $ convertLevMarError err - withWeights f g = maybeWithArray mWeights $ f ∘ g+ else do -- Converting results:+ info ← listToInfo <$> peekArray c'LM_INFO_SZ infoPtr+ let psV = VS.unsafeFromForeignPtr psFP 0 lenPs+ let covarM = reshape lenPs $+ VS.unsafeFromForeignPtr covarFP 0 covarLen -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+ return $ Right (psV, info, covarM)+ where+ lenPs = VS.length ps+ lenYs = VS.length ys+ covarLen = lenPs*lenPs+ (cMat, rhcVec) = fromJust mLinC -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+ -- Whether the parameters are constrained by a linear equation.+ linConstrained = isJust mLinC --- | 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])+ -- Whether the parameters are constrained by a bounding box.+ boxConstrained = isJust mLowBs ∨ isJust mUpBs + withBoxConstraints f g =+ maybeWithArray mLowBs $ \lBsPtr →+ maybeWithArray mUpBs $ \uBsPtr →+ f $ g lBsPtr uBsPtr + withLinConstraints f g =+ VS.unsafeWith (flatten cMat) $ \cMatPtr →+ VS.unsafeWith rhcVec $ \rhcVecPtr →+ f ∘ g cMatPtr rhcVecPtr $ rows cMat++ withWeights f g = maybeWithArray mWeights $ f ∘ g++maybeWithArray ∷ (Storable α) ⇒ Maybe (Vector α) → (Ptr α → IO β) → IO β+maybeWithArray Nothing f = f nullPtr+maybeWithArray (Just v) f = VS.unsafeWith v f++#if __GLASGOW_HASKELL__ >= 605+{-# INLINE fastMallocForeignPtrArray #-}+fastMallocForeignPtrArray ∷ ∀ α. Storable α ⇒ Int → IO (ForeignPtr α)+fastMallocForeignPtrArray n = mallocPlainForeignPtrBytes+ (n * sizeOf (undefined ∷ α))+#else+fastMallocForeignPtrArray ∷ ∀ α. Storable α ⇒ Int → IO (ForeignPtr α)+fastMallocForeignPtrArray = mallocForeignPtrArray+#endif++ -------------------------------------------------------------------------------- -- Minimization options. -------------------------------------------------------------------------------- -- | Minimization options data Options r =- Opts { optScaleInitMu ∷ r -- ^ Scale factor for initial /mu/.- , optStopNormInfJacTe ∷ r -- ^ Stopping thresholds for @||J^T e||_inf@.- , optStopNorm2Dp ∷ r -- ^ Stopping thresholds for @||Dp||_2@.- , optStopNorm2E ∷ r -- ^ Stopping thresholds for @||e||_2@.- , optDelta ∷ r -- ^ Step used in the difference- -- approximation to the Jacobian. If- -- @optDelta<0@, the Jacobian is approximated- -- with central differences which are more- -- accurate (but slower!) compared to the- -- forward differences employed by default.+ 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@@ -407,21 +445,30 @@ -- Constraints -------------------------------------------------------------------------------- +-- | Ensure that these vectors have the same length as the number of parameters. 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+ { lowerBounds ∷ !(Maybe (Vector r)) -- ^ Optional lower bounds+ , upperBounds ∷ !(Maybe (Vector r)) -- ^ Optional upper bounds+ , weights ∷ !(Maybe (Vector 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+-- | Linear constraints consisting of a constraints matrix, @k><m@ and+-- a right hand constraints vector, of length @k@ where @m@ is the number of+-- parameters and @k@ is the number of constraints.+type LinearConstraints r = (Matrix r, Vector r) -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+-- | * 'mempty' is defined as a 'Constraints' where all fields are 'Nothing'.+--+-- * 'mappend' merges two 'Constraints' by taking the first non-'Nothing' value+-- for each field.+instance Monoid (Constraints r) where+ mempty = Constraints Nothing Nothing Nothing Nothing+ mappend (Constraints lb1 ub1 w1 l1)+ (Constraints lb2 ub2 w2 l2) = Constraints (lb1 `mplus` lb2)+ (ub1 `mplus` ub2)+ (w1 `mplus` w2)+ (l1 `mplus` l2) --------------------------------------------------------------------------------@@ -430,26 +477,26 @@ -- | 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.+ { infNorm2initE ∷ !r -- ^ @||e||_2@ at initial parameters.+ , infNorm2E ∷ !r -- ^ @||e||_2@ at estimated parameters.+ , infNormInfJacTe ∷ !r -- ^ @||J^T e||_inf@ at estimated parameters.+ , infNorm2Dp ∷ !r -- ^ @||Dp||_2@ at estimated parameters.+ , infMuDivMax ∷ !r -- ^ @\mu/max[J^T J]_ii ]@ at estimated parameters.+ , infNumIter ∷ !Int -- ^ Number of iterations.+ , infStopReason ∷ !StopReason -- ^ Reason for terminating.+ , infNumFuncEvals ∷ !Int -- ^ Number of function evaluations.+ , infNumJacobEvals ∷ !Int -- ^ Number of jacobian evaluations.+ , infNumLinSysSolved ∷ !Int -- ^ Number of linear systems solved,+ -- i.e. attempts for reducing error. } deriving (Read, Show) -listToInfo ∷ (RealFrac cr, Fractional r) ⇒ [cr] → Info r+listToInfo ∷ (RealFrac r) ⇒ [r] → 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+ Info { infNorm2initE = a+ , infNorm2E = b+ , infNormInfJacTe = c+ , infNorm2Dp = d+ , infMuDivMax = e , infNumIter = floor f , infStopReason = toEnum $ floor g - 1 , infNumFuncEvals = floor h@@ -472,10 +519,7 @@ -- (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 --------------------------------------------------------------------------------@@ -504,7 +548,7 @@ -- Handy in case you want to thow a LevMarError as an exception: instance Exception LevMarError -levmarCErrorToLevMarError ∷ [(CInt, LevMarError)]+levmarCErrorToLevMarError ∷ [(Int, LevMarError)] levmarCErrorToLevMarError = [ (c'LM_ERROR, LevMarError) , (c'LM_ERROR_LAPACK_ERROR, LapackError)@@ -519,9 +563,9 @@ --, (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+convertLevMarError ∷ Int → LevMarError+convertLevMarError err = fromMaybe (error $ "Unknown levmar error: " ++ show err)+ (lookup err levmarCErrorToLevMarError) -- The End ---------------------------------------------------------------------
− Numeric/LevMar/Fitting.hs
@@ -1,142 +0,0 @@-{-# 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 ---------------------------------------------------------------------
+ README.markdown view
@@ -0,0 +1,26 @@+The Levenberg-Marquardt algorithm is an iterative technique that+finds a local minimum of a function that is expressed as the sum of+squares of nonlinear functions. It has become a standard technique+for nonlinear least-squares problems and can be thought of as a+combination of steepest descent and the Gauss-Newton method. When+the current solution is far from the correct one, the algorithm+behaves like a steepest descent method: slow, but guaranteed to+converge. When the current solution is close to the correct+solution, it becomes a Gauss-Newton method.++Optional box- and linear constraints can be given. Both single and+double precision floating point types are supported.++The actual algorithm is implemented in a [C library] which is bundled+with [bindings-levmar] which this package depends on.++License+=======++This library depends on [bindings-levmar] which is bundled together+with a [C library] which falls under the GPL. Please be aware of this+when distributing programs linked with this library. For details see+the description and license of [bindings-levmar].++[bindings-levmar]: http://hackage.haskell.org/package/bindings-levmar+[C library]: http://www.ics.forth.gr/~lourakis/levmar
levmar.cabal view
@@ -1,5 +1,5 @@ name: levmar-version: 0.3+version: 1.0 cabal-version: >= 1.6 build-type: Simple stability: experimental@@ -7,9 +7,11 @@ 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+copyright: (c) 2009 - 2011 Roel van Dijk & Bas van Dijk license: BSD3 license-file: LICENSE+homepage: https://github.com/basvandijk/levmar/+bug-reports: https://github.com/basvandijk/levmar/issues category: Numerical, Math synopsis: An implementation of the Levenberg-Marquardt algorithm description:@@ -27,30 +29,28 @@ 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:+ 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+ 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.+ the description and license of @bindings-levmar@. +extra-source-files: README.markdown+ source-repository head- Type: darcs- Location: http://code.haskell.org/levmar+ Type: git+ Location: git://github.com/basvandijk/levmar.git library build-depends: base >= 3 && < 4.4 , base-unicode-symbols >= 0.1.1 && < 0.3- , bindings-levmar >= 0.2 && < 0.3- exposed-modules: Numeric.LevMar, Numeric.LevMar.Fitting+ , bindings-levmar >= 1.0 && < 1.1+ , hmatrix >= 0.11 && < 0.12+ , vector >= 0.7 && < 0.8+ exposed-modules: Numeric.LevMar other-modules: Bindings.LevMar.CurryFriendly ghc-options: -Wall