{-# 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