Ipopt.AnyRF.hs

{-# LANGUAGE ConstraintKinds, RankNTypes, TypeFamilies, FlexibleInstances #-}
{- | 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.*!)

NLOpt.hs

{-# LANGUAGE TypeFamilies,DeriveDataTypeable, RankNTypes, ConstraintKinds, FlexibleContexts #-}
module NLOpt where

import Ipopt.AnyRF

import Foreign.C.Types
import Foreign.Ptr
import Foreign.Marshal
import Foreign.Storable
import Data.Int
import C2HS

import qualified Data.Vector.Storable.Mutable as VM
import qualified Data.Vector as V
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Generic as VG

import Control.Exception
import Data.Typeable
import Numeric.AD (grad, jacobian, hessian)
import Control.Monad

#include "nlopt.h"

-- * converting haskell functions

{- $conventions

NLOpt has three different types for functions 'FunPtrFunc' 'FunPtrMFunc' and 'FunPtrPrecond'.

[@AD@] means derivatives are calculated, and the haskell function does no IO.

[@G@] mean you are providing the derivatives

[@M@] means m functions are calculated at a time

-}

-- | an exact hessian calculated with AD. See 'toPrecondG'
-- XXX BFGS approx could also be done...
toPrecondAD :: (forall a. AnyRFCxt a => V.Vector a -> a) -> Precond
toPrecondAD f n xs vs r _usrData =
    toPrecondG (\x v -> return (hessian f x `mXv` v)) n xs vs r _usrData

-- | see <http://ab-initio.mit.edu/wiki/index.php/NLopt_Reference#Preconditioning_with_approximate_Hessians>
-- only applies to 'NLOPT_LD_CCSAQ'
toPrecondG :: (VG.Vector vec CDouble, v ~ vec CDouble)
    => (v -> v -> IO v) -- ^ given @x,v@ calculate  @H(x)v@ where H the hessian
    -> Precond
toPrecondG f n xs vs r _usrData = do
    xs' <- ptrToV n xs
    vs' <- ptrToV n vs
    copyInto n r =<< f xs' vs'

toFuncM :: VG.Vector v CDouble => (v CDouble -> IO (v CDouble) ) -> MFunc
toFuncM f m r n xs _g _usrData = do
    copyInto m r =<< f =<< ptrToV n xs

-- | @n@ and @m@ type variables indicate the vector size as
-- number of inputs and number of outputs respectively
toFuncMG :: (VG.Vector n CDouble, VG.Vector n (m CDouble), n ~ m)
    => (n CDouble -> IO (m CDouble) ) -- @f@
    -> (n CDouble -> IO (n (m CDouble))) -- @grad f@
    -> MFunc
toFuncMG f g m r n xs g' _usrData = do
    xs' <- ptrToV n xs
    copyInto m r =<< f xs'
    -- probably should do this better...
    copyInto (n*m) g' . VG.concat . VG.toList =<< g xs'

toFuncMAD :: (forall a. AnyRFCxt a => V.Vector a -> V.Vector a) -> MFunc
toFuncMAD f m r n xs g _usrData = do
    xs' <- ptrToV n xs
    copyInto m r (f xs')
    -- grad[i*n + j] should be satisfied...
    copyInto (n*m) g (V.concat (V.toList (jacobian f xs')))

toFunc :: (VG.Vector v CDouble) => (v CDouble -> IO CDouble) -> Func
toFunc f n xs _g _usrData = f =<< ptrToV n xs

-- | where the gradient happens via AD
toFuncAD :: (forall a. AnyRFCxt a => V.Vector a -> a) -> Func
toFuncAD f n xs g _usrData = do
    xs' <- ptrToV n xs
    copyInto n g (grad f xs')
    return (f xs')

toFuncG :: (VG.Vector v CDouble) => (v CDouble -> IO CDouble) -- ^ @f@
    -> (v CDouble -> IO (v CDouble)) -- ^ @grad(f)@
    -> Func
toFuncG f g n xs g' _usrData = do
    xs' <- ptrToV n xs
    copyInto n g' =<< g xs'
    f xs'

type Func = UnFunPtr FunPtrFunc
type MFunc = UnFunPtr FunPtrMFunc
type Precond = UnFunPtr FunPtrPrecond

type FunPtrFunc    = {# type nlopt_func #}
type FunPtrMFunc   = {# type nlopt_mfunc #}
type FunPtrPrecond = {# type nlopt_precond #}

{#pointer *nlopt_opt as NLOpt foreign newtype #}
{#enum nlopt_algorithm as ^ {} deriving (Show) #}

-- | negative (above NLOPT_SUCCESS) values of these are thrown as exceptions. The positive ones are
-- return values.
{#enum nlopt_result as ^ {} deriving (Show,Typeable) #}

instance Exception NloptResult

checkEC :: CInt -> IO NloptResult
checkEC n | n < 0 = throwIO e
          | otherwise = return e
    where e = fromCInt n :: NloptResult

{#fun nlopt_srand as ^ { `Int' } -> `()' #}
{#fun nlopt_srand_time as ^ { } -> `()' #}
{#fun nlopt_version as ^ {
        alloca- `Int' peekInt*,
        alloca- `Int' peekInt*,
        alloca- `Int' peekInt*} -> `()' #}
{#fun nlopt_create as ^ { toCInt `NloptAlgorithm', `Int' } -> `NLOpt' ptrToNLOpt * #}
{#fun nlopt_copy as ^ { withNLOpt_* `NLOpt' } -> `NLOpt' ptrToNLOpt * #}

-- | should not need to be called manually
{#fun nlopt_destroy as ^ { id `Ptr ()' } -> `()' #}

{#fun nlopt_optimize as ^ {
    withNLOpt_* `NLOpt', vmUnsafeWith* `Vec', alloca- `CDouble' peek* }
        -> `NloptResult' checkEC* #}

{#fun nlopt_set_min_objective as ^ {
    withNLOpt_* `NLOpt', withFunc* `Func', withNull- `()' }
        -> `NloptResult' checkEC* #}
{#fun nlopt_set_max_objective as ^ {
    withNLOpt_* `NLOpt', withFunc* `Func', withNull- `()' }
        -> `NloptResult' checkEC* #}

{#fun nlopt_set_precond_min_objective as ^ {
    withNLOpt_* `NLOpt',
    withFunc* `Func',
    withPrecond* `Precond',
    withNull- `()' }
    -> `NloptResult' checkEC*  #}

{#fun nlopt_set_precond_max_objective as ^ {
    withNLOpt_* `NLOpt',
    withFunc* `Func',
    withPrecond* `Precond',
    withNull- `()' } -> `NloptResult' checkEC* #}

{#fun nlopt_get_algorithm as ^ { withNLOpt_* `NLOpt' } -> `NloptAlgorithm' fromCInt #}
{#fun nlopt_get_dimension as ^ { withNLOpt_* `NLOpt' } -> `Int' #}

-- * constraints
{#fun nlopt_set_lower_bounds as ^ { withNLOpt_* `NLOpt', vmUnsafeWith* `Vec' }
    -> `NloptResult' checkEC* #}
{#fun nlopt_set_upper_bounds as ^ { withNLOpt_* `NLOpt', vmUnsafeWith* `Vec' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_get_lower_bounds as ^ { withNLOpt_* `NLOpt', vmUnsafeWith* `Vec' }
    -> `NloptResult' checkEC* #}
{#fun nlopt_get_upper_bounds as ^ { withNLOpt_* `NLOpt', vmUnsafeWith* `Vec' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_set_lower_bounds1 as ^ { withNLOpt_* `NLOpt', `Double' }
    -> `NloptResult' checkEC* #}
{#fun nlopt_set_upper_bounds1 as ^ { withNLOpt_* `NLOpt', `Double' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_remove_inequality_constraints as ^ { withNLOpt_* `NLOpt' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_add_inequality_constraint as ^
    { withNLOpt_* `NLOpt',
      withFunc* `Func',
      withNull- `()',
      `Double' } -> `NloptResult' checkEC* #}

{#fun nlopt_add_precond_inequality_constraint as ^
    { withNLOpt_* `NLOpt',
      withFunc* `Func',
      withPrecond* `Precond',
      withNull- `()',
      `Double' } -> `NloptResult' checkEC* #}

{#fun nlopt_add_inequality_mconstraint as ^
    { withNLOpt_* `NLOpt',
      `Int',
      withMFunc* `MFunc',
      withNull- `()',
      vmUnsafeWith* `Vec' } -> `NloptResult' checkEC* #}

{#fun nlopt_remove_equality_constraints as ^ { withNLOpt_* `NLOpt' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_add_equality_constraint as ^
    { withNLOpt_* `NLOpt',
      withFunc* `Func',
      withNull- `()',
      `Double' } -> `NloptResult' checkEC* #}

{#fun nlopt_add_precond_equality_constraint as ^
    { withNLOpt_* `NLOpt',
      withFunc* `Func',
      withPrecond* `Precond',
      withNull- `()',
      `Double' } -> `NloptResult' checkEC* #}

{#fun nlopt_add_equality_mconstraint as ^
    { withNLOpt_* `NLOpt',
      `Int',
      withMFunc* `MFunc',
      withNull- `()',
      vmUnsafeWith* `Vec' } -> `NloptResult' checkEC* #}

-- * stopping criteria
{#fun nlopt_set_stopval as ^ { withNLOpt_* `NLOpt', `Double' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_get_stopval as ^ { withNLOpt_* `NLOpt' } -> `Double' #}

{#fun nlopt_set_ftol_rel as ^ { withNLOpt_* `NLOpt', `Double' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_get_ftol_rel as ^ { withNLOpt_* `NLOpt' } -> `Double' #}

{#fun nlopt_set_ftol_abs as ^ { withNLOpt_* `NLOpt', `Double' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_get_ftol_abs as ^ { withNLOpt_* `NLOpt' } -> `Double' #}

{#fun nlopt_set_xtol_rel as ^ { withNLOpt_* `NLOpt', `Double' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_get_xtol_rel as ^ { withNLOpt_* `NLOpt' } -> `Double' #}

{#fun nlopt_set_xtol_abs1 as ^ { withNLOpt_* `NLOpt', `Double' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_set_xtol_abs as ^ { withNLOpt_* `NLOpt', vmUnsafeWith* `Vec' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_get_xtol_abs as ^ { withNLOpt_* `NLOpt', vmUnsafeWith* `Vec' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_set_maxeval as ^ { withNLOpt_* `NLOpt', `Int' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_get_maxeval as ^ { withNLOpt_* `NLOpt' } -> `Int'  #}

{#fun nlopt_set_maxtime as ^ { withNLOpt_* `NLOpt', `Double' }
    -> `NloptResult' checkEC* #}

{#fun nlopt_get_maxtime as ^ { withNLOpt_* `NLOpt' } -> `Double' #}

{#fun nlopt_force_stop as ^ { withNLOpt_* `NLOpt' } -> `NloptResult' checkEC* #}
{#fun nlopt_set_force_stop as ^ { withNLOpt_* `NLOpt', `Int' } -> `NloptResult' checkEC* #}
{#fun nlopt_get_force_stop as ^ { withNLOpt_* `NLOpt' } -> `Int' #}


-- * more algorithm-specific parameters
{#fun nlopt_set_local_optimizer as ^
    { withNLOpt_* `NLOpt', withNLOpt_* `NLOpt' } -> `NloptResult' checkEC* #}

{#fun nlopt_set_population as ^
    { withNLOpt_* `NLOpt', `Int' } -> `NloptResult' checkEC* #}

{#fun nlopt_get_population as ^ { withNLOpt_* `NLOpt' } -> `Int' #}

{#fun nlopt_set_vector_storage as ^
    { withNLOpt_* `NLOpt', `Int' } -> `NloptResult' checkEC* #}

{#fun nlopt_get_vector_storage as ^ { withNLOpt_* `NLOpt' } -> `Int' #}

{#fun nlopt_set_initial_step as ^ { withNLOpt_* `NLOpt',
    vmUnsafeWith* `Vec' } -> `NloptResult' checkEC* #}

{#fun nlopt_get_initial_step as ^ { withNLOpt_* `NLOpt',
    vmUnsafeWith* `Vec',
    vmUnsafeWith* `Vec' } -> `NloptResult' checkEC* #}

{#fun nlopt_set_initial_step1 as ^ { withNLOpt_* `NLOpt',
    `Double' } -> `NloptResult' checkEC* #}

-- * utils for ffi

-- $note
-- much is copied from Ipopt.Raw

withNLOpt_ x f = withNLOpt x (f . castPtr)

vmUnsafeWith = VM.unsafeWith

ptrToVS n p = do
    fp <- newForeignPtr_ p
    return (VM.unsafeFromForeignPtr0 fp (fromIntegral n))

ptrToV n p = fmap V.convert $ VS.unsafeFreeze =<< ptrToVS n p

copyInto n dest src = do
    to <- ptrToVS n dest
    VM.copy to =<< VS.unsafeThaw (V.convert src)

type family UnFunPtr a
type instance UnFunPtr (FunPtr a) = a
type Vec = VM.IOVector CDouble

toCInt x = fromIntegral (fromEnum x)
fromCInt x = toEnum (fromIntegral x)
peekInt ptr = fmap fromIntegral (peek ptr)


ptrToNLOpt :: Ptr () -> IO NLOpt
ptrToNLOpt p = do
    fp <- newForeignPtr nloptDestroyFP p
    return (NLOpt (castForeignPtr fp))

foreign import ccall "wrapper" mkNloptFinalizer :: (Ptr () -> IO ()) -> IO (FunPtr (Ptr () -> IO ()))

foreign import ccall unsafe "& nlopt_destroy" nloptDestroyFP :: FunPtr (Ptr () -> IO ())

foreign import ccall "wrapper" mkFunc :: Func -> IO FunPtrFunc
foreign import ccall "wrapper" mkPrecond :: Precond -> IO FunPtrPrecond
foreign import ccall "wrapper" mkMFunc :: MFunc -> IO FunPtrMFunc

withFunc :: Func -> (FunPtrFunc -> IO b) -> IO b
withFunc f g = do
    f' <- mkFunc f
    r <- g f'
    freeHaskellFunPtr f'
    return r

withPrecond :: Precond -> (FunPtrPrecond -> IO b) -> IO b
withPrecond f g = do
    f' <- mkPrecond f
    r <- g f'
    freeHaskellFunPtr f'
    return r

withMFunc :: MFunc -> (FunPtrMFunc -> IO b) -> IO b
withMFunc f g = do
    f' <- mkMFunc f
    r <- g f'
    freeHaskellFunPtr f'
    return r

withNull f = f nullPtr

-- | naive matrix × vector
mXv :: Num a => V.Vector (V.Vector a) -> V.Vector a -> V.Vector a
mXv m v = V.map (V.sum . V.zipWith (*) v) m

testnlopt.hs

import NLOpt
import qualified Data.Vector as V
import Data.Vector ((!))
import qualified Data.Vector.Storable as VS
main = do
    p <- nloptCreate NLOPT_LN_COBYLA 4
    nloptSetMinObjective p $ toFuncAD $ \x -> x!0 * x!3 * (x!0 + x!1 + x!2) + x!2

    let g1L = toFuncAD $ \x -> 2.0e19 - V.product x
        g1U = toFuncAD $ \x -> V.product x - 25
        g2 = toFuncAD $ \x -> V.sum (V.map (^2) x) - 40
    nloptAddInequalityConstraint p g1L 1e-2
    nloptAddInequalityConstraint p g1U 1e-2
    nloptAddEqualityConstraint p g2 1e-3

    x <- VS.thaw (VS.fromList [1,5,5,1])
    print =<< nloptOptimize p x
    x <- VS.freeze x
    print x

fixed by having copyTo check for a nullPtr