Monad.hs

{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{- |
The monadic interface to CLP allows to optimize
with respect to multiple objectives, successively.
-}
module Numeric.HiGHS.LP.Monad (
   -- * simple solver with warm restart
   T,
   run,
   solve,
   Direction(..),
   Method, Priv.simplex, Priv.choose, Priv.ipm,
   LPEnum.ModelStatus,
   Result,
   -- * solve with extra queries on the result
   solveWith,
   Query, Priv.getObjectiveValue,
   Priv.getOptimalVector, Priv.getSolutionVectors,
   Priv.getBasisStatus, Highs.BasisStatus,
   ) where

import qualified Numeric.HiGHS.LP.Enumeration as LPEnum
import qualified Numeric.HiGHS.LP.FFI as Highs
import qualified Numeric.HiGHS.LP.Debug as Debug
import qualified Numeric.HiGHS.LP.Private as Priv
import Numeric.HiGHS.LP.FFI (Highs)
import Numeric.HiGHS.LP.Private
         (Method, Result, Query, checkStatus, runContT, withBuffer,
          storeBounds, prepareRowBoundsArrays, prepareColumnBoundsArrays,
          storeConstraints, prepareConstraints,
          setMethod, objectiveSense, examineStatus)

import Numeric.LinearProgramming.Common
         (Bounds, Constraints, Direction(..), Objective)

import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape

import qualified Control.Monad.Trans.Cont as MC
import qualified Control.Monad.Trans.Reader as MR
import Control.Monad.IO.Class (liftIO)
import Control.Monad (when)
import Control.Exception (bracket)

import System.IO.Unsafe (unsafePerformIO)

import Foreign.Ptr (Ptr)


{- $setup
>>> :set -XTypeFamilies
>>> :set -XTypeOperators
>>> import qualified Numeric.HiGHS.LP.Monad as LP
>>> import qualified Numeric.HiGHS.LP as CLP
>>> import Test.Numeric.HiGHS.LP.Utility (traverse_Lag, traverseLag)
>>> import Test.Numeric.HiGHS.LP (TripletShape, tripletShape, forAllMethod)
>>> import Numeric.HiGHS.LP (Direction, (.*), (<=.))
>>>
>>> import qualified Numeric.LinearProgramming.Monad as LPMonad
>>> import qualified Numeric.LinearProgramming.Test as TestLP
>>> import Numeric.LinearProgramming.Common (Bounds, Objective)
>>>
>>> import qualified Data.Array.Comfort.Storable as Array
>>> import qualified Data.Array.Comfort.Shape as Shape
>>> import qualified Data.NonEmpty as NonEmpty
>>> import Data.Array.Comfort.Storable (Array)
>>> import Data.Traversable (Traversable)
>>> import Data.Foldable (Foldable)
>>>
>>> import qualified Control.Monad.Trans.Except as ME
>>>
>>> import qualified Data.List.HT as ListHT
>>> import Data.Tuple.HT (mapSnd)
>>>
>>> import Foreign.Storable (Storable)
>>>
>>> import qualified Test.QuickCheck as QC
>>>
>>>
>>> type Constraints ix = CLP.Constraints Double ix
>>>
>>>
>>> approxSuccession ::
>>>    (Shape.C sh, Show sh, Show a, Ord a, Num a, Storable a) =>
>>>    a ->
>>>    Either CLP.ModelStatus (NonEmpty.T [] (a, Array sh a)) ->
>>>    Either CLP.ModelStatus (NonEmpty.T [] (a, Array sh a)) ->
>>>    QC.Property
>>> approxSuccession tol x y =
>>>    QC.counterexample (show x) $
>>>    QC.counterexample (show y) $
>>>    case (x,y) of
>>>       (Left sx, Left sy) -> sx==sy
>>>       (Right (NonEmpty.Cons xh xs), Right (NonEmpty.Cons yh ys)) ->
>>>          let equalSol (optX, _) (optY, _) = TestLP.approxReal tol optX optY
>>>          in equalSol xh yh  &&  ListHT.equalWith equalSol xs ys
>>>       _ -> False
>>>
>>>
>>> runSuccessive ::
>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Foldable t) =>
>>>    CLP.Method ->
>>>    sh ->
>>>    Bounds ix ->
>>>    (Constraints ix, (Direction, Objective sh)) ->
>>>    t (Double -> Constraints ix, (Direction, Objective sh)) ->
>>>    Either CLP.ModelStatus ()
>>> runSuccessive method shape bounds (constrs,dirObj) objs =
>>>    let solve constrs_ dirObj_ = do
>>>          (status,result) <- LP.solve method constrs_ dirObj_
>>>          return $ maybe (Left status) Right result in
>>>    LP.run shape bounds $ ME.runExceptT $ do
>>>       (opt, _xs) <- ME.ExceptT $ solve constrs dirObj
>>>       traverse_Lag opt
>>>          (\prevResult (newConstr, dirObjI) -> do
>>>             (optI, _xs) <-
>>>                ME.ExceptT $ solve (newConstr prevResult) dirObjI
>>>             return optI)
>>>          objs
>>>
>>> solveSuccessiveWarm ::
>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>
>>>    CLP.Method ->
>>>    sh ->
>>>    Bounds ix ->
>>>    (Constraints ix, (Direction, Objective sh)) ->
>>>    t (Double -> Constraints ix, (Direction, Objective sh)) ->
>>>    Either CLP.ModelStatus (NonEmpty.T t (Double, Array sh Double))
>>> solveSuccessiveWarm method shape bounds (constrs,dirObj) objs =
>>>    let solve constrs_ dirObj_ = do
>>>          (status,result) <- LP.solve method constrs_ dirObj_
>>>          return $ maybe (Left status) Right result in
>>>    LP.run shape bounds $ ME.runExceptT $ do
>>>       result <- ME.ExceptT $ solve constrs dirObj
>>>       NonEmpty.Cons result <$>
>>>          traverseLag result
>>>             (\(prevOpt, _xs) (newConstr, dirObjI) ->
>>>                 ME.ExceptT $ solve (newConstr prevOpt) dirObjI)
>>>             objs
>>>
>>> solveSuccessiveGen ::
>>>    (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Traversable t) =>
>>>    CLP.Method ->
>>>    sh ->
>>>    Bounds ix ->
>>>    (Constraints ix, (Direction, Objective sh)) ->
>>>    t (Double -> Constraints ix, (Direction, Objective sh)) ->
>>>    Either CLP.ModelStatus (NonEmpty.T t (Double, Array sh Double))
>>> solveSuccessiveGen method shape bounds (constrs,dirObj) objs =
>>>    let solve bounds_ constrs_ dirObj_ =
>>>          case CLP.solve method bounds_ constrs_ dirObj_ of
>>>             (status,result) -> maybe (Left status) Right result in
>>>    LPMonad.run shape bounds $ ME.runExceptT $ do
>>>       result <- ME.ExceptT $ LPMonad.lift solve constrs dirObj
>>>       NonEmpty.Cons result <$>
>>>          traverseLag result
>>>             (\(prevOpt, _xs) (newConstr, dirObjI) ->
>>>                ME.ExceptT $ LPMonad.lift solve (newConstr prevOpt) dirObjI)
>>>             objs
-}


newtype T sh a = Cons (MR.ReaderT (sh, Ptr Highs) IO a)
   deriving (Functor, Applicative, Monad)

run ::
   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
   sh -> Bounds ix -> T sh a -> a
run shape bounds (Cons act) =
   unsafePerformIO $ runContT $ do
      model <- MC.ContT $ bracket Highs.create Highs.destroy
      liftIO $ Debug.initLog model
      startPtr <- withBuffer $ Array.vectorFromList [0]
      indexPtr <- withBuffer $ Array.vectorFromList []
      let numCols = Shape.size shape
      objPtr <- withBuffer $ Array.replicate (Shape.ZeroBased numCols) 0
      let emptyDoublePtr = objPtr
      (collbPtr,colubPtr) <-
         storeBounds $ prepareColumnBoundsArrays shape bounds
      liftIO $ checkStatus $ Highs.passLp model
         (fromIntegral numCols)
         0
         0
         Highs.matrixFormatRowwise
         Highs.objSenseMaximize
         0 objPtr
         collbPtr colubPtr
         emptyDoublePtr emptyDoublePtr
         startPtr indexPtr emptyDoublePtr
      liftIO $ MR.runReaderT act (shape, model)

{- |
Add new constraints to an existing problem
and run with a new direction and objective.

>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         fmap (mapSnd Array.toTuple) $ snd $
         LP.run tripletShape []
            (LP.solve LP.simplex
               [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20]
               (LP.Maximize, Array.fromTuple (4,-3,2)
                                 :: Array.Array TripletShape Double))
:}
Just (28.0,(5.0,0.0,4.0))

prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
   case (CLP.solve method bounds constrs (dir,obj),
         LP.run (Array.shape origin) bounds $
            LP.solve method constrs (dir,obj)) of
      ((_, Just (optA,_)), (_, Just (optB,_))) ->
         TestLP.approxReal 0.1 optA optB; _ -> False
:}

prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   TestLP.forAllObjectives origin $ \objs_ ->
   case TestLP.successiveObjectives origin 0.01 objs_ of
      (dirObj, objs) ->
         either
            (\msg -> QC.counterexample (show msg) False)
            (const $ QC.property True) $
         runSuccessive method (Array.shape origin) bounds (constrs,dirObj) objs
:}

prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   TestLP.forAllObjectives origin $ \objs_ ->
   let shape = Array.shape origin in
   case TestLP.successiveObjectives origin 0.01 objs_ of
      (dirObj, objs) ->
         approxSuccession 0.01
            (solveSuccessiveWarm method shape bounds (constrs,dirObj) objs)
            (solveSuccessiveGen method shape bounds (constrs,dirObj) objs)
:}
-}
solve ::
   (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>
   Method -> Constraints Double ix ->
   (Direction, Objective sh) -> T sh (Result sh)
solve = solveWith Priv.getResult

solveWith ::
   (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>
   Query sh result ->
   Method -> Constraints Double ix ->
   (Direction, Objective sh) -> T sh (LPEnum.ModelStatus, Maybe result)
solveWith query method constrs (dir,obj) = Cons $ do
   (shape, model) <- MR.ask
   when (shape /= Array.shape obj) $
      error "HiGHS.LP.Monad.solve: objective shape mismatch"

   let numRows = length constrs
   liftIO $ runContT $ do
      let (coefficients, indices, rowStarts) = prepareConstraints shape constrs
      (coefficientsPtr, indexPtr, startPtr)
         <- storeConstraints (coefficients, indices, rowStarts)
      (rowlbPtr,rowubPtr) <- storeBounds $ prepareRowBoundsArrays constrs
      objPtr <- withBuffer $ Array.map realToFrac obj
      liftIO $ do
         checkStatus $ Highs.addRows model (fromIntegral numRows)
            rowlbPtr rowubPtr
            (fromIntegral $ Shape.size $ Array.shape coefficients)
            startPtr indexPtr coefficientsPtr
         let numCols = Shape.size shape
         when (numCols>0) $ checkStatus $
            Highs.changeColsCostByRange model
               0 (fromIntegral numCols - 1) objPtr

   liftIO $ do
      setMethod model method
      checkStatus $ Highs.changeObjectiveSense model $ objectiveSense dir
      examineStatus query shape model =<< Highs.run model