Acceptance.hs

{-# LANGUAGE TemplateHaskell #-}

-- |
-- Module      :  Mcmc.Acceptance
-- Description :  Handle acceptance rates
-- Copyright   :  2021 Dominik Schrempf
-- License     :  GPL-3.0-or-later
--
-- Maintainer  :  dominik.schrempf@gmail.com
-- Stability   :  experimental
-- Portability :  portable
--
-- Creation date: Thu Jul  8 18:12:07 2021.
module Mcmc.Acceptance
  ( -- * Acceptance rates
    AcceptanceRate,
    AcceptanceCounts (..),
    AcceptanceRates (..),
    Acceptance,
    Acceptances (fromAcceptances),
    emptyA,
    pushAccept,
    pushReject,
    ResetAcceptance (..),
    resetA,
    transformKeysA,
    acceptanceRate,
    acceptanceRates,
  )
where

import Data.Aeson
import Data.Aeson.TH
import Data.Foldable
import qualified Data.Map.Strict as M

-- | Acceptance rate.
type AcceptanceRate = Double

-- | Number of accepted and rejected proposals.
data AcceptanceCounts = AcceptanceCounts
  { nAccepted :: !Int,
    nRejected :: !Int
  }
  deriving (Show, Eq, Ord)

$(deriveJSON defaultOptions ''AcceptanceCounts)

-- | Proposals based on Hamiltonian dynamics use expected acceptance rates, not counts.
data AcceptanceRates = AcceptanceRates
  { totalAcceptanceRate :: !Double,
    nAcceptanceRates :: !Int
  }
  deriving (Show, Eq)

$(deriveJSON defaultOptions ''AcceptanceRates)

-- | Stored actual acceptance counts and maybe expected acceptance rates.
data Acceptance = A AcceptanceCounts (Maybe AcceptanceRates)
  deriving (Show, Eq)

$(deriveJSON defaultOptions ''Acceptance)

addAccept :: Maybe AcceptanceRates -> Acceptance -> Acceptance
addAccept mr' (A (AcceptanceCounts a r) mr) = A (AcceptanceCounts (a + 1) r) (addAcceptanceRates mr' mr)

addReject :: Maybe AcceptanceRates -> Acceptance -> Acceptance
addReject mr' (A (AcceptanceCounts a r) mr) = A (AcceptanceCounts a (r + 1)) (addAcceptanceRates mr' mr)

addAcceptanceRates :: Maybe AcceptanceRates -> Maybe AcceptanceRates -> Maybe AcceptanceRates
addAcceptanceRates Nothing Nothing = Nothing
addAcceptanceRates (Just r) Nothing = Just r
addAcceptanceRates Nothing (Just r) = Just r
addAcceptanceRates (Just (AcceptanceRates al rl)) (Just (AcceptanceRates ar rr)) =
  Just $ AcceptanceRates (al + ar) (rl + rr)

-- | For each key @k@, store the number of accepted and rejected proposals.
newtype Acceptances k = Acceptances {fromAcceptances :: M.Map k Acceptance}
  deriving (Eq, Show)

instance (ToJSONKey k) => ToJSON (Acceptances k) where
  toJSON (Acceptances m) = toJSON m
  toEncoding (Acceptances m) = toEncoding m

instance (Ord k, FromJSONKey k) => FromJSON (Acceptances k) where
  parseJSON v = Acceptances <$> parseJSON v

-- | In the beginning there was the Word.
--
-- Initialize an empty storage of accepted/rejected values.
emptyA :: (Ord k) => [k] -> Acceptances k
emptyA ks = Acceptances $ M.fromList [(k, A noCounts Nothing) | k <- ks]
  where
    noCounts = AcceptanceCounts 0 0

-- | For key @k@, add an accept.
pushAccept :: (Ord k) => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k
pushAccept mr k = Acceptances . M.adjust (addAccept mr) k . fromAcceptances

-- | For key @k@, add a reject.
pushReject :: (Ord k) => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k
pushReject mr k = Acceptances . M.adjust (addReject mr) k . fromAcceptances

-- | Reset acceptance specification.
data ResetAcceptance
  = -- | Reset actual acceptance counts and expected acceptance rates.
    ResetEverything
  | -- | Only reset expected acceptance rates.
    ResetExpectedRatesOnly

-- | Reset acceptance counts.
resetA :: (Ord k) => ResetAcceptance -> Acceptances k -> Acceptances k
resetA ResetEverything = emptyA . M.keys . fromAcceptances
resetA ResetExpectedRatesOnly = Acceptances . M.map f . fromAcceptances
  where
    f (A cs _) = A cs Nothing

transformKeys :: (Ord k1, Ord k2) => [(k1, k2)] -> M.Map k1 v -> M.Map k2 v
transformKeys ks m = foldl' insrt M.empty ks
  where
    insrt m' (k1, k2) = M.insert k2 (m M.! k1) m'

-- | Transform keys using the given lists. Keys not provided will not be present
-- in the new 'Acceptance' variable.
transformKeysA :: (Ord k1, Ord k2) => [(k1, k2)] -> Acceptances k1 -> Acceptances k2
transformKeysA ks = Acceptances . transformKeys ks . fromAcceptances

-- | Compute acceptance counts, and actual and expected acceptances rates for a
-- specific proposal.
--
-- Return @Just (accepts, rejects, acceptance rate)@.
--
-- Return 'Nothing' if no proposals have been accepted or rejected (division by
-- zero).
acceptanceRate ::
  (Ord k) =>
  k ->
  Acceptances k ->
  -- | (nAccepts, nRejects, actualRate, expectedRate)
  (Int, Int, Maybe AcceptanceRate, Maybe AcceptanceRate)
acceptanceRate k a = case fromAcceptances a M.!? k of
  Just (A (AcceptanceCounts as rs) mrs) -> (as, rs, mar, mtr)
    where
      s = as + rs
      mar = if s <= 0 then Nothing else Just $ fromIntegral as / fromIntegral s
      mtr = case mrs of
        Nothing -> Nothing
        Just (AcceptanceRates xs n) -> Just $ xs / fromIntegral n
  Nothing -> error "acceptanceRate: Key not found in map."

-- | Compute actual acceptance rates for all proposals.
--
-- Set rate to 'Nothing' if no proposals have been accepted or rejected
-- (division by zero).
acceptanceRates :: Acceptances k -> M.Map k (Maybe AcceptanceRate)
acceptanceRates = M.map getRate . fromAcceptances
  where
    getRate (A (AcceptanceCounts as rs) _) =
      let s = as + rs
       in if s <= 0
            then Nothing
            else Just $ fromIntegral as / fromIntegral s