{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
module Numeric.COINOR.CLP (
   simplex,
   LP.Direction(..),
   PlusMinusOne(..),
   Term(..), (LP..*),
   Constraints,
   LP.free, (LP.<=.), (LP.>=.), (LP.==.), (LP.>=<.),
   Method, Priv.dual, Priv.primal,
   Priv.initialSolve, Priv.initialDualSolve, Priv.initialPrimalSolve,
   Priv.initialBarrierSolve, Priv.initialBarrierNoCrossSolve,
   FailureType(..),
   Result,
   ) where

import qualified Numeric.COINOR.CLP.FFI as FFI
import qualified Numeric.COINOR.CLP.Debug as Debug
import qualified Numeric.COINOR.CLP.Private as Priv
import Numeric.COINOR.CLP.Private
         (Method(runMethod), Result, FailureType(..),
          runContT, withBuffer, false,
          storeBounds, prepareRowBoundsArrays, prepareColumnBoundsArrays,
          storeConstraints, prepareConstraints,
          setOptimizationDirection, examineStatus)

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

import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape
import qualified Data.List.HT as ListHT

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

import System.IO.Unsafe (unsafePerformIO)

import Foreign.Ptr (Ptr, nullPtr)
import Foreign.C.Types (CDouble)


{- $setup
>>> import qualified Numeric.COINOR.CLP as LP
>>> import qualified Numeric.LinearProgramming.Test as TestLP
>>> import Numeric.COINOR.CLP
>>>    (PlusMinusOne(..), (.*), (==.), (<=.), (>=.), (>=<.))
>>>
>>> import qualified Data.Array.Comfort.Storable as Array
>>> import qualified Data.Array.Comfort.Shape as Shape
>>>
>>> import Data.Either.HT (mapRight)
>>> import Data.Tuple.HT (mapSnd)
>>>
>>> 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
>>>
>>> genMethod :: QC.Gen (String, LP.Method)
>>> genMethod = QC.elements $
>>>    ("dual", LP.dual) :
>>>    ("primal", LP.primal) :
>>>    ("initialSolve", LP.initialSolve) :
>>>    ("initialDualSolve", LP.initialDualSolve) :
>>>    ("initialPrimalSolve", LP.initialPrimalSolve) :
>>>    ("initialBarrierSolve", LP.initialBarrierSolve) :
>>>    -- let tests fail
>>>    -- ("initialBarrierNoCrossSolve", LP.initialBarrierNoCrossSolve) :
>>>    []
>>>
>>> forAllMethod ::
>>>    (QC.Testable prop) => (LP.Method -> prop) -> QC.Property
>>> forAllMethod prop = QC.forAllShow genMethod fst (prop . snd)
-}



data PlusMinusOne = MinusOne | PlusOne deriving (PlusMinusOne -> PlusMinusOne -> Bool
(PlusMinusOne -> PlusMinusOne -> Bool)
-> (PlusMinusOne -> PlusMinusOne -> Bool) -> Eq PlusMinusOne
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: PlusMinusOne -> PlusMinusOne -> Bool
== :: PlusMinusOne -> PlusMinusOne -> Bool
$c/= :: PlusMinusOne -> PlusMinusOne -> Bool
/= :: PlusMinusOne -> PlusMinusOne -> Bool
Eq, Int -> PlusMinusOne -> ShowS
[PlusMinusOne] -> ShowS
PlusMinusOne -> String
(Int -> PlusMinusOne -> ShowS)
-> (PlusMinusOne -> String)
-> ([PlusMinusOne] -> ShowS)
-> Show PlusMinusOne
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> PlusMinusOne -> ShowS
showsPrec :: Int -> PlusMinusOne -> ShowS
$cshow :: PlusMinusOne -> String
show :: PlusMinusOne -> String
$cshowList :: [PlusMinusOne] -> ShowS
showList :: [PlusMinusOne] -> ShowS
Show)


class Coefficient a where
   loadProblem ::
      (Shape.Indexed sh, Shape.Index sh ~ ix) =>
      sh ->
      Constraints a ix ->
      Ptr FFI.Simplex ->
      Ptr CDouble -> Ptr CDouble ->
      Ptr CDouble ->
      Ptr CDouble -> Ptr CDouble ->
      MC.ContT () IO ()

instance Coefficient Double where
   loadProblem :: forall sh ix.
(Indexed sh, Index sh ~ ix) =>
sh
-> Constraints Double ix
-> Ptr Simplex
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> ContT () IO ()
loadProblem sh
shape Constraints Double ix
constrs Ptr Simplex
lp Ptr CDouble
collbPtr Ptr CDouble
colubPtr Ptr CDouble
objPtr Ptr CDouble
rowlbPtr Ptr CDouble
rowubPtr = do
      let (Array ShapeInt CDouble
coefficients, Array ShapeInt CInt
indices, Array ShapeInt BigIndex
rowStarts) = sh
-> Constraints Double ix
-> (Array ShapeInt CDouble, Array ShapeInt CInt,
    Array ShapeInt BigIndex)
forall sh ix.
(Indexed sh, Index sh ~ ix) =>
sh
-> Constraints Double ix
-> (Array ShapeInt CDouble, Array ShapeInt CInt,
    Array ShapeInt BigIndex)
prepareConstraints sh
shape Constraints Double ix
constrs
      (Ptr CDouble
coefficientsPtr, Ptr CInt
indexPtr, Ptr BigIndex
startPtr)
         <- (Array ShapeInt CDouble, Array ShapeInt CInt,
 Array ShapeInt BigIndex)
-> ContT () IO (Ptr CDouble, Ptr CInt, Ptr BigIndex)
forall r.
(Array ShapeInt CDouble, Array ShapeInt CInt,
 Array ShapeInt BigIndex)
-> ContT r IO (Ptr CDouble, Ptr CInt, Ptr BigIndex)
storeConstraints (Array ShapeInt CDouble
coefficients, Array ShapeInt CInt
indices, Array ShapeInt BigIndex
rowStarts)
      let createMatrix :: IO (Ptr CoinPackedMatrix)
createMatrix =
            CBool
-> CInt
-> CInt
-> BigIndex
-> Ptr CDouble
-> Ptr CInt
-> Ptr BigIndex
-> Ptr CInt
-> IO (Ptr CoinPackedMatrix)
FFI.newCoinPackedMatrix
               CBool
false
               (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CInt) -> Int -> CInt
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
shape)
               (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CInt) -> Int -> CInt
forall a b. (a -> b) -> a -> b
$ Constraints Double ix -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Constraints Double ix
constrs)
               (Int -> BigIndex
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> BigIndex) -> Int -> BigIndex
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)
               Ptr CDouble
coefficientsPtr
               Ptr CInt
indexPtr
               Ptr BigIndex
startPtr
               Ptr CInt
forall a. Ptr a
nullPtr
      Ptr CoinPackedMatrix
matrix <- ((Ptr CoinPackedMatrix -> IO ()) -> IO ())
-> ContT () IO (Ptr CoinPackedMatrix)
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
MC.ContT (((Ptr CoinPackedMatrix -> IO ()) -> IO ())
 -> ContT () IO (Ptr CoinPackedMatrix))
-> ((Ptr CoinPackedMatrix -> IO ()) -> IO ())
-> ContT () IO (Ptr CoinPackedMatrix)
forall a b. (a -> b) -> a -> b
$ IO (Ptr CoinPackedMatrix)
-> (Ptr CoinPackedMatrix -> IO ())
-> (Ptr CoinPackedMatrix -> IO ())
-> IO ()
forall a b c. IO a -> (a -> IO b) -> (a -> IO c) -> IO c
bracket IO (Ptr CoinPackedMatrix)
createMatrix Ptr CoinPackedMatrix -> IO ()
FFI.deleteCoinPackedMatrix
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr Simplex
-> Ptr CoinPackedMatrix
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> IO ()
FFI.loadProblemFromCoinMatrix Ptr Simplex
lp Ptr CoinPackedMatrix
matrix
            Ptr CDouble
collbPtr Ptr CDouble
colubPtr
            Ptr CDouble
objPtr
            Ptr CDouble
rowlbPtr Ptr CDouble
rowubPtr
            Ptr CDouble
forall a. Ptr a
nullPtr

instance Coefficient PlusMinusOne where
   loadProblem :: forall sh ix.
(Indexed sh, Index sh ~ ix) =>
sh
-> Constraints PlusMinusOne ix
-> Ptr Simplex
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> ContT () IO ()
loadProblem sh
shape Constraints PlusMinusOne ix
constrs Ptr Simplex
lp Ptr CDouble
collbPtr Ptr CDouble
colubPtr Ptr CDouble
objPtr Ptr CDouble
rowlbPtr Ptr CDouble
rowubPtr = do
      let shapeOffset :: Index sh -> Int
shapeOffset = sh -> Index sh -> Int
forall sh. Indexed sh => sh -> Index sh -> Int
Shape.offset sh
shape
      let coefficients :: [([Term PlusMinusOne ix], [Term PlusMinusOne ix])]
coefficients =
            (Inequality [Term PlusMinusOne ix]
 -> ([Term PlusMinusOne ix], [Term PlusMinusOne ix]))
-> Constraints PlusMinusOne ix
-> [([Term PlusMinusOne ix], [Term PlusMinusOne ix])]
forall a b. (a -> b) -> [a] -> [b]
map
               (\(Inequality [Term PlusMinusOne ix]
terms Bound
_bnd) ->
                  (Term PlusMinusOne ix -> Bool)
-> [Term PlusMinusOne ix]
-> ([Term PlusMinusOne ix], [Term PlusMinusOne ix])
forall a. (a -> Bool) -> [a] -> ([a], [a])
ListHT.partition (\(Term PlusMinusOne
c ix
_) -> PlusMinusOne
c PlusMinusOne -> PlusMinusOne -> Bool
forall a. Eq a => a -> a -> Bool
== PlusMinusOne
PlusOne) [Term PlusMinusOne ix]
terms)
               Constraints PlusMinusOne ix
constrs
      Ptr CInt
indexPtr <-
         Array ShapeInt CInt -> ContT () IO (Ptr CInt)
forall sh a r. Array sh a -> ContT r IO (Ptr a)
withBuffer (Array ShapeInt CInt -> ContT () IO (Ptr CInt))
-> Array ShapeInt CInt -> ContT () IO (Ptr CInt)
forall a b. (a -> b) -> a -> b
$ [CInt] -> Array ShapeInt CInt
forall a. Storable a => [a] -> Array ShapeInt a
Array.vectorFromList ([CInt] -> Array ShapeInt CInt) -> [CInt] -> Array ShapeInt CInt
forall a b. (a -> b) -> a -> b
$
         (([Term PlusMinusOne ix], [Term PlusMinusOne ix]) -> [CInt])
-> [([Term PlusMinusOne ix], [Term PlusMinusOne ix])] -> [CInt]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap
            (\([Term PlusMinusOne ix]
positive,[Term PlusMinusOne ix]
negative) ->
               (Int -> CInt) -> [Int] -> [CInt]
forall a b. (a -> b) -> [a] -> [b]
map Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> [CInt]) -> [Int] -> [CInt]
forall a b. (a -> b) -> a -> b
$
                  (Term PlusMinusOne ix -> Int) -> [Term PlusMinusOne ix] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (\(Term PlusMinusOne
_ ix
ix) -> Index sh -> Int
shapeOffset ix
Index sh
ix) [Term PlusMinusOne ix]
positive
                  [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++
                  (Term PlusMinusOne ix -> Int) -> [Term PlusMinusOne ix] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (\(Term PlusMinusOne
_ ix
ix) -> Index sh -> Int
shapeOffset ix
Index sh
ix) [Term PlusMinusOne ix]
negative)
            [([Term PlusMinusOne ix], [Term PlusMinusOne ix])]
coefficients
      let rowStarts :: [Int]
rowStarts =
            (Int -> Int -> Int) -> Int -> [Int] -> [Int]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) Int
0 ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$
            (Inequality [Term PlusMinusOne ix] -> Int)
-> Constraints PlusMinusOne ix -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (\(Inequality [Term PlusMinusOne ix]
terms Bound
_bnd) -> [Term PlusMinusOne ix] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Term PlusMinusOne ix]
terms) Constraints PlusMinusOne ix
constrs
      Ptr BigIndex
startPositivePtr <-
         Array ShapeInt BigIndex -> ContT () IO (Ptr BigIndex)
forall sh a r. Array sh a -> ContT r IO (Ptr a)
withBuffer (Array ShapeInt BigIndex -> ContT () IO (Ptr BigIndex))
-> Array ShapeInt BigIndex -> ContT () IO (Ptr BigIndex)
forall a b. (a -> b) -> a -> b
$ [BigIndex] -> Array ShapeInt BigIndex
forall a. Storable a => [a] -> Array ShapeInt a
Array.vectorFromList ([BigIndex] -> Array ShapeInt BigIndex)
-> [BigIndex] -> Array ShapeInt BigIndex
forall a b. (a -> b) -> a -> b
$ (Int -> BigIndex) -> [Int] -> [BigIndex]
forall a b. (a -> b) -> [a] -> [b]
map Int -> BigIndex
forall a b. (Integral a, Num b) => a -> b
fromIntegral [Int]
rowStarts
      Ptr BigIndex
startNegativePtr <-
         Array ShapeInt BigIndex -> ContT () IO (Ptr BigIndex)
forall sh a r. Array sh a -> ContT r IO (Ptr a)
withBuffer (Array ShapeInt BigIndex -> ContT () IO (Ptr BigIndex))
-> Array ShapeInt BigIndex -> ContT () IO (Ptr BigIndex)
forall a b. (a -> b) -> a -> b
$ [BigIndex] -> Array ShapeInt BigIndex
forall a. Storable a => [a] -> Array ShapeInt a
Array.vectorFromList ([BigIndex] -> Array ShapeInt BigIndex)
-> [BigIndex] -> Array ShapeInt BigIndex
forall a b. (a -> b) -> a -> b
$
         (Int
 -> ([Term PlusMinusOne ix], [Term PlusMinusOne ix]) -> BigIndex)
-> [Int]
-> [([Term PlusMinusOne ix], [Term PlusMinusOne ix])]
-> [BigIndex]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
k ([Term PlusMinusOne ix]
pos,[Term PlusMinusOne ix]
_neg) -> Int -> BigIndex
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> BigIndex) -> Int -> BigIndex
forall a b. (a -> b) -> a -> b
$ Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Term PlusMinusOne ix] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Term PlusMinusOne ix]
pos)
            [Int]
rowStarts [([Term PlusMinusOne ix], [Term PlusMinusOne ix])]
coefficients
      let createMatrix :: IO (Ptr PlusMinusOneMatrix)
createMatrix =
            CInt
-> CInt
-> CBool
-> Ptr CInt
-> Ptr BigIndex
-> Ptr BigIndex
-> IO (Ptr PlusMinusOneMatrix)
FFI.newPlusMinusOneMatrix
               (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CInt) -> Int -> CInt
forall a b. (a -> b) -> a -> b
$ Constraints PlusMinusOne ix -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Constraints PlusMinusOne ix
constrs)
               (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CInt) -> Int -> CInt
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
shape)
               (Int -> CBool
forall a. Enum a => Int -> a
toEnum (Int -> CBool) -> Int -> CBool
forall a b. (a -> b) -> a -> b
$ Bool -> Int
forall a. Enum a => a -> Int
fromEnum Bool
False)
               Ptr CInt
indexPtr Ptr BigIndex
startPositivePtr Ptr BigIndex
startNegativePtr
      Ptr PlusMinusOneMatrix
matrix <- ((Ptr PlusMinusOneMatrix -> IO ()) -> IO ())
-> ContT () IO (Ptr PlusMinusOneMatrix)
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
MC.ContT (((Ptr PlusMinusOneMatrix -> IO ()) -> IO ())
 -> ContT () IO (Ptr PlusMinusOneMatrix))
-> ((Ptr PlusMinusOneMatrix -> IO ()) -> IO ())
-> ContT () IO (Ptr PlusMinusOneMatrix)
forall a b. (a -> b) -> a -> b
$ IO (Ptr PlusMinusOneMatrix)
-> (Ptr PlusMinusOneMatrix -> IO ())
-> (Ptr PlusMinusOneMatrix -> IO ())
-> IO ()
forall a b c. IO a -> (a -> IO b) -> (a -> IO c) -> IO c
bracket IO (Ptr PlusMinusOneMatrix)
createMatrix Ptr PlusMinusOneMatrix -> IO ()
FFI.deletePlusMinusOneMatrix
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr Simplex
-> Ptr PlusMinusOneMatrix
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> IO ()
forall matrix.
Ptr Simplex
-> Ptr matrix
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> IO ()
FFI.loadProblemFromMatrix Ptr Simplex
lp Ptr PlusMinusOneMatrix
matrix
            Ptr CDouble
collbPtr Ptr CDouble
colubPtr
            Ptr CDouble
objPtr
            Ptr CDouble
rowlbPtr Ptr CDouble
rowubPtr
            Ptr CDouble
forall a. Ptr a
nullPtr


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

>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         mapSnd Array.toTuple <$>
         LP.simplex LP.primal [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)
:}
Right (116.0,(22.0,12.0,32.0))

>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         mapSnd Array.toTuple <$>
         LP.simplex LP.primal [y >=<. (-12,12)]
            [[PlusOne .* x, MinusOne .* y] <=. 10,
             [MinusOne .* y, PlusOne .* z] <=. 20]
            (LP.Maximize, Array.fromTuple (4,-3,2)
               :: Array.Array TripletShape Double)
:}
Right (116.0,(22.0,12.0,32.0))

>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         mapSnd Array.toTuple <$>
         LP.simplex LP.primal [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)
:}
Left PrimalInfeasible

>>> :{
   case Shape.indexTupleFromShape tripletShape of
      (x,y,z) ->
         mapSnd Array.toTuple <$>
         LP.simplex LP.primal [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)
:}
Left DualInfeasible

prop> :{
   forAllMethod $ \method (QC.Positive posWeight) (QC.Positive negWeight) target ->
   case Shape.indexTupleFromShape pairShape of
      (pos,neg) ->
         case mapSnd Array.toTuple <$>
               LP.simplex method []
                  [[1.*pos, (-1::Double).*neg] ==. target]
                  (LP.Minimize, Array.fromTuple (posWeight,negWeight)
                     :: Array.Array PairShape Double) of
            Left _ -> QC.property False
            Right (absol,(posResult,negResult)) ->
               QC.property (absol>=0)
               .&&.
               (posResult === 0 .||. negResult === 0)
:}
prop> :{
   forAllMethod $ \method target ->
   case Shape.indexTupleFromShape pairShape of
      (pos,neg) ->
         case mapSnd Array.toTuple <$>
               LP.simplex method []
                  [[1.*pos, (-1::Double).*neg] ==. target]
                  (LP.Minimize, Array.fromTuple (1,1)
                     :: Array.Array PairShape Double) of
            Left _ -> QC.property False
            Right (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 LP.simplex method bounds constrs (dir,obj) of
      Left _ -> False
      Right _ -> True
:}
prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(dir,obj) ->
   case LP.simplex method bounds constrs (dir,obj) of
      Left _ -> QC.property False
      Right (_,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 LP.simplex method bounds constrs (dir,obj) of
      Left _ -> QC.property False
      Right (_,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 LP.simplex method bounds constrs (dir,obj) of
      Left _ -> QC.property False
      Right (opt,sol) ->
         QC.forAll (QC.choose (0,1)) $ \lambda ->
            let val = TestLP.scalarProduct obj $
                        TestLP.affineCombination lambda sol (Array.map fromIntegral origin)
            in 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 = LP.simplex method bounds constrs dirObjA in
   let solB = LP.simplex method bounds constrs dirObjB in
   QC.counterexample (show (mapRight fst solA, mapRight fst solB)) $
   case (solA, solB) of
      (Right _, Right _) -> True
      (Left _, Left _) -> True
      _ -> False
:}
prop> :{
   forAllMethod $ \method ->
   TestLP.forAllOrigin $ \origin ->
   TestLP.forAllProblem origin $ \bounds constrs ->
   QC.forAll (TestLP.genObjective origin) $ \(_dir,obj) ->
   case (LP.simplex method bounds constrs (LP.Minimize,obj),
         LP.simplex method bounds constrs (LP.Maximize,obj)) of
      (Right (optMin,_), Right (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 (LP.simplex method bounds allConstrs (dir,obj),
         LP.simplex method bounds someConstrs (dir,obj)) of
      (Right (optAll,_), Right (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 (LP.simplex methodA bounds constrs dirObj,
         LP.simplex methodB bounds constrs dirObj) of
      (Right (optA,_), Right (optB,_)) ->
         QC.counterexample (show (optA, optB)) $
         approxReal 0.01 optA optB
      _ -> QC.property False
:}
-}
simplex ::
   (Coefficient a, Shape.Indexed sh, Shape.Index sh ~ ix) =>
   Method -> Bounds ix -> Constraints a ix ->
   (Direction, Objective sh) -> Result sh
simplex :: forall a sh ix.
(Coefficient a, Indexed sh, Index sh ~ ix) =>
Method
-> Bounds ix
-> Constraints a ix
-> (Direction, Objective sh)
-> Result sh
simplex Method
method Bounds ix
bounds Constraints a ix
constrs (Direction
dir,Objective sh
obj) =
   IO (Result sh) -> Result sh
forall a. IO a -> a
unsafePerformIO (IO (Result sh) -> Result sh) -> IO (Result sh) -> Result sh
forall a b. (a -> b) -> a -> b
$
   IO (Ptr Simplex)
-> (Ptr Simplex -> IO ())
-> (Ptr Simplex -> IO (Result sh))
-> IO (Result sh)
forall a b c. IO a -> (a -> IO b) -> (a -> IO c) -> IO c
bracket IO (Ptr Simplex)
FFI.newModel Ptr Simplex -> IO ()
FFI.deleteModel ((Ptr Simplex -> IO (Result sh)) -> IO (Result sh))
-> (Ptr Simplex -> IO (Result sh)) -> IO (Result sh)
forall a b. (a -> b) -> a -> b
$ \Ptr Simplex
lp -> do

   Ptr Simplex -> IO ()
Debug.initLog Ptr Simplex
lp
   let shape :: sh
shape = Objective sh -> sh
forall sh a. Array sh a -> sh
Array.shape Objective sh
obj
   ContT () IO () -> IO ()
forall a. ContT a IO a -> IO a
runContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CDouble
objPtr <- Array sh CDouble -> ContT () IO (Ptr CDouble)
forall sh a r. Array sh a -> ContT r IO (Ptr a)
withBuffer (Array sh CDouble -> ContT () IO (Ptr CDouble))
-> Array sh CDouble -> ContT () 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 () 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 () IO (Ptr CDouble, Ptr CDouble))
-> (Array sh CDouble, Array sh CDouble)
-> ContT () 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 () 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 () IO (Ptr CDouble, Ptr CDouble))
-> (Array ShapeInt CDouble, Array ShapeInt CDouble)
-> ContT () IO (Ptr CDouble, Ptr CDouble)
forall a b. (a -> b) -> a -> b
$ Constraints a ix
-> (Array ShapeInt CDouble, Array ShapeInt CDouble)
forall ix.
Bounds ix -> (Array ShapeInt CDouble, Array ShapeInt CDouble)
prepareRowBoundsArrays Constraints a ix
constrs
      sh
-> Constraints a ix
-> Ptr Simplex
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> ContT () IO ()
forall a sh ix.
(Coefficient a, Indexed sh, Index sh ~ ix) =>
sh
-> Constraints a ix
-> Ptr Simplex
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> ContT () IO ()
forall sh ix.
(Indexed sh, Index sh ~ ix) =>
sh
-> Constraints a ix
-> Ptr Simplex
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> ContT () IO ()
loadProblem sh
shape Constraints a ix
constrs Ptr Simplex
lp Ptr CDouble
collbPtr Ptr CDouble
colubPtr Ptr CDouble
objPtr Ptr CDouble
rowlbPtr Ptr CDouble
rowubPtr
   Ptr Simplex -> Direction -> IO ()
setOptimizationDirection Ptr Simplex
lp Direction
dir
   Method -> Ptr Simplex -> IO ()
runMethod Method
method Ptr Simplex
lp
   sh -> Ptr Simplex -> IO (Result sh)
forall sh. C sh => sh -> Ptr Simplex -> IO (Result sh)
examineStatus sh
shape Ptr Simplex
lp