diff --git a/Numeric/Optimization/Algorithms/HagerZhang05.hsc b/Numeric/Optimization/Algorithms/HagerZhang05.hsc --- a/Numeric/Optimization/Algorithms/HagerZhang05.hsc +++ b/Numeric/Optimization/Algorithms/HagerZhang05.hsc @@ -53,10 +53,19 @@ import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Storable as S import qualified Data.Vector.Storable.Mutable as SM +import Control.Applicative import Control.Exception (bracket) import Control.Monad.Primitive (PrimMonad(..)) import Foreign import Foreign.C + +#ifdef DEBUG +import Debug.Trace (trace) +#else +trace :: String -> a -> a +trace _ x = x +#endif + #include "cg_user.h" -- $mainFunction @@ -95,14 +104,13 @@ -- Allocate everything. (ret, stats) <- SM.unsafeWith x $ \x_ptr -> - alloca $ \stats_ptr -> - alloca $ \param_ptr -> + allocaSet (Statistics 0 0 0 0 0) $ \stats_ptr -> + allocaSet params $ \param_ptr -> bracket (mkCFunction cf) freeHaskellFunPtr $ \cf_ptr -> bracket (mkCGradient cg) freeHaskellFunPtr $ \cg_ptr -> bracket (mkCCombined cc) freeHaskellFunPtr $ \cc_ptr -> allocateWorkSpace n $ \work_ptr -> do -- Go to C land. - poke param_ptr params ret <- cg_descent x_ptr (fromIntegral n) stats_ptr param_ptr grad_tol cf_ptr cg_ptr cc_ptr work_ptr @@ -113,6 +121,12 @@ x' <- G.unsafeFreeze x return $ ret `seq` (x', ret, stats) +-- | Allocates as 'alloca' and sets the memory area. +allocaSet :: Storable a => a -> (Ptr a -> IO b) -> IO b +allocaSet x f = alloca $ \x_ptr -> do + poke x_ptr x + f x_ptr + -- | Allocates enough work space for CG_DESCENT. If the number -- of dimensions is "small enough" then we allocate on the stack, -- otherwise we allocate via malloc. @@ -163,25 +177,46 @@ VFunction :: G.Vector v Double => (v Double -> Double) -> Function Simple - MFunction :: (forall m. PrimMonad m + MFunction :: (forall m. (PrimMonad m, Functor m) => PointMVector m -> m Double) -> Function Mutable + + +-- | Copies the input array from a mutable storable vector to any +-- pure vector. Used to convert pure functions into mutable +-- ones. +copyInput :: (PrimMonad m, G.Vector v Double) + => SM.MVector (PrimState m) Double + -> m (v Double) +copyInput mx = do + let s = trace " copyInput start" $ GM.length mx + mz <- GM.new s + let go i | i >= s = return () + | otherwise = GM.unsafeRead mx i >>= + GM.unsafeWrite mz i >> go (i+1) + go 0 + trace " stop" $ G.unsafeFreeze mz + +-- | Copies the output array from any pure vector to a mutable +-- storable array. Used to convert pure functions that return +-- the gradient into mutable ones. +copyOutput :: (PrimMonad m, G.Vector v Double) + => SM.MVector (PrimState m) Double + -> v Double + -> m () +copyOutput mret r = go $ trace " copyOutput start" $ 0 + where + s = min (GM.length mret) (G.length r) + go i | i >= s = trace " stop" $ return () + | otherwise = let !x = G.unsafeIndex r i + in GM.unsafeWrite mret i x >> go (i+1) + + + mutableF :: Function t -> Function Mutable -mutableF (VFunction f) = MFunction f' - where - f' mx = do - -- Copy the input to an immutable vector. - let s = GM.length mx - mz <- GM.new s - let go i | i > s = return () - | otherwise = GM.unsafeRead mx i >>= - GM.unsafeWrite mz i >> go (i+1) - go 0 - z <- G.unsafeFreeze mz - -- Run the user function. - return (f z) +mutableF (VFunction f) = MFunction (\mx -> f <$> copyInput mx) mutableF (MFunction f) = MFunction f prepareF :: Function Mutable -> CFunction @@ -189,13 +224,32 @@ \x_ptr n -> do let n' = fromIntegral n x_fptr <- newForeignPtr_ x_ptr - f (SM.unsafeFromForeignPtr x_fptr 0 n') + let x = SM.unsafeFromForeignPtr x_fptr 0 n' +#ifdef DEBUG + putStr $ unlines [ + "--> function:", + " x: " ++ showV x] +#endif + r <- f x +#ifdef DEBUG + putStrLn $ " r: " ++ show r +#endif + return r prepareF _ = error "HagerZhang05.prepareF: never here" +#ifdef DEBUG +showV :: SM.IOVector Double -> String +showV m = show $ go 0 (GM.length m) + where + go i n | i == n = [] + | otherwise = let !v = unsafePerformIO (GM.read m i) + in v : go (i+1) n +#endif + -- | Function calculating the value of the gradient of the -- objective function @f@ at a point @x@. -- @@ -206,7 +260,7 @@ VGradient :: G.Vector v Double => (v Double -> v Double) -> Gradient Simple - MGradient :: (forall m. PrimMonad m + MGradient :: (forall m. (PrimMonad m, Functor m) => PointMVector m -> GradientMVector m -> m ()) @@ -214,33 +268,27 @@ mutableG :: Gradient t -> Gradient Mutable mutableG (VGradient f) = MGradient f' where - f' mx mret = do - -- Copy the input to an immutable vector. - let s = GM.length mx - mz <- GM.new s - let go i | i > s = return () - | otherwise = GM.unsafeRead mx i >>= - GM.unsafeWrite mz i >> go (i+1) - go 0 - z <- G.unsafeFreeze mz - -- Run the user function. - let !r = f z - -- Copy the output to an immutable vector - let s' = min s (G.length r) - go' i | i > s' = return () - | otherwise = let !x = G.unsafeIndex r i - in GM.unsafeWrite mret i x >> go (i+1) - go' 0 + f' mx mret = f <$> copyInput mx >>= copyOutput mret mutableG (MGradient f) = MGradient f + prepareG :: Gradient Mutable -> CGradient prepareG (MGradient f) = \ret_ptr x_ptr n -> do let n' = fromIntegral n x_fptr <- newForeignPtr_ x_ptr ret_fptr <- newForeignPtr_ ret_ptr - f (SM.unsafeFromForeignPtr x_fptr 0 n') - (SM.unsafeFromForeignPtr ret_fptr 0 n') + let x = SM.unsafeFromForeignPtr x_fptr 0 n' + r = SM.unsafeFromForeignPtr ret_fptr 0 n' +#ifdef DEBUG + putStr $ unlines [ + "--> gradient:", + " x: " ++ showV x] +#endif + f x r +#ifdef DEBUG + putStrLn $ " r: " ++ showV r +#endif prepareG _ = error "HagerZhang05.prepareG: never here" @@ -257,7 +305,7 @@ VCombined :: G.Vector v Double => (v Double -> (Double, v Double)) -> Combined Simple - MCombined :: (forall m. PrimMonad m + MCombined :: (forall m. (PrimMonad m, Functor m) => PointMVector m -> GradientMVector m -> m Double) @@ -266,23 +314,8 @@ mutableC (VCombined f) = MCombined f' where f' mx mret = do - -- Copy the input to an immutable vector. - let s = GM.length mx - mz <- GM.new s - let go i | i > s = return () - | otherwise = GM.unsafeRead mx i >>= - GM.unsafeWrite mz i >> go (i+1) - go 0 - z <- G.unsafeFreeze mz - -- Run the user function. - let !(v,r) = f z - -- Copy the output to an immutable vector - let s' = min s (G.length r) - go' i | i > s' = return () - | otherwise = let !x = G.unsafeIndex r i - in GM.unsafeWrite mret i x >> go (i+1) - go' 0 - -- Return the value + (v,r) <- f <$> copyInput mx + copyOutput mret r return v mutableC (MCombined f) = MCombined f @@ -292,8 +325,18 @@ let n' = fromIntegral n x_fptr <- newForeignPtr_ x_ptr ret_fptr <- newForeignPtr_ ret_ptr - f (SM.unsafeFromForeignPtr x_fptr 0 n') - (SM.unsafeFromForeignPtr ret_fptr 0 n') + let x = SM.unsafeFromForeignPtr x_fptr 0 n' + r = SM.unsafeFromForeignPtr ret_fptr 0 n' +#ifdef DEBUG + putStr $ unlines [ + "--> combined:", + " x: " ++ showV x] +#endif + v <- f x r +#ifdef DEBUG + putStrLn $ " r: " ++ show v ++ ", " ++ showV r +#endif + return v prepareC _ = error "HagerZhang05.prepareC: never here" -- | Combine two separated functions into a single, combined one. diff --git a/nonlinear-optimization.cabal b/nonlinear-optimization.cabal --- a/nonlinear-optimization.cabal +++ b/nonlinear-optimization.cabal @@ -3,7 +3,7 @@ Tested-With: GHC Category: Math Name: nonlinear-optimization -Version: 0.3.1 +Version: 0.3.2 Stability: experimental License: GPL License-File: LICENSE @@ -27,13 +27,17 @@ matrices, although it would be possible to use something like @hmatrix@ easily. . - Currently only CM_DESCENT method is implemented. + Currently only CG_DESCENT method is implemented. Extra-Source-Files: CG_DESCENT-C-3.0/cg_descent.c, CG_DESCENT-C-3.0/cg_descent.h, CG_DESCENT-C-3.0/cg_user.h, CG_DESCENT-C-3.0/README +Flag Debug + Description: Enable some debug statements. + Default: False + Library Build-Depends: base >= 3 && < 5, vector >= 0.5 && < 0.6, @@ -56,4 +60,5 @@ Build-Tools: hsc2hs Extra-Libraries: m GHC-Options: -Wall - + if flag(Debug) + CPP-Options: -DDEBUG