diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,27 @@ +Copyright (c) Henning Thielemann 2023 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +SUCH DAMAGE. diff --git a/Makefile b/Makefile new file mode 100644 --- /dev/null +++ b/Makefile @@ -0,0 +1,12 @@ +run-test: update-test + runhaskell Setup configure --user --enable-tests + runhaskell Setup build + runhaskell Setup haddock + runhaskell Setup test highs-lp-test --show-details=streaming + + runhaskell Setup configure --user -fdebug + runhaskell Setup build + +update-test: + doctest-extract-0.1 -i src/ -o test/ --executable-main=Main.hs \ + Numeric.HiGHS.LP Numeric.HiGHS.LP.Monad diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/highs-lp.cabal b/highs-lp.cabal new file mode 100644 --- /dev/null +++ b/highs-lp.cabal @@ -0,0 +1,106 @@ +Cabal-Version: 2.2 +Name: highs-lp +Version: 0.0 +License: BSD-3-Clause +License-File: LICENSE +Author: Henning Thielemann +Maintainer: Henning Thielemann +Category: Math +Tested-With: GHC ==8.6.5 +Build-Type: Simple +Synopsis: Linear Programming using HiGHS and comfort-array +Description: + Simple interface to linear programming functions provided by HiGHS + using the flexible Array shape framework from @comfort-array@. + . + E.g. you can use @Shape.Tuple@ to convert safely + between nested tuples and arrays with the same number of elements. + . + > type X = Shape.Element + > type PairShape = Shape.NestedTuple Shape.TupleIndex (X,X) + > + > case Shape.indexTupleFromShape (Shape.static :: PairShape) of + > (posIx,negIx) -> + > case fmap (mapSnd Array.toTuple) $ snd $ + > LP.solve LP.choose [] [[1.*posIx, (-1).*negIx] ==. 314] + > (LP.Minimize, + > Array.fromTuple (23,42) :: Array PairShape Double) + > of + > Just (absol, (pos, neg)) -> + > printf "absol %f, pos %f, neg %f\n" absol pos neg + > _ -> fail "HiGHS solver failed" + . + Alternatives: @coinor-clp@, @comfort-glpk@, @hmatrix-glpk@, @glpk-hs@ + +Extra-Source-Files: + Makefile + +Flag debug + Description: Enable debug output + Default: False + Manual: True + +Source-Repository this + Tag: 0.0 + Type: darcs + Location: https://hub.darcs.net/thielema/highs-lp/ + +Source-Repository head + Type: darcs + Location: https://hub.darcs.net/thielema/highs-lp/ + +Library + Build-Depends: + linear-programming >=0.0 && <0.1, + comfort-array >=0.5.4 && <0.6, + guarded-allocation >=0.0.1 && <0.1, + QuickCheck >=2 && <3, + storable-record >=0.0.5 && <0.1, + deepseq >=1.3 && <1.6, + transformers >=0.3 && <0.7, + non-empty >=0.3.2 && <0.4, + utility-ht >=0.0.16 && <0.1, + base >=4.5 && <5 + + GHC-Options: -Wall + Default-Language: Haskell98 + Hs-Source-Dirs: src + If flag(debug) + Hs-Source-Dirs: src/debug-on + Else + Hs-Source-Dirs: src/debug-off + Exposed-Modules: + Numeric.HiGHS.LP + Numeric.HiGHS.LP.Monad + Numeric.HiGHS.LP.Example + Other-Modules: + Numeric.HiGHS.LP.Enumeration + Numeric.HiGHS.LP.FFI + Numeric.HiGHS.LP.Private + Numeric.HiGHS.LP.Debug + + PkgConfig-Depends: highs + +Test-Suite highs-lp-test + Type: exitcode-stdio-1.0 + Build-Depends: + highs-lp, + linear-programming >=0.0.1, + comfort-array, + transformers, + non-empty, + utility-ht, + doctest-exitcode-stdio >=0.0 && <0.1, + doctest-lib >=0.1 && <0.2, + QuickCheck >=2.1 && <3, + random >=1.0 && <1.3, + base >=4.5 && <5 + + GHC-Options: -Wall + Hs-Source-Dirs: test + Default-Language: Haskell98 + Main-Is: Main.hs + Other-Modules: + Test.Numeric.HiGHS.LP + Test.Numeric.HiGHS.LP.Monad + Test.Numeric.HiGHS.LP.Utility diff --git a/src/Numeric/HiGHS/LP.hs b/src/Numeric/HiGHS/LP.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/HiGHS/LP.hs @@ -0,0 +1,323 @@ +{-# 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 = solveWith 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 query method bounds constrs (dir,obj) = + unsafePerformIO $ + let shape = Array.shape obj in + let numCols = Shape.size shape in + let numRows = length constrs in + + runContT $ do + objPtr <- withBuffer $ Array.map realToFrac obj + (collbPtr,colubPtr) <- + storeBounds $ prepareColumnBoundsArrays shape bounds + (rowlbPtr,rowubPtr) <- storeBounds $ prepareRowBoundsArrays constrs + + let (coefficients, indices, rowStarts) = prepareConstraints shape constrs + (coefficientsPtr, indexPtr, startPtr) + <- storeConstraints (coefficients, indices, rowStarts) + + liftIO $ bracket Highs.create Highs.destroy $ \model -> do + Debug.initLog model + setMethod model method + checkStatus $ Highs.passLp model + (fromIntegral numCols) + (fromIntegral numRows) + (fromIntegral $ Shape.size $ Array.shape coefficients) + Highs.matrixFormatRowwise + (objectiveSense dir) + 0 objPtr + collbPtr colubPtr + rowlbPtr rowubPtr + startPtr indexPtr coefficientsPtr + + examineStatus query shape model =<< Highs.run model diff --git a/src/Numeric/HiGHS/LP/Enumeration.hsc b/src/Numeric/HiGHS/LP/Enumeration.hsc new file mode 100644 --- /dev/null +++ b/src/Numeric/HiGHS/LP/Enumeration.hsc @@ -0,0 +1,69 @@ +module Numeric.HiGHS.LP.Enumeration where + +import qualified Numeric.HiGHS.LP.FFI as Highs + +#include "interfaces/highs_c_api.h" + + +data ModelStatus = + ModelStatusNotset + | ModelStatusLoadError + | ModelStatusModelError + | ModelStatusPresolveError + | ModelStatusSolveError + | ModelStatusPostsolveError + | ModelStatusModelEmpty + | ModelStatusOptimal + | ModelStatusInfeasible + | ModelStatusUnboundedOrInfeasible + | ModelStatusUnbounded + | ModelStatusObjectiveBound + | ModelStatusObjectiveTarget + | ModelStatusTimeLimit + | ModelStatusIterationLimit + | ModelStatusUnknown + | ModelStatusSolutionLimit + | ModelStatusInterrupt + deriving (Eq, Show) + + +modelStatusFromC :: Highs.ModelStatus -> ModelStatus +modelStatusFromC (Highs.ModelStatus status) = + case status of + (#const kHighsModelStatusNotset) -> + ModelStatusNotset + (#const kHighsModelStatusLoadError) -> + ModelStatusLoadError + (#const kHighsModelStatusModelError) -> + ModelStatusModelError + (#const kHighsModelStatusPresolveError) -> + ModelStatusPresolveError + (#const kHighsModelStatusSolveError) -> + ModelStatusSolveError + (#const kHighsModelStatusPostsolveError) -> + ModelStatusPostsolveError + (#const kHighsModelStatusModelEmpty) -> + ModelStatusModelEmpty + (#const kHighsModelStatusOptimal) -> + ModelStatusOptimal + (#const kHighsModelStatusInfeasible) -> + ModelStatusInfeasible + (#const kHighsModelStatusUnboundedOrInfeasible) -> + ModelStatusUnboundedOrInfeasible + (#const kHighsModelStatusUnbounded) -> + ModelStatusUnbounded + (#const kHighsModelStatusObjectiveBound) -> + ModelStatusObjectiveBound + (#const kHighsModelStatusObjectiveTarget) -> + ModelStatusObjectiveTarget + (#const kHighsModelStatusTimeLimit) -> + ModelStatusTimeLimit + (#const kHighsModelStatusIterationLimit) -> + ModelStatusIterationLimit + (#const kHighsModelStatusUnknown) -> + ModelStatusUnknown + (#const kHighsModelStatusSolutionLimit) -> + ModelStatusSolutionLimit + (#const kHighsModelStatusInterrupt) -> + ModelStatusInterrupt + _ -> error $ "unknown model status " ++ show status diff --git a/src/Numeric/HiGHS/LP/Example.hs b/src/Numeric/HiGHS/LP/Example.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/HiGHS/LP/Example.hs @@ -0,0 +1,30 @@ +module Numeric.HiGHS.LP.Example + {-# WARNING "Only for documentation. Do not import." #-} + where + +import qualified Numeric.HiGHS.LP as LP +import Numeric.HiGHS.LP ((.*), (==.)) + +import qualified Data.Array.Comfort.Storable as Array +import qualified Data.Array.Comfort.Shape as Shape +import Data.Array.Comfort.Storable (Array) +import Data.Tuple.HT (mapSnd) + +import Text.Printf (printf) + + +type X = Shape.Element +type PairShape = Shape.NestedTuple Shape.TupleIndex (X,X) + +main :: IO () +main = + case Shape.indexTupleFromShape (Shape.static :: PairShape) of + (posIx,negIx) -> + case fmap (mapSnd Array.toTuple) $ snd $ + LP.solve LP.choose [] [[1.*posIx, (-1).*negIx] ==. 314] + (LP.Minimize, + Array.fromTuple (23,42) :: Array PairShape Double) + of + Just (absol, (pos, neg)) -> + printf "absol %f, pos %f, neg %f\n" absol pos neg + _ -> fail "HiGHS solver failed" diff --git a/src/Numeric/HiGHS/LP/FFI.hsc b/src/Numeric/HiGHS/LP/FFI.hsc new file mode 100644 --- /dev/null +++ b/src/Numeric/HiGHS/LP/FFI.hsc @@ -0,0 +1,249 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +module Numeric.HiGHS.LP.FFI where + +import qualified Foreign.Storable.Newtype as Store +import qualified Foreign.C.Types as C +import qualified Foreign.Storable as St +import Foreign.C.String (CString) +import Foreign.Storable (Storable) +import Foreign.Ptr (Ptr) + +import Data.Int (Int32) + +import Prelude hiding (Bool) + + +type CDouble = C.CDouble +type CInt = C.CInt + + +#include "interfaces/highs_c_api.h" + +type HighsInt = #{type HighsInt} + +newtype ObjSense = ObjSense HighsInt + +objSenseMinimize, objSenseMaximize :: ObjSense +objSenseMinimize = ObjSense (#const kHighsObjSenseMinimize) +objSenseMaximize = ObjSense (#const kHighsObjSenseMaximize) + +newtype MatrixFormat = MatrixFormat HighsInt + +matrixFormatColwise, matrixFormatRowwise :: MatrixFormat +matrixFormatColwise = MatrixFormat (#const kHighsMatrixFormatColwise) +matrixFormatRowwise = MatrixFormat (#const kHighsMatrixFormatRowwise) + +newtype Status = Status HighsInt + deriving (Eq, Show) + +statusError, statusOk, statusWarning :: Status +statusError = Status (#const kHighsStatusError) +statusOk = Status (#const kHighsStatusOk) +statusWarning = Status (#const kHighsStatusWarning) + +newtype ModelStatus = ModelStatus {deModelStatus :: HighsInt} + deriving (Eq, Show) + +modelStatusNotset :: ModelStatus +modelStatusLoadError :: ModelStatus +modelStatusModelError :: ModelStatus +modelStatusPresolveError :: ModelStatus +modelStatusSolveError :: ModelStatus +modelStatusPostsolveError :: ModelStatus +modelStatusModelEmpty :: ModelStatus +modelStatusOptimal :: ModelStatus +modelStatusInfeasible :: ModelStatus +modelStatusUnboundedOrInfeasible :: ModelStatus +modelStatusUnbounded :: ModelStatus +modelStatusObjectiveBound :: ModelStatus +modelStatusObjectiveTarget :: ModelStatus +modelStatusTimeLimit :: ModelStatus +modelStatusIterationLimit :: ModelStatus +modelStatusUnknown :: ModelStatus +modelStatusSolutionLimit :: ModelStatus +modelStatusInterrupt :: ModelStatus + +modelStatusNotset + = ModelStatus (#const kHighsModelStatusNotset) +modelStatusLoadError + = ModelStatus (#const kHighsModelStatusLoadError) +modelStatusModelError + = ModelStatus (#const kHighsModelStatusModelError) +modelStatusPresolveError + = ModelStatus (#const kHighsModelStatusPresolveError) +modelStatusSolveError + = ModelStatus (#const kHighsModelStatusSolveError) +modelStatusPostsolveError + = ModelStatus (#const kHighsModelStatusPostsolveError) +modelStatusModelEmpty + = ModelStatus (#const kHighsModelStatusModelEmpty) +modelStatusOptimal + = ModelStatus (#const kHighsModelStatusOptimal) +modelStatusInfeasible + = ModelStatus (#const kHighsModelStatusInfeasible) +modelStatusUnboundedOrInfeasible + = ModelStatus (#const kHighsModelStatusUnboundedOrInfeasible) +modelStatusUnbounded + = ModelStatus (#const kHighsModelStatusUnbounded) +modelStatusObjectiveBound + = ModelStatus (#const kHighsModelStatusObjectiveBound) +modelStatusObjectiveTarget + = ModelStatus (#const kHighsModelStatusObjectiveTarget) +modelStatusTimeLimit + = ModelStatus (#const kHighsModelStatusTimeLimit) +modelStatusIterationLimit + = ModelStatus (#const kHighsModelStatusIterationLimit) +modelStatusUnknown + = ModelStatus (#const kHighsModelStatusUnknown) +modelStatusSolutionLimit + = ModelStatus (#const kHighsModelStatusSolutionLimit) +modelStatusInterrupt + = ModelStatus (#const kHighsModelStatusInterrupt) + +instance Storable ModelStatus where + sizeOf = Store.sizeOf deModelStatus + alignment = Store.alignment deModelStatus + peek = Store.peek ModelStatus + poke = Store.poke $ deModelStatus + +newtype BasisStatus = BasisStatus {deBasisStatus :: HighsInt} + +basisStatusLower, basisStatusBasic, basisStatusUpper, + basisStatusZero, basisStatusNonbasic :: BasisStatus + +basisStatusLower = BasisStatus (#const kHighsBasisStatusLower) +basisStatusBasic = BasisStatus (#const kHighsBasisStatusBasic) +basisStatusUpper = BasisStatus (#const kHighsBasisStatusUpper) +basisStatusZero = BasisStatus (#const kHighsBasisStatusZero) +basisStatusNonbasic = BasisStatus (#const kHighsBasisStatusNonbasic) + +instance Show BasisStatus where + show (BasisStatus status) = + case status of + (#const kHighsBasisStatusLower) -> "Highs.basisStatusLower" + (#const kHighsBasisStatusBasic) -> "Highs.basisStatusBasic" + (#const kHighsBasisStatusUpper) -> "Highs.basisStatusUpper" + (#const kHighsBasisStatusZero) -> "Highs.basisStatusZero" + (#const kHighsBasisStatusNonbasic) -> "Highs.basisStatusNonbasic" + _ -> "(BasisStatus " ++ show status ++ ")" + +instance Storable BasisStatus where + sizeOf = Store.sizeOf deBasisStatus + alignment = Store.alignment deBasisStatus + peek = Store.peek BasisStatus + poke = Store.poke deBasisStatus + + +foreign import ccall "Highs_lpCall" + lpCall :: + HighsInt -> HighsInt -> + HighsInt -> MatrixFormat -> + ObjSense -> + CDouble -> Ptr CDouble -> + Ptr CDouble -> Ptr CDouble -> + Ptr CDouble -> Ptr CDouble -> + Ptr HighsInt -> Ptr HighsInt -> Ptr CDouble -> + Ptr CDouble -> Ptr CDouble -> + Ptr CDouble -> Ptr CDouble -> + Ptr BasisStatus -> + Ptr BasisStatus -> + Ptr ModelStatus -> + IO Status + + +data Highs = Highs + +foreign import ccall "Highs_create" + create :: IO (Ptr Highs) + +foreign import ccall "Highs_destroy" + destroy :: Ptr Highs -> IO () + +foreign import ccall "Highs_run" + run :: Ptr Highs -> IO Status + +foreign import ccall "Highs_passLp" + passLp :: + 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 + +foreign import ccall "Highs_getSolution" + getSolution :: + Ptr Highs -> + Ptr CDouble -> Ptr CDouble -> + Ptr CDouble -> Ptr CDouble -> + IO Status + +foreign import ccall "Highs_getBasis" + getBasis :: + Ptr Highs -> + Ptr BasisStatus -> + Ptr BasisStatus -> + IO Status + +foreign import ccall "Highs_getModelStatus" + getModelStatus :: Ptr Highs -> IO ModelStatus + +foreign import ccall "Highs_getObjectiveValue" + getObjectiveValue :: Ptr Highs -> IO CDouble + + +foreign import ccall "Highs_getNumCol" + getNumCol :: Ptr highs -> IO HighsInt + +foreign import ccall "Highs_getNumRow" + getNumRow :: Ptr highs -> IO HighsInt + + +foreign import ccall "Highs_addRows" + addRows :: + Ptr Highs -> + HighsInt -> + Ptr CDouble -> Ptr CDouble -> + HighsInt -> + Ptr HighsInt -> Ptr HighsInt -> + Ptr CDouble -> + IO Status + +foreign import ccall "Highs_changeObjectiveSense" + changeObjectiveSense :: Ptr Highs -> ObjSense -> IO Status + +foreign import ccall "Highs_changeObjectiveOffset" + changeObjectiveOffset :: Ptr Highs -> CDouble -> IO Status + +foreign import ccall "Highs_changeColsCostByRange" + changeColsCostByRange :: + Ptr Highs -> HighsInt -> HighsInt -> Ptr CDouble -> IO Status + + +foreign import ccall "Highs_readModel" + readModel :: Ptr Highs -> CString -> IO Status + +foreign import ccall "Highs_writeModel" + writeModel :: Ptr Highs -> CString -> IO Status + + +newtype Bool = Bool HighsInt + +false, true :: Bool +false = Bool 0 +true = Bool 1 + +foreign import ccall "Highs_setBoolOptionValue" + setBoolOptionValue :: Ptr Highs -> CString -> Bool -> IO Status + +foreign import ccall "Highs_setIntOptionValue" + setIntOptionValue :: Ptr Highs -> CString -> HighsInt -> IO Status + +foreign import ccall "Highs_setDoubleOptionValue" + setDoubleOptionValue :: Ptr Highs -> CString -> CDouble -> IO Status + +foreign import ccall "Highs_setStringOptionValue" + setStringOptionValue :: Ptr Highs -> CString -> CString -> IO Status diff --git a/src/Numeric/HiGHS/LP/Monad.hs b/src/Numeric/HiGHS/LP/Monad.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/HiGHS/LP/Monad.hs @@ -0,0 +1,284 @@ +{-# 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 diff --git a/src/Numeric/HiGHS/LP/Private.hs b/src/Numeric/HiGHS/LP/Private.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/HiGHS/LP/Private.hs @@ -0,0 +1,248 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +module Numeric.HiGHS.LP.Private where + +import qualified Numeric.HiGHS.LP.FFI as Highs +import Numeric.HiGHS.LP.Enumeration (ModelStatus, modelStatusFromC) +import Numeric.HiGHS.LP.FFI (Highs, HighsInt) +import Numeric.LinearProgramming.Common + (Term(..), Bound(..), Inequality(Inequality), + Bounds, Constraints, Direction(..)) + +import qualified Data.Array.Comfort.Boxed as BoxedArray +import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayMonadic +import qualified Data.Array.Comfort.Storable.Unchecked as ArrayUnchecked +import qualified Data.Array.Comfort.Storable as Array +import qualified Data.Array.Comfort.Shape as Shape +import qualified Data.Traversable as Trav +import Data.Array.Comfort.Storable (Array) +import Data.Foldable (for_) +import Data.Tuple.HT (mapPair) + +import qualified Control.Monad.Trans.Reader as MR +import qualified Control.Monad.Trans.Cont as MC +import qualified Control.Applicative.HT as AppHT +import qualified Control.Functor.HT as FuncHT +import Control.Monad (guard, when) +import Control.Applicative (liftA2, liftA3) + +import qualified Foreign.Marshal.Array.Guarded as ForeignArray +import Foreign.Storable (pokeElemOff, peekElemOff) +import Foreign.ForeignPtr (withForeignPtr) +import Foreign.Ptr (Ptr) +import Foreign.C.String (withCString) +import Foreign.C.Types (CDouble) + + + +withBuffer :: Array sh a -> MC.ContT r IO (Ptr a) +withBuffer arr = + MC.ContT $ withForeignPtr (ArrayUnchecked.buffer arr) + +runContT :: MC.ContT a IO a -> IO a +runContT act = MC.runContT act return + + + +positiveInfinity, negativeInfinity :: CDouble +positiveInfinity = 1/0 +negativeInfinity = -1/0 + +prepareBounds :: Inequality a -> (a, (CDouble, CDouble)) +prepareBounds (Inequality x bnd) = + (,) x $ + case bnd of + LessEqual up -> (negativeInfinity, realToFrac up) + GreaterEqual lo -> (realToFrac lo, positiveInfinity) + Between lo up -> (realToFrac lo, realToFrac up) + Equal y -> (realToFrac y, realToFrac y) + Free -> (negativeInfinity, positiveInfinity) + +prepareColumnBoundsArrays :: + (Shape.Indexed sh, Shape.Index sh ~ ix) => + sh -> Bounds ix -> (Array sh CDouble, Array sh CDouble) +prepareColumnBoundsArrays shape = + mapPair (Array.fromBoxed, Array.fromBoxed) . + FuncHT.unzip . + BoxedArray.fromAssociations (0, positiveInfinity) shape . + map prepareBounds + + +type ShapeInt = Shape.ZeroBased Int + +prepareRowBoundsArrays :: + Bounds ix -> (Array ShapeInt CDouble, Array ShapeInt CDouble) +prepareRowBoundsArrays constrs = + let shape = Shape.ZeroBased $ length constrs in + mapPair (Array.fromList shape, Array.fromList shape) $ + unzip $ map (snd . prepareBounds) constrs + +storeBounds :: + (Array sh CDouble, Array sh CDouble) -> + MC.ContT r IO (Ptr CDouble, Ptr CDouble) +storeBounds = AppHT.mapPair (withBuffer, withBuffer) + + +prepareConstraints :: + (Shape.Indexed sh, Shape.Index sh ~ ix) => + sh -> Constraints Double ix -> + (Array ShapeInt CDouble, Array ShapeInt HighsInt, Array ShapeInt HighsInt) +prepareConstraints shape constrs = + let {- + Highs.passLp returns Warning when there are zero coefficients. + I think zero coefficients are reasonably ok. + -} + constrsNonZero = map (fmap (filter (\(Term c _x) -> c/=0))) constrs + rowStarts = + Array.vectorFromList $ scanl (+) 0 $ + map (\(Inequality terms _bnd) -> fromIntegral $ length terms) + constrsNonZero + shapeOffset = Shape.offset shape + coefficients = + concatMap (\(Inequality terms _bnd) -> terms) constrsNonZero + indexArr = + Array.vectorFromList $ + map (\(Term _ ix) -> fromIntegral $ shapeOffset ix) coefficients + coefficientArr = + Array.vectorFromList $ + map (\(Term c _) -> realToFrac c) coefficients + in (coefficientArr, indexArr, rowStarts) + +storeConstraints :: + (Array ShapeInt CDouble, Array ShapeInt HighsInt, Array ShapeInt HighsInt) -> + MC.ContT r IO (Ptr CDouble, Ptr HighsInt, Ptr HighsInt) +storeConstraints (coefficients, indices, rowStarts) = + liftA3 (,,) + (withBuffer coefficients) + (withBuffer indices) + (withBuffer rowStarts) + + +objectiveSense :: Direction -> Highs.ObjSense +objectiveSense dir = + case dir of + Minimize -> Highs.objSenseMinimize + Maximize -> Highs.objSenseMaximize + + +setBoolOptionValue :: Ptr Highs -> String -> Highs.Bool -> IO () +setBoolOptionValue model key b = + checkStatus $ withCString key $ \cstr -> + Highs.setBoolOptionValue model cstr b + +setIntOptionValue :: Ptr Highs -> String -> HighsInt -> IO () +setIntOptionValue model key n = + checkStatus $ withCString key $ \cstr -> + Highs.setIntOptionValue model cstr n + +setDoubleOptionValue :: Ptr Highs -> String -> CDouble -> IO () +setDoubleOptionValue model key x = + checkStatus $ withCString key $ \cstr -> + Highs.setDoubleOptionValue model cstr x + +setStringOptionValue :: Ptr Highs -> String -> String -> IO () +setStringOptionValue model key value = + checkStatus $ + withCString key $ \keyPtr -> + withCString value $ \valuePtr -> + Highs.setStringOptionValue model keyPtr valuePtr + + + +newtype Method = Method {deMethod :: String} + deriving (Show) + +simplex, choose, ipm :: Method +simplex = Method "simplex" +choose = Method "choose" +ipm = Method "ipm" + +setMethod :: Ptr Highs -> Method -> IO () +setMethod model (Method method) = + setStringOptionValue model "solver" method + +checkStatus :: IO Highs.Status -> IO () +checkStatus act = do + status <- act + when (status == Highs.statusError) $ fail "Highs function failed" + + +newtype Query sh a = Query (MR.ReaderT (sh, Ptr Highs) IO a) + deriving (Functor, Applicative, Monad) + +getResult :: (Shape.C sh) => Query sh (Double, Array sh Double) +getResult = liftA2 (,) getObjectiveValue getOptimalVector + +getObjectiveValue :: Query sh Double +getObjectiveValue = + Query $ MR.ReaderT $ fmap realToFrac . Highs.getObjectiveValue . snd + +doubleFromCDoubleBuffer :: Int -> Ptr CDouble -> Ptr Double -> IO () +doubleFromCDoubleBuffer n srcPtr dstPtr = + for_ (take n [0..]) $ \k -> + pokeElemOff dstPtr k . realToFrac =<< peekElemOff srcPtr k + +getOptimalVector :: (Shape.C sh) => Query sh (Array sh Double) +getOptimalVector = + Query $ MR.ReaderT $ \(shape,model) -> + ArrayMonadic.unsafeCreateWithSize shape $ \numCols arrPtr -> + (fmap fromIntegral (Highs.getNumRow model) >>= ) $ \numRows -> + ForeignArray.alloca numCols $ \colValuePtr -> + ForeignArray.alloca numCols $ \colDualPtr -> + ForeignArray.alloca numRows $ \rowValuePtr -> + ForeignArray.alloca numRows $ \rowDualPtr -> do + + checkStatus $ + Highs.getSolution model colValuePtr colDualPtr rowValuePtr rowDualPtr + + doubleFromCDoubleBuffer numCols colValuePtr arrPtr + +getSolutionVectors :: + (Shape.C sh) => + Query sh + ((Array sh Double, Array sh Double), + (Array ShapeInt Double, Array ShapeInt Double)) +getSolutionVectors = + Query $ MR.ReaderT $ \(shape,model) -> + fmap (\(colPrimal,(colDual,row2)) -> ((colPrimal,colDual),row2)) $ + ArrayMonadic.unsafeCreateWithSizeAndResult shape $ \numCols colValueArr -> + ArrayMonadic.unsafeCreateWithSizeAndResult shape $ \_numCols colDualArr -> + (fmap fromIntegral (Highs.getNumRow model) >>= ) $ \numRows -> + let constraintsShape = Shape.ZeroBased numRows in + ArrayMonadic.unsafeCreateWithSizeAndResult constraintsShape $ + \_numRows rowValueArr -> + ArrayMonadic.unsafeCreate constraintsShape $ \rowDualArr -> + ForeignArray.alloca numCols $ \colValuePtr -> + ForeignArray.alloca numCols $ \colDualPtr -> + ForeignArray.alloca numRows $ \rowValuePtr -> + ForeignArray.alloca numRows $ \rowDualPtr -> do + + checkStatus $ + Highs.getSolution model colValuePtr colDualPtr rowValuePtr rowDualPtr + + doubleFromCDoubleBuffer numCols colValuePtr colValueArr + doubleFromCDoubleBuffer numCols colDualPtr colDualArr + doubleFromCDoubleBuffer numRows rowValuePtr rowValueArr + doubleFromCDoubleBuffer numRows rowDualPtr rowDualArr + +getBasisStatus :: + (Shape.C sh) => + Query sh (Array sh Highs.BasisStatus, Array ShapeInt Highs.BasisStatus) +getBasisStatus = + Query $ MR.ReaderT $ \(shape,model) -> + ArrayMonadic.unsafeCreateWithSizeAndResult shape $ \_numCols colStatusPtr -> + (fmap fromIntegral (Highs.getNumRow model) >>= ) $ \numRows -> + ArrayMonadic.unsafeCreate (Shape.ZeroBased numRows) $ \rowStatusPtr -> + checkStatus $ Highs.getBasis model colStatusPtr rowStatusPtr + + + +type Result sh = (ModelStatus, Maybe (Double, Array sh Double)) + +examineStatus :: (Shape.C sh) => + Query sh a -> sh -> Ptr Highs -> Highs.Status -> IO (ModelStatus, Maybe a) +examineStatus (Query query) shape model status = + liftA2 (,) (fmap modelStatusFromC $ Highs.getModelStatus model) $ + Trav.for (guard (status /= Highs.statusError)) $ \() -> + MR.runReaderT query (shape,model) diff --git a/src/debug-off/Numeric/HiGHS/LP/Debug.hs b/src/debug-off/Numeric/HiGHS/LP/Debug.hs new file mode 100644 --- /dev/null +++ b/src/debug-off/Numeric/HiGHS/LP/Debug.hs @@ -0,0 +1,9 @@ +module Numeric.HiGHS.LP.Debug where + +import qualified Numeric.HiGHS.LP.FFI as FFI +import Numeric.HiGHS.LP.Private (setBoolOptionValue) +import Foreign.Ptr (Ptr) + + +initLog :: Ptr FFI.Highs -> IO () +initLog model = setBoolOptionValue model "output_flag" FFI.false diff --git a/src/debug-on/Numeric/HiGHS/LP/Debug.hs b/src/debug-on/Numeric/HiGHS/LP/Debug.hs new file mode 100644 --- /dev/null +++ b/src/debug-on/Numeric/HiGHS/LP/Debug.hs @@ -0,0 +1,9 @@ +module Numeric.HiGHS.LP.Debug where + +import qualified Numeric.HiGHS.LP.FFI as FFI +import Numeric.HiGHS.LP.Private (setBoolOptionValue) +import Foreign.Ptr (Ptr) + + +initLog :: Ptr FFI.Highs -> IO () +initLog model = setBoolOptionValue model "output_flag" FFI.true diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -0,0 +1,12 @@ +-- Do not edit! Automatically created with doctest-extract. +module Main where + +import qualified Test.Numeric.HiGHS.LP +import qualified Test.Numeric.HiGHS.LP.Monad + +import qualified Test.DocTest.Driver as DocTest + +main :: IO () +main = DocTest.run $ do + Test.Numeric.HiGHS.LP.test + Test.Numeric.HiGHS.LP.Monad.test diff --git a/test/Test/Numeric/HiGHS/LP.hs b/test/Test/Numeric/HiGHS/LP.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/HiGHS/LP.hs @@ -0,0 +1,298 @@ +-- Do not edit! Automatically created with doctest-extract from src/Numeric/HiGHS/LP.hs +{-# LINE 46 "src/Numeric/HiGHS/LP.hs" #-} + +module Test.Numeric.HiGHS.LP where + +import Test.DocTest.Base +import qualified Test.DocTest.Driver as DocTest + +{-# LINE 47 "src/Numeric/HiGHS/LP.hs" #-} +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]) + +test :: DocTest.T () +test = do + DocTest.printPrefix "Numeric.HiGHS.LP:81: " +{-# LINE 81 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.example( +{-# LINE 81 "src/Numeric/HiGHS/LP.hs" #-} + + 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) + ) + [ExpectedLine [LineChunk "Just (28.0,(5.0,0.0,4.0))"]] + DocTest.printPrefix "Numeric.HiGHS.LP:92: " +{-# LINE 92 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.example( +{-# LINE 92 "src/Numeric/HiGHS/LP.hs" #-} + + 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) + ) + [ExpectedLine [LineChunk "Just (116.0,(22.0,12.0,32.0))"]] + DocTest.printPrefix "Numeric.HiGHS.LP:103: " +{-# LINE 103 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.example( +{-# LINE 103 "src/Numeric/HiGHS/LP.hs" #-} + + 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) + ) + [ExpectedLine [LineChunk "(ModelStatusInfeasible,",WildCardChunk,LineChunk ")"]] + DocTest.printPrefix "Numeric.HiGHS.LP:114: " +{-# LINE 114 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.example( +{-# LINE 114 "src/Numeric/HiGHS/LP.hs" #-} + + 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) + ) + [ExpectedLine [LineChunk "(ModelStatusUnbounded,",WildCardChunk,LineChunk ")"]] + DocTest.printPrefix "Numeric.HiGHS.LP:125: " +{-# LINE 125 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 125 "src/Numeric/HiGHS/LP.hs" #-} + + 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) + ) + DocTest.printPrefix "Numeric.HiGHS.LP:141: " +{-# LINE 141 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 141 "src/Numeric/HiGHS/LP.hs" #-} + + 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) + ) + DocTest.printPrefix "Numeric.HiGHS.LP:157: " +{-# LINE 157 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 157 "src/Numeric/HiGHS/LP.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP:166: " +{-# LINE 166 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 166 "src/Numeric/HiGHS/LP.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP:175: " +{-# LINE 175 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 175 "src/Numeric/HiGHS/LP.hs" #-} + + 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) + ) + DocTest.printPrefix "Numeric.HiGHS.LP:187: " +{-# LINE 187 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 187 "src/Numeric/HiGHS/LP.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP:204: " +{-# LINE 204 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 204 "src/Numeric/HiGHS/LP.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP:218: " +{-# LINE 218 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 218 "src/Numeric/HiGHS/LP.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP:229: " +{-# LINE 229 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 229 "src/Numeric/HiGHS/LP.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP:244: " +{-# LINE 244 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 244 "src/Numeric/HiGHS/LP.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP:265: " +{-# LINE 265 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.example( +{-# LINE 265 "src/Numeric/HiGHS/LP.hs" #-} + + 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) + ) + [ExpectedLine [LineChunk "Just (28.0,((Highs.basisStatusBasic,Highs.basisStatusLower,Highs.basisStatusBasic),[Highs.basisStatusUpper,Highs.basisStatusUpper]))"]] + DocTest.printPrefix "Numeric.HiGHS.LP:277: " +{-# LINE 277 "src/Numeric/HiGHS/LP.hs" #-} + DocTest.property( +{-# LINE 277 "src/Numeric/HiGHS/LP.hs" #-} + + 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 + ) diff --git a/test/Test/Numeric/HiGHS/LP/Monad.hs b/test/Test/Numeric/HiGHS/LP/Monad.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/HiGHS/LP/Monad.hs @@ -0,0 +1,183 @@ +-- Do not edit! Automatically created with doctest-extract from src/Numeric/HiGHS/LP/Monad.hs +{-# LINE 52 "src/Numeric/HiGHS/LP/Monad.hs" #-} + +{-# OPTIONS_GHC -XTypeFamilies #-} +{-# OPTIONS_GHC -XTypeOperators #-} +module Test.Numeric.HiGHS.LP.Monad where + +import Test.DocTest.Base +import qualified Test.DocTest.Driver as DocTest + +{-# LINE 55 "src/Numeric/HiGHS/LP/Monad.hs" #-} +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 + +test :: DocTest.T () +test = do + DocTest.printPrefix "Numeric.HiGHS.LP.Monad:198: " +{-# LINE 198 "src/Numeric/HiGHS/LP/Monad.hs" #-} + DocTest.example( +{-# LINE 198 "src/Numeric/HiGHS/LP/Monad.hs" #-} + + 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)) + ) + [ExpectedLine [LineChunk "Just (28.0,(5.0,0.0,4.0))"]] + DocTest.printPrefix "Numeric.HiGHS.LP.Monad:210: " +{-# LINE 210 "src/Numeric/HiGHS/LP/Monad.hs" #-} + DocTest.property( +{-# LINE 210 "src/Numeric/HiGHS/LP/Monad.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP.Monad:222: " +{-# LINE 222 "src/Numeric/HiGHS/LP/Monad.hs" #-} + DocTest.property( +{-# LINE 222 "src/Numeric/HiGHS/LP/Monad.hs" #-} + + 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 + ) + DocTest.printPrefix "Numeric.HiGHS.LP.Monad:235: " +{-# LINE 235 "src/Numeric/HiGHS/LP/Monad.hs" #-} + DocTest.property( +{-# LINE 235 "src/Numeric/HiGHS/LP/Monad.hs" #-} + + 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) + ) diff --git a/test/Test/Numeric/HiGHS/LP/Utility.hs b/test/Test/Numeric/HiGHS/LP/Utility.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/HiGHS/LP/Utility.hs @@ -0,0 +1,23 @@ +module Test.Numeric.HiGHS.LP.Utility where + +import Data.Traversable (Traversable, traverse) +import Data.Foldable (Foldable, traverse_) + +import qualified Control.Monad.Trans.State as MS + +import Data.Tuple.HT (double) + + +traverse_Lag :: + (Foldable t, Monad m) => + b -> (b -> a -> m b) -> t a -> m () +traverse_Lag b0 f = + flip MS.evalStateT b0 . + traverse_ (\a -> MS.StateT $ \b -> fmap double $ f b a) + +traverseLag :: + (Traversable t, Monad m) => + b -> (b -> a -> m b) -> t a -> m (t b) +traverseLag b0 f = + flip MS.evalStateT b0 . + traverse (\a -> MS.StateT $ \b -> fmap double $ f b a)