0 | module Numeric.Log
  1 |
  2 | ||| Numeric utility functions
  3 | export
  4 | isInf : (FromDouble a, Eq a, Fractional a) => a -> Bool
  5 | isInf = (1.0/0.0 ==)
  6 |
  7 | export
  8 | negInf : (Neg a, Fractional a) => a
  9 | negInf = -(1/0)
 10 |
 11 | %foreign "C:c_log1p,log-domain"
 12 | c_log1p : Double -> Double
 13 |
 14 | %foreign "C:c_expm1,log-domain"
 15 | c_expm1 : Double -> Double
 16 |
 17 | log1pexp : Double -> Double
 18 | log1pexp a = c_log1p (exp a)
 19 |
 20 | log1mexp : Double -> Double
 21 | log1mexp a = c_log1p (negate (exp a))
 22 |
 23 | ||| Log Domain
 24 | public export
 25 | record Log (a : Type) where
 26 |   constructor Exp
 27 |   ln : a
 28 |
 29 | public export
 30 | Show (Log Double) where
 31 |   show (Exp a) = show (exp a)
 32 |
 33 | public export
 34 | Functor Log where
 35 |   map f (Exp a) = Exp (f a)
 36 |
 37 | public export
 38 | Num (Log Double) where
 39 |   fromInteger = Exp . log . fromInteger
 40 |   Exp a * Exp b = Exp (a + b)
 41 |   Exp a + Exp b =
 42 |     if a == b && isInf a && isInf b then Exp a else
 43 |     if a >= b then Exp (a + log1pexp (b - a)) else
 44 |     Exp (b + log1pexp (a - b))
 45 |
 46 | public export
 47 | Abs (Log Double) where
 48 |   abs = id
 49 |
 50 | public export
 51 | Neg (Log Double) where
 52 |   negate (Exp a) = if isInf a && a < 0 then Exp negInf else Exp (0/0)
 53 |   Exp a - Exp b = if isInf a && isInf b && a < 0 && b < 0 then Exp negInf else Exp (a + log1mexp (b - a))
 54 |
 55 | public export
 56 | Fractional (Log Double) where
 57 |   Exp a / Exp b = Exp (a - b)
 58 |
 59 | ||| Note: Exp a == Exp b is not equivalent to exp a == exp b, and similarly for Exp a < Exp b.
 60 | ||| Checking that a value 'x' in (Exp x) is not zero in the non-log domain should be done via 'exp x /= 0', rather than as `Exp x /= Exp (log 0).
 61 | public export
 62 | Eq a => Eq (Log a) where
 63 |   Exp a == Exp b = a == b 
 64 |
 65 | public export
 66 | Ord a => Ord (Log a) where
 67 |   Exp a < Exp b = a < b
 68 |
 69 | ||| Explicitly convert a value to the log domain
 70 | public export
 71 | interface ToLogDomain a where
 72 |   toLogDomain   : a -> Log Double
 73 |
 74 | public export
 75 | ToLogDomain Double where
 76 |   toLogDomain = Exp . log
 77 |
 78 | public export
 79 | ToLogDomain Nat where
 80 |   toLogDomain = Exp . log . cast
 81 |
 82 | public export
 83 | ToLogDomain Int where
 84 |   toLogDomain = Exp . log . cast
 85 |
 86 | public export
 87 | ToLogDomain Integer where
 88 |   toLogDomain   = Exp . log . cast
 89 |
 90 | ||| Explicitly convert a double from the log domain
 91 | public export
 92 | fromLogDomain : Log Double -> Double
 93 | fromLogDomain = exp . ln 
 94 |
 95 | ||| Efficiently and accurately compute the sum of a set of log-domain numbers
 96 | data Acc a = MkAcc {-# UNPACK #-} Int64 a | None
 97 |
 98 | export
 99 | sum : (Foldable f) => f (Log Double) -> Log Double
100 | sum xs = Exp $ case foldl step1 None xs of
101 |   None => negInf
102 |   MkAcc nm1 a => 
103 |     if isInf a then a
104 |     else a + c_log1p (foldl (step2 a) 0 xs + cast nm1)
105 |   where
106 |     step1 : Acc Double -> Log Double -> Acc Double
107 |     step1 None      (Exp x) = MkAcc 0 x
108 |     step1 (MkAcc n y) (Exp x) = MkAcc (n + 1) (max x y)
109 |     step2 : Double -> Double -> Log Double -> Double
110 |     step2 a r (Exp x) = r + c_expm1 (x - a)