highs-lp (empty) → 0.0
raw patch · 16 files changed
+1885/−0 lines, 16 filesdep +QuickCheckdep +basedep +comfort-arraysetup-changed
Dependencies added: QuickCheck, base, comfort-array, deepseq, doctest-exitcode-stdio, doctest-lib, guarded-allocation, highs-lp, linear-programming, non-empty, random, storable-record, transformers, utility-ht
Files
- LICENSE +27/−0
- Makefile +12/−0
- Setup.lhs +3/−0
- highs-lp.cabal +106/−0
- src/Numeric/HiGHS/LP.hs +323/−0
- src/Numeric/HiGHS/LP/Enumeration.hsc +69/−0
- src/Numeric/HiGHS/LP/Example.hs +30/−0
- src/Numeric/HiGHS/LP/FFI.hsc +249/−0
- src/Numeric/HiGHS/LP/Monad.hs +284/−0
- src/Numeric/HiGHS/LP/Private.hs +248/−0
- src/debug-off/Numeric/HiGHS/LP/Debug.hs +9/−0
- src/debug-on/Numeric/HiGHS/LP/Debug.hs +9/−0
- test/Main.hs +12/−0
- test/Test/Numeric/HiGHS/LP.hs +298/−0
- test/Test/Numeric/HiGHS/LP/Monad.hs +183/−0
- test/Test/Numeric/HiGHS/LP/Utility.hs +23/−0
+ LICENSE view
@@ -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.
+ Makefile view
@@ -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
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ highs-lp.cabal view
@@ -0,0 +1,106 @@+Cabal-Version: 2.2+Name: highs-lp+Version: 0.0+License: BSD-3-Clause+License-File: LICENSE+Author: Henning Thielemann <haskell@henning-thielemann.de>+Maintainer: Henning Thielemann <haskell@henning-thielemann.de>+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
+ src/Numeric/HiGHS/LP.hs view
@@ -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
+ src/Numeric/HiGHS/LP/Enumeration.hsc view
@@ -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
+ src/Numeric/HiGHS/LP/Example.hs view
@@ -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"
+ src/Numeric/HiGHS/LP/FFI.hsc view
@@ -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
+ src/Numeric/HiGHS/LP/Monad.hs view
@@ -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
+ src/Numeric/HiGHS/LP/Private.hs view
@@ -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)
+ src/debug-off/Numeric/HiGHS/LP/Debug.hs view
@@ -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
+ src/debug-on/Numeric/HiGHS/LP/Debug.hs view
@@ -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
+ test/Main.hs view
@@ -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
+ test/Test/Numeric/HiGHS/LP.hs view
@@ -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+ )
+ test/Test/Numeric/HiGHS/LP/Monad.hs view
@@ -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)+ )
+ test/Test/Numeric/HiGHS/LP/Utility.hs view
@@ -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)