diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,24 @@ +Copyright (c) 2013, Gard Spreemann +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/examples/Example1.hs b/examples/Example1.hs new file mode 100644 --- /dev/null +++ b/examples/Example1.hs @@ -0,0 +1,42 @@ +module Main where + +import Numeric.LBFGSB +import Numeric.LBFGSB.Result +import qualified Data.Vector.Storable as V + +-- Example optimization problem taken from driver1.f in the tarball +-- for version 3.0 of L-BFGS-B itself. + +n :: Int +n = 25 + +f :: V.Vector Double -> Double +f x = 4* (V.foldl (\s i -> s + (x V.! i - (x V.! (i-1))^2)^2) (0.25* (x V.! 0 - 1)^2) (V.enumFromN 1 (n-1))) + +g :: V.Vector Double -> V.Vector Double +g x = V.generate n (\i -> 8*(t (i-1)) - 1.6e1*(x V.! i)*(t i)) + where + t i + | i == -1 = 0.25*(x V.! 0 - 1) + | i == n-1 = 0 + | otherwise = x V.! (i+1) - (x V.! i)^2 + +bounds :: [(Maybe Double, Maybe Double)] +bounds = map (\i -> if odd i then (Just 1e0, Just 1e2) else (Just (-1e2), Just 1e2)) [1..n] + +start :: V.Vector Double +start = V.replicate n 3.0 + +main :: IO () +main = putStrLn "Testing with function and parameters from driver1.f from L-BFGS-B 3.0 distribution archive." >> + let + res = minimize 5 1e7 1e-5 Nothing bounds start f g + in + putStrLn "Full results:" >> + print res >> + putStrLn "Solution point:" >> + print (solution res) >> + putStrLn "Steps needed:" >> + print (length (backtrace res)) >> + putStrLn "Function value at solution:" >> + print (f (solution res)) diff --git a/examples/Example2.hs b/examples/Example2.hs new file mode 100644 --- /dev/null +++ b/examples/Example2.hs @@ -0,0 +1,26 @@ +module Main where + +import Numeric.LBFGSB +import Numeric.LBFGSB.Convenience +import Numeric.LBFGSB.Result +import qualified Data.Vector.Storable as V + +-- 2D Rosenbrock function with approximate derivatives. + +rosenbrock :: V.Vector Double -> Double +rosenbrock x = (1 - x0)^2 + 100*(x1-x0^2)^2 + where + x0 = x V.! 0 + x1 = x V.! 1 + +main :: IO () +main = putStrLn "Testing with Rosenbrock function, both unbounded and bounded with minimum outside bounds." >> + let + resUnbounded = minimize 5 1e0 1e-12 (Just 100) [] (V.fromList [100, 100]) rosenbrock (approximateGradient 1e-6 rosenbrock) + resBounded = minimize 5 1e0 1e-12 (Just 100) [(Nothing, Nothing), (Nothing, Just 0.1)] (V.fromList [0, 0]) rosenbrock (approximateGradient 1e-6 rosenbrock) + in + putStrLn "Unbounded:" >> + print resUnbounded >> + putStrLn "-------------------" >> + putStrLn "Bounded to (-infty, 0.1) in the second direction only:" >> + print resBounded diff --git a/examples/build.sh b/examples/build.sh new file mode 100644 --- /dev/null +++ b/examples/build.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +ghc -threaded -i../src -O2 -fforce-recomp -o example1 --make Example1.hs -llbfgsb +ghc -threaded -i../src -O2 -fforce-recomp -o example2 --make Example2.hs -llbfgsb + + diff --git a/l-bfgs-b.cabal b/l-bfgs-b.cabal new file mode 100644 --- /dev/null +++ b/l-bfgs-b.cabal @@ -0,0 +1,55 @@ +Name: l-bfgs-b +Version: 0.1 +Synopsis: Bindings to L-BFGS-B, Fortran code for limited-memory quasi-Newton bound-constrained optimization +Homepage: http://nonempty.org/software/haskell-l-bfgs-b +License: BSD3 +License-file: LICENSE +Author: Gard Spreemann +Maintainer: Gard Spreemann +Copyright: 2013 Gard Spreemann +Category: Math +Build-type: Simple +Cabal-version: >=1.4 +Description: Bindings to L-BFGS-B, Fortran code for limited-memory quasi-Newton bound-constrained optimization. + . + L-BFGS-B is a Fortran library for limited-memory quasi-Newton bound-constrained optimization written + by Ciyou Zhu, Richard Byrd, Jorge Nocedal and Jose Luis Morales. More information can be found on its + homepage , or in [1]. + . + The L-BFGS-B Fortran code is not included in this package, as I consider it a dependency. This package expects to + be able to link against version 3.0 of the L-BFGS-B code, as built by a relatively recent version of gfortran. + Instructions on how to build L-BFGS-B as a shared library + can be found at . + . + The functions provided in this package wrap FFI calls in 'unsafePerformIO', which among other things means that + the called L-BFGS-B code should not output anything. The relevant @iprint@ flag is thus set negative to suppress + output as specified in the L-BFGS-B code. However, there are two places in said code where the flag is ignored + and output still occurs. If it bothers you that code exposed as pure prints things, see + for information on a simple patch for L-BFGS-B. The SciPy project + has described the same behavior at . + . + Example on usage can be found in the included @examples@ directiory. + . + The current version has only been lightly tested, and should not be trusted for serious work. Feedback is appreciated. + . + Changes in version 0.1: + . + * There has only been cursory testing, so do not trust these bindings yet. + . + * Initial release. + . + . + \[1] R. H. Byrd, P. Lu and J. Nocedal. A Limited Memory Algorithm for Bound Constrained Optimization, (1995), SIAM Journal on Scientific and Statistical Computing , 16, 5, pp. 1190-1208. + +Library + Exposed-modules: Numeric.LBFGSB, + Numeric.LBFGSB.Convenience, + Numeric.LBFGSB.Result + Other-modules: + hs-source-dirs: src + Build-depends: base >= 4 && <5, vector + c-sources: + cc-options: + extra-libraries: lbfgsb + include-dirs: + diff --git a/src/Numeric/LBFGSB.hs b/src/Numeric/LBFGSB.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LBFGSB.hs @@ -0,0 +1,330 @@ +{-# LANGUAGE ForeignFunctionInterface #-} + +-- | Minimize functions using the Fortran L-BFGS-B library for +-- limited-memory Broyden–Fletcher–Goldfarb–Shanno bound-constrained +-- minimization. More information on assumptions and function parameters can be found at the L-BFGS-B homepage and in its source code. +-- +-- A /bound-constrained/ domain is one that is a finite product of the +-- reals, closed intervals, and half-infinite intervals. We describe +-- the factors by @('Maybe' 'Double', 'Maybe' 'Double')@, with +-- @('Just' a, 'Just' b)@ describing the closed interval [a,b], and so +-- forth. +module Numeric.LBFGSB(minimize, minimize') where + +import Data.List +import Data.Maybe +import Control.Applicative +import Foreign.C.Types +import Foreign.C.String +import Foreign.Ptr +import Foreign.Storable +import Foreign.Marshal.Array +import Foreign.Marshal.Alloc +import Foreign.Marshal.Utils +import Foreign.ForeignPtr.Safe +import System.IO.Unsafe(unsafeDupablePerformIO) +import qualified Data.Vector.Storable as V +import qualified Numeric.LBFGSB.Result as R + + +-- | Minimization using L-BFGS-B. If you only require the solution +-- point, and not the full 'R.Result', see 'minimize''. If the +-- arguments do not satisfy the given requirements, behavior is +-- undefined and the program may even crash. +minimize :: Int -- ^ @m@: The maximum number of variable metric corrections used + -- to define the limited memory matrix. /Suggestion:/ @5@. + -> Double -- ^ @factr@: Iteration stops when the relative change in function value + -- is smaller than @factr*eps@, where @eps@ is a measure of machine precision + -- generated by the Fortran code. @1e12@ is low accuracy, @1e7@ is moderate, + -- and @1e1@ is extremely high. Must be @>=1@. /Suggestion:/ @1e7@. + -> Double -- ^ @pgtol@: Iteration stops when the largest component of the projected + -- gradient is smaller than @pgtol@. Must be @>=0@. /Suggestion:/ @1e-5@. + -> Maybe Int -- ^ @'Just' steps@ means the minimization is aborted if it has not converged after + -- @steps>0@ iterations. 'Nothing' signifies no limit. + -> [(Maybe Double, Maybe Double)] -- ^ Constraints, as described in the beginning of this module. If there are + -- fewer bounds than components in @x0@, the remaining dimensions are assumed + -- to be unbounded. @[]@ thus gives unbounded minimization. Moreover, + -- if there are more bounds than components in @x0@, only as many as needed are used. + -- @'repeat' ('Just' 0, 'Just' 1)@ thus specifies the unit cube of any dimension + -- as constraint. + -> V.Vector Double -- ^ @x0@: Starting point. The point /must/ be within the bounds. + -> (V.Vector Double -> Double) -- ^ @f@: Function to minimize. /Must/ take 'V.Vector's of precisely the same + -- length as @x0@. + -> (V.Vector Double -> V.Vector Double) -- ^ @g@: Gradient of @f@. /Must/ take and return 'V.Vector's of precisely the same + -- length as @x0@. "Numeric.LBFGSB.Convenience" provides a simple + -- approximation of the gradient if you do not have the real one. + -> R.Result +minimize m factr tol steps bounds x0 f g = unsafeDupablePerformIO (runDriver m factr tol steps bounds x0 f g) + +-- | If L-BFGS-B converges within the specified number of steps, the +-- solution point is returned as @'Just' solution@. Otherwise +-- 'Nothing' is returned. The arguments are the same as for +-- 'minimize'. +minimize' :: Int + -> Double + -> Double + -> Maybe Int + -> [(Maybe Double, Maybe Double)] + -> V.Vector Double + -> (V.Vector Double -> Double) + -> (V.Vector Double -> V.Vector Double) + -> Maybe (V.Vector Double) +minimize' m factr tol steps bounds x0 f g + | R.stopReason result == R.Converged = Just $ R.solution result + | otherwise = Nothing + where + result = minimize m factr tol steps bounds x0 f g + +data DriverContext = DriverContext { pn :: Ptr CInt + , pm :: Ptr CInt + , px :: Ptr Double + , pl :: Ptr Double + , pu :: Ptr Double + , pnbd :: Ptr CInt + , pf :: Ptr Double + , pg :: Ptr Double + , pfactr :: Ptr Double + , ppgtol :: Ptr Double + , pwa :: Ptr Double + , piwa :: Ptr CInt + , ptask :: Ptr CChar + , piprint :: Ptr CInt + , pcsave :: Ptr CChar + , plsave :: Ptr CInt + , pisave :: Ptr CInt + , pdsave :: Ptr Double } + +data TaskPrefix = FG | NewX | Start | Convergence | Other String + deriving (Eq, Show) + +stringToTaskPrefix :: String -> TaskPrefix +stringToTaskPrefix s + | "FG" `isPrefixOf` s = FG + | "NEW_X" `isPrefixOf` s = NewX + | "START" `isPrefixOf` s = Start + | "CONVERGENCE" `isPrefixOf` s = Convergence + | otherwise = Other s + +taskLength :: Int +taskLength = 60 + +csaveLength :: Int +csaveLength = 60 + +startString :: String +startString = take taskLength ("START" ++ repeat ' ') + +unzipBounds :: [(Maybe Double, Maybe Double)] -> ([Double], [Double], [CInt]) +unzipBounds bounds = unzip3 (map helper bounds) + where + helper (Nothing, Nothing) = (0, 0, 0) + helper (Just l, Nothing) = (l, 0, 1) + helper (Just l, Just u) = (l, u, 2) + helper (Nothing, Just u) = (0, u, 3) + +runDriver :: Int + -> Double + -> Double + -> Maybe Int + -> [(Maybe Double, Maybe Double)] + -> V.Vector Double + -> (V.Vector Double -> Double) + -> (V.Vector Double -> V.Vector Double) + -> IO R.Result +runDriver m factr tol steps bounds x0 f g + = let + n = V.length x0 + (ls, us, bds) = unzipBounds (take n (bounds ++ repeat (Nothing, Nothing))) + in + with (fromIntegral n) $ \pn -> + with (fromIntegral m) $ \pm -> + mallocForeignPtrArray n >>= \fpx -> + withForeignPtr fpx $ \px -> + withForeignPtr ((fst . V.unsafeToForeignPtr0) x0) (\px0 -> copyArray px px0 n) >>= \_ -> + withArray ls $ \pl -> + withArray us $ \pu -> + withArray bds $ \pnbd -> + alloca $ \pf -> + allocaArray n $ \pg -> + with factr $ \pfactr -> + with tol $ \ppgtol -> + allocaArray (2*m*n + 11*m*m + 5*n + 8*m) $ \pwa -> + allocaArray (3*n) $ \piwa -> + allocaArray taskLength $ \ptask -> + withCAString startString (\pstart -> copyArray ptask pstart taskLength) >>= \_ -> + with ((-1) :: CInt) $ \piprint -> + allocaArray csaveLength $ \pcsave -> + allocaArray 4 $ \plsave -> + allocaArray 44 $ \pisave -> + allocaArray 29 $ \pdsave -> + driver (DriverContext pn pm px pl pu pnbd pf pg pfactr ppgtol pwa piwa ptask piprint pcsave plsave pisave pdsave) steps [] f g >>= \(btrace, stopReason) -> + return (R.Result (V.unsafeFromForeignPtr0 fpx n) btrace stopReason) + +driver :: DriverContext + -> Maybe Int + -> [V.Vector Double] + -> (V.Vector Double -> Double) + -> (V.Vector Double -> V.Vector Double) + -> IO ([V.Vector Double], R.StopReason) +driver context stepsLeft backtrace f g + = readTaskPrefix context >>= \task -> + case task of + Convergence -> return (backtrace, R.Converged) + Start -> makeCall context >> + driver context stepsLeft backtrace f g + FG -> readX context >>= \x -> + update context f g x >> + makeCall context >> + driver context stepsLeft backtrace f g + NewX -> readX context >>= \x -> + makeCall context >> + if maybe True (> 0) stepsLeft + then driver context (Just (subtract 1) <*> stepsLeft) (x:backtrace) f g + else return (backtrace, R.StepCount) + Other s -> return (backtrace, R.Other s) + +update :: DriverContext -> (V.Vector Double -> Double) -> (V.Vector Double -> V.Vector Double) -> V.Vector Double -> IO () +update context f g x + = poke (pf context) (f x) >> + V.unsafeWith (g x) (\pgx -> copyArray (pg context) pgx (V.length x)) + +readX :: DriverContext -> IO (V.Vector Double) +readX context + = peek (pn context) >>= (return . fromIntegral) >>= \n -> + mallocForeignPtrArray n >>= \fpxCopy -> + withForeignPtr fpxCopy (\pxCopy -> copyArray pxCopy (px context) n) >> -- fpxCopy is a foreign pointer to a *copy* of the current x + return (V.unsafeFromForeignPtr0 fpxCopy n) + +readTaskPrefix :: DriverContext -> IO TaskPrefix +readTaskPrefix context = peekCStringLen (ptask context, taskLength) >>= (return . stringToTaskPrefix) + + +makeCall :: DriverContext -> IO () +makeCall context@(DriverContext pn pm px pl pu pnbd pf pg pfactr ppgtol pwa piwa ptask piprint pcsave plsave pisave pdsave) + = fortran_setulb pn pm px pl pu pnbd pf pg pfactr ppgtol pwa piwa ptask piprint pcsave plsave pisave pdsave --(fromIntegral taskLength) (fromIntegral csaveLength) + + + +foreign import ccall "setulb_" + fortran_setulb -- # Comments below are from the L-BFGS-B Fortran source code: + :: Ptr CInt --1 n is an integer variable. + -- On entry n is the dimension of the problem. + -- On exit n is unchanged. + -> Ptr CInt --2 m is an integer variable. + -- On entry m is the maximum number of variable metric corrections + -- used to define the limited memory matrix. + -- On exit m is unchanged. + -> Ptr Double --3 x is a double precision array of dimension n. + -- On entry x is an approximation to the solution. + -- On exit x is the current approximation. + -> Ptr Double --4 l is a double precision array of dimension n. + -- On entry l is the lower bound on x. + -- On exit l is unchanged. + -> Ptr Double --5 u is a double precision array of dimension n. + -- On entry u is the upper bound on x. + -- On exit u is unchanged. + -> Ptr CInt --6 nbd is an integer array of dimension n. + -- On entry nbd represents the type of bounds imposed on the + -- variables, and must be specified as follows: + -- nbd(i)=0 if x(i) is unbounded, + -- 1 if x(i) has only a lower bound, + -- 2 if x(i) has both lower and upper bounds, and + -- 3 if x(i) has only an upper bound. + -- On exit nbd is unchanged. + -> Ptr Double --7 f is a double precision variable. + -- On first entry f is unspecified. + -- On final exit f is the value of the function at x. + -> Ptr Double --8 g is a double precision array of dimension n. + -- On first entry g is unspecified. + -- On final exit g is the value of the gradient at x. + -> Ptr Double --9 factr is a double precision variable. + -- On entry factr >= 0 is specified by the user. The iteration + -- will stop when + -- (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch + -- where epsmch is the machine precision, which is automatically + -- generated by the code. Typical values for factr: 1.d+12 for + -- low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely + -- high accuracy. + -- On exit factr is unchanged. + -> Ptr Double --10 pgtol is a double precision variable. + -- On entry pgtol >= 0 is specified by the user. The iteration + -- will stop when + -- max{|proj g_i | i = 1, ..., n} <= pgtol + -- where pg_i is the ith component of the projected gradient. + -- On exit pgtol is unchanged. + -> Ptr Double --11 wa is a double precision working array of length -- In version 3.0: + -- (2mmax + 5)nmax + 12mmax^2 + 12mmax. -- 2*m*n + 11m*m + 5*n + 8*m + -> Ptr CInt --12 iwa is an integer working array of length 3nmax. + -> Ptr CChar --13 task is a working string of characters of length 60 indicating + -- the current job when entering and quitting this subroutine. + -> Ptr CInt --14 iprint is an integer variable that must be set by the user. + -- It controls the frequency and type of output generated: + -- iprint<0 no output is generated; + -- iprint=0 print only one line at the last iteration; + -- 0100 print details of every iteration including x and g; + -- When iprint > 0, the file iterate.dat will be created to + -- summarize the iteration. + -> Ptr CChar --15 csave is a working string of characters of length 60. + -> Ptr CInt --16 lsave is a logical working array of dimension 4. + -- On exit with 'task' = NEW_X, the following information is + -- available: + -- If lsave(1) = .true. then the initial X has been replaced by + -- its projection in the feasible set; + -- If lsave(2) = .true. then the problem is constrained; + -- If lsave(3) = .true. then each variable has upper and lower + -- bounds; + -> Ptr CInt --17 isave is an integer working array of dimension 44. + -- On exit with 'task' = NEW_X, the following information is + -- available: + -- isave(22) = the total number of intervals explored in the + -- search of Cauchy points; + -- isave(26) = the total number of skipped BFGS updates before + -- the current iteration; + -- isave(30) = the number of current iteration; + -- isave(31) = the total number of BFGS updates prior the current + -- iteration; + -- isave(33) = the number of intervals explored in the search of + -- Cauchy point in the current iteration; + -- isave(34) = the total number of function and gradient + -- evaluations; + -- isave(36) = the number of function value or gradient + -- evaluations in the current iteration; + -- if isave(37) = 0 then the subspace argmin is within the box; + -- if isave(37) = 1 then the subspace argmin is beyond the box; + -- isave(38) = the number of free variables in the current + -- iteration; + -- isave(39) = the number of active constraints in the current + -- iteration; + -- n + 1 - isave(40) = the number of variables leaving the set of + -- active constraints in the current iteration; + -- isave(41) = the number of variables entering the set of active + -- constraints in the current iteration. + -> Ptr Double --18 dsave is a double precision working array of dimension 29. + -- On exit with 'task' = NEW_X, the following information is + -- available: + -- dsave(1) = current 'theta' in the BFGS matrix; + -- dsave(2) = f(x) in the previous iteration; + -- dsave(3) = factr*epsmch; + -- dsave(4) = 2-norm of the line search direction vector; + -- dsave(5) = the machine precision epsmch generated by the code; + -- dsave(7) = the accumulated time spent on searching for + -- Cauchy points; + -- dsave(8) = the accumulated time spent on + -- subspace minimization; + -- dsave(9) = the accumulated time spent on line search; + -- dsave(11) = the slope of the line search function at + -- the current point of line search; + -- dsave(12) = the maximum relative step length imposed in + -- line search; + -- dsave(13) = the infinity norm of the projected gradient; + -- dsave(14) = the relative step length in the line search; + -- dsave(15) = the slope of the line search function at + -- the starting point of the line search; + -- dsave(16) = the square of the 2-norm of the line search + -- direction vector. +-- -> CInt -> CInt -- FORTRAN string lengths piled at the end. Fixme, is this right? + -> IO () + diff --git a/src/Numeric/LBFGSB/Convenience.hs b/src/Numeric/LBFGSB/Convenience.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LBFGSB/Convenience.hs @@ -0,0 +1,29 @@ +{-# LANGUAGE FlexibleContexts #-} + +-- | Some functions that can be useful together with L-BFGS-B. +module Numeric.LBFGSB.Convenience(approximateGradient, listFunction, vectorFunction) where + +import qualified Data.Vector.Storable as V +import qualified Data.Vector.Generic as GV + +-- | @'approximateGradient' h f x@ is an approximation of the gradient +-- of @f@ at @x@, computed using central differences with step size +-- @h@. +approximateGradient :: (GV.Vector v Double) => Double -> (v Double -> Double) -> (v Double -> v Double) +approximateGradient h f x = GV.generate (GV.length x) (\i -> central h (\y -> f (replaceAt i y x)) (x GV.! i)) + +-- | Turn a function on lists into one on 'V.Storable' 'V.Vector's. +listFunction :: ([Double] -> Double) -> (V.Vector Double -> Double) +listFunction f = f . V.toList + +-- | Turn a function on any generic 'GV.Vector's into one on 'V.Storable' 'V.Vector's. +vectorFunction :: (GV.Vector v Double) => (v Double -> Double) -> (V.Vector Double -> Double) +vectorFunction f = f . GV.convert + +replaceAt :: (GV.Vector v a) => Int -> a -> v a -> v a +replaceAt n x xs + | n > GV.length xs - 1 || n < 0 = xs + | otherwise = xs GV.// [(n, x)] + +central :: Double -> (Double -> Double) -> Double -> Double +central h f x = (f (x+h) - f (x-h))/(2*h) diff --git a/src/Numeric/LBFGSB/Result.hs b/src/Numeric/LBFGSB/Result.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LBFGSB/Result.hs @@ -0,0 +1,23 @@ +-- | The 'Result' data type encodes the minimization solution, as well +-- as auxiliary information about the minimization process. +module Numeric.LBFGSB.Result where + +import qualified Data.Vector.Storable as V + +-- | Stores the result of the minimization process. +data Result = Result { + solution :: V.Vector Double -- ^ Solution point /if the minimization completed successfully/. See 'stopReason'. + , backtrace :: [V.Vector Double] -- ^ The steps taken to reach the solution, in reverse order. Does not include the starting point. + , stopReason :: StopReason -- ^ The reason L-BFGS-B terminated. Only if this is + -- 'Converted' should you consider the solution correct! + } + deriving (Show) + +-- | The reason L-BFGS-B terminated. +data StopReason = + Converged -- ^ The solution converged. + | StepCount -- ^ The number of steps exceeded the user's request. + | Other String -- ^ Something else occured. In @'Other' s@, @s@ is the contents + -- of L-BFGS-B's @task@ variable on exit, as documented in the + -- source code of L-BFGS-B itself. + deriving (Eq, Show)