pcubature 0.1.0.0 → 0.2.0.0
raw patch · 5 files changed
+398/−339 lines, 5 filesdep +pcubaturedep +tastydep +tasty-hunitdep ~containersdep ~hspraydep ~vertexenum
Dependencies added: pcubature, tasty, tasty-hunit
Dependency ranges changed: containers, hspray, vertexenum
Files
- CHANGELOG.md +11/−5
- README.md +284/−287
- pcubature.cabal +62/−42
- src/Numeric/Integration/PolyhedralCubature.hs +7/−5
- tests/Main.hs +34/−0
CHANGELOG.md view
@@ -1,5 +1,11 @@-# Changelog for `pcubature`--## 0.1.0.0 - 2023-11-20--First release.+# Changelog for `pcubature` + + +## 0.1.0.0 - 2023-11-20 + +First release. + + +## 0.2.0.0 - 2024-05-XX + +The package does no longer depend on the 'hmatrix-glpk' package.
README.md view
@@ -1,288 +1,285 @@-# pcubature--<!-- badges: start -->-[](https://github.com/stla/pcubature/actions/workflows/Stack.yml)-<!-- badges: end -->--*Multiple integration over convex polytopes.*--***Warning:*** the package does not work in GHCi.--***Info:*** the package indirectly depends on the **hmatrix-glpk** package; -follow [this link](https://github.com/haskell-numerics/hmatrix/blob/master/INSTALL.md)-for installation instructions.--___--This package allows to evaluate a multiple integral over a convex polytope. -Let's consider for example the following integral:--$$\int_0^1\int_0^1\int_0^1 \exp(x+y+z)\,\text{d}z\,\text{d}y\,\text{d}x = {(e-1)}^3 \approx 5.07321411177285.$$--The domain of integration is the cube ${[0,1]}^3$. In order to use the package, -one has to provide the vertices of this cube:--```haskell-integrateOnPolytope'- :: (Vector Double -> Double) -- ^ integrand- -> [[Double]] -- ^ vertices of the polytope- -> Int -- ^ maximum number of evaluations- -> Double -- ^ desired absolute error- -> Double -- ^ desired relative error- -> Int -- ^ integration rule: 1, 2, 3 or 4- -> IO Result -- ^ value, error estimate, evaluations, success-```--Let's go:--```haskell-module Main - where-import Numeric.Integration.PolyhedralCubature-import Data.Vector.Unboxed as V--f :: Vector Double -> Double-f v = exp (V.sum v)--cube :: [[Double]]-cube = [- [0, 0, 0]- , [0, 0, 1]- , [0, 1, 0]- , [0, 1, 1]- , [1, 0, 0]- , [1, 0, 1]- , [1, 1, 0]- , [1, 1, 1]- ]--integral :: IO Result-integral = integrateOnPolytope' f cube 10000 0 1e-6 3--main :: IO ()-main = do - i <- integral- print i--- Result {--- value = 5.073214090351428--- , errorEstimate = 2.8421152805879766e-6--- , evaluations = 710--- , success = True--- }-```--This cube is axis-aligned. So it may be better to use the **adaptive-cubature** -package here. The **pcubature** package allows to evaluate multiple integrals -whose bounds are (roughly speaking) linear combinations of the variables, -such as:--$$\int_{-5}^4\int_{-5}^{3-x}\int_{-10}^{6-2x-y} f(x, y, z)\,\text{d}z\,\text{d}y\,\text{d}x.$$--Here, the domain of integration is given by the set of linear inequalities:--$$\left\{\begin{matrix} -5 & \leq & x & \leq & 4 \\\ -5 & \leq & y & \leq & 3-x \\\ -10 & \leq & z & \leq & 6-2x-y \end{matrix}\right.$$--Each of these linear inequalities defines a halfspace of $\mathbb{R}^3$, and -the intersection of these six halfspaces is a convex polytope (a polyhedron).--But it is not easy to get the vertices of this polytope. This is why the -**pcubature** package depends on the **vertexenum** package, whose purpose is -to enumerate the vertices of a polytope given as above, with linear -inequalities. Let's take as example the function $f(x,y,z) = x(x+1) - yz^2$:--```haskell-module Main- where-import Numeric.Integration.PolyhedralCubature-import Geometry.VertexEnum-import Data.VectorSpace ( - AdditiveGroup((^+^), (^-^))- , VectorSpace((*^)) - )-import Data.Vector.Unboxed as V--f :: Vector Double -> Double-f v = x * (x+1) - y * z * z- where- x = v ! 0- y = v ! 1- z = v ! 2--polytope :: [Constraint Double]-polytope = [- 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--integral :: IO Result-integral = integrateOnPolytope' f polytope 10000 0 1e-6 3--main :: IO ()-main = do - i <- integral- print i--- Result {--- value = 74321.77499999988--- , errorEstimate = 1.0533262499999988e-7--- , evaluations = 330--- , success = True--- }-```--The exact value of this integral is $74321.775$, as we shall see later.--The function $f$ of this example is polynomial. So we can use the function -`integratePolynomialOnPolytope` to integrate it. This requires to define -the polynomial with the help of the **hspray** package; we also import some -modules of the **numeric-prelude** package, which allows to define a **hspray** -polynomial more conveniently:--```haskell-module Main- where-import Numeric.Integration.PolyhedralCubature-import Geometry.VertexEnum-import Data.VectorSpace ( - AdditiveGroup((^+^), (^-^))- , VectorSpace((*^)) - )-import Math.Algebra.Hspray ( Spray, lone, (^**^) )-import Prelude hiding ( (*), (+), (-) )-import qualified Prelude as P-import Algebra.Additive -import Algebra.Module -import Algebra.Ring--p :: Spray Double-p = x * (x + one) - (y * z^**^2) - where- x = lone 1 :: Spray Double- y = lone 2 :: Spray Double- z = lone 3 :: Spray Double--polytope :: [Constraint Double]-polytope = [- 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--integral :: IO Double-integral = integratePolynomialOnPolytope' p polytope--main :: IO ()-main = do - i <- integral- print i--- 74321.77499999967-```--The function `integratePolynomialOnSimplex` implements an exact procedure. -However we didn't get the exact result. That's because of (small) -numerical errors. The first numerical errors occur in the vertex enumeration -performed by the **vertexenum** package:--```haskell-module Main- where-import Geometry.VertexEnum-import Data.VectorSpace ( - AdditiveGroup((^+^), (^-^))- , VectorSpace((*^)) - )--polytope :: [Constraint Double]-polytope = [- x .>= (-5) - , x .<= 4- , y .>= (-5)- , y .<=. cst 3 ^-^ x - , z .>= (-10)- , z .<=. cst 6 ^-^ 2*^x ^-^ y - ]- where- x = newVar 1- y = newVar 2- z = newVar 3--vertices :: IO [[Double]]-vertices = vertexenum polytope Nothing--main :: IO ()-main = do - vs <- vertices- print vs--- [--- [-5.000000000000003, 8.000000000000004, 8.000000000000004]--- , [-4.999999999999998, -4.999999999999996, 20.999999999999993]--- , [3.999999999999999, -0.9999999999999997, -1.0]--- , [3.999999999999999, -5.0, 3.0000000000000004]--- , [-5.0, -5.0, -10.0]--- , [-5.0, 8.000000000000002, -10.0]--- , [4.0, -0.9999999999999999, -10.0]--- , [4.0, -5.0, -10.0]--- ]-```--Since all coefficients of the linear inequalities are rational (they even are -integral), the vertices should be rational as well. -Unfortunately, **vertexenum** only allows to get vertices with double -coordinates. So if we want to use `Rational`, we have to manually enter -the vertices:--```haskell-module Main- where-import Numeric.Integration.PolyhedralCubature-import Math.Algebra.Hspray ( Spray, lone, (^**^) )-import Prelude hiding ( (*), (+), (-) )-import qualified Prelude as P-import Algebra.Additive -import Algebra.Module -import Algebra.Ring--p :: Spray Rational-p = x * (x + one) - (y * z^**^2) - where- x = lone 1 :: Spray Rational- y = lone 2 :: Spray Rational- z = lone 3 :: Spray Rational--polytope :: [[Rational]]-polytope = [- [-5, 8, 8]- , [-5, -5, 21]- , [4, -1, -1]- , [4, -5, 3]- , [-5, -5, -10]- , [-5, 8, -10]- , [4, -1, -10]- , [4, -5, -10]- ]--integral :: IO Rational-integral = integratePolynomialOnPolytope p polytope--main :: IO ()-main = do - i <- integral- print i--- 2972871 % 40-```-+# pcubature + +<!-- badges: start --> +[](https://github.com/stla/pcubature/actions/workflows/Stack-lts.yml) +[](https://github.com/stla/pcubature/actions/workflows/Stack-nightly.yml) +<!-- badges: end --> + +***Multiple integration over convex polytopes.*** + +***Warning:*** the package does not work in GHCi. + +___ + +This package allows to evaluate a multiple integral over a convex polytope. +Let's consider for example the following integral: + +$$\int_0^1\int_0^1\int_0^1 \exp(x+y+z)\,\text{d}z\,\text{d}y\,\text{d}x = {(e-1)}^3 \approx 5.07321411177285.$$ + +The domain of integration is the cube ${[0,1]}^3$. In order to use the package, +one has to provide the vertices of this cube: + +```haskell +integrateOnPolytope' + :: (Vector Double -> Double) -- ^ integrand + -> [[Double]] -- ^ vertices of the polytope + -> Int -- ^ maximum number of evaluations + -> Double -- ^ desired absolute error + -> Double -- ^ desired relative error + -> Int -- ^ integration rule: 1, 2, 3 or 4 + -> IO Result -- ^ value, error estimate, evaluations, success +``` + +Let's go: + +```haskell +module Main + where +import Numeric.Integration.PolyhedralCubature +import Data.Vector.Unboxed as V + +f :: Vector Double -> Double +f v = exp (V.sum v) + +cube :: [[Double]] +cube = [ + [0, 0, 0] + , [0, 0, 1] + , [0, 1, 0] + , [0, 1, 1] + , [1, 0, 0] + , [1, 0, 1] + , [1, 1, 0] + , [1, 1, 1] + ] + +integral :: IO Result +integral = integrateOnPolytope' f cube 10000 0 1e-6 3 + +main :: IO () +main = do + i <- integral + print i +-- Result { +-- value = 5.073214090351428 +-- , errorEstimate = 2.8421152805879766e-6 +-- , evaluations = 710 +-- , success = True +-- } +``` + +This cube is axis-aligned. So it may be better to use the **adaptive-cubature** +package here. The **pcubature** package allows to evaluate multiple integrals +whose bounds are (roughly speaking) linear combinations of the variables, +such as: + +$$\int_{-5}^4\int_{-5}^{3-x}\int_{-10}^{6-2x-y} f(x, y, z)\,\text{d}z\,\text{d}y\,\text{d}x.$$ + +Here, the domain of integration is given by the set of linear inequalities: + +$$\left\{\begin{matrix} -5 & \leq & x & \leq & 4 \\\ -5 & \leq & y & \leq & 3-x \\\ -10 & \leq & z & \leq & 6-2x-y \end{matrix}\right.$$ + +Each of these linear inequalities defines a halfspace of $\mathbb{R}^3$, and +the intersection of these six halfspaces is a convex polytope (a polyhedron). + +But it is not easy to get the vertices of this polytope. This is why the +**pcubature** package depends on the **vertexenum** package, whose purpose is +to enumerate the vertices of a polytope given as above, with linear +inequalities. Let's take as example the function $f(x,y,z) = x(x+1) - yz^2$: + +```haskell +module Main + where +import Numeric.Integration.PolyhedralCubature +import Geometry.VertexEnum +import Data.VectorSpace ( + AdditiveGroup((^+^), (^-^)) + , VectorSpace((*^)) + ) +import Data.Vector.Unboxed as V + +f :: Vector Double -> Double +f v = x * (x+1) - y * z * z + where + x = v ! 0 + y = v ! 1 + z = v ! 2 + +polytope :: [Constraint Double] +polytope = [ + 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 + +integral :: IO Result +integral = integrateOnPolytope' f polytope 10000 0 1e-6 3 + +main :: IO () +main = do + i <- integral + print i +-- Result { +-- value = 74321.77499999988 +-- , errorEstimate = 1.0533262499999988e-7 +-- , evaluations = 330 +-- , success = True +-- } +``` + +The exact value of this integral is $74321.775$, as we shall see later. + +The function $f$ of this example is polynomial. So we can use the function +`integratePolynomialOnPolytope` to integrate it. This requires to define +the polynomial with the help of the **hspray** package; we also import some +modules of the **numeric-prelude** package, which allows to define a **hspray** +polynomial more conveniently: + +```haskell +module Main + where +import Numeric.Integration.PolyhedralCubature +import Geometry.VertexEnum +import Data.VectorSpace ( + AdditiveGroup((^+^), (^-^)) + , VectorSpace((*^)) + ) +import Math.Algebra.Hspray ( Spray, lone, (^**^) ) +import Prelude hiding ( (*), (+), (-) ) +import qualified Prelude as P +import Algebra.Additive +import Algebra.Module +import Algebra.Ring + +p :: Spray Double +p = x * (x + one) - (y * z^**^2) + where + x = lone 1 :: Spray Double + y = lone 2 :: Spray Double + z = lone 3 :: Spray Double + +polytope :: [Constraint Double] +polytope = [ + 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 + +integral :: IO Double +integral = integratePolynomialOnPolytope' p polytope + +main :: IO () +main = do + i <- integral + print i +-- 74321.77499999967 +``` + +The function `integratePolynomialOnPolytope` implements an exact procedure. +However we didn't get the exact result. That's because of (small) +numerical errors. The first numerical errors occur in the vertex enumeration +performed by the **vertexenum** package: + +```haskell +module Main + where +import Geometry.VertexEnum +import Data.VectorSpace ( + AdditiveGroup((^+^), (^-^)) + , VectorSpace((*^)) + ) + +polytope :: [Constraint Double] +polytope = [ + x .>= (-5) + , x .<= 4 + , y .>= (-5) + , y .<=. cst 3 ^-^ x + , z .>= (-10) + , z .<=. cst 6 ^-^ 2*^x ^-^ y + ] + where + x = newVar 1 + y = newVar 2 + z = newVar 3 + +vertices :: IO [[Double]] +vertices = vertexenum polytope Nothing + +main :: IO () +main = do + vs <- vertices + print vs +-- [ +-- [-5.000000000000003, 8.000000000000004, 8.000000000000004] +-- , [-4.999999999999998, -4.999999999999996, 20.999999999999993] +-- , [3.999999999999999, -0.9999999999999997, -1.0] +-- , [3.999999999999999, -5.0, 3.0000000000000004] +-- , [-5.0, -5.0, -10.0] +-- , [-5.0, 8.000000000000002, -10.0] +-- , [4.0, -0.9999999999999999, -10.0] +-- , [4.0, -5.0, -10.0] +-- ] +``` + +Since all coefficients of the linear inequalities are rational (they even are +integral), the vertices should be rational as well. +Unfortunately, **vertexenum** only allows to get vertices with double +coordinates. So if we want to use `Rational`, we have to manually enter +the vertices: + +```haskell +module Main + where +import Numeric.Integration.PolyhedralCubature +import Math.Algebra.Hspray ( Spray, lone, (^**^) ) +import Prelude hiding ( (*), (+), (-) ) +import qualified Prelude as P +import Algebra.Additive +import Algebra.Module +import Algebra.Ring + +p :: Spray Rational +p = x * (x + one) - (y * z^**^2) + where + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + +polytope :: [[Rational]] +polytope = [ + [-5, 8, 8] + , [-5, -5, 21] + , [4, -1, -1] + , [4, -5, 3] + , [-5, -5, -10] + , [-5, 8, -10] + , [4, -1, -10] + , [4, -5, -10] + ] + +integral :: IO Rational +integral = integratePolynomialOnPolytope p polytope + +main :: IO () +main = do + i <- integral + print i +-- 2972871 % 40 +``` + We get it, the exact value $74321.775$, as promised.
pcubature.cabal view
@@ -1,42 +1,62 @@-cabal-version: 2.2--name: pcubature-version: 0.1.0.0-synopsis: Integration over convex polytopes-description: Multiple integration over convex polytopes.-homepage: https://github.com/stla/pcubature#readme-license: GPL-3.0-only-license-file: LICENSE-author: Stéphane Laurent-maintainer: laurent_step@outlook.fr-copyright: 2023 Stéphane Laurent-category: Numeric, Integration-build-type: Simple-extra-source-files: README.md- CHANGELOG.md--library- hs-source-dirs: src- exposed-modules: Numeric.Integration.PolyhedralCubature- build-depends: base >= 4.7 && < 5- , containers >= 0.6.2.1 && < 0.8- , delaunayNd >= 0.1.0.2 && < 0.2- , hspray >= 0.1.0.0 && < 0.2- , numeric-prelude >= 0.4.4 && < 0.5- , scubature >= 1.1.0.0 && < 1.2- , vector >= 0.12.3 && < 0.14- , vertexenum >= 0.1.1.0 && < 0.2- 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- location: https://github.com/stla/pcubature+cabal-version: 2.2 + +name: pcubature +version: 0.2.0.0 +synopsis: Integration over convex polytopes +description: Multiple integration over convex polytopes. +homepage: https://github.com/stla/pcubature#readme +license: GPL-3.0-only +license-file: LICENSE +author: Stéphane Laurent +maintainer: laurent_step@outlook.fr +copyright: 2023-2024 Stéphane Laurent +category: Numeric, Integration +build-type: Simple +extra-source-files: README.md +extra-doc-files: CHANGELOG.md + +library + hs-source-dirs: src + exposed-modules: Numeric.Integration.PolyhedralCubature + build-depends: base >= 4.7 && < 5 + , containers >= 0.6.5.1 && < 0.7 + , delaunayNd >= 0.1.0.2 && < 0.2 + , hspray >= 0.1.0.0 && < 0.5.3 + , numeric-prelude >= 0.4.4 && < 0.5 + , scubature >= 1.1.0.0 && < 1.2 + , vector >= 0.12.3 && < 0.14 + , vertexenum >= 1.0.0.0 && < 1.1 + default-language: Haskell2010 + 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/ + Build-Depends: base >= 4.7 && < 5 + , tasty >= 1.4 && < 1.5 + , tasty-hunit >= 0.10 && < 0.11 + , hspray >= 0.1.0.0 && < 0.5.3 + , pcubature + 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 + location: https://github.com/stla/pcubature
src/Numeric/Integration/PolyhedralCubature.hs view
@@ -1,20 +1,22 @@ {-| Module : Numeric.Integration.PolyhedralCubature Description : Multiple integration over convex polytopes. -Copyright : (c) Stéphane Laurent, 2023 +Copyright : (c) Stéphane Laurent, 2023-2024 License : GPL-3 Maintainer : laurent_step@outlook.fr Evaluation of integrals over a convex polytope. See README for examples. -} module Numeric.Integration.PolyhedralCubature - ( integrateOnPolytopeN + ( + Result(..) + , Results(..) + , Constraint(..) + , VectorD + , integrateOnPolytopeN , integrateOnPolytope , integrateOnPolytopeN' , integrateOnPolytope' - , Result(..) - , Results(..) - , Constraint(..) , integratePolynomialOnPolytope , integratePolynomialOnPolytope' )
+ tests/Main.hs view
@@ -0,0 +1,34 @@+module Main ( main ) where +import Data.Ratio ( (%) ) +import Numeric.Integration.PolyhedralCubature ( integratePolynomialOnPolytope ) +import Math.Algebra.Hspray ( Spray, lone, (^**^), (^*^), (^+^), (^-^), unitSpray ) +import Test.Tasty ( defaultMain, testGroup ) +import Test.Tasty.HUnit ( testCase, assertEqual ) + + +main :: IO () +main = defaultMain $ + testGroup "Tests" + [ + + testCase "exact integral of a polynomial" $ do + let + x = lone 1 :: Spray Rational + y = lone 2 :: Spray Rational + z = lone 3 :: Spray Rational + p = x ^*^ (x ^+^ unitSpray) ^-^ (y ^*^ z^**^2) + polytope :: [[Rational]] + polytope = [ + [-5, 8, 8] + , [-5, -5, 21] + , [4, -1, -1] + , [4, -5, 3] + , [-5, -5, -10] + , [-5, 8, -10] + , [4, -1, -10] + , [4, -5, -10] + ] + integral <- integratePolynomialOnPolytope p polytope + assertEqual "" integral (2972871 % 40) + + ]