module Ipopt.Raw (
createIpoptProblemAD,
createIpoptProblemADSparse,
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 qualified Foreign.C.String as C2HSImp
import qualified Foreign.C.Types as C2HSImp
import qualified Foreign.ForeignPtr as C2HSImp
import qualified Foreign.Marshal.Utils as C2HSImp
import qualified Foreign.Ptr as C2HSImp
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 qualified Numeric.AD.Rank1.Sparse as Sparse
import qualified Numeric.AD.Internal.Sparse as Sparse
import Numeric.AD.Rank1.Sparse (Sparse)
import Ipopt.AnyRF
import Data.VectorSpace (VectorSpace,Scalar)
type IpNumber = (C2HSImp.CDouble)
type IpIndex = (C2HSImp.CInt)
type IpInt = (C2HSImp.CInt)
type IpBool = (C2HSImp.CInt)
type IpF = ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt))))))))
type IpGradF = ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt))))))))
type IpG = ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt)))))))))
type IpJacG = ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CInt -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CInt) -> ((C2HSImp.Ptr C2HSImp.CInt) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt))))))))))))
type IpH = ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CDouble -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CInt) -> ((C2HSImp.Ptr C2HSImp.CInt) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt)))))))))))))))
type IpIntermediateCB = ((C2HSImp.FunPtr (C2HSImp.CInt -> (C2HSImp.CInt -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CInt -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt)))))))))))))))
type IntermediateCB = CInt
-> CInt
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> CInt
-> Ptr ()
-> IO IpBool
data ApplicationReturnStatus = InternalError
| InsufficientMemory
| NonipoptExceptionThrown
| UnrecoverableException
| InvalidNumberDetected
| InvalidOption
| InvalidProblemDefinition
| NotEnoughDegreesOfFreedom
| MaximumCputimeExceeded
| ErrorInStepComputation
| RestorationFailed
| MaximumIterationsExceeded
| SolveSucceeded
| SolvedToAcceptableLevel
| InfeasibleProblemDetected
| SearchDirectionBecomesTooSmall
| DivergingIterates
| UserRequestedStop
| FeasiblePointFound
deriving (Show)
instance Enum ApplicationReturnStatus where
succ InternalError = InsufficientMemory
succ InsufficientMemory = NonipoptExceptionThrown
succ NonipoptExceptionThrown = UnrecoverableException
succ UnrecoverableException = InvalidNumberDetected
succ InvalidNumberDetected = InvalidOption
succ InvalidOption = InvalidProblemDefinition
succ InvalidProblemDefinition = NotEnoughDegreesOfFreedom
succ NotEnoughDegreesOfFreedom = MaximumCputimeExceeded
succ MaximumCputimeExceeded = ErrorInStepComputation
succ ErrorInStepComputation = RestorationFailed
succ RestorationFailed = MaximumIterationsExceeded
succ MaximumIterationsExceeded = SolveSucceeded
succ SolveSucceeded = SolvedToAcceptableLevel
succ SolvedToAcceptableLevel = InfeasibleProblemDetected
succ InfeasibleProblemDetected = SearchDirectionBecomesTooSmall
succ SearchDirectionBecomesTooSmall = DivergingIterates
succ DivergingIterates = UserRequestedStop
succ UserRequestedStop = FeasiblePointFound
succ FeasiblePointFound = error "ApplicationReturnStatus.succ: FeasiblePointFound has no successor"
pred InsufficientMemory = InternalError
pred NonipoptExceptionThrown = InsufficientMemory
pred UnrecoverableException = NonipoptExceptionThrown
pred InvalidNumberDetected = UnrecoverableException
pred InvalidOption = InvalidNumberDetected
pred InvalidProblemDefinition = InvalidOption
pred NotEnoughDegreesOfFreedom = InvalidProblemDefinition
pred MaximumCputimeExceeded = NotEnoughDegreesOfFreedom
pred ErrorInStepComputation = MaximumCputimeExceeded
pred RestorationFailed = ErrorInStepComputation
pred MaximumIterationsExceeded = RestorationFailed
pred SolveSucceeded = MaximumIterationsExceeded
pred SolvedToAcceptableLevel = SolveSucceeded
pred InfeasibleProblemDetected = SolvedToAcceptableLevel
pred SearchDirectionBecomesTooSmall = InfeasibleProblemDetected
pred DivergingIterates = SearchDirectionBecomesTooSmall
pred UserRequestedStop = DivergingIterates
pred FeasiblePointFound = UserRequestedStop
pred InternalError = error "ApplicationReturnStatus.pred: InternalError has no predecessor"
enumFromTo from to = go from
where
end = fromEnum to
go v = case compare (fromEnum v) end of
LT -> v : go (succ v)
EQ -> [v]
GT -> []
enumFrom from = enumFromTo from FeasiblePointFound
fromEnum InternalError = (199)
fromEnum InsufficientMemory = (102)
fromEnum NonipoptExceptionThrown = (101)
fromEnum UnrecoverableException = (100)
fromEnum InvalidNumberDetected = (13)
fromEnum InvalidOption = (12)
fromEnum InvalidProblemDefinition = (11)
fromEnum NotEnoughDegreesOfFreedom = (10)
fromEnum MaximumCputimeExceeded = (4)
fromEnum ErrorInStepComputation = (3)
fromEnum RestorationFailed = (2)
fromEnum MaximumIterationsExceeded = (1)
fromEnum SolveSucceeded = 0
fromEnum SolvedToAcceptableLevel = 1
fromEnum InfeasibleProblemDetected = 2
fromEnum SearchDirectionBecomesTooSmall = 3
fromEnum DivergingIterates = 4
fromEnum UserRequestedStop = 5
fromEnum FeasiblePointFound = 6
toEnum (199) = InternalError
toEnum (102) = InsufficientMemory
toEnum (101) = NonipoptExceptionThrown
toEnum (100) = UnrecoverableException
toEnum (13) = InvalidNumberDetected
toEnum (12) = InvalidOption
toEnum (11) = InvalidProblemDefinition
toEnum (10) = NotEnoughDegreesOfFreedom
toEnum (4) = MaximumCputimeExceeded
toEnum (3) = ErrorInStepComputation
toEnum (2) = RestorationFailed
toEnum (1) = MaximumIterationsExceeded
toEnum 0 = SolveSucceeded
toEnum 1 = SolvedToAcceptableLevel
toEnum 2 = InfeasibleProblemDetected
toEnum 3 = SearchDirectionBecomesTooSmall
toEnum 4 = DivergingIterates
toEnum 5 = UserRequestedStop
toEnum 6 = FeasiblePointFound
toEnum unmatched = error ("ApplicationReturnStatus.toEnum: Cannot match " ++ show unmatched)
data AlgorithmMode = Regularmode
| Restorationphasemode
deriving (Show)
instance Enum AlgorithmMode where
succ Regularmode = Restorationphasemode
succ Restorationphasemode = error "AlgorithmMode.succ: Restorationphasemode has no successor"
pred Restorationphasemode = Regularmode
pred Regularmode = error "AlgorithmMode.pred: Regularmode has no predecessor"
enumFromTo from to = go from
where
end = fromEnum to
go v = case compare (fromEnum v) end of
LT -> v : go (succ v)
EQ -> [v]
GT -> []
enumFrom from = enumFromTo from Restorationphasemode
fromEnum Regularmode = 0
fromEnum Restorationphasemode = 1
toEnum 0 = Regularmode
toEnum 1 = Restorationphasemode
toEnum unmatched = error ("AlgorithmMode.toEnum: Cannot match " ++ show unmatched)
newtype IpProblem = IpProblem (C2HSImp.ForeignPtr (IpProblem))
withIpProblem :: IpProblem -> (C2HSImp.Ptr IpProblem -> IO b) -> IO b
withIpProblem (IpProblem fptr) = C2HSImp.withForeignPtr fptr
ipp x = withIpProblem x
type family UnFunPtr a
type instance UnFunPtr (FunPtr a) = a
data IpBoolT = IpFalse
| IpTrue
deriving (Eq,Ord,Show)
instance Enum IpBoolT where
succ IpFalse = IpTrue
succ IpTrue = error "IpBoolT.succ: IpTrue has no successor"
pred IpTrue = IpFalse
pred IpFalse = error "IpBoolT.pred: IpFalse has no predecessor"
enumFromTo from to = go from
where
end = fromEnum to
go v = case compare (fromEnum v) end of
LT -> v : go (succ v)
EQ -> [v]
GT -> []
enumFrom from = enumFromTo from IpTrue
fromEnum IpFalse = 0
fromEnum IpTrue = 1
toEnum 0 = IpFalse
toEnum 1 = IpTrue
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' ->
C2HSImp.withCString a2 $ \a2' ->
let {a3' = realToFrac a3} in
addIpoptNumOption'_ a1' a2' a3' >>= \res ->
let {res' = C2HSImp.toBool res} in
return (res')
addIpoptStrOption :: (IpProblem) -> (String) -> (String) -> IO ((Bool))
addIpoptStrOption a1 a2 a3 =
ipp a1 $ \a1' ->
C2HSImp.withCString a2 $ \a2' ->
C2HSImp.withCString a3 $ \a3' ->
addIpoptStrOption'_ a1' a2' a3' >>= \res ->
let {res' = C2HSImp.toBool res} in
return (res')
addIpoptIntOption :: (IpProblem) -> (String) -> (Int) -> IO ((Bool))
addIpoptIntOption a1 a2 a3 =
ipp a1 $ \a1' ->
C2HSImp.withCString a2 $ \a2' ->
let {a3' = fromIntegral a3} in
addIpoptIntOption'_ a1' a2' a3' >>= \res ->
let {res' = C2HSImp.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' ->
C2HSImp.withCString a2 $ \a2' ->
let {a3' = fromIntegral a3} in
openIpoptOutputFile'_ a1' a2' a3' >>= \res ->
let {res' = C2HSImp.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' = C2HSImp.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' = C2HSImp.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
createIpoptProblemADSparse
:: Vec
-> Vec
-> Vec
-> Vec
-> (V.Vector (Sparse CDouble) -> Sparse CDouble)
-> (V.Vector (Sparse CDouble) -> V.Vector (Sparse CDouble))
-> IO IpProblem
createIpoptProblemADSparse 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) <- mkFsSparse 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)
sparsePrimal :: Num a => Sparse a -> a
sparsePrimal (Sparse.Sparse a _) = a
sparsePrimal Sparse.Zero = 0
mkFsSparse ::
Int
-> Int
-> (V.Vector (Sparse CDouble) -> Sparse CDouble)
-> (V.Vector (Sparse CDouble) -> V.Vector (Sparse CDouble))
-> IO (IpF, IpGradF, IpG, IpJacG, IpH)
mkFsSparse n m f g = do
ipF <- wrapIpF $ \x -> do
x <- VG.convert `fmap` VS.unsafeFreeze x
return $ sparsePrimal $ f (Sparse.vars x)
ipGradF <- wrapIpGradF $ \x -> do
x <- VG.convert `fmap` VS.unsafeFreeze x
VS.unsafeThaw $ VG.convert (Sparse.grad f x)
ipG <- wrapIpG $ \x -> do
x <- VG.convert `fmap` VS.unsafeFreeze x
VS.unsafeThaw $ VG.convert $ V.map sparsePrimal $ g $ Sparse.vars 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 $ Sparse.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 $ Sparse.hessian f x
gj = V.zipWith (\l v -> V.map (l*) v) lambda (V.map tri (Sparse.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'_ :: (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CInt -> (C2HSImp.CInt -> ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt))))))) -> ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt)))))))) -> ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt))))))) -> ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CInt -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CInt) -> ((C2HSImp.Ptr C2HSImp.CInt) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt))))))))))) -> ((C2HSImp.FunPtr (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CDouble -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (C2HSImp.CInt -> (C2HSImp.CInt -> ((C2HSImp.Ptr C2HSImp.CInt) -> ((C2HSImp.Ptr C2HSImp.CInt) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt)))))))))))))) -> (IO (C2HSImp.Ptr (IpProblem)))))))))))))))))
foreign import ccall safe "Ipopt/Raw.chs.h AddIpoptNumOption"
addIpoptNumOption'_ :: ((C2HSImp.Ptr (IpProblem)) -> ((C2HSImp.Ptr C2HSImp.CChar) -> (C2HSImp.CDouble -> (IO C2HSImp.CInt))))
foreign import ccall safe "Ipopt/Raw.chs.h AddIpoptStrOption"
addIpoptStrOption'_ :: ((C2HSImp.Ptr (IpProblem)) -> ((C2HSImp.Ptr C2HSImp.CChar) -> ((C2HSImp.Ptr C2HSImp.CChar) -> (IO C2HSImp.CInt))))
foreign import ccall safe "Ipopt/Raw.chs.h AddIpoptIntOption"
addIpoptIntOption'_ :: ((C2HSImp.Ptr (IpProblem)) -> ((C2HSImp.Ptr C2HSImp.CChar) -> (C2HSImp.CInt -> (IO C2HSImp.CInt))))
foreign import ccall safe "Ipopt/Raw.chs.h OpenIpoptOutputFile"
openIpoptOutputFile'_ :: ((C2HSImp.Ptr (IpProblem)) -> ((C2HSImp.Ptr C2HSImp.CChar) -> (C2HSImp.CInt -> (IO C2HSImp.CInt))))
foreign import ccall safe "Ipopt/Raw.chs.h SetIpoptProblemScaling"
setIpoptProblemScaling'_ :: ((C2HSImp.Ptr (IpProblem)) -> (C2HSImp.CDouble -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> (IO C2HSImp.CInt)))))
foreign import ccall safe "Ipopt/Raw.chs.h SetIntermediateCallback"
setIntermediateCallback1'_ :: ((C2HSImp.Ptr (IpProblem)) -> ((C2HSImp.FunPtr (C2HSImp.CInt -> (C2HSImp.CInt -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CDouble -> (C2HSImp.CInt -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt)))))))))))))) -> (IO C2HSImp.CInt)))
foreign import ccall safe "Ipopt/Raw.chs.h IpoptSolve"
ipoptSolve2'_ :: ((C2HSImp.Ptr (IpProblem)) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr C2HSImp.CDouble) -> ((C2HSImp.Ptr ()) -> (IO C2HSImp.CInt)))))))))