packages feed

numeric-optimization (empty) → 0.1.0.0

raw patch · 9 files changed

+1155/−0 lines, 9 filesdep +HUnitdep +basedep +constraintssetup-changed

Dependencies added: HUnit, base, constraints, containers, data-default-class, hmatrix, hspec, lbfgs, nonlinear-optimization, numeric-optimization, primitive, vector

Files

+ CHANGELOG.md view
@@ -0,0 +1,11 @@+# Changelog for `numeric-optimization`++All notable changes to this project will be documented in this file.++The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),+and this project adheres to the+[Haskell Package Versioning Policy](https://pvp.haskell.org/).++## Unreleased++## 0.1.0.0 - YYYY-MM-DD
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Masahiro Sakai (c) 2023++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 Masahiro Sakai nor the names of other+      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 COPYRIGHT+OWNER OR CONTRIBUTORS 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.
+ README.md view
@@ -0,0 +1,57 @@+# numeric-optimization++Unified interface to various numerical optimization algorithms.++Note that the package name is numeric-optimization and not numeri**cal**-optimization.+The name `numeric-optimization` comes from the module name `Numeric.Optimization`.+++## Example Usage++```haskell+{-# LANGUAGE OverloadedLists #-}+import Data.Vector.Storable (Vector)+import Numeric.Optimization++main :: IO ()+main = do+  result <- minimize LBFGS def (WithGrad rosenbrock rosenbrock') [-3,-4]+  print (resultSuccess result)  -- True+  print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]+  print (resultValue result)  -- 1.8129771632403013e-18++-- https://en.wikipedia.org/wiki/Rosenbrock_function+rosenbrock :: Vector Double -> Double+rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)++rosenbrock' :: Vector Double -> Vector Double+rosenbrock' [x,y] =+  [ 2 * (1 - x) * (-1) + 100 * 2 * (y - sq x) * (-2) * x+  , 100 * 2 * (y - sq x)+  ]++sq :: Floating a => a -> a+sq x = x ** 2+```++## Supported Algorithms++|Algorithm|Solver implemention|Haskell binding| |+|---------|-------------------|---------------|-|+|CG\_DESCENT|[CG_DESCENT-C](https://www.math.lsu.edu/~hozhang/SoftArchive/CG_DESCENT-C-3.0.tar.gz)|[nonlinear-optimization](https://hackage.haskell.org/package/nonlinear-optimization)|Requires `with-cg-descent` flag|+|Limited memory BFGS (L-BFGS)|[liblbfgs](https://github.com/chokkan/liblbfgs)|[lbfgs](https://hackage.haskell.org/package/lbfgs)|+|Newton's method|Pure Haskell implementation using [HMatrix](https://hackage.haskell.org/package/hmatrix)|-|+++## Related Packages++* Packages for using with automatic differentiation:+  * [numerical-optimization-ad](https://hackage.haskell.org/package/numerical-optimization-ad) for using with [ad](https://hackage.haskell.org/package/ad) package+  * [numerical-optimization-backprop](https://hackage.haskell.org/package/numerical-optimization-backprop) for using with [backprop](https://hackage.haskell.org/package/backprop) package+* [MIP](https://hackage.haskell.org/package/MIP) for solving linear programming and mixed-integer linear programming problems++## LICENSE++The code in thie packaged is licensed under [BSD-3-Clause](LIENSE).++If you enable `with-cg-descent` flag, it uses GPL-licensed packages and the resulting binary should be distributed under GPL.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ examples/rosenbrock.hs view
@@ -0,0 +1,23 @@+{-# LANGUAGE OverloadedLists #-}+import Data.Vector.Storable (Vector)+import Numeric.Optimization++main :: IO ()+main = do+  result <- minimize LBFGS def (WithGrad rosenbrock rosenbrock') [-3,-4]+  print (resultSuccess result)  -- True+  print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]+  print (resultValue result)  -- 1.8129771632403013e-18++-- https://en.wikipedia.org/wiki/Rosenbrock_function+rosenbrock :: Vector Double -> Double+rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)++rosenbrock' :: Vector Double -> Vector Double+rosenbrock' [x,y] =+  [ 2 * (1 - x) * (-1) + 100 * 2 * (y - sq x) * (-2) * x+  , 100 * 2 * (y - sq x)+  ]++sq :: Floating a => a -> a+sq x = x ** 2
+ numeric-optimization.cabal view
@@ -0,0 +1,100 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.35.1.+--+-- see: https://github.com/sol/hpack++name:           numeric-optimization+version:        0.1.0.0+synopsis:       Unified interface to various numerical optimization algorithms+description:    Please see the README on GitHub at <https://github.com/msakai/nonlinear-optimization-ad/tree/master/numeric-optimization#readme>+category:       Math, Algorithms, Optimisation, Optimization+homepage:       https://github.com/msakai/numeric-optimization#readme+bug-reports:    https://github.com/msakai/numeric-optimization/issues+author:         Masahiro Sakai+maintainer:     masahiro.sakai@gmail.com+copyright:      Masahiro Sakai &lt;masahiro.sakai@gmail.com&gt;+license:        BSD3+license-file:   LICENSE+build-type:     Simple+tested-with:+    GHC == 9.4.5+  , GHC == 9.2.7+  , GHC == 9.0.2+  , GHC == 8.10.7+  , GHC == 8.8.4+  , GHC == 8.6.5+extra-source-files:+    README.md+    CHANGELOG.md++source-repository head+  type: git+  location: https://github.com/msakai/numeric-optimization++flag build-examples+  description: Build example programs+  manual: True+  default: False++flag with-cg-descent+  description: Enable CGDescent optimization algorithm provided by nonlinear-optimization package and CG_DESCENT-C library. Since they are licensed under GPL, setting this flag True implies that resulting binary is also under GPL.+  manual: True+  default: False++library+  exposed-modules:+      Numeric.Optimization+  other-modules:+      Paths_numeric_optimization+  hs-source-dirs:+      src+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints+  build-depends:+      base >=4.12 && <5+    , constraints+    , data-default-class >=0.1.2.0 && <0.2+    , hmatrix >=0.20.0.0+    , lbfgs ==0.1.*+    , primitive >=0.6.4.0+    , vector >=0.12.0.2 && <0.14+  default-language: Haskell2010+  if flag(with-cg-descent)+    cpp-options: -DWITH_CG_DESCENT+    build-depends:+        nonlinear-optimization >=0.3.7 && <0.4+  else+    cpp-options:  ++executable rosenbrock+  main-is: rosenbrock.hs+  other-modules:+      Paths_numeric_optimization+  hs-source-dirs:+      examples+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      base >=4.12 && <5+    , data-default-class >=0.1.2.0 && <0.2+    , numeric-optimization+    , vector >=0.12.0.2 && <0.14+  default-language: Haskell2010++test-suite numeric-optimization-test+  type: exitcode-stdio-1.0+  main-is: Spec.hs+  other-modules:+      IsClose+      Paths_numeric_optimization+  hs-source-dirs:+      test+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      HUnit >=1.6.0.0 && <1.7+    , base >=4.12 && <5+    , containers >=0.6.0.1 && <0.7+    , data-default-class >=0.1.2.0 && <0.2+    , hspec >=2.7.1 && <3.0+    , numeric-optimization+    , vector >=0.12.0.2 && <0.14+  default-language: Haskell2010
+ src/Numeric/Optimization.hs view
@@ -0,0 +1,749 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.Optimization+-- Copyright   :  (c) Masahiro Sakai 2023+-- License     :  BSD-style+--+-- Maintainer  :  masahiro.sakai@gmail.com+-- Stability   :  provisional+-- Portability :  non-portable+--+-- This module aims to provides unifined interface to various numerical+-- optimization, like [scipy.optimize](https://docs.scipy.org/doc/scipy/reference/optimize.html) in Python.+--+-- In this module, you need to explicitly provide the function to calculate the+-- gradient, -- but you can use @numeric-optimization-ad@ or+-- @numeric-optimization-backprop@ to define it using automatic differentiation.+--+-----------------------------------------------------------------------------+module Numeric.Optimization+  (++  -- * Main function+    minimize++  -- * Problem specification+  --+  -- $problemDefinition+  , IsProblem (..)+  , HasGrad (..)+  , HasHessian (..)+  , Constraint (..)+  , boundsUnconstrained+  , isUnconstainedBounds+  -- ** Wrapper types+  , WithGrad (..)+  , WithHessian (..)+  , WithBounds (..)+  , WithConstraints (..)++  -- * Algorithm selection+  , Method (..)+  , isSupportedMethod+  , Params (..)++  -- * Result+  , Result (..)+  , Statistics (..)+  , OptimizationException (..)++  -- * Utilities and Re-export+  , Default (..)+  , Optionally (..)+  , hasOptionalDict+  ) where++import Control.Exception+import Control.Monad.Primitive+import Control.Monad.ST+import Data.Coerce+import Data.Constraint (Dict (..))+import Data.Default.Class+import Data.Functor.Contravariant+import Data.IORef+import Data.Maybe+import qualified Data.Vector as V+import Data.Vector.Storable (Vector)+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import qualified Data.Vector.Storable.Mutable as VSM+import Foreign.C+import qualified Numeric.LBFGS.Vector as LBFGS+#ifdef WITH_CG_DESCENT+import qualified Numeric.Optimization.Algorithms.HagerZhang05 as CG+#endif+import Numeric.LinearAlgebra (Matrix)+import qualified Numeric.LinearAlgebra as LA+++-- | Selection of numerical optimization algorithms+data Method+  = CGDescent+    -- ^ Conjugate gradient method based on Hager and Zhang [1].+    --+    -- The implementation is provided by nonlinear-optimization package [3]+    -- which is a binding library of [2].+    --+    -- This method requires gradient but does not require hessian.+    --+    -- * [1] Hager, W. W. and Zhang, H.  /A new conjugate gradient/+    --   /method with guaranteed descent and an efficient line/+    --   /search./ Society of Industrial and Applied Mathematics+    --   Journal on Optimization, 16 (2005), 170-192.+    --+    -- * [2] <https://www.math.lsu.edu/~hozhang/SoftArchive/CG_DESCENT-C-3.0.tar.gz>+    --+    -- * [3] <https://hackage.haskell.org/package/nonlinear-optimization>+  | LBFGS+    -- ^ Limited memory BFGS (L-BFGS) algorithm [1]+    --+    -- The implementtion is provided by lbfgs package [2]+    -- which is a binding of liblbfgs [3].+    --+    -- This method requires gradient but does not require hessian.+    --+    -- * [1] <https://en.wikipedia.org/wiki/Limited-memory_BFGS>+    --+    -- * [2] <https://hackage.haskell.org/package/lbfgs>+    --+    -- * [3] <https://github.com/chokkan/liblbfgs>+  | Newton+    -- ^ Native implementation of Newton method+    --+    -- This method requires both gradient and hessian.+  deriving (Eq, Ord, Enum, Show, Bounded)+++-- | Whether a 'Method' is supported under the current environment.+isSupportedMethod :: Method -> Bool+isSupportedMethod LBFGS = True+#ifdef WITH_CG_DESCENT+isSupportedMethod CGDescent = True+#else+isSupportedMethod CGDescent = False+#endif+isSupportedMethod Newton = True+++-- | Parameters for optimization algorithms+--+-- TODO:+--+-- * How to pass algorithm specific parameters?+--+-- * Separate 'callback' from other more concrete serializeable parameters?+data Params a+  = Params+  { paramsCallback :: Maybe (a -> IO Bool)+    -- ^ If callback function returns @True@, the algorithm execution is terminated.+  , paramsTol :: Maybe Double+    -- ^ Tolerance for termination. When 'tol' is specified, the selected algorithm sets+    -- some relevant solver-specific tolerance(s) equal to 'tol'.+  }++instance Default (Params a) where+  def =+    Params+    { paramsCallback = Nothing+    , paramsTol = Nothing+    }++instance Contravariant Params where+  contramap f params =+    params+    { paramsCallback = fmap ((. f)) (paramsCallback params)+    }+++-- | Optimization result+data Result a+  = Result+  { resultSuccess :: Bool+    -- ^ Whether or not the optimizer exited successfully.+  , resultMessage :: String+    -- ^ Description of the cause of the termination.+  , resultSolution :: a+    -- ^ Solution+  , resultValue :: Double+    -- ^ Value of the function at the solution.+  , resultGrad :: Maybe a+    -- ^ Gradient at the solution+  , resultHessian :: Maybe (Matrix Double)+    -- ^ Hessian at the solution; may be an approximation.+  , resultHessianInv :: Maybe (Matrix Double)+    -- ^ Inverse of Hessian at the solution; may be an approximation.+  , resultStatistics :: Statistics+    -- ^ Statistics of optimizaion process+  }++instance Functor Result where+  fmap f result =+    result+    { resultSolution = f (resultSolution result)+    , resultGrad = fmap f (resultGrad result)+    }+++-- | Statistics of optimizaion process+data Statistics+  = Statistics+  { totalIters :: Int+    -- ^ Total number of iterations.+  , funcEvals :: Int+    -- ^ Total number of function evaluations.+  , gradEvals :: Int+    -- ^ Total number of gradient evaluations.+  , hessEvals :: Int+    -- ^ Total number of hessian evaluations.+  }+++-- | The bad things that can happen when you use the library.+data OptimizationException+  = UnsupportedProblem String+  | UnsupportedMethod Method+  | GradUnavailable+  | HessianUnavailable+  deriving (Show)++instance Exception OptimizationException++++-- $problemDefinition+--+-- Problems are specified by types of 'IsProblem' type class.+--+-- In the simplest case, @'VS.Vector' Double -> Double@ is a instance+-- of 'IsProblem' class. It is enough if your problem does not have+-- constraints and the selected algorithm does not further information+-- (e.g. gradients and hessians),+--+-- You can equip a problem with other information using wrapper types:+--+-- * 'WithBounds'+--+-- * 'WithConstraints'+--+-- * 'WithGrad'+--+-- * 'WithHessian'+--+-- If you need further flexibility or efficient implementation, you can+-- define instance of 'IsProblem' by yourself.++-- | Optimization problems+class IsProblem prob where+  -- | Objective function+  --+  -- It is called @fun@ in @scipy.optimize.minimize@.+  func :: prob -> Vector Double -> Double++  -- | Bounds+  --+  bounds :: prob -> Maybe (V.Vector (Double, Double))+  bounds _ = Nothing++  -- | Constraints+  constraints :: prob -> [Constraint]+  constraints _ = []++  {-# MINIMAL func #-}+++-- | Optimization problem equipped with gradient information+class IsProblem prob => HasGrad prob where+  -- | Gradient of a function computed by 'func'+  --+  -- It is called @jac@ in @scipy.optimize.minimize@.+  grad :: prob -> Vector Double -> Vector Double+  grad prob = snd . grad' prob++  -- | Pair of 'func' and 'grad'+  grad' :: prob -> Vector Double -> (Double, Vector Double)+  grad' prob x = runST $ do+    gret <- VGM.new (VG.length x)+    y <- grad'M prob x gret+    g <- VG.unsafeFreeze gret+    return (y, g)++  -- | Similar to 'grad'' but destination passing style is used for gradient vector+  grad'M :: PrimMonad m => prob -> Vector Double -> VSM.MVector (PrimState m) Double -> m Double+  grad'M prob x gvec = do+    let y = func prob x+    VG.imapM_ (VGM.write gvec) (grad prob x)+    return y++  {-# MINIMAL grad | grad' | grad'M #-}+++-- | Optimization problem equipped with hessian information+class IsProblem prob => HasHessian prob where+  -- | Hessian of a function computed by 'func'+  --+  -- It is called @hess@ in @scipy.optimize.minimize@.+  hessian :: prob -> Vector Double -> Matrix Double++  -- | The product of the hessian @H@ of a function @f@ at @x@ with a vector @x@.+  --+  -- It is called @hessp@ in @scipy.optimize.minimize@.+  --+  -- See also <https://hackage.haskell.org/package/ad-4.5.4/docs/Numeric-AD.html#v:hessianProduct>.+  hessianProduct :: prob -> Vector Double -> Vector Double -> Vector Double+  hessianProduct prob x v = hessian prob x LA.#> v++  {-# MINIMAL hessian #-}+++-- | Optional constraint+class Optionally c where+  optionalDict :: Maybe (Dict c)+++-- | Utility function to define 'Optionally' instances+hasOptionalDict :: c => Maybe (Dict c)+hasOptionalDict = Just Dict+++-- | Type of constraint+--+-- Currently, no constraints are supported.+data Constraint++-- | Bounds for unconstrained problems, i.e. (-∞,+∞).+boundsUnconstrained :: Int -> V.Vector (Double, Double)+boundsUnconstrained n = V.replicate n (-1/0, 1/0)++-- | Whether all lower bounds are -∞ and all upper bounds are +∞.+isUnconstainedBounds :: V.Vector (Double, Double) -> Bool+isUnconstainedBounds = V.all p+  where+    p (lb, ub) = isInfinite lb && lb < 0 && isInfinite ub && ub > 0+++-- | Minimization of scalar function of one or more variables.+--+-- This function is intended to provide functionality similar to Python's @scipy.optimize.minimize@.+--+-- Example:+--+-- > {-# LANGUAGE OverloadedLists #-}+-- >+-- > import Data.Vector.Storable (Vector)+-- > import Numeric.Optimization+-- >+-- > main :: IO ()+-- > main = do+-- >   (x, result, stat) <- minimize LBFGS def (WithGrad rosenbrock rosenbrock') [-3,-4]+-- >   print (resultSuccess result)  -- True+-- >   print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]+-- >   print (resultValue result)  -- 1.8129771632403013e-18+-- >+-- > -- https://en.wikipedia.org/wiki/Rosenbrock_function+-- > rosenbrock :: Vector Double -> Double+-- > rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)+-- >+-- > rosenbrock' :: Vector Double -> Vector Double+-- > rosenbrock' [x,y] =+-- >   [ 2 * (1 - x) * (-1) + 100 * 2 * (y - sq x) * (-2) * x+-- >   , 100 * 2 * (y - sq x)+-- >   ]+-- >+-- > sq :: Floating a => a -> a+-- > sq x = x ** 2+minimize+  :: forall prob. (IsProblem prob, Optionally (HasGrad prob), Optionally (HasHessian prob))+  => Method  -- ^ Numerical optimization algorithm to use+  -> Params (Vector Double) -- ^ Parameters for optimization algorithms. Use 'def' as a default.+  -> prob  -- ^ Optimization problem to solve+  -> Vector Double  -- ^ Initial value+  -> IO (Result (Vector Double))+#ifdef WITH_CG_DESCENT+minimize CGDescent =+  case optionalDict @(HasGrad prob) of+    Just Dict -> minimize_CGDescent+    Nothing -> \_ _ _ -> throwIO GradUnavailable+#endif+minimize LBFGS =+  case optionalDict @(HasGrad prob) of+    Just Dict -> minimize_LBFGS+    Nothing -> \_ _ _ -> throwIO GradUnavailable+minimize Newton =+  case optionalDict @(HasGrad prob) of+    Nothing -> \_ _ _ -> throwIO GradUnavailable+    Just Dict ->+      case optionalDict @(HasHessian prob) of+        Nothing -> \_ _ _ -> throwIO HessianUnavailable+        Just Dict -> minimize_Newton+minimize method = \_ _ _ -> throwIO (UnsupportedMethod method)+++#ifdef WITH_CG_DESCENT++minimize_CGDescent :: HasGrad prob => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double))+minimize_CGDescent _params prob _ | not (isNothing (bounds prob)) = throwIO (UnsupportedProblem "CGDescent does not support bounds")+minimize_CGDescent _params prob _ | not (null (constraints prob)) = throwIO (UnsupportedProblem "CGDescent does not support constraints")+minimize_CGDescent params prob x0 = do+  let grad_tol = fromMaybe 1e-6 $ paramsTol params++      cg_params =+        CG.defaultParameters+        { CG.printFinal = False+        }++      mf :: forall m. PrimMonad m => CG.PointMVector m -> m Double+      mf mx = do+        x <- VG.unsafeFreeze mx+        return $ func prob x++      mg :: forall m. PrimMonad m => CG.PointMVector m -> CG.GradientMVector m -> m ()+      mg mx mret = do+        x <- VG.unsafeFreeze mx+        _ <- grad'M prob x mret+        return ()++      mc :: forall m. PrimMonad m => CG.PointMVector m -> CG.GradientMVector m -> m Double+      mc mx mret = do+        x <- VG.unsafeFreeze mx+        grad'M prob x mret++  (x, result, stat) <-+    CG.optimize+      cg_params+      grad_tol+      x0+      (CG.MFunction mf)+      (CG.MGradient mg)+      (Just (CG.MCombined mc))++  let (success, msg) =+        case result of+          CG.ToleranceStatisfied      -> (True, "convergence tolerance satisfied")+          CG.FunctionChange           -> (True, "change in func <= feps*|f|")+          CG.MaxTotalIter             -> (False, "total iterations exceeded maxit")+          CG.NegativeSlope            -> (False, "slope always negative in line search")+          CG.MaxSecantIter            -> (False, "number secant iterations exceed nsecant")+          CG.NotDescent               -> (False, "search direction not a descent direction")+          CG.LineSearchFailsInitial   -> (False, "line search fails in initial interval")+          CG.LineSearchFailsBisection -> (False, "line search fails during bisection")+          CG.LineSearchFailsUpdate    -> (False, "line search fails during interval update")+          CG.DebugTol                 -> (False, "debugger is on and the function value increases")+          CG.FunctionValueNaN         -> (False, "function value became nan")+          CG.StartFunctionValueNaN    -> (False, "starting function value is nan")++  return $+    Result+    { resultSuccess = success+    , resultMessage = msg+    , resultSolution = x+    , resultValue = CG.finalValue stat+    , resultGrad = Nothing+    , resultHessian = Nothing+    , resultHessianInv = Nothing+    , resultStatistics =+        Statistics+        { totalIters = fromIntegral $ CG.totalIters stat+        , funcEvals = fromIntegral $ CG.funcEvals stat+        , gradEvals = fromIntegral $ CG.gradEvals stat+        , hessEvals = 0+        }+    }++#endif+++minimize_LBFGS :: HasGrad prob => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double))+minimize_LBFGS _params prob _ | not (isNothing (bounds prob)) = throwIO (UnsupportedProblem "LBFGS does not support bounds")+minimize_LBFGS _params prob _ | not (null (constraints prob)) = throwIO (UnsupportedProblem "LBFGS does not support constraints")+minimize_LBFGS params prob x0 = do+  evalCounter <- newIORef (0::Int)+  iterRef <- newIORef (0::Int)++  let lbfgsParams =+        LBFGS.LBFGSParameters+        { LBFGS.lbfgsPast = Nothing+        , LBFGS.lbfgsDelta = fromMaybe 0 $ paramsTol params+        , LBFGS.lbfgsLineSearch = LBFGS.DefaultLineSearch+        , LBFGS.lbfgsL1NormCoefficient = Nothing+        }++      instanceData :: ()+      instanceData = ()++      evalFun :: () -> VSM.IOVector CDouble -> VSM.IOVector CDouble -> CInt -> CDouble -> IO CDouble+      evalFun _inst xvec gvec _n _step = do+        modifyIORef' evalCounter (+1)+#if MIN_VERSION_vector(0,13,0)+        x <- VG.unsafeFreeze (VSM.unsafeCoerceMVector xvec :: VSM.IOVector Double)+        y <- grad'M prob x (VSM.unsafeCoerceMVector gvec :: VSM.IOVector Double)+#else+        x <- VG.unsafeFreeze (coerce xvec :: VSM.IOVector Double)+        y <- grad'M prob x (coerce gvec :: VSM.IOVector Double)+#endif+        return (coerce y)++      progressFun :: () -> VSM.IOVector CDouble -> VSM.IOVector CDouble -> CDouble -> CDouble -> CDouble -> CDouble -> CInt -> CInt -> CInt -> IO CInt+      progressFun _inst xvec _gvec _fx _xnorm _gnorm _step _n iter _nev = do+        writeIORef iterRef $! fromIntegral iter+        shouldStop <-+          case paramsCallback params of+            Nothing -> return False+            Just callback -> do+#if MIN_VERSION_vector(0,13,0)+              x <- VG.freeze (VSM.unsafeCoerceMVector xvec :: VSM.IOVector Double)+#else+              x <- VG.freeze (coerce xvec :: VSM.IOVector Double)+#endif+              callback x+        return $ if shouldStop then 1 else 0++  (result, x_) <- LBFGS.lbfgs lbfgsParams evalFun progressFun instanceData (VG.toList x0)+  let x = VG.fromList x_+      (success, msg) =+        case result of+          LBFGS.Success                -> (True,  "Success")+          LBFGS.Stop                   -> (True,  "Stop")+          LBFGS.AlreadyMinimized       -> (True,  "The initial variables already minimize the objective function.")+          LBFGS.UnknownError           -> (False, "Unknown error.")+          LBFGS.LogicError             -> (False, "Logic error.")+          LBFGS.OutOfMemory            -> (False, "Insufficient memory.")+          LBFGS.Canceled               -> (False, "The minimization process has been canceled.")+          LBFGS.InvalidN               -> (False, "Invalid number of variables specified.")+          LBFGS.InvalidNSSE            -> (False, "Invalid number of variables (for SSE) specified.")+          LBFGS.InvalidXSSE            -> (False, "The array x must be aligned to 16 (for SSE).")+          LBFGS.InvalidEpsilon         -> (False, "Invalid parameter lbfgs_parameter_t::epsilon specified.")+          LBFGS.InvalidTestPeriod      -> (False, "Invalid parameter lbfgs_parameter_t::past specified.")+          LBFGS.InvalidDelta           -> (False, "Invalid parameter lbfgs_parameter_t::delta specified.")+          LBFGS.InvalidLineSearch      -> (False, "Invalid parameter lbfgs_parameter_t::linesearch specified.")+          LBFGS.InvalidMinStep         -> (False, "Invalid parameter lbfgs_parameter_t::max_step specified.")+          LBFGS.InvalidMaxStep         -> (False, "Invalid parameter lbfgs_parameter_t::max_step specified.")+          LBFGS.InvalidFtol            -> (False, "Invalid parameter lbfgs_parameter_t::ftol specified.")+          LBFGS.InvalidWolfe           -> (False, "Invalid parameter lbfgs_parameter_t::wolfe specified.")+          LBFGS.InvalidGtol            -> (False, "Invalid parameter lbfgs_parameter_t::gtol specified.")+          LBFGS.InvalidXtol            -> (False, "Invalid parameter lbfgs_parameter_t::xtol specified.")+          LBFGS.InvalidMaxLineSearch   -> (False, "Invalid parameter lbfgs_parameter_t::max_linesearch specified.")+          LBFGS.InvalidOrthantwise     -> (False, "Invalid parameter lbfgs_parameter_t::orthantwise_c specified.")+          LBFGS.InvalidOrthantwiseStart-> (False, "Invalid parameter lbfgs_parameter_t::orthantwise_start specified.")+          LBFGS.InvalidOrthantwiseEnd  -> (False, "Invalid parameter lbfgs_parameter_t::orthantwise_end specified.")+          LBFGS.OutOfInterval          -> (False, "The line-search step went out of the interval of uncertainty.")+          LBFGS.IncorrectTMinMax       -> (False, "A logic error occurred; alternatively, the interval of uncertainty became too small.")+          LBFGS.RoundingError          -> (False, "A rounding error occurred; alternatively, no line-search step satisfies the sufficient decrease and curvature conditions.")+          LBFGS.MinimumStep            -> (False, "The line-search step became smaller than lbfgs_parameter_t::min_step.")+          LBFGS.MaximumStep            -> (False, "The line-search step became larger than lbfgs_parameter_t::max_step.")+          LBFGS.MaximumLineSearch      -> (False, "The line-search routine reaches the maximum number of evaluations.")+          LBFGS.MaximumIteration       -> (False, "The algorithm routine reaches the maximum number of iterations.")+          LBFGS.WidthTooSmall          -> (False, "Relative width of the interval of uncertainty is at most lbfgs_parameter_t::xtol.")+          LBFGS.InvalidParameters      -> (False, "A logic error (negative line-search step) occurred.")+          LBFGS.IncreaseGradient       -> (False, "The current search direction increases the objective function value.")++  nEvals <- readIORef evalCounter++  return $+    Result+    { resultSuccess = success+    , resultMessage = msg+    , resultSolution = x+    , resultValue = func prob x+    , resultGrad = Nothing+    , resultHessian = Nothing+    , resultHessianInv = Nothing+    , resultStatistics =+        Statistics+        { totalIters = undefined+        , funcEvals = nEvals + 1+        , gradEvals = nEvals + 1+        , hessEvals = 0+        }+    }+++minimize_Newton :: (HasGrad prob, HasHessian prob) => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double))+minimize_Newton _params prob _ | not (isNothing (bounds prob)) = throwIO (UnsupportedProblem "Newton does not support bounds")+minimize_Newton _params prob _ | not (null (constraints prob)) = throwIO (UnsupportedProblem "Newton does not support constraints")+minimize_Newton params prob x0 = do+  let tol = fromMaybe 1e-6 (paramsTol params)+      loop !x !y !g !h !n = do+        shouldStop <-+          case paramsCallback params of+            Just callback -> callback x+            Nothing -> return False+        if shouldStop then do+          return $+            Result+            { resultSuccess = False+            , resultMessage = "The minimization process has been canceled."+            , resultSolution = x+            , resultValue = y+            , resultGrad = Just g+            , resultHessian = Just h+            , resultHessianInv = Nothing+            , resultStatistics =+                Statistics+                { totalIters = n+                , funcEvals = n+                , gradEvals = n+                , hessEvals = n+                }+            }+        else do+          let p = h LA.<\> g+              x' = VG.zipWith (-) x p+          if LA.norm_Inf (VG.zipWith (-) x' x) > tol then do+            let (y', g') = grad' prob x'+                h' = hessian prob x'+            loop x' y' g' h' (n+1)+          else do+            return $+              Result+              { resultSuccess = True+              , resultMessage = "success"+              , resultSolution = x+              , resultValue = y+              , resultGrad = Just g+              , resultHessian = Just h+              , resultHessianInv = Nothing+              , resultStatistics =+                  Statistics+                  { totalIters = n+                  , funcEvals = n+                  , gradEvals = n+                  , hessEvals = n+                  }+              }+  let (y0, g0) = grad' prob x0+      h0 = hessian prob x0+  loop x0 y0 g0 h0 1++-- ------------------------------------------------------------------------++instance IsProblem (Vector Double -> Double) where+  func f = f++instance Optionally (HasGrad (Vector Double -> Double)) where+  optionalDict = Nothing++instance Optionally (HasHessian (Vector Double -> Double)) where+  optionalDict = Nothing++-- ------------------------------------------------------------------------++-- | Wrapper type for adding gradient function to a problem+data WithGrad prob = WithGrad prob (Vector Double -> Vector Double)++instance IsProblem prob => IsProblem (WithGrad prob) where+  func (WithGrad prob _g) = func prob+  bounds (WithGrad prob _g) = bounds prob+  constraints (WithGrad prob _g) = constraints prob++instance IsProblem prob => HasGrad (WithGrad prob) where+  grad (WithGrad _prob g) = g++instance HasHessian prob => HasHessian (WithGrad prob) where+  hessian (WithGrad prob _g) = hessian prob+  hessianProduct (WithGrad prob _g) = hessianProduct prob++instance IsProblem prob => Optionally (HasGrad (WithGrad prob)) where+  optionalDict = hasOptionalDict++instance Optionally (HasHessian prob) => Optionally (HasHessian (WithGrad prob)) where+  optionalDict =+    case optionalDict @(HasHessian prob) of+      Just Dict -> hasOptionalDict+      Nothing -> Nothing++-- ------------------------------------------------------------------------++-- | Wrapper type for adding hessian to a problem+data WithHessian prob = WithHessian prob (Vector Double -> Matrix Double)++instance IsProblem prob => IsProblem (WithHessian prob) where+  func (WithHessian prob _hess) = func prob+  bounds (WithHessian prob _hess) = bounds prob+  constraints (WithHessian prob _hess) = constraints prob++instance HasGrad prob => HasGrad (WithHessian prob) where+  grad (WithHessian prob _) = grad prob++instance IsProblem prob => HasHessian (WithHessian prob) where+  hessian (WithHessian _prob hess) = hess++instance Optionally (HasGrad prob) => Optionally (HasGrad (WithHessian prob)) where+  optionalDict =+    case optionalDict @(HasGrad prob) of+      Just Dict -> hasOptionalDict+      Nothing -> Nothing++instance IsProblem prob => Optionally (HasHessian (WithHessian prob)) where+  optionalDict = hasOptionalDict++-- ------------------------------------------------------------------------++-- | Wrapper type for adding bounds to a problem+data WithBounds prob = WithBounds prob (V.Vector (Double, Double))++instance IsProblem prob => IsProblem (WithBounds prob) where+  func (WithBounds prob _bounds) = func prob+  bounds (WithBounds _prob bounds) = Just bounds+  constraints (WithBounds prob _bounds) = constraints prob++instance HasGrad prob => HasGrad (WithBounds prob) where+  grad (WithBounds prob _bounds) = grad prob+  grad' (WithBounds prob _bounds) = grad' prob+  grad'M (WithBounds prob _bounds) = grad'M prob++instance HasHessian prob => HasHessian (WithBounds prob) where+  hessian (WithBounds prob _bounds) = hessian prob+  hessianProduct (WithBounds prob _bounds) = hessianProduct prob++instance Optionally (HasGrad prob) => Optionally (HasGrad (WithBounds prob)) where+  optionalDict =+    case optionalDict @(HasGrad prob) of+      Just Dict -> hasOptionalDict+      Nothing -> Nothing++instance Optionally (HasHessian prob) => Optionally (HasHessian (WithBounds prob)) where+  optionalDict =+    case optionalDict @(HasHessian prob) of+      Just Dict -> hasOptionalDict+      Nothing -> Nothing++-- ------------------------------------------------------------------------++-- | Wrapper type for adding constraints to a problem+data WithConstraints prob = WithConstraints prob [Constraint]++instance IsProblem prob => IsProblem (WithConstraints prob) where+  func (WithConstraints prob _constraints) = func prob+  bounds (WithConstraints prob _constraints) = bounds prob+  constraints (WithConstraints _prob constraints) = constraints++instance HasGrad prob => HasGrad (WithConstraints prob) where+  grad (WithConstraints prob _constraints) = grad prob+  grad' (WithConstraints prob _constraints) = grad' prob+  grad'M (WithConstraints prob _constraints) = grad'M prob++instance HasHessian prob => HasHessian (WithConstraints prob) where+  hessian (WithConstraints prob _constraints) = hessian prob+  hessianProduct (WithConstraints prob _constraints) = hessianProduct prob++instance Optionally (HasGrad prob) => Optionally (HasGrad (WithConstraints prob)) where+  optionalDict =+    case optionalDict @(HasGrad prob) of+      Just Dict -> hasOptionalDict+      Nothing -> Nothing++instance Optionally (HasHessian prob) => Optionally (HasHessian (WithConstraints prob)) where+  optionalDict =+    case optionalDict @(HasHessian prob) of+      Just Dict -> hasOptionalDict+      Nothing -> Nothing++-- ------------------------------------------------------------------------
+ test/IsClose.hs view
@@ -0,0 +1,153 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module IsClose+  (+  -- Tolerance type+    Tol (..)++  -- AllClose class+  , AllClose (..)+  , allCloseRawUnit+  , allCloseRawRealFrac+  , allCloseRawRealFloat++  -- * Re-exports+  , Default (..)++  -- * HUnit+  , assertAllClose+  ) where++import Data.Default.Class+import Data.List.NonEmpty (NonEmpty (..))+import Data.Map (Map)+import qualified Data.Map as Map+import Data.Monoid+import Data.Semigroup+import qualified Data.Vector as V+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Unboxed as VU+import GHC.Stack (HasCallStack)+import Test.HUnit+import Text.Printf++-- ------------------------------------------------------------------------++-- | Tolerance+--+-- Values @a@ and @b@ are considered /close/ if @abs (a - b) <= atol + rtol * abs b@.+data Tol a+  = Tol+  { rtol :: a -- ^ The relative tolerance parameter (default: @1e-05@)+  , atol :: a -- ^ The absolute tolerance parameter (default: @1e-08@)+  , equalNan :: Bool -- ^ Whether to compare NaN’s as equal (default: @False@)+  } deriving (Show)++instance RealFrac a => Default (Tol a) where+  def = Tol+    { rtol = 1e-05+    , atol = 1e-08+    , equalNan = False+    }++-- ------------------------------------------------------------------------++class Real r => AllClose r a where+  -- | Returns number of mismatches, number of elements, maximal absolute difference, and maximal relative difference.+  -- Returns @'Ap' 'Nothing'@ if given values are incomparable.+  allCloseRaw :: Tol r -> a -> a -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)++  -- | Returns 'True' if the two arrays are equal within the given tolerance; 'False' otherwise.+  allClose :: Tol r -> a -> a -> Bool+  allClose tol x y =+    case getAp (allCloseRaw tol x y) of+      Nothing -> False+      Just (Sum numMismatched, _, _, _) -> numMismatched == 0++allCloseRawRealFrac :: RealFrac r => Tol r -> r -> r -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)+allCloseRawRealFrac t a b = Ap $ Just $+  ( Sum $ if abs (a - b) <= atol t + rtol t * abs b then 0 else 1+  , Sum 1+  , Max (abs (a - b))+  , Max (abs (a - b) / abs b)+  )++allCloseRawRealFloat :: RealFloat r => Tol r -> r -> r -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)+allCloseRawRealFloat t a b+  | isNaN a /= isNaN b = Ap Nothing+  | otherwise = Ap $ Just $+      ( Sum $ if (equalNan t && isNaN a && isNaN b) || a == b || abs (a - b) <= atol t + rtol t * abs b then 0 else 1+      , Sum 1+      , Max (abs (a - b))+      , Max (abs (a - b) / abs b)+      )++allCloseRawUnit :: Num r => Ap Maybe (Sum Int, Sum Int, Max r, Max r)+allCloseRawUnit = Ap (Just (Sum 0, Sum 0, Max 0, Max 0))++instance AllClose Rational Rational where+  allCloseRaw = allCloseRawRealFrac++instance AllClose Double Double where+  allCloseRaw = allCloseRawRealFloat++instance (AllClose r a) => AllClose r (Maybe a) where+  allCloseRaw tol (Just a) (Just b) = allCloseRaw tol a b+  allCloseRaw _ Nothing Nothing = allCloseRawUnit+  allCloseRaw _ _ _ = Ap Nothing++instance (AllClose r v) => AllClose r [v] where+  allCloseRaw tol xs ys+    | length xs == length ys = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip xs ys])+    | otherwise = Ap Nothing++instance (Ord k, AllClose r v) => AllClose r (Map k v) where+  allCloseRaw tol m1 m2+    | Map.keys m1 == Map.keys m2 = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip (Map.elems m1) (Map.elems m2)])+    | otherwise = Ap Nothing++instance (AllClose r v) => AllClose r (V.Vector v) where+  allCloseRaw tol xs ys+    | VG.length xs == VG.length ys = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip (VG.toList xs) (VG.toList ys)])+    | otherwise = Ap Nothing++instance (AllClose r v, VS.Storable v) => AllClose r (VS.Vector v) where+  allCloseRaw tol xs ys+    | VG.length xs == VG.length ys = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip (VG.toList xs) (VG.toList ys)])+    | otherwise = Ap Nothing++instance (AllClose r v, VU.Unbox v) => AllClose r (VU.Vector v) where+  allCloseRaw tol xs ys+    | VG.length xs == VG.length ys = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip (VG.toList xs) (VG.toList ys)])+    | otherwise = Ap Nothing++-- ------------------------------------------------------------------------++-- | Assert that two objects are equal up to desired tolerance.+assertAllClose+  :: (HasCallStack, AllClose r a, Show r, Show a)+  => Tol r+  -> a -- ^ actual+  -> a -- ^ desired+  -> Assertion+assertAllClose tol a b =+  case getAp (allCloseRaw tol a b) of+    Nothing ->+      assertString $ unlines $ header ++ ["x and y nan location mismatch:"] ++ footer+    Just (Sum numMismatch, Sum numTotal, Max absDiff, Max relDiff)+      | numMismatch == 0 -> return ()+      | otherwise ->+          assertString $ unlines $+            header +++            [ printf "Mismatched elements: %d / %d (%f%%)" numMismatch numTotal (fromIntegral numMismatch * 100 / fromIntegral numTotal :: Double)+            , " Max absolute difference: " ++ show absDiff+            , " Max relative difference: " ++ show relDiff+            ] ++ footer+   where+     header, footer :: [String]+     header = [printf "Not equal to tolerance rtol=%s, atol=%s" (show (rtol tol)) (show (atol tol)), ""]+     footer = [" x: " ++ show a, " y: " ++ show b]++-- ------------------------------------------------------------------------
+ test/Spec.hs view
@@ -0,0 +1,30 @@+{-# LANGUAGE OverloadedLists #-}+import Test.Hspec++import Data.Vector.Storable (Vector)+import Numeric.Optimization+import IsClose+++main :: IO ()+main = hspec $ do+  describe "minimize" $ do+    context "when given rosenbrock function" $+      it "returns the global optimum" $ do+        result <- minimize LBFGS def (WithGrad rosenbrock rosenbrock') [-3,-4]+        resultSuccess result `shouldBe` True+        assertAllClose (def :: Tol Double) (resultSolution result) [1,1]+++-- https://en.wikipedia.org/wiki/Rosenbrock_function+rosenbrock :: Vector Double -> Double+rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)++rosenbrock' :: Vector Double -> Vector Double+rosenbrock' [x,y] =+  [ 2 * (1 - x) * (-1) + 100 * 2 * (y - sq x) * (-2) * x+  , 100 * 2 * (y - sq x)+  ]++sq :: Floating a => a -> a+sq x = x ** 2