{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
module Numeric.HiGHS.LP (
   -- * simple solver
   solve,
   LP.Direction(..),
   Term(..), (LP..*),
   Constraints,
   LP.free, (LP.<=.), (LP.>=.), (LP.==.), (LP.>=<.),
   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

-- ToDo: extend .* to functions of indices using a type class
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.Private
         (Method, Result, Query, checkStatus, runContT, withBuffer,
          storeBounds, prepareRowBoundsArrays, prepareColumnBoundsArrays,
          storeConstraints, prepareConstraints,
          setMethod, objectiveSense, examineStatus)

import qualified Numeric.LinearProgramming.Common as LP
import Numeric.LinearProgramming.Common
         (Bounds, Term(Term), Constraints, Direction(..), Objective)

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

import Control.Monad.IO.Class (liftIO)
import Control.Exception (bracket)

import System.IO.Unsafe (unsafePerformIO)



{- $setup
>>> import qualified Numeric.HiGHS.LP as LP
>>> import qualified Numeric.LinearProgramming.Test as TestLP
>>> import Numeric.HiGHS.LP ((.*), (==.), (<=.), (>=.), (>=<.))
>>>
>>> import qualified Data.Array.Comfort.Storable as Array
>>> import qualified Data.Array.Comfort.Shape as Shape
>>> import Data.Tuple.HT (mapPair, mapSnd)
>>>
>>> import Control.Applicative (liftA2)
>>>
>>> import qualified Test.QuickCheck as QC
>>> import Test.QuickCheck ((===), (.&&.), (.||.))
>>>
>>> type X = Shape.Element
>>> type PairShape = Shape.NestedTuple Shape.TupleIndex (X,X)
>>> type TripletShape = Shape.NestedTuple Shape.TupleIndex (X,X,X)
>>>
>>> pairShape :: PairShape
>>> pairShape = Shape.static
>>>
>>> tripletShape :: TripletShape
>>> tripletShape = Shape.static
>>>
>>> approxReal :: (Ord a, Num a) => a -> a -> a -> Bool
>>> approxReal tol x y = abs (x-y) <= tol
>>>
>>> forAllMethod ::
>>>    (QC.Testable prop) => (LP.Method -> prop) -> QC.Property
>>> forAllMethod = QC.forAll (QC.elements [LP.simplex, LP.choose, LP.ipm])
-}



{- |
>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         fmap (mapSnd Array.toTuple) $ snd $
         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))

>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         fmap (mapSnd Array.toTuple) $ snd $
         LP.solve LP.choose [y >=<. (-12,12)]
            [[1.*x, (-1).*y] <=. 10, [(-1).*y, (1::Double).*z] <=. 20]
            (LP.Maximize, Array.fromTuple (4,-3,2)
               :: Array.Array TripletShape Double)
:}
Just (116.0,(22.0,12.0,32.0))

>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         mapSnd (fmap (mapSnd Array.toTuple)) $
         LP.solve LP.choose [y >=<. (-12,12)]
            [[1.*x, 1.*y] <=. 10, [1.*y, (-1::Double).*z] >=. 20]
            (LP.Maximize, Array.fromTuple (4,3,2)
               :: Array.Array TripletShape Double)
:}
(ModelStatusInfeasible,...)

>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         mapSnd (fmap (mapSnd Array.toTuple)) $
         LP.solve LP.choose [y >=<. (-12,12)]
            [[1.*x, 1.*y] <=. 10, [1.*y, (1::Double).*z] >=. 20]
            (LP.Maximize, Array.fromTuple (4,3,2)
               :: Array.Array TripletShape Double)
:}
(ModelStatusUnbounded,...)

prop> :{
   forAllMethod $ \method
      (QC.Positive posWeight) (QC.Positive negWeight) target ->
   case Shape.indexTupleFromShape pairShape of
      (pos,neg) ->
         case fmap (mapSnd Array.toTuple) $ snd $ LP.solve method []
                  [[1.*pos, (-1::Double).*neg] ==. target]
                  (LP.Minimize,
                   Array.fromTuple (posWeight,negWeight)
                     :: Array.Array PairShape Double) of
            Nothing -> QC.property False
            Just (absol,(posResult,negResult)) ->
               QC.property (absol>=0)
               .&&.
               (posResult === 0 .||. negResult === 0)
:}
prop> :{
   forAllMethod $ \method target ->
   case Shape.indexTupleFromShape pairShape of
      (pos,neg) ->
         case fmap (mapSnd Array.toTuple) $ snd $ LP.solve method []
                  [[1.*pos, (-1::Double).*neg] ==. target]
                  (LP.Minimize, Array.fromTuple (1,1)
                     :: Array.Array PairShape Double) of
            Nothing -> QC.property False
            Just (absol,(posResult,negResult)) ->
               QC.counterexample (show(absol,(posResult,negResult))) $
               QC.property (approxReal 0.001 absol (abs target))
               .&&.
               (posResult === 0 .||. negResult === 0)
:}

prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
   case snd $ LP.solve method bounds constrs (dir,obj) of
      Nothing -> False
      Just _ -> True
:}
prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
   case snd $ LP.solve method bounds constrs (dir,obj) of
      Nothing -> QC.property False
      Just (_,sol) -> TestLP.checkFeasibility 0.1 bounds constrs sol
:}
prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
   case snd $ LP.solve method bounds constrs (dir,obj) of
      Nothing -> QC.property False
      Just (_,sol) ->
         QC.forAll (QC.choose (0,1)) $ \lambda ->
         TestLP.checkFeasibility 0.1 bounds constrs $
         TestLP.affineCombination lambda sol (Array.map fromIntegral origin)
:}
prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
   case snd $ LP.solve method bounds constrs (dir,obj) of
      Nothing -> QC.property False
      Just (opt,sol) ->
         QC.forAll (QC.choose (0,1)) $ \lambda ->
            let val = TestLP.scalarProduct obj $
                        TestLP.affineCombination lambda sol $
                        Array.map fromIntegral origin
            in QC.counterexample (show (dir,opt,val)) $
               case dir of
                  LP.Minimize -> opt-0.01 <= val
                  LP.Maximize -> opt+0.01 >= val
:}
prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllBoundedProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \dirObjA ->
   QC.forAll (TestLP.genObjective origin) $ \dirObjB ->
   let solA = snd $ LP.solve method bounds constrs dirObjA in
   let solB = snd $ LP.solve method bounds constrs dirObjB in
   QC.counterexample (show (fmap fst solA, fmap fst solB)) $
   case (solA, solB) of
      (Just _, Just _) -> True
      (Nothing, Nothing) -> True
      _ -> False
:}
prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(_dir,obj) ->
   case (snd $ LP.solve method bounds constrs (LP.Minimize,obj),
         snd $ LP.solve method bounds constrs (LP.Maximize,obj)) of
      (Just (optMin,_), Just (optMax,_)) ->
         QC.counterexample (show (optMin, optMax)) $ optMin <= optMax + 0.01
      _ -> QC.property False
:}
prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds allConstrs ->
   QC.forAll (QC.sublistOf allConstrs) $ \someConstrs ->
   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
   case (snd $ LP.solve method bounds allConstrs (dir,obj),
         snd $ LP.solve method bounds someConstrs (dir,obj)) of
      (Just (optAll,_), Just (optSome,_)) ->
         QC.counterexample (show (optAll, optSome)) $
         case dir of
            LP.Minimize -> optAll >= optSome-0.01
            LP.Maximize -> optAll <= optSome+0.01
      _ -> QC.property False
:}
prop> :{
   forAllMethod $ \methodA ->
   forAllMethod $ \methodB ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \dirObj ->
   case (snd $ LP.solve methodA bounds constrs dirObj,
         snd $ LP.solve methodB bounds constrs dirObj) of
      (Just (optA,_), Just (optB,_)) ->
         QC.counterexample (show (optA, optB)) $
         approxReal 0.01 optA optB
      _ -> QC.property False
:}
-}
solve ::
   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
   Method -> Bounds ix -> Constraints Double ix ->
   (Direction, Objective sh) -> Result sh
solve :: forall sh ix.
(Indexed sh, Index sh ~ ix) =>
Method
-> Bounds ix
-> Constraints Double ix
-> (Direction, Objective sh)
-> Result sh
solve = Query sh (Double, Array sh Double)
-> Method
-> Bounds ix
-> Constraints Double ix
-> (Direction, Array sh Double)
-> (ModelStatus, Maybe (Double, Array sh Double))
forall sh ix result.
(Indexed sh, Index sh ~ ix) =>
Query sh result
-> Method
-> Bounds ix
-> Constraints Double ix
-> (Direction, Objective sh)
-> (ModelStatus, Maybe result)
solveWith Query sh (Double, Array sh Double)
forall sh. C sh => Query sh (Double, Array sh Double)
Priv.getResult

{- |
>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         fmap (mapSnd (mapPair (Array.toTuple, Array.toList))) $ snd $
         LP.solveWith
            (liftA2 (,) LP.getObjectiveValue LP.getBasisStatus)
            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,((Highs.basisStatusBasic,Highs.basisStatusLower,Highs.basisStatusBasic),[Highs.basisStatusUpper,Highs.basisStatusUpper]))

prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
   case (snd $ LP.solve method bounds constrs (dir,obj),
         snd $ LP.solveWith LP.getSolutionVectors
                  method bounds constrs (dir,obj)) of
      (Just (_,sol0), Just ((sol1,_),_)) -> sol0 == sol1
      _ -> False
:}
-}
solveWith ::
   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
   Query sh result -> Method -> Bounds ix -> Constraints Double ix ->
   (Direction, Objective sh) -> (LPEnum.ModelStatus, Maybe result)
solveWith :: forall sh ix result.
(Indexed sh, Index sh ~ ix) =>
Query sh result
-> Method
-> Bounds ix
-> Constraints Double ix
-> (Direction, Objective sh)
-> (ModelStatus, Maybe result)
solveWith Query sh result
query Method
method Bounds ix
bounds Constraints Double ix
constrs (Direction
dir,Objective sh
obj) =
   IO (ModelStatus, Maybe result) -> (ModelStatus, Maybe result)
forall a. IO a -> a
unsafePerformIO (IO (ModelStatus, Maybe result) -> (ModelStatus, Maybe result))
-> IO (ModelStatus, Maybe result) -> (ModelStatus, Maybe result)
forall a b. (a -> b) -> a -> b
$
   let shape :: sh
shape = Objective sh -> sh
forall sh a. Array sh a -> sh
Array.shape Objective sh
obj in
   let numCols :: Int
numCols = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
shape in
   let numRows :: Int
numRows = Constraints Double ix -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Constraints Double ix
constrs in

   ContT (ModelStatus, Maybe result) IO (ModelStatus, Maybe result)
-> IO (ModelStatus, Maybe result)
forall a. ContT a IO a -> IO a
runContT (ContT (ModelStatus, Maybe result) IO (ModelStatus, Maybe result)
 -> IO (ModelStatus, Maybe result))
-> ContT (ModelStatus, Maybe result) IO (ModelStatus, Maybe result)
-> IO (ModelStatus, Maybe result)
forall a b. (a -> b) -> a -> b
$ do
   Ptr CDouble
objPtr <- Array sh CDouble
-> ContT (ModelStatus, Maybe result) IO (Ptr CDouble)
forall sh a r. Array sh a -> ContT r IO (Ptr a)
withBuffer (Array sh CDouble
 -> ContT (ModelStatus, Maybe result) IO (Ptr CDouble))
-> Array sh CDouble
-> ContT (ModelStatus, Maybe result) IO (Ptr CDouble)
forall a b. (a -> b) -> a -> b
$ (Double -> CDouble) -> Objective sh -> Array sh CDouble
forall sh a b.
(C sh, Storable a, Storable b) =>
(a -> b) -> Array sh a -> Array sh b
Array.map Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac Objective sh
obj
   (Ptr CDouble
collbPtr,Ptr CDouble
colubPtr) <-
      (Array sh CDouble, Array sh CDouble)
-> ContT (ModelStatus, Maybe result) IO (Ptr CDouble, Ptr CDouble)
forall sh r.
(Array sh CDouble, Array sh CDouble)
-> ContT r IO (Ptr CDouble, Ptr CDouble)
storeBounds ((Array sh CDouble, Array sh CDouble)
 -> ContT (ModelStatus, Maybe result) IO (Ptr CDouble, Ptr CDouble))
-> (Array sh CDouble, Array sh CDouble)
-> ContT (ModelStatus, Maybe result) IO (Ptr CDouble, Ptr CDouble)
forall a b. (a -> b) -> a -> b
$ sh -> Bounds ix -> (Array sh CDouble, Array sh CDouble)
forall sh ix.
(Indexed sh, Index sh ~ ix) =>
sh -> Bounds ix -> (Array sh CDouble, Array sh CDouble)
prepareColumnBoundsArrays sh
shape Bounds ix
bounds
   (Ptr CDouble
rowlbPtr,Ptr CDouble
rowubPtr) <- (Array ShapeInt CDouble, Array ShapeInt CDouble)
-> ContT (ModelStatus, Maybe result) IO (Ptr CDouble, Ptr CDouble)
forall sh r.
(Array sh CDouble, Array sh CDouble)
-> ContT r IO (Ptr CDouble, Ptr CDouble)
storeBounds ((Array ShapeInt CDouble, Array ShapeInt CDouble)
 -> ContT (ModelStatus, Maybe result) IO (Ptr CDouble, Ptr CDouble))
-> (Array ShapeInt CDouble, Array ShapeInt CDouble)
-> ContT (ModelStatus, Maybe result) IO (Ptr CDouble, Ptr CDouble)
forall a b. (a -> b) -> a -> b
$ Constraints Double ix
-> (Array ShapeInt CDouble, Array ShapeInt CDouble)
forall ix.
Bounds ix -> (Array ShapeInt CDouble, Array ShapeInt CDouble)
prepareRowBoundsArrays Constraints Double ix
constrs

   let (Array ShapeInt CDouble
coefficients, Array ShapeInt HighsInt
indices, Array ShapeInt HighsInt
rowStarts) = sh
-> Constraints Double ix
-> (Array ShapeInt CDouble, Array ShapeInt HighsInt,
    Array ShapeInt HighsInt)
forall sh ix.
(Indexed sh, Index sh ~ ix) =>
sh
-> Constraints Double ix
-> (Array ShapeInt CDouble, Array ShapeInt HighsInt,
    Array ShapeInt HighsInt)
prepareConstraints sh
shape Constraints Double ix
constrs
   (Ptr CDouble
coefficientsPtr, Ptr HighsInt
indexPtr, Ptr HighsInt
startPtr)
      <- (Array ShapeInt CDouble, Array ShapeInt HighsInt,
 Array ShapeInt HighsInt)
-> ContT
     (ModelStatus, Maybe result)
     IO
     (Ptr CDouble, Ptr HighsInt, Ptr HighsInt)
forall r.
(Array ShapeInt CDouble, Array ShapeInt HighsInt,
 Array ShapeInt HighsInt)
-> ContT r IO (Ptr CDouble, Ptr HighsInt, Ptr HighsInt)
storeConstraints (Array ShapeInt CDouble
coefficients, Array ShapeInt HighsInt
indices, Array ShapeInt HighsInt
rowStarts)

   IO (ModelStatus, Maybe result)
-> ContT (ModelStatus, Maybe result) IO (ModelStatus, Maybe result)
forall a. IO a -> ContT (ModelStatus, Maybe result) IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (ModelStatus, Maybe result)
 -> ContT
      (ModelStatus, Maybe result) IO (ModelStatus, Maybe result))
-> IO (ModelStatus, Maybe result)
-> ContT (ModelStatus, Maybe result) IO (ModelStatus, Maybe result)
forall a b. (a -> b) -> a -> b
$ IO (Ptr Highs)
-> (Ptr Highs -> IO ())
-> (Ptr Highs -> IO (ModelStatus, Maybe result))
-> IO (ModelStatus, Maybe result)
forall a b c. IO a -> (a -> IO b) -> (a -> IO c) -> IO c
bracket IO (Ptr Highs)
Highs.create Ptr Highs -> IO ()
Highs.destroy ((Ptr Highs -> IO (ModelStatus, Maybe result))
 -> IO (ModelStatus, Maybe result))
-> (Ptr Highs -> IO (ModelStatus, Maybe result))
-> IO (ModelStatus, Maybe result)
forall a b. (a -> b) -> a -> b
$ \Ptr Highs
model -> do
      Ptr Highs -> IO ()
Debug.initLog Ptr Highs
model
      Ptr Highs -> Method -> IO ()
setMethod Ptr Highs
model Method
method
      IO Status -> IO ()
checkStatus (IO Status -> IO ()) -> IO Status -> IO ()
forall a b. (a -> b) -> a -> b
$ Ptr Highs
-> HighsInt
-> HighsInt
-> HighsInt
-> MatrixFormat
-> ObjSense
-> CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr HighsInt
-> Ptr HighsInt
-> Ptr CDouble
-> IO Status
Highs.passLp Ptr Highs
model
         (Int -> HighsInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
numCols)
         (Int -> HighsInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
numRows)
         (Int -> HighsInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> HighsInt) -> Int -> HighsInt
forall a b. (a -> b) -> a -> b
$ ShapeInt -> Int
forall sh. C sh => sh -> Int
Shape.size (ShapeInt -> Int) -> ShapeInt -> Int
forall a b. (a -> b) -> a -> b
$ Array ShapeInt CDouble -> ShapeInt
forall sh a. Array sh a -> sh
Array.shape Array ShapeInt CDouble
coefficients)
         MatrixFormat
Highs.matrixFormatRowwise
         (Direction -> ObjSense
objectiveSense Direction
dir)
         CDouble
0 Ptr CDouble
objPtr
         Ptr CDouble
collbPtr Ptr CDouble
colubPtr
         Ptr CDouble
rowlbPtr Ptr CDouble
rowubPtr
         Ptr HighsInt
startPtr Ptr HighsInt
indexPtr Ptr CDouble
coefficientsPtr

      Query sh result
-> sh -> Ptr Highs -> Status -> IO (ModelStatus, Maybe result)
forall sh a.
C sh =>
Query sh a
-> sh -> Ptr Highs -> Status -> IO (ModelStatus, Maybe a)
examineStatus Query sh result
query sh
shape Ptr Highs
model (Status -> IO (ModelStatus, Maybe result))
-> IO Status -> IO (ModelStatus, Maybe result)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr Highs -> IO Status
Highs.run Ptr Highs
model