diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,1 @@ +Probability distributions, represented by piecewise polynomials. diff --git a/benchmark/Main.hs b/benchmark/Main.hs new file mode 100644 --- /dev/null +++ b/benchmark/Main.hs @@ -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 + ] + ] diff --git a/probability-polynomial.cabal b/probability-polynomial.cabal new file mode 100644 --- /dev/null +++ b/probability-polynomial.cabal @@ -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 diff --git a/src/Data/Function/Class.hs b/src/Data/Function/Class.hs new file mode 100644 --- /dev/null +++ b/src/Data/Function/Class.hs @@ -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 diff --git a/src/Numeric/Function/Piecewise.hs b/src/Numeric/Function/Piecewise.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Function/Piecewise.hs @@ -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" diff --git a/src/Numeric/Measure/Discrete.hs b/src/Numeric/Measure/Discrete.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Measure/Discrete.hs @@ -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 + ] diff --git a/src/Numeric/Measure/Finite/Mixed.hs b/src/Numeric/Measure/Finite/Mixed.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Measure/Finite/Mixed.hs @@ -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 diff --git a/src/Numeric/Measure/Probability.hs b/src/Numeric/Measure/Probability.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Measure/Probability.hs @@ -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 diff --git a/src/Numeric/Polynomial/Simple.hs b/src/Numeric/Polynomial/Simple.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Polynomial/Simple.hs @@ -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) diff --git a/src/Numeric/Probability/Moments.hs b/src/Numeric/Probability/Moments.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Probability/Moments.hs @@ -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 +-} diff --git a/test/Numeric/Function/PiecewiseSpec.hs b/test/Numeric/Function/PiecewiseSpec.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/Function/PiecewiseSpec.hs @@ -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 diff --git a/test/Numeric/Measure/DiscreteSpec.hs b/test/Numeric/Measure/DiscreteSpec.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/Measure/DiscreteSpec.hs @@ -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 diff --git a/test/Numeric/Measure/Finite/MixedSpec.hs b/test/Numeric/Measure/Finite/MixedSpec.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/Measure/Finite/MixedSpec.hs @@ -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 diff --git a/test/Numeric/Measure/ProbabilitySpec.hs b/test/Numeric/Measure/ProbabilitySpec.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/Measure/ProbabilitySpec.hs @@ -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] diff --git a/test/Numeric/Polynomial/SimpleSpec.hs b/test/Numeric/Polynomial/SimpleSpec.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/Polynomial/SimpleSpec.hs @@ -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 diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,1 @@ +{-# OPTIONS_GHC -F -pgmF hspec-discover #-}