vertexenum 0.1.1.0 → 1.0.0.0
raw patch · 9 files changed
+210/−110 lines, 9 filesdep +extradep +monad-loggerdep +simplex-methoddep −hmatrix-glpkdep ~containers
Dependencies added: extra, monad-logger, simplex-method
Dependencies removed: hmatrix-glpk
Dependency ranges changed: containers
Files
- CHANGELOG.md +6/−0
- README.md +10/−11
- src/Geometry/VertexEnum.hs +7/−5
- src/Geometry/VertexEnum/Constraint.hs +7/−2
- src/Geometry/VertexEnum/Internal.hs +101/−50
- src/Geometry/VertexEnum/LinearCombination.hs +6/−1
- src/Geometry/VertexEnum/VertexEnum.hs +47/−27
- tests/Main.hs +10/−9
- vertexenum.cabal +16/−5
CHANGELOG.md view
@@ -1,5 +1,11 @@ # Changelog for `vertexenum` + +## 1.0.0.0 - 2024-05-08 + +The package does no longer depend on the 'hmatrix-glpk' package. + + ## 0.1.1.0 - 2023-11-20 The types `LinearCombination` and `Constraint` are parametric now.
README.md view
@@ -1,17 +1,14 @@ # vertexenum <!-- badges: start --> -[](https://github.com/stla/vertexenum/actions/workflows/Stack.yml) +[](https://github.com/stla/vertexenum/actions/workflows/Stack-lts.yml) +[](https://github.com/stla/vertexenum/actions/workflows/Stack-nightly.yml) <!-- badges: end --> -*Get the vertices of an intersection of halfspaces.* +***Get the vertices of an intersection of halfspaces.*** ____ -This package depends on the packages **hmatrix** and **hmatrix-glpk**; follow -[this link](https://github.com/haskell-numerics/hmatrix/blob/master/INSTALL.md) -for installation instructions. - Consider the following system of linear inequalities: $$\left\{\begin{matrix} -5 & \leqslant & x & \leqslant & 4 \\ -5 & \leqslant & y & \leqslant & 3-x \\ -10 & \leqslant & z & \leqslant & 6-2x-y \end{matrix}.\right.$$ @@ -21,12 +18,14 @@ polytope: ```haskell -import Data.VectorSpace ( AdditiveGroup((^+^), (^-^)) - , VectorSpace((*^)) ) +import Data.VectorSpace ( + AdditiveGroup( (^+^), (^-^) ) + , VectorSpace( (*^) ) + ) import Geometry.VertexEnum -constraints :: [Constraint Double] -constraints = +inequalities :: [Constraint Rational] +inequalities = [ x .>= (-5) -- shortcut for `x .>=. cst (-5)` , x .<= 4 , y .>= (-5) @@ -46,4 +45,4 @@ interior to the polytope. If this argument is `Nothing`, an interior point is automatically calculated. You can get it with the `interiorPoint` function. It is easy to mentally get an interior point for the above example, but in -general this is not an easy problem.+general this is not an easy problem.
src/Geometry/VertexEnum.hs view
@@ -1,21 +1,23 @@ {-| Module : Geometry.VertexEnum Description : Vertex enumeration of convex polytopes. -Copyright : (c) Stéphane Laurent, 2023 +Copyright : (c) Stéphane Laurent, 2023-2024 License : GPL-3 Maintainer : laurent_step@outlook.fr -See README for an example. +Enumeration of the vertices of a convex polytope given by linear +inequalities. See README for an example. -} module Geometry.VertexEnum ( module X ) where -import Geometry.VertexEnum.Constraint as X ( Constraint(..) - , Sense(..) +import Geometry.VertexEnum.Constraint as X ( Constraint (..) + , Sense (..) , (.>=.), (.<=.), (.>=), (.<=) ) import Geometry.VertexEnum.LinearCombination as X ( VarIndex - , LinearCombination(..) + , LinearCombination (..) + , Var , newVar , linearCombination , constant
src/Geometry/VertexEnum/Constraint.hs view
@@ -2,12 +2,13 @@ module Geometry.VertexEnum.Constraint ( Sense (..) , Constraint (..) + , toRationalConstraint , (.>=.) , (.<=.) , (.>=) , (.<=) ) where -import Geometry.VertexEnum.LinearCombination ( LinearCombination, constant ) +import Geometry.VertexEnum.LinearCombination ( LinearCombination, constant, toRationalLinearCombination ) data Sense = Gt | Lt deriving Eq @@ -19,6 +20,10 @@ data Constraint a = Constraint (LinearCombination a) Sense (LinearCombination a) +toRationalConstraint :: Real a => Constraint a -> Constraint Rational +toRationalConstraint (Constraint lhs sense rhs) = + Constraint (toRationalLinearCombination lhs) sense (toRationalLinearCombination rhs) + instance Show a => Show (Constraint a) where show :: Constraint a -> String show (Constraint lhs sense rhs) = show lhs ++ " " ++ show sense ++ " " ++ show rhs @@ -35,4 +40,4 @@ (.<=) :: LinearCombination a -> a -> Constraint a (.<=) lhs x = (.<=.) lhs (constant x) -infix 4 .<=., .>=. +infix 4 .<=., .>=., .<=, .>=
src/Geometry/VertexEnum/Internal.hs view
@@ -1,24 +1,35 @@ module Geometry.VertexEnum.Internal ( normalizeConstraints , varsOfConstraint + , feasiblePoint + , findSigns , iPoint ) where +import Prelude hiding ( EQ ) +import Control.Monad.Logger ( + runStdoutLoggingT + , filterLogger + ) import Data.IntMap.Strict ( IntMap, mergeWithKey ) import qualified Data.IntMap.Strict as IM +import qualified Data.Map.Strict as DM +import Data.Maybe ( fromJust, isJust ) import Data.List ( nub, union ) +import Data.List.Extra ( unsnoc ) import Geometry.VertexEnum.Constraint ( Constraint (..), Sense (..) ) import Geometry.VertexEnum.LinearCombination ( LinearCombination (..), VarIndex ) -import Numeric.LinearProgramming ( simplex, - Bound(Free, (:<=:)), - Constraints(Dense), - Optimization(Maximize), - Solution( - Undefined - , Feasible - , Infeasible - , NoFeasible - , Optimal - , Unbounded) ) +import Linear.Simplex.Solver.TwoPhase ( + twoPhaseSimplex + , findFeasibleSolution + ) +import Linear.Simplex.Types ( + Result ( .. ) + , PolyConstraint ( .. ) + , ObjectiveFunction ( .. ) + ) +import Linear.Simplex.Util ( + simplifySystem + ) normalizeLinearCombination :: Num a => [VarIndex] -> LinearCombination a -> IntMap a @@ -29,58 +40,98 @@ varsOfLinearCombo (LinearCombination imap) = IM.keys imap varsOfConstraint :: Constraint a -> [VarIndex] -varsOfConstraint (Constraint lhs _ rhs) = - varsOfLinearCombo lhs `union` varsOfLinearCombo rhs +varsOfConstraint (Constraint left _ right) = + varsOfLinearCombo left `union` varsOfLinearCombo right -normalizeConstraint :: Real a => [VarIndex] -> Constraint a -> [Double] -normalizeConstraint vars (Constraint lhs sense rhs) = +normalizeConstraint :: Real a => [VarIndex] -> Constraint a -> [a] +normalizeConstraint vars (Constraint left sense right) = if sense == Lt then xs ++ [x] else map negate xs ++ [-x] where - lhs' = normalizeLinearCombination vars lhs - rhs' = normalizeLinearCombination vars rhs + lhs' = normalizeLinearCombination vars left + rhs' = normalizeLinearCombination vars right coefs = IM.elems $ mergeWithKey (\_ a b -> Just (a-b)) id id lhs' rhs' - coefs' :: [Double] - coefs' = map realToFrac coefs - (x, xs) = case coefs' of - (xx:xxs) -> (xx, xxs) - [] -> (0, []) - -- let (x:xs) = map realToFrac $ - -- IM.elems $ mergeWithKey (\_ a b -> Just (a-b)) id id lhs' rhs' - -- in - -- if sense == Lt - -- then xs ++ [x] - -- else map negate xs ++ [-x] - -- where lhs' = normalizeLinearCombination vars lhs - -- rhs' = normalizeLinearCombination vars rhs + (x, xs) = case coefs of + (xx:xxs) -> (xx, xxs) + [] -> (0, []) -normalizeConstraints :: Real a => [Constraint a] -> [[Double]] -- for qhalf +normalizeConstraints :: Real a => [Constraint a] -> [[a]] normalizeConstraints constraints = map (normalizeConstraint vars) constraints where vars = nub $ concatMap varsOfConstraint constraints -inequality :: [Double] -> Bound [Double] -inequality row = (coeffs ++ [1.0]) :<=: bound +negateIf :: Bool -> Rational -> Rational +negateIf test x = if test then -x else x + +inequality :: [Bool] -> [Rational] -> PolyConstraint +inequality toNegate row = + LEQ { + lhs = DM.fromList (zip [0 ..] (1 : coeffs')), rhs = -bound + } where - coeffs = init row - bound = -(last row) + (coeffs, bound) = fromJust $ unsnoc row + coeffs' = zipWith negateIf toNegate coeffs -inequalities :: [[Double]] -> Constraints -inequalities normConstraints = Dense (map inequality normConstraints) +inequalities :: [[Rational]] -> [Bool] -> [PolyConstraint] +inequalities normConstraints toNegate = + simplifySystem $ map (inequality toNegate) normConstraints -iPoint :: [[Double]] -> [Double] -iPoint halfspacesMatrix = case solution of - Optimal (_, point) -> init point - Undefined -> error "Failed to find interior point (undefined)." - Feasible (_, _) -> error "Failed to find interior point (feasible)." - Infeasible (_, _) -> error "Failed to find interior point (infeasible)." - NoFeasible -> error "Failed to find interior point (no feasible)." - Unbounded -> error "Failed to find interior point (unbounded)." +-- iPoint does not necessarily return the optimal interior point, because +-- this point possibly corresponds to another [Bool] combination; it just +-- returns a feasible point +iPoint :: [[Rational]] -> [Bool] -> IO [Double] +iPoint halfspacesMatrix toNegate = do + maybeResult <- runStdoutLoggingT $ filterLogger (\_ _ -> False) $ + twoPhaseSimplex objFunc polyConstraints + return $ case maybeResult of + Just (Result var varLitMap) -> + let sol = DM.delete 0 $ DM.delete var varLitMap + nvars = length toNegate + sol' = DM.union sol (DM.fromList (zip [1 .. nvars] (repeat 0))) + in + map fromRational + ( + zipWith negateIf toNegate (DM.elems sol') + ) + Nothing -> error "iPoint: should not happen." where - constraints' = inequalities halfspacesMatrix - n = length (head halfspacesMatrix) - objective = Maximize (replicate (n-1) 0 ++ [1]) - bounds = map Free [1 .. (n-1)] - solution = simplex objective constraints' bounds + polyConstraints = inequalities halfspacesMatrix toNegate + objFunc = Max { + objective = DM.singleton 0 1 + } + +feasiblePoint :: [[Rational]] -> [Bool] -> IO Bool +feasiblePoint halfspacesMatrix toNegate = do + maybeFS <- runStdoutLoggingT $ filterLogger (\_ _ -> False) $ + findFeasibleSolution polyConstraints + return $ isJust maybeFS + where + polyConstraints = simplifySystem $ map ineq halfspacesMatrix + ineq row = + LEQ { + lhs = DM.fromList (zip [1 ..] coeffs'), rhs = -bound + } + where + (coeffs, bound) = fromJust $ unsnoc row + coeffs' = zipWith negateIf toNegate coeffs + +findSigns :: [[Rational]] -> IO [Bool] +findSigns halfspacesMatrix = do + go 0 + where + nvars = length (halfspacesMatrix !! 0) - 1 + combinations = sequence $ replicate nvars [False, True] + ncombinations = length combinations + go i + | i == ncombinations = do + return [] + | otherwise = do + let combo = combinations !! i + test <- feasiblePoint halfspacesMatrix combo + if test + then do + return $ combo + else do + go (i+1)
src/Geometry/VertexEnum/LinearCombination.hs view
@@ -2,11 +2,13 @@ {-# LANGUAGE InstanceSigs #-} module Geometry.VertexEnum.LinearCombination ( LinearCombination (..) + , Var , newVar , VarIndex , linearCombination , constant , cst + , toRationalLinearCombination ) where import Data.AdditiveGroup ( AdditiveGroup(zeroV, negateV, (^+^)) ) @@ -18,6 +20,9 @@ newtype LinearCombination a = LinearCombination (IntMap a) +toRationalLinearCombination :: Real a => LinearCombination a -> LinearCombination Rational +toRationalLinearCombination (LinearCombination imap) = LinearCombination (IM.map toRational imap) + instance (Eq a) => Eq (LinearCombination a) where (==) :: LinearCombination a -> LinearCombination a -> Bool (==) (LinearCombination x) (LinearCombination y) = x == y @@ -55,7 +60,7 @@ newVar :: Num a => VarIndex -> Var a newVar i = if i >= 0 then LinearCombination (IM.singleton i 1) - else error "negative index" + else error "newVar: negative index" -- | linear combination from list of terms linearCombination :: Num a => [(a, Var a)] -> LinearCombination a
src/Geometry/VertexEnum/VertexEnum.hs view
@@ -4,13 +4,14 @@ , interiorPoint ) where import Control.Monad ( unless, when, (<$!>) ) +import Data.Maybe ( isJust, fromJust ) import Foreign.C.Types ( CDouble, CUInt ) import Foreign.Marshal.Alloc ( free, mallocBytes ) import Foreign.Marshal.Array ( peekArray, pokeArray ) import Foreign.Storable ( peek, sizeOf ) import Geometry.VertexEnum.CVertexEnum ( c_intersections ) -import Geometry.VertexEnum.Constraint ( Constraint ) -import Geometry.VertexEnum.Internal ( iPoint, normalizeConstraints ) +import Geometry.VertexEnum.Constraint ( Constraint, toRationalConstraint ) +import Geometry.VertexEnum.Internal ( iPoint, normalizeConstraints, findSigns ) hsintersections :: [[Double]] -- halfspaces -> [Double] -- interior point @@ -20,11 +21,11 @@ let n = length halfspaces dim = length ipoint unless (all ((== dim+1) . length) halfspaces) $ - error "the points must have the same dimension" + error "the length of the point does not match the number of variables." when (dim < 2) $ - error "dimension must be at least 2" + error "dimension must be at least 2." when (n <= dim) $ - error "insufficient number of halfspaces" + error "insufficient number of halfspaces." hsPtr <- mallocBytes (n * (dim+1) * sizeOf (undefined :: CDouble)) pokeArray hsPtr (concatMap (map realToFrac) halfspaces) ipointPtr <- mallocBytes (dim * sizeOf (undefined :: CDouble)) @@ -41,7 +42,7 @@ then do free resultPtr free nintersectionsPtr - error $ "qhull returned an error (code " ++ show exitcode ++ ")" + error $ "qhull returned an error (code " ++ show exitcode ++ ")." else do nintersections <- (<$!>) fromIntegral (peek nintersectionsPtr) result <- (<$!>) (map (map realToFrac)) @@ -52,36 +53,55 @@ return result -- | Vertex enumeration -vertexenum :: Real a => [Constraint a] -- ^ list of inequalities - -> Maybe [Double] -- ^ point in the interior of the polytope - -> IO [[Double]] +vertexenum :: + Real a + => [Constraint a] -- ^ linear inequalities + -> Maybe [Double] -- ^ point satisfying the inequalities, @Nothing@ for automatic point + -> IO [[Double]] -- ^ vertices of the polytope defined by the inequalities vertexenum constraints point = do - let halfspacesMatrix = normalizeConstraints constraints - ipoint = case point of - Just x -> x - Nothing -> iPoint halfspacesMatrix - hsintersections halfspacesMatrix ipoint False + let halfspacesMatrix = + map (map realToFrac) (normalizeConstraints constraints) + if isJust point + then do + let check = checkConstraints constraints (fromJust point) + when (not $ all snd check) $ + error "vertexenum: the provided point does not fulfill the inequalities." + hsintersections halfspacesMatrix (fromJust point) False + else do + ipoint <- interiorPoint constraints + hsintersections halfspacesMatrix ipoint False --- | Check whether a point fulfills some constraints; returns the +-- | Checks whether a point fulfills some inequalities; returns the -- difference between the upper member and the lower member for each --- constraint, which is positive in case if the constraint is fulfilled -checkConstraints :: Real a => [Constraint a] -- ^ list of inequalities - -> [Double] -- ^ point to be tested - -> [(Double, Bool)] -- ^ difference and status for each constraint +-- inequality, which is positive in case if the inequality is fulfilled. +checkConstraints :: + Real a + => [Constraint a] -- ^ linear inequalities + -> [Double] -- ^ point to be tested + -> [(Double, Bool)] -- ^ difference and status for each constraint checkConstraints constraints point = if nvars == length point + 1 then zip differences (map (>= 0) differences) else - error "The length of the point does not match the number of variables." + error "checkConstraints: the length of the point does not match the number of variables." where - halfspacesMatrix = normalizeConstraints constraints - nvars = length (head halfspacesMatrix) + halfspacesMatrix = + map (map realToFrac) (normalizeConstraints constraints) + nvars = length (halfspacesMatrix !! 0) checkRow pt row = - sum (zipWith (*) row (pt ++ [1])) differences = map (checkRow point) halfspacesMatrix --- | Return a point fulfilling a list of constraints -interiorPoint :: Real a => [Constraint a] -> [Double] -interiorPoint constraints = iPoint halfspacesMatrix - where - halfspacesMatrix = normalizeConstraints constraints +-- | Returns a point fulfilling a list of inequalities +interiorPoint :: + Real a + => [Constraint a] -- ^ linear inequalities + -> IO [Double] -- ^ point fulfilling the inequaities +interiorPoint constraints = do + let + constraints' = map toRationalConstraint constraints + halfspacesMatrix = normalizeConstraints constraints' + signs <- findSigns halfspacesMatrix + when (null signs) $ + error "interiorPoint: no feasible point." + iPoint halfspacesMatrix signs
tests/Main.hs view
@@ -1,11 +1,11 @@-module Main where -import Approx ( assertApproxZero ) +module Main ( main ) where +-- import Approx ( assertApproxZero ) import Geometry.VertexEnum ( (.<=), (.>=), Constraint, newVar , vertexenum, interiorPoint, checkConstraints ) import Test.Tasty ( defaultMain, testGroup ) import Test.Tasty.HUnit ( testCase, assertEqual, assertBool ) -cubeConstraints :: [Constraint Double] +cubeConstraints :: [Constraint Rational] cubeConstraints = [ x .<= 1 , x .>= (-1) @@ -27,14 +27,15 @@ let check = checkConstraints cubeConstraints [0, 0, 0] assertBool "" (all snd check) - , testCase "cube vertices" $ do + , testCase "there are 8 cube vertices" $ do vertices <- vertexenum cubeConstraints (Just [0, 0, 0]) assertEqual "" (length vertices) 8 - , testCase "interior point of the cube" $ do - let ipoint = interiorPoint cubeConstraints - norm = (ipoint !! 0) * (ipoint !! 0) + - (ipoint !! 1) * (ipoint !! 1) + (ipoint !! 2) * (ipoint !! 2) - assertApproxZero "" 6 norm + , testCase "interior point of the cube is [0, 0, 0]" $ do + ipoint <- interiorPoint cubeConstraints + -- let norm = (ipoint !! 0) * (ipoint !! 0) + + -- (ipoint !! 1) * (ipoint !! 1) + (ipoint !! 2) * (ipoint !! 2) + -- assertApproxZero "" 6 norm + assertEqual "" ipoint [0, 0, 0] ]
vertexenum.cabal view
@@ -1,7 +1,7 @@ cabal-version: 2.2 name: vertexenum -version: 0.1.1.0 +version: 1.0.0.0 synopsis: Vertex enumeration description: Vertex enumeration of convex polytopes given by linear inequalities. homepage: https://github.com/stla/vertexenum#readme @@ -9,11 +9,11 @@ license-file: LICENSE author: Stéphane Laurent maintainer: laurent_step@outlook.fr -copyright: 2023 Stéphane Laurent +copyright: 2023-2024 Stéphane Laurent category: Math, Geometry build-type: Simple extra-source-files: README.md - CHANGELOG.md +extra-doc-files: CHANGELOG.md library hs-source-dirs: src @@ -24,9 +24,11 @@ , Geometry.VertexEnum.CVertexEnum , Geometry.VertexEnum.VertexEnum build-depends: base >= 4.7 && < 5 - , containers >= 0.6.2.1 && < 0.8 - , hmatrix-glpk >= 0.19.0.0 && < 0.20 + , simplex-method >= 0.2.0.0 && < 0.3 + , containers >= 0.6.5.1 && < 0.7 + , monad-logger >= 0.3.40 && < 0.4 , vector-space >= 0.15 && < 0.17 + , extra >= 1.7 && < 1.8 other-extensions: ForeignFunctionInterface , TypeFamilies , InstanceSigs @@ -81,6 +83,15 @@ , tasty-hunit >= 0.10 && < 0.11 , vertexenum Default-Language: Haskell2010 + ghc-options: -Wall + -Wcompat + -Widentities + -Wincomplete-record-updates + -Wincomplete-uni-patterns + -Wmissing-export-lists + -Wmissing-home-modules + -Wpartial-fields + -Wredundant-constraints source-repository head type: git