ipopt-hs 0.2.0.0 → 0.3.0.0
raw patch · 11 files changed
+410/−357 lines, 11 filesdep +Rlang-QQdep +ansi-wl-pprintdep +splinesdep ~addep ~basedep ~containersnew-component:exe:ipopt-hs_Tests
Dependencies added: Rlang-QQ, ansi-wl-pprint, splines, vector-space
Dependency ranges changed: ad, base, containers, lens
Files
- Ipopt.hs +15/−11
- Ipopt/AnyRF.hs +127/−0
- Ipopt/NLP.hs +92/−100
- Ipopt/PP.hs +64/−0
- Ipopt/Raw.chs +38/−22
- Ipopt/Sparsity.hs +13/−10
- examples/AllTests.hs +23/−0
- examples/Test1.hs +0/−123
- examples/Test2.hs +0/−52
- examples/Test3.hs +0/−20
- ipopt-hs.cabal +38/−19
Ipopt.hs view
@@ -4,33 +4,37 @@ -- 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(..),+ AnyRF(..), Identity(..), -- ** functions addG, addF,- -- * low-level bits still needed- IpOptSolved(..),- ApplicationReturnStatus(..),-- -- ** solver options+ -- ** running the solver+ ppSoln, + NLPT, nlpstate0,+ module Control.Monad.State,+ solveNLP',+ -- *** solver options -- $solverOptRef addIpoptNumOption, addIpoptStrOption, addIpoptIntOption, setIpoptProblemScaling, openIpoptOutputFile,-- -- *** types+ + -- * low-level bits still needed+ IpOptSolved(..),+ ApplicationReturnStatus(..),+ -- ** types Vec, IpNumber, ) where +import Ipopt.PP import Ipopt.NLP import Ipopt.Raw+import Ipopt.AnyRF import Control.Monad.State+import Control.Monad.Identity (Identity(..)) -- $solverOptRef -- see <http://www.coin-or.org/Ipopt/documentation/node39.html ipopt's options reference>
+ Ipopt/AnyRF.hs view
@@ -0,0 +1,127 @@+{- | Description: representing functions that can be differentiated++The 'AnyRF' 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'.++Values of type @AnyRF Identity@ can be generated using functions +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.+-}+module Ipopt.AnyRF where++import Data.Sequence (Seq)+import Data.Vector (Vector)+import Data.Monoid+import Control.Monad.Identity+import qualified Data.VectorSpace as VectorSpace+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++-- | @AnyRF cb@ is a function that uses variables from the nonlinear+-- program in a way supported by 'AnyRFCxt'. The @cb@ is+-- usually 'Identity'+data AnyRF cb = AnyRF (forall a. AnyRFCxt a => Vector a -> cb a)++-- | RealFloat gives most numerical operations,+-- 'VectorSpace' is involved to allow using definitions from the splines package+type AnyRFCxt a = (RealFloat a, VectorSpace a, Scalar a ~ a)++-- *** helpers for defining instances+liftOp0 :: (forall a. AnyRFCxt a => a) -> AnyRF Identity+liftOp0 op = AnyRF $ \x -> Identity op++liftOp1 :: (forall a. AnyRFCxt a => a -> a) -> AnyRF Identity -> AnyRF Identity+liftOp1 op (AnyRF a) = AnyRF $ \x -> Identity (op (runIdentity (a x)))++liftOp2 :: (forall a. AnyRFCxt 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) where+ compare _ = error "anyRF compare"+ max = liftOp2 max+ min = liftOp2 min+instance Eq (AnyRF Identity) where+ (==) = error "anyRF =="+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 VectorSpace.VectorSpace (AnyRF Identity) where+ type Scalar (AnyRF Identity) = Double+ x *^ v = realToFrac x*v++instance VectorSpace.AdditiveGroup (AnyRF Identity) where+ zeroV = liftOp0 0+ (^+^) = (+)+ negateV = negate++-- * orphan instances +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.*!)
Ipopt/NLP.hs view
@@ -1,9 +1,14 @@-{-# LANGUAGE FlexibleInstances, GADTs, RankNTypes, TemplateHaskell #-}+{-# LANGUAGE TemplateHaskell #-} -- | Description : a EDSL for describing nonlinear programs -- -- see usage in @examples/Test3.hs@+--+-- IPOPT does support naming variables if you use c+++-- (by overriding a @virtual void finalize_metadata@), but+-- it's not clear that we can set that from c/haskell module Ipopt.NLP where +import Ipopt.AnyRF import Control.Applicative import Control.Lens import Control.Monad@@ -23,8 +28,13 @@ import qualified Data.Sequence as Seq import qualified Data.Set as S import qualified Data.Vector as V+import qualified Data.Vector.Generic as VG+import Data.Vector ((!)) import qualified Data.Vector.Storable as VS+import Data.Maybe +import qualified Data.VectorSpace as VectorSpace+ import Ipopt.Raw -- * state@@ -40,10 +50,12 @@ -- | what namespace (see 'inEnv') _currentEnv :: [String], _variables :: M.Map String Ix,+ _variablesInv :: IxMap String, -- invert variables map -- | human-readable descriptions for the constraint, objective and -- variables- _constraintLabels, _objLabels, _varLabels :: IntMap String,- _varEnv :: IntMap (S.Set [String]),+ _constraintLabels, _objLabels :: IntMap String,+ _varLabels :: IxMap String,+ _varEnv :: IxMap (S.Set [String]), _constraintEnv, _objEnv :: IntMap [String], _nlpfun :: NLPFun, _defaultBounds :: (Double,Double),@@ -51,14 +63,19 @@ _initX :: Vector Double } deriving (Show) --- | solver deals with arrays. This type is for indexes into the array+newtype IxMap a = IxMap (IntMap a)+ deriving (Show, Monoid, Functor)++-- | the 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 mempty -- env (mempty :: NLPFun) (-1/0, 1/0)@@ -67,90 +84,7 @@ 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@@ -159,7 +93,13 @@ makeLenses ''NLPFun makeLenses ''Ix makeLenses ''NLPState+makeIso ''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)+ix_ (Ix i) = from ixMap . ix i+at_ (Ix i) = from ixMap . at i+ -- | @to@ is one of 'varEnv', 'constraintEnv', 'objEnv' copyEnv to n = do ev <- use currentEnv@@ -174,9 +114,9 @@ -- | calls 'createIpoptProblemAD' and 'ipoptSolve'. To be used at the -- end of a do-block.-solveNLP' :: MonadIO m =>+solveNLP' :: (VG.Vector v Double, MonadIO m) => (IpProblem -> IO ()) -- ^ set ipopt options (using functions from "Ipopt.Raw")- -> NLPT m IpOptSolved+ -> NLPT m (IpOptSolved v) solveNLP' setOpts = do (xl,xu) <- join $ uses (nlpfun . boundX) seqToVecs (gl,gu) <- join $ uses (nlpfun . boundG) seqToVecs@@ -223,13 +163,12 @@ 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'+ -> Maybe String -- ^ optional longer description -> 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+ -> NLPT m Ix -- ^ the index (into the rawvector of variables that the solver sees)+var' bs longDescription s = do ev <- use currentEnv m <- use variables let s' = intercalate "." (reverse (s:ev))@@ -240,14 +179,24 @@ nlpfun . boundX %= (Seq.|> db) n' <- use nMax variables %= M.insert s' n'+ variablesInv . at_ n' .= Just s'+ F.for_ longDescription $ \d -> varLabels . at_ n' %= (<> Just d)++ -- try to find a sane initial value+ initX %= let x0 = fromMaybe 0 $+ bs >>= \(a,b) -> find valid [(a+b)/2, a, b]+ valid x = not $ isInfinite x || isNaN x+ in (`V.snoc` x0) return n' Just n -> return n- varEnv %= IM.insert (view varIx n) (S.singleton ev)+ varEnv . at_ n %= (<> Just (S.singleton ev))+ F.for_ longDescription $ \str -> varLabels . at_ n %= (<> Just str)+ F.traverse_ (narrowBounds n) bs- return (AnyRF $ \x -> Identity $ x V.! view varIx n, n)+ return n --- | 'var'' without the usually unnecessary 'Ix'-var bs s = fmap fst (var' bs s)+-- | a combination of 'var'' and from'Ix'+var bs s = ixToVar <$> var' bs Nothing s {- | 'var', except this causes the solver to get a new variable, so that you can use:@@ -258,10 +207,20 @@ 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)+varFresh' :: (Monad m, Functor m) => Maybe (Double,Double) -> String -> NLPT m Ix+varFresh' bs s = do+ existing <- gets (^? variables . ix s)+ case existing of+ Just a -> return a+ Nothing -> do+ m <- use variables+ let n = M.size m + 1+ -- get the first of "x", "x1", "x1_", "x1__", "x1___"+ Just sUniq <- return $ find (`M.notMember` m) $ s : iterate (++"_") (s ++ show n)+ var' bs Nothing sUniq +varFresh bs s = fmap ixToVar $ varFresh' bs s+ -- *** namespace {- $namespace @@ -296,6 +255,39 @@ n : ns <- use currentEnv currentEnv .= ns return n++-- *** piecewise linear++-- $piecewise+-- see for example chapter 20 of <http://www.ampl.com/BOOK/download.html>+-- and use of the splines package in @examples\/Test4@ and @examples\/Test5@++{- | splits a variable @x@ into two positive variables such that+@x = x^+ - x^-@ the new variables represent the positive and negative+parts of @x - b@++> (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+-}+splitVar :: (Monad m, Functor m)+ => Double -- ^ @b@+ -> Ix -- ^ index for @x@+ -> NLPT m (AnyRF Identity, AnyRF Identity) -- ^ @(b-x)_+, (x-b)_+@+splitVar b i = do+ -- need to have a variable name... + Just s <- gets (^? variablesInv . ix_ i)+ let x = ixToVar i+ xPlus <- varFresh' (Just (0, 1/0)) (s ++ "+")+ xMinus <- varFresh' (Just (0, 1/0)) (s ++ "-")+ x0 <- uses initX (! view varIx i)+ initX %= (V.// [(view varIx xPlus, max 0 (x0 - b)), (view varIx xMinus, max 0 ( b - x0) )])+ addG (Just ("splitVar: " ++ s)) (b, b) (x + ixToVar xMinus - ixToVar xPlus)+ return (ixToVar xMinus, ixToVar xPlus)++ixToVar :: Ix -> AnyRF Identity+ixToVar (Ix i) = AnyRF (\v -> Identity (v V.! i)) -- *** bounds
+ Ipopt/PP.hs view
@@ -0,0 +1,64 @@+{-# LANGUAGE TemplateHaskell, NoMonomorphismRestriction, OverloadedStrings #-}+module Ipopt.PP where++import Data.List+import Data.Ord+import Ipopt.Raw+import Ipopt.NLP+import Ipopt.AnyRF+import Text.PrettyPrint.ANSI.Leijen hiding ((<>))++import Text.Printf+import qualified Data.IntMap as IM+import qualified Data.Vector.Storable as VS+import Data.Vector (Vector)+import qualified Data.Vector as V+import Control.Monad.State+import Control.Monad.Writer+import Foreign.C.Types(CDouble(..))+import Control.Lens+import Data.List+import Data.Foldable (toList,for_, foldMap)++-- * lenses for IpOptSolved+makeLensesWith (lensRules & lensField .~ stripPrefix "_ipOptSolved_") ''IpOptSolved++-- * pretty printing+ppSoln :: (MonadIO m, Functor m) =>+ NLPState -> NLPT m (IpOptSolved Vector) -> m (IpOptSolved Vector, Doc)+ppSoln state0 problem = flip evalStateT state0 $ runWriterT $ do+ s <- lift problem+ st <- get++ tell (dullgreen "status: " <> (s^.status.to colorStatus))+ br+ tell $ "obj_tot" <> double (s^.objective)+ join $ uses (nlpfun . funF) $ \(AnyRF f) -> case sortBy (comparing fst) $+ toList (f (s^.x)) `zip` [1 .. ] of+ [_] -> return ()+ [] -> return ()+ xs -> for_ xs $ \(x,i) -> tell $ "obj" <> int i <> colon <$> double x++ for_ (st ^. variablesInv . from ixMap . to IM.toList) $ \(k,desc) -> do+ br+ tell $ string desc <> "(" <> int k <> ")" <> "=" <> string (printf "%.3g" (s ^?! x . ix k))+ br+ tell $ "g" <$> foldMap (\e -> mempty <$$> double e) (s ^. g)++ return s++-- * internal++statusOk :: ApplicationReturnStatus -> Bool+statusOk x = case x of+ SolveSucceeded -> True+ SolvedToAcceptableLevel -> True+ UserRequestedStop -> True+ FeasiblePointFound -> True+ _ -> False++colorStatus x = (if statusOk x then id else black . onred) (string (show x))++br = tell (mempty <$$> mempty)++
Ipopt/Raw.chs view
@@ -1,4 +1,4 @@-{-# LANGUAGE ForeignFunctionInterface, PatternGuards, RankNTypes, TypeFamilies #-}+{-# LANGUAGE ForeignFunctionInterface #-} {- | Copyright: (C) 2013 Adam Vogt@@ -84,6 +84,9 @@ import qualified Data.Vector.Storable as VS import qualified Data.Vector.Storable.Mutable as VM +import Ipopt.AnyRF+import Data.VectorSpace (VectorSpace,Scalar)+ #include "IpStdCInterface.h" type IpNumber = {# type Number #}@@ -165,6 +168,12 @@ -- | Vector of numbers type Vec = VM.IOVector IpNumber +fromVec :: VG.Vector v Double => Vec -> IO (v Double)+fromVec mv = do+ v <- VS.freeze mv + return (VG.convert (VS.map (\(CDouble a) -> a) v))++ 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@@ -201,19 +210,19 @@ 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 }+-- | lenses are in "Ipopt.PP"+data IpOptSolved v = IpOptSolved+ { _ipOptSolved_status :: ApplicationReturnStatus,+ _ipOptSolved_objective :: Double,+ _ipOptSolved_x,+ _ipOptSolved_g,+ _ipOptSolved_mult_g,+ _ipOptSolved_mult_x_L,+ _ipOptSolved_mult_x_U :: v Double } -ipoptSolve :: IpProblem+ipoptSolve :: VG.Vector v Double => IpProblem -> Vec -- ^ starting point @x@. Note that the value is overwritten with the final @x@.- -> IO IpOptSolved+ -> IO (IpOptSolved v) ipoptSolve problem x = do g <- VM.new (VM.length x) mult_g <- VM.new (VM.length x)@@ -228,14 +237,21 @@ mult_x_L mult_x_U nullPtr+ + x' <- fromVec x+ g' <- fromVec g+ mult_g' <- fromVec mult_g + mult_x_L' <- fromVec mult_x_L+ mult_x_U' <- fromVec mult_x_U+ return $ IpOptSolved (fst out) (snd out)- x- g- mult_g- mult_x_L- mult_x_U+ x'+ g'+ mult_g'+ mult_x_L'+ mult_x_U' _ = {#fun IpoptSolve as ipoptSolve2 { unIpProblem `IpProblem',@@ -274,15 +290,15 @@ -> 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@+ -> (forall a. AnyRFCxt a => V.Vector a -> a) -- ^ objective function @f : R^n -> R@+ -> (forall a. AnyRFCxt 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+ (eval_f, eval_grad_f, eval_g, eval_jac_g, eval_h) <- mkFs (False,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 @@ -290,8 +306,8 @@ (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@+ -> (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 ipF <- wrapIpF $ \x -> do
Ipopt/Sparsity.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE PatternGuards, RankNTypes, TupleSections #-} {- | Description: find out derivatives are unnecessary@@ -6,7 +5,9 @@ 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.+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) @@ -47,6 +48,8 @@ import qualified Data.IntSet as S import Data.Monoid +import Ipopt.AnyRF+ {- $setup >>> let g x = x!0 + x!1 * x!1 * x!2 @@ -62,7 +65,7 @@ -} nanPropagateG1 :: Int -- ^ size of input vector- -> (forall a. RealFloat a => V.Vector a -> a)+ -> (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@@ -80,7 +83,7 @@ -} nanPropagateG :: Int -- ^ size of input vector- -> (forall a. RealFloat a => V.Vector a -> a)+ -> (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)))@@ -92,7 +95,7 @@ variable 3 isn't even used. -} nanPropagate1 :: Int -- ^ size of input vector- -> (forall a. RealFloat a => V.Vector a -> a)+ -> (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],@@ -106,17 +109,17 @@ -} nanPropagateH :: Int- -> (forall a. RealFloat a => V.Vector a -> a)+ -> (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. RealFloat a => V.Vector a -> V.Vector a)+ -> (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. RealFloat a => V.Vector a -> V.Vector a)+ -> (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))@@ -138,8 +141,8 @@ -jacobianSS :: RealFloat a => V.Vector (Int,Int)- -> (forall a. RealFloat a => V.Vector a -> V.Vector a)+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
+ examples/AllTests.hs view
@@ -0,0 +1,23 @@+import Test1+import Test2+import Test3+import Test4+import Test5+import System.Environment++-- probably should get autogenerated if lots of examples get generated...+main = do+ as <- getArgs+ case as of+ [] -> error "run like: ipopt-hs_Tests 1\nipopt-hs_Tests all"+ ["1"] -> main1+ ["2"] -> main2+ ["3"] -> main3+ ["4"] -> main4+ ["5"] -> main5+ ["all"] -> do+ main1+ main2+ main3+ main4+ main5
− examples/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.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
@@ -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.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
@@ -1,20 +0,0 @@-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,5 +1,5 @@ name: ipopt-hs-version: 0.2.0.0+version: 0.3.0.0 synopsis: haskell binding to ipopt including automatic differentiation description: a haskell binding to the nonlinear programming solver ipopt <http://projects.coin-or.org/Ipopt>@@ -38,38 +38,57 @@ license-file: LICENSE author: Adam Vogt <vogt.adam@gmail.com> maintainer: Adam Vogt <vogt.adam@gmail.com>-category: Optimisation Math Numeric+category: Optimisation, Math, Numeric build-type: Simple-cabal-version: >=1.8+cabal-version: >=1.10 +flag build_examples+ description: build executable from examples/+ default: False+ source-repository head type: darcs location: http://code.haskell.org/~aavogt/ipopt-hs library- exposed-modules: Ipopt, Ipopt.Raw, Ipopt.Sparsity, Ipopt.NLP+ exposed-modules: Ipopt,+ Ipopt.AnyRF,+ Ipopt.Raw,+ Ipopt.Sparsity,+ Ipopt.NLP,+ Ipopt.PP other-modules: C2HS- build-depends: base ==4.6.*,+ build-depends: base < 5, vector ==0.10.*, ad ==3.4.*,- containers,+ containers == 0.5.*, mtl == 2.*,- lens+ lens >= 3.10,+ ansi-wl-pprint >= 0.6.7,+ vector-space >= 0.8.6+ default-language: Haskell2010+ default-extensions: ConstraintKinds,+ FlexibleContexts,+ FlexibleInstances,+ GeneralizedNewtypeDeriving,+ RankNTypes,+ TupleSections,+ TypeFamilies+ other-extensions: ForeignFunctionInterface,+ TemplateHaskell pkgconfig-depends: ipopt build-tools: c2hs -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+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+ hs-source-dirs: examples+ default-language: Haskell2010+ other-modules: Paths_ipopt_hs+ if !flag(build_examples)+ buildable: False