{-# 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, 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.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 (Eq, 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 shape constrs lp collbPtr colubPtr objPtr rowlbPtr rowubPtr = do (coefficientsShape, coefficientsPtr, indexPtr, startPtr) <- storeConstraints shape constrs let createMatrix = FFI.newCoinPackedMatrix false (fromIntegral $ Shape.size shape) (fromIntegral $ length constrs) (fromIntegral $ Shape.size coefficientsShape) coefficientsPtr indexPtr startPtr nullPtr matrix <- MC.ContT $ bracket createMatrix FFI.deleteCoinPackedMatrix liftIO $ FFI.loadProblemFromCoinMatrix lp matrix collbPtr colubPtr objPtr rowlbPtr rowubPtr nullPtr instance Coefficient PlusMinusOne where loadProblem shape constrs lp collbPtr colubPtr objPtr rowlbPtr rowubPtr = do let shapeOffset = Shape.offset shape let coefficients = map (\(Inequality terms _bnd) -> ListHT.partition (\(Term c _) -> c == PlusOne) terms) constrs indexPtr <- withBuffer $ Array.vectorFromList $ concatMap (\(positive,negative) -> map fromIntegral $ map (\(Term _ ix) -> shapeOffset ix) positive ++ map (\(Term _ ix) -> shapeOffset ix) negative) coefficients let rowStarts = scanl (+) 0 $ map (\(Inequality terms _bnd) -> length terms) constrs startPositivePtr <- withBuffer $ Array.vectorFromList $ map fromIntegral rowStarts startNegativePtr <- withBuffer $ Array.vectorFromList $ zipWith (\k (pos,_neg) -> fromIntegral $ k + length pos) rowStarts coefficients let createMatrix = FFI.newPlusMinusOneMatrix (fromIntegral $ length constrs) (fromIntegral $ Shape.size shape) (toEnum $ fromEnum False) indexPtr startPositivePtr startNegativePtr matrix <- MC.ContT $ bracket createMatrix FFI.deletePlusMinusOneMatrix liftIO $ FFI.loadProblemFromMatrix lp matrix collbPtr colubPtr objPtr rowlbPtr rowubPtr 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)) 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 (Right (absol,(posResult,negResult))) -> QC.property (absol>=0) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False 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 (Right (absol,(posResult,negResult))) -> QC.counterexample (show(absol,(posResult,negResult))) $ QC.property (approxReal 0.001 absol (abs target)) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property 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 (dir,obj) of Right _ -> 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 (dir,obj) of Right (_,sol) -> TestLP.checkFeasibility 0.1 bounds constrs sol; _ -> QC.property 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 (dir,obj) of Right (_,sol) -> QC.forAll (QC.choose (0,1)) $ \lambda -> TestLP.checkFeasibility 0.1 bounds constrs $ TestLP.affineCombination lambda sol (Array.map fromIntegral origin); _ -> QC.property 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 (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 -} simplex :: (Coefficient a, Shape.Indexed sh, Shape.Index sh ~ ix) => Method -> Bounds ix -> Constraints a ix -> (Direction, Objective sh) -> Result sh simplex method bounds constrs (dir,obj) = unsafePerformIO $ bracket FFI.newModel FFI.deleteModel $ \lp -> do Debug.initLog lp let shape = Array.shape obj runContT $ do objPtr <- withBuffer $ Array.map realToFrac obj (collbPtr,colubPtr) <- storeBounds $ prepareColumnBoundsArrays shape bounds (rowlbPtr,rowubPtr) <- storeBounds $ prepareRowBoundsArrays constrs loadProblem shape constrs lp collbPtr colubPtr objPtr rowlbPtr rowubPtr setOptimizationDirection lp dir runMethod method lp examineStatus shape lp