diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. 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/rosenbrock.hs b/examples/rosenbrock.hs new file mode 100644 --- /dev/null +++ b/examples/rosenbrock.hs @@ -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 diff --git a/numeric-optimization.cabal b/numeric-optimization.cabal new file mode 100644 --- /dev/null +++ b/numeric-optimization.cabal @@ -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 +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 <masahiro.sakai@gmail.com> +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 diff --git a/src/Numeric/Optimization.hs b/src/Numeric/Optimization.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization.hs @@ -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] + -- + -- * [3] + | 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] + -- + -- * [2] + -- + -- * [3] + | 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 . + 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 + +-- ------------------------------------------------------------------------ diff --git a/test/IsClose.hs b/test/IsClose.hs new file mode 100644 --- /dev/null +++ b/test/IsClose.hs @@ -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] + +-- ------------------------------------------------------------------------ diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -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