diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2009-2011 Roel van Dijk, Bas van Dijk +Copyright (c) 2009-2012 Roel van Dijk, Bas van Dijk All rights reserved. diff --git a/Numeric/LevMar.hs b/Numeric/LevMar.hs --- a/Numeric/LevMar.hs +++ b/Numeric/LevMar.hs @@ -11,7 +11,7 @@ -------------------------------------------------------------------------------- -- | -- Module: Numeric.LevMar --- Copyright: (c) 2009 - 2011 Roel van Dijk & Bas van Dijk +-- Copyright: (c) 2009 - 2012 Roel van Dijk & Bas van Dijk -- License: BSD-style (see the file LICENSE) -- Maintainer: Roel van Dijk -- Bas van Dijk @@ -25,24 +25,20 @@ module Numeric.LevMar ( -- * Model & Jacobian. - Model - , Jacobian + Params, Samples + , Model, Jacobian -- * Levenberg-Marquardt algorithm. , LevMarable(levmar) -- * Minimization options. - , Options(..) - , defaultOpts + , Options(..), defaultOpts -- * Constraints - , Constraints(..) - , LinearConstraints + , Constraints(..), LinearConstraints -- * Output - , Info(..) - , StopReason(..) - , LevMarError(..) + , Info(..), StopReason(..), LevMarError(..) ) where @@ -53,10 +49,11 @@ -- from base: import Control.Monad ( return, mplus ) import Control.Exception ( Exception ) +import Data.Bool ( otherwise ) import Data.Data ( Data ) import Data.Typeable ( Typeable ) import Data.Either ( Either(Left, Right) ) -import Data.Eq ( Eq ) +import Data.Eq ( Eq, (==) ) import Data.Function ( ($) ) import Data.Functor ( (<$>) ) import Data.Int ( Int ) @@ -64,12 +61,13 @@ import Data.Maybe ( Maybe(Nothing, Just), isJust, fromJust, fromMaybe ) import Data.Monoid ( Monoid, mempty, mappend ) import Data.Ord ( Ord, (<) ) +import Foreign.C.Types ( CInt ) import Foreign.Marshal.Array ( allocaArray, withArray, peekArray, copyArray ) import Foreign.Ptr ( Ptr, nullPtr ) import Foreign.ForeignPtr ( ForeignPtr, newForeignPtr_, withForeignPtr ) import Foreign.Storable ( Storable ) import Prelude ( Enum, Fractional, RealFrac, Float, Double - , toEnum, (-), (*), error, floor + , fromIntegral, toEnum, (-), (*), error, floor ) import System.IO ( IO ) import System.IO.Unsafe ( unsafePerformIO ) @@ -148,31 +146,43 @@ -- Model & Jacobian. -------------------------------------------------------------------------------- +-- | Parameter vector of length @m@. +-- +-- Ensure that @m <= n@ where @n@ is the length of the 'Samples' vector. +type Params r = Vector r + +-- | Sample vector of length @n@. +-- +-- Ensure that @n >= m@ where @m@ is the length of the 'Params' vector. +type Samples r = Vector r + {-| A functional relation describing measurements represented as a function -from a vector of parameters to a vector of expected measurements. +from a vector of parameters to a vector of expected samples. - * Ensure that the length of the parameters vector equals the length of the - initial parameters vector in 'levmar'. + * Ensure that the length @m@ of the parameter vector equals the length of the + initial parameter vector in 'levmar'. - * Ensure that the length of the ouput vector equals the length of the samples - vector in 'levmar'. + * Ensure that the length @n@ of the output sample vector equals the length of + the sample vector in 'levmar'. + + * Ensure that the length @n@ of the output sample vector vector is bigger than or + equal to the length @m@ of the parameter vector. -} -type Model r = Vector r → Vector r +type Model r = Params r → Samples r {-| 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. +of parameters to a matrix which for each expected sample describes the partial +derivatives of the parameters. See: - * Ensure that the length of the parameter vector equals the length of the initial - parameter vector in 'levmar'. + * Ensure that the length @m@ of the parameter vector equals the length of the + initial parameter vector in 'levmar'. * Ensure that the output matrix has the dimension @n>= m@. levmar ∷ Model r -- ^ Model → Maybe (Jacobian r) -- ^ Optional jacobian - → Vector r -- ^ Initial parameters - → Vector r -- ^ Samples + → Params r -- ^ Initial parameters of length @m@ + → Samples r -- ^ Sample vector of length @n@ → Int -- ^ Maximum iterations → Options r -- ^ Minimization options → Constraints r -- ^ Constraints - → Either LevMarError (Vector r, Info r, Matrix r) + → Either LevMarError (Params r, Info r, Matrix r) instance LevMarable Float where levmar = gen_levmar slevmar_der @@ -241,12 +253,12 @@ → Model r -- ^ Model → Maybe (Jacobian r) -- ^ Optional jacobian - → Vector r -- ^ Initial parameters - → Vector r -- ^ Samples + → Params r -- ^ Initial parameters + → Samples r -- ^ Samples → Int -- ^ Maximum iterations → Options r -- ^ Options → Constraints r -- ^ Constraints - → Either LevMarError (Vector r, Info r, Matrix r) + → Either LevMarError (Params r, Info r, Matrix r) gen_levmar f_der f_dif f_bc_der @@ -255,7 +267,9 @@ f_lec_dif f_blec_der f_blec_dif - model mJac ps ys itMax opts (Constraints mLowBs mUpBs mWeights mLinC) = + model mJac ps ys itMax opts (Constraints mLowBs mUpBs mWeights mLinC) + | m == 0 = Left LevMarError -- LAPACK will crash otherwise! + | otherwise = -- All effects are contained, so we can safely perform: unsafePerformIO $ do @@ -268,10 +282,10 @@ -- 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 + psFP ← fastMallocForeignPtrArray m withForeignPtr psFP $ \psPtr → do VS.unsafeWith ps $ \psPtrInp → - copyArray psPtr psPtrInp lenPs + copyArray psPtr psPtrInp m -- Retrieve the (read-only) pointer 'ysPtr' to the samples vector 'ys' -- so we can pass it to the C function: @@ -292,7 +306,7 @@ -- Like the parameters array the matrix -- needs to be returned from this function. -- So we also allocate it using a ForeignPtr: - covarFP ← fastMallocForeignPtrArray covarLen + covarFP ← fastMallocForeignPtrArray mm withForeignPtr covarFP $ \covarPtr → -- 'cmodel' is the low-level model function which is converted @@ -306,20 +320,20 @@ let cmodel ∷ Bindings.LevMar.Model r cmodel parPtr hxPtr _ _ _ = do parFP ← newForeignPtr_ parPtr - let psV = VS.unsafeFromForeignPtr parFP 0 lenPs + let psV = VS.unsafeFromForeignPtr parFP 0 m vector = model psV VS.unsafeWith vector $ \p → copyArray hxPtr p (VS.length vector) in withModel cmodel $ \modelFunPtr → do -- All the low-level C functions share a common set of arguments. -- 'runDif' applies these arguments to the given C function 'f': - let runDif ∷ LevMarDif r → IO Int + let runDif ∷ LevMarDif r → IO CInt runDif f = f modelFunPtr psPtr ysPtr - lenPs - lenYs - itMax + (fromIntegral m) + (fromIntegral n) + (fromIntegral itMax) optsPtr infoPtr nullPtr @@ -341,14 +355,14 @@ let cjacobian ∷ Bindings.LevMar.Jacobian r cjacobian parPtr jPtr _ _ _ = do parFP ← newForeignPtr_ parPtr - let psV = VS.unsafeFromForeignPtr parFP 0 lenPs + let psV = VS.unsafeFromForeignPtr parFP 0 m matrix = jac psV vector = flatten matrix VS.unsafeWith vector $ \p → copyArray jPtr p (VS.length vector) in withJacobian cjacobian $ \jacobPtr → - let runDer ∷ LevMarDer r → IO Int + let runDer ∷ LevMarDer r → IO CInt runDer f = runDif $ f jacobPtr in if boxConstrained then if linConstrained @@ -368,20 +382,19 @@ then return $ Left $ convertLevMarError err 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 + info ← listToInfo <$> peekArray c'LM_INFO_SZ infoPtr + let psV = VS.unsafeFromForeignPtr psFP 0 m + let covarM = reshape m $ VS.unsafeFromForeignPtr covarFP 0 mm - return $ Right (psV, info, covarM) + return $ Right (psV, info, covarM) where - lenPs = VS.length ps - lenYs = VS.length ys - covarLen = lenPs*lenPs - (cMat, rhcVec) = fromJust mLinC + m = VS.length ps + n = VS.length ys + mm = m*m -- Whether the parameters are constrained by a linear equation. - linConstrained = isJust mLinC + linConstrained = isJust mLinC + (cMat, rhcVec) = fromJust mLinC -- Whether the parameters are constrained by a bounding box. boxConstrained = isJust mLowBs ∨ isJust mUpBs @@ -394,7 +407,7 @@ withLinConstraints f g = VS.unsafeWith (flatten cMat) $ \cMatPtr → VS.unsafeWith rhcVec $ \rhcVecPtr → - f ∘ g cMatPtr rhcVecPtr $ rows cMat + f ∘ g cMatPtr rhcVecPtr $ fromIntegral $ rows cMat withWeights f g = maybeWithArray mWeights $ f ∘ g @@ -451,9 +464,9 @@ -- | Ensure that these vectors have the same length as the number of parameters. data Constraints r = Constraints - { lowerBounds ∷ !(Maybe (Vector r)) -- ^ Optional lower bounds - , upperBounds ∷ !(Maybe (Vector r)) -- ^ Optional upper bounds - , weights ∷ !(Maybe (Vector r)) -- ^ Optional weights + { lowerBounds ∷ !(Maybe (Params r)) -- ^ Optional lower bounds + , upperBounds ∷ !(Maybe (Params r)) -- ^ Optional upper bounds + , weights ∷ !(Maybe (Params r)) -- ^ Optional weights , linearConstraints ∷ !(Maybe (LinearConstraints r)) -- ^ Optional linear constraints } deriving (Read, Show, Typeable) @@ -554,7 +567,7 @@ -- Handy in case you want to thow a LevMarError as an exception: instance Exception LevMarError -levmarCErrorToLevMarError ∷ [(Int, LevMarError)] +levmarCErrorToLevMarError ∷ [(CInt, LevMarError)] levmarCErrorToLevMarError = [ (c'LM_ERROR, LevMarError) , (c'LM_ERROR_LAPACK_ERROR, LapackError) @@ -569,6 +582,6 @@ --, (c'LM_ERROR_SUM_OF_SQUARES_NOT_FINITE, we don't treat this as an error) ] -convertLevMarError ∷ Int → LevMarError +convertLevMarError ∷ CInt → LevMarError convertLevMarError err = fromMaybe (error $ "Unknown levmar error: " ++ show err) (lookup err levmarCErrorToLevMarError) diff --git a/levmar.cabal b/levmar.cabal --- a/levmar.cabal +++ b/levmar.cabal @@ -1,5 +1,5 @@ name: levmar -version: 1.2.0.1 +version: 1.2.1 cabal-version: >= 1.6 build-type: Simple stability: experimental @@ -7,7 +7,7 @@ Bas van Dijk maintainer: Roel van Dijk Bas van Dijk -copyright: (c) 2009 - 2011 Roel van Dijk & Bas van Dijk +copyright: (c) 2009 - 2012 Roel van Dijk & Bas van Dijk license: BSD3 license-file: LICENSE homepage: https://github.com/basvandijk/levmar @@ -48,7 +48,7 @@ library build-depends: base >= 3 && < 4.6 , base-unicode-symbols >= 0.1.1 && < 0.3 - , bindings-levmar >= 1.0 && < 1.1 + , bindings-levmar >= 1.1 && < 1.2 , hmatrix >= 0.12 && < 0.14 , vector >= 0.8 && < 0.10 exposed-modules: Numeric.LevMar