diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,46 +1,49 @@ -# vertexenum - - -[![Stack](https://github.com/stla/vertexenum/actions/workflows/Stack.yml/badge.svg)](https://github.com/stla/vertexenum/actions/workflows/Stack.yml) - - -*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 + + +[![Stack](https://github.com/stla/vertexenum/actions/workflows/Stack.yml/badge.svg)](https://github.com/stla/vertexenum/actions/workflows/Stack.yml) + + +*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. diff --git a/src/Geometry/VertexEnum/Constraint.hs b/src/Geometry/VertexEnum/Constraint.hs --- a/src/Geometry/VertexEnum/Constraint.hs +++ b/src/Geometry/VertexEnum/Constraint.hs @@ -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 .<=., .>=. diff --git a/src/Geometry/VertexEnum/Internal.hs b/src/Geometry/VertexEnum/Internal.hs --- a/src/Geometry/VertexEnum/Internal.hs +++ b/src/Geometry/VertexEnum/Internal.hs @@ -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 diff --git a/src/Geometry/VertexEnum/LinearCombination.hs b/src/Geometry/VertexEnum/LinearCombination.hs --- a/src/Geometry/VertexEnum/LinearCombination.hs +++ b/src/Geometry/VertexEnum/LinearCombination.hs @@ -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 diff --git a/src/Geometry/VertexEnum/VertexEnum.hs b/src/Geometry/VertexEnum/VertexEnum.hs --- a/src/Geometry/VertexEnum/VertexEnum.hs +++ b/src/Geometry/VertexEnum/VertexEnum.hs @@ -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 diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -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) diff --git a/vertexenum.cabal b/vertexenum.cabal --- a/vertexenum.cabal +++ b/vertexenum.cabal @@ -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