module Ipopt.Raw (
createIpoptProblemAD,
ipoptSolve,
IpOptSolved(..),
addIpoptNumOption,
addIpoptStrOption,
addIpoptIntOption,
openIpoptOutputFile,
Vec,
IpNumber(..),
IpIndex(..),
IpInt(..),
IpBool(..),
IpF(..),
IpGradF(..),
IpG(..),
IpJacG(..),
IpH(..),
IpProblem(..),
IntermediateCB,
ApplicationReturnStatus(..),
createIpoptProblem,
freeIpoptProblem,
setIpoptProblemScaling,
setIntermediateCallback,
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
import qualified Numeric.AD.Internal.Identity as AD
import Ipopt.AnyRF
import Data.VectorSpace (VectorSpace,Scalar)
type IpNumber = (CDouble)
type IpIndex = (CInt)
type IpInt = (CInt)
type IpBool = (CInt)
type IpF = ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt))))))))
type IpGradF = ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt))))))))
type IpG = ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> (CInt -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt)))))))))
type IpJacG = ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> (CInt -> (CInt -> ((Ptr CInt) -> ((Ptr CInt) -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt))))))))))))
type IpH = ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> (CDouble -> (CInt -> ((Ptr CDouble) -> (CInt -> (CInt -> ((Ptr CInt) -> ((Ptr CInt) -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt)))))))))))))))
type IpIntermediateCB = ((FunPtr (CInt -> (CInt -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CInt -> ((Ptr ()) -> (IO CInt)))))))))))))))
type IntermediateCB = CInt
-> CInt
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> CInt
-> Ptr ()
-> IO IpBool
data ApplicationReturnStatus = SolveSucceeded
| SolvedToAcceptableLevel
| InfeasibleProblemDetected
| SearchDirectionBecomesTooSmall
| DivergingIterates
| UserRequestedStop
| FeasiblePointFound
| MaximumIterationsExceeded
| RestorationFailed
| ErrorInStepComputation
| MaximumCputimeExceeded
| NotEnoughDegreesOfFreedom
| InvalidProblemDefinition
| InvalidOption
| InvalidNumberDetected
| UnrecoverableException
| NonipoptExceptionThrown
| InsufficientMemory
| InternalError
deriving (Show)
instance Enum ApplicationReturnStatus where
fromEnum SolveSucceeded = 0
fromEnum SolvedToAcceptableLevel = 1
fromEnum InfeasibleProblemDetected = 2
fromEnum SearchDirectionBecomesTooSmall = 3
fromEnum DivergingIterates = 4
fromEnum UserRequestedStop = 5
fromEnum FeasiblePointFound = 6
fromEnum MaximumIterationsExceeded = (1)
fromEnum RestorationFailed = (2)
fromEnum ErrorInStepComputation = (3)
fromEnum MaximumCputimeExceeded = (4)
fromEnum NotEnoughDegreesOfFreedom = (10)
fromEnum InvalidProblemDefinition = (11)
fromEnum InvalidOption = (12)
fromEnum InvalidNumberDetected = (13)
fromEnum UnrecoverableException = (100)
fromEnum NonipoptExceptionThrown = (101)
fromEnum InsufficientMemory = (102)
fromEnum InternalError = (199)
toEnum 0 = SolveSucceeded
toEnum 1 = SolvedToAcceptableLevel
toEnum 2 = InfeasibleProblemDetected
toEnum 3 = SearchDirectionBecomesTooSmall
toEnum 4 = DivergingIterates
toEnum 5 = UserRequestedStop
toEnum 6 = FeasiblePointFound
toEnum (1) = MaximumIterationsExceeded
toEnum (2) = RestorationFailed
toEnum (3) = ErrorInStepComputation
toEnum (4) = MaximumCputimeExceeded
toEnum (10) = NotEnoughDegreesOfFreedom
toEnum (11) = InvalidProblemDefinition
toEnum (12) = InvalidOption
toEnum (13) = InvalidNumberDetected
toEnum (100) = UnrecoverableException
toEnum (101) = NonipoptExceptionThrown
toEnum (102) = InsufficientMemory
toEnum (199) = InternalError
toEnum unmatched = error ("ApplicationReturnStatus.toEnum: Cannot match " ++ show unmatched)
data AlgorithmMode = Regularmode
| Restorationphasemode
deriving (Show)
instance Enum AlgorithmMode where
fromEnum Regularmode = 0
fromEnum Restorationphasemode = 1
toEnum 0 = Regularmode
toEnum 1 = Restorationphasemode
toEnum unmatched = error ("AlgorithmMode.toEnum: Cannot match " ++ show unmatched)
newtype IpProblem = IpProblem (ForeignPtr (IpProblem))
withIpProblem :: IpProblem -> (Ptr IpProblem -> IO b) -> IO b
withIpProblem (IpProblem fptr) = withForeignPtr fptr
ipp x = withIpProblem x
type family UnFunPtr a
type instance UnFunPtr (FunPtr a) = a
data IpBoolT = IpTrue
| IpFalse
deriving (Eq,Ord,Show)
instance Enum IpBoolT where
fromEnum IpTrue = 1
fromEnum IpFalse = 0
toEnum 1 = IpTrue
toEnum 0 = IpFalse
toEnum unmatched = error ("IpBoolT.toEnum: Cannot match " ++ show unmatched)
ipTrue = toEnum (fromEnum IpTrue) :: IpBool
ipFalse = toEnum (fromEnum IpFalse) :: IpBool
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
foreign import ccall "wrapper" wrapIpIntermediateCB :: IntermediateCB -> IO IpIntermediateCB
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 :: Vec -> (Ptr CDouble -> IO r) -> IO r
vmUnsafeWith v f = VM.unsafeWith v (f . castPtr)
type Vec = VM.IOVector Double
_ = CDouble (5 :: Double) :: IpNumber
fromVec :: VG.Vector v Double => Vec -> IO (v Double)
fromVec mv = do
v <- VS.freeze mv
return (VG.convert 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
| lx <- VM.length xL,
lx == VM.length xU,
lg <- VM.length gL,
lg == VM.length gU = do
p <- createIpoptProblem3 lx xL xU lg gL gU nJac nHess 0 f g gradF jacG hess
p' <- newForeignPtr freeIpoptProblem (castPtr p)
return (IpProblem (castForeignPtr p'))
| otherwise = error "dimensions wrong!"
createIpoptProblem3 :: (Int) -> (Vec) -> (Vec) -> (Int) -> (Vec) -> (Vec) -> (Int) -> (Int) -> (Int) -> (IpF) -> (IpG) -> (IpGradF) -> (IpJacG) -> (IpH) -> IO ((Ptr IpProblem))
createIpoptProblem3 a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 =
let {a1' = fromIntegral a1} in
vmUnsafeWith a2 $ \a2' ->
vmUnsafeWith a3 $ \a3' ->
let {a4' = fromIntegral a4} in
vmUnsafeWith a5 $ \a5' ->
vmUnsafeWith a6 $ \a6' ->
let {a7' = fromIntegral a7} in
let {a8' = fromIntegral a8} in
let {a9' = fromIntegral a9} in
let {a10' = id a10} in
let {a11' = id a11} in
let {a12' = id a12} in
let {a13' = id a13} in
let {a14' = id a14} in
createIpoptProblem3'_ a1' a2' a3' a4' a5' a6' a7' a8' a9' a10' a11' a12' a13' a14' >>= \res ->
let {res' = id res} in
return (res')
addIpoptNumOption :: (IpProblem) -> (String) -> (Double) -> IO ((Bool))
addIpoptNumOption a1 a2 a3 =
ipp a1 $ \a1' ->
withCString a2 $ \a2' ->
let {a3' = realToFrac a3} in
addIpoptNumOption'_ a1' a2' a3' >>= \res ->
let {res' = toBool res} in
return (res')
addIpoptStrOption :: (IpProblem) -> (String) -> (String) -> IO ((Bool))
addIpoptStrOption a1 a2 a3 =
ipp a1 $ \a1' ->
withCString a2 $ \a2' ->
withCString a3 $ \a3' ->
addIpoptStrOption'_ a1' a2' a3' >>= \res ->
let {res' = toBool res} in
return (res')
addIpoptIntOption :: (IpProblem) -> (String) -> (Int) -> IO ((Bool))
addIpoptIntOption a1 a2 a3 =
ipp a1 $ \a1' ->
withCString a2 $ \a2' ->
let {a3' = fromIntegral a3} in
addIpoptIntOption'_ a1' a2' a3' >>= \res ->
let {res' = toBool res} in
return (res')
foreign import ccall unsafe "&FreeIpoptProblem"
freeIpoptProblem :: FunPtr (Ptr () -> IO ())
openIpoptOutputFile :: (IpProblem) -> (String) -> (Int) -> IO ((Bool))
openIpoptOutputFile a1 a2 a3 =
ipp a1 $ \a1' ->
withCString a2 $ \a2' ->
let {a3' = fromIntegral a3} in
openIpoptOutputFile'_ a1' a2' a3' >>= \res ->
let {res' = toBool res} in
return (res')
setIpoptProblemScaling :: (IpProblem) -> (Double) -> (Vec) -> (Vec) -> IO ((Bool))
setIpoptProblemScaling a1 a2 a3 a4 =
ipp a1 $ \a1' ->
let {a2' = realToFrac a2} in
vmUnsafeWith a3 $ \a3' ->
vmUnsafeWith a4 $ \a4' ->
setIpoptProblemScaling'_ a1' a2' a3' a4' >>= \res ->
let {res' = toBool res} in
return (res')
setIntermediateCallback1 :: (IpProblem) -> (IpIntermediateCB) -> IO ((Bool))
setIntermediateCallback1 a1 a2 =
ipp a1 $ \a1' ->
let {a2' = id a2} in
setIntermediateCallback1'_ a1' a2' >>= \res ->
let {res' = toBool res} in
return (res')
setIntermediateCallback pp cb = do
cb' <- wrapIpIntermediateCB cb
setIntermediateCallback1 pp cb'
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 :: VG.Vector v Double => IpProblem
-> Vec
-> IO (IpOptSolved v)
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
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'
ipoptSolve2 :: (IpProblem) -> (Vec) -> (Vec) -> (Vec) -> (Vec) -> (Vec) -> (Ptr ()) -> IO ((ApplicationReturnStatus), (Double))
ipoptSolve2 a1 a2 a3 a5 a6 a7 a8 =
ipp a1 $ \a1' ->
vmUnsafeWith a2 $ \a2' ->
vmUnsafeWith a3 $ \a3' ->
alloca $ \a4' ->
vmUnsafeWith a5 $ \a5' ->
vmUnsafeWith a6 $ \a6' ->
vmUnsafeWith a7 $ \a7' ->
let {a8' = id a8} in
ipoptSolve2'_ a1' a2' a3' a4' a5' a6' a7' a8' >>= \res ->
let {res' = cToEnum res} in
peekFloatConv a4'>>= \a4'' ->
return (res', a4'')
createIpoptProblemAD
:: Vec
-> Vec
-> Vec
-> Vec
-> (forall a. AnyRFCxt a => V.Vector a -> a)
-> (forall a. AnyRFCxt a => V.Vector a -> V.Vector a)
-> 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
-> Int
-> (forall a. AnyRFCxt a => V.Vector a -> a)
-> (forall a. AnyRFCxt a => V.Vector a -> V.Vector a)
-> 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)
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))
denseIJh n m iRow jCol = do
i <- newIORef 0
forM_ [0 .. fromIntegral n1] $ \ row ->
forM_ [ 0 .. row ] $ \col -> do
ii <- readIORef i
VM.write iRow ii row
VM.write jCol ii col
writeIORef i (ii+1)
foreign import ccall safe "Ipopt/Raw.chs.h CreateIpoptProblem"
createIpoptProblem3'_ :: (CInt -> ((Ptr CDouble) -> ((Ptr CDouble) -> (CInt -> ((Ptr CDouble) -> ((Ptr CDouble) -> (CInt -> (CInt -> (CInt -> ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt))))))) -> ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> (CInt -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt)))))))) -> ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt))))))) -> ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> (CInt -> (CInt -> ((Ptr CInt) -> ((Ptr CInt) -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt))))))))))) -> ((FunPtr (CInt -> ((Ptr CDouble) -> (CInt -> (CDouble -> (CInt -> ((Ptr CDouble) -> (CInt -> (CInt -> ((Ptr CInt) -> ((Ptr CInt) -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt)))))))))))))) -> (IO (Ptr (IpProblem)))))))))))))))))
foreign import ccall safe "Ipopt/Raw.chs.h AddIpoptNumOption"
addIpoptNumOption'_ :: ((Ptr (IpProblem)) -> ((Ptr CChar) -> (CDouble -> (IO CInt))))
foreign import ccall safe "Ipopt/Raw.chs.h AddIpoptStrOption"
addIpoptStrOption'_ :: ((Ptr (IpProblem)) -> ((Ptr CChar) -> ((Ptr CChar) -> (IO CInt))))
foreign import ccall safe "Ipopt/Raw.chs.h AddIpoptIntOption"
addIpoptIntOption'_ :: ((Ptr (IpProblem)) -> ((Ptr CChar) -> (CInt -> (IO CInt))))
foreign import ccall safe "Ipopt/Raw.chs.h OpenIpoptOutputFile"
openIpoptOutputFile'_ :: ((Ptr (IpProblem)) -> ((Ptr CChar) -> (CInt -> (IO CInt))))
foreign import ccall safe "Ipopt/Raw.chs.h SetIpoptProblemScaling"
setIpoptProblemScaling'_ :: ((Ptr (IpProblem)) -> (CDouble -> ((Ptr CDouble) -> ((Ptr CDouble) -> (IO CInt)))))
foreign import ccall safe "Ipopt/Raw.chs.h SetIntermediateCallback"
setIntermediateCallback1'_ :: ((Ptr (IpProblem)) -> ((FunPtr (CInt -> (CInt -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CDouble -> (CInt -> ((Ptr ()) -> (IO CInt)))))))))))))) -> (IO CInt)))
foreign import ccall safe "Ipopt/Raw.chs.h IpoptSolve"
ipoptSolve2'_ :: ((Ptr (IpProblem)) -> ((Ptr CDouble) -> ((Ptr CDouble) -> ((Ptr CDouble) -> ((Ptr CDouble) -> ((Ptr CDouble) -> ((Ptr CDouble) -> ((Ptr ()) -> (IO CInt)))))))))