diff --git a/Ipopt.hs b/Ipopt.hs --- a/Ipopt.hs +++ b/Ipopt.hs @@ -18,6 +18,8 @@ ipopts, setIpoptProblemScaling, openIpoptOutputFile, + setIntermediateCallback, + IntermediateCB, -- * low-level bits still needed IpOptSolved(..), diff --git a/Ipopt/AnyRF.hs b/Ipopt/AnyRF.hs --- a/Ipopt/AnyRF.hs +++ b/Ipopt/AnyRF.hs @@ -1,4 +1,5 @@ -{-# LANGUAGE ConstraintKinds, RankNTypes, TypeFamilies, FlexibleInstances #-} +{-# LANGUAGE ConstraintKinds, RankNTypes, TypeFamilies, FlexibleInstances, + UndecidableInstances #-} {- | Description: representing functions that can be differentiated The 'AnyRF' wrapper holds functions that can be used @@ -13,6 +14,18 @@ defined in "Ipopt.NLP" (also exported by "Ipopt"). Directly using the constructor is another option: @AnyRF $ Identity . V.sum@, calculates the sum of all variables in the problem. + +convergence can sometimes be improved by exposing additional +variables/derivatives to the solver. IE. instead of maximizing f(u) where +we internally calculate x(u), maximize f(u,y) with another constraint that +y = x(u). Albersmeyer 2010 SIAM. Also generally known in process optimization +(ie. "equation-oriented" which has many variables works, while "sequential-modular" +modes that have very few variables left does not work) + +AnyRF could be used to generate some C and feed it to something like casadi, +though this has to happen a runtime unless there's a clean way to handle constant +parameters such as (fromIntegral 3) + -} module Ipopt.AnyRF where @@ -24,9 +37,15 @@ import Data.VectorSpace (VectorSpace, Scalar) import qualified Numeric.AD as AD -import qualified Numeric.AD.Types as AD -import qualified Numeric.AD.Internal.Classes as AD +import qualified Numeric.AD.Mode as AD +import qualified Numeric.AD.Internal.Forward as AD +import qualified Numeric.AD.Internal.Identity as AD +import qualified Numeric.AD.Internal.Kahn as AD +import qualified Numeric.AD.Internal.On as AD +import qualified Numeric.AD.Internal.Reverse as AD +import qualified Numeric.AD.Internal.Sparse as AD + -- | @AnyRF cb@ is a function that uses variables from the nonlinear -- program in a way supported by 'AnyRFCxt'. The @cb@ is -- usually 'Identity' @@ -35,7 +54,7 @@ -- | RealFloat gives most numerical operations, -- 'VectorSpace' is involved to allow using definitions from the -- package -type AnyRFCxt a = (RealFloat a, VectorSpace a, Scalar a ~ a) +type AnyRFCxt a = (VectorSpace a, RealFloat a, VectorSpace.Scalar a ~ a) -- *** helpers for defining instances liftOp0 :: (forall a. AnyRFCxt a => a) -> AnyRF Identity @@ -124,10 +143,34 @@ -- $orphans -- these belong somewhere between the @ad@ package and @vector-space@ -instance (Num a, AD.Mode f) => VectorSpace.AdditiveGroup (AD.AD f a) where - zeroV = AD.zero - (^+^) = (AD.<+>) - negateV = AD.negate1 -instance (Num a, AD.Mode f) => VectorSpace.VectorSpace (AD.AD f a) where - type Scalar (AD.AD f a) = AD.AD f a - (*^) = (AD.*!) + +instance (AD.Mode a) => VectorSpace.AdditiveGroup (AD.AD s a) where + zeroV = AD.zero + (^+^) = (+) + negateV = negate + +instance (AD.Mode a) => VectorSpace.VectorSpace (AD.AD s a) where + type Scalar (AD.AD s a) = AD.AD s a + (*^) = (*) + +instance (AD.Mode (AD.On a)) => VectorSpace.AdditiveGroup (AD.On a) where + zeroV = AD.zero + (^+^) = (+) + negateV = negate + +instance (AD.Mode (AD.On a)) => VectorSpace.VectorSpace (AD.On a) where + type Scalar (AD.On a) = AD.On a + (*^) = (*) + + + +instance (Num a, AD.Mode (AD.Reverse s a)) => VectorSpace.AdditiveGroup (AD.Reverse s a) where + zeroV = AD.zero + (^+^) = (+) + negateV = negate + +instance (Num a, AD.Mode (AD.Reverse s a)) => VectorSpace.VectorSpace (AD.Reverse s a) where + type Scalar (AD.Reverse a s) = AD.Reverse a s + (*^) = (*) + + diff --git a/Ipopt/NLP.hs b/Ipopt/NLP.hs --- a/Ipopt/NLP.hs +++ b/Ipopt/NLP.hs @@ -279,7 +279,17 @@ > (xMinus, xPlus) <- splitVar b x Using @max (x-b) 0@ instead of xPlus (ie. not telling the solver that @b@ is -a special point) seems to work just as well +a special point) seems to work just as well: additional special treatment +is needed. For example see chapter 11 of + +> Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes +> Lorenz T. Biegler +> SIAM 2010 + +which discusses several ways to reformulate the problem so that +an ordinary NLP solver will not have trouble with the fact that one of +the pair of constraints (@x+ = 0 | x- = 0@) is tight at an optimum. + -} splitVar :: (Monad m, Functor m) => Double -- ^ @b@ diff --git a/Ipopt/Raw.chs b/Ipopt/Raw.chs --- a/Ipopt/Raw.chs +++ b/Ipopt/Raw.chs @@ -38,6 +38,8 @@ IpProblem(..), + IntermediateCB, + ApplicationReturnStatus(..), -- * lower-level parts of the binding @@ -45,6 +47,7 @@ freeIpoptProblem, setIpoptProblemScaling, + setIntermediateCallback, -- ** marshalling functions @@ -68,8 +71,6 @@ ) where -import Ipopt.Sparsity - import C2HS import Control.Exception import Control.Monad @@ -84,6 +85,8 @@ import qualified Data.Vector.Storable as VS import qualified Data.Vector.Storable.Mutable as VM +import qualified Numeric.AD.Internal.Identity as AD + import Ipopt.AnyRF import Data.VectorSpace (VectorSpace,Scalar) @@ -100,6 +103,22 @@ type IpJacG = {# type Eval_Jac_G_CB #} type IpH = {# type Eval_H_CB #} +type IpIntermediateCB = {# type Intermediate_CB #} + +type IntermediateCB = CInt -- ^ alg_mod (0 regular, 1 is resto) + -> CInt -- ^ iter count + -> Double -- ^ obj value + -> Double -- ^ inf_pr + -> Double -- ^ inf_du + -> Double -- ^ mu + -> Double -- ^ d_norm + -> Double -- ^ regularization_size + -> Double -- ^ alpha_du + -> Double -- ^ alpha_pr + -> CInt -- ^ ls_trials + -> Ptr () -- ^ user_data (usually null) + -> IO IpBool + {#enum ApplicationReturnStatus as ^ {underscoreToCase} deriving (Show) #} {#enum AlgorithmMode as ^ {underscoreToCase} deriving (Show) #} @@ -110,9 +129,12 @@ type family UnFunPtr a type instance UnFunPtr (FunPtr a) = a -ipTrue = 1 :: IpBool -ipFalse = 0 :: IpBool +{#enum define IpBoolT { TRUE as IpTrue, FALSE as IpFalse } deriving (Eq,Ord,Show) #} +ipTrue = toEnum (fromEnum IpTrue) :: IpBool +ipFalse = toEnum (fromEnum IpFalse) :: IpBool + + -- | likely an unsafe method for getting a "Data.Vector.Storable.Mutable" out of a 'Ptr' ptrToVS n p = do fp <- newForeignPtr_ p @@ -124,7 +146,9 @@ foreign import ccall "wrapper" wrapIpJacG1 :: UnFunPtr IpJacG -> IO IpJacG foreign import ccall "wrapper" wrapIpH1 :: UnFunPtr IpH -> IO IpH +foreign import ccall "wrapper" wrapIpIntermediateCB :: IntermediateCB -> IO IpIntermediateCB + toB x = either (\ e@SomeException {} -> print e >> return ipFalse) (\ _ -> return ipTrue ) =<< try x @@ -208,7 +232,6 @@ foreign import ccall unsafe "&FreeIpoptProblem" freeIpoptProblem :: FunPtr (Ptr () -> IO ()) --- {#fun FreeIpoptProblem as ^ { ipp* `IpProblem' } -> `()' #} {#fun OpenIpoptOutputFile as ^ { ipp* `IpProblem', `String', `Int' } -> `Bool' #} @@ -219,6 +242,14 @@ vmUnsafeWith* `Vec' } -> `Bool' #} +{#fun SetIntermediateCallback as setIntermediateCallback1 + { ipp* `IpProblem', + id `IpIntermediateCB' } -> `Bool' #} + +setIntermediateCallback pp cb = do + cb' <- wrapIpIntermediateCB cb + setIntermediateCallback1 pp cb' + -- | lenses are in "Ipopt.PP" data IpOptSolved v = IpOptSolved { _ipOptSolved_status :: ApplicationReturnStatus, @@ -307,18 +338,17 @@ n == VM.length xU, m <- VM.length gL, m == VM.length gU = do - (eval_f, eval_grad_f, eval_g, eval_jac_g, eval_h) <- mkFs (False,True) n m f g + (eval_f, eval_grad_f, eval_g, eval_jac_g, eval_h) <- mkFs n m f g createIpoptProblem xL xU gL gU (n*m) (((n+1)*n) `div` 2) eval_f eval_g eval_grad_f eval_jac_g eval_h mkFs :: - (Bool,Bool) -- ^ grad, hessian are sparse? - -> Int -- ^ @n@ number of variables + Int -- ^ @n@ number of variables -> Int -- ^ @m@ number of constraints -> (forall a. AnyRFCxt a => V.Vector a -> a) -- ^ objective function @R^n -> R@ -> (forall a. AnyRFCxt a => V.Vector a -> V.Vector a) -- ^ constraint functions @R^n -> R^m@ -> IO (IpF, IpGradF, IpG, IpJacG, IpH) -mkFs (jacGs,hessGs) n m f g = do +mkFs n m f g = do ipF <- wrapIpF $ \x -> do x <- VG.convert `fmap` VS.unsafeFreeze x return $ f x @@ -331,13 +361,7 @@ x <- VG.convert `fmap` VS.unsafeFreeze x VS.unsafeThaw $ VG.convert (g x) - let ijJ = nanPropagateJ n g - ipJacG <- - if jacGs then wrapIpJacG (sparseIJ ijJ) $ \x y -> do - x <- VG.convert `fmap` VS.unsafeFreeze x - jac <- VS.unsafeThaw $ VG.convert $ jacobianSS ijJ g x - VM.copy y jac - else wrapIpJacG (denseIJ n m) $ \x y -> do + ipJacG <- wrapIpJacG (denseIJ n m) $ \x y -> do x <- VG.convert `fmap` VS.unsafeFreeze x jac <- VS.unsafeThaw $ VG.convert $ VG.concat $ VG.toList $ jacobian g x VM.copy y jac @@ -358,12 +382,6 @@ return (ipF, ipGradF, ipG, ipJacG, ipH) - -sparseIJ ijJ iRow jCol = do - let (iJ,jJ) = V.unzip ijJ - conv = VS.unsafeThaw . V.convert . V.map fromIntegral - VM.copy iRow =<< conv iJ - VM.copy jCol =<< conv jJ -- | indexes the same as http://www.coin-or.org/Ipopt/documentation/node40.html denseIJ n m iRow jCol = do diff --git a/Ipopt/Sparsity.hs b/Ipopt/Sparsity.hs deleted file mode 100644 --- a/Ipopt/Sparsity.hs +++ /dev/null @@ -1,158 +0,0 @@ -{- | - -Description: find out derivatives are unnecessary - -gradients and hessians tend to be sparse. This can be probed by -calculating these with NaNs as inputs. http://www.gpops2.com/ uses this -strategy, but perhaps there are relatively common cases (calls to -blas/lapack for example) that do not propagate NaNs correctly: perhaps -Ipopt.NLP provides some information about the structure of the problem, -provided that all variables used are lexically scoped? - -All functions provide indices that are affected (and should thus be included) - -functions ending in 1 set one variable at a time to NaN, and additionally -provide a hint as to which variable to change. - -functions not ending in 1 set all input variables to NaN - -it seems that adding a 1 achieves the same thing as doing one more derivative. -In other words, nanPropagateG1 tells the same thing as nanPropagateH. - --} -module Ipopt.Sparsity ( - -- $setup - -- * gradient sparsity - nanPropagate1, - nanPropagateG, - - -- * jacobian - nanPropagateJ, - -- * hessian sparsity - nanPropagateH, - nanPropagateG1, - - nanPropagateHF, - - -- * sparse derivatives (faked for now) - jacobianSS, - ) where - -import Control.Applicative -import Control.Monad -import qualified Data.Vector as V -import Data.Vector ((!)) -import Numeric.AD - -import qualified Data.IntMap as M -import qualified Data.IntSet as S -import Data.Monoid - -import Ipopt.AnyRF - -{- $setup ->>> let g x = x!0 + x!1 * x!1 * x!2 - --} - -{- | a nonzero gradient when inputs are NaN ==> no need to include -that row/column in the hessian, since it will be zero - ->>> nanPropagateG1 4 g -[(1,fromList [1,2]),(2,fromList [1])] - - - --} -nanPropagateG1 :: Int -- ^ size of input vector - -> (forall a. AnyRFCxt a => V.Vector a -> a) - -> [(Int, V.Vector Int)] -- ^ @(i,js)@ - -- shows that a NaN at index @i@ produces NaNs in gradient indexes - -- @js@... so the hessians only need to include - -- @i@ and @js@ -nanPropagateG1 nx f | x0:_ <- dropWhile (isNaN . f) (trialV nx [0,0.5,1]) - = [ (i,j) | i <- [0 .. nx-1], - let g = grad f (x0 V.// [(i,(0/0))]) - j = V.findIndices isNaN g, - not (V.null j) ] - -{- | - ->>> nanPropagateG 4 g -fromList [0,1,2] - --} -nanPropagateG :: Int -- ^ size of input vector - -> (forall a. AnyRFCxt a => V.Vector a -> a) - -> V.Vector Int -nanPropagateG nx f = V.findIndices (\x -> x /= 0) - (grad f (V.replicate nx (0/0))) - -{- | ->>> nanPropagate1 4 g -[0,1,2] - -variable 3 isn't even used. --} -nanPropagate1 :: Int -- ^ size of input vector - -> (forall a. AnyRFCxt a => V.Vector a -> a) - -> [Int] -- ^ inputs that don't become NaN -nanPropagate1 nx f | x0:_ <- dropWhile (isNaN . f) (trialV nx [0,0.5,1]) - = [ i | i <- [0 .. nx-1], - let g = f (x0 V.// [(i,(0/0))]), - isNaN g ] - - -{- | ->>> nanPropagateH 4 g -fromList [(1,1),(1,2),(2,1)] - --} -nanPropagateH :: Int - -> (forall a. AnyRFCxt a => V.Vector a -> a) - -> V.Vector (Int, Int) -- ^ (i,j) indexes -nanPropagateH nx f = nonzeroIxs $ hessian f (V.replicate nx (0/0::Double)) - -nanPropagateJ :: Int - -> (forall a. AnyRFCxt a => V.Vector a -> V.Vector a) - -> V.Vector (Int,Int) -nanPropagateJ nx f = nonzeroIxs $ jacobian f (V.replicate nx (0/0::Double)) - -nanPropagateHF :: Int - -> (forall a. AnyRFCxt a => V.Vector a -> V.Vector a) - -> V.Vector (V.Vector (Int,Int)) -nanPropagateHF nx f = - V.map nonzeroIxs $ hessianF f (V.replicate nx (0/0::Double)) - -{- -complement :: (Int,Int) -> V.Vector (Int,Int) -> V.Vector (Int,Int) -complement (mx,my) = M.fromList . M.toList --} - -collapse :: V.Vector (V.Vector (Int,Int)) -> V.Vector (Int,Int) -collapse = V.fromList . concatMap (\(a,bs) -> map (a,) (S.toList bs)) . M.toList . - M.mapWithKey (\k s -> fst (S.split (k+1) s)) . - M.fromListWith (<>) . map (\(a,b) -> (a, S.singleton b)) . concatMap V.toList . V.toList - -nonzeroIxs = V.concatMap id . V.imap (\i -> V.map (i,) . V.findIndices (/= 0)) - -trialV :: Int -> [Double] -> [V.Vector Double] -trialV nx cand = V.fromList <$> replicateM nx cand - - - -jacobianSS :: AnyRFCxt a => V.Vector (Int,Int) - -> (forall a. AnyRFCxt a => V.Vector a -> V.Vector a) - -> V.Vector a - -> V.Vector a -jacobianSS ijs f x = let v = jacobian f x - in V.map (\(i,j) -> v V.! i V.! j) ijs - -{- -hessianFSS :: RealFloat a => V.Vector (Int,Int) - -> (forall a. RealFloat a => V.Vector a -> V.Vector a) - -> V.Vector a - -> V.Vector a -hessianFSS = undefined --} - diff --git a/Nlopt/Raw.chs b/Nlopt/Raw.chs --- a/Nlopt/Raw.chs +++ b/Nlopt/Raw.chs @@ -114,7 +114,7 @@ instance Exception NloptResult checkEC :: CInt -> IO NloptResult -checkEC n | n < 0 = throwIO e +checkEC n | n `elem` [-1,-2,-3] = throwIO e | otherwise = return e where e = fromCInt n :: NloptResult diff --git a/ipopt-hs.cabal b/ipopt-hs.cabal --- a/ipopt-hs.cabal +++ b/ipopt-hs.cabal @@ -1,5 +1,5 @@ name: ipopt-hs -version: 0.4.0.1 +version: 0.4.2.0 synopsis: haskell binding to ipopt and nlopt including automatic differentiation description: a haskell binding to the nonlinear programming solver . Bindings to @@ -10,15 +10,14 @@ needs the c library ipopt installed. Also by default a binding to nlopt is included. Nlopt by default does not include shared libraries, which seems to be needed to run things from ghci (ie. you need to - @./configure --enable-shared@ when buildin nlopt). + @./configure --enable-shared@ when building nlopt). . A embedded language, similar to the one provided by glpk-hs, is defined in "Ipopt.NLP". The goal is to define problems at a level similar to other "algebraic modeling languages", but retain some of the safety and flexibility available in haskell. There is some overhead - but at least on the small 4-variable constrained optimization - problem. + but perhaps it is negligible for your uses. . Current limitations include: . @@ -29,7 +28,7 @@ . * sparseness of derivatives isn't used . - * no binding to SetIntermediateCallback + * no binding to sensitivity parts as-implemented in ipopt license: BSD3 license-file: LICENSE author: Adam Vogt @@ -54,7 +53,6 @@ exposed-modules: Ipopt, Ipopt.AnyRF, Ipopt.Raw, - Ipopt.Sparsity, Ipopt.NLP, Ipopt.Options, Ipopt.PP @@ -66,7 +64,7 @@ other-modules: C2HS build-depends: base < 5, - ad >=3.4, + ad >=4.2, ansi-wl-pprint >= 0.6.7, containers < 0.6, lens >= 3.10 && < 5, @@ -95,7 +93,7 @@ ipopt-hs, lens, mtl, ansi-wl-pprint, Rlang-QQ, vector-space, splines, ad, criterion, random-shuffle, - optimization >= 0.1.3, linear + linear hs-source-dirs: examples default-language: Haskell2010