packages feed

l-bfgs-b (empty) → 0.1

raw patch · 9 files changed

+537/−0 lines, 9 filesdep +basedep +vectorsetup-changed

Dependencies added: base, vector

Files

+ LICENSE view
@@ -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 <organization> 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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ examples/Example1.hs view
@@ -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))
+ examples/Example2.hs view
@@ -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
+ examples/build.sh view
@@ -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++
+ l-bfgs-b.cabal view
@@ -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 <gspreemann@gmail.com>+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 <http://users.eecs.northwestern.edu/~nocedal/lbfgsb.html>, 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 <http://nonempty.org/software/haskell-l-bfgs-b>.+		     .+                     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 +                     <http://nonempty.org/software/haskell-l-bfgs-b> for information on a simple patch for L-BFGS-B. The SciPy project+                     has described the same behavior at <http://projects.scipy.org/scipy/ticket/1742>.+                     .+		     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:            +
+ src/Numeric/LBFGSB.hs view
@@ -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 <http://users.eecs.northwestern.edu/~nocedal/lbfgsb.html> 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;+                           --        0<iprint<99 print also f and |proj g| every iprint iterations;+                           --        iprint=99   print details of every iteration except n-vectors;+                           --        iprint=100  print also the changes of active set and final x;+                           --        iprint>100  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 ()+
+ src/Numeric/LBFGSB/Convenience.hs view
@@ -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)
+ src/Numeric/LBFGSB/Result.hs view
@@ -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)