SciBaseTypes 0.0.0.1 → 0.1.0.0

raw patch · 9 files changed

+525/−343 lines, 9 filesdep +semiringsdep ~basedep ~log-domainPVP ok

version bump matches the API change (PVP)

Dependencies added: semirings

Dependency ranges changed: base, log-domain

API changes (from Hackage documentation)

Files

• Algebra/Structure/SemiRing.hs +0/−194
• Algebra/Structure/Semiring.hs +204/−0
• Numeric/Discretized.hs +115/−24
• Numeric/LogDomain.hs +16/−2
• SciBaseTypes.cabal +26/−33
• StatisticalMechanics/Ensemble.hs +27/−3
• Statistics/Odds.hs +25/−11
• Statistics/Probability.hs +107/−75
• changelog.md +5/−1

@@ -1,194 +0,0 @@---- | A set with two binary operations, one for addition (@srplus@), one for--- multiplication (@srmul@). Together with a neutral element for @srplus@,--- named @srzero@, and one for @srmul@, named @srone@.--module Algebra.Structure.SemiRing where--import Control.DeepSeq (NFData(..))-import Data.Coerce-import Data.Monoid hiding ((<>))-import Data.Semigroup-import Data.Vector.Unboxed.Deriving-import Data.Vector.Unboxed (Unbox)-import GHC.Generics-import Unsafe.Coerce--import Numeric.Limits------ * The 'SemiRing' type class.---- | The semiring operations and neutral elements.--class SemiRing a where-  srplus  ∷ a → a → a-  srmul   ∷ a → a → a-  srzero  ∷ a-  srone   ∷ a---- | Unicode variant of @srplus@.--infixl 6 ⊕-infixl 6 `srplus`-(⊕) ∷ SemiRing a ⇒ a → a → a-(⊕) = srplus-{-# Inline (⊕) #-}---- | Unicode variant of @srmul@.--infixl 7 ⊗-infixl 7 `srmul`-(⊗) ∷ SemiRing a ⇒ a → a → a-(⊗) = srmul-{-# Inline (⊗) #-}------ * Newtype wrappers for 'SemiRing' that make the semiring to use explicit.--- This is important, because several types, say Prob(ability) have multiple--- useful semiring instances.------ 'Data.Monoid' in @base@ provides a number of newtype wrappers (@Sum@,--- @Product@, etc) for monoids, which have one binary operation and identity.--- There is, obviously, overlap with the structures constructed here.---- | The Viterbi SemiRing. It maximizes over the product.--newtype Viterbi x = Viterbi { getViterbi ∷ x }-  deriving (Eq, Ord, Read, Show, Bounded, Generic, Generic1, Num)--derivingUnbox "Viterbi"-  [t| forall x . Unbox x ⇒ Viterbi x → x |]  [| getViterbi |]  [| Viterbi |]--instance NFData x ⇒ NFData (Viterbi x) where-  rnf (Viterbi x) = rnf x-  {-# Inline rnf #-}---- |------ TODO Shall we have generic instances, or specific ones like @SemiRing--- (Viterbi Prob)@?------ TODO Consider either a constraint @ProbLike x@ or the above.--instance (Ord x, Num x) ⇒ SemiRing (Viterbi x) where-  srplus (Viterbi x) (Viterbi y) = Viterbi $ max x y-  srmul  (Viterbi x) (Viterbi y) = Viterbi $ x * y-  srzero = Viterbi 0-  srone  = Viterbi 1-  {-# Inline srplus #-}-  {-# Inline srmul  #-}-  {-# Inline srzero #-}-  {-# Inline srone  #-}---- | The tropical MinPlus SemiRing. It minimizes over the sum.--newtype MinPlus x = MinPlus { getMinPlus ∷ x }-  deriving (Eq, Ord, Read, Show, Bounded, Generic, Generic1, Num)--derivingUnbox "MinPlus"-  [t| forall x . Unbox x ⇒ MinPlus x → x |]  [| getMinPlus |]  [| MinPlus |]--instance NFData x ⇒ NFData (MinPlus x) where-  rnf (MinPlus x) = rnf x-  {-# Inline rnf #-}---- |------ TODO Shall we have generic instances, or specific ones like @SemiRing--- (Viterbi Prob)@?------ TODO Consider either a constraint @ProbLike x@ or the above.--instance (Ord x, Num x, NumericLimits x) ⇒ SemiRing (MinPlus x) where-  srplus (MinPlus x) (MinPlus y) = MinPlus $ min x y-  srmul  (MinPlus x) (MinPlus y) = MinPlus $ x + y-  srzero = MinPlus maxFinite-  srone  = 0-  {-# Inline srplus #-}-  {-# Inline srmul  #-}-  {-# Inline srzero #-}-  {-# Inline srone  #-}------ | The tropical MaxPlus SemiRing. It maximizes over the sum.--newtype MaxPlus x = MaxPlus { getMaxPlus ∷ x }-  deriving (Eq, Ord, Read, Show, Bounded, Generic, Generic1, Num)--derivingUnbox "MaxPlus"-  [t| forall x . Unbox x ⇒ MaxPlus x → x |]  [| getMaxPlus |]  [| MaxPlus |]--instance NFData x ⇒ NFData (MaxPlus x) where-  rnf (MaxPlus x) = rnf x-  {-# Inline rnf #-}--instance NumericLimits x ⇒ NumericLimits (MaxPlus x) where-  minFinite = MaxPlus minFinite-  maxFinite = MaxPlus maxFinite---- |------ TODO Shall we have generic instances, or specific ones like @SemiRing--- (Viterbi Prob)@?------ TODO Consider either a constraint @ProbLike x@ or the above.--instance (Ord x, Num x, NumericLimits x) ⇒ SemiRing (MaxPlus x) where-  srplus (MaxPlus x) (MaxPlus y) = MaxPlus $ max x y-  srmul  (MaxPlus x) (MaxPlus y) = MaxPlus $ x + y-  srzero = MaxPlus minFinite-  srone  = 0-  {-# Inline srplus #-}-  {-# Inline srmul  #-}-  {-# Inline srzero #-}-  {-# Inline srone  #-}------ * Generic semiring structure encoding.---- | The generic semiring, defined over two 'Semigroup' and 'Monoid'--- constructions.------ It can be used like this:--- @--- srzero ∷ GSemiRing Min Sum Int  == maxBound--- srone  ∷ GSemiRing Min Sum Int  == 0--- @------ It is generally useful to still provide explicit instances, since @Min@--- requires a @Bounded@ instance.--newtype GSemiRing (zeroMonoid ∷ * → *) (oneMonoid ∷ * → *) (x ∷ *) = GSemiRing { getSemiRing ∷ x }-  deriving (Eq, Ord, Read, Show, Generic)--instance-  forall zeroMonoid oneMonoid x-  . ( Semigroup (zeroMonoid x)-    , Monoid    (zeroMonoid x)-    , Semigroup ( oneMonoid x)-    , Monoid    ( oneMonoid x)-    )-  ⇒ SemiRing (GSemiRing zeroMonoid oneMonoid x) where-  srplus (GSemiRing x) (GSemiRing y) =-    let x' ∷ zeroMonoid x = unsafeCoerce x-        y' ∷ zeroMonoid x = unsafeCoerce y-    in  unsafeCoerce $ x' <> y'-  srmul (GSemiRing x) (GSemiRing y) =-    let x' ∷ oneMonoid x = unsafeCoerce x-        y' ∷ oneMonoid x = unsafeCoerce y-    in  unsafeCoerce $ x' <> y'-  srzero = unsafeCoerce (mempty ∷ zeroMonoid x)-  srone  = unsafeCoerce (mempty ∷  oneMonoid x)-  {-# Inline srplus #-}-  {-# Inline srmul  #-}-  {-# Inline srzero #-}-  {-# Inline srone  #-}---- ** Variants of 'Semigroup' structures, that use @NumericLimits@ instead of--- @Bounded@.-

@@ -0,0 +1,204 @@++-- | A set with two binary operations, one for addition (@srplus@), one for+-- multiplication (@srmul@). Together with a neutral element for @srplus@,+-- named @srzero@, and one for @srmul@, named @srone@.++module Algebra.Structure.Semiring+  ( module Algebra.Structure.Semiring+  , Data.Semiring.Semiring (..)+  ) where++import Control.DeepSeq (NFData(..))+import Data.Coerce+import Data.Monoid hiding ((<>))+import Data.Semigroup+import Data.Semiring (Semiring(..))+import Data.Vector.Unboxed.Deriving+import Data.Vector.Unboxed (Unbox)+import GHC.Generics+import Numeric.Log+import Unsafe.Coerce++import Numeric.Limits++++-- | Unicode variant of @srplus@.++infixl 6 ⊕+(⊕) ∷ Semiring a ⇒ a → a → a+(⊕) = plus+{-# Inline (⊕) #-}++-- | Unicode variant of @srmul@.++infixl 7 ⊗+(⊗) ∷ Semiring a ⇒ a → a → a+(⊗) = times+{-# Inline (⊗) #-}++++-- * Newtype wrappers for 'SemiRing' that make the semiring to use explicit.+-- This is important, because several types, say Prob(ability) have multiple+-- useful semiring instances.+--+-- 'Data.Monoid' in @base@ provides a number of newtype wrappers (@Sum@,+-- @Product@, etc) for monoids, which have one binary operation and identity.+-- There is, obviously, overlap with the structures constructed here.++-- | The Viterbi SemiRing. It maximizes over the product.++newtype Viterbi x = Viterbi { getViterbi ∷ x }+  deriving (Eq, Ord, Read, Show, Bounded, Generic, Generic1, Num)++derivingUnbox "Viterbi"+  [t| forall x . Unbox x ⇒ Viterbi x → x |]  [| getViterbi |]  [| Viterbi |]++instance NFData x ⇒ NFData (Viterbi x) where+  rnf (Viterbi x) = rnf x+  {-# Inline rnf #-}++-- |+--+-- TODO Shall we have generic instances, or specific ones like @SemiRing+-- (Viterbi Prob)@?+--+-- TODO Consider either a constraint @ProbLike x@ or the above.++instance (Ord x, Semiring x) ⇒ Semiring (Viterbi x) where+  plus  (Viterbi x) (Viterbi y) = Viterbi $ max x y+  times (Viterbi x) (Viterbi y) = Viterbi $ x `times` y+  zero = Viterbi zero+  one  = Viterbi one+  {-# Inline plus  #-}+  {-# Inline times #-}+  {-# Inline zero  #-}+  {-# Inline one   #-}++-- | The tropical MinPlus SemiRing. It minimizes over the sum.++newtype MinPlus x = MinPlus { getMinPlus ∷ x }+  deriving (Eq, Ord, Read, Show, Bounded, Generic, Generic1, Num)++derivingUnbox "MinPlus"+  [t| forall x . Unbox x ⇒ MinPlus x → x |]  [| getMinPlus |]  [| MinPlus |]++instance NFData x ⇒ NFData (MinPlus x) where+  rnf (MinPlus x) = rnf x+  {-# Inline rnf #-}++instance NumericLimits x ⇒ NumericLimits (MinPlus x) where+  minFinite = MinPlus minFinite+  maxFinite = MinPlus maxFinite++-- |+--+-- Be careful, if the numeric limits are hits, underflows, etc will happen.++instance (Ord x, Semiring x, NumericLimits x) ⇒ Semiring (MinPlus x) where+  plus  (MinPlus x) (MinPlus y) = MinPlus $ min x y+  times (MinPlus x) (MinPlus y) = MinPlus $ x `plus` y+  zero = MinPlus maxFinite+  one  = MinPlus zero+  {-# Inline plus  #-}+  {-# Inline times #-}+  {-# Inline zero  #-}+  {-# Inline one   #-}++++-- | The tropical MaxPlus SemiRing. It maximizes over the sum.++newtype MaxPlus x = MaxPlus { getMaxPlus ∷ x }+  deriving (Eq, Ord, Read, Show, Bounded, Generic, Generic1, Num)++derivingUnbox "MaxPlus"+  [t| forall x . Unbox x ⇒ MaxPlus x → x |]  [| getMaxPlus |]  [| MaxPlus |]++instance NFData x ⇒ NFData (MaxPlus x) where+  rnf (MaxPlus x) = rnf x+  {-# Inline rnf #-}++instance NumericLimits x ⇒ NumericLimits (MaxPlus x) where+  minFinite = MaxPlus minFinite+  maxFinite = MaxPlus maxFinite++-- |+--+-- TODO Shall we have generic instances, or specific ones like @SemiRing+-- (Viterbi Prob)@?+--+-- TODO Consider either a constraint @ProbLike x@ or the above.++instance (Ord x, Semiring x, NumericLimits x) ⇒ Semiring (MaxPlus x) where+  plus  (MaxPlus x) (MaxPlus y) = MaxPlus $ max x y+  times (MaxPlus x) (MaxPlus y) = MaxPlus $ x `plus` y+  zero = MaxPlus minFinite+  one  = MaxPlus zero+  {-# Inline plus  #-}+  {-# Inline times #-}+  {-# Inline zero  #-}+  {-# Inline one   #-}++++-- * Generic semiring structure encoding.++-- | The generic semiring, defined over two 'Semigroup' and 'Monoid'+-- constructions.+--+-- It can be used like this:+-- @+-- zero ∷ GSemiring Min Sum Int  == maxBound+-- one  ∷ GSemiring Min Sum Int  == 0+-- @+--+-- It is generally useful to still provide explicit instances, since @Min@+-- requires a @Bounded@ instance.++newtype GSemiring (zeroMonoid ∷ * → *) (oneMonoid ∷ * → *) (x ∷ *) = GSemiring { getSemiring ∷ x }+  deriving (Eq, Ord, Read, Show, Generic)++instance+  forall zeroMonoid oneMonoid x+  . ( Semigroup (zeroMonoid x)+    , Monoid    (zeroMonoid x)+    , Semigroup ( oneMonoid x)+    , Monoid    ( oneMonoid x)+    , Coercible (zeroMonoid x) (GSemiring zeroMonoid oneMonoid x)+    , Coercible (oneMonoid x) (GSemiring zeroMonoid oneMonoid x)+    )+  ⇒ Semiring (GSemiring zeroMonoid oneMonoid x) where+  plus (GSemiring x) (GSemiring y) =+    let x' ∷ zeroMonoid x = coerce x+        y' ∷ zeroMonoid x = coerce y+    in  coerce $ x' <> y'+  times (GSemiring x) (GSemiring y) =+    let x' ∷ oneMonoid x = coerce x+        y' ∷ oneMonoid x = coerce y+    in  coerce $ x' <> y'+  zero = coerce (mempty ∷ zeroMonoid x)+  one  = coerce (mempty ∷  oneMonoid x)+  {-# Inline plus  #-}+  {-# Inline times #-}+  {-# Inline zero  #-}+  {-# Inline one   #-}++++-- * Semiring on 'Numeric.Log'. This is an orphan instance, but it can't be+-- helped much, unless we want to wrap into yet another newtype.++instance (Precise a, RealFloat a) ⇒ Semiring (Log a) where+  plus  = (+)+  times = (*)+  zero  = 0+  one   = 1+  {-# Inline plus  #-}+  {-# Inline times #-}+  {-# Inline zero  #-}+  {-# Inline one   #-}+++

@@ -5,20 +5,64 @@ module Numeric.Discretized where  import Control.Applicative+import Control.DeepSeq (NFData(..))+import Data.Aeson (FromJSON,ToJSON)+import Data.Binary (Binary)+import Data.Hashable (Hashable) import Data.Proxy import Data.Ratio+import Data.Serialize (Serialize)+import Data.Vector.Unboxed.Deriving import Debug.Trace import GHC.Generics-import GHC.TypeLits import GHC.Real (Ratio(..))+import GHC.TypeLits +import Algebra.Structure.Semiring+import Numeric.Limits  --- | A discretized value takes a floating point number @n@ and produces @n *--- fromIntegral l / fromIntegral u@ where both @u@ and @l@ are given as--- @TypeLits@. I.e. a scaling factor of @ (u / l) = (1 / 100)@ does all--- calculations in subdivisions of 100.++-- | Some discretizations are of the type @ln 2 / 2@ (@PAM@ matrices in Blast+-- for example). Using this type, we can annotate as follows: @Discretized+-- (RTyLn 2 :% RTyId 2)@. --+-- One may use @Unknown@ if the scale is not known. For example, the blast+-- matrices use different scales internally and one needs to read the header to+-- get the scale.++data RatioTy a = RTyExp a | RTyId a | RTyLn a | RTyPlus (RatioTy a) (RatioTy a) | RTyTimes (RatioTy a) (RatioTy a) | Unknown++class RatioTyConstant a where+  ratioTyConstant ∷ Proxy a → Ratio Integer++instance (KnownNat k) ⇒ RatioTyConstant (RTyExp (k∷Nat)) where+  {-# Inline ratioTyConstant #-}+  ratioTyConstant Proxy = let n = natVal @k Proxy in toRational (exp $ fromInteger n)++instance (KnownNat k) ⇒ RatioTyConstant (RTyId (k∷Nat)) where+  {-# Inline ratioTyConstant #-}+  ratioTyConstant Proxy = let n = natVal @k Proxy in toRational n++instance (KnownNat k) ⇒ RatioTyConstant (RTyLn (k∷Nat)) where+  {-# Inline ratioTyConstant #-}+  ratioTyConstant Proxy = let n = natVal @k Proxy in toRational (log $ fromInteger n)++instance (RatioTyConstant a, RatioTyConstant b) ⇒ RatioTyConstant (RTyPlus (a∷RatioTy k) (b∷RatioTy k)) where+  {-# Inline ratioTyConstant #-}+  ratioTyConstant Proxy = ratioTyConstant @a Proxy + ratioTyConstant @b Proxy++instance (RatioTyConstant a, RatioTyConstant b) ⇒ RatioTyConstant (RTyTimes (a∷RatioTy k) (b∷RatioTy k)) where+  {-# Inline ratioTyConstant #-}+  ratioTyConstant Proxy = ratioTyConstant @a Proxy * ratioTyConstant @b Proxy++-- | A discretized value takes a floating point number @n@ and produces a+-- discretized value. The actual discretization formula is given on the type+-- level, freeing us from having to carry around some scaling function.+--+-- Typically, one might use types likes @100@, @(100 :% 1)@, or @(RTyLn 2 :%+-- RTyId 2)@.+-- -- The main use of a 'Discretized' value is to enable calculations with 'Int' -- while somewhat pretending to use floating point values. --@@ -36,65 +80,112 @@ -- some thought on in which direction to wrap. Maybe, we want to log-domain -- Discretized values, which probably just works. -newtype Discretized (u ∷ Nat) (l ∷ Nat) = Discretized { getDiscretized ∷ Int }+newtype Discretized (b ∷ k) = Discretized { getDiscretized ∷ Int }   deriving (Eq,Ord,Generic,Show,Read) -instance (KnownNat u, KnownNat l) ⇒ Num (Discretized u l) where+derivingUnbox "Discretized"+  [t| forall t . Discretized t → Int |]  [| getDiscretized |]  [| Discretized |]++instance NFData (Discretized t) where+  rnf (Discretized k) = rnf k+  {-# Inline rnf #-}++instance Binary    (Discretized t)+instance Serialize (Discretized t)+instance FromJSON  (Discretized t)+instance ToJSON    (Discretized t)+instance Hashable  (Discretized t)++instance Num (Discretized Unknown) where+  Discretized x + Discretized y = Discretized $ x+y+  Discretized x - Discretized y = Discretized $ x-y+  (*) = error "Discretized Unknown does not admit (*)"+  abs (Discretized x) = Discretized $ abs x+  signum (Discretized x) = Discretized $ signum x+  fromInteger = error "Discretized Unknown does not admit fromInteger"+  {-# Inline (+) #-}+  {-# Inline (-) #-}+  {-# Inline abs #-}+  {-# Inline signum #-}++instance (KnownNat u, KnownNat l) ⇒ Num (Discretized ((u∷Nat) :% (l∷Nat))) where+  {-# Inline (+) #-}   Discretized x + Discretized y = Discretized (x+y)+  {-# Inline (-) #-}   Discretized x - Discretized y = Discretized (x-y)+  -- TODO it should be possible to generalize this over arbitrary value, or+  -- replace @KnownNat@ with the above @ratioTyConstant@.+  {-# Inline (*) #-}   Discretized x * Discretized y =     let u = fromInteger $ natVal @u Proxy         l = fromInteger $ natVal @l Proxy     in  Discretized $ (x*y*u) `div` l-  abs (Discretized x) = Discretized (abs x)-  signum (Discretized x) = Discretized $ signum x-  fromInteger = Discretized . fromInteger-  {-# Inline (+) #-}-  {-# Inline (-) #-}-  {-# Inline (*) #-}   {-# Inline abs #-}+  abs (Discretized x) = Discretized (abs x)   {-# Inline signum #-}+  signum (Discretized x) = Discretized $ signum x   {-# Inline fromInteger #-}+  fromInteger x =+    let u = fromInteger $ natVal @u Proxy+        l = fromInteger $ natVal @l Proxy+    in  Discretized $ (fromInteger x*u) `div` l -instance Enum (Discretized u l) where+instance Enum (Discretized b) where   toEnum = Discretized   {-# Inline toEnum #-}   fromEnum = getDiscretized   {-# Inline fromEnum #-} -instance (Enum (Discretized u l), KnownNat u, KnownNat l) ⇒ Integral (Discretized u l) where+-- instance (Enum (Discretized b), KnownNat u, KnownNat l) ⇒ Integral (Discretized u l) where -instance (KnownNat u, KnownNat l) ⇒ Fractional (Discretized u l) where+instance (KnownNat u, KnownNat l) ⇒ Fractional (Discretized ((u∷Nat) :% (l∷Nat))) where+  {-# Inline (/) #-}   Discretized x / Discretized y =     let u = fromInteger $ natVal @u Proxy         l = fromInteger $ natVal @l Proxy     in  Discretized $ (x * l) `div` (y * u)-  {-# Inline (/) #-}+  {-# Inline recip #-}   recip (Discretized x) =     let u = fromInteger $ natVal @u Proxy         l = fromInteger $ natVal @l Proxy     in  error "need to find approximately ok transformation"-  {-# Inline recip #-}+  {-# Inline fromRational #-}   fromRational (a :% b) =     let u = natVal @u Proxy         l = natVal @l Proxy     in  Discretized . fromInteger $ (a * l) `div` (b * u) -instance (KnownNat u, KnownNat l) ⇒ Real (Discretized u l) where+instance (KnownNat u, KnownNat l) ⇒ Real (Discretized ((u∷Nat) :% (l∷Nat))) where+  {-# Inline toRational #-}   toRational (Discretized d) =     let u = natVal @u Proxy         l = natVal @l Proxy     in  (fromIntegral d * u) % l-  {-# Inline toRational #-} +instance (Num (Discretized k)) ⇒ Semiring (Discretized k) where+  plus = (+)+  times = (*)+  zero = 0+  one = 1+  {-# Inline plus  #-}+  {-# Inline times #-}+  {-# Inline zero  #-}+  {-# Inline one   #-}++instance (NumericLimits (Discretized t)) where+  minFinite = Discretized minFinite+  {-# Inline minFinite #-}+  maxFinite = Discretized maxFinite+  {-# Inline maxFinite #-}+ -- | Discretizes any @Real a@ into the @Discretized@ value. This conversion--- is /lossy/!+-- is /lossy/ and uses a type-level rational of @u :% l@! -discretize ∷ forall a u l . (Real a, KnownNat u, KnownNat l) ⇒ a → Discretized u l-discretize a =+discretizeRatio ∷ forall a u l . (Real a, KnownNat u, KnownNat l) ⇒ a → Discretized ((u∷Nat) :% (l∷Nat))+{-# Inline discretizeRatio #-}+discretizeRatio a =   let u = natVal @u Proxy       l = natVal @l Proxy       k = toRational a   in  Discretized . fromIntegral $ numerator k * l `div` (denominator k * u)-{-# Inline discretize #-} 

@@ -7,15 +7,16 @@ module Numeric.LogDomain where  import Control.Monad.Except+import Numeric.Log    -- | Instances for @LogDomain x@ should be for specific types.  class LogDomain x where-  -- | The data family to connect a type @x@ with the type @Ln x@ in the+  -- | The type family to connect a type @x@ with the type @Ln x@ in the   -- log-domain.-  data Ln x ∷ *+  type Ln x ∷ *   -- | Transport a value in @x@ into the log-domain. @logdom@ should throw an   -- exception if @log x@ is not valid.   logdom ∷ (MonadError String m) ⇒ x → m (Ln x)@@ -23,4 +24,17 @@   unsafelogdom ∷ x → Ln x   -- | Transport a value @Ln x@ back into the linear domain @x@.   lindom ∷ Ln x → x++++instance LogDomain Double where+  type Ln Double = Log Double+  {-# Inline logdom #-}+  logdom x+    | x < 0     = throwError "log of negative number"+    | otherwise = return $ unsafelogdom x+  {-# Inline unsafelogdom #-}+  unsafelogdom = Exp . log+  {-# Inline lindom #-}+  lindom = exp . ln 

@@ -1,16 +1,16 @@+Cabal-version:  2.2 Name:           SciBaseTypes-Version:        0.0.0.1-License:        BSD3+Version:        0.1.0.0+License:        BSD-3-Clause License-file:   LICENSE Author:         Christian Hoener zu Siederdissen Maintainer:     choener@bioinf.uni-leipzig.de-Copyright:      Christian Hoener zu Siederdissen, 2018+Copyright:      Christian Hoener zu Siederdissen, 2018-2019 homepage:       https://github.com/choener/SciBaseTypes bug-reports:    https://github.com/choener/SciBaseTypes/issues Stability:      Experimental Category:       Data Build-type:     Simple-Cabal-version:  >= 1.10 tested-with:    GHC == 8.4.4 Synopsis:       Base types and classes for statistics, sciences and humanities Description:@@ -25,16 +25,7 @@   -library-  exposed-modules:-    Algebra.Structure.SemiRing-    Numeric.Discretized-    Numeric.Limits-    Numeric.LogDomain-    StatisticalMechanics.Ensemble-    Statistics.Odds-    Statistics.Probability-+common deps   build-depends: base                     >= 4.7      &&  < 5.0                , aeson                    >= 1.0                , binary                   >= 0.7@@ -44,18 +35,20 @@                , lens                     >= 4.0                , log-domain               >= 0.12                , mtl                      >= 2.0+               , semirings                >= 0.3                , vector                   >= 0.10                , vector-th-unbox          >= 0.2+  default-language:+    Haskell2010   ghc-options:     -O2     -funbox-strict-fields-  default-language:-    Haskell2010   default-extensions: BangPatterns                     , ConstraintKinds                     , DataKinds                     , DeriveGeneric                     , FlexibleContexts+                    , FlexibleInstances                     , GeneralizedNewtypeDeriving                     , MultiParamTypeClasses                     , PolyKinds@@ -66,51 +59,51 @@                     , TupleSections                     , TypeApplications                     , TypeFamilies+                    , TypeOperators                     , UndecidableInstances                     , UnicodeSyntax   +library+  import: deps+  exposed-modules:+    Algebra.Structure.Semiring+    Numeric.Discretized+    Numeric.Limits+    Numeric.LogDomain+    StatisticalMechanics.Ensemble+    Statistics.Odds+    Statistics.Probability+++ benchmark BenchmarkSciBaseTypes+  import: deps   build-depends: base                , SciBaseTypes   hs-source-dirs:     tests   main-is:     Benchmark.hs-  default-language:-    Haskell2010   type:     exitcode-stdio-1.0-  default-extensions: BangPatterns-                    , FlexibleContexts-                    , ScopedTypeVariables-                    , TypeFamilies-                    , UnicodeSyntax   ghc-options:-    -O2-    -funbox-strict-fields     -funfolding-use-threshold1000     -funfolding-keeness-factor1000    test-suite properties+  import: deps   type:     exitcode-stdio-1.0   main-is:     properties.hs   ghc-options:-    -threaded -rtsopts -with-rtsopts=-N+    -O0 -threaded -rtsopts -with-rtsopts=-N   hs-source-dirs:     tests-  default-language:-    Haskell2010-  default-extensions: BangPatterns-                    , ScopedTypeVariables-                    , TemplateHaskell-                    , TypeFamilies-                    , UnicodeSyntax   build-depends: base                , SciBaseTypes                , QuickCheck                   >= 2.7

@@ -4,8 +4,12 @@  module StatisticalMechanics.Ensemble where +import Numeric.Log +import Statistics.Probability ++ -- | The state probability functions provide conversion from some types @a@ -- into non-normalized probabilities. For "real" applications, using the -- @logProbability@ function is preferred. This functions allows for easy@@ -17,11 +21,31 @@ -- feasible in a number of problems. -- -- TODO replace @()@ with temperature and results with non-normalized @P@ or--- @LogP@, depending.+-- @LogP@, depending. At some point we want to have type-level physical+-- quantities, hence the need for the second type.  class StateProbability a where   -- | Given a temperature and a state "energy", return the corresponding   -- non-normalized probability.-  stateProbability    ∷ () → a → ()-  stateLogProbability ∷ () → a → ()+  stateProbability+    ∷ Double+    -- ^ this is @k*T@+    → a+    -- ^ the energy (or discretized energy)+    → Probability NotNormalized Double+    -- ^ probability of being in state @a@, but only proportional up to @1/Z@.+  stateLogProbability+    ∷ Double+    -- ^ this is @1/(k * T)@+    → a+    -- ^ the energy (or discretized energy)+    → Log (Probability NotNormalized Double)+    -- ^ resulting probability++instance StateProbability Double where+  stateProbability kT x = Prob . exp . negate $ x/kT+  {-# Inline stateProbability #-}+  --stateLogProbability kT x = Exp . log . Prob . exp . negate $ x/kT+  stateLogProbability kT x = Exp . Prob . negate $ x/kT+  {-# Inline stateLogProbability #-} 

@@ -1,5 +1,8 @@  -- | Provides newtypes for odds, log-odds, and discretized versions.+--+-- TODO This is currently quite ad-hoc and needs better formalization. In+-- particular in terms of wrapping and usage of @Num@ and @Semiring@.  module Statistics.Odds where @@ -11,6 +14,8 @@ import Data.Vector.Unboxed.Deriving import GHC.Generics (Generic) +import Algebra.Structure.Semiring+import Numeric.Discretized import Numeric.Limits  @@ -20,30 +25,39 @@ newtype Odds = Odds { getOdds ∷ Double }   deriving (Generic,Eq,Ord,Show,Read,Num) +deriving instance Semiring Odds+++ -- | Encodes log-odds that have been rounded or clamped to integral numbers. -- One advantage this provides is more efficient "maximum/minimum" calculations -- compared to using @Double@s. -- -- Note that these are "explicit" log-odds. Each numeric operation uses the--- underlying operation on @Int@.+-- underlying operation on @Int@. If you want automatic handling, choose @Log+-- Odds@. -newtype DiscLogOdds = DiscLogOdds { getDiscLogOdds ∷ Int }-  deriving (Generic,Eq,Ord,Show,Read,Num)+newtype DiscLogOdds (t∷k) = DiscLogOdds { getDiscLogOdds ∷ Discretized t }+  deriving (Generic,Eq,Ord,Show,Read) +deriving instance (Num (Discretized (t∷k))) ⇒ Num (DiscLogOdds t)++deriving instance (Semiring (Discretized (t∷k))) ⇒ Semiring (DiscLogOdds t)+ derivingUnbox "DiscretizedLogOdds"-  [t| DiscLogOdds → Int |]  [| getDiscLogOdds |]  [| DiscLogOdds |]+  [t| forall t . DiscLogOdds t → Int |]  [| getDiscretized . getDiscLogOdds |]  [| DiscLogOdds . Discretized |] -instance Binary    DiscLogOdds-instance Serialize DiscLogOdds-instance FromJSON  DiscLogOdds-instance ToJSON    DiscLogOdds-instance Hashable  DiscLogOdds+instance Binary    (DiscLogOdds t)+instance Serialize (DiscLogOdds t)+instance FromJSON  (DiscLogOdds t)+instance ToJSON    (DiscLogOdds t)+instance Hashable  (DiscLogOdds t) -instance NFData DiscLogOdds where+instance (NFData (Discretized t)) ⇒ NFData (DiscLogOdds t) where   rnf (DiscLogOdds k) = rnf k   {-# Inline rnf #-} -instance NumericLimits DiscLogOdds where+instance (NumericLimits (Discretized t)) ⇒ NumericLimits (DiscLogOdds t) where   minFinite = DiscLogOdds minFinite   {-# Inline minFinite #-}   maxFinite = DiscLogOdds maxFinite

@@ -7,11 +7,12 @@ module Statistics.Probability where  import Control.Lens-import Numeric.Log+import Data.Char (chr) import Data.Vector.Unboxed.Deriving import Data.Vector.Unboxed (Unbox)+import Numeric.Log -import Algebra.Structure.SemiRing+import Algebra.Structure.Semiring import Numeric.LogDomain import Numeric.Limits @@ -27,97 +28,128 @@ -- @Normalized@, the contained values are in the range @[0,...,1]@, otherwise -- they are in the range @[0,...,∞]@. -newtype Prob (n ∷ IsNormalized) x = Prob { getProb ∷ x }+newtype Probability (n ∷ IsNormalized) x = Prob { getProb ∷ x }   deriving (Eq,Ord,Show,Read) -derivingUnbox "Prob"-  [t| forall n x. Unbox x ⇒ Prob n x → x |]  [| getProb |]  [| Prob |]+derivingUnbox "Probability"+  [t| forall n x. Unbox x ⇒ Probability n x → x |]  [| getProb |]  [| Prob |] -deriving instance (Enum       x) ⇒ Enum       (Prob n x)-deriving instance (Num        x) ⇒ Num        (Prob n x)-deriving instance (Fractional x) ⇒ Fractional (Prob n x)-deriving instance (Floating   x) ⇒ Floating   (Prob n x)-deriving instance (Real       x) ⇒ Real       (Prob n x)-deriving instance (RealFrac   x) ⇒ RealFrac   (Prob n x)-deriving instance (RealFloat  x) ⇒ RealFloat  (Prob n x)+deriving instance (Enum       x) ⇒ Enum       (Probability n x)+deriving instance (Num        x) ⇒ Num        (Probability n x)+deriving instance (Fractional x) ⇒ Fractional (Probability n x)+deriving instance (Floating   x) ⇒ Floating   (Probability n x)+deriving instance (Real       x) ⇒ Real       (Probability n x)+deriving instance (RealFrac   x) ⇒ RealFrac   (Probability n x)+deriving instance (RealFloat  x) ⇒ RealFloat  (Probability n x)+deriving instance (Precise    x) ⇒ Precise    (Probability n x) -instance (Num r) ⇒ SemiRing (Prob n r) where-  srplus = (+)-  srmul  = (*)-  srzero = 0-  srone  = 1+instance (Num r) ⇒ Semiring (Probability n r) where+  plus  = (+)+  times = (*)+  zero = 0+  one  = 1+  {-# Inline plus  #-}+  {-# Inline times #-}+  {-# Inline zero  #-}+  {-# Inline one   #-}  -- | Turns a value into a normalized probability. @error@ if the value is not -- in the range @[0,...,1]@. -prob ∷ (Ord x, Num x, Show x) ⇒ x → Prob Normalized x+prob ∷ (Ord x, Num x, Show x) ⇒ x → Probability Normalized x prob x   | x >= 0 && x <= 1 = Prob x   | otherwise        = error $ show x ++ " not in range of [0,...,1]" {-# Inline prob #-} --- | Simple wrapper around @Prob@ that fixes non-normalization.+-- | Simple wrapper around @Probability@ that fixes non-normalization. -prob' ∷ (Ord x, Num x, Show x) ⇒ x → Prob NotNormalized x+prob' ∷ (Ord x, Num x, Show x) ⇒ x → Probability NotNormalized x prob' = Prob {-# Inline prob' #-} ----- * Probability in log space. A number of operations internally cast to @Log@--- from @log-domain@, but the values themselves are *not* newtype-wrapped @Log--- x@ values. This is because we want to be *explicit* that these are--- log-probabilities.------ @Log@ numbers in cases like @fromIntegral 1 :: Log Double@ are treated as--- not being in the log-domain, hence @fromIntegral performs a @log@--- operations.--newtype LogProb (n ∷ IsNormalized) x = LogProb { getLogProb ∷ x }-  deriving (Eq,Ord,Show)--derivingUnbox "LogProb"-  [t| forall n x. Unbox x ⇒ LogProb n x → x |]  [| getLogProb |]  [| LogProb |]--instance (Precise x, RealFloat x) ⇒ Num (LogProb n x) where-  (+) = withLog2 (+)-  (*) = withLog2 (*)-  abs = withLog1 abs-  signum = withLog1 signum-  fromInteger = LogProb . fromInteger-  negate = withLog1 negate-  (-) = withLog2 (-)--instance (Num d, Fractional d) ⇒ NumericLimits (LogProb n d) where-  minFinite = LogProb 0-  maxFinite = LogProb (1/0)--withLog1 ∷ (Log x → Log y) → LogProb n x → LogProb n y-withLog1 op (LogProb x) = LogProb . ln $ op (Exp x)-{-# Inline withLog1 #-}+-- | This simple function represents probabilities with characters between '0'+-- @0.0 -- 0.05@ up to '9' @0.85 -- 0.95@ and finally '*' for @>0.95@. -withLog2 ∷ (Log x → Log y → Log z) → LogProb n x → LogProb n y → LogProb n z-withLog2 op (LogProb x) (LogProb y) = LogProb . ln $ op (Exp x) (Exp y)-{-# Inline withLog2 #-}+probabilityToChar ∷ (Num k, RealFrac k) ⇒ Probability Normalized k → Char+probabilityToChar (Prob p')+  | i >= 10 = '*'+  | otherwise = chr $ 48 + i+  where p = max 0.0 $ min p' 1.0+        i = round $ p * 10  --- * Conversion between probability in linear and log-space.---- | Turn probability into log-probability.--p2lp ∷ (Floating x) ⇒ Prob n x → LogProb n x-p2lp (Prob x) = LogProb $ log x-{-# Inline p2lp #-}---- | Turn log-probability into probability.--lp2p ∷ (Floating x) ⇒ LogProb n x → Prob n x-lp2p (LogProb x) = Prob $ exp x-{-# Inline lp2p #-}---- | An isomorphism between @Prob@ and @LogProb@.--aslp ∷ (Floating x) ⇒ Iso' (Prob n x) (LogProb n x)-aslp = iso p2lp lp2p-{-# Inline aslp #-}+-- -- * Probability in log space. A number of operations internally cast to @Log@+-- -- from @log-domain@, but the values themselves are *not* newtype-wrapped @Log+-- -- x@ values. This is because we want to be *explicit* that these are+-- -- log-probabilities.+-- --+-- -- @Log@ numbers in cases like @fromIntegral 1 :: Log Double@ are treated as+-- -- not being in the log-domain, hence @fromIntegral performs a @log@+-- -- operations.+-- +-- newtype LogProb (n ∷ IsNormalized) x = LogProb { getLogProb ∷ x }+--   deriving (Eq,Ord,Show)+-- +-- derivingUnbox "LogProb"+--   [t| forall n x. Unbox x ⇒ LogProb n x → x |]  [| getLogProb |]  [| LogProb |]+-- +-- instance (Precise x, RealFloat x) ⇒ Num (LogProb n x) where+--   (+) = withLog2 (+)+--   {-# Inline (+) #-}+--   (*) = withLog2 (*)+--   {-# Inline (*) #-}+--   abs = withLog1 abs+--   {-# Inline abs #-}+--   signum = withLog1 signum+--   {-# Inline signum #-}+--   fromInteger = LogProb . fromInteger+--   {-# Inline fromInteger #-}+--   negate = withLog1 negate+--   {-# Inline negate #-}+--   (-) = withLog2 (-)+--   {-# Inline (-) #-}+-- +-- instance (Precise r, RealFloat r, Num r) ⇒ SemiRing (LogProb n r) where+--   srplus = (+)+--   {-# Inline srplus #-}+--   srmul  = (*)+--   {-# Inline srmul #-}+--   srzero = 0+--   {-# Inline srzero #-}+--   srone  = 1+--   {-# Inline srone #-}+-- +-- instance (Num d, Fractional d) ⇒ NumericLimits (LogProb n d) where+--   minFinite = LogProb 0+--   maxFinite = LogProb (1/0)+-- +-- withLog1 ∷ (Log x → Log y) → LogProb n x → LogProb n y+-- withLog1 op (LogProb x) = LogProb . ln $ op (Exp x)+-- {-# Inline withLog1 #-}+-- +-- withLog2 ∷ (Log x → Log y → Log z) → LogProb n x → LogProb n y → LogProb n z+-- withLog2 op (LogProb x) (LogProb y) = LogProb . ln $ op (Exp x) (Exp y)+-- {-# Inline withLog2 #-}+-- +-- +-- -- * Conversion between probability in linear and log-space.+-- +-- -- | Turn probability into log-probability.+-- +-- p2lp ∷ (Floating x) ⇒ Prob n x → LogProb n x+-- p2lp (Prob x) = LogProb $ log x+-- {-# Inline p2lp #-}+-- +-- -- | Turn log-probability into probability.+-- +-- lp2p ∷ (Floating x) ⇒ LogProb n x → Prob n x+-- lp2p (LogProb x) = Prob $ exp x+-- {-# Inline lp2p #-}+-- +-- -- | An isomorphism between @Prob@ and @LogProb@.+-- +-- aslp ∷ (Floating x) ⇒ Iso' (Prob n x) (LogProb n x)+-- aslp = iso p2lp lp2p+-- {-# Inline aslp #-} 

@@ -1,4 +1,8 @@-0.0.0.1+0.1.0.0 -------+- depending on @semirings@ now, instead of our own type class. breaks existing+  code, but harmless to move to. +0.0.0.1+------- - initial creation