ipopt-hs 0.4.0.1 → 0.4.2.0
raw patch · 7 files changed
+114/−201 lines, 7 filesdep −optimizationdep ~ad
Dependencies removed: optimization
Dependency ranges changed: ad
Files
- Ipopt.hs +2/−0
- Ipopt/AnyRF.hs +54/−11
- Ipopt/NLP.hs +11/−1
- Ipopt/Raw.chs +40/−22
- Ipopt/Sparsity.hs +0/−158
- Nlopt/Raw.chs +1/−1
- ipopt-hs.cabal +6/−8
Ipopt.hs view
@@ -18,6 +18,8 @@ ipopts, setIpoptProblemScaling, openIpoptOutputFile,+ setIntermediateCallback,+ IntermediateCB, -- * low-level bits still needed IpOptSolved(..),
Ipopt/AnyRF.hs view
@@ -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 -- <http://hackage.haskell.org/package/splines splines> 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+ (*^) = (*)++
Ipopt/NLP.hs view
@@ -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@
Ipopt/Raw.chs view
@@ -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
− Ipopt/Sparsity.hs
@@ -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--}-
Nlopt/Raw.chs view
@@ -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
ipopt-hs.cabal view
@@ -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 <http://projects.coin-or.org/Ipopt Ipopt>. 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 <http://code.haskell.org/~aavogt/ipopt-hs/examples/bench.html>- 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 <vogt.adam@gmail.com>@@ -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