scubature 1.0.0.1 → 1.1.0.0

raw patch · 5 files changed

+129/−5 lines, 5 filesdep +hspraydep +numeric-preludePVP ok

version bump matches the API change (PVP)

Dependencies added: hspray, numeric-prelude

API changes (from Hackage documentation)

Files

• CHANGELOG.md +5/−0
• README.md +75/−2
• scubature.cabal +5/−2
• src/Numeric/Integration/IntegratePolynomialOnSimplex.hs +43/−0
• src/Numeric/Integration/SimplexCubature.hs +1/−1

@@ -5,3 +5,8 @@ 1.0.0.1 ------- * removed the upper bounds of the dependencies++1.1.0.0+-------+* new function `integratePolynomialOnSimplex`, returning the exact value of +  the integral of a polynomial over a simplex

@@ -2,6 +2,18 @@  Pure Haskell implementation of simplicial cubature (integration over a simplex). +This library is a port of a part of the R package **SimplicalCubature**, +written by John P. Nolan, and which contains R translations of +some Matlab and Fortran code written by Alan Genz. +It is also a port of a part of the R package **SphericalCubature**, also +written by John P. Nolan. +In addition it provides a function for the exact computation of the integral +of a polynomial over a simplex.++___++## Integral of a function on a simplex+ ```haskell integrateOnSimplex     :: (VectorD -> VectorD)   -- integrand@@ -14,7 +26,7 @@     -> IO Results             -- values, error estimates, evaluations, success ``` -## Example+### Example  ![equation](http://latex.codecogs.com/gif.latex?%5Cint_0%5E1%5Cint_0%5Ex%5Cint_0%5Ey%5Cexp%28x+y+z%29%5C,%5Cmathrm%7Bd%7Dz%5C,%5Cmathrm%7Bd%7Dy%5C,%5Cmathrm%7Bd%7Dx=%5Cfrac%7B1%7D%7B6%7D%28e-1%29%5E3%5Capprox%20.8455356853) @@ -28,7 +40,7 @@ :} ``` -Define the simplex:+Define the simplex (tetrahedron in dimension 3) by the list of its vertices:  ```haskell simplex = [[0, 0, 0], [1, 1, 1], [0, 1, 1], [0, 0, 1]]@@ -65,11 +77,61 @@ --        , success       = True } ``` ++## Exact integral of a polynomial on a simplex++```haskell+integratePolynomialOnSimplex+  :: (C a, Fractional a, Ord a) -- `C a` means that `a` must be a ring+  => Spray a -- ^ polynomial to be integrated+  -> [[a]]   -- ^ simplex to integrate over+  -> a+```++### Example++We take as an example the rational numbers for `a`. Thus we must take a +polynomial with rational coefficients and a simplex whose vertices +have rational coordinates. Then the integral will be a rational number.++Our polynomial is ++![equation](https://latex.codecogs.com/gif.image?\dpi{110}P(x,&space;y,&space;z)&space;=&space;x^4&space;+&space;y&space;+&space;2xy^2&space;-&space;3z.)++It must be defined in Haskell with the +[**hspray**](https://github.com/stla/hspray) library.++```haskell+import Numeric.Integration.IntegratePolynomialOnSimplex+import Data.Ratio+import Math.Algebra.Hspray ++:{+simplex :: [[Rational]]+simplex = [[1, 1, 1], [2, 2, 3], [3, 4, 5], [3, 2, 1]]+:}++x = lone 1 :: Spray Rational+y = lone 2 :: Spray Rational+z = lone 3 :: Spray Rational++:{+poly :: Spray Rational+poly = x^**^4 ^+^ y ^+^ 2.^(x ^*^ y^**^2) ^-^ 3.^z+:}++integratePolynomialOnSimplex poly simplex+-- 1387 % 42+```++ ## Integration on a spherical triangle  The library also allows to evaluate an integral on a spherical simplex on the unit sphere (in dimension 3, a spherical triangle). +### Example+ For example take the first orthant in dimension 3:  ```haskell@@ -99,3 +161,14 @@ --        , evaluations   = 17065 --        , success       = True } ```+++## References++- A. Genz and R. Cools. +  *An adaptive numerical cubature algorithm for simplices.* +  ACM Trans. Math. Software 29, 297-308 (2003).++- Jean B. Lasserre.+  *Simple formula for the integration of polynomials on a simplex.* +  BIT Numerical Mathematics 61, 523-533 (2021).

@@ -1,6 +1,6 @@ cabal-version:      >=1.10 name:               scubature-version:            1.0.0.1+version:            1.1.0.0 license:            GPL-3 license-file:       LICENSE copyright:          2022 Stéphane Laurent@@ -25,6 +25,7 @@     exposed-modules:         Numeric.Integration.SimplexCubature         Numeric.Integration.SphericalSimplexCubature+        Numeric.Integration.IntegratePolynomialOnSimplex      hs-source-dirs:   src     other-modules:@@ -41,4 +42,6 @@         vector >=0.12.3.1,         matrix >=0.3.6.1,         containers >=0.6.4.1,-        ilist >=0.4.0.1+        ilist >=0.4.0.1,+        hspray >=0.1.1.0,+        numeric-prelude >=0.4.4

@@ -0,0 +1,43 @@+module Numeric.Integration.IntegratePolynomialOnSimplex where+import           Algebra.Ring                   ( C )+import           Data.List                      ( transpose )+import           Data.Matrix                    ( detLU+                                                , fromLists+                                                )+import           Math.Algebra.Hspray            ( (*^)+                                                , Spray+                                                , (^+^)+                                                , bombieriSpray+                                                , composeSpray+                                                , constantSpray+                                                , lone+                                                , toList+                                                )++-- | Exact integral of a polynomial over a simplex+integratePolynomialOnSimplex+  :: (C a, Fractional a, Ord a) +  => Spray a -- ^ polynomial to be integrated+  -> [[a]]   -- ^ simplex to integrate over+  -> a+integratePolynomialOnSimplex p simplex =+  s * abs (detLU $ fromLists b) / (fromIntegral $ product [2 .. n])+ where+  v            = last simplex+  n            = length v+  b            = map (\column -> zipWith (-) column v) (take n simplex)+  vb           = zip v (transpose b)+  variables    = map lone [1 .. n]+  newvariables = map+    (\(vi, bi) ->+      (constantSpray vi) ^+^ foldl1 (^+^) (zipWith (*^) bi variables)+    )+    vb+  q      = composeSpray p newvariables+  qterms = toList $ bombieriSpray q+  s      = sum $ map f qterms+   where+    f (exponents, coef) = if d == 0+      then coef+      else coef / (fromIntegral $ product [n + 1 .. n + d])+      where d = sum exponents

@@ -6,7 +6,7 @@ import           Data.Array.Unsafe                   (unsafeThaw) import qualified Data.Vector.Unboxed                 as UV import           Numeric.Integration.Simplex.Simplex ( isValidSimplices, Simplices )-import Numeric.Integration.SimplexCubature.Internal  ( VectorD, IO3dArray, adsimp )+import           Numeric.Integration.SimplexCubature.Internal  ( VectorD, IO3dArray, adsimp )  data Results = Results   { values         :: [Double]