src/Numeric/Probability.hs

   1 -- For `ProbabilityScale`
   2 {-# LANGUAGE DataKinds #-}
   3
   4 module Numeric.Probability (
   5   Probability,
   6   ProbabilityScale,
   7   ProbabilityBounded (..),
   8   probability,
   9   safeProbability,
  10   runProbability,
  11   assertProbability,
  12 )
  13 where
  14
  15 import Control.Monad (Monad (..))
  16 import Data.Eq (Eq)
  17 import Data.Function (id)
  18 import Data.Maybe (Maybe (..), fromJust)
  19 import Data.Monoid (Monoid (..))
  20 import Data.Ord (Ord (..))
  21 import Data.Proxy (Proxy (..))
  22 import Data.Semigroup (Semigroup (..))
  23 import Data.Validity (Validity (..), declare)
  24 import Data.Word (Word64)
  25 import GHC.Generics (Generic)
  26 import GHC.TypeNats (Natural, natVal)
  27 import Logic
  28 import Logic.Theory.Bool (type (&&))
  29 import Logic.Theory.Ord (type (<=))
  30 import Numeric.Decimal (Decimal (..), MonadThrow (..))
  31 import Numeric.Decimal qualified as Decimal
  32 import System.Random (Random)
  33 import Text.Show (Show (show))
  34 import Prelude (Bounded (..), Enum, Fractional (..), Integral, Num (..), Rational, Real, error, (^))
  35
  36 type Probability = Decimal Decimal.RoundHalfEven ProbabilityScale ProbabilityBounded
  37 instance Validity Probability where
  38   validate (Decimal wb) = validate wb
  39
  40 probability :: MonadThrow m => Rational -> m Probability
  41 probability = Decimal.fromRationalDecimalBoundedWithRounding
  42 {-# INLINE probability #-}
  43
  44 safeProbability :: r ::: Rational / r <= 0 && r <= 1 -> Probability
  45 safeProbability (Named r) = fromJust (probability r)
  46
  47 runProbability :: Probability -> Rational
  48 runProbability = Decimal.toRationalDecimal
  49 {-# INLINE runProbability #-}
  50
  51 assertProbability :: Rational -> Probability
  52 assertProbability r = case probability r of
  53   Just p -> p
  54   Nothing -> error ("assertProbability: " <> show r)
  55
  56 instance Num (Decimal.Arith Probability) where
  57   (+) = Decimal.bindM2 Decimal.plusDecimalBounded
  58   (-) = Decimal.bindM2 Decimal.minusDecimalBounded
  59   (*) = Decimal.bindM2 Decimal.timesDecimalBoundedWithRounding
  60   abs = id
  61   signum m = m >>= Decimal.signumDecimalBounded
  62   fromInteger = Decimal.fromIntegerDecimalBoundedIntegral
  63
  64 instance Fractional (Decimal.Arith Probability) where
  65   (/) = Decimal.bindM2 Decimal.divideDecimalBoundedWithRounding
  66   fromRational = probability
  67
  68 -- >>> 10^19 <= (fromIntegral (maxBound :: Word64) :: Integer
  69 -- True
  70 -- >>> 10^20 <= (fromIntegral (maxBound :: Word64) :: Integer
  71 -- False
  72 type ProbabilityScale = 19
  73
  74 newtype ProbabilityBounded = ProbabilityBounded {unProbabilityBounded :: Word64}
  75   deriving (Show, Eq, Ord, Num, Real, Integral, Enum, Random, Generic)
  76 instance Bounded ProbabilityBounded where
  77   minBound = ProbabilityBounded 0
  78   maxBound = ProbabilityBounded (10 ^ (natVal (Proxy @ProbabilityScale)))
  79 instance Validity ProbabilityBounded where
  80   validate (ProbabilityBounded w) =
  81     mconcat
  82       [ declare ("The contained word is smaller or equal to 10 ^ ProbabilityScale = " <> show (10 ^ n :: Natural)) (w <= 10 ^ n)
  83       ]
  84     where
  85       n :: Natural
  86       n = natVal (Proxy @ProbabilityScale)