nonlinear-optimization 0.3.1 → 0.3.2
raw patch · 2 files changed
+109/−61 lines, 2 files
Files
Numeric/Optimization/Algorithms/HagerZhang05.hsc view
@@ -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.
nonlinear-optimization.cabal view
@@ -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