diff --git a/C2HS.hs b/C2HS.hs --- a/C2HS.hs +++ b/C2HS.hs @@ -97,21 +97,6 @@ -- Everything else below is deprecated. -- These functions are not used by code generated by c2hs. -{-# DEPRECATED withBool "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED peekBool "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED withEnum "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED peekEnum "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED nothingIf "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED nothingIfNull "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED combineBitMasks "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED containsBitMask "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED extractBitMasks "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED cIntConv "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED cFloatConv "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED cFromBool "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED cToBool "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED cToEnum "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -{-# DEPRECATED cFromEnum "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} -- Passing Booleans by reference diff --git a/Ipopt.hs b/Ipopt.hs --- a/Ipopt.hs +++ b/Ipopt.hs @@ -1,4 +1,4 @@ --- | Description: exports most things you should need +-- | Description: exports most things you should need -- -- This module exports most things you should need. -- Also take a look at "Ipopt.NLP" and "Ipopt.Raw" and @examples/@ @@ -10,7 +10,7 @@ -- ** functions addG, addF, -- ** running the solver - ppSoln, + ppSoln, NLPT, nlpstate0, module Control.Monad.State, solveNLP', @@ -20,7 +20,7 @@ openIpoptOutputFile, setIntermediateCallback, IntermediateCB, - + -- * low-level bits still needed IpOptSolved(..), ApplicationReturnStatus(..), diff --git a/Ipopt/AnyRF.hs b/Ipopt/AnyRF.hs --- a/Ipopt/AnyRF.hs +++ b/Ipopt/AnyRF.hs @@ -174,3 +174,12 @@ (*^) = (*) + +instance Num a => VectorSpace.VectorSpace (AD.Sparse a) where + type Scalar (AD.Sparse a) = AD.Sparse a + (*^) = (*) + +instance (Num a, AD.Mode (AD.Sparse a)) => VectorSpace.AdditiveGroup (AD.Sparse a) where + zeroV = AD.zero + (^+^) = (+) + negateV = negate diff --git a/Ipopt/NLP.hs b/Ipopt/NLP.hs --- a/Ipopt/NLP.hs +++ b/Ipopt/NLP.hs @@ -102,7 +102,7 @@ makeLenses ''NLPFun makeLenses ''Ix makeLenses ''NLPState -makeIso ''IxMap +ixMap x = iso IxMap (\(IxMap a) -> a) x -- could makeLenses ''IxMap -- | should be a way to write an instance of At that'll make the normal -- at/ix work with the IxMap / Ix (as opposed to IntMap/Int) @@ -124,7 +124,7 @@ -- | calls 'createIpoptProblemAD' and 'ipoptSolve'. To be used at the -- end of a do-block. solveNLP' :: (VG.Vector v Double, MonadIO m) => - (IpProblem -> IO ()) -- ^ set ipopt options (using functions from "Ipopt.Raw") + (IpProblem -> IO ()) -- ^ set ipopt options (using functions from "Ipopt.Raw") or the 'ipopts' quasiquoter -> NLPT m (IpOptSolved v) solveNLP' setOpts = do (xl,xu) <- join $ uses (nlpfun . boundX) seqToVecs @@ -137,6 +137,23 @@ x0 <- uses initX V.convert r <- liftIO (ipoptSolve p =<< VS.thaw x0) return r + +-- | a slower version of 'solveNLP'' that uses 'createIpoptProblemADSparse' +solveNLP'sparse :: (VG.Vector v Double, MonadIO m) => + (IpProblem -> IO ()) -- ^ set ipopt options (using functions from "Ipopt.Raw") or the 'ipopts' quasiquoter + -> NLPT m (IpOptSolved v) +solveNLP'sparse setOpts = do + (xl,xu) <- join $ uses (nlpfun . boundX) seqToVecs + (gl,gu) <- join $ uses (nlpfun . boundG) seqToVecs + AnyRF fs <- use (nlpfun . funF) + AnyRF gs <- use (nlpfun . funG) + + p <- liftIO (createIpoptProblemADSparse xl xu gl gu (F.sum . fs) (V.fromList . toList . gs)) + liftIO (setOpts p) + x0 <- uses initX V.convert + r <- liftIO (ipoptSolve p =<< VS.thaw x0) + return r + -- | add a constraint addG :: Monad m diff --git a/Ipopt/PP.hs b/Ipopt/PP.hs --- a/Ipopt/PP.hs +++ b/Ipopt/PP.hs @@ -22,7 +22,13 @@ import Data.Foldable (toList,for_, foldMap) -- * lenses for IpOptSolved -makeLensesWith (lensRules & lensField .~ stripPrefix "_ipOptSolved_") ''IpOptSolved +status f s = (\b -> s { _ipOptSolved_status = b }) `fmap` f (_ipOptSolved_status s) +objective f s = (\b -> s { _ipOptSolved_objective = b }) `fmap` f (_ipOptSolved_objective s) +x f s = (\b -> s { _ipOptSolved_x = b }) `fmap` f (_ipOptSolved_x s) +g f s = (\b -> s { _ipOptSolved_g = b }) `fmap` f (_ipOptSolved_g s) +mult_g f s = (\b -> s { _ipOptSolved_mult_g = b }) `fmap` f (_ipOptSolved_mult_g s) +mult_x_L f s = (\b -> s { _ipOptSolved_mult_x_L = b }) `fmap` f (_ipOptSolved_mult_x_L s) +mult_x_U f s = (\b -> s { _ipOptSolved_mult_x_U = b }) `fmap` f (_ipOptSolved_mult_x_U s) -- * pretty printing ppSoln state0 problem = flip evalStateT state0 $ runWriterT $ do diff --git a/Ipopt/Raw.chs b/Ipopt/Raw.chs --- a/Ipopt/Raw.chs +++ b/Ipopt/Raw.chs @@ -11,6 +11,7 @@ module Ipopt.Raw ( -- * specifying problem createIpoptProblemAD, + createIpoptProblemADSparse, -- ** solve ipoptSolve, @@ -87,6 +88,10 @@ import qualified Numeric.AD.Internal.Identity as AD +import qualified Numeric.AD.Rank1.Sparse as Sparse +import qualified Numeric.AD.Internal.Sparse as Sparse +import Numeric.AD.Rank1.Sparse (Sparse) + import Ipopt.AnyRF import Data.VectorSpace (VectorSpace,Scalar) @@ -323,7 +328,7 @@ > ipoptSolve nlp x0 -Refer to the example @Test2.hs@ for details of setting up the vectors supplied. +Refer to @examples/HS71ad.hs@ for details of setting up the vectors supplied. -} createIpoptProblemAD :: Vec -- ^ @xL@ 'VM.Vector' of lower bounds for decision variables with length @n@ @@ -342,6 +347,31 @@ createIpoptProblem xL xU gL gU (n*m) (((n+1)*n) `div` 2) eval_f eval_g eval_grad_f eval_jac_g eval_h + +{- | this is 50% slower than 'createIpoptProblemAD' in one instance +). +But the benefit is that no RankN types are used (so it is possible to implement +more functions without having to modify 'AnyRFCxt') +-} +createIpoptProblemADSparse + + :: Vec -- ^ @xL@ 'VM.Vector' of lower bounds for decision variables with length @n@ + -> Vec -- ^ @xU@ 'VM.Vector' of upper bounds for decision variables + -> Vec -- ^ @gL@ 'VM.Vector' of lower bounds for constraint functions @g(x)@ with length @m@ + -> Vec -- ^ @gU@ 'VM.Vector' of upper bounds for constraint functions @g(x)@ + -> (V.Vector (Sparse CDouble) -> Sparse CDouble) -- ^ objective function @f : R^n -> R@ + -> (V.Vector (Sparse CDouble) -> V.Vector (Sparse CDouble)) -- ^ constraint functions @g : R^n -> R^m@ + -> IO IpProblem +createIpoptProblemADSparse xL xU gL gU f g + | n <- VM.length xL, + 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) <- mkFsSparse 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 :: Int -- ^ @n@ number of variables -> Int -- ^ @m@ number of constraints @@ -382,6 +412,49 @@ return (ipF, ipGradF, ipG, ipJacG, ipH) +sparsePrimal :: Num a => Sparse a -> a +sparsePrimal (Sparse.Sparse a _) = a +sparsePrimal Sparse.Zero = 0 + +mkFsSparse :: + Int -- ^ @n@ number of variables + -> Int -- ^ @m@ number of constraints + -> (V.Vector (Sparse CDouble) -> Sparse CDouble) -- ^ objective function @R^n -> R@ + -> (V.Vector (Sparse CDouble) -> V.Vector (Sparse CDouble)) -- ^ constraint functions @R^n -> R^m@ + -> IO (IpF, IpGradF, IpG, IpJacG, IpH) +mkFsSparse n m f g = do + ipF <- wrapIpF $ \x -> do + x <- VG.convert `fmap` VS.unsafeFreeze x + return $ sparsePrimal $ f (Sparse.vars x) + + ipGradF <- wrapIpGradF $ \x -> do + x <- VG.convert `fmap` VS.unsafeFreeze x + VS.unsafeThaw $ VG.convert (Sparse.grad f x) + + ipG <- wrapIpG $ \x -> do + x <- VG.convert `fmap` VS.unsafeFreeze x + VS.unsafeThaw $ VG.convert $ V.map sparsePrimal $ g $ Sparse.vars x + + ipJacG <- wrapIpJacG (denseIJ n m) $ \x y -> do + x <- VG.convert `fmap` VS.unsafeFreeze x + jac <- VS.unsafeThaw $ VG.convert $ VG.concat $ VG.toList $ Sparse.jacobian g x + VM.copy y jac + + ipH <- wrapIpH (denseIJh n m) + ( \ obj_factor lambda x values -> do + + x <- VG.convert `fmap` VS.unsafeFreeze x + lambda <- VG.convert `fmap` VS.unsafeFreeze lambda + + let tri = VG.concat . VG.toList . V.imap (\n -> V.take (n+1)) + obj = V.map (*obj_factor) $ tri $ Sparse.hessian f x + gj = V.zipWith (\l v -> V.map (l*) v) lambda (V.map tri (Sparse.hessianF g x)) + lagrangian = V.foldl (V.zipWith (+)) obj gj + + VM.copy values =<< VS.unsafeThaw (VG.convert lagrangian) + ) + + return (ipF, ipGradF, ipG, ipJacG, ipH) -- | indexes the same as http://www.coin-or.org/Ipopt/documentation/node40.html denseIJ n m iRow jCol = do 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.2.0 +version: 0.5.0.0 synopsis: haskell binding to ipopt and nlopt including automatic differentiation description: a haskell binding to the nonlinear programming solver . Bindings to @@ -24,7 +24,9 @@ * copying in every iteration happens between between "Data.Vector.Storable" and "Data.Vector" might be avoidable somehow. Currently it is done because AD needs a Traversable - structure, but Storable vectors are not traversable. + structure, but Storable vectors are not traversable. Note that + there is a "Numeric.AD.Mode.Forward.Double" but as yet none + for Vectors of Double . * sparseness of derivatives isn't used . @@ -67,7 +69,7 @@ ad >=4.2, ansi-wl-pprint >= 0.6.7, containers < 0.6, - lens >= 3.10 && < 5, + lens >= 3.7 && < 5, mtl == 2.*, template-haskell, uu-parsinglib >= 2.8, @@ -85,15 +87,18 @@ TemplateHaskell pkgconfig-depends: ipopt build-tools: c2hs + -- for ghci-7.8 + cc-options: -fPIC executable ipopt-hs_Tests main-is: AllTests.hs - build-depends: base <= 5, vector ==0.10.*, - ipopt-hs, lens, mtl, ansi-wl-pprint, - Rlang-QQ, vector-space, splines, ad, - criterion, random-shuffle, - linear + if (flag(build_examples)) + build-depends: base <= 5, vector ==0.10.*, + ipopt-hs, lens, mtl, ansi-wl-pprint, + Rlang-QQ, vector-space, splines, ad, + criterion, random-shuffle, + linear hs-source-dirs: examples default-language: Haskell2010