diff --git a/C2HS.hs b/C2HS.hs new file mode 100644 --- /dev/null +++ b/C2HS.hs @@ -0,0 +1,238 @@ +-- C->Haskell Compiler: Marshalling library +-- +-- Copyright (c) [1999...2005] Manuel M T Chakravarty +-- +-- Redistribution and use in source and binary forms, with or without +-- modification, are permitted provided that the following conditions are met: +-- +-- 1. Redistributions of source code must retain the above copyright notice, +-- this list of conditions and the following disclaimer. +-- 2. Redistributions in binary form must reproduce the above copyright +-- notice, this list of conditions and the following disclaimer in the +-- documentation and/or other materials provided with the distribution. +-- 3. The name of the author may not be used to endorse or promote products +-- derived from this software without specific prior written permission. +-- +-- THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR +-- IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +-- OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN +-- NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +-- SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +-- TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +-- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +-- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +-- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +-- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +-- +--- Description --------------------------------------------------------------- +-- +-- Language: Haskell 98 +-- +-- This module provides the marshaling routines for Haskell files produced by +-- C->Haskell for binding to C library interfaces. It exports all of the +-- low-level FFI (language-independent plus the C-specific parts) together +-- with the C->HS-specific higher-level marshalling routines. +-- + +module C2HS ( + + -- * Re-export the language-independent component of the FFI + module Foreign, + + -- * Re-export the C language component of the FFI + module Foreign.C, + + -- * Composite marshalling functions + withCStringLenIntConv, peekCStringLenIntConv, withIntConv, withFloatConv, + peekIntConv, peekFloatConv, withBool, peekBool, withEnum, peekEnum, + + -- * Conditional results using 'Maybe' + nothingIf, nothingIfNull, + + -- * Bit masks + combineBitMasks, containsBitMask, extractBitMasks, + + -- * Conversion between C and Haskell types + cIntConv, cFloatConv, cToBool, cFromBool, cToEnum, cFromEnum +) where + + +import Foreign +import Foreign.C + +import Control.Monad (liftM) + + +-- Composite marshalling functions +-- ------------------------------- + +-- Strings with explicit length +-- +withCStringLenIntConv :: Num n => String -> ((CString, n) -> IO a) -> IO a +withCStringLenIntConv s f = withCStringLen s $ \(p, n) -> f (p, fromIntegral n) + +peekCStringLenIntConv :: Integral n => (CString, n) -> IO String +peekCStringLenIntConv (s, n) = peekCStringLen (s, fromIntegral n) + +-- Marshalling of numerals +-- + +withIntConv :: (Storable b, Integral a, Integral b) + => a -> (Ptr b -> IO c) -> IO c +withIntConv = with . fromIntegral + +withFloatConv :: (Storable b, RealFloat a, RealFloat b) + => a -> (Ptr b -> IO c) -> IO c +withFloatConv = with . realToFrac + +peekIntConv :: (Storable a, Integral a, Integral b) + => Ptr a -> IO b +peekIntConv = liftM fromIntegral . peek + +peekFloatConv :: (Storable a, RealFloat a, RealFloat b) + => Ptr a -> IO b +peekFloatConv = liftM realToFrac . peek + + +-- Everything else below is deprecated. +-- These functions are not used by code generated by c2hs. + +{-# DEPRECATED withBool "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED peekBool "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED withEnum "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED peekEnum "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED nothingIf "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED nothingIfNull "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED combineBitMasks "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED containsBitMask "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED extractBitMasks "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED cIntConv "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED cFloatConv "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED cFromBool "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED cToBool "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED cToEnum "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} +{-# DEPRECATED cFromEnum "The C2HS module will soon stop providing unnecessary\nutility functions. Please use standard FFI library functions instead." #-} + + +-- Passing Booleans by reference +-- + +withBool :: (Integral a, Storable a) => Bool -> (Ptr a -> IO b) -> IO b +withBool = with . fromBool + +peekBool :: (Integral a, Storable a) => Ptr a -> IO Bool +peekBool = liftM toBool . peek + + +-- Passing enums by reference +-- + +withEnum :: (Enum a, Integral b, Storable b) => a -> (Ptr b -> IO c) -> IO c +withEnum = with . cFromEnum + +peekEnum :: (Enum a, Integral b, Storable b) => Ptr b -> IO a +peekEnum = liftM cToEnum . peek + + +-- Storing of 'Maybe' values +-- ------------------------- + +--TODO: kill off this orphan instance! + +instance Storable a => Storable (Maybe a) where + sizeOf _ = sizeOf (undefined :: Ptr ()) + alignment _ = alignment (undefined :: Ptr ()) + + peek p = do + ptr <- peek (castPtr p) + if ptr == nullPtr + then return Nothing + else liftM Just $ peek ptr + + poke p v = do + ptr <- case v of + Nothing -> return nullPtr + Just v' -> new v' + poke (castPtr p) ptr + + +-- Conditional results using 'Maybe' +-- --------------------------------- + +-- Wrap the result into a 'Maybe' type. +-- +-- * the predicate determines when the result is considered to be non-existing, +-- ie, it is represented by `Nothing' +-- +-- * the second argument allows to map a result wrapped into `Just' to some +-- other domain +-- +nothingIf :: (a -> Bool) -> (a -> b) -> a -> Maybe b +nothingIf p f x = if p x then Nothing else Just $ f x + +-- |Instance for special casing null pointers. +-- +nothingIfNull :: (Ptr a -> b) -> Ptr a -> Maybe b +nothingIfNull = nothingIf (== nullPtr) + + +-- Support for bit masks +-- --------------------- + +-- Given a list of enumeration values that represent bit masks, combine these +-- masks using bitwise disjunction. +-- +combineBitMasks :: (Num b, Enum a, Bits b) => [a] -> b +combineBitMasks = foldl (.|.) 0 . map (fromIntegral . fromEnum) + +-- Tests whether the given bit mask is contained in the given bit pattern +-- (i.e., all bits set in the mask are also set in the pattern). +-- +containsBitMask :: (Num a, Bits a, Enum b) => a -> b -> Bool +bits `containsBitMask` bm = let bm' = fromIntegral . fromEnum $ bm + in + bm' .&. bits == bm' + +-- |Given a bit pattern, yield all bit masks that it contains. +-- +-- * This does *not* attempt to compute a minimal set of bit masks that when +-- combined yield the bit pattern, instead all contained bit masks are +-- produced. +-- +extractBitMasks :: (Num a, Bits a, Enum b, Bounded b) => a -> [b] +extractBitMasks bits = + [bm | bm <- [minBound..maxBound], bits `containsBitMask` bm] + + +-- Conversion routines +-- ------------------- + +-- |Integral conversion +-- +cIntConv :: (Integral a, Integral b) => a -> b +cIntConv = fromIntegral + +-- |Floating conversion +-- +cFloatConv :: (RealFloat a, RealFloat b) => a -> b +cFloatConv = realToFrac + +-- |Obtain C value from Haskell 'Bool'. +-- +cFromBool :: Num a => Bool -> a +cFromBool = fromBool + +-- |Obtain Haskell 'Bool' from C value. +-- +cToBool :: (Eq a, Num a) => a -> Bool +cToBool = toBool + +-- |Convert a C enumeration to Haskell. +-- +cToEnum :: (Integral i, Enum e) => i -> e +cToEnum = toEnum . fromIntegral + +-- |Convert a Haskell enumeration to C. +-- +cFromEnum :: (Enum e, Integral i) => e -> i +cFromEnum = fromIntegral . fromEnum diff --git a/Ipopt.chs b/Ipopt.chs new file mode 100644 --- /dev/null +++ b/Ipopt.chs @@ -0,0 +1,359 @@ +{-# LANGUAGE ForeignFunctionInterface, PatternGuards, RankNTypes, TypeFamilies #-} +{- | + +Copyright: (C) 2013 Adam Vogt +Maintainer: Adam Vogt +Stability: unstable + +Binding to 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) + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2013, Adam Vogt + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Adam Vogt nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/Test1.hs b/Test1.hs new file mode 100644 --- /dev/null +++ b/Test1.hs @@ -0,0 +1,123 @@ +{- | Translation of the "number 71 from the Hock-Schittkowsky test suite" + 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 + diff --git a/Test2.hs b/Test2.hs new file mode 100644 --- /dev/null +++ b/Test2.hs @@ -0,0 +1,52 @@ +{-# 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 + ( 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/ipopt-hs.cabal b/ipopt-hs.cabal new file mode 100644 --- /dev/null +++ b/ipopt-hs.cabal @@ -0,0 +1,25 @@ +name: ipopt-hs +version: 0.0.0.0 +synopsis: haskell binding to ipopt including automatic differentiation +description: +license: BSD3 +license-file: LICENSE +author: Adam Vogt +maintainer: Adam Vogt +category: Data +build-type: Simple +cabal-version: >=1.8 + +library + exposed-modules: Ipopt + other-modules: C2HS + build-depends: base ==4.6.*, vector ==0.10.*, ad ==3.4.* + pkgconfig-depends: ipopt + +executable ipopt-hs_Test1 + main-is: Test1.hs + build-depends: base ==4.6.*, vector ==0.10.*, ipopt-hs + +executable ipopt-hs_Test2 + main-is: Test2.hs + build-depends: base ==4.6.*, vector ==0.10.*, ipopt-hs