probability-polynomial (empty) → 1.0.0.0
raw patch · 18 files changed
+3246/−0 lines, 18 filesdep +QuickCheckdep +basedep +containers
Dependencies added: QuickCheck, base, containers, criterion, deepseq, exact-combinatorics, hspec, probability-polynomial
Files
- CHANGELOG.md +8/−0
- LICENSE +28/−0
- README.md +1/−0
- benchmark/Main.hs +54/−0
- probability-polynomial.cabal +98/−0
- src/Data/Function/Class.hs +58/−0
- src/Numeric/Function/Piecewise.hs +268/−0
- src/Numeric/Measure/Discrete.hs +149/−0
- src/Numeric/Measure/Finite/Mixed.hs +381/−0
- src/Numeric/Measure/Probability.hs +192/−0
- src/Numeric/Polynomial/Simple.hs +618/−0
- src/Numeric/Probability/Moments.hs +91/−0
- test/Numeric/Function/PiecewiseSpec.hs +290/−0
- test/Numeric/Measure/DiscreteSpec.hs +140/−0
- test/Numeric/Measure/Finite/MixedSpec.hs +216/−0
- test/Numeric/Measure/ProbabilitySpec.hs +252/−0
- test/Numeric/Polynomial/SimpleSpec.hs +401/−0
- test/Spec.hs +1/−0
+ CHANGELOG.md view
@@ -0,0 +1,8 @@+# Revision history for `probability-polynomial`++## 1.0.0.0 — 2024-12-23++* Initial release+ * Polynomials+ * Finite, signed measures on the number line+ * Probability measures
+ LICENSE view
@@ -0,0 +1,28 @@+BSD 3-Clause License++Copyright (c) 2020-2024, Predictable Network Solutions Ltd.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++1. Redistributions of source code must retain the above copyright notice, this+ list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright notice,+ this list of conditions and the following disclaimer in the documentation+ and/or other materials provided with the distribution.++3. Neither the name of the copyright holder nor the names of its+ contributors may be used to endorse or promote products derived from+ this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,1 @@+Probability distributions, represented by piecewise polynomials.
+ benchmark/Main.hs view
@@ -0,0 +1,54 @@+{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE DeriveGeneric #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+-}+module Main (main) where++import Criterion.Main+ ( bench+ , bgroup+ , defaultMain+ , nf+ )+import Numeric.Polynomial.Simple+ ( Poly+ )++import qualified Numeric.Polynomial.Simple as Poly++longPoly :: (Integral b, Floating a, Eq a) => b -> Poly a+longPoly m = Poly.fromCoefficients $ replicate (2 ^ m) pi++mulLongPolys :: Int -> Poly Double+mulLongPolys n = longPoly n * longPoly n++addLongPolys :: Int -> Poly Double+addLongPolys n = longPoly n + longPoly n++convPolys :: Int -> [(Double, Poly Double)]+convPolys n = Poly.convolve (0, 1, Poly.constant 1) (0, 1, longPoly n)++main :: IO ()+main =+ defaultMain+ [ bgroup+ "con"+ [ bench "a1" $ nf addLongPolys 1+ , bench "a5" $ nf addLongPolys 5+ , bench "a10" $ nf addLongPolys 10+ , bench "a15" $ nf addLongPolys 15+ , bench "a20" $ nf addLongPolys 20+ , bench "m1" $ nf mulLongPolys 1+ , bench "m3" $ nf mulLongPolys 3+ , bench "m5" $ nf mulLongPolys 5+ , bench "m7" $ nf mulLongPolys 7+ , bench "m9" $ nf mulLongPolys 9+ , bench "c1" $ nf convPolys 1+ , bench "c2" $ nf convPolys 2+ , bench "c3" $ nf convPolys 3+ , bench "c4" $ nf convPolys 4+ ]+ ]
+ probability-polynomial.cabal view
@@ -0,0 +1,98 @@+cabal-version: 3.0+name: probability-polynomial++-- Package Versioning Policy: https://pvp.haskell.org+-- PVP summary: +-+------- breaking API changes+-- | | +----- non-breaking API additions+-- | | | +--- code changes with no API change+version: 1.0.0.0+synopsis: Probability distributions via piecewise polynomials+description:+ Package for manipulating finite probability distributions.++ Both discrete, continuous and mixed probability distributions are supported.+ Continuous probability distributions are represented+ in terms of piecewise polynomials.++ Also includes an implementation of polynomials in one variable.++category: Probability, Math, Numeric, DeltaQ+homepage: https://github.com/DeltaQ-SD/deltaq+license: BSD-3-Clause+license-file: LICENSE+copyright: Predictable Network Solutions Ltd., 2020-2024+author: Peter W. Thompson, Heinrich Apfelmus+maintainer: peter.thompson@pnsol.com++extra-doc-files:+ CHANGELOG.md+ README.md++tested-with:+ , GHC == 9.10.1++common warnings+ ghc-options: -Wall++source-repository head+ type: git+ location: git://github.com/DeltaQ-SD/deltaq.git+ subdir: lib/probability-polynomial++library+ import: warnings+ hs-source-dirs: src+ default-language: Haskell2010++ build-depends:+ , base >= 4.14.3.0 && < 5+ , containers >= 0.6 && < 0.8+ , deepseq >= 1.4.4.0 && < 1.6+ , exact-combinatorics >= 0.2 && < 0.3++ exposed-modules:+ Data.Function.Class+ Numeric.Function.Piecewise+ Numeric.Measure.Discrete+ Numeric.Measure.Finite.Mixed+ Numeric.Measure.Probability+ Numeric.Polynomial.Simple+ Numeric.Probability.Moments++test-suite test+ import: warnings+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ default-language: Haskell2010++ build-tool-depends: hspec-discover:hspec-discover+ + build-depends:+ , base+ , containers+ , probability-polynomial+ , hspec >= 2.11.0 && < 2.12+ , QuickCheck >= 2.14 && < 2.16+ + main-is:+ Spec.hs+ + other-modules:+ Numeric.Function.PiecewiseSpec+ Numeric.Measure.DiscreteSpec+ Numeric.Measure.Finite.MixedSpec+ Numeric.Measure.ProbabilitySpec+ Numeric.Polynomial.SimpleSpec++benchmark probability-polynomial-benchmark+ import: warnings+ type: exitcode-stdio-1.0+ hs-source-dirs: benchmark+ default-language: Haskell2010+ main-is: Main.hs++ build-depends:+ , base+ , probability-polynomial+ , criterion >= 1.6 && < 1.7+ , deepseq
+ src/Data/Function/Class.hs view
@@ -0,0 +1,58 @@+{-# LANGUAGE TypeFamilies #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+Description : Type class for functions, e.g. polynomials.+-}+module Data.Function.Class+ ( Function (..)+ ) where++import qualified Data.Map as Map+import qualified Data.Set as Set++-- | An instance of 'Function' is a type that represents functions.+-- Function can be evaluated at points in their 'Domain'.+--+-- Examples: Polynomials, trigonometric polynomials, piecewise polynomials, …+class Function f where+ -- | The __domain__ of definition of the function.+ type Domain f+ -- | The __codomain__ of a function is the set of potential function values,+ -- i.e. function values never lie outside this set.+ --+ -- In contrast, the set of actual function values+ -- is called the __image__ and+ -- is typically a strict subset of the codomain.+ type Codomain f++ -- | Evaluate a function at a point in its 'Domain'.+ eval :: f -> Domain f -> Codomain f++-- | Functions are 'Function'.+instance Function (a -> b) where+ type Domain (a -> b) = a+ type Codomain (a -> b) = b++ eval = id++-- | @'Map.Map' k v@ represents a function @k -> Maybe v@.+--+-- > Domain (Map k v) = k+-- > Codomain (Map k v) = Maybe v+instance Ord k => Function (Map.Map k v) where+ type instance Domain (Map.Map k v) = k+ type instance Codomain (Map.Map k v) = Maybe v++ eval = flip Map.lookup++-- | @'Set.Set' v@ represents a function @v -> Bool@.+--+-- > Domain (Set v) = v+-- > Codomain (Set v) = Bool+instance Ord v => Function (Set.Set v) where+ type Domain (Set.Set v) = v+ type Codomain (Set.Set v) = Bool++ eval = flip Set.member
+ src/Numeric/Function/Piecewise.hs view
@@ -0,0 +1,268 @@+{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE UndecidableInstances #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+Description : Piecewise functions on the number line. +-}+module Numeric.Function.Piecewise+ ( -- * Type+ Piecewise++ -- * Basic operations+ , zero+ , fromInterval+ , fromAscPieces+ , toAscPieces+ , intervals++ -- * Structure+ , mapPieces+ , mergeBy+ , trim++ -- * Numerical+ , evaluate+ , translateWith++ -- * Zip+ , zipPointwise+ ) where++import Control.DeepSeq+ ( NFData+ )+import GHC.Generics+ ( Generic+ )++import qualified Data.Function.Class as Fun++{-----------------------------------------------------------------------------+ Type+------------------------------------------------------------------------------}+-- | Internal representation of a single piece,+-- starting at a basepoint of type @a@+-- and containing an object of type @o@.+data Piece a o = Piece+ { basepoint :: a+ , object :: o+ }+ deriving (Eq, Show, Generic, NFData)++{- | A function defined piecewise on numerical intervals.+ +* @o@ = type of function on every piece+ e.g. polynomials or other specialized representations of functions+* @'Fun.Domain' o@ = numerical type for the number line, e.g. 'Rational' or 'Double'++A value @f :: Piecewise o@ represents a function++> eval f x = { 0 if -∞ < x < x1+> { eval o1 x if x1 <= x < x2+> { eval o2 x if x2 <= x < x3+> { …+> { eval on x if xn <= x < +∞++where @x1, …, xn@ are points on the real number line+(in strictly increasing order)+and where @o1, …, on@ are specialized representations functions,+e.g. polynomials.++In other words, the value @f@ represents a function that+is defined piecewise on half-open intervals.++The function 'intervals' returns the half-open intervals in the middle:++> intervals f = [(x1,x2), (x2,x3), …, (xn-1, xn)]++No attempt is made to merge intervals if the piecewise objects are equal,+e.g. the situation @o1 == o2@ may occur.++-}+data Piecewise o+ = Pieces [Piece (Fun.Domain o) o]+ deriving (Generic)++deriving instance (Show (Fun.Domain o), Show o) => Show (Piecewise o)+deriving instance (NFData (Fun.Domain o), NFData o) => NFData (Piecewise o)++{-$Piecewise Invariants++* The empty list represents the zero function.+* The 'basepoint's are in strictly increasing order.+* The internal representation of the function mentioned in the definition is++ > f = Pieces [Piece x1 o1, Piece x2 o2, …, Piece xn on]+-}++{-----------------------------------------------------------------------------+ Operations+------------------------------------------------------------------------------}+-- | The function which is zero everywhere.+zero :: Piecewise o+zero = Pieces []++-- | @fromInterval (x1,x2) o@ creates a 'Piecewise' function+-- from a single function @o@ by restricting it to the+-- to half-open interval @x1 <= x < x2@.+-- The result is zero outside this interval.+fromInterval+ :: (Ord (Fun.Domain o), Num o)+ => (Fun.Domain o, Fun.Domain o) -> o -> Piecewise o+fromInterval (x,y) o = Pieces [Piece start o, Piece end 0]+ where+ start = min x y+ end = max x y++-- | Build a piecewise function from an ascending list of contiguous pieces.+--+-- /The precondition (`map fst` of input list is ascending) is not checked./+fromAscPieces :: Ord (Fun.Domain o) => [(Fun.Domain o, o)] -> Piecewise o+fromAscPieces = Pieces . map (uncurry Piece)++-- | Convert the piecewise function to a list of contiguous pieces+-- where the starting points of the pieces are in ascending order.+toAscPieces :: Ord (Fun.Domain o) => Piecewise o -> [(Fun.Domain o, o)]+toAscPieces (Pieces xos) = [ (x, o) | Piece x o <- xos ]++-- | Intervals on which the piecewise function is defined, in sequence.+-- The last half-open interval, @xn <= x < +∞@, is omitted.+intervals :: Piecewise o -> [(Fun.Domain o, Fun.Domain o)]+intervals (Pieces ys) =+ zip (map basepoint ys) (drop 1 $ map basepoint ys)++{-----------------------------------------------------------------------------+ Operations+ Structure+------------------------------------------------------------------------------}+-- | Map the objects of pieces.+mapPieces+ :: Fun.Domain o ~ Fun.Domain o'+ => (o -> o') -> Piecewise o -> Piecewise o'+mapPieces f (Pieces ps) = Pieces [ Piece x (f o) | Piece x o <- ps ]++-- | Merge all adjacent pieces whose functions are considered+-- equal by the given predicate.+mergeBy :: Num o => (o -> o -> Bool) -> Piecewise o -> Piecewise o+mergeBy eq (Pieces pieces) = Pieces $ go 0 pieces+ where+ go _ [] = []+ go before (p : ps)+ | before `eq` object p = go before ps+ | otherwise = p : go (object p) ps++-- | Merge all adjacent pieces whose functions are equal according to '(==)'.+trim :: (Eq o, Num o) => Piecewise o -> Piecewise o+trim = mergeBy (==)++{-----------------------------------------------------------------------------+ Operations+ Evaluation+------------------------------------------------------------------------------}+{-|+Evaluate a piecewise function at a point.++* @'Fun.Domain' ('Piecewise' o) = 'Fun.Domain' o@+* @'Fun.Codomain' ('Piecewise' o) = 'Fun.Codomain' o@+-}+instance (Fun.Function o, Num o, Ord (Fun.Domain o), Num (Fun.Codomain o))+ => Fun.Function (Piecewise o)+ where+ type instance Domain (Piecewise o) = Fun.Domain o+ type instance Codomain (Piecewise o) = Fun.Codomain o+ eval = evaluate++-- | Evaluate the piecewise function at a point.+-- See 'Piecewise' for the semantics.+evaluate+ :: (Fun.Function o, Num o, Ord (Fun.Domain o), Num (Fun.Codomain o))+ => Piecewise o -> Fun.Domain o -> Fun.Codomain o+evaluate (Pieces pieces) x = go 0 pieces+ where+ go before [] = Fun.eval before x+ go before (p:ps)+ | basepoint p <= x = go (object p) ps+ | otherwise = Fun.eval before x++-- | Translate a piecewise function,+-- given a way to translate each piece.+--+-- > eval (translate' y o) = eval o (x - y)+-- > implies+-- > eval (translateWith translate' y p) = eval p (x - y)+translateWith+ :: (Ord (Fun.Domain o), Num (Fun.Domain o), Num o)+ => (Fun.Domain o -> o -> o)+ -> Fun.Domain o -> Piecewise o -> Piecewise o+translateWith trans y (Pieces pieces) =+ Pieces [ Piece (x + y) (trans y o) | Piece x o <- pieces ]++{-----------------------------------------------------------------------------+ Operations+ Zip+------------------------------------------------------------------------------}+-- | Combine two piecewise functions by combining the pieces+-- with a pointwise operation that preserves @0@.+--+-- For example, `(+)` and `(*)` are pointwise operations on functions,+-- but convolution is not a pointwise operation.+--+-- Preconditions on the argument @f@:+--+-- * @f 0 0 = 0@+-- * @f@ is a pointwise operations on functions,+-- e.g. commutes with pointwise evaluation.+--+-- /The preconditions are not checked!/+zipPointwise+ :: (Ord (Fun.Domain o), Num o)+ => (o -> o -> o)+ -- ^ @f@+ -> Piecewise o -> Piecewise o -> Piecewise o+zipPointwise f (Pieces xs') (Pieces ys') =+ Pieces $ go 0 xs' 0 ys'+ where+ -- We split the intervals and combine the pieces in a single pass.+ --+ -- The algorithm is similar to mergesort:+ -- We walk both lists in parallel and generate a new piece by+ -- * taking the basepoint of the nearest piece+ -- * and combining it with the object that was overhanging from+ -- the previous piece (`xhang`, `yhang`)+ go _ [] _ [] = []+ go _ (Piece x ox : xstail) yhang [] =+ Piece x (f ox yhang) : go ox xstail yhang []+ go xhang [] _ (Piece y oy : ystail) =+ Piece y (f xhang oy) : go xhang [] oy ystail+ go xhang xs@(Piece x ox : xstail) yhang ys@(Piece y oy : ystail) =+ case compare x y of+ LT -> Piece x (f ox yhang) : go ox xstail yhang ys+ EQ -> Piece x (f ox oy ) : go ox xstail oy ystail+ GT -> Piece y (f xhang oy ) : go xhang xs oy ystail++{-----------------------------------------------------------------------------+ Operations+ Numeric+------------------------------------------------------------------------------}+{-| Algebraic operations '(+)', '(*)' and 'negate' on piecewise functions.++The functions 'abs' and 'signum' are defined using 'abs' and 'signum'+for every piece.++TODO: 'fromInteger' is __undefined__+-}+instance (Ord (Fun.Domain o), Num o) => Num (Piecewise o) where+ (+) = zipPointwise (+)+ (*) = zipPointwise (*)+ negate = mapPieces negate+ abs = mapPieces abs+ signum = mapPieces signum+ fromInteger 0 = zero+ fromInteger _ = error "TODO: fromInteger not implemented"
+ src/Numeric/Measure/Discrete.hs view
@@ -0,0 +1,149 @@+{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+Description : Discrete, finite signed measures on the number line.+-}+module Numeric.Measure.Discrete+ ( -- * Type+ Discrete+ , fromMap+ , toMap+ , zero+ , dirac+ , distribution++ -- * Observations+ , total+ , integrate++ -- * Operations, numerical+ , add+ , scale+ , translate+ , convolve+ ) where++import Data.List+ ( scanl'+ )+import Data.Map+ ( Map+ )+import Numeric.Function.Piecewise+ ( Piecewise+ )+import Numeric.Polynomial.Simple+ ( Poly+ )++import qualified Data.Map.Strict as Map+import qualified Numeric.Function.Piecewise as Piecewise+import qualified Numeric.Polynomial.Simple as Poly++{-----------------------------------------------------------------------------+ Type+------------------------------------------------------------------------------}+-- | A discrete, finite+-- [signed measure](https://en.wikipedia.org/wiki/Signed_measure)+-- on the number line.+newtype Discrete a = Discrete (Map a a)+ -- INVARIANT: All values are non-zero.+ deriving (Show)++-- | Internal.+-- Remove all zero values.+trim :: (Ord a, Num a) => Map a a -> Map a a+trim m = Map.filter (/= 0) m++-- | Two measures are equal if they yield the same measures on every set.+--+-- > mx == my+-- > implies+-- > forall t. eval (distribution mx) t = eval (distribution my) t+instance (Ord a, Num a) => Eq (Discrete a) where+ Discrete mx == Discrete my = mx == my++{-----------------------------------------------------------------------------+ Operations+------------------------------------------------------------------------------}+-- | The measure that assigns @0@ to every set.+zero :: Num a => Discrete a+zero = Discrete Map.empty++-- | A+-- [Dirac measure](https://en.wikipedia.org/wiki/Dirac_measure)+-- at the given point @x@.+--+-- > total (dirac x) = 1+dirac :: (Ord a, Num a) => a -> Discrete a+dirac x = Discrete (Map.singleton x 1)++-- | Construct a discrete measure+-- from a collection of points and their measures.+fromMap :: (Ord a, Num a) => Map a a -> Discrete a+fromMap = Discrete . trim++-- | Decompose the discrete measure into a collection of points+-- and their measures.+toMap :: Num a => Discrete a -> Map a a+toMap (Discrete m) = m++-- | The total of the measure applied to the set of real numbers.+total :: Num a => Discrete a -> a+total (Discrete m) = sum m++-- | Integrate a function @f@ with respect to the given measure @m@,+-- \( \int f(x) dm(x) \).+integrate :: (Ord a, Num a) => (a -> a) -> Discrete a -> a+integrate f (Discrete m) = sum $ Map.mapWithKey (\x w -> f x * w) m++-- | @eval (distribution m) x@ is the measure of the interval \( (-∞, x] \).+--+-- This is known as the [distribution function+-- ](https://en.wikipedia.org/wiki/Distribution_function_%28measure_theory%29).+distribution :: (Ord a, Num a) => Discrete a -> Piecewise (Poly a)+distribution (Discrete m) =+ Piecewise.fromAscPieces+ $ zipWith (\(x,_) s -> (x,Poly.constant s)) diracs steps+ where+ diracs = Map.toAscList m+ steps = tail $ scanl' (+) 0 $ map snd diracs++-- | Add two measures.+--+-- > total (add mx my) = total mx + total my+add :: (Ord a, Num a) => Discrete a -> Discrete a -> Discrete a+add (Discrete mx) (Discrete my) =+ Discrete $ trim $ Map.unionWith (+) mx my++-- | Scale a measure by a constant.+--+-- > total (scale a mx) = a * total mx+scale :: (Ord a, Num a) => a -> Discrete a -> Discrete a+scale 0 (Discrete _) = Discrete Map.empty+scale s (Discrete m) = Discrete $ Map.map (s *) m++-- | Translate a measure along the number line.+--+-- > eval (distribution (translate y m)) x+-- > = eval (distribution m) (x - y)+translate :: (Ord a, Num a) => a -> Discrete a -> Discrete a+translate y (Discrete m) = Discrete $ Map.mapKeys (y +) m++-- | Additive convolution of two measures.+--+-- Properties:+--+-- > convolve (dirac x) (dirac y) = dirac (x + y)+convolve :: (Ord a, Num a) => Discrete a -> Discrete a -> Discrete a+-- >+-- > convolve mx my = convolve my mx+-- > convolve (add mx my) mz = add (convolve mx mz) (convolve my mz)+-- > translate z (convolve mx my) = convolve (translate z mx) my+-- > total (convolve mx my) = total mx * total myconvolve :: (Ord a, Num a) => Discrete a -> Discrete a -> Discrete a+convolve (Discrete mx) (Discrete my) =+ Discrete $ trim $ Map.fromListWith (+)+ [ (x + y, wx * wy)+ | (x,wx) <- Map.toList mx+ , (y,wy) <- Map.toList my+ ]
+ src/Numeric/Measure/Finite/Mixed.hs view
@@ -0,0 +1,381 @@+{-# LANGUAGE GeneralizedNewtypeDeriving #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+Description : Finite signed measures on the number line.+-}+module Numeric.Measure.Finite.Mixed+ ( -- * Type+ Measure+ , zero+ , dirac+ , uniform+ , distribution+ , fromDistribution++ -- * Observations+ , total+ , support+ , isPositive+ , integrate++ -- * Operations, numerical+ , add+ , scale+ , translate+ , convolve+ ) where++import Data.Function.Class+ ( Function (..)+ )+import Data.List+ ( scanl'+ )+import Control.DeepSeq+ ( NFData+ )+import Numeric.Function.Piecewise+ ( Piecewise+ )+import Numeric.Polynomial.Simple+ ( Poly+ )++import qualified Data.Map.Strict as Map+import qualified Numeric.Function.Piecewise as Piecewise+import qualified Numeric.Measure.Discrete as D+import qualified Numeric.Polynomial.Simple as Poly++{-----------------------------------------------------------------------------+ Type+------------------------------------------------------------------------------}+-- | A finite+-- [signed measure](https://en.wikipedia.org/wiki/Signed_measure)+-- on the number line.+newtype Measure a = Measure (Piecewise (Poly a))+ -- INVARIANT: Adjacent pieces contain distinct objects.+ -- INVARIANT: The last piece is a constant polynomial,+ -- so that the measure is finite.+ deriving (Show, NFData)++-- | @eval (distribution m) x@ is the measure of the interval \( (-∞, x] \).+--+-- This is known as the [distribution function+-- ](https://en.wikipedia.org/wiki/Distribution_function_%28measure_theory%29).+distribution :: (Ord a, Num a) => Measure a -> Piecewise (Poly a)+distribution (Measure p) = p++-- | Construct a signed measure from its+-- [distribution function+-- ](https://en.wikipedia.org/wiki/Distribution_function_%28measure_theory%29).+--+-- Return 'Nothing' if the measure is not finite,+-- that is if the last piece of the piecewise function is not constant.+fromDistribution+ :: (Ord a, Num a)+ => Piecewise (Poly a) -> Maybe (Measure a)+fromDistribution pieces+ | isEventuallyConstant pieces = Just $ Measure $ trim pieces+ | otherwise = Nothing++-- | Test whether a piecewise polynomial is consant as x -> ∞.+isEventuallyConstant :: (Ord a, Num a) => Piecewise (Poly a) -> Bool+isEventuallyConstant pieces+ | null xpolys = True+ | otherwise = (<= 0) . Poly.degree . snd $ last xpolys+ where+ xpolys = Piecewise.toAscPieces pieces++-- | Internal.+-- Join all intervals whose polynomials are equal.+trim :: (Ord a, Num a) => Piecewise (Poly a) -> Piecewise (Poly a)+trim = Piecewise.trim++-- | Two measures are equal if they yield the same measures on every set.+--+-- > mx == my+-- > implies+-- > forall t. eval (distribution mx) t = eval (distribution my) t+instance (Ord a, Num a) => Eq (Measure a) where+ Measure mx == Measure my =+ Piecewise.toAscPieces mx == Piecewise.toAscPieces my++{-----------------------------------------------------------------------------+ Operations+------------------------------------------------------------------------------}+-- | The measure that assigns @0@ to every set.+zero :: Num a => Measure a+zero = Measure Piecewise.zero++-- | A+-- [Dirac measure](https://en.wikipedia.org/wiki/Dirac_measure)+-- at the given point @x@.+--+-- > total (dirac x) = 1+dirac :: (Ord a, Num a) => a -> Measure a+dirac x = Measure $ Piecewise.fromAscPieces [(x, Poly.constant 1)]++-- | The probability measure of a uniform probability distribution+-- in the interval \( [x,y) \).+--+-- > total (uniform x y) = 1+uniform :: (Ord a, Num a, Fractional a) => a -> a -> Measure a+uniform x y = Measure $ case compare x y of+ EQ -> Piecewise.fromAscPieces [(x, 1)]+ _ -> Piecewise.fromAscPieces [(low, poly), (high, 1)]+ where+ low = min x y+ high = max x y+ poly = Poly.lineFromTo (low, 0) (high, 1)++-- | The total of the measure applied to the set of real numbers.+total :: (Ord a, Num a) => Measure a -> a+total (Measure p) =+ case Piecewise.toAscPieces p of+ [] -> 0+ ps -> eval (snd (last ps)) 0++-- | The 'support' is the smallest closed, contiguous interval \( [x,y] \)+-- outside of which the measure is zero.+--+-- Returns 'Nothing' if the interval is empty.+support :: (Ord a, Num a) => Measure a -> Maybe (a, a)+support (Measure pieces) =+ case Piecewise.toAscPieces pieces of+ [] -> Nothing+ ps -> Just (fst $ head ps, fst $ last ps)++-- | Check whether a signed measure is positive.+--+-- A signed measure is /positive/ if the measure of any set+-- is nonnegative. In other words a positive signed measure+-- is just a measure in the ordinary sense.+--+-- This test is nontrivial, as we have to check that the distribution+-- function is monotonically increasing.+isPositive :: (Ord a, Num a, Fractional a) => Measure a -> Bool+isPositive (Measure m) = go 0 $ Piecewise.toAscPieces m+ where+ go _ [] =+ True+ go before ((x, o) : []) =+ eval before x <= eval o x+ go before ((x1, o) : xos@((x2, _) : _)) =+ (eval before x1 <= eval o x1)+ && Poly.isMonotonicallyIncreasingOn o (x1,x2)+ && go o xos++{-----------------------------------------------------------------------------+ Operations+ Numerical+------------------------------------------------------------------------------}+-- | Add two measures.+--+-- > total (add mx my) = total mx + total my+add :: (Ord a, Num a) => Measure a -> Measure a -> Measure a+add (Measure mx) (Measure my) =+ Measure $ trim $ Piecewise.zipPointwise (+) mx my++-- | Scale a measure by a constant.+--+-- > total (scale a mx) = a * total mx+scale :: (Ord a, Num a) => a -> Measure a -> Measure a+scale 0 (Measure _) = zero+scale x (Measure m) = Measure $ Piecewise.mapPieces (Poly.scale x) m++-- | Translate a measure along the number line.+--+-- > eval (distribution (translate y m)) x+-- > = eval (distribution m) (x - y)+translate :: (Ord a, Num a, Fractional a) => a -> Measure a -> Measure a+translate y (Measure m) =+ Measure $ Piecewise.translateWith Poly.translate y m++{-----------------------------------------------------------------------------+ Operations+ Decomposition into continuous and discrete measures,+ needed for convolution.+------------------------------------------------------------------------------}+-- | Measure that is absolutely continuous+-- with respect to the Lebesgue measure,+-- Represented via its distribution function.+newtype Continuous a = Continuous { unContinuous :: Piecewise (Poly a) }+ -- INVARIANT: The last piece is @Poly.constant p@ for some @p :: a@.++-- | Density function (Radon–Nikodym derivative) of an absolutely+-- continuous measure.+newtype Density a = Density (Piecewise (Poly a))+ -- INVARIANT: The last piece is @Poly.constant 0@.++-- | Density function of an absolutely continuous measure.+toDensity+ :: (Ord a, Num a, Fractional a)+ => Continuous a -> Density a+toDensity = Density . Piecewise.mapPieces Poly.differentiate . unContinuous++-- | Decompose a mixed measure into+-- a continuous measure and a discrete measure.+-- See also [Lebesgue's decomposition theorem+-- ](https://en.wikipedia.org/wiki/Lebesgue%27s_decomposition_theorem)+decompose+ :: (Ord a, Num a, Fractional a)+ => Measure a -> (Continuous a, D.Discrete a)+decompose (Measure m) =+ ( Continuous $ trim $ Piecewise.fromAscPieces withoutJumps+ , D.fromMap $ Map.fromList jumps+ )+ where+ pieces = Piecewise.toAscPieces m++ withoutJumps =+ zipWith (\(x,o) j -> (x, o - Poly.constant j)) pieces totalJumps+ totalJumps = tail $ scanl' (+) 0 $ map snd jumps++ jumps = go 0 pieces+ where+ go _ [] = []+ go prev ((x,o) : xos) =+ (x, Poly.eval o x - Poly.eval prev x) : go o xos++{-----------------------------------------------------------------------------+ Observations+ Integration+------------------------------------------------------------------------------}+-- | Integrate a polynomial @f@ with respect to the given measure @m@,+-- \( \int f(x) dm(x) \).+integrate :: (Ord a, Num a, Fractional a) => Poly a -> Measure a -> a+integrate f m =+ integrateContinuous f continuous+ + D.integrate (eval f) discrete+ where+ (continuous, discrete) = decompose m++-- | Integrate a polynomial over an absolutely continuous measure.+integrateContinuous+ :: (Ord a, Num a, Fractional a)+ => Poly a -> Continuous a -> a+integrateContinuous f gg+ | null gpieces = 0+ | otherwise = sum $ map integrateOverInterval $ integrands+ where+ Density g = toDensity gg+ gpieces = Piecewise.toAscPieces g++ -- Pieces on the bounded intervals+ boundedPieces xos =+ zipWith (\(x1,o) (x2,_) -> ((x1, x2), o)) xos (drop 1 xos)++ integrands = [ (x12, f * o) | (x12, o) <- boundedPieces gpieces ]++ integrateOverInterval ((x1, x2), p) =+ eval pp x2 - eval pp x1+ where+ pp = Poly.integrate p++{-----------------------------------------------------------------------------+ Operations+ Convolution+------------------------------------------------------------------------------}+{-$ NOTE [Convolution]++In order to compute a convolution,+we convolve a density with the distribution function.++Let $f$ denote a density, which can be continuous or a Dirac delta.+Let $G$ denote a distribution function.+Let $H = f * G$ be the result of the convolution.+It can be shown that this is the distribution function of the+convolution of the densities, $h = f * g$.++The formula for convolution is++$ H(y) = ∫ f(y - x) G(x) dx = ∫ f (x) G(y - x) dx$.++When $f$ is a sum of delta functions, $f = Σ w_j delta_{x_j}(x)$,+this integral becomes ($y - x = x_j$ => $x = y - x_j$)++$ H(y) = Σ w_j G(y - x_j) $.++When $f$ is a piecewise polynomial, we can convolve the pieces.++When convolving with a distribution function, the final piece+will be a constant $g_n$ on the interval $[x_n,∞)$.+In this case, the convolution is given by++\[+H(y)+ = ∫ f (x) G(y - x) dx+ = ∫_{ -∞}^{y-x_n} f(x) g_n dx+ = g_n F(y-x_n)+\]++where $F$ is the distribution function of the density $f$.+-}++-- | Convolve a discrete measure with a mixed measure.+--+-- See NOTE [Convolution].+convolveDiscrete+ :: (Ord a, Num a, Fractional a)+ => D.Discrete a -> Measure a -> Measure a+convolveDiscrete f gg =+ foldr add zero+ [ scale w (translate x gg)+ | (x, w) <- Map.toAscList $ D.toMap f+ ]++-- | Convolve an absolutely continuous measure with a mixed measure.+--+-- See NOTE [Convolution].+convolveContinuous+ :: (Ord a, Num a, Fractional a)+ => Continuous a -> Measure a -> Measure a+convolveContinuous (Continuous ff) (Measure gg)+ | null ffpieces = zero+ | null ggpieces = zero+ | otherwise = Measure $ trim $ boundedConvolutions + lastConvolution+ where+ ffpieces = Piecewise.toAscPieces ff+ ggpieces = Piecewise.toAscPieces gg++ Density f = toDensity (Continuous ff)+ fpieces = Piecewise.toAscPieces f++ -- Pieces on the bounded intervals+ boundedPieces xos =+ zipWith (\(x,o) (y,_) -> (x, y, o)) xos (drop 1 xos)++ boundedConvolutions =+ sum $+ [ Piecewise.fromAscPieces (Poly.convolve fo ggo)+ | fo <- boundedPieces fpieces+ , ggo <- boundedPieces ggpieces+ ]++ (xlast, plast) = last ggpieces+ glast = case Poly.toCoefficients plast of+ [] -> 0+ (a0:_) -> a0+ lastConvolution =+ Piecewise.mapPieces (Poly.scale glast)+ $ Piecewise.translateWith Poly.translate xlast ff++-- | Additive convolution of two measures.+--+-- Properties:+--+-- > convolve (dirac x) (dirac y) = dirac (x + y)+-- >+-- > convolve mx my = convolve my mx+-- > convolve (add mx my) mz = add (convolve mx mz) (convolve my mz)+-- > translate z (convolve mx my) = convolve (translate z mx) my+-- > total (convolve mx my) = total mx * total my+convolve+ :: (Ord a, Num a, Fractional a)+ => Measure a -> Measure a -> Measure a+convolve mx my =+ add (convolveContinuous contx my) (convolveDiscrete deltax my)+ where+ (contx, deltax) = decompose mx
+ src/Numeric/Measure/Probability.hs view
@@ -0,0 +1,192 @@+{-# LANGUAGE GeneralizedNewtypeDeriving #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+Description : Probability measures on the number line.+-}+module Numeric.Measure.Probability+ ( -- * Type+ Prob+ , dirac+ , uniform+ , distribution+ , fromDistribution+ , measure+ , fromMeasure+ , unsafeFromMeasure++ -- * Observations+ , support+ , expectation+ , moments++ -- * Operations, numerical+ , choice+ , translate+ , convolve+ ) where++import Control.DeepSeq+ ( NFData+ )+import Numeric.Function.Piecewise+ ( Piecewise+ )+import Numeric.Measure.Finite.Mixed+ ( Measure+ )+import Numeric.Polynomial.Simple+ ( Poly+ )+import Numeric.Probability.Moments+ ( Moments (..)+ , fromExpectedPowers+ )++import qualified Numeric.Measure.Finite.Mixed as M+import qualified Numeric.Polynomial.Simple as Poly++{-----------------------------------------------------------------------------+ Type+------------------------------------------------------------------------------}+-- | A+-- [probability measure](https://en.wikipedia.org/wiki/Probability_measure)+-- on the number line.+--+-- A probability measure is a 'M.Measure' whose 'M.total' is @1@.+newtype Prob a = Prob (Measure a)+ -- INVARIANT: 'M.isPositive' equals 'True'.+ -- INVARIANT: 'M.total' equals 1+ deriving (Show, NFData)++-- | View the probability measure as a 'M.Measure'.+measure :: (Ord a, Num a) => Prob a -> Measure a+measure (Prob m) = m++-- | View a 'M.Measure' as a probability distribution.+--+-- The measure @m@ must be positive, with total weight @1@, that is+--+-- > isPositive m == True+-- > total m == 1+--+-- These preconditions are checked and the function returns 'Nothing'+-- if they fail. +fromMeasure :: (Ord a, Num a, Fractional a) => Measure a -> Maybe (Prob a)+fromMeasure m+ | M.isPositive m && M.total m == 1 = Just $ Prob m+ | otherwise = Nothing++-- | View a 'M.Measure' as a probability distribution.+--+-- Variant of 'fromMeasure' where /the precondition are not checked!/+unsafeFromMeasure :: Measure a -> Prob a+unsafeFromMeasure = Prob++-- | @eval (distribution m) x@ is the probability of picking a number @<= x@.+--+-- This is known as the+-- [cumulative distribution function+-- ](https://en.wikipedia.org/wiki/Cumulative_distribution_function).+distribution :: (Ord a, Num a) => Prob a -> Piecewise (Poly a)+distribution (Prob m) = M.distribution m++-- | Construct a probability distribution from its+-- [cumulative distribution function+-- ](https://en.wikipedia.org/wiki/Cumulative_distribution_function).+--+-- Return 'Nothing' if+-- * the cumulative distribution function is not monotonicall increasing+-- * the last piece of the piecewise function is not a constant+-- equal to @1@.+fromDistribution+ :: (Ord a, Num a, Fractional a)+ => Piecewise (Poly a) -> Maybe (Prob a)+fromDistribution pieces+ | Just m <- M.fromDistribution pieces = fromMeasure m+ | otherwise = Nothing++-- | Two probability measures are equal if they have the same cumulative+-- distribution functions.+--+-- > px == py+-- > implies+-- > forall t. eval (distribution px) t = eval (distribution py) t+instance (Ord a, Num a) => Eq (Prob a) where+ Prob mx == Prob my = mx == my++{-----------------------------------------------------------------------------+ Construction+------------------------------------------------------------------------------}+-- | A+-- [Dirac measure](https://en.wikipedia.org/wiki/Dirac_measure)+-- at the given point @x@.+--+-- @dirac x@ is the probability distribution where @x@ occurs with certainty.+dirac :: (Ord a, Num a) => a -> Prob a+dirac = Prob . M.dirac++-- | The uniform probability distribution on the interval \( [x,y) \).+uniform :: (Ord a, Num a, Fractional a) => a -> a -> Prob a+uniform x y = Prob $ M.uniform x y++{-----------------------------------------------------------------------------+ Construction+------------------------------------------------------------------------------}+-- | The 'support' is the smallest closed, contiguous interval \( [x,y] \)+-- outside of which the probability is zero.+--+-- Returns 'Nothing' if the interval is empty.+support :: (Ord a, Num a) => Prob a -> Maybe (a, a)+support (Prob m) = M.support m++-- | Compute the+-- [expected value](https://en.wikipedia.org/wiki/Expected_value)+-- of a polynomial @f@ with respect to the given probability distribution,+-- \( E[f(X)] \).+expectation :: (Ord a, Num a, Fractional a) => Poly a -> Prob a -> a+expectation f (Prob m) = M.integrate f m++-- | Compute the first four+-- commonly used moments of a probability distribution.+moments :: (Ord a, Num a, Fractional a) => Prob a -> Moments a+moments m =+ fromExpectedPowers (ex 1, ex 2, ex 3, ex 4)+ where+ ex n = expectation (Poly.monomial n 1) m++{-----------------------------------------------------------------------------+ Operations+------------------------------------------------------------------------------}+-- | Left-biased random choice.+--+-- @choice p@ is a probability distribution where+-- events from the left argument are chosen with probablity @p@+-- and events from the right argument are chosen with probability @(1-p)@.+--+-- > eval (distribution (choice p mx my)) z+-- > = p * eval (distribution mx) z + (1-p) * eval (distribution my) z+choice :: (Ord a, Num a, Fractional a) => a -> Prob a -> Prob a -> Prob a+choice p (Prob mx) (Prob my) = Prob $+ M.add (M.scale p mx) (M.scale (1 - p) my)++-- | Translate a probability distribution along the number line.+--+-- > eval (distribution (translate y m)) x+-- > = eval (distribution m) (x - y)+translate :: (Ord a, Num a, Fractional a) => a -> Prob a -> Prob a+translate y (Prob m) = Prob $ M.translate y m++-- | Additive convolution of two probability measures.+--+-- Properties:+--+-- > convolve (dirac x) (dirac y) = dirac (x + y)+-- >+-- > convolve mx my = convolve my mx+-- > translate z (convolve mx my) = convolve (translate z mx) my+convolve+ :: (Ord a, Num a, Fractional a)+ => Prob a -> Prob a -> Prob a+convolve (Prob mx) (Prob my) = Prob $ M.convolve mx my
+ src/Numeric/Polynomial/Simple.hs view
@@ -0,0 +1,618 @@+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+Description : Polynomials and computations with them.+-}+module Numeric.Polynomial.Simple+ ( -- * Basic operations+ Poly+ , eval+ , degree+ , constant+ , zero+ , monomial+ , fromCoefficients+ , toCoefficients+ , scale+ , scaleX++ -- * Advanced operations++ -- ** Convenience+ , display+ , lineFromTo++ -- ** Algebraic+ , translate+ , integrate+ , differentiate+ , euclidianDivision+ , convolve++ -- ** Numerical+ , compareToZero+ , countRoots+ , isMonotonicallyIncreasingOn+ , root+ ) where++import Control.DeepSeq+ ( NFData+ , NFData1+ )+import GHC.Generics+ ( Generic+ , Generic1+ )+import Math.Combinatorics.Exact.Binomial -- needed to automatically derive NFData+ ( choose+ )++import qualified Data.Function.Class as Fun++{-----------------------------------------------------------------------------+ Basic operations+------------------------------------------------------------------------------}++-- | Polynomial with coefficients in @a@.+newtype Poly a = Poly [a]+ -- INVARIANT: List of coefficients from lowest to highest degree.+ -- INVARIANT: The empty list is not allowed,+ -- the zero polynomial is represented as [0].+ deriving (Show, Generic, Generic1)++instance NFData a => NFData (Poly a)+instance NFData1 Poly++instance (Eq a, Num a) => Eq (Poly a) where+ x == y =+ toCoefficients (trimPoly x) == toCoefficients (trimPoly y)++{-| The constant polynomial.++> eval (constant a) = const a+-}+constant :: a -> Poly a+constant x = Poly [x]++-- | The zero polynomial.+zero :: Num a => Poly a+zero = constant 0++{-| Degree of a polynomial.++The degree of a constant polynomial is @0@, but+the degree of the zero polynomial is @-1@ for Euclidean division.+-}+degree :: (Eq a, Num a) => Poly a -> Int+degree x = case trimPoly x of+ Poly [0] -> -1+ Poly xs -> length xs - 1++-- | remove top zeroes+trimPoly :: (Eq a, Num a) => Poly a -> Poly a+trimPoly (Poly as) = Poly (reverse $ goTrim $ reverse as)+ where+ goTrim [] = error "Empty polynomial"+ goTrim xss@[_] = xss -- can't use dropWhile as it would remove the last zero+ goTrim xss@(x : xs) = if x == 0 then goTrim xs else xss++-- | @monomial n a@ is the polynomial @a * x^n@.+monomial :: (Eq a, Num a) => Int -> a -> Poly a+monomial n x = if x == 0 then zero else Poly (reverse (x : replicate n 0))++{-| Construct a polynomial @a0 + a1·x + …@ from+its list of coefficients @[a0, a1, …]@.+-}+fromCoefficients :: (Eq a, Num a) => [a] -> Poly a+fromCoefficients [] = zero+fromCoefficients as = trimPoly $ Poly as++{-| List the coefficients @[a0, a1, …]@+of a polynomial @a0 + a1·x + …@.+-}+toCoefficients :: Poly a -> [a]+toCoefficients (Poly as) = as++{-| Multiply the polynomial by the unknown @x@.++> eval (scaleX p) x = x * eval p x+> degree (scaleX p) = 1 + degree p if degree p >= 0+-}+scaleX :: (Eq a, Num a) => Poly a -> Poly a+scaleX (Poly xs)+ | xs == [0] = Poly xs -- don't shift up zero+ | otherwise = Poly (0 : xs)++{-| Scale a polynomial by a scalar.+More efficient than multiplying by a constant polynomial.++> eval (scale a p) x = a * eval p x+-}+scale :: Num a => a -> Poly a -> Poly a+scale x (Poly xs) = Poly (map (* x) xs)++-- Does not agree with naming conventions in `Data.Poly`.++{-|+ Add polynomials by simply adding their coefficients as long as both lists continue.+ When one list runs out we take the tail of the longer list (this prevents us from just using zipWith!).+ Addtion might cancel out the highest order terms, so need to trim just in case.+-}+addPolys :: (Eq a, Num a) => Poly a -> Poly a -> Poly a+addPolys (Poly as) (Poly bs) = trimPoly (Poly (go as bs))+ where+ go [] ys = ys+ go xs [] = xs+ go (x : xs) (y : ys) = (x + y) : go xs ys++{-|+ multiply term-wise and then add (very simple - FFTs might be faster, but not for today)+ (a0 + a1x + a2x^2 + ...) * (b0 + b1x + b2x^2 ...)+ = a0 * (b0 + b1x + b2x^2 +...) + a1x * (b0 + b1x + ...)+ = (a0*b0) + (a0*b1x) + ...+ + (a1*b0x) ++ + ...+ (may be an optimisation to be done by getting the shortest poly in the right place)+-}+mulPolys :: (Eq a, Num a) => Poly a -> Poly a -> Poly a+mulPolys as bs = sum (intermediateSums as bs)+ where+ intermediateSums :: (Eq a, Num a) => Poly a -> Poly a -> [Poly a]+ intermediateSums _ (Poly []) = error "Second polynomial was empty"+ intermediateSums (Poly []) _ = [] -- stop when we exhaust the first list+ -- as we consume the coeffecients of the first list, we shift up the second list to increase the power under consideration+ intermediateSums (Poly (x : xs)) ys =+ scale x ys : intermediateSums (Poly xs) (scaleX ys)++{-| Algebraic operations '(+)', '(*)' and 'negate' on polynomials.++The functions 'abs' and 'signum' are undefined.+-}+instance (Eq a, Num a) => Num (Poly a) where+ (+) = addPolys+ (*) = mulPolys+ negate (Poly a) = Poly (map negate a)+ abs = undefined+ signum = undefined+ fromInteger n = Poly [Prelude.fromInteger n]++{-|+Evaluate a polynomial at a point.++> eval :: Poly a -> a -> a+-}+instance Num a => Fun.Function (Poly a) where+ type instance Domain (Poly a) = a+ type instance Codomain (Poly a) = a+ eval = eval++{-|+Evaluate a polynomial at a point.++> eval :: Poly a -> a -> a++Uses Horner's method to minimise the number of multiplications.++@+a0 + a1·x + a2·x^2 + ... + a{n-1}·x^{n-1} + an·x^n+ = a0 + x·(a1 + x·(a2 + x·(… + x·(a{n-1} + x·an)) ))+@+-}+eval :: Num a => Poly a -> a -> a+eval (Poly as) x = foldr (\ai result -> x * result + ai) 0 as++{-----------------------------------------------------------------------------+ Advanced operations+ Convenience+------------------------------------------------------------------------------}++{-|+Return a list of pairs @(x, eval p x)@ from the graph of the polynomial.+The values @x@ are from the range @(l, u)@ with uniform spacing @s@.++Specifically,++> map fst (display p (l, u) s)+> = [l, l+s, l + 2·s, … , u'] ++ if u' == l then [] else [l]++where @u'@ is the largest number of the form @u' = l + s·k@, @k@ natural,+that still satisfies @u' < l@.+We always display the last point as well.+-}+display :: (Ord a, Eq a, Num a) => Poly a -> (a, a) -> a -> [(a, a)]+display p (l, u) s+ | s == 0 = map evalPoint [l, u]+ | otherwise = map evalPoint (l : go (l + s))+ where+ evalPoint x = (x, eval p x)+ go x+ | x >= u = [u] -- always include the last point+ | otherwise = x : go (x + s)++{-| Linear polymonial connecting the points @(x1, y1)@ and @(x2, y2)@,+assuming that @x1 ≠ x2@.++If the points are equal, we return a constant polynomial.++> let p = lineFromTo (x1, y1) (x2, y2)+>+> degree p <= 1+> eval p x1 = y1+> eval p x2 = y2+-}+lineFromTo :: (Eq a, Fractional a) => (a, a) -> (a, a) -> Poly a+lineFromTo (x1, y1) (x2, y2)+ | x1 == x2 = constant y1+ | slope == 0 = constant y1+ | otherwise = fromCoefficients [shift, slope]+ where+ -- slope of the linear function+ slope = (y2 - y1) / (x2 - x1)+ -- the constant shift is fixed by+ -- the fact that the line needs to pass through (x1,y1)+ shift = y1 - x1 * slope++{-----------------------------------------------------------------------------+ Advanced operations+ Algebraic+------------------------------------------------------------------------------}++{-| Indefinite integral of a polynomial with constant term zero.++The integral of @x^n@ is @1/(n+1)·x^(n+1)@.++> eval (integrate p) 0 = 0+> integrate (differentiate p) = p - constant (eval p 0)+-}+integrate :: (Eq a, Fractional a) => Poly a -> Poly a+integrate (Poly as) =+ -- Integrate by puting a zero constant term at the bottom and+ -- converting a x^n into a/(n+1) x^(n+1).+ -- 0 -> 0x is the first non-constant term, so we start at 1.+ -- When integrating a zero polynomial with a zero constant+ -- we get [0,0] so need to trim+ trimPoly (Poly (0 : zipWith (/) as (iterate (+ 1) 1)))++{-| Differentiate a polynomial.++We have @dx^n/dx = n·x^(n-1)@.++> differentiate (integrate p) = p+> differentiate (p * q) = (differentiate p) * q + p * (differentiate q)+-}+differentiate :: Num a => Poly a -> Poly a+differentiate (Poly []) = error "Polynomial was empty"+differentiate (Poly [_]) = zero -- constant differentiates to zero+differentiate (Poly (_ : as)) =+ -- discard the constant term, everything else noves down one+ Poly (zipWith (*) as (iterate (+ 1) 1))++{-| Convolution of two polynomials defined on bounded intervals.+Produces three contiguous pieces as a result.+-}+convolve+ :: (Fractional a, Eq a, Ord a) => (a, a, Poly a) -> (a, a, Poly a) -> [(a, Poly a)]+convolve (lf, uf, Poly fs) (lg, ug, Poly gs)+ | (lf < 0) || (lg < 0) = error "Interval bounds cannot be negative"+ | (lf >= uf) || (lg >= ug) = error "Invalid interval" -- upper bounds should be strictly greater than lower bounds+ | (ug - lg) > (uf - lf) = convolve (lg, ug, Poly gs) (lf, uf, Poly fs) -- if g is wider than f, swap the terms+ | otherwise -- we know g is narrower than f+ =+ let+ -- sum a set of terms depending on an iterator k (assumed to go down to 0), where each term is a k-dependent+ -- polynomial with a k-dependent multiplier+ sumSeries k mulFactor poly = sum [mulFactor n `scale` poly n | n <- [0 .. k]]++ -- the inner summation has a similar structure each time+ innerSum m n term k = sumSeries (m + k + 1) innerMult (\j -> monomial (m + n + 1 - j) (term j))+ where+ innerMult j =+ fromIntegral+ (if even j then (m + k + 1) `choose` j else negate ((m + k + 1) `choose` j))++ convolveMonomials m n innerTerm = sumSeries n (multiplier m n) (innerTerm m n)+ where+ multiplier p q k =+ fromIntegral (if even k then q `choose` k else negate (q `choose` k))+ / fromIntegral (p + k + 1)++ {-+ For each term, clock through the powers of each polynomial to give convolutions of monomials, which we sum.+ We extract each coefficient of each polynomial, together with an integer recording their position (i.e. power of x),+ and multiply the coefficients together with the new polynomial generated by convolving the monomials.+ -}+ makeTerm f =+ sum+ [ (a * b) `scale` convolveMonomials m n f+ | (m, a) <- zip [0 ..] fs+ , (n, b) <- zip [0 ..] gs+ ]++ firstTerm =+ makeTerm (\m n k -> innerSum m n (lg ^) k - monomial (n - k) (lf ^ (m + k + 1)))++ secondTerm = makeTerm (\m n -> innerSum m n (\k -> lg ^ k - ug ^ k))++ thirdTerm =+ makeTerm (\m n k -> monomial (n - k) (uf ^ (m + k + 1)) - innerSum m n (ug ^) k)+ in+ {-+ When convolving distributions, both distributions will start at 0 and so there will always be a pair of intervals+ with lg = lf = 0, so we don't need to add an initial zero piece.+ We must have lf + lg < lf + ug due to initial interval validity check. However, it's possible that lf + ug = uf + lg, so+ we need to test for a redundant middle interval+ -}+ if lf + ug == uf + lg+ then [(lf + lg, firstTerm), (uf + lg, thirdTerm), (uf + ug, zero)]+ else+ [ (lf + lg, firstTerm)+ , (lf + ug, secondTerm)+ , (uf + lg, thirdTerm)+ , (uf + ug, zero)+ ]++{-| Translate the argument of a polynomial by summing binomial expansions.++> eval (translate y p) x = eval p (x - y)+-}+translate :: forall a. (Fractional a, Eq a, Num a) => a -> Poly a -> Poly a+translate y (Poly ps) =+ sum+ [ b `scale` binomialExpansion n+ | (n, b) <- zip [0 ..] ps+ ]+ where+ -- binomialTerm n k = coefficient of x^k in the expensation of (x - y)^n+ binomialTerm :: Integer -> Integer -> a+ binomialTerm n k = fromInteger (n `choose` k) * (-y) ^ (n - k)++ -- binomialExpansion n = (x - y)^n expanded as a polyonial in x+ binomialExpansion :: Integer -> Poly a+ binomialExpansion n = Poly (map (binomialTerm n) [0 .. n])++{-|+[Euclidian division of polynomials+](https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclidean_division)+takes two polynomials @a@ and @b ≠ 0@,+and returns two polynomials, the quotient @q@ and the remainder @r@,+such that++> a = q * b + r+> degree r < degree b+-}+euclidianDivision+ :: forall a. (Fractional a, Eq a, Ord a)+ => Poly a -> Poly a -> (Poly a, Poly a)+euclidianDivision pa pb+ | pb == zero = error "Division by zero polynomial"+ | otherwise = goDivide (zero, pa)+ where+ degB = degree pb++ -- Coefficient of the highest power term+ leadingCoefficient :: Poly a -> a+ leadingCoefficient (Poly x) = last x++ lcB = leadingCoefficient pb++ goDivide :: (Poly a, Poly a) -> (Poly a, Poly a)+ goDivide (q, r)+ | degree r < degB = (q, r)+ | otherwise = goDivide (q + s, r - s * pb)+ where+ s = monomial (degree r - degB) (leadingCoefficient r / lcB)++{-----------------------------------------------------------------------------+ Advanced operations+ Numerical+------------------------------------------------------------------------------}+{-|+@'countRoots' (x1, x2, p)@ returns the number of /distinct/ real roots+of the polynomial on the open interval \( (x_1, x_2) \).++(Roots with higher multiplicity are each counted as a single distinct root.)++This function uses [Sturm's theorem+](https://en.wikipedia.org/wiki/Sturm%27s_theorem),+with special provisions for roots on the boundary of the interval.+-}+countRoots :: (Fractional a, Ord a) => (a, a, Poly a) -> Int+countRoots (l, r, p) =+ countRoots' $ (p `factorOutRoot` l) `factorOutRoot` r+ where+ -- we can now assume that the polynomial has no roots at the boundary+ countRoots' q = case degree q of+ -- q is the zero polynomial, so it doesn't *cross* zero+ -1 -> 0+ -- q is a non-zero constant polynomial - no root+ 0 -> 0+ -- q is a linear polynomial,+ 1 -> if eval q l * eval q r < 0 then 1 else 0+ -- q has degree 2 or more so we can construct the Sturm sequence+ _ -> countRootsSturm (l, r, q)++-- | Given a polynomial \( p(x) \) and a value \( a \),+-- this functions factors out the polynomial \( (x-a)^m \),+-- where \( m \) is the highest power where this polynomial+-- divides \( p(x) \) without remainder.+--+-- * If the value \( a \) is a root of the polynomial,+-- then \( m \) is the multiplicity of the root.+-- * If the value \( a \) is not a root, then+-- \( m = 0 \) and the function returns \( p (x) \).+--+-- In other words, this function returns a polynomial \( q (x) \)+-- such that+--+-- \( p(x) = q(x)·(x - a)^m \)+--+-- where \( q(a) ≠ 0 \).+-- If the polynomial \( p(x) \) is identically 'zero',+-- we return 'zero' as well.+factorOutRoot :: (Fractional a, Ord a) => Poly a -> a -> Poly a+factorOutRoot p0 x0+ | p0 == zero = zero+ | otherwise = go p0+ where+ go p+ | eval p x0 == 0 = factorOutRoot pDividedByXMinusX0 x0+ | otherwise = p+ where+ xMinusX0 = monomial 1 1 - constant x0+ (pDividedByXMinusX0, _) = p `euclidianDivision` xMinusX0++{-|+@'countRootsSturm' (x1, x2, p)@ returns the number of /distinct/ real roots+of the polynomial @p@ on the half-open interval \( (x_1, x_2] \),+under the following assumptions:++* @'degree' p >= 2@+* neither \( x_1 \) nor \( x_2 \) are multiple roots of \( p(x) \).++This function is an implementation of [Sturm's theorem+](https://en.wikipedia.org/wiki/Sturm%27s_theorem).+-}+countRootsSturm :: (Fractional a, Eq a, Ord a) => (a, a, Poly a) -> Int+countRootsSturm (l, r, p) =+ -- p has degree 2 or more so we can construct the Sturm sequence+ signVariations psl - signVariations psr+ where+ ps = reversedSturmSequence p+ psl = map (flip eval l) ps+ psr = map (flip eval r) ps++{-| Number of sign variations in a list of real numbers.++Given a list @c0, c1, c2, . . . ck@,+then a sign variation (or sign change) in the sequence+is a pair of indices @i < j@ such that @ci*cj < 0@,+and either @j = i + 1@ or @ck = 0@ for all @@ such that @i < k < j@.+-}+signVariations :: (Fractional a, Ord a) => [a] -> Int+signVariations xs =+ length (filter (< 0) pairsMultiplied)+ where+ -- we simply remove zero elements to implement the clause+ -- "ck = 0 for all k such that i < k < j"+ zeroesRemoved = filter (/= 0) xs+ pairsMultiplied = zipWith (*) zeroesRemoved (drop 1 zeroesRemoved)++{-|+Construct the [Sturm sequence+](https://en.wikipedia.org/wiki/Sturm%27s_theorem)+of a given polynomial @p@. The Sturm sequence is given by the polynomials++> p0 = p+> p1 = differentiate p+> p{i+1} = - rem(p{i-1}, pi)++where @rem@ denotes the remainder under 'euclidianDivision'.+We truncate the list when one of the @pi = 0@.++For ease of implementation, we++* construct the 'reverse' of the Sturm sequence.+ This does not affect the number of sign variations that the usage site+ will be interested in.++* assume that the @degree p >= 1@.+-}+reversedSturmSequence :: (Fractional a, Ord a) => Poly a -> [Poly a]+reversedSturmSequence p =+ go [differentiate p, p]+ where+ -- Note that this is called with a list of length 2 and grows the list,+ -- so we don't need to match all cases.+ go ps@(pI : pIminusOne : _)+ | remainder == zero = ps+ | otherwise = go (negate remainder : ps)+ where+ remainder = snd $ euclidianDivision pIminusOne pI+ go _ = error "reversedSturmSequence: impossible"++-- | Check whether a polynomial is monotonically increasing on+-- a given interval.+isMonotonicallyIncreasingOn+ :: (Fractional a, Eq a, Ord a) => Poly a -> (a,a) -> Bool+isMonotonicallyIncreasingOn p (x1,x2) =+ eval p x1 <= eval p x2+ && countRoots (x1, x2, differentiate p) == 0++{-|+Measure whether or not a polynomial is consistently above or below zero,+or equals zero.++Need to consider special cases where there is a root at a boundary point.+-}+compareToZero :: (Fractional a, Eq a, Ord a) => (a, a, Poly a) -> Maybe Ordering+compareToZero (l, u, p)+ | l >= u = error "Invalid interval"+ | p == zero = Just EQ+ | lower * upper < 0 = Nothing -- quick test to eliminate simple cases+ | countRoots (l, u, p) > 0 = Nothing -- polynomial crosses zero+ -- since the polynomial has no roots, the comparison is detmined by the boundary values+ | lower == 0 = Just (compare upper lower)+ | upper == 0 = Just (compare lower upper)+ | lower > 0 = Just GT -- upper must also be > 0 due to the lack of roots+ | otherwise = Just LT -- upper and lower both < 0 due to the lack of roots+ where+ lower = eval p l+ upper = eval p u++{-|+Find the root of a polynomial in a given interval,+assuming that there is exactly one root in the given interval.+This precondition has to be checked through other means,+e.g. 'countRoots'.++We find the root by repeatedly halving the interval in which the root must lie+until its width is less than the specified precision.+Constant and linear polynomials, @degree p <= 1@, are treated as special cases.+-}+findRoot+ :: (Fractional a, Eq a, Num a, Ord a) => a -> (a, a) -> Poly a -> Maybe a+findRoot precision (l, u) p+ -- if the polynomial is zero, the whole interval is a root, so return the basepoint+ | degp < 0 = Just l+ -- if the poly is a non-zero constant, no root is present+ | degp == 0 = Nothing+ -- if the polynomial has degree 1, can calculate the root exactly+ | degp == 1 = Just (-(head ps / last ps)) -- p0 + p1x = 0 => x = -p0/p1+ | precision <= 0 = error "Invalid precision value"+ | otherwise = Just (halveInterval precision l u pl pu)+ where+ Poly ps = p+ degp = degree p+ pu = eval p u+ pl = eval p l+ halveInterval eps x y px py+ -- when the interval is small enough, stop:+ -- the root is in this interval, so take the mid point+ | width <= eps = mid+ -- choose the lower half,+ -- as the polynomial has different signs at the ends+ | px * pmid < 0 = halveInterval eps x mid px pmid+ -- choose the upper half+ | otherwise = halveInterval eps mid y pmid py+ where+ width = y - x+ mid = x + width / 2+ pmid = eval p mid++{-| Otherwise we have a polynomial:+subtract the value we are looking for so that we seek a zero crossing+-}+root+ :: (Ord a, Num a, Eq a, Fractional a)+ => a+ -> a+ -> (a, a)+ -> Poly a+ -> Maybe a+root e x (l, u) p = findRoot e (l, u) (p - constant x)
+ src/Numeric/Probability/Moments.hs view
@@ -0,0 +1,91 @@+{-# LANGUAGE NamedFieldPuns #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+Description : Moments of probability distributions.+-}+module Numeric.Probability.Moments+ ( Moments (..)+ , fromExpectedPowers+ ) where++{-----------------------------------------------------------------------------+ Test+------------------------------------------------------------------------------}++-- | The first four commonly used moments of a probability distribution.+data Moments a = Moments+ { mean :: a+ -- ^ [Mean or Expected Value](https://en.wikipedia.org/wiki/Expected_value)+ -- \( \mu \).+ -- Defined as \( \mu = E[X] \).+ , variance :: a+ -- ^ [Variance](https://en.wikipedia.org/wiki/Variance) \( \sigma^2 \).+ -- Defined as \( \sigma^2 = E[(X - \mu)^2] \).+ -- Equal to \( \sigma^2 = E[X^2] - \mu^2 \).+ , skewness :: a+ -- ^ [Skewness](https://en.wikipedia.org/wiki/Skewness) \( \gamma_1 \).+ -- Defined as+ -- \( \gamma_1 = E\left[\left(\frac{(X - \mu)}{\sigma}\right)^3 \right] \).+ , kurtosis :: a+ -- ^ [Kurtosis](https://en.wikipedia.org/wiki/Kurtosis) \( \kappa \).+ -- Defined as+ -- \( \kappa = E\left[\left(\frac{(X - \mu)}{\sigma}\right)^4 \right] \).+ --+ -- The kurtosis is bounded below: \( \kappa \geq \gamma_1^2 + 1 \).+ }+ deriving (Eq, Show)++-- | Compute the 'Moments' of a probability distribution given+-- the expectation values of the first four powers \( m_k = E[X^k] \).+--+-- > fromExpectedPowers (m1,m2,m3,m4)+fromExpectedPowers+ :: (Ord a, Num a, Fractional a)+ => (a, a, a, a) -> Moments a+fromExpectedPowers (mean, m2, m3, m4)+ | variance == 0 =+ Moments{mean, variance, skewness = 0, kurtosis = 1}+ | otherwise =+ Moments{mean, variance, skewness, kurtosis}+ where+ meanSq = mean * mean++ variance = m2 - meanSq+ sigma = squareRoot variance++ skewness =+ (m3 - 3 * mean * variance - mean * meanSq+ ) / (sigma * variance)++ kurtosis =+ (m4+ - 4 * mean * skewness * sigma * variance+ - 6 * meanSq * variance+ - meanSq * meanSq+ ) / (variance * variance)++-- | Helper function to approximate the square root.+-- Precision: 1e-4 of the given value.+--+-- Uses Heron's iterative method.+squareRoot :: (Ord a, Num a, Fractional a) => a -> a+squareRoot x+ | x < 0 = error "Negative square root input"+ | x == 0 = 0+ | otherwise = goRoot x0+ where+ precision = x / 10000+ x0 = x/2 -- initial guess+ goRoot xi+ | abs (x - xi * xi) <= precision = xi+ | otherwise = goRoot ((xi + x / xi)/2)++{-sqRoot :: a -> a+sqRoot x = + let+ y :: Double+ y = toRational x+ in fromRational . toRational . sqrt y+-}
+ test/Numeric/Function/PiecewiseSpec.hs view
@@ -0,0 +1,290 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# OPTIONS_GHC -Wno-orphans #-}+{-# OPTIONS_GHC -Wno-missing-methods #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+-}+module Numeric.Function.PiecewiseSpec+ ( spec+ , genInterval+ , genPiecewise+ ) where++import Prelude++import Data.Function.Class+ ( eval+ )+import Numeric.Function.Piecewise+ ( Piecewise+ , fromAscPieces+ , fromInterval+ , intervals+ , toAscPieces+ , translateWith+ , trim+ , zipPointwise+ )+import Test.Hspec+ ( Spec+ , describe+ , it+ )+import Test.QuickCheck+ ( Arbitrary+ , Gen+ , Positive (..)+ , (===)+ , (.&&.)+ , arbitrary+ , frequency+ , listOf+ , property+ )++import qualified Data.Function.Class as Fun++{-----------------------------------------------------------------------------+ Tests+------------------------------------------------------------------------------}+spec :: Spec+spec = do+ describe "Test consistency" $ do+ describe "Linear" $ do+ it "eval . translate" $ property $+ \p y x ->+ evalLinear (translateLinear y p) x+ === evalLinear p (x - y)++ describe "Interval" $ do+ it "member intersect" $ property $+ \x y z ->+ member z (intersect x y) === (member z x && member z y)++ describe "fromInterval" $ do+ it "intervals" $ property $+ \(x :: Rational) (Positive d) (o :: Constant) ->+ let y = x + d+ in intervals (fromInterval (x,y) o) === [(x,y)]++ it "eval" $ property $+ \(x :: Rational) (Positive d) (o :: Linear) z ->+ let y = x + d+ p = fromInterval (x, y) o+ in + eval p z+ === (if x <= z && z < y then eval o z else 0)++ describe "mergeBy" $ do+ it "(p + negate p) trims to 0" $ property $+ \(p :: Piecewise Linear) ->+ let z = trim (p + negate p)+ in null (toAscPieces z) === True+ .&&. eval z 0 === 0++ describe "translateWith" $ do+ it "eval . translate" $ property $+ \(p :: Piecewise Linear) x y ->+ eval (translateWith translateLinear y p) x+ === eval p (x - y)++ describe "zipPointwise" $ do+ it "intersects intervals" $ property $+ \p (q :: Piecewise Constant) ->+ allIntervals (zipPointwise (+) p q)+ === [ i+ | ip <- allIntervals p+ , iq <- allIntervals q+ , let i = intersect ip iq+ , i /= Empty+ ]++ it "eval, +" $ property $+ \p (q :: Piecewise Linear) x ->+ eval (zipPointwise (+) p q) x+ === (eval p x + eval q x)++ it "eval, *" $ property $+ \p (q :: Piecewise Constant) x ->+ eval (zipPointwise (*) p q) x+ === (eval p x * eval q x)++ describe "instance Num (Piecewise Q Constant)" $ do+ it "(+)" $ property $+ \p (q :: Piecewise Constant) x ->+ eval (p + q) x+ === (eval p x + eval q x)++ it "(*)" $ property $+ \p (q :: Piecewise Constant) x ->+ eval (p * q) x+ === (eval p x * eval q x)++ it "negate" $ property $+ \(p :: Piecewise Constant) x ->+ eval (negate p) x+ === negate (eval p x)++ it "abs" $ property $+ \(p :: Piecewise Constant) x ->+ eval (abs p) x+ === abs (eval p x)++ it "signum" $ property $+ \(p :: Piecewise Constant) x ->+ eval (signum p) x+ === signum (eval p x)++{-----------------------------------------------------------------------------+ Helper types+ Constant and linear functions+------------------------------------------------------------------------------}+type Q = Rational++-- | Constant function+newtype Constant = Constant Q+ deriving (Eq, Show)++instance Num Constant where+ Constant a1 + Constant a2 = Constant (a1 + a2)+ Constant a1 * Constant a2 = Constant (a1 * a2)+ negate (Constant a) = Constant (negate a)+ abs (Constant a) = Constant (abs a)+ signum (Constant a) = Constant (signum a)+ fromInteger n = Constant (fromInteger n)++instance Fun.Function Constant where+ type instance Domain Constant = Q+ type instance Codomain Constant = Q+ eval (Constant a) _ = a++-- | Linear function with a constant and a slope+data Linear = Linear Q Q+ deriving (Eq, Show)++instance Num Linear where+ Linear a1 b1 + Linear a2 b2 = Linear (a1 + a2) (b1 + b2)+ negate (Linear a b) = Linear (negate a) (negate b)+ fromInteger n = Linear 0 (fromInteger n)++instance Fun.Function Linear where+ type instance Domain Linear = Q+ type instance Codomain Linear = Q+ eval = evalLinear++translateLinear :: Q -> Linear -> Linear+translateLinear y (Linear a b) = Linear a (b - a*y)++evalLinear :: Linear -> Q -> Q+evalLinear (Linear a b) x = a*x + b++{-----------------------------------------------------------------------------+ Helper types+ Intervals+------------------------------------------------------------------------------}+-- | Interval on the real number line.+-- This type does not represent all interval types,+-- only those that are relevant to our purposes here.+data Interval+ = All+ | Empty+ | Before Q -- exclusive+ | After Q -- inclusive+ | FromTo Q Q+ deriving (Eq, Show)++-- | Definition of membership.+member :: Q -> Interval -> Bool+member _ All = True+member _ Empty = False+member z (Before y) = z < y+member z (After x) = x <= z+member z (FromTo x y) = x <= z && z < y++-- | The intersection of two 'Interval' is again an 'Interval'.+intersect :: Interval -> Interval -> Interval+intersect All x = x+intersect x All = x+intersect Empty _ = Empty+intersect _ Empty = Empty+intersect (Before y1) (Before y2) = Before (min y1 y2)+intersect (Before y1) (After x2) = mkFromTo x2 y1+intersect (Before y1) (FromTo x2 y2) = mkFromTo x2 (min y1 y2)+intersect (After x1) (After x2) = After (max x1 x2)+intersect (After x1) (Before y2) = mkFromTo x1 y2+intersect (After x1) (FromTo x2 y2) = mkFromTo (max x1 x2) y2+intersect (FromTo x1 y1) (Before y2) = mkFromTo x1 (min y1 y2)+intersect (FromTo x1 y1) (After x2) = mkFromTo (max x1 x2) y1+intersect (FromTo x1 y1) (FromTo x2 y2) = mkFromTo (max x1 x2) (min y1 y2)++-- | Smart constructor,+-- returns 'Empty' if the endpoint does not come after the starting point.+mkFromTo :: Q -> Q -> Interval+mkFromTo x y = if x < y then FromTo x y else Empty++-- | Return all intervals, +allIntervals :: Fun.Domain o ~ Q => Piecewise o -> [Interval]+allIntervals pieces+ | null xs = [All]+ | otherwise = [Before xmin] <> map (uncurry FromTo) is <> [After xmax]+ where+ xs = map fst (toAscPieces pieces)+ is = zip xs (drop 1 xs)+ xmin = minimum xs+ xmax = maximum xs++{-----------------------------------------------------------------------------+ Random generators+------------------------------------------------------------------------------}+instance Arbitrary Constant where+ arbitrary = Constant <$> arbitrary++instance Arbitrary Linear where+ arbitrary = Linear <$> arbitrary <*> arbitrary++genInterval :: Gen (Q,Q)+genInterval = do+ x <- arbitrary+ Positive d <- arbitrary+ pure (x, x + d)++genFromTo :: Gen Interval+genFromTo = uncurry FromTo <$> genInterval++instance Arbitrary Interval where+ arbitrary = frequency+ [ (1, pure All)+ , (1, pure Empty)+ , (3, Before <$> arbitrary)+ , (3, After <$> arbitrary)+ , (20, genFromTo)+ ]++-- | A list of disjoint and sorted elements.+newtype DisjointSorted a = DisjointSorted [a]+ deriving (Eq, Show)++genDisjointSorted :: Gen (DisjointSorted Rational)+genDisjointSorted =+ DisjointSorted . drop 1 . scanl (\s (Positive d) -> s + d) 0+ <$> listOf arbitrary++instance Arbitrary (DisjointSorted Rational) where+ arbitrary = genDisjointSorted++genPiecewise :: Fun.Domain o ~ Rational => Gen o -> Gen (Piecewise o)+genPiecewise gen = do+ DisjointSorted xs <- genDisjointSorted+ os <- mapM (const gen) xs+ pure $ fromAscPieces $ zip xs os++instance+ (Fun.Domain o ~ Rational, Arbitrary o)+ => Arbitrary (Piecewise o)+ where+ arbitrary = genPiecewise arbitrary
+ test/Numeric/Measure/DiscreteSpec.hs view
@@ -0,0 +1,140 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# OPTIONS_GHC -Wno-orphans #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+-}+module Numeric.Measure.DiscreteSpec+ ( spec+ ) where++import Prelude++import Data.Function.Class+ ( eval+ )+import Numeric.Measure.Discrete+ ( Discrete+ , add+ , convolve+ , dirac+ , distribution+ , fromMap+ , integrate+ , scale+ , toMap+ , total+ , translate+ , zero+ )+import Test.Hspec+ ( Spec+ , describe+ , it+ )+import Test.QuickCheck+ ( Arbitrary+ , Positive (..)+ , (===)+ , (==>)+ , arbitrary+ , cover+ , property+ )++import qualified Data.Map.Strict as Map++{-----------------------------------------------------------------------------+ Tests+------------------------------------------------------------------------------}+spec :: Spec+spec = do+ describe "instance Eq" $ do+ it "add m (scale (-1) m) == zero" $ property $+ \(m :: Discrete Rational) ->+ cover 80 (total m /= 0) "nontrivial"+ $ add m (scale (-1) m) === zero++ it "dirac x /= dirac y" $ property $+ \(x :: Rational) (y :: Rational) ->+ x /= y ==> dirac x /= dirac y++ describe "distribution" $ do+ it "eval and total" $ property $+ \(m :: Discrete Rational) ->+ let xlast = maybe 0 fst $ Map.lookupMax $ toMap m+ in total m+ === eval (distribution m) xlast++ it "eval and scale" $ property $+ \(m :: Discrete Rational) x s->+ eval (distribution (scale s m)) x+ === s * eval (distribution m) x++ describe "integrate" $ do+ it "total" $ property $+ \(m :: Discrete Rational) ->+ integrate (const 1) m+ === total m++ it "linearity, function (+)" $ property $+ \(mx :: Discrete Rational) ->+ let f = id+ in integrate (\x -> f x + f x) mx+ === integrate f mx + integrate f mx ++ it "linearity, measure add" $ property $+ \(mx :: Discrete Rational) my ->+ let f = id+ in integrate f (add mx my)+ === integrate f mx + integrate f my ++ it "linearity, measure scale" $ property $+ \(mx :: Discrete Rational) a ->+ let f = id+ in integrate f (scale a mx)+ === a * integrate f mx++ describe "translate" $ do+ it "distribution" $ property $+ \(m :: Discrete Rational) y x ->+ eval (distribution (translate y m)) x+ === eval (distribution m) (x - y)++ describe "convolve" $ do+ it "dirac" $ property $+ \(x :: Rational) y ->+ convolve (dirac x) (dirac y)+ === dirac (x + y)++ it "total" $ property $+ \mx (my :: Discrete Rational) ->+ total (convolve mx my)+ === total mx * total my++ it "symmetric" $ property $+ \mx (my :: Discrete Rational) ->+ convolve mx my+ === convolve my mx++ it "distributive, left" $ property $+ \mx my (mz :: Discrete Rational) ->+ convolve (add mx my) mz+ === add (convolve mx mz) (convolve my mz) ++ it "distributive, right" $ property $+ \mx my (mz :: Discrete Rational) ->+ convolve mx (add my mz)+ === add (convolve mx my) (convolve mx mz) ++ it "translate, left" $ property $+ \mx (my :: Discrete Rational) (Positive z) ->+ translate z (convolve mx my)+ === convolve (translate z mx) my++{-----------------------------------------------------------------------------+ Random generators+------------------------------------------------------------------------------}+instance (Ord a, Num a, Arbitrary a) => Arbitrary (Discrete a) where+ arbitrary = fromMap . Map.fromList <$> arbitrary
+ test/Numeric/Measure/Finite/MixedSpec.hs view
@@ -0,0 +1,216 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# OPTIONS_GHC -Wno-orphans #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+-}+module Numeric.Measure.Finite.MixedSpec+ ( spec+ ) where++import Prelude++import Data.Function.Class+ ( eval+ )+import Data.Maybe+ ( fromJust+ )+import Numeric.Measure.Finite.Mixed+ ( Measure+ , add+ , convolve+ , dirac+ , distribution+ , fromDistribution+ , integrate+ , isPositive+ , scale+ , support+ , total+ , translate+ , uniform+ , zero+ )+import Numeric.Function.PiecewiseSpec+ ( genPiecewise+ )+import Numeric.Polynomial.SimpleSpec+ ( genPoly+ )+import Test.Hspec+ ( Spec+ , describe+ , it+ )+import Test.QuickCheck+ ( Arbitrary+ , Gen+ , Positive (..)+ , (===)+ , (==>)+ , arbitrary+ , conjoin+ , counterexample+ , cover+ , mapSize+ , once+ , property+ )++import qualified Numeric.Function.Piecewise as Piecewise+import qualified Numeric.Polynomial.Simple as Poly++{-----------------------------------------------------------------------------+ Tests+------------------------------------------------------------------------------}+spec :: Spec+spec = do+ describe "dirac" $ do+ it "total" $ property $+ \(x :: Rational) ->+ total (dirac x) === 1++ describe "uniform" $ do+ it "total" $ property $+ \(x :: Rational) y ->+ total (uniform x y) === 1++ it "support" $ property $+ \(x :: Rational) y ->+ support (uniform x y) === Just (min x y, max x y)++ it "distribution at midpoint" $ property $+ \(x :: Rational) (y :: Rational) ->+ x /= y ==>+ eval (distribution (uniform x y)) ((x + y) / 2) === 1/2++ describe "instance Eq" $ do+ it "add m (scale (-1) m) == zero" $ property $+ \(m :: Measure Rational) ->+ cover 80 (total m /= 0) "nontrivial"+ $ add m (scale (-1) m) === zero+ + it "dirac x /= dirac y" $ property $+ \(x :: Rational) (y :: Rational) ->+ x /= y ==> dirac x /= dirac y++ describe "add" $ do+ it "total" $ property $+ \(mx :: Measure Rational) my ->+ total (add mx my) === total mx + total my++ describe "translate" $ do+ it "distribution" $ property $+ \(m :: Measure Rational) y x ->+ eval (distribution (translate y m)) x+ === eval (distribution m) (x - y)++ describe "convolve" $ do+ it "dirac dirac" $ property $+ \(x :: Rational) y ->+ convolve (dirac x) (dirac y)+ === dirac (x + y)++ it "total" $ property $ mapSize (`div` 10) $+ \mx (my :: Measure Rational) ->+ total (convolve mx my)+ === total mx * total my++ it "dirac translate, left" $ property $ mapSize (`div` 10) $+ \(mx :: Measure Rational) (y :: Rational) ->+ convolve mx (dirac y)+ === translate y mx++ it "dirac translate, right" $ property $ mapSize (`div` 10) $+ \(x :: Rational) (my :: Measure Rational) ->+ convolve (dirac x) my+ === translate x my++ it "symmetric" $ property $ mapSize (`div` 10) $+ \mx (my :: Measure Rational) ->+ convolve mx my+ === convolve my mx++ it "distributive, left" $ property $ mapSize (`div` 12) $+ \mx my (mz :: Measure Rational) ->+ convolve (add mx my) mz+ === add (convolve mx mz) (convolve my mz) ++ it "distributive, right" $ property $ mapSize (`div` 12) $+ \mx my (mz :: Measure Rational) ->+ convolve mx (add my mz)+ === add (convolve mx my) (convolve mx mz) ++ it "translate, left" $ property $ mapSize (`div` 10) $+ \mx (my :: Measure Rational) (Positive z) ->+ translate z (convolve mx my)+ === convolve (translate z mx) my++ describe "isPositive" $ do+ it "scale dirac" $ property $+ \(x :: Rational) w ->+ isPositive (scale w (dirac x))+ === (w >= 0)++ it "sum of positive dirac" $ property $+ \(ws :: [Positive Rational]) ->+ let mkDirac i (Positive w) = scale w (dirac i)+ diracs = zipWith mkDirac [1..] ws+ in isPositive (foldr add zero diracs)+ === True++ it "nfold convolution of uniform" $ once $+ let convolutions :: [Measure Rational]+ convolutions =+ iterate (convolve (uniform 0 1)) (dirac 0)+ prop_isPositive m =+ counterexample (show m)+ $ isPositive m === True+ in conjoin+ $ take 20+ $ map prop_isPositive convolutions++ describe "integrate" $ do+ it "total" $ mapSize (`div` 10) $ property $+ \(m :: Measure Rational) ->+ integrate (Poly.constant 1) m+ === total m++ it "linearity, function (+)" $ mapSize (`div` 10) $ property $+ \f g (mx :: Measure Rational) ->+ integrate (f + g) mx+ === integrate f mx + integrate g mx ++ it "linearity, measure add" $ mapSize (`div` 10) $ property $+ \(mx :: Measure Rational) my ->+ let f = Poly.fromCoefficients [0,1]+ in integrate f (add mx my)+ === integrate f mx + integrate f my ++ it "linearity, measure scale" $ mapSize (`div` 10) $ property $+ \(mx :: Measure Rational) a ->+ let f = Poly.fromCoefficients [0,1]+ in integrate f (scale a mx)+ === a * integrate f mx++{-----------------------------------------------------------------------------+ Random generators+------------------------------------------------------------------------------}+genMeasure :: Gen (Measure Rational)+genMeasure =+ fromJust . fromDistribution . setLastPieceConstant <$> genPiecewise genPoly+ where+ setLastPieceConstant =+ Piecewise.fromAscPieces+ . setLastPieceConstant'+ . Piecewise.toAscPieces++ setLastPieceConstant' [] = []+ setLastPieceConstant' [(x, o)] = [(x, Poly.constant (eval o x))]+ setLastPieceConstant' (p : ps) = p : setLastPieceConstant' ps++instance Arbitrary (Measure Rational) where+ arbitrary = genMeasure
+ test/Numeric/Measure/ProbabilitySpec.hs view
@@ -0,0 +1,252 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# OPTIONS_GHC -Wno-orphans #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+-}+module Numeric.Measure.ProbabilitySpec+ ( spec+ ) where++import Prelude++import Data.Function.Class+ ( eval+ )+import Data.Ratio+ ( (%)+ )+import Numeric.Polynomial.SimpleSpec+ ( genPositivePoly+ )+import Numeric.Measure.Probability+ ( Prob+ , choice+ , convolve+ , dirac+ , distribution+ , expectation+ , fromDistribution+ , fromMeasure+ , unsafeFromMeasure+ , measure+ , moments+ , support+ , translate+ , uniform+ )+import Numeric.Probability.Moments+ ( Moments (..)+ )+import Test.Hspec+ ( Spec+ , describe+ , it+ )+import Test.QuickCheck+ ( Arbitrary+ , Gen+ , NonNegative (..)+ , Positive (..)+ , (===)+ , (==>)+ , arbitrary+ , choose+ , chooseInteger+ , frequency+ , getSize+ , mapSize+ , oneof+ , property+ , scale+ , vectorOf+ )++import qualified Numeric.Measure.Finite.Mixed as M+import qualified Numeric.Polynomial.Simple as Poly++{-----------------------------------------------------------------------------+ Tests+------------------------------------------------------------------------------}+spec :: Spec+spec = do+ describe "uniform" $ do+ it "support" $ property $+ \(x :: Rational) y ->+ support (uniform x y) === Just (min x y, max x y)++ it "distribution at midpoint" $ property $+ \(x :: Rational) (y :: Rational) ->+ x /= y ==>+ eval (distribution (uniform x y)) ((x + y) / 2) === 1/2++ describe "instance Eq" $ do + it "dirac x /= dirac y" $ property $+ \(x :: Rational) (y :: Rational) ->+ x /= y ==> dirac x /= dirac y++ describe "elimination . introduction" $ do + it "unsafe fromMeasure . measure" $ property $+ \(m :: Prob Rational) ->+ m === (unsafeFromMeasure . measure) m++ it "fromMeasure . measure" $ property $+ \(m :: Prob Rational) ->+ Just m === (fromMeasure . measure) m++ it "unsafe fromDistribution . distribution" $ property $+ \(m :: Prob Rational) ->+ Just m ===+ (fmap unsafeFromMeasure . M.fromDistribution . distribution) m++ it "fromDistribution . distribution" $ property $+ \(m :: Prob Rational) ->+ Just m ===+ (fromDistribution . distribution) m++ describe "expectation" $ do+ it "unit" $ property $+ \(m :: Prob Rational) ->+ expectation (Poly.constant 1) m+ === 1++ it "positivity" $ mapSize (`div` 2) $ property $+ \(m :: Prob Rational) (PositivePoly p) ->+ expectation p m+ >= 0++ describe "moments" $ do+ it "mean is additive" $ mapSize (`div` 10) $ property $+ \(mx :: Prob Rational) my ->+ let mean' = mean . moments+ in mean' (convolve mx my)+ === mean' mx + mean' my++ it "variance is additive" $ mapSize (`div` 10) $ property $+ \(mx :: Prob Rational) my ->+ let variance' = variance . moments+ in variance' (convolve mx my)+ === variance' mx + variance' my++ it "skewness absorbs translate" $ property $+ \(m :: Prob Rational) y ->+ let skewness' = skewness . moments+ in skewness' (translate y m)+ === skewness' m++ it "kurtosis absorbs translate" $ property $+ \(m :: Prob Rational) y ->+ let kurtosis' = kurtosis . moments+ in kurtosis' (translate y m)+ === kurtosis' m++ it "kurtosis bounded below" $ property $+ \(m :: Prob Rational) ->+ let ms = moments m+ in kurtosis ms+ >= (skewness ms)^(2 :: Int) + 1++ describe "choice" $ do+ it "distribution" $ property $+ \(Probability p) (mx :: Prob Rational) my z ->+ eval (distribution (choice p mx my)) z+ === p * eval (distribution mx) z+ + (1-p) * eval (distribution my) z++ describe "translate" $ do+ it "distribution" $ property $+ \(m :: Prob Rational) y x ->+ eval (distribution (translate y m)) x+ === eval (distribution m) (x - y)++ describe "convolve" $ do+ it "dirac dirac" $ property $+ \(x :: Rational) y ->+ convolve (dirac x) (dirac y)+ === dirac (x + y)++ it "dirac translate, left" $ property $ mapSize (`div` 10) $+ \(mx :: Prob Rational) (y :: Rational) ->+ convolve mx (dirac y)+ === translate y mx++ it "dirac translate, right" $ property $ mapSize (`div` 10) $+ \(x :: Rational) (my :: Prob Rational) ->+ convolve (dirac x) my+ === translate x my++ it "symmetric" $ property $ mapSize (`div` 10) $+ \mx (my :: Prob Rational) ->+ convolve mx my+ === convolve my mx++ it "translate, left" $ property $ mapSize (`div` 10) $+ \mx (my :: Prob Rational) (Positive z) ->+ translate z (convolve mx my)+ === convolve (translate z mx) my++{-----------------------------------------------------------------------------+ Random generators+------------------------------------------------------------------------------}+newtype PositivePoly = PositivePoly (Poly.Poly Rational)+ deriving (Eq, Show)++instance Arbitrary PositivePoly where+ arbitrary = PositivePoly <$> genPositivePoly++newtype Probability = Probability Rational+ deriving (Eq, Show)++instance Arbitrary Probability where+ arbitrary = Probability <$> genProbability++instance Arbitrary (Prob Rational) where+ arbitrary = scale (`div` 15) genProb++-- | Generate a random 'Prob' by generating a random expression.+genProb :: Gen (Prob Rational)+genProb = do+ size <- getSize+ genProbFromList =<< vectorOf size genSimpleProb++-- | Generate a 'uniform'.+genUniform :: Gen (Prob Rational)+genUniform = do+ NonNegative a <- arbitrary+ Positive d <- arbitrary+ pure $ uniform a (a + d)++-- | Generate a 'dirac'.+genDirac :: Gen (Prob Rational)+genDirac = do+ NonNegative a <- arbitrary+ pure $ dirac a++-- | Generate a simple probability measure — one of 'uniform', 'dirac'.+genSimpleProb :: Gen (Prob Rational)+genSimpleProb =+ frequency [(20, genUniform), (4, genDirac)]++-- | Generate a random probability in the interval (0,1).+genProbability :: Gen Rational+genProbability = do+ denominator <- chooseInteger (1,2^(20 :: Int))+ numerator <- chooseInteger (0, denominator)+ pure (numerator % denominator)++-- | Generate a random 'Prob' by combining a given list+-- of 'Prob' with random operations.+genProbFromList :: [Prob Rational] -> Gen (Prob Rational)+genProbFromList [] = pure $ dirac 0+genProbFromList [x] = pure x+genProbFromList xs = do+ n <- choose (1, length xs - 1)+ let (ys, zs) = splitAt n xs+ genOp <*> genProbFromList ys <*> genProbFromList zs+ where+ genChoice = do+ p <- genProbability+ pure $ choice p+ genOp = oneof [pure convolve, genChoice]
+ test/Numeric/Polynomial/SimpleSpec.hs view
@@ -0,0 +1,401 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# OPTIONS_GHC -Wno-orphans #-}+{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}++{-|+Copyright : Predictable Network Solutions Ltd., 2020-2024+License : BSD-3-Clause+-}+module Numeric.Polynomial.SimpleSpec+ ( spec+ , genPoly+ , genPositivePoly+ ) where++import Prelude++import Data.List+ ( nub+ )+import Data.Traversable+ ( for+ )+import Numeric.Polynomial.Simple+ ( Poly+ , compareToZero+ , constant+ , convolve+ , countRoots+ , degree+ , differentiate+ , display+ , euclidianDivision+ , eval+ , fromCoefficients+ , integrate+ , isMonotonicallyIncreasingOn+ , lineFromTo+ , monomial+ , root+ , scale+ , scaleX+ , translate+ , zero+ )+import Test.Hspec+ ( Spec+ , before_+ , describe+ , it+ , pendingWith+ )+import Test.QuickCheck+ ( Arbitrary+ , Gen+ , NonNegative (..)+ , Positive (..)+ , Property+ , (===)+ , (==>)+ , (.&&.)+ , arbitrary+ , counterexample+ , forAll+ , frequency+ , listOf+ , mapSize+ , property+ , withMaxSuccess+ )++import qualified Test.QuickCheck as QC++{-----------------------------------------------------------------------------+ Tests+------------------------------------------------------------------------------}+xit' :: String -> String -> Property -> Spec+xit' reason label = before_ (pendingWith reason) . it label++spec :: Spec+spec = do+ describe "constant" $ do+ it "eval" $ property $+ \c (x :: Rational) ->+ eval (constant c) x === c++ describe "scale" $ do+ it "eval" $ property $+ \a p (x :: Rational) ->+ eval (scale a p) x === a * eval p x++ describe "scaleX" $ do+ it "degree" $ property $+ \(p :: Poly Rational) ->+ (degree p >= 0)+ ==> (degree (scaleX p) === 1 + degree p)++ it "eval" $ property $+ \p (x :: Rational) ->+ eval (scaleX p) x === x * eval p x++ it "zero" $ withMaxSuccess 1 $ property $+ scaleX zero == (zero :: Poly Rational)++ describe "(+)" $ do+ it "eval" $ property $+ \p q (x :: Rational) ->+ eval (p + q) x === eval p x + eval q x++ describe "(*)" $ do+ it "eval" $ property $+ \p q (x :: Rational) ->+ eval (p * q) x === eval p x * eval q x++ describe "display" $ do+ it "step == 0" $ property $+ \(l :: Rational) (Positive d) ->+ let u = l + d+ in display zero (l, u) 0+ === zip [l, u] (repeat 0)++ it "zero" $ property $+ \(l :: Rational) (Positive d) (Positive (n :: Integer)) ->+ let u = l + d+ s = (u - l) / fromIntegral (min 100 n)+ in display zero (l, u) s+ === zip (nub ([l, l+s .. u] <> [u])) (repeat 0)++ describe "lineFromTo" $ do+ it "degree" $ property $+ \x1 (x2 :: Rational) y1 y2 ->+ let p = lineFromTo (x1, y1) (x2, y2)+ in degree p <= 1++ it "eval" $ property $+ \x1 (x2 :: Rational) y1 y2 ->+ let p = lineFromTo (x1, y1) (x2, y2)+ in x1 /= x2+ ==> (eval p x1 === y1 .&&. eval p x2 == y2)+++ describe "integrate" $ do+ it "eval" $ property $+ \(p :: Poly Rational) ->+ eval (integrate p) 0 === 0++ it "integrate . differentiate" $ property $+ \(p :: Poly Rational) ->+ integrate (differentiate p) === p - constant (eval p 0)++ describe "differentiate" $ do+ it "differentiate . integrate" $ property $+ \(p :: Poly Rational) ->+ differentiate (integrate p) === p++ it "Leibniz rule" $ property $+ \(p :: Poly Rational) q ->+ differentiate (p * q)+ === differentiate p * q + p * differentiate q++ describe "translate" $ do+ it "eval" $ property $+ \p y (x :: Rational) ->+ eval (translate y p) x === eval p (x - y)++ it "differentiate" $ property $+ \p (y :: Rational) ->+ differentiate (translate y p)+ === translate y (differentiate p)++ describe "euclidianDivision" $ do+ it "a = q * b + r, and degree r < degree b" $ property $+ \a (b :: Poly Rational) ->+ let (q, r) = euclidianDivision a b in+ b /= zero ==>+ (a === q*b + r .&&. degree r < degree b)++ describe "convolve" $ do+ it "product of integrals" $ property $ mapSize (`div` 6) $+ \(NonNegative x1) (Positive d1) (NonNegative x2) (Positive d2)+ p (q :: Poly Rational) ->+ let p1 = (x1, x1 + d1, p)+ q1 = (x2, x2 + d2, q)+ in+ integrateInterval p1 * integrateInterval q1+ === integratePieces (convolve p1 q1)++ describe "countRoots" $ do+ it "counts distinct roots in open interval" $ property $+ \(PolyWithRealRoots p roots) (x1 :: Rational) (Positive d) ->+ let x2 = x1 + d in+ countRoots (x1, x2, p)+ === countRoots' (x1, x2) roots++ it "handles roots at boundary" $ mapSize (`div` 2) $ property $+ \(PolyWithRealRoots p _) (x1 :: Rational) (Positive d) ->+ let x2 = x1 + d+ xx = monomial 1 1+ rootCount = countRoots (x1, x2, p)+ in countRoots (x1, x2, p * (xx - constant x1))+ === rootCount+ .&&.+ countRoots (x1, x2, p * (xx - constant x2))+ === rootCount++ describe "root" $ do+ it "cubic polynomial" $ property $ mapSize (`div` 5) $+ \(x1 :: Rational) (Positive dx3) ->+ let xx = scaleX (constant 1) :: Poly Rational+ x2 = 0.6 * x1 + 0.4 * x3+ x3 = x1 + dx3+ p = (xx - constant x1) * (xx - constant x2) * (xx - constant x3)+ l = x1 + 100 * epsilon+ u = x3 - 100 * epsilon+ epsilon = (x3-x1)/(1000*1000*50)+ Just x2' = root epsilon 0 (l, u) p+ in+ property $ abs (x2' - x2) <= epsilon++ xit' "bug" "cubic polynomial, midpoint" $ property $ mapSize (`div` 5) $+ \(x1 :: Rational) (Positive dx3) ->+ let xx = scaleX (constant 1) :: Poly Rational+ x2 = (x1 + x3) / 2+ x3 = x1 + dx3+ p = (xx - constant x1) * (xx - constant x2) * (xx - constant x3)+ l = x1 + 100 * epsilon+ u = x3 - 100 * epsilon+ epsilon = (x3-x1)/(1000*1000*50)+ Just x2' = root epsilon 0 (l, u) p+ in+ id+ $ counterexample ("interval = " <> show (l,u))+ $ counterexample ("countRoots = " <> show (countRoots (l, u, p)))+ $ counterexample ("expected root = " <> show x2)+ $ counterexample ("eval polynomial at expected root = " <> show (eval p x2))+ $ counterexample ("epsilon = " <> show epsilon)+ $ counterexample ("found root = " <> show x2')+ $ counterexample ("root within range of other root " <> show (abs (x2' - x3) <= 20*epsilon))+ $ property $ abs (x2' - x2) <= epsilon++ describe "isMonotonicallyIncreasingOn" $+ it "quadratic polynomial" $ property $+ \(x1 :: Rational) (Positive d) ->+ let xx = scaleX (constant 1)+ p = negate ((xx - constant x1) * (xx - constant x2))+ x2 = x1 + d+ xmid = (x1 + x2) / 2+ in+ isMonotonicallyIncreasingOn p (x1,xmid) === True++ describe "compareToZero" $ do+ it "lineFromTo" $ property $+ \(x1 :: Rational) (Positive dx) y1 (Positive dy) ->+ let x2 = x1 + dx+ y2 = y1 + dy+ p = lineFromTo (x1, y1) (x2, y2)+ result+ | y1 == 0 && y2 == 0 = Just EQ+ | y1 >= 0 = Just GT+ | y2 <= 0 = Just LT+ | otherwise = Nothing+ in+ compareToZero (x1, x2, p)+ === result++ it "quadratic polynomial with two roots" $ property $+ \(x1 :: Rational) (Positive d) ->+ let xx = scaleX (constant 1)+ p = (xx - constant x1 + 1) * (xx - constant x2 - 1)+ x2 = x1 + d+ in+ compareToZero (x1, x2, p) === Just LT++ it "quadratic polynomial + a0" $ property $+ \(x1 :: Rational) a0 ->+ let xx = scaleX (constant 1)+ p = (xx - constant x1)^(2 :: Int) + constant a0+ in+ compareToZero (x1 - abs a0 - 1, x1 + abs a0 + 1, p)+ === + if a0 > 0+ then Just GT+ else Nothing++ describe "genPositivePoly" $+ it "eval" $ property $+ \(x :: Rational) ->+ forAll genPositivePoly $ \p ->+ eval p x > 0++ describe "genPolyWithRealRoots" $+ it "eval" $ property $+ \(PolyWithRealRoots (p :: Poly Rational) (Roots roots)) ->+ all (\x -> eval p x == 0) $ map fst roots++{-----------------------------------------------------------------------------+ Helper functions+------------------------------------------------------------------------------}+-- | Definite integral of a polynomial over an interval.+integrateInterval+ :: (Eq a, Num a, Fractional a) => (a, a, Poly a) -> a+integrateInterval (x, y, p) = eval pp y - eval pp x+ where pp = integrate p++-- | Definite integral of a sequence of polynomials over pieces.+integratePieces+ :: (Eq a, Num a, Fractional a) => [(a, Poly a)] -> a+integratePieces = sum . map integrateInterval . intervals+ where+ intervals pieces =+ [ (x, y, p)+ | ((x, p), y) <- zip pieces $ drop 1 $ map fst pieces+ ]++-- | Multiplicity of a root.+type Multiplicity = Int++-- | A list of roots with multiplicity.+newtype Roots a = Roots [(a, Multiplicity)]+ deriving (Eq, Show)++-- | Use [Vieta's theorem+-- ](https://en.wikipedia.org/wiki/Vieta%27s_formulas)+-- to convert a list of roots with mulitiplicities into+-- a polynomial with exactly those roots.+fromRoots :: (Ord a, Num a) => Roots a -> Poly a+fromRoots (Roots xms) =+ product $ map (\(r,m) -> (xx - constant r) ^ m) xms+ where+ xx = monomial 1 1++-- | Count the distinct number of real roots+-- that lie in the given, open interval.+countRoots' :: Ord a => (a, a) -> Roots a -> Int+countRoots' (xl, xr) (Roots xs) =+ length . filter (\x -> xl < x && x < xr) $ map fst xs++{-----------------------------------------------------------------------------+ Random generators+------------------------------------------------------------------------------}+-- | Generate an arbitrary polynomial.+genPoly :: Gen (Poly Rational)+genPoly = fromCoefficients <$> listOf arbitrary++instance Arbitrary (Poly Rational) where+ arbitrary = genPoly++-- | Generate a quadratic polynomial that is positive,+-- i.e. has no real roots and is always larger than zero.+genQuadraticPositivePoly :: Gen (Poly Rational)+genQuadraticPositivePoly = do+ let xx = fromCoefficients [0, 1]+ x0 <- constant <$> arbitrary+ Positive b <- arbitrary+ pure $ (xx - x0) * (xx - x0) + constant b++-- | Generate a positive polynomial, i.e. @eval p x > 0@ for all @x@.+genPositivePoly :: Gen (Poly Rational)+genPositivePoly =+ QC.scale (`div` 3) $ product <$> listOf genQuadraticPositivePoly++-- | A list of disjoint and sorted elements.+newtype DisjointSorted a = DisjointSorted [a]+ deriving (Eq, Show)++genDisjointSorted :: Gen (DisjointSorted Rational)+genDisjointSorted = + DisjointSorted . drop 1 . scanl (\s (Positive d) -> s + d) 0+ <$> listOf arbitrary++instance Arbitrary (DisjointSorted Rational) where+ arbitrary = genDisjointSorted++genMultiplicity :: Gen Multiplicity+genMultiplicity =+ frequency [(20, pure 1), (2, pure 2), (2, pure 3), (1, pure 7)]++genRoots :: Gen (Roots Rational)+genRoots = do+ DisjointSorted xs <- arbitrary+ xms <- for xs $ \x -> do+ multiplicity <- genMultiplicity+ pure $ (x, multiplicity)+ pure $ Roots xms++instance Arbitrary (Roots Rational) where+ arbitrary = genRoots++-- | A polynomial with known real roots.+-- The polynomial may have additional complex roots.+data PolyWithRealRoots a = PolyWithRealRoots (Poly a) (Roots a)+ deriving (Eq, Show)++genPolyWithRealRoots :: Gen (PolyWithRealRoots Rational)+genPolyWithRealRoots = do+ roots <- QC.scale (`div` 7) $ arbitrary+ q <- QC.scale (`div` 11) $ genPositivePoly+ pure $ PolyWithRealRoots (fromRoots roots * q) roots++instance Arbitrary (PolyWithRealRoots Rational) where+ arbitrary = genPolyWithRealRoots
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}