vertexenum 0.1.0.0 → 0.1.1.0
raw patch · 8 files changed
+193/−188 lines, 8 files
Files
- CHANGELOG.md +10/−5
- README.md +48/−45
- src/Geometry/VertexEnum/Constraint.hs +11/−8
- src/Geometry/VertexEnum/Internal.hs +8/−9
- src/Geometry/VertexEnum/LinearCombination.hs +23/−28
- src/Geometry/VertexEnum/VertexEnum.hs +5/−5
- tests/Main.hs +1/−1
- vertexenum.cabal +87/−87
CHANGELOG.md view
@@ -1,5 +1,10 @@-# Changelog for `vertexenum`--## 0.1.0.0 - 2023-11-18--First release.+# Changelog for `vertexenum` + +## 0.1.1.0 - 2023-11-20 + +The types `LinearCombination` and `Constraint` are parametric now. + + +## 0.1.0.0 - 2023-11-18 + +First release.
README.md view
@@ -1,46 +1,49 @@-# vertexenum--<!-- badges: start -->-[](https://github.com/stla/vertexenum/actions/workflows/Stack.yml)-<!-- badges: end -->--*Get the vertices of an intersection of halfspaces.*--____--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.$$--Each inequality defines a halfspace. The intersection of the six halfspaces is-a convex polytope. The `vertexenum` function can calculate the vertices of this -polytope:--```haskell-import Data.Ratio ( (%) )-import Data.VectorSpace ( AdditiveGroup((^+^), (^-^))- , VectorSpace((*^)) )-import Geometry.VertexEnum--constraints :: [Constraint]-constraints =- [ x .>= (-5) -- shortcut for `x .>=. cst (-5)`- , x .<= 4- , y .>= (-5)- , y .<=. cst 3 ^-^ x -- we need `cst` here- , z .>= (-10)- , z .<=. cst 6 ^-^ 2*^x ^-^ y ]- where- x = newVar 1- y = newVar 2- z = newVar 3--vertexenum constraints Nothing-```--The type of the second argument of `vertexenum` is `Maybe [Double]`. If this -argument is `Just point`, then `point` must be the coordinates of a point -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 +# vertexenum + +<!-- badges: start --> +[](https://github.com/stla/vertexenum/actions/workflows/Stack.yml) +<!-- badges: end --> + +*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.$$ + +Each inequality defines a halfspace. The intersection of the six halfspaces is +a convex polytope. The `vertexenum` function can calculate the vertices of this +polytope: + +```haskell +import Data.VectorSpace ( AdditiveGroup((^+^), (^-^)) + , VectorSpace((*^)) ) +import Geometry.VertexEnum + +constraints :: [Constraint Double] +constraints = + [ x .>= (-5) -- shortcut for `x .>=. cst (-5)` + , x .<= 4 + , y .>= (-5) + , y .<=. cst 3 ^-^ x -- we need `cst` here + , z .>= (-10) + , z .<=. cst 6 ^-^ 2*^x ^-^ y ] + where + x = newVar 1 + y = newVar 2 + z = newVar 3 + +vertexenum constraints Nothing +``` + +The type of the second argument of `vertexenum` is `Maybe [Double]`. If this +argument is `Just point`, then `point` must be the coordinates of a point +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.
src/Geometry/VertexEnum/Constraint.hs view
@@ -17,19 +17,22 @@ show Gt = ">=" show Lt = "<=" -data Constraint = Constraint LinearCombination Sense LinearCombination - deriving (Eq, Show) +data Constraint a = Constraint (LinearCombination a) Sense (LinearCombination a) -(.>=.) :: LinearCombination -> LinearCombination -> Constraint -(.>=.) lhs rhs = Constraint lhs Gt rhs +instance Show a => Show (Constraint a) where + show :: Constraint a -> String + show (Constraint lhs sense rhs) = show lhs ++ " " ++ show sense ++ " " ++ show rhs -(.<=.) :: LinearCombination -> LinearCombination -> Constraint -(.<=.) lhs rhs = Constraint lhs Lt rhs +(.>=.) :: LinearCombination a -> LinearCombination a -> Constraint a +(.>=.) lhs = Constraint lhs Gt -(.>=) :: LinearCombination -> Rational -> Constraint +(.<=.) :: LinearCombination a -> LinearCombination a -> Constraint a +(.<=.) lhs = Constraint lhs Lt + +(.>=) :: LinearCombination a -> a -> Constraint a (.>=) lhs x = (.>=.) lhs (constant x) -(.<=) :: LinearCombination -> Rational -> Constraint +(.<=) :: LinearCombination a -> a -> Constraint a (.<=) lhs x = (.<=.) lhs (constant x) infix 4 .<=., .>=.
src/Geometry/VertexEnum/Internal.hs view
@@ -6,7 +6,6 @@ import Data.IntMap.Strict ( IntMap, mergeWithKey ) import qualified Data.IntMap.Strict as IM import Data.List ( nub, union ) -import Data.Ratio ( (%), numerator, denominator ) import Geometry.VertexEnum.Constraint ( Constraint (..), Sense (..) ) import Geometry.VertexEnum.LinearCombination ( LinearCombination (..), VarIndex ) import Numeric.LinearProgramming ( simplex, @@ -22,18 +21,18 @@ , Unbounded) ) normalizeLinearCombination :: - [VarIndex] -> LinearCombination -> IntMap Rational + Num a => [VarIndex] -> LinearCombination a -> IntMap a normalizeLinearCombination vars (LinearCombination lc) = IM.union lc (IM.fromList [(i,0) | i <- vars `union` [0]]) -varsOfLinearCombo :: LinearCombination -> [VarIndex] +varsOfLinearCombo :: LinearCombination a -> [VarIndex] varsOfLinearCombo (LinearCombination imap) = IM.keys imap -varsOfConstraint :: Constraint -> [VarIndex] +varsOfConstraint :: Constraint a -> [VarIndex] varsOfConstraint (Constraint lhs _ rhs) = varsOfLinearCombo lhs `union` varsOfLinearCombo rhs -normalizeConstraint :: [VarIndex] -> Constraint -> [Double] +normalizeConstraint :: Real a => [VarIndex] -> Constraint a -> [Double] normalizeConstraint vars (Constraint lhs sense rhs) = if sense == Lt then xs ++ [x] @@ -42,9 +41,9 @@ lhs' = normalizeLinearCombination vars lhs rhs' = normalizeLinearCombination vars rhs coefs = IM.elems $ mergeWithKey (\_ a b -> Just (a-b)) id id lhs' rhs' - denominators = map denominator coefs - ppcm = foldr lcm 1 denominators % 1 - (x, xs) = case map (realToFrac . numerator . (*ppcm)) coefs of + coefs' :: [Double] + coefs' = map realToFrac coefs + (x, xs) = case coefs' of (xx:xxs) -> (xx, xxs) [] -> (0, []) -- let (x:xs) = map realToFrac $ @@ -56,7 +55,7 @@ -- where lhs' = normalizeLinearCombination vars lhs -- rhs' = normalizeLinearCombination vars rhs -normalizeConstraints :: [Constraint] -> [[Double]] -- for qhalf +normalizeConstraints :: Real a => [Constraint a] -> [[Double]] -- for qhalf normalizeConstraints constraints = map (normalizeConstraint vars) constraints where
src/Geometry/VertexEnum/LinearCombination.hs view
@@ -9,68 +9,63 @@ , cst ) where -import Data.AdditiveGroup ( AdditiveGroup(zeroV, negateV, (^+^)) ) +import Data.AdditiveGroup ( AdditiveGroup(zeroV, negateV, (^+^)) ) import Data.IntMap.Strict ( IntMap, mergeWithKey ) import qualified Data.IntMap.Strict as IM import Data.List ( intercalate ) -import Data.Ratio ( numerator, denominator ) import Data.Tuple ( swap ) import Data.VectorSpace ( linearCombo, VectorSpace(..) ) -newtype LinearCombination = LinearCombination (IntMap Rational) - deriving Eq +newtype LinearCombination a = LinearCombination (IntMap a) -instance Show LinearCombination where - show :: LinearCombination -> String +instance (Eq a) => Eq (LinearCombination a) where + (==) :: LinearCombination a -> LinearCombination a -> Bool + (==) (LinearCombination x) (LinearCombination y) = x == y + +instance (Show a) => Show (LinearCombination a) where + show :: LinearCombination a -> String show (LinearCombination x) = intercalate " + " $ map (\(i, r) -> if i == 0 - then showRational r - else if r == 1 - then "x" ++ show i - else showRational r ++ "*x" ++ show i) + then show r + else show r ++ "*x" ++ show i + ) (IM.toAscList x) - where - showRational :: Rational -> String - showRational r = if q == 1 then show p else show p ++ "/" ++ show q - where - p = numerator r - q = denominator r -instance AdditiveGroup LinearCombination where - zeroV :: LinearCombination +instance Num a => AdditiveGroup (LinearCombination a) where + zeroV :: LinearCombination a zeroV = LinearCombination (IM.singleton 0 0) - (^+^) :: LinearCombination -> LinearCombination -> LinearCombination + (^+^) :: LinearCombination a -> LinearCombination a -> LinearCombination a (^+^) (LinearCombination imap1) (LinearCombination imap2) = LinearCombination (mergeWithKey (\_ x y -> Just (x+y)) id id imap1 imap2) - negateV :: LinearCombination -> LinearCombination + negateV :: LinearCombination a -> LinearCombination a negateV (LinearCombination imap) = LinearCombination (IM.map negate imap) -instance VectorSpace LinearCombination where - type Scalar LinearCombination = Rational - (*^) :: Scalar LinearCombination -> LinearCombination -> LinearCombination +instance Num a => VectorSpace (LinearCombination a) where + type Scalar (LinearCombination a) = a + (*^) :: Scalar (LinearCombination a) -> LinearCombination a -> LinearCombination a (*^) lambda (LinearCombination imap) = LinearCombination (IM.map (*lambda) imap) -type Var = LinearCombination +type Var a = LinearCombination a type VarIndex = Int -- | new variable -newVar :: VarIndex -> Var +newVar :: Num a => VarIndex -> Var a newVar i = if i >= 0 then LinearCombination (IM.singleton i 1) else error "negative index" -- | linear combination from list of terms -linearCombination :: [(Rational, Var)] -> LinearCombination +linearCombination :: Num a => [(a, Var a)] -> LinearCombination a linearCombination terms = linearCombo (map swap terms) -- LinearCombination (IM.fromListWith (+) (map swap terms)) -- | constant linear combination -constant :: Rational -> LinearCombination +constant :: a -> LinearCombination a constant x = LinearCombination (IM.singleton 0 x) -- | alias for `constant` -cst :: Rational -> LinearCombination +cst :: a -> LinearCombination a cst = constant
src/Geometry/VertexEnum/VertexEnum.hs view
@@ -13,8 +13,8 @@ import Geometry.VertexEnum.Internal ( iPoint, normalizeConstraints ) hsintersections :: [[Double]] -- halfspaces - -> [Double] -- interior point - -> Bool -- print to stdout + -> [Double] -- interior point + -> Bool -- print to stdout -> IO [[Double]] hsintersections halfspaces ipoint stdout = do let n = length halfspaces @@ -52,7 +52,7 @@ return result -- | Vertex enumeration -vertexenum :: [Constraint] -- ^ list of inequalities +vertexenum :: Real a => [Constraint a] -- ^ list of inequalities -> Maybe [Double] -- ^ point in the interior of the polytope -> IO [[Double]] vertexenum constraints point = do @@ -65,7 +65,7 @@ -- | Check whether a point fulfills some constraints; 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 :: [Constraint] -- ^ list of inequalities +checkConstraints :: Real a => [Constraint a] -- ^ list of inequalities -> [Double] -- ^ point to be tested -> [(Double, Bool)] -- ^ difference and status for each constraint checkConstraints constraints point = @@ -81,7 +81,7 @@ differences = map (checkRow point) halfspacesMatrix -- | Return a point fulfilling a list of constraints -interiorPoint :: [Constraint] -> [Double] +interiorPoint :: Real a => [Constraint a] -> [Double] interiorPoint constraints = iPoint halfspacesMatrix where halfspacesMatrix = normalizeConstraints constraints
tests/Main.hs view
@@ -5,7 +5,7 @@ import Test.Tasty ( defaultMain, testGroup ) import Test.Tasty.HUnit ( testCase, assertEqual, assertBool ) -cubeConstraints :: [Constraint] +cubeConstraints :: [Constraint Double] cubeConstraints = [ x .<= 1 , x .>= (-1)
vertexenum.cabal view
@@ -1,87 +1,87 @@-cabal-version: 2.2--name: vertexenum-version: 0.1.0.0-synopsis: Vertex enumeration-description: Vertex enumeration of convex polytopes.-homepage: https://github.com/stla/vertexenum#readme-license: GPL-3.0-only-license-file: LICENSE-author: Stéphane Laurent-maintainer: laurent_step@outlook.fr-copyright: 2023 Stéphane Laurent-category: Math, Geometry-build-type: Simple-extra-source-files: README.md- CHANGELOG.md--library- hs-source-dirs: src- exposed-modules: Geometry.VertexEnum- other-modules: Geometry.VertexEnum.Constraint- , Geometry.VertexEnum.Internal- , Geometry.VertexEnum.LinearCombination- , 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- , vector-space >= 0.15 && < 0.17- other-extensions: ForeignFunctionInterface- , TypeFamilies- , InstanceSigs- default-language: Haskell2010- include-dirs: C- C-sources: C/libqhull_r.c- , C/geom_r.c- , C/geom2_r.c- , C/global_r.c- , C/io_r.c- , C/mem_r.c- , C/merge_r.c- , C/poly_r.c- , C/poly2_r.c- , C/qset_r.c- , C/random_r.c- , C/usermem_r.c- , C/userprintf_r.c- , C/user_r.c- , C/stat_r.c- , C/halfspaces.c- , C/utils.c- install-includes: C/libqhull_r.h- , C/geom_r.h- , C/io_r.h- , C/mem_r.h- , C/merge_r.h- , C/poly_r.h- , C/qhull_ra.h- , C/qset_r.h- , C/random_r.h- , C/user_r.h- , C/stat_r.h- , C/utils.h- ghc-options: -Wall- -Wcompat- -Widentities- -Wincomplete-record-updates- -Wincomplete-uni-patterns- -Wmissing-export-lists- -Wmissing-home-modules- -Wpartial-fields- -Wredundant-constraints--test-suite unit-tests- type: exitcode-stdio-1.0- main-is: Main.hs- hs-source-dirs: tests/- other-modules: Approx- Build-Depends: base >= 4.7 && < 5- , tasty >= 1.4 && < 1.5- , tasty-hunit >= 0.10 && < 0.11- , vertexenum- Default-Language: Haskell2010--source-repository head- type: git- location: https://github.com/stla/vertexenum+cabal-version: 2.2 + +name: vertexenum +version: 0.1.1.0 +synopsis: Vertex enumeration +description: Vertex enumeration of convex polytopes given by linear inequalities. +homepage: https://github.com/stla/vertexenum#readme +license: GPL-3.0-only +license-file: LICENSE +author: Stéphane Laurent +maintainer: laurent_step@outlook.fr +copyright: 2023 Stéphane Laurent +category: Math, Geometry +build-type: Simple +extra-source-files: README.md + CHANGELOG.md + +library + hs-source-dirs: src + exposed-modules: Geometry.VertexEnum + other-modules: Geometry.VertexEnum.Constraint + , Geometry.VertexEnum.Internal + , Geometry.VertexEnum.LinearCombination + , 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 + , vector-space >= 0.15 && < 0.17 + other-extensions: ForeignFunctionInterface + , TypeFamilies + , InstanceSigs + default-language: Haskell2010 + include-dirs: C + C-sources: C/libqhull_r.c + , C/geom_r.c + , C/geom2_r.c + , C/global_r.c + , C/io_r.c + , C/mem_r.c + , C/merge_r.c + , C/poly_r.c + , C/poly2_r.c + , C/qset_r.c + , C/random_r.c + , C/usermem_r.c + , C/userprintf_r.c + , C/user_r.c + , C/stat_r.c + , C/halfspaces.c + , C/utils.c + install-includes: C/libqhull_r.h + , C/geom_r.h + , C/io_r.h + , C/mem_r.h + , C/merge_r.h + , C/poly_r.h + , C/qhull_ra.h + , C/qset_r.h + , C/random_r.h + , C/user_r.h + , C/stat_r.h + , C/utils.h + ghc-options: -Wall + -Wcompat + -Widentities + -Wincomplete-record-updates + -Wincomplete-uni-patterns + -Wmissing-export-lists + -Wmissing-home-modules + -Wpartial-fields + -Wredundant-constraints + +test-suite unit-tests + type: exitcode-stdio-1.0 + main-is: Main.hs + hs-source-dirs: tests/ + other-modules: Approx + Build-Depends: base >= 4.7 && < 5 + , tasty >= 1.4 && < 1.5 + , tasty-hunit >= 0.10 && < 0.11 + , vertexenum + Default-Language: Haskell2010 + +source-repository head + type: git + location: https://github.com/stla/vertexenum