{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RankNTypes #-}

-- |
-- Module      :  Mcmc.Proposal
-- Description :  Proposals and cycles
-- Copyright   :  (c) Dominik Schrempf 2020
-- License     :  GPL-3.0-or-later
--
-- Maintainer  :  dominik.schrempf@gmail.com
-- Stability   :  unstable
-- Portability :  portable
--
-- Creation date: Wed May 20 13:42:53 2020.
module Mcmc.Proposal
  ( -- * Proposal
    PName (..),
    PDescription (..),
    PWeight (..),
    Proposal (..),
    (@~),
    ProposalSimple,
    Tuner (tParam, tFunc),
    Tune (..),
    createProposal,
    tune,

    -- * Cycle
    Order (..),
    Cycle (ccProposals),
    fromList,
    setOrder,
    getNIterations,
    tuneCycle,
    autotuneCycle,
    summarizeCycle,

    -- * Acceptance
    Acceptance (fromAcceptance),
    emptyA,
    pushA,
    resetA,
    transformKeysA,
    acceptanceRatios,
  )
where

import Data.Aeson
import Data.Bifunctor
import qualified Data.ByteString.Builder as BB
import qualified Data.ByteString.Lazy.Char8 as BL
import Data.Default
import qualified Data.Double.Conversion.ByteString as BC
import Data.Function
import Data.List
import Data.Map.Strict (Map)
import qualified Data.Map.Strict as M
import Data.Maybe
import Lens.Micro
import Mcmc.Internal.ByteString
import Mcmc.Internal.Shuffle
import Numeric.Log hiding (sum)
import System.Random.MWC

-- | Proposal name.
newtype PName = PName {PName -> String
fromPName :: String}
  deriving (Int -> PName -> ShowS
[PName] -> ShowS
PName -> String
(Int -> PName -> ShowS)
-> (PName -> String) -> ([PName] -> ShowS) -> Show PName
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [PName] -> ShowS
$cshowList :: [PName] -> ShowS
show :: PName -> String
$cshow :: PName -> String
showsPrec :: Int -> PName -> ShowS
$cshowsPrec :: Int -> PName -> ShowS
Show, PName -> PName -> Bool
(PName -> PName -> Bool) -> (PName -> PName -> Bool) -> Eq PName
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: PName -> PName -> Bool
$c/= :: PName -> PName -> Bool
== :: PName -> PName -> Bool
$c== :: PName -> PName -> Bool
Eq, Eq PName
Eq PName
-> (PName -> PName -> Ordering)
-> (PName -> PName -> Bool)
-> (PName -> PName -> Bool)
-> (PName -> PName -> Bool)
-> (PName -> PName -> Bool)
-> (PName -> PName -> PName)
-> (PName -> PName -> PName)
-> Ord PName
PName -> PName -> Bool
PName -> PName -> Ordering
PName -> PName -> PName
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: PName -> PName -> PName
$cmin :: PName -> PName -> PName
max :: PName -> PName -> PName
$cmax :: PName -> PName -> PName
>= :: PName -> PName -> Bool
$c>= :: PName -> PName -> Bool
> :: PName -> PName -> Bool
$c> :: PName -> PName -> Bool
<= :: PName -> PName -> Bool
$c<= :: PName -> PName -> Bool
< :: PName -> PName -> Bool
$c< :: PName -> PName -> Bool
compare :: PName -> PName -> Ordering
$ccompare :: PName -> PName -> Ordering
$cp1Ord :: Eq PName
Ord)

-- | Proposal description.
newtype PDescription = PDescription {PDescription -> String
fromPDescription :: String}
  deriving (Int -> PDescription -> ShowS
[PDescription] -> ShowS
PDescription -> String
(Int -> PDescription -> ShowS)
-> (PDescription -> String)
-> ([PDescription] -> ShowS)
-> Show PDescription
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [PDescription] -> ShowS
$cshowList :: [PDescription] -> ShowS
show :: PDescription -> String
$cshow :: PDescription -> String
showsPrec :: Int -> PDescription -> ShowS
$cshowsPrec :: Int -> PDescription -> ShowS
Show, PDescription -> PDescription -> Bool
(PDescription -> PDescription -> Bool)
-> (PDescription -> PDescription -> Bool) -> Eq PDescription
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: PDescription -> PDescription -> Bool
$c/= :: PDescription -> PDescription -> Bool
== :: PDescription -> PDescription -> Bool
$c== :: PDescription -> PDescription -> Bool
Eq, Eq PDescription
Eq PDescription
-> (PDescription -> PDescription -> Ordering)
-> (PDescription -> PDescription -> Bool)
-> (PDescription -> PDescription -> Bool)
-> (PDescription -> PDescription -> Bool)
-> (PDescription -> PDescription -> Bool)
-> (PDescription -> PDescription -> PDescription)
-> (PDescription -> PDescription -> PDescription)
-> Ord PDescription
PDescription -> PDescription -> Bool
PDescription -> PDescription -> Ordering
PDescription -> PDescription -> PDescription
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: PDescription -> PDescription -> PDescription
$cmin :: PDescription -> PDescription -> PDescription
max :: PDescription -> PDescription -> PDescription
$cmax :: PDescription -> PDescription -> PDescription
>= :: PDescription -> PDescription -> Bool
$c>= :: PDescription -> PDescription -> Bool
> :: PDescription -> PDescription -> Bool
$c> :: PDescription -> PDescription -> Bool
<= :: PDescription -> PDescription -> Bool
$c<= :: PDescription -> PDescription -> Bool
< :: PDescription -> PDescription -> Bool
$c< :: PDescription -> PDescription -> Bool
compare :: PDescription -> PDescription -> Ordering
$ccompare :: PDescription -> PDescription -> Ordering
$cp1Ord :: Eq PDescription
Ord)

-- | The weight determines how often a 'Proposal' is executed per iteration of
-- the Markov chain.
newtype PWeight = PWeight {PWeight -> Int
fromPWeight :: Int}
  deriving (Int -> PWeight -> ShowS
[PWeight] -> ShowS
PWeight -> String
(Int -> PWeight -> ShowS)
-> (PWeight -> String) -> ([PWeight] -> ShowS) -> Show PWeight
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [PWeight] -> ShowS
$cshowList :: [PWeight] -> ShowS
show :: PWeight -> String
$cshow :: PWeight -> String
showsPrec :: Int -> PWeight -> ShowS
$cshowsPrec :: Int -> PWeight -> ShowS
Show, PWeight -> PWeight -> Bool
(PWeight -> PWeight -> Bool)
-> (PWeight -> PWeight -> Bool) -> Eq PWeight
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: PWeight -> PWeight -> Bool
$c/= :: PWeight -> PWeight -> Bool
== :: PWeight -> PWeight -> Bool
$c== :: PWeight -> PWeight -> Bool
Eq, Eq PWeight
Eq PWeight
-> (PWeight -> PWeight -> Ordering)
-> (PWeight -> PWeight -> Bool)
-> (PWeight -> PWeight -> Bool)
-> (PWeight -> PWeight -> Bool)
-> (PWeight -> PWeight -> Bool)
-> (PWeight -> PWeight -> PWeight)
-> (PWeight -> PWeight -> PWeight)
-> Ord PWeight
PWeight -> PWeight -> Bool
PWeight -> PWeight -> Ordering
PWeight -> PWeight -> PWeight
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: PWeight -> PWeight -> PWeight
$cmin :: PWeight -> PWeight -> PWeight
max :: PWeight -> PWeight -> PWeight
$cmax :: PWeight -> PWeight -> PWeight
>= :: PWeight -> PWeight -> Bool
$c>= :: PWeight -> PWeight -> Bool
> :: PWeight -> PWeight -> Bool
$c> :: PWeight -> PWeight -> Bool
<= :: PWeight -> PWeight -> Bool
$c<= :: PWeight -> PWeight -> Bool
< :: PWeight -> PWeight -> Bool
$c< :: PWeight -> PWeight -> Bool
compare :: PWeight -> PWeight -> Ordering
$ccompare :: PWeight -> PWeight -> Ordering
$cp1Ord :: Eq PWeight
Ord)

-- | A 'Proposal' is an instruction about how the Markov chain will traverse the
-- state space @a@. Essentially, it is a probability mass or probability density
-- conditioned on the current state (i.e., a kernel).
--
-- A 'Proposal' may be tuneable in that it contains information about how to enlarge
-- or shrink the step size to tune the acceptance ratio.
--
-- No proposals with the same name and description are allowed in a 'Cycle'.
data Proposal a = Proposal
  { -- | Name of the affected variable.
    Proposal a -> PName
pName :: PName,
    -- | Description of the proposal type and parameters.
    Proposal a -> PDescription
pDescription :: PDescription,
    -- | The weight determines how often a 'Proposal' is executed per iteration of
    -- the Markov chain.
    Proposal a -> PWeight
pWeight :: PWeight,
    -- | Simple proposal without name, weight, and tuning information.
    Proposal a -> ProposalSimple a
pSimple :: ProposalSimple a,
    -- | Tuning is disabled if set to 'Nothing'.
    Proposal a -> Maybe (Tuner a)
pTuner :: Maybe (Tuner a)
  }

-- XXX: This should be removed.
instance Show (Proposal a) where
  show :: Proposal a -> String
show Proposal a
m = PName -> String
fromPName (Proposal a -> PName
forall a. Proposal a -> PName
pName Proposal a
m) String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
" " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> PDescription -> String
fromPDescription (Proposal a -> PDescription
forall a. Proposal a -> PDescription
pDescription Proposal a
m) String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
", weight " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Int -> String
forall a. Show a => a -> String
show (PWeight -> Int
fromPWeight (PWeight -> Int) -> PWeight -> Int
forall a b. (a -> b) -> a -> b
$ Proposal a -> PWeight
forall a. Proposal a -> PWeight
pWeight Proposal a
m)

instance Eq (Proposal a) where
  Proposal a
m == :: Proposal a -> Proposal a -> Bool
== Proposal a
n = Proposal a -> PName
forall a. Proposal a -> PName
pName Proposal a
m PName -> PName -> Bool
forall a. Eq a => a -> a -> Bool
== Proposal a -> PName
forall a. Proposal a -> PName
pName Proposal a
n Bool -> Bool -> Bool
&& Proposal a -> PDescription
forall a. Proposal a -> PDescription
pDescription Proposal a
m PDescription -> PDescription -> Bool
forall a. Eq a => a -> a -> Bool
== Proposal a -> PDescription
forall a. Proposal a -> PDescription
pDescription Proposal a
n

instance Ord (Proposal a) where
  compare :: Proposal a -> Proposal a -> Ordering
compare = (PDescription, PName, PWeight)
-> (PDescription, PName, PWeight) -> Ordering
forall a. Ord a => a -> a -> Ordering
compare ((PDescription, PName, PWeight)
 -> (PDescription, PName, PWeight) -> Ordering)
-> (Proposal a -> (PDescription, PName, PWeight))
-> Proposal a
-> Proposal a
-> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` (\Proposal a
p -> (Proposal a -> PDescription
forall a. Proposal a -> PDescription
pDescription Proposal a
p, Proposal a -> PName
forall a. Proposal a -> PName
pName Proposal a
p, Proposal a -> PWeight
forall a. Proposal a -> PWeight
pWeight Proposal a
p))

-- | Convert a proposal from one data type to another using a lens.
--
-- For example:
--
-- @
-- scaleFirstEntryOfTuple = _1 @~ scale
-- @
(@~) :: Lens' b a -> Proposal a -> Proposal b
@~ :: Lens' b a -> Proposal a -> Proposal b
(@~) Lens' b a
l (Proposal PName
n PDescription
d PWeight
w ProposalSimple a
s Maybe (Tuner a)
t) = PName
-> PDescription
-> PWeight
-> ProposalSimple b
-> Maybe (Tuner b)
-> Proposal b
forall a.
PName
-> PDescription
-> PWeight
-> ProposalSimple a
-> Maybe (Tuner a)
-> Proposal a
Proposal PName
n PDescription
d PWeight
w (Lens' b a -> ProposalSimple a -> ProposalSimple b
forall b a. Lens' b a -> ProposalSimple a -> ProposalSimple b
convertS Lens' b a
l ProposalSimple a
s) (Lens' b a -> Tuner a -> Tuner b
forall b a. Lens' b a -> Tuner a -> Tuner b
convertT Lens' b a
l (Tuner a -> Tuner b) -> Maybe (Tuner a) -> Maybe (Tuner b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Maybe (Tuner a)
t)

-- | Simple proposal without tuning information.
--
-- Instruction about randomly moving from the current state to a new state,
-- given some source of randomness.
--
-- In order to calculate the Metropolis-Hastings ratio, we need to know the
-- ratio of the backward to forward kernels (i.e., the probability masses or
-- probability densities) and the absolute value of the determinant of the
-- Jacobian matrix.
--
-- For unbiased proposals, these values are 1.0 such that
--
-- @
-- proposalSimpleUnbiased x g = return (x', 1.0, 1.0)
-- @
--
-- For biased proposals, the kernel ratio is qYX / qXY, where qXY is the
-- probability density to move from X to Y, and the absolute value of the
-- determinant of the Jacobian matrix differs from 1.0.
type ProposalSimple a = a -> GenIO -> IO (a, Log Double, Log Double)

convertS :: Lens' b a -> ProposalSimple a -> ProposalSimple b
convertS :: Lens' b a -> ProposalSimple a -> ProposalSimple b
convertS Lens' b a
l ProposalSimple a
s = b -> Gen RealWorld -> IO (b, Log Double, Log Double)
ProposalSimple b
s'
  where
    s' :: b -> Gen RealWorld -> IO (b, Log Double, Log Double)
s' b
v Gen RealWorld
g = do
      (a
x', Log Double
r, Log Double
j) <- ProposalSimple a
s (b
v b -> Getting a b a -> a
forall s a. s -> Getting a s a -> a
^. Getting a b a
Lens' b a
l) Gen RealWorld
Gen (PrimState IO)
g
      (b, Log Double, Log Double) -> IO (b, Log Double, Log Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (ASetter b b a a -> a -> b -> b
forall s t a b. ASetter s t a b -> b -> s -> t
set ASetter b b a a
Lens' b a
l a
x' b
v, Log Double
r, Log Double
j)

-- | Tune the acceptance ratio of a 'Proposal'; see 'tune', or 'autotuneCycle'.
data Tuner a = Tuner
  { Tuner a -> Double
tParam :: Double,
    Tuner a -> Double -> ProposalSimple a
tFunc :: Double -> ProposalSimple a
  }

convertT :: Lens' b a -> Tuner a -> Tuner b
convertT :: Lens' b a -> Tuner a -> Tuner b
convertT Lens' b a
l (Tuner Double
p Double -> ProposalSimple a
f) = Double -> (Double -> ProposalSimple b) -> Tuner b
forall a. Double -> (Double -> ProposalSimple a) -> Tuner a
Tuner Double
p Double -> b -> Gen RealWorld -> IO (b, Log Double, Log Double)
Double -> ProposalSimple b
f'
  where
    f' :: Double -> ProposalSimple b
f' Double
x = Lens' b a -> ProposalSimple a -> ProposalSimple b
forall b a. Lens' b a -> ProposalSimple a -> ProposalSimple b
convertS Lens' b a
l (ProposalSimple a -> ProposalSimple b)
-> ProposalSimple a -> ProposalSimple b
forall a b. (a -> b) -> a -> b
$ Double -> ProposalSimple a
f Double
x

-- | Tune the proposal?
data Tune = Tune | NoTune
  deriving (Int -> Tune -> ShowS
[Tune] -> ShowS
Tune -> String
(Int -> Tune -> ShowS)
-> (Tune -> String) -> ([Tune] -> ShowS) -> Show Tune
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Tune] -> ShowS
$cshowList :: [Tune] -> ShowS
show :: Tune -> String
$cshow :: Tune -> String
showsPrec :: Int -> Tune -> ShowS
$cshowsPrec :: Int -> Tune -> ShowS
Show, Tune -> Tune -> Bool
(Tune -> Tune -> Bool) -> (Tune -> Tune -> Bool) -> Eq Tune
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Tune -> Tune -> Bool
$c/= :: Tune -> Tune -> Bool
== :: Tune -> Tune -> Bool
$c== :: Tune -> Tune -> Bool
Eq)

-- | Create a possibly tuneable proposal.
createProposal ::
  -- | Description of the proposal type and parameters.
  PDescription ->
  -- | Function creating a simple proposal for a given tuning parameter. The
  -- larger the tuning parameter, the larger the proposal (and the lower the
  -- expected acceptance ratio), and vice versa.
  (Double -> ProposalSimple a) ->
  -- | Name.
  PName ->
  -- | PWeight.
  PWeight ->
  -- | Activate tuning?
  Tune ->
  Proposal a
createProposal :: PDescription
-> (Double -> ProposalSimple a)
-> PName
-> PWeight
-> Tune
-> Proposal a
createProposal PDescription
d Double -> ProposalSimple a
f PName
n PWeight
w Tune
Tune = PName
-> PDescription
-> PWeight
-> ProposalSimple a
-> Maybe (Tuner a)
-> Proposal a
forall a.
PName
-> PDescription
-> PWeight
-> ProposalSimple a
-> Maybe (Tuner a)
-> Proposal a
Proposal PName
n PDescription
d PWeight
w (Double -> ProposalSimple a
f Double
1.0) (Tuner a -> Maybe (Tuner a)
forall a. a -> Maybe a
Just (Tuner a -> Maybe (Tuner a)) -> Tuner a -> Maybe (Tuner a)
forall a b. (a -> b) -> a -> b
$ Double -> (Double -> ProposalSimple a) -> Tuner a
forall a. Double -> (Double -> ProposalSimple a) -> Tuner a
Tuner Double
1.0 Double -> ProposalSimple a
f)
createProposal PDescription
d Double -> ProposalSimple a
f PName
n PWeight
w Tune
NoTune = PName
-> PDescription
-> PWeight
-> ProposalSimple a
-> Maybe (Tuner a)
-> Proposal a
forall a.
PName
-> PDescription
-> PWeight
-> ProposalSimple a
-> Maybe (Tuner a)
-> Proposal a
Proposal PName
n PDescription
d PWeight
w (Double -> ProposalSimple a
f Double
1.0) Maybe (Tuner a)
forall a. Maybe a
Nothing

-- Minimal tuning parameter; subject to change.
tuningParamMin :: Double
tuningParamMin :: Double
tuningParamMin = Double
1e-12

-- | Tune a 'Proposal'. Return 'Nothing' if 'Proposal' is not tuneable. If the parameter
--   @dt@ is larger than 1.0, the 'Proposal' is enlarged, if @0<dt<1.0@, it is
--   shrunk. Negative tuning parameters are not allowed.
tune :: Double -> Proposal a -> Maybe (Proposal a)
tune :: Double -> Proposal a -> Maybe (Proposal a)
tune Double
dt Proposal a
m
  | Double
dt Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 = String -> Maybe (Proposal a)
forall a. HasCallStack => String -> a
error (String -> Maybe (Proposal a)) -> String -> Maybe (Proposal a)
forall a b. (a -> b) -> a -> b
$ String
"tune: Tuning parameter not positive: " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Double -> String
forall a. Show a => a -> String
show Double
dt String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
"."
  | Bool
otherwise = do
    (Tuner Double
t Double -> ProposalSimple a
f) <- Proposal a -> Maybe (Tuner a)
forall a. Proposal a -> Maybe (Tuner a)
pTuner Proposal a
m
    -- Ensure that the tuning parameter is not too small.
    let t' :: Double
t' = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
tuningParamMin (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dt)
    Proposal a -> Maybe (Proposal a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Proposal a -> Maybe (Proposal a))
-> Proposal a -> Maybe (Proposal a)
forall a b. (a -> b) -> a -> b
$ Proposal a
m {pSimple :: ProposalSimple a
pSimple = Double -> ProposalSimple a
f Double
t', pTuner :: Maybe (Tuner a)
pTuner = Tuner a -> Maybe (Tuner a)
forall a. a -> Maybe a
Just (Tuner a -> Maybe (Tuner a)) -> Tuner a -> Maybe (Tuner a)
forall a b. (a -> b) -> a -> b
$ Double -> (Double -> ProposalSimple a) -> Tuner a
forall a. Double -> (Double -> ProposalSimple a) -> Tuner a
Tuner Double
t' Double -> ProposalSimple a
f}

-- The desired acceptance ratio 0.44 is optimal for one-dimensional proposals;
-- one could also store the affected number of dimensions with the proposal and
-- tune towards an acceptance ratio accounting for the number of dimensions.
--
-- The optimal ratios seem to be:
-- - One dimension: 0.44 (numerical result).
-- - Five and more dimensions: 0.234 seems to be a good value (numerical result).
-- - Infinite dimensions: 0.234 (theorem for specific target distributions).
-- See Handbook of Markov chain Monte Carlo, chapter 4.
ratioOpt :: Double
ratioOpt :: Double
ratioOpt = Double
0.44

-- Warn if acceptance ratio is lower.
ratioMin :: Double
ratioMin :: Double
ratioMin = Double
0.1

-- Warn if acceptance ratio is larger.
ratioMax :: Double
ratioMax :: Double
ratioMax = Double
0.9

-- | Define the order in which 'Proposal's are executed in a 'Cycle'. The total
-- number of 'Proposal's per 'Cycle' may differ between 'Order's (e.g., compare
-- 'RandomO' and 'RandomReversibleO').
data Order
  = -- | Shuffle the 'Proposal's in the 'Cycle'. The 'Proposal's are replicated
    -- according to their weights and executed in random order. If a 'Proposal' has
    -- weight @w@, it is executed exactly @w@ times per iteration.
    RandomO
  | -- | The 'Proposal's are executed sequentially, in the order they appear in the
    -- 'Cycle'. 'Proposal's with weight @w>1@ are repeated immediately @w@ times
    -- (and not appended to the end of the list).
    SequentialO
  | -- | Similar to 'RandomO'. However, a reversed copy of the list of
    --  shuffled 'Proposal's is appended such that the resulting Markov chain is
    --  reversible.
    --  Note: the total number of 'Proposal's executed per cycle is twice the number
    --  of 'RandomO'.
    RandomReversibleO
  | -- | Similar to 'SequentialO'. However, a reversed copy of the list of
    -- sequentially ordered 'Proposal's is appended such that the resulting Markov
    -- chain is reversible.
    SequentialReversibleO
  deriving (Order -> Order -> Bool
(Order -> Order -> Bool) -> (Order -> Order -> Bool) -> Eq Order
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Order -> Order -> Bool
$c/= :: Order -> Order -> Bool
== :: Order -> Order -> Bool
$c== :: Order -> Order -> Bool
Eq, Int -> Order -> ShowS
[Order] -> ShowS
Order -> String
(Int -> Order -> ShowS)
-> (Order -> String) -> ([Order] -> ShowS) -> Show Order
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Order] -> ShowS
$cshowList :: [Order] -> ShowS
show :: Order -> String
$cshow :: Order -> String
showsPrec :: Int -> Order -> ShowS
$cshowsPrec :: Int -> Order -> ShowS
Show)

instance Default Order where def :: Order
def = Order
RandomO

-- | In brief, a 'Cycle' is a list of proposals.
--
-- The state of the Markov chain will be logged only after all 'Proposal's in
-- the 'Cycle' have been completed, and the iteration counter will be increased
-- by one. The order in which the 'Proposal's are executed is specified by
-- 'Order'. The default is 'RandomO'.
--
-- __Proposals must have unique names__, so that they can be identified.
data Cycle a = Cycle
  { Cycle a -> [Proposal a]
ccProposals :: [Proposal a],
    Cycle a -> Order
ccOrder :: Order
  }

-- | Create a 'Cycle' from a list of 'Proposal's.
fromList :: [Proposal a] -> Cycle a
fromList :: [Proposal a] -> Cycle a
fromList [] =
  String -> Cycle a
forall a. HasCallStack => String -> a
error String
"fromList: Received an empty list but cannot create an empty Cycle."
fromList [Proposal a]
xs =
  if [Proposal a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Proposal a] -> [Proposal a]
forall a. Eq a => [a] -> [a]
nub [Proposal a]
xs) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Proposal a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Proposal a]
xs
    then [Proposal a] -> Order -> Cycle a
forall a. [Proposal a] -> Order -> Cycle a
Cycle [Proposal a]
xs Order
forall a. Default a => a
def
    else String -> Cycle a
forall a. HasCallStack => String -> a
error String
"fromList: Proposals are not unique."

-- | Set the order of 'Proposal's in a 'Cycle'.
setOrder :: Order -> Cycle a -> Cycle a
setOrder :: Order -> Cycle a -> Cycle a
setOrder Order
o Cycle a
c = Cycle a
c {ccOrder :: Order
ccOrder = Order
o}

-- | Replicate 'Proposal's according to their weights and possibly shuffle them.
getNIterations :: Cycle a -> Int -> GenIO -> IO [[Proposal a]]
getNIterations :: Cycle a -> Int -> Gen (PrimState IO) -> IO [[Proposal a]]
getNIterations (Cycle [Proposal a]
xs Order
o) Int
n Gen (PrimState IO)
g = case Order
o of
  Order
RandomO -> [Proposal a] -> Int -> Gen (PrimState IO) -> IO [[Proposal a]]
forall a. [a] -> Int -> Gen (PrimState IO) -> IO [[a]]
shuffleN [Proposal a]
ps Int
n Gen (PrimState IO)
g
  Order
SequentialO -> [[Proposal a]] -> IO [[Proposal a]]
forall (m :: * -> *) a. Monad m => a -> m a
return ([[Proposal a]] -> IO [[Proposal a]])
-> [[Proposal a]] -> IO [[Proposal a]]
forall a b. (a -> b) -> a -> b
$ Int -> [Proposal a] -> [[Proposal a]]
forall a. Int -> a -> [a]
replicate Int
n [Proposal a]
ps
  Order
RandomReversibleO -> do
    [[Proposal a]]
psRs <- [Proposal a] -> Int -> Gen (PrimState IO) -> IO [[Proposal a]]
forall a. [a] -> Int -> Gen (PrimState IO) -> IO [[a]]
shuffleN [Proposal a]
ps Int
n Gen (PrimState IO)
g
    [[Proposal a]] -> IO [[Proposal a]]
forall (m :: * -> *) a. Monad m => a -> m a
return [[Proposal a]
psR [Proposal a] -> [Proposal a] -> [Proposal a]
forall a. [a] -> [a] -> [a]
++ [Proposal a] -> [Proposal a]
forall a. [a] -> [a]
reverse [Proposal a]
psR | [Proposal a]
psR <- [[Proposal a]]
psRs]
  Order
SequentialReversibleO -> [[Proposal a]] -> IO [[Proposal a]]
forall (m :: * -> *) a. Monad m => a -> m a
return ([[Proposal a]] -> IO [[Proposal a]])
-> [[Proposal a]] -> IO [[Proposal a]]
forall a b. (a -> b) -> a -> b
$ Int -> [Proposal a] -> [[Proposal a]]
forall a. Int -> a -> [a]
replicate Int
n ([Proposal a] -> [[Proposal a]]) -> [Proposal a] -> [[Proposal a]]
forall a b. (a -> b) -> a -> b
$ [Proposal a]
ps [Proposal a] -> [Proposal a] -> [Proposal a]
forall a. [a] -> [a] -> [a]
++ [Proposal a] -> [Proposal a]
forall a. [a] -> [a]
reverse [Proposal a]
ps
  where
    !ps :: [Proposal a]
ps = [[Proposal a]] -> [Proposal a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [Int -> Proposal a -> [Proposal a]
forall a. Int -> a -> [a]
replicate (PWeight -> Int
fromPWeight (PWeight -> Int) -> PWeight -> Int
forall a b. (a -> b) -> a -> b
$ Proposal a -> PWeight
forall a. Proposal a -> PWeight
pWeight Proposal a
m) Proposal a
m | Proposal a
m <- [Proposal a]
xs]

-- | Tune 'Proposal's in the 'Cycle'. See 'tune'.
tuneCycle :: Map (Proposal a) Double -> Cycle a -> Cycle a
tuneCycle :: Map (Proposal a) Double -> Cycle a -> Cycle a
tuneCycle Map (Proposal a) Double
m Cycle a
c =
  if [Proposal a] -> [Proposal a]
forall a. Ord a => [a] -> [a]
sort (Map (Proposal a) Double -> [Proposal a]
forall k a. Map k a -> [k]
M.keys Map (Proposal a) Double
m) [Proposal a] -> [Proposal a] -> Bool
forall a. Eq a => a -> a -> Bool
== [Proposal a] -> [Proposal a]
forall a. Ord a => [a] -> [a]
sort [Proposal a]
ps
    then Cycle a
c {ccProposals :: [Proposal a]
ccProposals = (Proposal a -> Proposal a) -> [Proposal a] -> [Proposal a]
forall a b. (a -> b) -> [a] -> [b]
map Proposal a -> Proposal a
tuneF [Proposal a]
ps}
    else String -> Cycle a
forall a. HasCallStack => String -> a
error String
"tuneCycle: Map contains proposals that are not in the cycle."
  where
    ps :: [Proposal a]
ps = Cycle a -> [Proposal a]
forall a. Cycle a -> [Proposal a]
ccProposals Cycle a
c
    tuneF :: Proposal a -> Proposal a
tuneF Proposal a
p = case Map (Proposal a) Double
m Map (Proposal a) Double -> Proposal a -> Maybe Double
forall k a. Ord k => Map k a -> k -> Maybe a
M.!? Proposal a
p of
      Maybe Double
Nothing -> Proposal a
p
      Just Double
x -> Proposal a -> Maybe (Proposal a) -> Proposal a
forall a. a -> Maybe a -> a
fromMaybe Proposal a
p (Double -> Proposal a -> Maybe (Proposal a)
forall a. Double -> Proposal a -> Maybe (Proposal a)
tune Double
x Proposal a
p)

-- | Calculate acceptance ratios and auto tune the 'Proposal's in the 'Cycle'. For
-- now, a 'Proposal' is enlarged when the acceptance ratio is above 0.44, and
-- shrunk otherwise. Do not change 'Proposal's that are not tuneable.
autotuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a
autotuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a
autotuneCycle Acceptance (Proposal a)
a = Map (Proposal a) Double -> Cycle a -> Cycle a
forall a. Map (Proposal a) Double -> Cycle a -> Cycle a
tuneCycle ((Double -> Double)
-> Map (Proposal a) Double -> Map (Proposal a) Double
forall a b k. (a -> b) -> Map k a -> Map k b
M.map (\Double
x -> Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ratioOpt) (Map (Proposal a) Double -> Map (Proposal a) Double)
-> Map (Proposal a) Double -> Map (Proposal a) Double
forall a b. (a -> b) -> a -> b
$ Acceptance (Proposal a) -> Map (Proposal a) Double
forall k. Acceptance k -> Map k Double
acceptanceRatios Acceptance (Proposal a)
a)

renderRow ::
  BL.ByteString ->
  BL.ByteString ->
  BL.ByteString ->
  BL.ByteString ->
  BL.ByteString ->
  BL.ByteString ->
  BL.ByteString ->
  BL.ByteString ->
  BL.ByteString
renderRow :: ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
renderRow ByteString
name ByteString
ptype ByteString
weight ByteString
nAccept ByteString
nReject ByteString
acceptRatio ByteString
tuneParam ByteString
manualAdjustment = ByteString
"   " ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
nm ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
pt ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
wt ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
na ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
nr ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
ra ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
tp ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
mt
  where
    nm :: ByteString
nm = Int -> ByteString -> ByteString
alignLeft Int
30 ByteString
name
    pt :: ByteString
pt = Int -> ByteString -> ByteString
alignLeft Int
50 ByteString
ptype
    wt :: ByteString
wt = Int -> ByteString -> ByteString
alignRight Int
8 ByteString
weight
    na :: ByteString
na = Int -> ByteString -> ByteString
alignRight Int
15 ByteString
nAccept
    nr :: ByteString
nr = Int -> ByteString -> ByteString
alignRight Int
15 ByteString
nReject
    ra :: ByteString
ra = Int -> ByteString -> ByteString
alignRight Int
15 ByteString
acceptRatio
    tp :: ByteString
tp = Int -> ByteString -> ByteString
alignRight Int
20 ByteString
tuneParam
    mt :: ByteString
mt = Int -> ByteString -> ByteString
alignRight Int
30 ByteString
manualAdjustment

proposalHeader :: BL.ByteString
proposalHeader :: ByteString
proposalHeader =
  ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
renderRow ByteString
"Name" ByteString
"Description" ByteString
"Weight" ByteString
"Accepted" ByteString
"Rejected" ByteString
"Ratio" ByteString
"Tuning parameter" ByteString
"Consider manual adjustment"

summarizeProposal :: Proposal a -> Maybe (Int, Int, Double) -> BL.ByteString
summarizeProposal :: Proposal a -> Maybe (Int, Int, Double) -> ByteString
summarizeProposal Proposal a
m Maybe (Int, Int, Double)
r =
  ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
renderRow
    (String -> ByteString
BL.pack (String -> ByteString) -> String -> ByteString
forall a b. (a -> b) -> a -> b
$ PName -> String
fromPName (PName -> String) -> PName -> String
forall a b. (a -> b) -> a -> b
$ Proposal a -> PName
forall a. Proposal a -> PName
pName Proposal a
m)
    (String -> ByteString
BL.pack (String -> ByteString) -> String -> ByteString
forall a b. (a -> b) -> a -> b
$ PDescription -> String
fromPDescription (PDescription -> String) -> PDescription -> String
forall a b. (a -> b) -> a -> b
$ Proposal a -> PDescription
forall a. Proposal a -> PDescription
pDescription Proposal a
m)
    ByteString
weight
    ByteString
nAccept
    ByteString
nReject
    ByteString
acceptRatio
    ByteString
tuneParamStr
    ByteString
manualAdjustmentStr
  where
    weight :: ByteString
weight = Builder -> ByteString
BB.toLazyByteString (Builder -> ByteString) -> Builder -> ByteString
forall a b. (a -> b) -> a -> b
$ Int -> Builder
BB.intDec (Int -> Builder) -> Int -> Builder
forall a b. (a -> b) -> a -> b
$ PWeight -> Int
fromPWeight (PWeight -> Int) -> PWeight -> Int
forall a b. (a -> b) -> a -> b
$ Proposal a -> PWeight
forall a. Proposal a -> PWeight
pWeight Proposal a
m
    nAccept :: ByteString
nAccept = Builder -> ByteString
BB.toLazyByteString (Builder -> ByteString) -> Builder -> ByteString
forall a b. (a -> b) -> a -> b
$ Builder
-> ((Int, Int, Double) -> Builder)
-> Maybe (Int, Int, Double)
-> Builder
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Builder
"" (Int -> Builder
BB.intDec (Int -> Builder)
-> ((Int, Int, Double) -> Int) -> (Int, Int, Double) -> Builder
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int, Double) -> Getting Int (Int, Int, Double) Int -> Int
forall s a. s -> Getting a s a -> a
^. Getting Int (Int, Int, Double) Int
forall s t a b. Field1 s t a b => Lens s t a b
_1)) Maybe (Int, Int, Double)
r
    nReject :: ByteString
nReject = Builder -> ByteString
BB.toLazyByteString (Builder -> ByteString) -> Builder -> ByteString
forall a b. (a -> b) -> a -> b
$ Builder
-> ((Int, Int, Double) -> Builder)
-> Maybe (Int, Int, Double)
-> Builder
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Builder
"" (Int -> Builder
BB.intDec (Int -> Builder)
-> ((Int, Int, Double) -> Int) -> (Int, Int, Double) -> Builder
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int, Double) -> Getting Int (Int, Int, Double) Int -> Int
forall s a. s -> Getting a s a -> a
^. Getting Int (Int, Int, Double) Int
forall s t a b. Field2 s t a b => Lens s t a b
_2)) Maybe (Int, Int, Double)
r
    acceptRatio :: ByteString
acceptRatio = ByteString -> ByteString
BL.fromStrict (ByteString -> ByteString) -> ByteString -> ByteString
forall a b. (a -> b) -> a -> b
$ ByteString
-> ((Int, Int, Double) -> ByteString)
-> Maybe (Int, Int, Double)
-> ByteString
forall b a. b -> (a -> b) -> Maybe a -> b
maybe ByteString
"" (Int -> Double -> ByteString
BC.toFixed Int
3 (Double -> ByteString)
-> ((Int, Int, Double) -> Double)
-> (Int, Int, Double)
-> ByteString
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int, Double)
-> Getting Double (Int, Int, Double) Double -> Double
forall s a. s -> Getting a s a -> a
^. Getting Double (Int, Int, Double) Double
forall s t a b. Field3 s t a b => Lens s t a b
_3)) Maybe (Int, Int, Double)
r
    tuneParamStr :: ByteString
tuneParamStr = ByteString -> ByteString
BL.fromStrict (ByteString -> ByteString) -> ByteString -> ByteString
forall a b. (a -> b) -> a -> b
$ ByteString -> (Double -> ByteString) -> Maybe Double -> ByteString
forall b a. b -> (a -> b) -> Maybe a -> b
maybe ByteString
"" (Int -> Double -> ByteString
BC.toFixed Int
3) (Tuner a -> Double
forall a. Tuner a -> Double
tParam (Tuner a -> Double) -> Maybe (Tuner a) -> Maybe Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Proposal a -> Maybe (Tuner a)
forall a. Proposal a -> Maybe (Tuner a)
pTuner Proposal a
m)
    check :: Double -> p
check Double
v
      | Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
ratioMin = p
"ratio too low"
      | Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
ratioMax = p
"ratio too high"
      | Bool
otherwise = p
""
    manualAdjustmentStr :: ByteString
manualAdjustmentStr = ByteString -> ByteString
BL.fromStrict (ByteString -> ByteString) -> ByteString -> ByteString
forall a b. (a -> b) -> a -> b
$ ByteString
-> ((Int, Int, Double) -> ByteString)
-> Maybe (Int, Int, Double)
-> ByteString
forall b a. b -> (a -> b) -> Maybe a -> b
maybe ByteString
"" (Double -> ByteString
forall p. IsString p => Double -> p
check (Double -> ByteString)
-> ((Int, Int, Double) -> Double)
-> (Int, Int, Double)
-> ByteString
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int, Double)
-> Getting Double (Int, Int, Double) Double -> Double
forall s a. s -> Getting a s a -> a
^. Getting Double (Int, Int, Double) Double
forall s t a b. Field3 s t a b => Lens s t a b
_3)) Maybe (Int, Int, Double)
r

hLine :: BL.ByteString -> BL.ByteString
hLine :: ByteString -> ByteString
hLine ByteString
s = ByteString
"   " ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> Int64 -> Char -> ByteString
BL.replicate (ByteString -> Int64
BL.length ByteString
s Int64 -> Int64 -> Int64
forall a. Num a => a -> a -> a
- Int64
3) Char
'-'

-- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance ratios.
summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> BL.ByteString
summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> ByteString
summarizeCycle Acceptance (Proposal a)
a Cycle a
c =
  ByteString -> [ByteString] -> ByteString
BL.intercalate ByteString
"\n" ([ByteString] -> ByteString) -> [ByteString] -> ByteString
forall a b. (a -> b) -> a -> b
$
    [ ByteString
"Summary of proposal(s) in cycle. " ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
mpi ByteString -> ByteString -> ByteString
forall a. Semigroup a => a -> a -> a
<> ByteString
" proposal(s) per iteration.",
      ByteString
proposalHeader,
      ByteString -> ByteString
hLine ByteString
proposalHeader
    ]
      [ByteString] -> [ByteString] -> [ByteString]
forall a. [a] -> [a] -> [a]
++ [Proposal a -> Maybe (Int, Int, Double) -> ByteString
forall a. Proposal a -> Maybe (Int, Int, Double) -> ByteString
summarizeProposal Proposal a
m (Proposal a -> Maybe (Int, Int, Double)
ar Proposal a
m) | Proposal a
m <- [Proposal a]
ps]
      [ByteString] -> [ByteString] -> [ByteString]
forall a. [a] -> [a] -> [a]
++ [ByteString -> ByteString
hLine ByteString
proposalHeader]
  where
    ps :: [Proposal a]
ps = Cycle a -> [Proposal a]
forall a. Cycle a -> [Proposal a]
ccProposals Cycle a
c
    mpi :: ByteString
mpi = Builder -> ByteString
BB.toLazyByteString (Builder -> ByteString) -> Builder -> ByteString
forall a b. (a -> b) -> a -> b
$ Int -> Builder
BB.intDec (Int -> Builder) -> Int -> Builder
forall a b. (a -> b) -> a -> b
$ [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ (Proposal a -> Int) -> [Proposal a] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (PWeight -> Int
fromPWeight (PWeight -> Int) -> (Proposal a -> PWeight) -> Proposal a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Proposal a -> PWeight
forall a. Proposal a -> PWeight
pWeight) [Proposal a]
ps
    ar :: Proposal a -> Maybe (Int, Int, Double)
ar Proposal a
m = Proposal a -> Acceptance (Proposal a) -> Maybe (Int, Int, Double)
forall k.
(Show k, Ord k) =>
k -> Acceptance k -> Maybe (Int, Int, Double)
acceptanceRatio Proposal a
m Acceptance (Proposal a)
a

-- | For each key @k@, store the number of accepted and rejected proposals.
newtype Acceptance k = Acceptance {Acceptance k -> Map k (Int, Int)
fromAcceptance :: Map k (Int, Int)}

instance ToJSONKey k => ToJSON (Acceptance k) where
  toJSON :: Acceptance k -> Value
toJSON (Acceptance Map k (Int, Int)
m) = Map k (Int, Int) -> Value
forall a. ToJSON a => a -> Value
toJSON Map k (Int, Int)
m
  toEncoding :: Acceptance k -> Encoding
toEncoding (Acceptance Map k (Int, Int)
m) = Map k (Int, Int) -> Encoding
forall a. ToJSON a => a -> Encoding
toEncoding Map k (Int, Int)
m

instance (Ord k, FromJSONKey k) => FromJSON (Acceptance k) where
  parseJSON :: Value -> Parser (Acceptance k)
parseJSON Value
v = Map k (Int, Int) -> Acceptance k
forall k. Map k (Int, Int) -> Acceptance k
Acceptance (Map k (Int, Int) -> Acceptance k)
-> Parser (Map k (Int, Int)) -> Parser (Acceptance k)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Value -> Parser (Map k (Int, Int))
forall a. FromJSON a => Value -> Parser a
parseJSON Value
v

-- | In the beginning there was the Word.
--
-- Initialize an empty storage of accepted/rejected values.
emptyA :: Ord k => [k] -> Acceptance k
emptyA :: [k] -> Acceptance k
emptyA [k]
ks = Map k (Int, Int) -> Acceptance k
forall k. Map k (Int, Int) -> Acceptance k
Acceptance (Map k (Int, Int) -> Acceptance k)
-> Map k (Int, Int) -> Acceptance k
forall a b. (a -> b) -> a -> b
$ [(k, (Int, Int))] -> Map k (Int, Int)
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [(k
k, (Int
0, Int
0)) | k
k <- [k]
ks]

-- | For key @k@, prepend an accepted (True) or rejected (False) proposal.
pushA :: (Ord k, Show k) => k -> Bool -> Acceptance k -> Acceptance k
pushA :: k -> Bool -> Acceptance k -> Acceptance k
pushA k
k Bool
True = Map k (Int, Int) -> Acceptance k
forall k. Map k (Int, Int) -> Acceptance k
Acceptance (Map k (Int, Int) -> Acceptance k)
-> (Acceptance k -> Map k (Int, Int))
-> Acceptance k
-> Acceptance k
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int) -> (Int, Int))
-> k -> Map k (Int, Int) -> Map k (Int, Int)
forall k a. Ord k => (a -> a) -> k -> Map k a -> Map k a
M.adjust ((Int -> Int) -> (Int, Int) -> (Int, Int)
forall (p :: * -> * -> *) a b c.
Bifunctor p =>
(a -> b) -> p a c -> p b c
first Int -> Int
forall a. Enum a => a -> a
succ) k
k (Map k (Int, Int) -> Map k (Int, Int))
-> (Acceptance k -> Map k (Int, Int))
-> Acceptance k
-> Map k (Int, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Acceptance k -> Map k (Int, Int)
forall k. Acceptance k -> Map k (Int, Int)
fromAcceptance
pushA k
k Bool
False = Map k (Int, Int) -> Acceptance k
forall k. Map k (Int, Int) -> Acceptance k
Acceptance (Map k (Int, Int) -> Acceptance k)
-> (Acceptance k -> Map k (Int, Int))
-> Acceptance k
-> Acceptance k
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Int) -> (Int, Int))
-> k -> Map k (Int, Int) -> Map k (Int, Int)
forall k a. Ord k => (a -> a) -> k -> Map k a -> Map k a
M.adjust ((Int -> Int) -> (Int, Int) -> (Int, Int)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second Int -> Int
forall a. Enum a => a -> a
succ) k
k (Map k (Int, Int) -> Map k (Int, Int))
-> (Acceptance k -> Map k (Int, Int))
-> Acceptance k
-> Map k (Int, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Acceptance k -> Map k (Int, Int)
forall k. Acceptance k -> Map k (Int, Int)
fromAcceptance
{-# INLINEABLE pushA #-}

-- | Reset acceptance storage.
resetA :: Ord k => Acceptance k -> Acceptance k
resetA :: Acceptance k -> Acceptance k
resetA = [k] -> Acceptance k
forall k. Ord k => [k] -> Acceptance k
emptyA ([k] -> Acceptance k)
-> (Acceptance k -> [k]) -> Acceptance k -> Acceptance k
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Map k (Int, Int) -> [k]
forall k a. Map k a -> [k]
M.keys (Map k (Int, Int) -> [k])
-> (Acceptance k -> Map k (Int, Int)) -> Acceptance k -> [k]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Acceptance k -> Map k (Int, Int)
forall k. Acceptance k -> Map k (Int, Int)
fromAcceptance

transformKeys :: (Ord k1, Ord k2) => [k1] -> [k2] -> Map k1 v -> Map k2 v
transformKeys :: [k1] -> [k2] -> Map k1 v -> Map k2 v
transformKeys [k1]
ks1 [k2]
ks2 Map k1 v
m = (Map k2 v -> (k1, k2) -> Map k2 v)
-> Map k2 v -> [(k1, k2)] -> Map k2 v
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Map k2 v -> (k1, k2) -> Map k2 v
forall k. Ord k => Map k v -> (k1, k) -> Map k v
insrt Map k2 v
forall k a. Map k a
M.empty ([(k1, k2)] -> Map k2 v) -> [(k1, k2)] -> Map k2 v
forall a b. (a -> b) -> a -> b
$ [k1] -> [k2] -> [(k1, k2)]
forall a b. [a] -> [b] -> [(a, b)]
zip [k1]
ks1 [k2]
ks2
  where
    insrt :: Map k v -> (k1, k) -> Map k v
insrt Map k v
m' (k1
k1, k
k2) = k -> v -> Map k v -> Map k v
forall k a. Ord k => k -> a -> Map k a -> Map k a
M.insert k
k2 (Map k1 v
m Map k1 v -> k1 -> v
forall k a. Ord k => Map k a -> k -> a
M.! k1
k1) Map k v
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] -> Acceptance k1 -> Acceptance k2
transformKeysA :: [k1] -> [k2] -> Acceptance k1 -> Acceptance k2
transformKeysA [k1]
ks1 [k2]
ks2 = Map k2 (Int, Int) -> Acceptance k2
forall k. Map k (Int, Int) -> Acceptance k
Acceptance (Map k2 (Int, Int) -> Acceptance k2)
-> (Acceptance k1 -> Map k2 (Int, Int))
-> Acceptance k1
-> Acceptance k2
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [k1] -> [k2] -> Map k1 (Int, Int) -> Map k2 (Int, Int)
forall k1 k2 v.
(Ord k1, Ord k2) =>
[k1] -> [k2] -> Map k1 v -> Map k2 v
transformKeys [k1]
ks1 [k2]
ks2 (Map k1 (Int, Int) -> Map k2 (Int, Int))
-> (Acceptance k1 -> Map k1 (Int, Int))
-> Acceptance k1
-> Map k2 (Int, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Acceptance k1 -> Map k1 (Int, Int)
forall k. Acceptance k -> Map k (Int, Int)
fromAcceptance

-- | Acceptance counts and ratio for a specific proposal.
acceptanceRatio :: (Show k, Ord k) => k -> Acceptance k -> Maybe (Int, Int, Double)
acceptanceRatio :: k -> Acceptance k -> Maybe (Int, Int, Double)
acceptanceRatio k
k Acceptance k
a = case Acceptance k -> Map k (Int, Int)
forall k. Acceptance k -> Map k (Int, Int)
fromAcceptance Acceptance k
a Map k (Int, Int) -> k -> Maybe (Int, Int)
forall k a. Ord k => Map k a -> k -> Maybe a
M.!? k
k of
  Just (Int
0, Int
0) -> Maybe (Int, Int, Double)
forall a. Maybe a
Nothing
  Just (Int
as, Int
rs) -> (Int, Int, Double) -> Maybe (Int, Int, Double)
forall a. a -> Maybe a
Just (Int
as, Int
rs, Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
as Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
as Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
rs))
  Maybe (Int, Int)
Nothing -> String -> Maybe (Int, Int, Double)
forall a. HasCallStack => String -> a
error (String -> Maybe (Int, Int, Double))
-> String -> Maybe (Int, Int, Double)
forall a b. (a -> b) -> a -> b
$ String
"acceptanceRatio: Key not found in map: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ k -> String
forall a. Show a => a -> String
show k
k String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"."

-- | Acceptance ratios for all proposals.
acceptanceRatios :: Acceptance k -> Map k Double
acceptanceRatios :: Acceptance k -> Map k Double
acceptanceRatios = ((Int, Int) -> Double) -> Map k (Int, Int) -> Map k Double
forall a b k. (a -> b) -> Map k a -> Map k b
M.map (\(Int
as, Int
rs) -> Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
as Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
as Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
rs)) (Map k (Int, Int) -> Map k Double)
-> (Acceptance k -> Map k (Int, Int))
-> Acceptance k
-> Map k Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Acceptance k -> Map k (Int, Int)
forall k. Acceptance k -> Map k (Int, Int)
fromAcceptance