{-# 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 provide unified 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](https://hackage.haskell.org/package/numeric-optimization-ad) or
-- [numeric-optimization-backprop](https://hackage.haskell.org/package/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.Applicative
import Control.Exception
import Control.Monad
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
#ifdef WITH_LBFGS
import qualified Numeric.LBFGS.Vector as LBFGS
import qualified Numeric.LBFGS.Raw as LBFGS (unCLBFGSResult, lbfgserrCanceled)
#endif
#ifdef WITH_CG_DESCENT
import qualified Numeric.Optimization.Algorithms.HagerZhang05 as CG
#endif
#ifdef WITH_LBFGSB
import qualified Numeric.LBFGSB as LBFGSB
import qualified Numeric.LBFGSB.Result as LBFGSB
#endif
import Numeric.Limits
import Numeric.LinearAlgebra (Matrix)
import qualified Numeric.LinearAlgebra as LA
import System.IO.Unsafe


-- | 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>
  | LBFGSB
    -- ^ Limited memory BFGS algorithm with bound constraints (L-BFGS-B) [1][2][3]
    --
    -- The implementation is provided by l-bfgs-b package [4]
    -- which is a bindign to L-BFGS-B fortran code [5].
    --
    -- * [1] R. H. Byrd, P. Lu and J. Nocedal. [A Limited Memory Algorithm for Bound Constrained Optimization](http://www.ece.northwestern.edu/~nocedal/PSfiles/limited.ps.gz), (1995), SIAM Journal on Scientific and Statistical Computing , 16, 5, pp. 1190-1208.
    --
    -- * [2] C. Zhu, R. H. Byrd and J. Nocedal. [L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization](http://www.ece.northwestern.edu/~nocedal/PSfiles/lbfgsb.ps.gz) (1997), ACM Transactions on Mathematical Software, Vol 23, Num. 4, pp. 550-560.
    --
    -- * [3] J. L. Morales and J. Nocedal. [L-BFGS-B: Remark on Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization](http://www.ece.northwestern.edu/~morales/PSfiles/acm-remark.pdf) (2011), ACM Transactions on Mathematical Software, Vol 38, Num. 7, pp. 1–4
    --
    -- * [4] <https://hackage.haskell.org/package/l-bfgs-b>
    --
    -- * [5] <http://users.iems.northwestern.edu/~nocedal/lbfgsb.html>
    --
    -- @since 0.1.1.0
  | Newton
    -- ^ Naïve implementation of Newton method in Haskell
    --
    -- This method requires both gradient and hessian.
  deriving (Method -> Method -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Method -> Method -> Bool
$c/= :: Method -> Method -> Bool
== :: Method -> Method -> Bool
$c== :: Method -> Method -> Bool
Eq, Eq Method
Method -> Method -> Bool
Method -> Method -> Ordering
Method -> Method -> Method
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Method -> Method -> Method
$cmin :: Method -> Method -> Method
max :: Method -> Method -> Method
$cmax :: Method -> Method -> Method
>= :: Method -> Method -> Bool
$c>= :: Method -> Method -> Bool
> :: Method -> Method -> Bool
$c> :: Method -> Method -> Bool
<= :: Method -> Method -> Bool
$c<= :: Method -> Method -> Bool
< :: Method -> Method -> Bool
$c< :: Method -> Method -> Bool
compare :: Method -> Method -> Ordering
$ccompare :: Method -> Method -> Ordering
Ord, Int -> Method
Method -> Int
Method -> [Method]
Method -> Method
Method -> Method -> [Method]
Method -> Method -> Method -> [Method]
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: Method -> Method -> Method -> [Method]
$cenumFromThenTo :: Method -> Method -> Method -> [Method]
enumFromTo :: Method -> Method -> [Method]
$cenumFromTo :: Method -> Method -> [Method]
enumFromThen :: Method -> Method -> [Method]
$cenumFromThen :: Method -> Method -> [Method]
enumFrom :: Method -> [Method]
$cenumFrom :: Method -> [Method]
fromEnum :: Method -> Int
$cfromEnum :: Method -> Int
toEnum :: Int -> Method
$ctoEnum :: Int -> Method
pred :: Method -> Method
$cpred :: Method -> Method
succ :: Method -> Method
$csucc :: Method -> Method
Enum, Int -> Method -> ShowS
[Method] -> ShowS
Method -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Method] -> ShowS
$cshowList :: [Method] -> ShowS
show :: Method -> String
$cshow :: Method -> String
showsPrec :: Int -> Method -> ShowS
$cshowsPrec :: Int -> Method -> ShowS
Show, Method
forall a. a -> a -> Bounded a
maxBound :: Method
$cmaxBound :: Method
minBound :: Method
$cminBound :: Method
Bounded)


-- | Whether a 'Method' is supported under the current environment.
isSupportedMethod :: Method -> Bool
#ifdef WITH_LBFGS
isSupportedMethod :: Method -> Bool
isSupportedMethod Method
LBFGS = Bool
True
#else
isSupportedMethod LBFGS = False
#endif
#ifdef WITH_CG_DESCENT
isSupportedMethod CGDescent = True
#else
isSupportedMethod Method
CGDescent = Bool
False
#endif
#ifdef WITH_LBFGSB
isSupportedMethod LBFGSB = True
#else
isSupportedMethod Method
LBFGSB = Bool
False
#endif
isSupportedMethod Method
Newton = Bool
True


-- | Parameters for optimization algorithms
--
-- TODO:
--
-- * Better way to pass algorithm specific parameters?
--
-- * Separate 'paramsCallback' from other more concrete serializeable parameters?
data Params a
  = Params
  { forall a. Params a -> Maybe (a -> IO Bool)
paramsCallback :: Maybe (a -> IO Bool)
    -- ^ If callback function returns @True@, the algorithm execution is terminated.
  , forall a. Params a -> Maybe Double
paramsTol :: Maybe Double
    -- ^ Tolerance for termination. When @tol@ is specified, the selected algorithm sets
    -- some relevant solver-specific tolerance(s) equal to @tol@.
    --
    -- If specified, this value is used as defaults for 'paramsFTol' and 'paramsGTol'.
  , forall a. Params a -> Maybe Double
paramsFTol :: Maybe Double
    -- ^ 'LBFGS' stops iteration when delta-based convergence test
    -- (see 'paramsPast') is enabled and the following condition is
    -- met:
    --
    -- \[
    --     \left|\frac{f' - f}{f}\right| < \mathrm{ftol},
    -- \]
    --
    -- where @f'@ is the objective value of @past@ ('paramsPast') iterations ago,
    -- and @f@ is the objective value of the current iteration.
    -- The default value is @1e-5@.
    --
    -- 'LBFGSB' stops iteration when the following condition is met:
    --
    -- \[
    --     \frac{f^k - f^{k+1}}{\mathrm{max}\{|f^k|,|f^{k+1}|,1\}} \le \mathrm{ftol}.
    -- \]
    --
    -- The default value is @1e7 * ('epsilon' :: Double) = 2.220446049250313e-9@.
    --
    -- @since 0.1.1.0
  , forall a. Params a -> Maybe Double
paramsGTol :: Maybe Double
    -- ^ 'LBFGSB' stops iteration when \(\mathrm{max}\{|\mathrm{pg}_i| \mid i = 1, \ldots, n\} \le \mathrm{gtol}\)
    -- where \(\mathrm{pg}_i\) is the i-th component of the projected gradient.
    --
    -- @since 0.1.1.0
  , forall a. Params a -> Maybe Int
paramsMaxIters :: Maybe Int
    -- ^ Maximum number of iterations.
    --
    -- Currently only 'LBFGSB', 'CGDescent', and 'Newton' uses this.
    --
    -- @since 0.1.1.0
  , forall a. Params a -> Maybe Int
paramsPast :: Maybe Int
    -- ^ Distance for delta-based convergence test in 'LBFGS'
    --
    -- This parameter determines the distance, in iterations, to compute
    -- the rate of decrease of the objective function. If the value of this
    -- parameter is @Nothing@, the library does not perform the delta-based
    -- convergence test. The default value is @Nothing@.
    --
    -- @since 0.1.1.0
  , forall a. Params a -> Maybe Int
paramsMaxCorrections :: Maybe Int
    -- ^ The maximum number of variable metric corrections used in 'LBFGSB'
    -- to define the limited memory matrix.
    --
    -- @since 0.1.1.0
  }

instance Default (Params a) where
  def :: Params a
def =
    Params
    { paramsCallback :: Maybe (a -> IO Bool)
paramsCallback = forall a. Maybe a
Nothing
    , paramsTol :: Maybe Double
paramsTol = forall a. Maybe a
Nothing
    , paramsFTol :: Maybe Double
paramsFTol = forall a. Maybe a
Nothing
    , paramsGTol :: Maybe Double
paramsGTol = forall a. Maybe a
Nothing
    , paramsMaxIters :: Maybe Int
paramsMaxIters = forall a. Maybe a
Nothing
    , paramsPast :: Maybe Int
paramsPast = forall a. Maybe a
Nothing
    , paramsMaxCorrections :: Maybe Int
paramsMaxCorrections = forall a. Maybe a
Nothing
    }

instance Contravariant Params where
  contramap :: forall a' a. (a' -> a) -> Params a -> Params a'
contramap a' -> a
f Params a
params =
    Params a
params
    { paramsCallback :: Maybe (a' -> IO Bool)
paramsCallback = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((forall b c a. (b -> c) -> (a -> b) -> a -> c
. a' -> a
f)) (forall a. Params a -> Maybe (a -> IO Bool)
paramsCallback Params a
params)
    }


-- | Optimization result
data Result a
  = Result
  { forall a. Result a -> Bool
resultSuccess :: Bool
    -- ^ Whether or not the optimizer exited successfully.
  , forall a. Result a -> String
resultMessage :: String
    -- ^ Description of the cause of the termination.
  , forall a. Result a -> a
resultSolution :: a
    -- ^ Solution
  , forall a. Result a -> Double
resultValue :: Double
    -- ^ Value of the function at the solution.
  , forall a. Result a -> Maybe a
resultGrad :: Maybe a
    -- ^ Gradient at the solution
  , forall a. Result a -> Maybe (Matrix Double)
resultHessian :: Maybe (Matrix Double)
    -- ^ Hessian at the solution; may be an approximation.
  , forall a. Result a -> Maybe (Matrix Double)
resultHessianInv :: Maybe (Matrix Double)
    -- ^ Inverse of Hessian at the solution; may be an approximation.
  , forall a. Result a -> Statistics
resultStatistics :: Statistics
    -- ^ Statistics of optimizaion process
  }
  deriving (Int -> Result a -> ShowS
forall a. Show a => Int -> Result a -> ShowS
forall a. Show a => [Result a] -> ShowS
forall a. Show a => Result a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Result a] -> ShowS
$cshowList :: forall a. Show a => [Result a] -> ShowS
show :: Result a -> String
$cshow :: forall a. Show a => Result a -> String
showsPrec :: Int -> Result a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Result a -> ShowS
Show)

instance Functor Result where
  fmap :: forall a b. (a -> b) -> Result a -> Result b
fmap a -> b
f Result a
result =
    Result a
result
    { resultSolution :: b
resultSolution = a -> b
f (forall a. Result a -> a
resultSolution Result a
result)
    , resultGrad :: Maybe b
resultGrad = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> b
f (forall a. Result a -> Maybe a
resultGrad Result a
result)
    }


-- | Statistics of optimizaion process
data Statistics
  = Statistics
  { Statistics -> Int
totalIters :: Int
    -- ^ Total number of iterations.
  , Statistics -> Int
funcEvals :: Int
    -- ^ Total number of function evaluations.
  , Statistics -> Int
gradEvals :: Int
    -- ^ Total number of gradient evaluations.
  , Statistics -> Int
hessianEvals :: Int
    -- ^ Total number of hessian evaluations.
  , Statistics -> Int
hessEvals :: Int
    -- ^ Total number of hessian evaluations.
  }
  deriving (Int -> Statistics -> ShowS
[Statistics] -> ShowS
Statistics -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Statistics] -> ShowS
$cshowList :: [Statistics] -> ShowS
show :: Statistics -> String
$cshow :: Statistics -> String
showsPrec :: Int -> Statistics -> ShowS
$cshowsPrec :: Int -> Statistics -> ShowS
Show)

{-# DEPRECATED hessEvals "Use hessianEvals instead" #-}


-- | The bad things that can happen when you use the library.
data OptimizationException
  = UnsupportedProblem String
  | UnsupportedMethod Method
  | GradUnavailable
  | HessianUnavailable
  deriving (Int -> OptimizationException -> ShowS
[OptimizationException] -> ShowS
OptimizationException -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [OptimizationException] -> ShowS
$cshowList :: [OptimizationException] -> ShowS
show :: OptimizationException -> String
$cshow :: OptimizationException -> String
showsPrec :: Int -> OptimizationException -> ShowS
$cshowsPrec :: Int -> OptimizationException -> ShowS
Show, OptimizationException -> OptimizationException -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: OptimizationException -> OptimizationException -> Bool
$c/= :: OptimizationException -> OptimizationException -> Bool
== :: OptimizationException -> OptimizationException -> Bool
$c== :: OptimizationException -> OptimizationException -> Bool
Eq)

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 require 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 prob
_ = forall a. Maybe a
Nothing

  -- | Constraints
  constraints :: prob -> [Constraint]
  constraints prob
_ = []

  {-# 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
prob = forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall prob.
HasGrad prob =>
prob -> Vector Double -> (Double, Vector Double)
grad' prob
prob

  -- | Pair of 'func' and 'grad'
  grad' :: prob -> Vector Double -> (Double, Vector Double)
  grad' prob
prob Vector Double
x = forall a. (forall s. ST s a) -> a
runST forall a b. (a -> b) -> a -> b
$ do
    MVector s Double
gret <- forall (m :: * -> *) (v :: * -> * -> *) a.
(HasCallStack, PrimMonad m, MVector v a) =>
Int -> m (v (PrimState m) a)
VGM.new (forall (v :: * -> *) a. Vector v a => v a -> Int
VG.length Vector Double
x)
    Double
y <- forall prob (m :: * -> *).
(HasGrad prob, PrimMonad m) =>
prob -> Vector Double -> MVector (PrimState m) Double -> m Double
grad'M prob
prob Vector Double
x MVector s Double
gret
    Vector Double
g <- forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
Mutable v (PrimState m) a -> m (v a)
VG.unsafeFreeze MVector s Double
gret
    forall (m :: * -> *) a. Monad m => a -> m a
return (Double
y, Vector Double
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
prob Vector Double
x MVector (PrimState m) Double
gvec = do
    let y :: Double
y = forall prob. IsProblem prob => prob -> Vector Double -> Double
func prob
prob Vector Double
x
    forall (m :: * -> *) (v :: * -> *) a b.
(Monad m, Vector v a) =>
(Int -> a -> m b) -> v a -> m ()
VG.imapM_ (forall (m :: * -> *) (v :: * -> * -> *) a.
(HasCallStack, PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> a -> m ()
VGM.write MVector (PrimState m) Double
gvec) (forall prob. HasGrad prob => prob -> Vector Double -> Vector Double
grad prob
prob Vector Double
x)
    forall (m :: * -> *) a. Monad m => a -> m a
return Double
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
prob Vector Double
x Vector Double
v = forall prob.
HasHessian prob =>
prob -> Vector Double -> Matrix Double
hessian prob
prob Vector Double
x forall t. Numeric t => Matrix t -> Vector t -> Vector t
LA.#> Vector Double
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 :: forall (c :: Constraint). c => Maybe (Dict c)
hasOptionalDict = forall a. a -> Maybe a
Just forall (a :: Constraint). a => Dict a
Dict


-- | Type of constraint
--
-- Currently, no constraints are supported.
data Constraint

-- | Bounds for unconstrained problems, i.e. (-∞,+∞).
boundsUnconstrained :: Int -> V.Vector (Double, Double)
boundsUnconstrained :: Int -> Vector (Double, Double)
boundsUnconstrained Int
n = forall a. Int -> a -> Vector a
V.replicate Int
n (-forall a. RealFloat a => a
infinity, forall a. RealFloat a => a
infinity)

-- | Whether all lower bounds are -∞ and all upper bounds are +∞.
isUnconstainedBounds :: V.Vector (Double, Double) -> Bool
isUnconstainedBounds :: Vector (Double, Double) -> Bool
isUnconstainedBounds = forall a. (a -> Bool) -> Vector a -> Bool
V.all forall {a} {a}. (RealFloat a, RealFloat a) => (a, a) -> Bool
p
  where
    p :: (a, a) -> Bool
p (a
lb, a
ub) = forall a. RealFloat a => a -> Bool
isInfinite a
lb Bool -> Bool -> Bool
&& a
lb forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
&& forall a. RealFloat a => a -> Bool
isInfinite a
ub Bool -> Bool -> Bool
&& a
ub forall a. Ord a => a -> a -> Bool
> a
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
#ifdef WITH_LBFGS
minimize :: forall prob.
(IsProblem prob, Optionally (HasGrad prob),
 Optionally (HasHessian prob)) =>
Method
-> Params (Vector Double)
-> prob
-> Vector Double
-> IO (Result (Vector Double))
minimize Method
LBFGS =
  case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasGrad prob) of
    Just Dict (HasGrad prob)
Dict -> forall prob.
HasGrad prob =>
Params (Vector Double)
-> prob -> Vector Double -> IO (Result (Vector Double))
minimize_LBFGS
    Maybe (Dict (HasGrad prob))
Nothing -> \Params (Vector Double)
_ prob
_ Vector Double
_ -> forall e a. Exception e => e -> IO a
throwIO OptimizationException
GradUnavailable
#endif
#ifdef WITH_LBFGSB
minimize LBFGSB =
  case optionalDict @(HasGrad prob) of
    Just Dict -> minimize_LBFGSB
    Nothing -> \_ _ _ -> throwIO GradUnavailable
#endif
minimize Method
Newton =
  case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasGrad prob) of
    Maybe (Dict (HasGrad prob))
Nothing -> \Params (Vector Double)
_ prob
_ Vector Double
_ -> forall e a. Exception e => e -> IO a
throwIO OptimizationException
GradUnavailable
    Just Dict (HasGrad prob)
Dict ->
      case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasHessian prob) of
        Maybe (Dict (HasHessian prob))
Nothing -> \Params (Vector Double)
_ prob
_ Vector Double
_ -> forall e a. Exception e => e -> IO a
throwIO OptimizationException
HessianUnavailable
        Just Dict (HasHessian prob)
Dict -> forall prob.
(HasGrad prob, HasHessian prob) =>
Params (Vector Double)
-> prob -> Vector Double -> IO (Result (Vector Double))
minimize_Newton
minimize Method
method = \Params (Vector Double)
_ prob
_ Vector Double
_ -> forall e a. Exception e => e -> IO a
throwIO (Method -> OptimizationException
UnsupportedMethod Method
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
        , CG.maxItersFac =
            case paramsMaxIters params of
              Nothing -> CG.maxItersFac CG.defaultParameters
              Just m -> fromIntegral m / fromIntegral (VG.length x0)
        }

      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
        , hessianEvals = 0
        }
    }

#endif


#ifdef WITH_LBFGS

minimize_LBFGS :: HasGrad prob => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double))
minimize_LBFGS :: forall prob.
HasGrad prob =>
Params (Vector Double)
-> prob -> Vector Double -> IO (Result (Vector Double))
minimize_LBFGS Params (Vector Double)
_params prob
prob Vector Double
_ | Bool -> Bool
not (forall a. Maybe a -> Bool
isNothing (forall prob.
IsProblem prob =>
prob -> Maybe (Vector (Double, Double))
bounds prob
prob)) = forall e a. Exception e => e -> IO a
throwIO (String -> OptimizationException
UnsupportedProblem String
"LBFGS does not support bounds")
minimize_LBFGS Params (Vector Double)
_params prob
prob Vector Double
_ | Bool -> Bool
not (forall (t :: * -> *) a. Foldable t => t a -> Bool
null (forall prob. IsProblem prob => prob -> [Constraint]
constraints prob
prob)) = forall e a. Exception e => e -> IO a
throwIO (String -> OptimizationException
UnsupportedProblem String
"LBFGS does not support constraints")
minimize_LBFGS Params (Vector Double)
params prob
prob Vector Double
x0 = do
  IORef Int
evalCounter <- forall a. a -> IO (IORef a)
newIORef (Int
0::Int)
  IORef Int
iterRef <- forall a. a -> IO (IORef a)
newIORef (Int
0::Int)

  let lbfgsParams :: LBFGSParameters
lbfgsParams =
        LBFGS.LBFGSParameters
        { lbfgsPast :: Maybe Int
LBFGS.lbfgsPast = forall a. Params a -> Maybe Int
paramsPast Params (Vector Double)
params
        , lbfgsDelta :: Double
LBFGS.lbfgsDelta = forall a. a -> Maybe a -> a
fromMaybe Double
1e-5 forall a b. (a -> b) -> a -> b
$ forall a. Params a -> Maybe Double
paramsFTol Params (Vector Double)
params forall (f :: * -> *) a. Alternative f => f a -> f a -> f a
<|> forall a. Params a -> Maybe Double
paramsTol Params (Vector Double)
params
        , lbfgsLineSearch :: LineSearchAlgorithm
LBFGS.lbfgsLineSearch = LineSearchAlgorithm
LBFGS.DefaultLineSearch
        , lbfgsL1NormCoefficient :: Maybe Double
LBFGS.lbfgsL1NormCoefficient = forall a. Maybe a
Nothing
        }

      instanceData :: ()
      instanceData :: ()
instanceData = ()

      evalFun :: () -> VSM.IOVector CDouble -> VSM.IOVector CDouble -> CInt -> CDouble -> IO CDouble
      evalFun :: ()
-> IOVector CDouble
-> IOVector CDouble
-> CInt
-> CDouble
-> IO CDouble
evalFun ()
_inst IOVector CDouble
xvec IOVector CDouble
gvec CInt
_n CDouble
_step = do
        forall a. IORef a -> (a -> a) -> IO ()
modifyIORef' IORef Int
evalCounter (forall a. Num a => a -> a -> a
+Int
1)
#if MIN_VERSION_vector(0,13,0)
        Vector Double
x <- forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
Mutable v (PrimState m) a -> m (v a)
VG.unsafeFreeze (forall a b s. Coercible a b => MVector s a -> MVector s b
VSM.unsafeCoerceMVector IOVector CDouble
xvec :: VSM.IOVector Double)
        Double
y <- forall prob (m :: * -> *).
(HasGrad prob, PrimMonad m) =>
prob -> Vector Double -> MVector (PrimState m) Double -> m Double
grad'M prob
prob Vector Double
x (forall a b s. Coercible a b => MVector s a -> MVector s b
VSM.unsafeCoerceMVector IOVector CDouble
gvec :: VSM.IOVector Double)
#else
        x <- VG.unsafeFreeze (coerce xvec :: VSM.IOVector Double)
        y <- grad'M prob x (coerce gvec :: VSM.IOVector Double)
#endif
        forall (m :: * -> *) a. Monad m => a -> m a
return (coerce :: forall a b. Coercible a b => a -> b
coerce Double
y)

      progressFun :: () -> VSM.IOVector CDouble -> VSM.IOVector CDouble -> CDouble -> CDouble -> CDouble -> CDouble -> CInt -> CInt -> CInt -> IO CInt
      progressFun :: ()
-> IOVector CDouble
-> IOVector CDouble
-> CDouble
-> CDouble
-> CDouble
-> CDouble
-> CInt
-> CInt
-> CInt
-> IO CInt
progressFun ()
_inst IOVector CDouble
xvec IOVector CDouble
_gvec CDouble
_fx CDouble
_xnorm CDouble
_gnorm CDouble
_step CInt
_n CInt
iter CInt
_nev = do
        forall a. IORef a -> a -> IO ()
writeIORef IORef Int
iterRef forall a b. (a -> b) -> a -> b
$! forall a b. (Integral a, Num b) => a -> b
fromIntegral CInt
iter
        Bool
shouldStop <-
          case forall a. Params a -> Maybe (a -> IO Bool)
paramsCallback Params (Vector Double)
params of
            Maybe (Vector Double -> IO Bool)
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return Bool
False
            Just Vector Double -> IO Bool
callback -> do
#if MIN_VERSION_vector(0,13,0)
              Vector Double
x <- forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
Mutable v (PrimState m) a -> m (v a)
VG.freeze (forall a b s. Coercible a b => MVector s a -> MVector s b
VSM.unsafeCoerceMVector IOVector CDouble
xvec :: VSM.IOVector Double)
#else
              x <- VG.freeze (coerce xvec :: VSM.IOVector Double)
#endif
              Vector Double -> IO Bool
callback Vector Double
x
        forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ if Bool
shouldStop then forall a b. (Integral a, Num b) => a -> b
fromIntegral (CLBFGSResult -> CInt
LBFGS.unCLBFGSResult CLBFGSResult
LBFGS.lbfgserrCanceled) else CInt
0

  (LBFGSResult
result, [Double]
x_) <- forall a.
LBFGSParameters
-> EvaluateFun a
-> ProgressFun a
-> a
-> [Double]
-> IO (LBFGSResult, [Double])
LBFGS.lbfgs LBFGSParameters
lbfgsParams ()
-> IOVector CDouble
-> IOVector CDouble
-> CInt
-> CDouble
-> IO CDouble
evalFun ()
-> IOVector CDouble
-> IOVector CDouble
-> CDouble
-> CDouble
-> CDouble
-> CDouble
-> CInt
-> CInt
-> CInt
-> IO CInt
progressFun ()
instanceData (forall (v :: * -> *) a. Vector v a => v a -> [a]
VG.toList Vector Double
x0)
  let x :: Vector Double
x = forall (v :: * -> *) a. Vector v a => [a] -> v a
VG.fromList [Double]
x_
      (Bool
success, String
msg) =
        case LBFGSResult
result of
          LBFGSResult
LBFGS.Success                -> (Bool
True,  String
"Success")
          LBFGSResult
LBFGS.Stop                   -> (Bool
True,  String
"Stop")
          LBFGSResult
LBFGS.AlreadyMinimized       -> (Bool
True,  String
"The initial variables already minimize the objective function.")
          LBFGSResult
LBFGS.UnknownError           -> (Bool
False, String
"Unknown error.")
          LBFGSResult
LBFGS.LogicError             -> (Bool
False, String
"Logic error.")
          LBFGSResult
LBFGS.OutOfMemory            -> (Bool
False, String
"Insufficient memory.")
          LBFGSResult
LBFGS.Canceled               -> (Bool
False, String
"The minimization process has been canceled.")
          LBFGSResult
LBFGS.InvalidN               -> (Bool
False, String
"Invalid number of variables specified.")
          LBFGSResult
LBFGS.InvalidNSSE            -> (Bool
False, String
"Invalid number of variables (for SSE) specified.")
          LBFGSResult
LBFGS.InvalidXSSE            -> (Bool
False, String
"The array x must be aligned to 16 (for SSE).")
          LBFGSResult
LBFGS.InvalidEpsilon         -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::epsilon specified.")
          LBFGSResult
LBFGS.InvalidTestPeriod      -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::past specified.")
          LBFGSResult
LBFGS.InvalidDelta           -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::delta specified.")
          LBFGSResult
LBFGS.InvalidLineSearch      -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::linesearch specified.")
          LBFGSResult
LBFGS.InvalidMinStep         -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::max_step specified.")
          LBFGSResult
LBFGS.InvalidMaxStep         -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::max_step specified.")
          LBFGSResult
LBFGS.InvalidFtol            -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::ftol specified.")
          LBFGSResult
LBFGS.InvalidWolfe           -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::wolfe specified.")
          LBFGSResult
LBFGS.InvalidGtol            -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::gtol specified.")
          LBFGSResult
LBFGS.InvalidXtol            -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::xtol specified.")
          LBFGSResult
LBFGS.InvalidMaxLineSearch   -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::max_linesearch specified.")
          LBFGSResult
LBFGS.InvalidOrthantwise     -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::orthantwise_c specified.")
          LBFGSResult
LBFGS.InvalidOrthantwiseStart-> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::orthantwise_start specified.")
          LBFGSResult
LBFGS.InvalidOrthantwiseEnd  -> (Bool
False, String
"Invalid parameter lbfgs_parameter_t::orthantwise_end specified.")
          LBFGSResult
LBFGS.OutOfInterval          -> (Bool
False, String
"The line-search step went out of the interval of uncertainty.")
          LBFGSResult
LBFGS.IncorrectTMinMax       -> (Bool
False, String
"A logic error occurred; alternatively, the interval of uncertainty became too small.")
          LBFGSResult
LBFGS.RoundingError          -> (Bool
False, String
"A rounding error occurred; alternatively, no line-search step satisfies the sufficient decrease and curvature conditions.")
          LBFGSResult
LBFGS.MinimumStep            -> (Bool
False, String
"The line-search step became smaller than lbfgs_parameter_t::min_step.")
          LBFGSResult
LBFGS.MaximumStep            -> (Bool
False, String
"The line-search step became larger than lbfgs_parameter_t::max_step.")
          LBFGSResult
LBFGS.MaximumLineSearch      -> (Bool
False, String
"The line-search routine reaches the maximum number of evaluations.")
          LBFGSResult
LBFGS.MaximumIteration       -> (Bool
False, String
"The algorithm routine reaches the maximum number of iterations.")
          LBFGSResult
LBFGS.WidthTooSmall          -> (Bool
False, String
"Relative width of the interval of uncertainty is at most lbfgs_parameter_t::xtol.")
          LBFGSResult
LBFGS.InvalidParameters      -> (Bool
False, String
"A logic error (negative line-search step) occurred.")
          LBFGSResult
LBFGS.IncreaseGradient       -> (Bool
False, String
"The current search direction increases the objective function value.")

  Int
iters <- forall a. IORef a -> IO a
readIORef IORef Int
iterRef
  Int
nEvals <- forall a. IORef a -> IO a
readIORef IORef Int
evalCounter

  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$
    Result
    { resultSuccess :: Bool
resultSuccess = Bool
success
    , resultMessage :: String
resultMessage = String
msg
    , resultSolution :: Vector Double
resultSolution = Vector Double
x
    , resultValue :: Double
resultValue = forall prob. IsProblem prob => prob -> Vector Double -> Double
func prob
prob Vector Double
x
    , resultGrad :: Maybe (Vector Double)
resultGrad = forall a. Maybe a
Nothing
    , resultHessian :: Maybe (Matrix Double)
resultHessian = forall a. Maybe a
Nothing
    , resultHessianInv :: Maybe (Matrix Double)
resultHessianInv = forall a. Maybe a
Nothing
    , resultStatistics :: Statistics
resultStatistics =
        Statistics
        { totalIters :: Int
totalIters = Int
iters
        , funcEvals :: Int
funcEvals = Int
nEvals forall a. Num a => a -> a -> a
+ Int
1  -- +1 is for computing 'resultValue'
        , gradEvals :: Int
gradEvals = Int
nEvals
        , hessEvals :: Int
hessEvals = Int
0
        , hessianEvals :: Int
hessianEvals = Int
0
        }
    }

#endif


#ifdef WITH_LBFGSB

minimize_LBFGSB :: HasGrad prob => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double))
minimize_LBFGSB _params prob _ | not (null (constraints prob)) = throwIO (UnsupportedProblem "LBFGSB does not support constraints")
minimize_LBFGSB params prob x0 = do
  funcEvalRef <- newIORef (0::Int)
  gradEvalRef <- newIORef (0::Int)

  let bounds' =
        case bounds prob of
          Nothing -> []
          Just vec -> map convertB (VG.toList vec)
      convertB (lb, ub) =
        ( if isInfinite lb && lb < 0
          then Nothing
          else Just lb
        , if isInfinite ub && ub > 0
          then Nothing
          else Just ub
        )
      func' x = unsafePerformIO $ do
        modifyIORef' funcEvalRef (+1)
        evaluate (func prob x)
      grad' x = unsafePerformIO $ do
        modifyIORef' gradEvalRef (+1)
        evaluate (grad prob x)

  let -- | @m@: The maximum number of variable metric corrections used
      -- to define the limited memory matrix. /Suggestion:/ @5@.
      m = fromMaybe 5 (paramsMaxCorrections params)

      -- | @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@.
      factr = fromMaybe 1e7 $ (/ epsilon) <$> (paramsFTol params <|> paramsTol params)

      -- ^ @pgtol@: Iteration stops when the largest component of the projected
      -- gradient is smaller than @pgtol@. Must be @>=0@. /Suggestion:/ @1e-5@.
      pgtol = fromMaybe 1e-5 $ paramsGTol params <|> paramsTol params

      -- | @'Just' steps@ means the minimization is aborted if it has not converged after
      -- @steps>0@ iterations. 'Nothing' signifies no limit.
      steps = paramsMaxIters params

  result <- evaluate $ LBFGSB.minimize m factr pgtol steps bounds' x0 func' grad'

  let x = LBFGSB.solution result
      (success, msg) =
         case LBFGSB.stopReason result of
           LBFGSB.Converged -> (True, "The solution converged.")
           LBFGSB.StepCount -> (False, "The number of steps exceeded the user's request.")
           LBFGSB.Other msg -> (False, msg)

  funcEvals <- readIORef funcEvalRef
  gradEvals <- readIORef gradEvalRef

  return $
    Result
    { resultSuccess = success
    , resultMessage = msg
    , resultSolution = x
    , resultValue = func prob x
    , resultGrad = Nothing
    , resultHessian = Nothing
    , resultHessianInv = Nothing
    , resultStatistics =
        Statistics
        { totalIters = length (LBFGSB.backtrace result)
        , funcEvals = funcEvals
        , gradEvals = gradEvals
        , hessEvals = 0
        , hessianEvals = 0
        }
    }

#endif


minimize_Newton :: (HasGrad prob, HasHessian prob) => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double))
minimize_Newton :: forall prob.
(HasGrad prob, HasHessian prob) =>
Params (Vector Double)
-> prob -> Vector Double -> IO (Result (Vector Double))
minimize_Newton Params (Vector Double)
_params prob
prob Vector Double
_ | Bool -> Bool
not (forall a. Maybe a -> Bool
isNothing (forall prob.
IsProblem prob =>
prob -> Maybe (Vector (Double, Double))
bounds prob
prob)) = forall e a. Exception e => e -> IO a
throwIO (String -> OptimizationException
UnsupportedProblem String
"Newton does not support bounds")
minimize_Newton Params (Vector Double)
_params prob
prob Vector Double
_ | Bool -> Bool
not (forall (t :: * -> *) a. Foldable t => t a -> Bool
null (forall prob. IsProblem prob => prob -> [Constraint]
constraints prob
prob)) = forall e a. Exception e => e -> IO a
throwIO (String -> OptimizationException
UnsupportedProblem String
"Newton does not support constraints")
minimize_Newton Params (Vector Double)
params prob
prob Vector Double
x0 = do
  let tol :: Double
tol = forall a. a -> Maybe a -> a
fromMaybe Double
1e-6 (forall a. Params a -> Maybe Double
paramsTol Params (Vector Double)
params)

      loop :: Vector Double
-> Double
-> Vector Double
-> Matrix Double
-> Int
-> IO (Result (Vector Double))
loop !Vector Double
x !Double
y !Vector Double
g !Matrix Double
h !Int
iter = do
        Maybe String
shouldStop <- forall (t :: * -> *) (m :: * -> *) a.
(Foldable t, MonadPlus m) =>
t (m a) -> m a
msum forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence
          [ forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a b. (a -> b) -> a -> b
$ case forall a. Params a -> Maybe Int
paramsMaxIters Params (Vector Double)
params of
              Just Int
maxIter | Int
maxIter forall a. Ord a => a -> a -> Bool
<= Int
iter -> forall a. a -> Maybe a
Just String
"maximum number of iterations reached"
              Maybe Int
_ -> forall a. Maybe a
Nothing
          , case forall a. Params a -> Maybe (a -> IO Bool)
paramsCallback Params (Vector Double)
params of
              Maybe (Vector Double -> IO Bool)
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Maybe a
Nothing
              Just Vector Double -> IO Bool
callback -> do
                Bool
flag <- Vector Double -> IO Bool
callback Vector Double
x
                forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ if Bool
flag then forall a. a -> Maybe a
Just String
"The minimization process has been canceled." else forall a. Maybe a
Nothing
          ]
        case Maybe String
shouldStop of
          Just String
reason ->
            forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$
              Result
              { resultSuccess :: Bool
resultSuccess = Bool
False
              , resultMessage :: String
resultMessage = String
reason
              , resultSolution :: Vector Double
resultSolution = Vector Double
x
              , resultValue :: Double
resultValue = Double
y
              , resultGrad :: Maybe (Vector Double)
resultGrad = forall a. a -> Maybe a
Just Vector Double
g
              , resultHessian :: Maybe (Matrix Double)
resultHessian = forall a. a -> Maybe a
Just Matrix Double
h
              , resultHessianInv :: Maybe (Matrix Double)
resultHessianInv = forall a. Maybe a
Nothing
              , resultStatistics :: Statistics
resultStatistics =
                  Statistics
                  { totalIters :: Int
totalIters = Int
iter
                  , funcEvals :: Int
funcEvals = Int
iter forall a. Num a => a -> a -> a
+ Int
1
                  , gradEvals :: Int
gradEvals = Int
iter forall a. Num a => a -> a -> a
+ Int
1
                  , hessEvals :: Int
hessEvals = Int
iter forall a. Num a => a -> a -> a
+ Int
1
                  , hessianEvals :: Int
hessianEvals = Int
iter forall a. Num a => a -> a -> a
+ Int
1
                  }
              }
          Maybe String
Nothing -> do
            let p :: Vector Double
p = Matrix Double
h forall (c :: * -> *) t.
(LSDiv c, Field t) =>
Matrix t -> c t -> c t
LA.<\> Vector Double
g
                x' :: Vector Double
x' = forall (v :: * -> *) a b c.
(Vector v a, Vector v b, Vector v c) =>
(a -> b -> c) -> v a -> v b -> v c
VG.zipWith (-) Vector Double
x Vector Double
p
            if forall a. Normed a => a -> Double
LA.norm_Inf (forall (v :: * -> *) a b c.
(Vector v a, Vector v b, Vector v c) =>
(a -> b -> c) -> v a -> v b -> v c
VG.zipWith (-) Vector Double
x' Vector Double
x) forall a. Ord a => a -> a -> Bool
> Double
tol then do
              let (Double
y', Vector Double
g') = forall prob.
HasGrad prob =>
prob -> Vector Double -> (Double, Vector Double)
grad' prob
prob Vector Double
x'
                  h' :: Matrix Double
h' = forall prob.
HasHessian prob =>
prob -> Vector Double -> Matrix Double
hessian prob
prob Vector Double
x'
              Vector Double
-> Double
-> Vector Double
-> Matrix Double
-> Int
-> IO (Result (Vector Double))
loop Vector Double
x' Double
y' Vector Double
g' Matrix Double
h' (Int
iter forall a. Num a => a -> a -> a
+ Int
1)
            else do
              forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$
                Result
                { resultSuccess :: Bool
resultSuccess = Bool
True
                , resultMessage :: String
resultMessage = String
"success"
                , resultSolution :: Vector Double
resultSolution = Vector Double
x
                , resultValue :: Double
resultValue = Double
y
                , resultGrad :: Maybe (Vector Double)
resultGrad = forall a. a -> Maybe a
Just Vector Double
g
                , resultHessian :: Maybe (Matrix Double)
resultHessian = forall a. a -> Maybe a
Just Matrix Double
h
                , resultHessianInv :: Maybe (Matrix Double)
resultHessianInv = forall a. Maybe a
Nothing
                , resultStatistics :: Statistics
resultStatistics =
                    Statistics
                    { totalIters :: Int
totalIters = Int
iter
                    , funcEvals :: Int
funcEvals = Int
iter forall a. Num a => a -> a -> a
+ Int
1
                    , gradEvals :: Int
gradEvals = Int
iter forall a. Num a => a -> a -> a
+ Int
1
                    , hessEvals :: Int
hessEvals = Int
iter forall a. Num a => a -> a -> a
+ Int
1
                    , hessianEvals :: Int
hessianEvals = Int
iter forall a. Num a => a -> a -> a
+ Int
1
                    }
                }

  let (Double
y0, Vector Double
g0) = forall prob.
HasGrad prob =>
prob -> Vector Double -> (Double, Vector Double)
grad' prob
prob Vector Double
x0
      h0 :: Matrix Double
h0 = forall prob.
HasHessian prob =>
prob -> Vector Double -> Matrix Double
hessian prob
prob Vector Double
x0
  Vector Double
-> Double
-> Vector Double
-> Matrix Double
-> Int
-> IO (Result (Vector Double))
loop Vector Double
x0 Double
y0 Vector Double
g0 Matrix Double
h0 Int
0

-- ------------------------------------------------------------------------

instance IsProblem (Vector Double -> Double) where
  func :: (Vector Double -> Double) -> Vector Double -> Double
func Vector Double -> Double
f = Vector Double -> Double
f

instance Optionally (HasGrad (Vector Double -> Double)) where
  optionalDict :: Maybe (Dict (HasGrad (Vector Double -> Double)))
optionalDict = forall a. Maybe a
Nothing

instance Optionally (HasHessian (Vector Double -> Double)) where
  optionalDict :: Maybe (Dict (HasHessian (Vector Double -> Double)))
optionalDict = forall a. Maybe a
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 -> Vector Double -> Double
func (WithGrad prob
prob Vector Double -> Vector Double
_g) = forall prob. IsProblem prob => prob -> Vector Double -> Double
func prob
prob
  bounds :: WithGrad prob -> Maybe (Vector (Double, Double))
bounds (WithGrad prob
prob Vector Double -> Vector Double
_g) = forall prob.
IsProblem prob =>
prob -> Maybe (Vector (Double, Double))
bounds prob
prob
  constraints :: WithGrad prob -> [Constraint]
constraints (WithGrad prob
prob Vector Double -> Vector Double
_g) = forall prob. IsProblem prob => prob -> [Constraint]
constraints prob
prob

instance IsProblem prob => HasGrad (WithGrad prob) where
  grad :: WithGrad prob -> Vector Double -> Vector Double
grad (WithGrad prob
_prob Vector Double -> Vector Double
g) = Vector Double -> Vector Double
g

instance HasHessian prob => HasHessian (WithGrad prob) where
  hessian :: WithGrad prob -> Vector Double -> Matrix Double
hessian (WithGrad prob
prob Vector Double -> Vector Double
_g) = forall prob.
HasHessian prob =>
prob -> Vector Double -> Matrix Double
hessian prob
prob
  hessianProduct :: WithGrad prob -> Vector Double -> Vector Double -> Vector Double
hessianProduct (WithGrad prob
prob Vector Double -> Vector Double
_g) = forall prob.
HasHessian prob =>
prob -> Vector Double -> Vector Double -> Vector Double
hessianProduct prob
prob

instance IsProblem prob => Optionally (HasGrad (WithGrad prob)) where
  optionalDict :: Maybe (Dict (HasGrad (WithGrad prob)))
optionalDict = forall (c :: Constraint). c => Maybe (Dict c)
hasOptionalDict

instance Optionally (HasHessian prob) => Optionally (HasHessian (WithGrad prob)) where
  optionalDict :: Maybe (Dict (HasHessian (WithGrad prob)))
optionalDict =
    case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasHessian prob) of
      Just Dict (HasHessian prob)
Dict -> forall (c :: Constraint). c => Maybe (Dict c)
hasOptionalDict
      Maybe (Dict (HasHessian prob))
Nothing -> forall a. Maybe a
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 -> Vector Double -> Double
func (WithHessian prob
prob Vector Double -> Matrix Double
_hess) = forall prob. IsProblem prob => prob -> Vector Double -> Double
func prob
prob
  bounds :: WithHessian prob -> Maybe (Vector (Double, Double))
bounds (WithHessian prob
prob Vector Double -> Matrix Double
_hess) = forall prob.
IsProblem prob =>
prob -> Maybe (Vector (Double, Double))
bounds prob
prob
  constraints :: WithHessian prob -> [Constraint]
constraints (WithHessian prob
prob Vector Double -> Matrix Double
_hess) = forall prob. IsProblem prob => prob -> [Constraint]
constraints prob
prob

instance HasGrad prob => HasGrad (WithHessian prob) where
  grad :: WithHessian prob -> Vector Double -> Vector Double
grad (WithHessian prob
prob Vector Double -> Matrix Double
_) = forall prob. HasGrad prob => prob -> Vector Double -> Vector Double
grad prob
prob

instance IsProblem prob => HasHessian (WithHessian prob) where
  hessian :: WithHessian prob -> Vector Double -> Matrix Double
hessian (WithHessian prob
_prob Vector Double -> Matrix Double
hess) = Vector Double -> Matrix Double
hess

instance Optionally (HasGrad prob) => Optionally (HasGrad (WithHessian prob)) where
  optionalDict :: Maybe (Dict (HasGrad (WithHessian prob)))
optionalDict =
    case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasGrad prob) of
      Just Dict (HasGrad prob)
Dict -> forall (c :: Constraint). c => Maybe (Dict c)
hasOptionalDict
      Maybe (Dict (HasGrad prob))
Nothing -> forall a. Maybe a
Nothing

instance IsProblem prob => Optionally (HasHessian (WithHessian prob)) where
  optionalDict :: Maybe (Dict (HasHessian (WithHessian prob)))
optionalDict = forall (c :: Constraint). c => Maybe (Dict c)
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 -> Vector Double -> Double
func (WithBounds prob
prob Vector (Double, Double)
_bounds) = forall prob. IsProblem prob => prob -> Vector Double -> Double
func prob
prob
  bounds :: WithBounds prob -> Maybe (Vector (Double, Double))
bounds (WithBounds prob
_prob Vector (Double, Double)
bounds) = forall a. a -> Maybe a
Just Vector (Double, Double)
bounds
  constraints :: WithBounds prob -> [Constraint]
constraints (WithBounds prob
prob Vector (Double, Double)
_bounds) = forall prob. IsProblem prob => prob -> [Constraint]
constraints prob
prob

instance HasGrad prob => HasGrad (WithBounds prob) where
  grad :: WithBounds prob -> Vector Double -> Vector Double
grad (WithBounds prob
prob Vector (Double, Double)
_bounds) = forall prob. HasGrad prob => prob -> Vector Double -> Vector Double
grad prob
prob
  grad' :: WithBounds prob -> Vector Double -> (Double, Vector Double)
grad' (WithBounds prob
prob Vector (Double, Double)
_bounds) = forall prob.
HasGrad prob =>
prob -> Vector Double -> (Double, Vector Double)
grad' prob
prob
  grad'M :: forall (m :: * -> *).
PrimMonad m =>
WithBounds prob
-> Vector Double -> MVector (PrimState m) Double -> m Double
grad'M (WithBounds prob
prob Vector (Double, Double)
_bounds) = forall prob (m :: * -> *).
(HasGrad prob, PrimMonad m) =>
prob -> Vector Double -> MVector (PrimState m) Double -> m Double
grad'M prob
prob

instance HasHessian prob => HasHessian (WithBounds prob) where
  hessian :: WithBounds prob -> Vector Double -> Matrix Double
hessian (WithBounds prob
prob Vector (Double, Double)
_bounds) = forall prob.
HasHessian prob =>
prob -> Vector Double -> Matrix Double
hessian prob
prob
  hessianProduct :: WithBounds prob -> Vector Double -> Vector Double -> Vector Double
hessianProduct (WithBounds prob
prob Vector (Double, Double)
_bounds) = forall prob.
HasHessian prob =>
prob -> Vector Double -> Vector Double -> Vector Double
hessianProduct prob
prob

instance Optionally (HasGrad prob) => Optionally (HasGrad (WithBounds prob)) where
  optionalDict :: Maybe (Dict (HasGrad (WithBounds prob)))
optionalDict =
    case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasGrad prob) of
      Just Dict (HasGrad prob)
Dict -> forall (c :: Constraint). c => Maybe (Dict c)
hasOptionalDict
      Maybe (Dict (HasGrad prob))
Nothing -> forall a. Maybe a
Nothing

instance Optionally (HasHessian prob) => Optionally (HasHessian (WithBounds prob)) where
  optionalDict :: Maybe (Dict (HasHessian (WithBounds prob)))
optionalDict =
    case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasHessian prob) of
      Just Dict (HasHessian prob)
Dict -> forall (c :: Constraint). c => Maybe (Dict c)
hasOptionalDict
      Maybe (Dict (HasHessian prob))
Nothing -> forall a. Maybe a
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 -> Vector Double -> Double
func (WithConstraints prob
prob [Constraint]
_constraints) = forall prob. IsProblem prob => prob -> Vector Double -> Double
func prob
prob
  bounds :: WithConstraints prob -> Maybe (Vector (Double, Double))
bounds (WithConstraints prob
prob [Constraint]
_constraints) = forall prob.
IsProblem prob =>
prob -> Maybe (Vector (Double, Double))
bounds prob
prob
  constraints :: WithConstraints prob -> [Constraint]
constraints (WithConstraints prob
_prob [Constraint]
constraints) = [Constraint]
constraints

instance HasGrad prob => HasGrad (WithConstraints prob) where
  grad :: WithConstraints prob -> Vector Double -> Vector Double
grad (WithConstraints prob
prob [Constraint]
_constraints) = forall prob. HasGrad prob => prob -> Vector Double -> Vector Double
grad prob
prob
  grad' :: WithConstraints prob -> Vector Double -> (Double, Vector Double)
grad' (WithConstraints prob
prob [Constraint]
_constraints) = forall prob.
HasGrad prob =>
prob -> Vector Double -> (Double, Vector Double)
grad' prob
prob
  grad'M :: forall (m :: * -> *).
PrimMonad m =>
WithConstraints prob
-> Vector Double -> MVector (PrimState m) Double -> m Double
grad'M (WithConstraints prob
prob [Constraint]
_constraints) = forall prob (m :: * -> *).
(HasGrad prob, PrimMonad m) =>
prob -> Vector Double -> MVector (PrimState m) Double -> m Double
grad'M prob
prob

instance HasHessian prob => HasHessian (WithConstraints prob) where
  hessian :: WithConstraints prob -> Vector Double -> Matrix Double
hessian (WithConstraints prob
prob [Constraint]
_constraints) = forall prob.
HasHessian prob =>
prob -> Vector Double -> Matrix Double
hessian prob
prob
  hessianProduct :: WithConstraints prob
-> Vector Double -> Vector Double -> Vector Double
hessianProduct (WithConstraints prob
prob [Constraint]
_constraints) = forall prob.
HasHessian prob =>
prob -> Vector Double -> Vector Double -> Vector Double
hessianProduct prob
prob

instance Optionally (HasGrad prob) => Optionally (HasGrad (WithConstraints prob)) where
  optionalDict :: Maybe (Dict (HasGrad (WithConstraints prob)))
optionalDict =
    case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasGrad prob) of
      Just Dict (HasGrad prob)
Dict -> forall (c :: Constraint). c => Maybe (Dict c)
hasOptionalDict
      Maybe (Dict (HasGrad prob))
Nothing -> forall a. Maybe a
Nothing

instance Optionally (HasHessian prob) => Optionally (HasHessian (WithConstraints prob)) where
  optionalDict :: Maybe (Dict (HasHessian (WithConstraints prob)))
optionalDict =
    case forall (c :: Constraint). Optionally c => Maybe (Dict c)
optionalDict @(HasHessian prob) of
      Just Dict (HasHessian prob)
Dict -> forall (c :: Constraint). c => Maybe (Dict c)
hasOptionalDict
      Maybe (Dict (HasHessian prob))
Nothing -> forall a. Maybe a
Nothing

-- ------------------------------------------------------------------------