ipopt-hs 0.0.0.0 → 0.2.0.0
raw patch · 11 files changed
+1114/−540 lines, 11 filesdep +containersdep +lensdep +mtlnew-component:exe:ipopt-hs_Test3
Dependencies added: containers, lens, mtl
Files
- Ipopt.chs +0/−359
- Ipopt.hs +36/−0
- Ipopt/NLP.hs +315/−0
- Ipopt/Raw.chs +357/−0
- Ipopt/Sparsity.hs +155/−0
- Test1.hs +0/−123
- Test2.hs +0/−52
- examples/Test1.hs +123/−0
- examples/Test2.hs +52/−0
- examples/Test3.hs +20/−0
- ipopt-hs.cabal +56/−6
− Ipopt.chs
@@ -1,359 +0,0 @@-{-# LANGUAGE ForeignFunctionInterface, PatternGuards, RankNTypes, TypeFamilies #-}-{- |--Copyright: (C) 2013 Adam Vogt-Maintainer: Adam Vogt <vogt.adam@gmail.com>-Stability: unstable--Binding to ipopt <http://projects.coin-or.org/Ipopt>. Uses "Numeric.AD" to compute-derivatives required by ipopt. Current limitations include:--* derivatives are computed and stored without taking advantage of sparsity--* copying is done in converting between "Data.Vector.Storable" and "Data.Vector" might be unnecessary. Currently it is done because AD needs a Traversable structure, but Storable vectors are not traversable.--* probably doesn't work if @coin\/IpStdCInterface.h@ has Number =/= 'CDouble'--* no binding to SetIntermediateCallback--* garbage collection of 'IpProblem' won't free C-side resources--* specifying problems might be made easier by following (extending) the approach taken by glpk-hs--Refer to @Test1.hs@ for an example where the derivatives are computed by hand,-and @Test2.hs@ for the use of 'createIpoptProblemAD'.---}-module Ipopt (- -- * specifying problem- createIpoptProblemAD,-- -- ** solve- ipoptSolve,- IpOptSolved(..),-- -- ** solver options- addIpoptNumOption,- addIpoptStrOption,- addIpoptIntOption,- openIpoptOutputFile,-- -- * types- Vec,-- IpNumber(..),- IpIndex(..),- IpInt(..),- IpBool(..),-- IpF(..),- IpGradF(..),- IpG(..),- IpJacG(..),- IpH(..),-- IpProblem(..),-- ApplicationReturnStatus(..),-- -- * lower-level parts of the binding- createIpoptProblem,-- freeIpoptProblem,- setIpoptProblemScaling,--- -- ** marshalling functions- wrapIpF,- wrapIpGradF,- wrapIpG,- wrapIpJacG,- wrapIpH,-- wrapIpF1,- wrapIpGradF1,- wrapIpG1,- wrapIpJacG1,- wrapIpH1,-- wrapIpF2,- wrapIpGradF2,- wrapIpG2,- wrapIpJacG2,- wrapIpH2,-- ) where--import C2HS-import Control.Exception-import Control.Monad-import Data.IORef-import Foreign.C-import Foreign.ForeignPtr-import Foreign.Ptr-import Foreign.Storable-import Numeric.AD-import qualified Data.Vector as V-import qualified Data.Vector.Generic as VG-import qualified Data.Vector.Storable as VS-import qualified Data.Vector.Storable.Mutable as VM----#include "coin/IpStdCInterface.h"--type IpNumber = {# type Number #}-type IpIndex = {# type Index #}-type IpInt = {# type Int #}-type IpBool = {# type Bool #}--type IpF = {# type Eval_F_CB #}-type IpGradF = {# type Eval_Grad_F_CB #}-type IpG = {# type Eval_G_CB #}-type IpJacG = {# type Eval_Jac_G_CB #}-type IpH = {# type Eval_H_CB #}--{#enum ApplicationReturnStatus as ^ {underscoreToCase} deriving (Show) #}-{#enum AlgorithmMode as ^ {underscoreToCase} deriving (Show) #}---newtype IpProblem = IpProblem { unIpProblem :: Ptr ()}--type family UnFunPtr a-type instance UnFunPtr (FunPtr a) = a--ipTrue = 1 :: IpBool-ipFalse = 0 :: IpBool---- | likely an unsafe method for getting a "Data.Vector.Storable.Mutable" out of a 'Ptr'-ptrToVS n p = do- fp <- newForeignPtr_ p- return (VM.unsafeFromForeignPtr0 fp (fromIntegral n))--foreign import ccall "wrapper" wrapIpF1 :: UnFunPtr IpF -> IO IpF-foreign import ccall "wrapper" wrapIpG1 :: UnFunPtr IpG -> IO IpG-foreign import ccall "wrapper" wrapIpGradF1 :: UnFunPtr IpGradF -> IO IpGradF-foreign import ccall "wrapper" wrapIpJacG1 :: UnFunPtr IpJacG -> IO IpJacG-foreign import ccall "wrapper" wrapIpH1 :: UnFunPtr IpH -> IO IpH---toB x = either (\ e@SomeException {} -> print e >> return ipFalse)- (\ _ -> return ipTrue ) =<< try x--wrapIpF2' fun n xin new_x obj_val _userData = do- toB $ poke obj_val =<< fun =<< ptrToVS n xin--wrapIpF2 fun n xin new_x obj_val _userData = do- toB $ poke obj_val =<< fun =<< ptrToVS n xin--wrapIpG2 fun n xin new_x m gout _userData = do- toB $ join $ liftM2 VM.copy (ptrToVS m gout) (fun =<< ptrToVS n xin)--wrapIpGradF2 fun n x new_x grad_f _userData = do- toB $ join $ liftM2 VM.copy (ptrToVS n grad_f) (fun =<< ptrToVS n x)--wrapIpJacG2 fun1 fun2 n x new_x m nj iRow jCol jacs _userData- | jacs == nullPtr = do- toB $ join $ liftM2 fun1 (ptrToVS nj iRow) (ptrToVS nj jCol)- | otherwise = do- toB $ join $ liftM2 fun2 (ptrToVS n x) (ptrToVS nj jacs)--wrapIpH2 funSparsity funEval n x new_x obj_factor m lambda new_lambda nHess iRow jCol values _userData- | iRow == nullPtr = do- toB $ join $ liftM3 (funEval obj_factor)- (ptrToVS m lambda)- (ptrToVS n x)- (ptrToVS nHess values)- | otherwise = do- toB $ join $ liftM2 funSparsity- (ptrToVS nHess iRow)- (ptrToVS nHess jCol)--wrapIpF f = wrapIpF1 (wrapIpF2 f)-wrapIpG f = wrapIpG1 (wrapIpG2 f)-wrapIpGradF f = wrapIpGradF1 (wrapIpGradF2 f)-wrapIpJacG f1 f2 = wrapIpJacG1 (wrapIpJacG2 f1 f2)-wrapIpH fSparsity fEval = wrapIpH1 (wrapIpH2 fSparsity fEval)---vmUnsafeWith = VM.unsafeWith---- | Vector of numbers-type Vec = VM.IOVector IpNumber--createIpoptProblem :: Vec -> Vec -> Vec -> Vec- -> Int -> Int -> IpF -> IpG -> IpGradF -> IpJacG -> IpH -> IO IpProblem-createIpoptProblem xL xU gL gU nJac nHess f g gradF jacG hess- | lx <- VM.length xL,- lx == VM.length xU,- lg <- VM.length gL,- lg == VM.length gU = createIpoptProblem3 lx xL xU lg gL gU nJac nHess 0 f g gradF jacG hess- | otherwise = error "dimensions wrong!"--{#fun CreateIpoptProblem as createIpoptProblem3- { `Int', vmUnsafeWith* `Vec', vmUnsafeWith* `Vec',- `Int', vmUnsafeWith* `Vec', vmUnsafeWith* `Vec',- `Int', `Int', `Int', id `IpF', id `IpG', id `IpGradF',- id `IpJacG', id `IpH' } -> `IpProblem' IpProblem #}-_ = {#fun AddIpoptNumOption as ^- { unIpProblem `IpProblem', `String', `Double' } -> `Bool' #}--_ = {#fun AddIpoptStrOption as ^- { unIpProblem `IpProblem', `String', `String' } -> `Bool' #}--_ = {#fun AddIpoptIntOption as ^- { unIpProblem `IpProblem', `String', `Int' } -> `Bool' #}--_ = {#fun FreeIpoptProblem as ^- { unIpProblem `IpProblem' } -> `()' #}--_ = {#fun OpenIpoptOutputFile as ^- { unIpProblem `IpProblem', `String', `Int' } -> `Bool' #}--_ = {#fun SetIpoptProblemScaling as ^- { unIpProblem `IpProblem',- `Double',- vmUnsafeWith* `Vec',- vmUnsafeWith* `Vec'- } -> `Bool' #}---data IpOptSolved = IpOptSolved- { ipOptSolved_status :: ApplicationReturnStatus,- ipOptSolved_objective :: Double,- ipOptSolved_x,- ipOptSolved_g,- ipOptSolved_mult_g,- ipOptSolved_mult_x_L,- ipOptSolved_mult_x_U :: Vec }--ipoptSolve :: IpProblem- -> Vec -- ^ starting point @x@. Note that the value is overwritte with the final @x@.- -> IO IpOptSolved-ipoptSolve problem x = do- g <- VM.new (VM.length x)- mult_g <- VM.new (VM.length x)- mult_x_L <- VM.new (VM.length x)- mult_x_U <- VM.new (VM.length x)-- out <- ipoptSolve2- problem- x- g- mult_g- mult_x_L- mult_x_U- nullPtr- return $ IpOptSolved- (fst out)- (snd out)- x- g- mult_g- mult_x_L- mult_x_U--_ = {#fun IpoptSolve as ipoptSolve2- { unIpProblem `IpProblem',- vmUnsafeWith* `Vec',- vmUnsafeWith* `Vec',- alloca- `Double' peekFloatConv*,- vmUnsafeWith* `Vec',- vmUnsafeWith* `Vec',- vmUnsafeWith* `Vec',- id `Ptr ()' } -> `ApplicationReturnStatus' cToEnum #}---{- | Set-up an 'IpProblem' to be solved later. Only objective function (@f@)-and constraint functions (@g@) need to be specified. Derivatives needed by ipopt-are computed by "Numeric.AD".--To solve the optimization problem:--> min f(x)-> such that-> xL <= x <= xU-> gL <= g(x) <= gU--First create an opaque 'IpProblem' object (nlp):--> nlp <- createIpOptProblemAD xL xU gL gU f g--Then pass it off to 'ipoptSolve'.--> ipoptSolve nlp x0--Refer to the example @Test2.hs@ for details of setting up the vectors supplied.--}-createIpoptProblemAD- :: 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)@- -> (forall a. Num a => V.Vector a -> a) -- ^ objective function @f : R^n -> R@- -> (forall a. Num a => V.Vector a -> V.Vector a) -- ^ constraint functions @g : R^n -> R^m@- -> IO IpProblem-createIpoptProblemAD 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) <- 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 :: Int -- ^ @n@ number of variables- -> Int -- ^ @m@ number of constraints- -> (forall a. Num a => V.Vector a -> a) -- ^ objective function @R^n -> R@- -> (forall a. Num a => V.Vector a -> V.Vector a) -- ^ constraint functions @R^n -> R^m@- -> IO (IpF, IpGradF, IpG, IpJacG, IpH)-mkFs n m f g = do- ipF <- wrapIpF $ \x -> do- x <- VG.convert `fmap` VS.unsafeFreeze x- return $ f x-- ipGradF <- wrapIpGradF $ \x -> do- x <- VG.convert `fmap` VS.unsafeFreeze x- VS.unsafeThaw $ VG.convert (grad f x)-- ipG <- wrapIpG $ \x -> do- x <- VG.convert `fmap` VS.unsafeFreeze x- VS.unsafeThaw $ VG.convert (g 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 $ 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 $ hessian f x- gj = V.zipWith (\l v -> V.map (l*) v) lambda (V.map tri (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- VM.copy iRow =<< VS.unsafeThaw (VS.generate (n*m) (\x -> fromIntegral $ x `div` n))- VM.copy jCol =<< VS.unsafeThaw (VS.generate (n*m) (\x -> fromIntegral $ x `mod` n))---- | indexes the same as http://www.coin-or.org/Ipopt/documentation/node41.html-denseIJh n m iRow jCol = do- i <- newIORef 0- forM_ [0 .. fromIntegral n-1] $ \ row ->- forM_ [ 0 .. row ] $ \col -> do- ii <- readIORef i- VM.write iRow ii row- VM.write jCol ii col- writeIORef i (ii+1)-
+ Ipopt.hs view
@@ -0,0 +1,36 @@+-- | 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/@+module Ipopt (+ -- * high-level+ NLPT, nlpstate0,+ module Control.Monad.State,+ solveNLP',+ -- ** variables+ var', var, varFresh,+ AnyRF(..),+ -- ** functions+ addG, addF,+ -- * low-level bits still needed+ IpOptSolved(..),+ ApplicationReturnStatus(..),++ -- ** solver options+ -- $solverOptRef+ addIpoptNumOption,+ addIpoptStrOption, + addIpoptIntOption,+ setIpoptProblemScaling,+ openIpoptOutputFile,++ -- *** types+ Vec, IpNumber,+) where++import Ipopt.NLP+import Ipopt.Raw+import Control.Monad.State++-- $solverOptRef+-- see <http://www.coin-or.org/Ipopt/documentation/node39.html ipopt's options reference>
+ Ipopt/NLP.hs view
@@ -0,0 +1,315 @@+{-# LANGUAGE FlexibleInstances, GADTs, RankNTypes, TemplateHaskell #-}+-- | Description : a EDSL for describing nonlinear programs+--+-- see usage in @examples/Test3.hs@+module Ipopt.NLP where++import Control.Applicative+import Control.Lens+import Control.Monad+import Control.Monad.Identity+import Control.Monad.State+import Data.Foldable (toList)+import Data.IntMap (IntMap)+import Data.List+import Data.Monoid+import Data.Monoid (First)+import Data.Sequence (Seq)+import Data.Vector (Vector)+import Foreign.C.Types (CDouble(..))+import qualified Data.Foldable as F+import qualified Data.IntMap as IM+import qualified Data.Map as M+import qualified Data.Sequence as Seq+import qualified Data.Set as S+import qualified Data.Vector as V+import qualified Data.Vector.Storable as VS++import Ipopt.Raw++-- * state+data NLPFun = NLPFun+ { _funF, _funG :: AnyRF Seq,+ _boundX, _boundG :: Seq (Double,Double) }+instance Show NLPFun where+ show (NLPFun f g a b) = "NLPFun <f> <g> {" ++ show a ++ "}{" ++ show b ++ "}"++data NLPState = NLPState+ { -- | current maximum index+ _nMax :: Ix,+ -- | what namespace (see 'inEnv')+ _currentEnv :: [String],+ _variables :: M.Map String Ix,+ -- | human-readable descriptions for the constraint, objective and+ -- variables+ _constraintLabels, _objLabels, _varLabels :: IntMap String,+ _varEnv :: IntMap (S.Set [String]),+ _constraintEnv, _objEnv :: IntMap [String],+ _nlpfun :: NLPFun,+ _defaultBounds :: (Double,Double),+ -- | inital state variable for the solver+ _initX :: Vector Double }+ deriving (Show)++-- | solver deals with arrays. This type is for indexes into the array+-- for the current variables that the solver is trying to find.+newtype Ix = Ix { _varIx :: Int } deriving Show++-- | the initial state to use when you actually have to get to IO+-- with the solution+nlpstate0 = NLPState (Ix (-1)) mempty mempty+ mempty mempty mempty -- labels+ mempty mempty mempty -- env+ (mempty :: NLPFun)+ (-1/0, 1/0)+ V.empty++type NLPT = StateT NLPState+type NLP = NLPT IO+ +-- ** representing functions++{- | this wrapper holds functions that can be used+for the objective (`f`) or for constraints (`g`). Many functions+in the instances provided are partial: this seems to be unavoidable+because the input variables haven't been decided yet, so you should+not be allowed to use 'compare' on these. But for now just use the+standard Prelude classes, and unimplementable functions (which+would not produce an 'AnyRF') are calls to 'error'++generate these using 'var', or perhaps by directly using the constructor:+@AnyRF $ Identity . V.sum@, would for example give the sum of all variables.+-}+data AnyRF cb = AnyRF (forall a. RealFloat a => Vector a -> cb a)++-- *** helpers for defining instances+liftOp0 :: (forall a. RealFloat a => a) -> AnyRF Identity+liftOp0 op = AnyRF $ \x -> Identity op++liftOp1 :: (forall a. RealFloat a => a -> a) -> AnyRF Identity -> AnyRF Identity+liftOp1 op (AnyRF a) = AnyRF $ \x -> Identity (op (runIdentity (a x)))++liftOp2 :: (forall a. RealFloat a => a -> a -> a) -> AnyRF Identity -> AnyRF Identity -> AnyRF Identity+liftOp2 op (AnyRF a) (AnyRF b) = AnyRF $ \x -> Identity (runIdentity (a x) `op` runIdentity (b x))++instance Num (AnyRF Identity) where+ (+) = liftOp2 (+)+ (*) = liftOp2 (*)+ (-) = liftOp2 (-)+ abs = liftOp1 abs+ signum = liftOp1 signum+ fromInteger n = liftOp0 (fromInteger n)++instance Fractional (AnyRF Identity) where+ (/) = liftOp2 (/)+ recip = liftOp1 recip+ fromRational n = liftOp0 (fromRational n)++instance Floating (AnyRF Identity) where+ pi = liftOp0 pi+ exp = liftOp1 exp+ sqrt = liftOp1 sqrt+ log = liftOp1 log+ sin = liftOp1 sin+ tan = liftOp1 tan+ cos = liftOp1 cos+ asin = liftOp1 asin+ atan = liftOp1 atan+ acos = liftOp1 acos+ sinh = liftOp1 sinh+ tanh = liftOp1 tanh+ cosh = liftOp1 cosh+ asinh = liftOp1 asinh+ atanh = liftOp1 atanh+ acosh = liftOp1 acosh+ (**) = liftOp2 (**)+ logBase = liftOp2 logBase++instance Real (AnyRF Identity) where+ toRational _ = error "Real AnyRF Identity"++instance Ord (AnyRF Identity)+instance Eq (AnyRF Identity)+instance RealFrac (AnyRF Identity) where+ properFraction = error "properFraction AnyRF"+instance RealFloat (AnyRF Identity) where+ isInfinite = error "isInfinite AnyRF"+ isNaN = error "isNaN AnyRF"+ decodeFloat = error "decodeFloat AnyRF"+ floatRange = error "floatRange AnyRF"+ isNegativeZero = error "isNegativeZero AnyRF"+ isIEEE = error "isIEEE AnyRF"+ isDenormalized = error "isDenormalized AnyRF"+ floatDigits _ = floatDigits (error "RealFrac AnyRF Identity floatDigits" :: Double)+ floatRadix _ = floatRadix (error "RealFrac AnyRF Identity floatRadix" :: Double)+ atan2 = liftOp2 atan2+ significand = liftOp1 significand+ scaleFloat n = liftOp1 (scaleFloat n)+ encodeFloat a b = liftOp0 (encodeFloat a b)++instance Monoid (AnyRF Seq) where+ AnyRF f `mappend` AnyRF g = AnyRF (f `mappend` g)+ mempty = AnyRF mempty++instance Monoid NLPFun where+ NLPFun f g bx bg `mappend` NLPFun f' g' bx' bg' = NLPFun (f <> f') (g <> g') (bx <> bx') (bg <> bg')+ mempty = NLPFun mempty mempty mempty mempty++-- ** low level lenses to NLPState+makeLenses ''NLPFun+makeLenses ''Ix+makeLenses ''NLPState++-- | @to@ is one of 'varEnv', 'constraintEnv', 'objEnv'+copyEnv to n = do+ ev <- use currentEnv+ cloneLens to %= IM.insert n ev++-- | @to@ should be 'constraintLabels', 'objLabels', 'varLabels'+addDesc to Nothing n = return ()+addDesc to (Just x) n = do+ to %= IM.insert n x++-- * high-level functions++-- | calls 'createIpoptProblemAD' and 'ipoptSolve'. To be used at the+-- end of a do-block.+solveNLP' :: MonadIO m =>+ (IpProblem -> IO ()) -- ^ set ipopt options (using functions from "Ipopt.Raw")+ -> NLPT m IpOptSolved+solveNLP' 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 (createIpoptProblemAD xl xu gl gu (F.sum . fs) (V.fromList . toList . gs))+ liftIO (setOpts p)+ x0 <- uses initX (V.convert . V.map CDouble)+ r <- liftIO (ipoptSolve p =<< VS.thaw x0)+ liftIO $ freeIpoptProblem p+ return r++-- | add a constraint+addG :: Monad m+ => Maybe String -- ^ optional description+ -> (Double,Double) -- ^ bounds @(gl,gu)@ for the single inequality @gl_i <= g_i(x) <= gu_i@+ -> AnyRF Identity -- ^ @g_i(x)@+ -> NLPT m ()+addG d b (AnyRF f) = do+ nlpfun . boundG %= (Seq.|> b)+ nlpfun . funG %= \(AnyRF fs) -> AnyRF $ \x -> fs x Seq.|> runIdentity (f x)+ n <- use (nlpfun . boundG . to Seq.length)+ copyEnv constraintEnv n+ addDesc constraintLabels d n++{- | add a piece of the objective function, which is added in the form+`f_1 + f_2 + ...`, to make it easier to understand (at some point)+which components are responsible for the majority of the cost, and+which are irrelevant.+-}+addF :: Monad m+ => Maybe String -- ^ description+ -> AnyRF Identity -- ^ `f_i(x)`+ -> NLPT m ()+addF d (AnyRF f) = do+ nlpfun . funF %= \(AnyRF fs) -> AnyRF $ \x -> fs x Seq.|> runIdentity (f x)+ n <- use (objEnv . to ((+1) . IM.size))+ copyEnv objEnv n+ addDesc objLabels d n++-- | add a variable, or get a reference to the the same variable if it has+-- already been used+var' :: (Monad m, Functor m)+ => Maybe (Double,Double) -- ^ bounds @(xl,xu)@ to request that @xl <= x <= xu@.+ -- if Nothing, you get whatever is in 'defaultBounds'+ -> String -- ^ variable name (namespace from the 'pushEnv' / 'popEnv' can+ -- make an @"x"@ you request here different from one you+ -- previously requested+ -> NLPT m (AnyRF Identity, Ix) -- ^ the value, and index (into the raw+ -- vector of variables that the solver+ -- sees)+var' bs s = do+ ev <- use currentEnv+ m <- use variables+ let s' = intercalate "." (reverse (s:ev))+ n <- case M.lookup s' m of+ Nothing -> do+ nMax %= over varIx (+1)+ db <- use defaultBounds+ nlpfun . boundX %= (Seq.|> db)+ n' <- use nMax+ variables %= M.insert s' n'+ return n'+ Just n -> return n+ varEnv %= IM.insert (view varIx n) (S.singleton ev)+ F.traverse_ (narrowBounds n) bs+ return (AnyRF $ \x -> Identity $ x V.! view varIx n, n)++-- | 'var'' without the usually unnecessary 'Ix'+var bs s = fmap fst (var' bs s)++{- | 'var', except this causes the solver to get a new variable,+so that you can use:++> [a,b,c,d,e] <- replicateM 5 (varFresh (Just (0, 10)) "x")++and the different letters can take different values (between 0 and 10)+in the optimal solution (depending on what you do with @a@ and similar+in the objective function).+-}+varFresh bs s = do+ n <- uses variables ((+1) . M.size)+ var bs (s ++ show n)++-- *** namespace+{- $namespace++When you build up an optimization problem, it may be composed of pieces.+Functions in this section help to ease the pain of making unique variables.+To illustrate:++> m <- inEnv "A" (var b "x")+> n <- var b "A.x" ++@m@ and @n@ above should refer to the same variable. In some sense this+is \"better\" that using 'varFresh' all the time, since perhaps you would+like to add dependencies between components (say the size of a header pipe,+refridgeration unit, foundation etc. has to satisfy sizes of individual+components).++-}++-- | combination of 'pushEnv' and 'popEnv'+inEnv :: Monad m => String -> NLPT m a -> NLPT m a+inEnv s action = do+ pushEnv s+ r <- action+ popEnv+ return r++pushEnv :: Monad m => String -> NLPT m ()+pushEnv s = currentEnv %= (s:)++popEnv :: Monad m => NLPT m String+popEnv = do+ n : ns <- use currentEnv+ currentEnv .= ns+ return n++-- *** bounds++-- | override bounds. Should be unnecessary given 'var' takes bounds.+setBounds :: Monad m => Ix -> (Double,Double) -> NLPT m ()+setBounds (Ix i) bs = nlpfun . boundX %= Seq.update i bs++-- | shrink the interval in which that variable is allowed.+narrowBounds :: Monad m => Ix -> (Double,Double) -> NLPT m ()+narrowBounds (Ix i) (a,b) = nlpfun . boundX . ix i %= \(a',b') -> (max a a', min b b')++-- * internal+seqToVecs :: MonadIO m => Seq (Double,Double) -> m (Vec,Vec)+seqToVecs x = let (a,b) = unzip (toList x) in liftIO $ do+ a' <- VS.thaw (VS.map CDouble (VS.fromList a))+ b' <- VS.thaw (VS.map CDouble (VS.fromList b))+ return (a',b')
+ Ipopt/Raw.chs view
@@ -0,0 +1,357 @@+{-# LANGUAGE ForeignFunctionInterface, PatternGuards, RankNTypes, TypeFamilies #-}+{- |++Copyright: (C) 2013 Adam Vogt+Maintainer: Adam Vogt <vogt.adam@gmail.com>+Stability: unstable+Description: lowest-level parts of the binding+++-}+module Ipopt.Raw (+ -- * specifying problem+ createIpoptProblemAD,++ -- ** solve+ ipoptSolve,+ IpOptSolved(..),++ -- ** solver options+ addIpoptNumOption,+ addIpoptStrOption,+ addIpoptIntOption,+ openIpoptOutputFile,++ -- * types+ Vec,++ IpNumber(..),+ IpIndex(..),+ IpInt(..),+ IpBool(..),++ IpF(..),+ IpGradF(..),+ IpG(..),+ IpJacG(..),+ IpH(..),++ IpProblem(..),++ ApplicationReturnStatus(..),++ -- * lower-level parts of the binding+ createIpoptProblem,++ freeIpoptProblem,+ setIpoptProblemScaling,+++ -- ** marshalling functions+ wrapIpF,+ wrapIpGradF,+ wrapIpG,+ wrapIpJacG,+ wrapIpH,++ wrapIpF1,+ wrapIpGradF1,+ wrapIpG1,+ wrapIpJacG1,+ wrapIpH1,++ wrapIpF2,+ wrapIpGradF2,+ wrapIpG2,+ wrapIpJacG2,+ wrapIpH2,++ ) where++import Ipopt.Sparsity++import C2HS+import Control.Exception+import Control.Monad+import Data.IORef+import Foreign.C+import Foreign.ForeignPtr+import Foreign.Ptr+import Foreign.Storable+import Numeric.AD+import qualified Data.Vector as V+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VM++#include "IpStdCInterface.h"++type IpNumber = {# type Number #}+type IpIndex = {# type Index #}+type IpInt = {# type Int #}+type IpBool = {# type Bool #}++type IpF = {# type Eval_F_CB #}+type IpGradF = {# type Eval_Grad_F_CB #}+type IpG = {# type Eval_G_CB #}+type IpJacG = {# type Eval_Jac_G_CB #}+type IpH = {# type Eval_H_CB #}++{#enum ApplicationReturnStatus as ^ {underscoreToCase} deriving (Show) #}+{#enum AlgorithmMode as ^ {underscoreToCase} deriving (Show) #}+++newtype IpProblem = IpProblem { unIpProblem :: Ptr ()}++type family UnFunPtr a+type instance UnFunPtr (FunPtr a) = a++ipTrue = 1 :: IpBool+ipFalse = 0 :: IpBool++-- | likely an unsafe method for getting a "Data.Vector.Storable.Mutable" out of a 'Ptr'+ptrToVS n p = do+ fp <- newForeignPtr_ p+ return (VM.unsafeFromForeignPtr0 fp (fromIntegral n))++foreign import ccall "wrapper" wrapIpF1 :: UnFunPtr IpF -> IO IpF+foreign import ccall "wrapper" wrapIpG1 :: UnFunPtr IpG -> IO IpG+foreign import ccall "wrapper" wrapIpGradF1 :: UnFunPtr IpGradF -> IO IpGradF+foreign import ccall "wrapper" wrapIpJacG1 :: UnFunPtr IpJacG -> IO IpJacG+foreign import ccall "wrapper" wrapIpH1 :: UnFunPtr IpH -> IO IpH+++toB x = either (\ e@SomeException {} -> print e >> return ipFalse)+ (\ _ -> return ipTrue ) =<< try x++wrapIpF2' fun n xin new_x obj_val _userData = do+ toB $ poke obj_val =<< fun =<< ptrToVS n xin++wrapIpF2 fun n xin new_x obj_val _userData = do+ toB $ poke obj_val =<< fun =<< ptrToVS n xin++wrapIpG2 fun n xin new_x m gout _userData = do+ toB $ join $ liftM2 VM.copy (ptrToVS m gout) (fun =<< ptrToVS n xin)++wrapIpGradF2 fun n x new_x grad_f _userData = do+ toB $ join $ liftM2 VM.copy (ptrToVS n grad_f) (fun =<< ptrToVS n x)++wrapIpJacG2 fun1 fun2 n x new_x m nj iRow jCol jacs _userData+ | jacs == nullPtr = do+ toB $ join $ liftM2 fun1 (ptrToVS nj iRow) (ptrToVS nj jCol)+ | otherwise = do+ toB $ join $ liftM2 fun2 (ptrToVS n x) (ptrToVS nj jacs)++wrapIpH2 funSparsity funEval n x new_x obj_factor m lambda new_lambda nHess iRow jCol values _userData+ | iRow == nullPtr = do+ toB $ join $ liftM3 (funEval obj_factor)+ (ptrToVS m lambda)+ (ptrToVS n x)+ (ptrToVS nHess values)+ | otherwise = do+ toB $ join $ liftM2 funSparsity+ (ptrToVS nHess iRow)+ (ptrToVS nHess jCol)++wrapIpF f = wrapIpF1 (wrapIpF2 f)+wrapIpG f = wrapIpG1 (wrapIpG2 f)+wrapIpGradF f = wrapIpGradF1 (wrapIpGradF2 f)+wrapIpJacG f1 f2 = wrapIpJacG1 (wrapIpJacG2 f1 f2)+wrapIpH fSparsity fEval = wrapIpH1 (wrapIpH2 fSparsity fEval)+++vmUnsafeWith = VM.unsafeWith++-- | Vector of numbers+type Vec = VM.IOVector IpNumber++createIpoptProblem :: Vec -> Vec -> Vec -> Vec+ -> Int -> Int -> IpF -> IpG -> IpGradF -> IpJacG -> IpH -> IO IpProblem+createIpoptProblem xL xU gL gU nJac nHess f g gradF jacG hess+ | lx <- VM.length xL,+ lx == VM.length xU,+ lg <- VM.length gL,+ lg == VM.length gU = createIpoptProblem3 lx xL xU lg gL gU nJac nHess 0 f g gradF jacG hess+ | otherwise = error "dimensions wrong!"++{#fun CreateIpoptProblem as createIpoptProblem3+ { `Int', vmUnsafeWith* `Vec', vmUnsafeWith* `Vec',+ `Int', vmUnsafeWith* `Vec', vmUnsafeWith* `Vec',+ `Int', `Int', `Int', id `IpF', id `IpG', id `IpGradF',+ id `IpJacG', id `IpH' } -> `IpProblem' IpProblem #}+_ = {#fun AddIpoptNumOption as ^+ { unIpProblem `IpProblem', `String', `Double' } -> `Bool' #}++_ = {#fun AddIpoptStrOption as ^+ { unIpProblem `IpProblem', `String', `String' } -> `Bool' #}++_ = {#fun AddIpoptIntOption as ^+ { unIpProblem `IpProblem', `String', `Int' } -> `Bool' #}++_ = {#fun FreeIpoptProblem as ^+ { unIpProblem `IpProblem' } -> `()' #}++_ = {#fun OpenIpoptOutputFile as ^+ { unIpProblem `IpProblem', `String', `Int' } -> `Bool' #}++_ = {#fun SetIpoptProblemScaling as ^+ { unIpProblem `IpProblem',+ `Double',+ vmUnsafeWith* `Vec',+ vmUnsafeWith* `Vec'+ } -> `Bool' #}+++data IpOptSolved = IpOptSolved+ { ipOptSolved_status :: ApplicationReturnStatus,+ ipOptSolved_objective :: Double,+ ipOptSolved_x,+ ipOptSolved_g,+ ipOptSolved_mult_g,+ ipOptSolved_mult_x_L,+ ipOptSolved_mult_x_U :: Vec }++ipoptSolve :: IpProblem+ -> Vec -- ^ starting point @x@. Note that the value is overwritten with the final @x@.+ -> IO IpOptSolved+ipoptSolve problem x = do+ g <- VM.new (VM.length x)+ mult_g <- VM.new (VM.length x)+ mult_x_L <- VM.new (VM.length x)+ mult_x_U <- VM.new (VM.length x)++ out <- ipoptSolve2+ problem+ x+ g+ mult_g+ mult_x_L+ mult_x_U+ nullPtr+ return $ IpOptSolved+ (fst out)+ (snd out)+ x+ g+ mult_g+ mult_x_L+ mult_x_U++_ = {#fun IpoptSolve as ipoptSolve2+ { unIpProblem `IpProblem',+ vmUnsafeWith* `Vec',+ vmUnsafeWith* `Vec',+ alloca- `Double' peekFloatConv*,+ vmUnsafeWith* `Vec',+ vmUnsafeWith* `Vec',+ vmUnsafeWith* `Vec',+ id `Ptr ()' } -> `ApplicationReturnStatus' cToEnum #}+++{- | Set-up an 'IpProblem' to be solved later. Only objective function (@f@)+and constraint functions (@g@) need to be specified. Derivatives needed by ipopt+are computed by "Numeric.AD".++To solve the optimization problem:++> min f(x)+> such that+> xL <= x <= xU+> gL <= g(x) <= gU++First create an opaque 'IpProblem' object (nlp):++> nlp <- createIpOptProblemAD xL xU gL gU f g++Then pass it off to 'ipoptSolve'.++> ipoptSolve nlp x0++Refer to the example @Test2.hs@ for details of setting up the vectors supplied.+-}+createIpoptProblemAD+ :: 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)@+ -> (forall a. RealFloat a => V.Vector a -> a) -- ^ objective function @f : R^n -> R@+ -> (forall a. RealFloat a => V.Vector a -> V.Vector a) -- ^ constraint functions @g : R^n -> R^m@+ -> IO IpProblem+createIpoptProblemAD 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) <- mkFs (True,True) 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 -- ^ @m@ number of constraints+ -> (forall a. RealFloat a => V.Vector a -> a) -- ^ objective function @R^n -> R@+ -> (forall a. RealFloat 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+ ipF <- wrapIpF $ \x -> do+ x <- VG.convert `fmap` VS.unsafeFreeze x+ return $ f x++ ipGradF <- wrapIpGradF $ \x -> do+ x <- VG.convert `fmap` VS.unsafeFreeze x+ VS.unsafeThaw $ VG.convert (grad f x)++ ipG <- wrapIpG $ \x -> do+ 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+ x <- VG.convert `fmap` VS.unsafeFreeze x+ jac <- VS.unsafeThaw $ VG.convert $ VG.concat $ VG.toList $ 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 $ hessian f x+ gj = V.zipWith (\l v -> V.map (l*) v) lambda (V.map tri (hessianF g x))+ lagrangian = V.foldl (V.zipWith (+)) obj gj++ VM.copy values =<< VS.unsafeThaw (VG.convert lagrangian)+ )++ 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+ VM.copy iRow =<< VS.unsafeThaw (VS.generate (n*m) (\x -> fromIntegral $ x `div` n))+ VM.copy jCol =<< VS.unsafeThaw (VS.generate (n*m) (\x -> fromIntegral $ x `mod` n))++-- | indexes the same as http://www.coin-or.org/Ipopt/documentation/node41.html+denseIJh n m iRow jCol = do+ i <- newIORef 0+ forM_ [0 .. fromIntegral n-1] $ \ row ->+ forM_ [ 0 .. row ] $ \col -> do+ ii <- readIORef i+ VM.write iRow ii row+ VM.write jCol ii col+ writeIORef i (ii+1)+
+ Ipopt/Sparsity.hs view
@@ -0,0 +1,155 @@+{-# LANGUAGE PatternGuards, RankNTypes, TupleSections #-}+{- |++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.++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++{- $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. RealFloat 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. RealFloat 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. RealFloat 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. RealFloat 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. RealFloat 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. RealFloat 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 :: RealFloat a => V.Vector (Int,Int)+ -> (forall a. RealFloat 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+-}+
− Test1.hs
@@ -1,123 +0,0 @@-{- | Translation of the "number 71 from the Hock-Schittkowsky test suite"-<http://www.coin-or.org/Ipopt/documentation/node28.html> to check that the binding works.----}-module Main where---import Ipopt- ( wrapIpF,- wrapIpG,- wrapIpGradF,- wrapIpJacG,- wrapIpH,- createIpoptProblem,- addIpoptNumOption,- addIpoptStrOption,- ipoptSolve )-import qualified Data.Vector.Storable as VS-import qualified Data.Vector.Storable.Mutable as VM-import Data.Vector.Storable ( (!) )-import Control.Monad ( forM_, forM )-import Data.IORef ( writeIORef, readIORef, newIORef )-import Text.Printf ( printf )--main = do- let n = 4- m = 2-- fromList = VS.unsafeThaw . VS.fromList- xL <- fromList (replicate n 1)- xU <- fromList (replicate n 5)-- gL <- fromList [25, 40]- gU <- fromList [2.0e19,40]-- eval_f <- wrapIpF $ \y -> do- x <- VS.unsafeFreeze y- return (x!0 * x!3 * (x!0 + x!1 + x!2) + x!2)-- eval_grad_f <- wrapIpGradF $ \y -> do- x <- VS.unsafeFreeze y- fromList- [ x!0 * x!3 + x!3 * (x!0 + x!1 + x!2),- x!0 * x!3,- x!0 * x!3 + 1,- x!0 * (x!0 + x!1 + x!2)]-- eval_g <- wrapIpG $ \y -> do- x <- VS.unsafeFreeze y- fromList [x!0 * x!1 * x!2 * x!3,- x!0*x!0 + x!1*x!1 + x!2*x!2 + x!3*x!3]-- eval_jac_g <- wrapIpJacG- (\ iRow jCol -> do- print (VM.length iRow)- VM.set (VM.slice 0 4 iRow) 0- VM.set (VM.slice 4 4 iRow) 1- forM_ (zip [0 .. 7] $ cycle [0 .. 3]) $ \(a,b) ->- VM.write jCol a b- )- (\ x values -> do- x <- VS.unsafeFreeze x- values_new <- fromList [- x!1*x!2*x!3,- x!0*x!2*x!3,- x!0*x!1*x!3,- x!0*x!1*x!2,-- 2*x!0,- 2*x!1,- 2*x!2,- 2*x!3 ]- VM.copy values values_new- )--- eval_h <- wrapIpH- ( \ iRow jCol -> do- i <- newIORef 0- forM [0 .. 3] $ \ row ->- forM [ 0 .. row ] $ \col -> do- ii <- readIORef i- VM.write iRow ii row- VM.write jCol ii col- writeIORef i (ii+1)- )-- ( \ obj_factor lambda x values -> do-- lambda <- VS.unsafeFreeze lambda- x <- VS.unsafeFreeze x-- VM.copy values =<< fromList- [ obj_factor * (2*x!3) + lambda!1 * 2,- obj_factor * (x!3) + lambda!0 * (x!2 * x!3),- lambda!1 * 2,- obj_factor * (x!3) + lambda!0 * (x!1 * x!3),- lambda!0 * (x!0 * x!3),- lambda!1 * 2,- obj_factor * (2*x!0 + x!1 + x!2) + lambda!0 * (x!1 * x!2),- obj_factor * (x!0) + lambda!0 * (x!0 * x!2),- obj_factor * (x!0) + lambda!0 * (x!0 * x!1),- lambda!1 * 2 ]- )-- nlp <- createIpoptProblem xL xU gL gU 8 10 eval_f eval_g eval_grad_f eval_jac_g eval_h-- addIpoptNumOption nlp "tol" 1.0e-9- addIpoptStrOption nlp "mu_strategy" "adaptive"-- x <- fromList [1,5, 5,1]-- ans <- ipoptSolve nlp x-- let x_official = VS.fromList [1, 4.74299964, 3.82114998, 1.37940829]- x <- VS.freeze x-- let sse :: Double- sse = realToFrac $ VS.sum $ VS.zipWith (\a b -> (a-b)^2) x x_official-- printf "\n\n||x - x_official||_2 = %f\n" sse-
− Test2.hs
@@ -1,52 +0,0 @@-{-# LANGUAGE RankNTypes #-}-{- | Translation of the "number 71 from the Hock-Schittkowsky test suite"-<http://www.coin-or.org/Ipopt/documentation/node28.html> testing AD---}-module Main where--import qualified Data.Vector.Storable as VS-import Text.Printf ( printf )-import Data.Vector ((!))-import qualified Data.Vector as V-import qualified Data.Vector.Storable as VS-import Ipopt- ( addIpoptNumOption,- addIpoptStrOption,- ipoptSolve,- createIpoptProblemAD )--main = do- let n = 4 -- number of variables- fromList = VS.unsafeThaw . VS.fromList-- xL <- fromList (replicate n 1)- xU <- fromList (replicate n 5)-- gL <- fromList [25, 40]- gU <- fromList [2.0e19,40]-- let f x = x!0 * x!3 * (x!0 + x!1 + x!2) + x!2- g x = V.fromList- [x!0 * x!1 * x!2 * x!3,- x!0*x!0 + x!1*x!1 + x!2*x!2 + x!3*x!3]-- nlp <- createIpoptProblemAD xL xU gL gU f g-- addIpoptNumOption nlp "tol" 1.0e-9- addIpoptStrOption nlp "mu_strategy" "adaptive"-- x <- fromList [1,5, 5,1]-- ans <- ipoptSolve nlp x-- let x_official = VS.fromList [1, 4.74299964, 3.82114998, 1.37940829]-- x <- VS.freeze x-- let sse :: Double- sse = realToFrac $ VS.sum $ VS.zipWith (\a b -> (a-b)^2) x x_official-- printf "\n\n||x - x_official||_2 = %f\n" sse--
+ examples/Test1.hs view
@@ -0,0 +1,123 @@+{- | Translation of the "number 71 from the Hock-Schittkowsky test suite"+<http://www.coin-or.org/Ipopt/documentation/node28.html> to check that the binding works.+++-}+module Main where+++import Ipopt.Raw+ ( wrapIpF,+ wrapIpG,+ wrapIpGradF,+ wrapIpJacG,+ wrapIpH,+ createIpoptProblem,+ addIpoptNumOption,+ addIpoptStrOption,+ ipoptSolve )+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VM+import Data.Vector.Storable ( (!) )+import Control.Monad ( forM_, forM )+import Data.IORef ( writeIORef, readIORef, newIORef )+import Text.Printf ( printf )++main = do+ let n = 4+ m = 2++ fromList = VS.unsafeThaw . VS.fromList+ xL <- fromList (replicate n 1)+ xU <- fromList (replicate n 5)++ gL <- fromList [25, 40]+ gU <- fromList [2.0e19,40]++ eval_f <- wrapIpF $ \y -> do+ x <- VS.unsafeFreeze y+ return (x!0 * x!3 * (x!0 + x!1 + x!2) + x!2)++ eval_grad_f <- wrapIpGradF $ \y -> do+ x <- VS.unsafeFreeze y+ fromList+ [ x!0 * x!3 + x!3 * (x!0 + x!1 + x!2),+ x!0 * x!3,+ x!0 * x!3 + 1,+ x!0 * (x!0 + x!1 + x!2)]++ eval_g <- wrapIpG $ \y -> do+ x <- VS.unsafeFreeze y+ fromList [x!0 * x!1 * x!2 * x!3,+ x!0*x!0 + x!1*x!1 + x!2*x!2 + x!3*x!3]++ eval_jac_g <- wrapIpJacG+ (\ iRow jCol -> do+ print (VM.length iRow)+ VM.set (VM.slice 0 4 iRow) 0+ VM.set (VM.slice 4 4 iRow) 1+ forM_ (zip [0 .. 7] $ cycle [0 .. 3]) $ \(a,b) ->+ VM.write jCol a b+ )+ (\ x values -> do+ x <- VS.unsafeFreeze x+ values_new <- fromList [+ x!1*x!2*x!3,+ x!0*x!2*x!3,+ x!0*x!1*x!3,+ x!0*x!1*x!2,++ 2*x!0,+ 2*x!1,+ 2*x!2,+ 2*x!3 ]+ VM.copy values values_new+ )+++ eval_h <- wrapIpH+ ( \ iRow jCol -> do+ i <- newIORef 0+ forM [0 .. 3] $ \ row ->+ forM [ 0 .. row ] $ \col -> do+ ii <- readIORef i+ VM.write iRow ii row+ VM.write jCol ii col+ writeIORef i (ii+1)+ )++ ( \ obj_factor lambda x values -> do++ lambda <- VS.unsafeFreeze lambda+ x <- VS.unsafeFreeze x++ VM.copy values =<< fromList+ [ obj_factor * (2*x!3) + lambda!1 * 2,+ obj_factor * (x!3) + lambda!0 * (x!2 * x!3),+ lambda!1 * 2,+ obj_factor * (x!3) + lambda!0 * (x!1 * x!3),+ lambda!0 * (x!0 * x!3),+ lambda!1 * 2,+ obj_factor * (2*x!0 + x!1 + x!2) + lambda!0 * (x!1 * x!2),+ obj_factor * (x!0) + lambda!0 * (x!0 * x!2),+ obj_factor * (x!0) + lambda!0 * (x!0 * x!1),+ lambda!1 * 2 ]+ )++ nlp <- createIpoptProblem xL xU gL gU 8 10 eval_f eval_g eval_grad_f eval_jac_g eval_h++ addIpoptNumOption nlp "tol" 1.0e-9+ addIpoptStrOption nlp "mu_strategy" "adaptive"++ x <- fromList [1,5, 5,1]++ ans <- ipoptSolve nlp x++ let x_official = VS.fromList [1, 4.74299964, 3.82114998, 1.37940829]+ x <- VS.freeze x++ let sse :: Double+ sse = realToFrac $ VS.sum $ VS.zipWith (\a b -> (a-b)^2) x x_official++ printf "\n\n||x - x_official||_2 = %f\n" sse+
+ examples/Test2.hs view
@@ -0,0 +1,52 @@+{-# LANGUAGE RankNTypes #-}+{- | Translation of the "number 71 from the Hock-Schittkowsky test suite"+<http://www.coin-or.org/Ipopt/documentation/node28.html> testing AD++-}+module Main where++import qualified Data.Vector.Storable as VS+import Text.Printf ( printf )+import Data.Vector ((!))+import qualified Data.Vector as V+import qualified Data.Vector.Storable as VS+import Ipopt.Raw+ ( addIpoptNumOption,+ addIpoptStrOption,+ ipoptSolve,+ createIpoptProblemAD )++main = do+ let n = 4 -- number of variables+ fromList = VS.unsafeThaw . VS.fromList++ xL <- fromList (replicate n 1)+ xU <- fromList (replicate n 5)++ gL <- fromList [25, 40]+ gU <- fromList [2.0e19,40]++ let f x = x!0 * x!3 * (x!0 + x!1 + x!2) + x!2+ g x = V.fromList+ [x!0 * x!1 * x!2 * x!3,+ x!0*x!0 + x!1*x!1 + x!2*x!2 + x!3*x!3]++ nlp <- createIpoptProblemAD xL xU gL gU f g++ addIpoptNumOption nlp "tol" 1.0e-9+ addIpoptStrOption nlp "mu_strategy" "adaptive"++ x <- fromList [1,5, 5,1]++ ans <- ipoptSolve nlp x++ let x_official = VS.fromList [1, 4.74299964, 3.82114998, 1.37940829]++ x <- VS.freeze x++ let sse :: Double+ sse = realToFrac $ VS.sum $ VS.zipWith (\a b -> (a-b)^2) x x_official++ printf "\n\n||x - x_official||_2 = %f\n" sse++
+ examples/Test3.hs view
@@ -0,0 +1,20 @@+module Main where+import Ipopt+import qualified Data.Vector as V+import Control.Lens+import Control.Monad+import Control.Monad.State++-- | the same problem as solved in Test1 and Test2, except+-- showing how bits from Ipopt.NLP help+test22 = do+ [x0, x1, x2, x3] <- replicateM 4 (varFresh (Just (1,5)) "x")+ addF (Just "test22") $ x0*x3*(x0 + x1 + x2) + x2+ addG Nothing (25,2.0e19) $ x0*x1*x2*x3+ addG Nothing (40,40) $ x0*x0 + x1*x1 + x2*x2 + x3*x3+ initX .= V.fromList [1,5,5,1]+ solveNLP' (const (return ()))+ `runStateT` nlpstate0++main = test22+
ipopt-hs.cabal view
@@ -1,25 +1,75 @@ name: ipopt-hs-version: 0.0.0.0+version: 0.2.0.0 synopsis: haskell binding to ipopt including automatic differentiation-description: +description: a haskell binding to the nonlinear programming solver ipopt+ <http://projects.coin-or.org/Ipopt>+ .+ [@installation@]+ needs ipopt installed: For example if you have a+ `/usr/include/coin/IpStdCInterface.h` from your ipopt installation,+ use:+ .+ > cabal install ipopt-hs --extra-include-dirs=/usr/include/coin+ .+ 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.+ .+ Refer to @examples/Test1.hs@ for an example where the derivatives+ are computed by hand, @Test2.hs@ for the use of+ 'createIpoptProblemAD' and @Test3.hs@ for the highest level.+ .+ Current limitations include:+ .+ * 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.+ .+ * sparseness of derivatives isn't used to decide which way to calculate (forward vs. backward mode)+ .+ * probably doesn't work if @IpStdCInterface.h@ has Number =/= 'CDouble'+ .+ * no binding to SetIntermediateCallback+ .+ * garbage collection of 'IpProblem' won't free C-side resources license: BSD3 license-file: LICENSE author: Adam Vogt <vogt.adam@gmail.com> maintainer: Adam Vogt <vogt.adam@gmail.com>-category: Data+category: Optimisation Math Numeric build-type: Simple cabal-version: >=1.8 +source-repository head+ type: darcs+ location: http://code.haskell.org/~aavogt/ipopt-hs+ library- exposed-modules: Ipopt+ exposed-modules: Ipopt, Ipopt.Raw, Ipopt.Sparsity, Ipopt.NLP other-modules: C2HS- build-depends: base ==4.6.*, vector ==0.10.*, ad ==3.4.*+ build-depends: base ==4.6.*,+ vector ==0.10.*,+ ad ==3.4.*,+ containers,+ mtl == 2.*,+ lens pkgconfig-depends: ipopt+ build-tools: c2hs -executable ipopt-hs_Test1 +executable ipopt-hs_Test1 main-is: Test1.hs build-depends: base ==4.6.*, vector ==0.10.*, ipopt-hs+ hs-source-dirs: examples+ other-modules: Paths_ipopt_hs executable ipopt-hs_Test2 main-is: Test2.hs build-depends: base ==4.6.*, vector ==0.10.*, ipopt-hs+ hs-source-dirs: examples++executable ipopt-hs_Test3+ main-is: Test3.hs+ build-depends: base ==4.6.*, vector ==0.10.*, ipopt-hs, lens, mtl+ hs-source-dirs: examples