diff --git a/Ipopt.hs b/Ipopt.hs --- a/Ipopt.hs +++ b/Ipopt.hs @@ -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 diff --git a/Ipopt/AnyRF.hs b/Ipopt/AnyRF.hs new file mode 100644 --- /dev/null +++ b/Ipopt/AnyRF.hs @@ -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.*!) diff --git a/Ipopt/NLP.hs b/Ipopt/NLP.hs --- a/Ipopt/NLP.hs +++ b/Ipopt/NLP.hs @@ -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 +-- 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 diff --git a/Ipopt/PP.hs b/Ipopt/PP.hs new file mode 100644 --- /dev/null +++ b/Ipopt/PP.hs @@ -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) + + diff --git a/Ipopt/Raw.chs b/Ipopt/Raw.chs --- a/Ipopt/Raw.chs +++ b/Ipopt/Raw.chs @@ -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 diff --git a/Ipopt/Sparsity.hs b/Ipopt/Sparsity.hs --- a/Ipopt/Sparsity.hs +++ b/Ipopt/Sparsity.hs @@ -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 diff --git a/examples/AllTests.hs b/examples/AllTests.hs new file mode 100644 --- /dev/null +++ b/examples/AllTests.hs @@ -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 diff --git a/examples/Test1.hs b/examples/Test1.hs deleted file mode 100644 --- a/examples/Test1.hs +++ /dev/null @@ -1,123 +0,0 @@ -{- | Translation of the "number 71 from the Hock-Schittkowsky test suite" - 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 - diff --git a/examples/Test2.hs b/examples/Test2.hs deleted file mode 100644 --- a/examples/Test2.hs +++ /dev/null @@ -1,52 +0,0 @@ -{-# LANGUAGE RankNTypes #-} -{- | Translation of the "number 71 from the Hock-Schittkowsky test suite" - 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 - - diff --git a/examples/Test3.hs b/examples/Test3.hs deleted file mode 100644 --- a/examples/Test3.hs +++ /dev/null @@ -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 - 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.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 @@ -38,38 +38,57 @@ license-file: LICENSE author: Adam Vogt maintainer: Adam Vogt -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