src/Worksheets/Utils/Probability.hs

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