{-# LANGUAGE DerivingVia #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RankNTypes #-}

-- |
-- Module      :  Mcmc.Proposal
-- Description :  Proposals are instruction to move around the state space
-- Copyright   :  2021 Dominik Schrempf
-- 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
  ( -- * Proposals
    PName (..),
    PDescription (..),
    PWeight (fromPWeight),
    pWeight,
    PDimension (..),
    PSpeed (..),
    Proposal (..),
    KernelRatio,
    PResult (..),
    Jacobian,
    JacobianFunction,
    (@~),
    liftProposal,
    liftProposalWith,
    PFunction,
    createProposal,

    -- * Tuners
    Tuner (..),
    Tune (..),
    TuningParameter,
    TuningType (..),
    TuningFunction,
    AuxiliaryTuningParameters,
    tuningFunction,
    tuningFunctionWithAux,
    tuningFunctionOnlyAux,
    tuningParameterMin,
    tuningParameterMax,
    tuneWithTuningParameters,
    getOptimalRate,

    -- * Output
    proposalHeader,
    summarizeProposal,
  )
where

import Data.Bifunctor
import qualified Data.ByteString.Builder as BB
import qualified Data.ByteString.Lazy.Char8 as BL
import qualified Data.Double.Conversion.ByteString as BC
import Data.Function
import qualified Data.Vector as VB
import qualified Data.Vector.Unboxed as VU
import Lens.Micro
import Mcmc.Acceptance
import Mcmc.Internal.ByteString
import Mcmc.Jacobian
import Numeric.Log hiding (sum)
import System.Random.Stateful

-- | 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)
  deriving (Semigroup PName
PName
Semigroup PName
-> PName
-> (PName -> PName -> PName)
-> ([PName] -> PName)
-> Monoid PName
[PName] -> PName
PName -> PName -> PName
forall a.
Semigroup a -> a -> (a -> a -> a) -> ([a] -> a) -> Monoid a
mconcat :: [PName] -> PName
$cmconcat :: [PName] -> PName
mappend :: PName -> PName -> PName
$cmappend :: PName -> PName -> PName
mempty :: PName
$cmempty :: PName
$cp1Monoid :: Semigroup PName
Monoid, b -> PName -> PName
NonEmpty PName -> PName
PName -> PName -> PName
(PName -> PName -> PName)
-> (NonEmpty PName -> PName)
-> (forall b. Integral b => b -> PName -> PName)
-> Semigroup PName
forall b. Integral b => b -> PName -> PName
forall a.
(a -> a -> a)
-> (NonEmpty a -> a)
-> (forall b. Integral b => b -> a -> a)
-> Semigroup a
stimes :: b -> PName -> PName
$cstimes :: forall b. Integral b => b -> PName -> PName
sconcat :: NonEmpty PName -> PName
$csconcat :: NonEmpty PName -> PName
<> :: PName -> PName -> PName
$c<> :: PName -> PName -> PName
Semigroup) via String

-- | 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 positive weight determines how often a 'Proposal' is executed per
-- iteration of the Markov chain. Abstract data type; for construction, see
-- 'pWeight'.
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)

-- | Check if the weight is positive.
--
-- Call 'error' if weight is zero or negative.
pWeight :: Int -> PWeight
pWeight :: Int -> PWeight
pWeight Int
n
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = String -> PWeight
forall a. HasCallStack => String -> a
error (String -> PWeight) -> String -> PWeight
forall a b. (a -> b) -> a -> b
$ String
"pWeight: Proposal weight is zero or negative: " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Int -> String
forall a. Show a => a -> String
show Int
n String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
"."
  | Bool
otherwise = Int -> PWeight
PWeight Int
n

-- | Proposal dimension.
--
-- The number of affected, independent parameters.
--
-- The dimension is used to calculate the optimal acceptance rate, and does not
-- have to be exact.
--
-- Usually, the optimal acceptance rate of low dimensional proposals is higher
-- than for high dimensional ones. However, this is not always true (see below).
--
-- Further, optimal acceptance rates are still subject to controversies. To my
-- knowledge, research has focused on random walk proposals with multivariate
-- normal distributions of dimension @d@. In this case, the following acceptance
-- rates are desired:
--
-- - one dimension: 0.44 (numerical results);
--
-- - five and more dimensions: 0.234 (numerical results);
--
-- - infinite dimensions: 0.234 (theorem for specific target distributions).
--
-- See Handbook of Markov chain Monte Carlo, chapter 4.
--
-- Of course, many proposals may not be classical random walk proposals. For
-- example, the beta proposal on a simplex ('Mcmc.Proposal.Simplex.beta')
-- samples one new variable of the simplex from a beta distribution while
-- rescaling all other variables. What is the dimension of this proposal? Here,
-- the dimension is set to 2. The reason is that if the dimension of the simplex
-- is 2, two variables are changed. If the dimension of the simplex is high, one
-- variable is changed substantially, while all others are changed marginally.
--
-- Further, if a proposal changes a number of variables in the same way (and not
-- independently like in a random walk proposal), the dimension of the proposal
-- is set to the number of variables changed.
--
-- Moreover, proposals of unknown dimension are assumed to have high dimension,
-- and the optimal acceptance rate 0.234 is used.
--
-- Finally, special proposals may have completely different desired acceptance
-- rates. For example. the Hamiltonian Monte Carlo proposal (see
-- Mcmc.Proposal.Hamiltonian.hmc) has a desired acceptance rate of 0.65.
-- Specific acceptance rates can be set with 'PSpecial'.
data PDimension
  = PDimension Int
  | PDimensionUnknown
  | -- | Provide dimension ('Int') and desired acceptance rate ('Double').
    PSpecial Int Double

-- | Rough indication whether a proposal is fast or slow.
--
-- Useful during burn in. Slow proposals are not executed during fast auto
-- tuning periods.
--
-- See 'Mcmc.Settings.BurnInSettings'.
data PSpeed = PFast | PSlow
  deriving (PSpeed -> PSpeed -> Bool
(PSpeed -> PSpeed -> Bool)
-> (PSpeed -> PSpeed -> Bool) -> Eq PSpeed
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: PSpeed -> PSpeed -> Bool
$c/= :: PSpeed -> PSpeed -> Bool
== :: PSpeed -> PSpeed -> Bool
$c== :: PSpeed -> PSpeed -> Bool
Eq)

-- | 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 Markov kernel).
--
-- A 'Proposal' may be tuneable in that it contains information about how to
-- enlarge or shrink the proposal size to decrease or increase the acceptance
-- rate.
--
-- Predefined proposals are provided. To create custom proposals, one may use
-- the convenience function 'createProposal'.
data Proposal a = Proposal
  { -- | Name of the affected variable.
    Proposal a -> PName
prName :: PName,
    -- | Description of the proposal type and parameters.
    Proposal a -> PDescription
prDescription :: PDescription,
    Proposal a -> PSpeed
prSpeed :: PSpeed,
    Proposal a -> PDimension
prDimension :: PDimension,
    -- | The weight determines how often a 'Proposal' is executed per iteration of
    -- the Markov chain.
    Proposal a -> PWeight
prWeight :: PWeight,
    -- | Simple proposal function without name, weight, and tuning information.
    Proposal a -> PFunction a
prFunction :: PFunction a,
    -- | Tuning is disabled if set to 'Nothing'.
    Proposal a -> Maybe (Tuner a)
prTuner :: Maybe (Tuner a)
  }

instance Eq (Proposal a) where
  Proposal a
m == :: Proposal a -> Proposal a -> Bool
== Proposal a
n = Proposal a -> PName
forall a. Proposal a -> PName
prName Proposal a
m PName -> PName -> Bool
forall a. Eq a => a -> a -> Bool
== Proposal a -> PName
forall a. Proposal a -> PName
prName Proposal a
n Bool -> Bool -> Bool
&& Proposal a -> PDescription
forall a. Proposal a -> PDescription
prDescription Proposal a
m PDescription -> PDescription -> Bool
forall a. Eq a => a -> a -> Bool
== Proposal a -> PDescription
forall a. Proposal a -> PDescription
prDescription 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
prDescription Proposal a
p, Proposal a -> PName
forall a. Proposal a -> PName
prName Proposal a
p, Proposal a -> PWeight
forall a. Proposal a -> PWeight
prWeight Proposal a
p))

-- | Ratio of the proposal kernels.
--
-- For unbiased, volume preserving proposals, the values is 1.0.
--
-- For biased proposals, the kernel ratio is qYX / qXY, where qAB is the
-- probability density to move from A to B.
type KernelRatio = Log Double

-- | Proposal result.
data PResult a
  = -- | Accept the new value regardless of the prior, likelihood or Jacobian.
    ForceAccept !a
  | -- | Reject the proposal regardless of the prior, likelihood or Jacobian.
    ForceReject
  | -- | Propose a new value.
    --
    -- In order to calculate the Metropolis-Hastings-Green ratio, we need to know
    -- the ratio of the backward to forward kernels (the 'KernelRatio' or the
    -- probability masses or probability densities) and the 'Jacobian'.
    --
    -- NOTE: Actually the 'Jacobian' should be part of the 'KernelRatio'. However,
    -- it is more declarative to have them separate. Like so, we are constantly
    -- reminded: Is the Jacobian modifier different from 1.0?
    Propose !a !KernelRatio !Jacobian
  deriving (Int -> PResult a -> ShowS
[PResult a] -> ShowS
PResult a -> String
(Int -> PResult a -> ShowS)
-> (PResult a -> String)
-> ([PResult a] -> ShowS)
-> Show (PResult a)
forall a. Show a => Int -> PResult a -> ShowS
forall a. Show a => [PResult a] -> ShowS
forall a. Show a => PResult a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [PResult a] -> ShowS
$cshowList :: forall a. Show a => [PResult a] -> ShowS
show :: PResult a -> String
$cshow :: forall a. Show a => PResult a -> String
showsPrec :: Int -> PResult a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> PResult a -> ShowS
Show, PResult a -> PResult a -> Bool
(PResult a -> PResult a -> Bool)
-> (PResult a -> PResult a -> Bool) -> Eq (PResult a)
forall a. Eq a => PResult a -> PResult a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: PResult a -> PResult a -> Bool
$c/= :: forall a. Eq a => PResult a -> PResult a -> Bool
== :: PResult a -> PResult a -> Bool
$c== :: forall a. Eq a => PResult a -> PResult a -> Bool
Eq)

-- | Lift a proposal from one data type to another.
--
-- Assume the Jacobian is 1.0.
--
-- For example:
--
-- @
-- scaleFirstEntryOfTuple = _1 @~ scale
-- @
--
-- See also 'liftProposalWith'.
infixl 7 @~

(@~) :: Lens' b a -> Proposal a -> Proposal b
@~ :: Lens' b a -> Proposal a -> Proposal b
(@~) = Lens' b a -> Proposal a -> Proposal b
forall b a. Lens' b a -> Proposal a -> Proposal b
liftProposal

-- | See '(@~)'.
liftProposal :: Lens' b a -> Proposal a -> Proposal b
liftProposal :: Lens' b a -> Proposal a -> Proposal b
liftProposal = JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b
forall b a.
JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b
liftProposalWith (JacobianG Double -> JacobianFunction b
forall a b. a -> b -> a
const JacobianG Double
1.0)

-- | Lift a proposal from one data type to another.
--
-- A function to calculate the Jacobian has to be provided. If the Jabobian is
-- 1.0, use 'liftProposal'.
--
-- For further reference, please see the [example
-- @Pair@](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Pair/Pair.hs).
liftProposalWith :: JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b
liftProposalWith :: JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b
liftProposalWith JacobianFunction b
jf Lens' b a
l (Proposal PName
n PDescription
r PSpeed
d PDimension
p PWeight
w PFunction a
s Maybe (Tuner a)
t) =
  PName
-> PDescription
-> PSpeed
-> PDimension
-> PWeight
-> PFunction b
-> Maybe (Tuner b)
-> Proposal b
forall a.
PName
-> PDescription
-> PSpeed
-> PDimension
-> PWeight
-> PFunction a
-> Maybe (Tuner a)
-> Proposal a
Proposal PName
n PDescription
r PSpeed
d PDimension
p PWeight
w (JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b
forall b a.
JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b
liftPFunctionWith JacobianFunction b
jf Lens' b a
l PFunction a
s) (JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b
forall b a. JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b
liftTunerWith JacobianFunction b
jf 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)

-- TODO @Dominik (high, feature): Proposals should be aware of: Is this Burn in, or not?

-- | Simple proposal function without tuning information.
--
-- Instruction about randomly moving from the current state to a new state,
-- given some source of randomness.
--
-- Maybe report acceptance counts internal to the proposal (e.g., used by
-- proposals based on Hamiltonian dynamics).
type PFunction a = a -> IOGenM StdGen -> IO (PResult a, Maybe AcceptanceCounts)

-- Lift a proposal function from one data type to another.
liftPFunctionWith :: JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b
liftPFunctionWith :: JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b
liftPFunctionWith JacobianFunction b
jf Lens' b a
l PFunction a
s = PFunction b
s'
  where
    s' :: PFunction b
s' b
y IOGenM StdGen
g = do
      (PResult a
pr, Maybe AcceptanceCounts
ac) <- PFunction a
s (b
y b -> Getting a b a -> a
forall s a. s -> Getting a s a -> a
^. Getting a b a
Lens' b a
l) IOGenM StdGen
g
      let pr' :: PResult b
pr' = case PResult a
pr of
            ForceAccept a
x' -> b -> PResult b
forall a. a -> PResult a
ForceAccept (b -> PResult b) -> b -> PResult b
forall a b. (a -> b) -> a -> b
$ 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
y
            PResult a
ForceReject -> PResult b
forall a. PResult a
ForceReject
            Propose a
x' JacobianG Double
r JacobianG Double
j ->
              let y' :: b
y' = 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
y
                  jxy :: JacobianG Double
jxy = JacobianFunction b
jf b
y
                  jyx :: JacobianG Double
jyx = JacobianFunction b
jf b
y'
                  j' :: JacobianG Double
j' = JacobianG Double
j JacobianG Double -> JacobianG Double -> JacobianG Double
forall a. Num a => a -> a -> a
* JacobianG Double
jyx JacobianG Double -> JacobianG Double -> JacobianG Double
forall a. Fractional a => a -> a -> a
/ JacobianG Double
jxy
               in b -> JacobianG Double -> JacobianG Double -> PResult b
forall a. a -> JacobianG Double -> JacobianG Double -> PResult a
Propose b
y' JacobianG Double
r JacobianG Double
j'
      (PResult b, Maybe AcceptanceCounts)
-> IO (PResult b, Maybe AcceptanceCounts)
forall (f :: * -> *) a. Applicative f => a -> f a
pure (PResult b
pr', Maybe AcceptanceCounts
ac)

-- | Create a proposal with a single tuning parameter.
--
-- Proposals with auxiliary tuning parameters have to be created manually. See
-- 'Tuner' for more information, and 'Mcmc.Proposal.Hamiltonian' for an example.
createProposal ::
  -- | Description of the proposal type and parameters.
  PDescription ->
  -- | Function creating a simple proposal function for a given tuning parameter.
  (TuningParameter -> PFunction a) ->
  -- | Speed.
  PSpeed ->
  -- | Dimension.
  PDimension ->
  -- | Name.
  PName ->
  -- | Weight.
  PWeight ->
  -- | Activate tuning?
  Tune ->
  Proposal a
createProposal :: PDescription
-> (Double -> PFunction a)
-> PSpeed
-> PDimension
-> PName
-> PWeight
-> Tune
-> Proposal a
createProposal PDescription
r Double -> PFunction a
f PSpeed
s PDimension
d PName
n PWeight
w Tune
Tune =
  PName
-> PDescription
-> PSpeed
-> PDimension
-> PWeight
-> PFunction a
-> Maybe (Tuner a)
-> Proposal a
forall a.
PName
-> PDescription
-> PSpeed
-> PDimension
-> PWeight
-> PFunction a
-> Maybe (Tuner a)
-> Proposal a
Proposal PName
n PDescription
r PSpeed
s PDimension
d PWeight
w (Double -> PFunction a
f Double
1.0) (Tuner a -> Maybe (Tuner a)
forall a. a -> Maybe a
Just Tuner a
tuner)
  where
    fT :: TuningFunction a
fT = TuningFunction a
forall a. TuningFunction a
tuningFunction
    g :: Double -> p -> Either a (PFunction a)
g Double
t p
_ = PFunction a -> Either a (PFunction a)
forall a b. b -> Either a b
Right (PFunction a -> Either a (PFunction a))
-> PFunction a -> Either a (PFunction a)
forall a b. (a -> b) -> a -> b
$ Double -> PFunction a
f Double
t
    tuner :: Tuner a
tuner = Double
-> AuxiliaryTuningParameters
-> Bool
-> TuningFunction a
-> (Double
    -> AuxiliaryTuningParameters -> Either String (PFunction a))
-> Tuner a
forall a.
Double
-> AuxiliaryTuningParameters
-> Bool
-> TuningFunction a
-> (Double
    -> AuxiliaryTuningParameters -> Either String (PFunction a))
-> Tuner a
Tuner Double
1.0 AuxiliaryTuningParameters
forall a. Unbox a => Vector a
VU.empty Bool
False TuningFunction a
forall a. TuningFunction a
fT Double -> AuxiliaryTuningParameters -> Either String (PFunction a)
forall p a. Double -> p -> Either a (PFunction a)
g
createProposal PDescription
r Double -> PFunction a
f PSpeed
s PDimension
d PName
n PWeight
w Tune
NoTune =
  PName
-> PDescription
-> PSpeed
-> PDimension
-> PWeight
-> PFunction a
-> Maybe (Tuner a)
-> Proposal a
forall a.
PName
-> PDescription
-> PSpeed
-> PDimension
-> PWeight
-> PFunction a
-> Maybe (Tuner a)
-> Proposal a
Proposal PName
n PDescription
r PSpeed
s PDimension
d PWeight
w (Double -> PFunction a
f Double
1.0) Maybe (Tuner a)
forall a. Maybe a
Nothing

-- | Required information to tune 'Proposal's.
data Tuner a = Tuner
  { Tuner a -> Double
tTuningParameter :: TuningParameter,
    Tuner a -> AuxiliaryTuningParameters
tAuxiliaryTuningParameters :: AuxiliaryTuningParameters,
    -- | Does the tuner require the trace over the last tuning period?
    Tuner a -> Bool
tRequireTrace :: Bool,
    Tuner a -> TuningFunction a
tTuningFunction :: TuningFunction a,
    -- | Given the tuning parameter, and the auxiliary tuning parameters, get
    -- the tuned propose function.
    --
    -- Should return 'Left' if the vector of auxiliary tuning parameters is
    -- invalid.
    Tuner a
-> Double
-> AuxiliaryTuningParameters
-> Either String (PFunction a)
tGetPFunction ::
      TuningParameter ->
      AuxiliaryTuningParameters ->
      Either String (PFunction a)
  }

-- Lift tuner from one data type to another.
liftTunerWith :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b
liftTunerWith :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b
liftTunerWith JacobianFunction b
jf Lens' b a
l (Tuner Double
p AuxiliaryTuningParameters
ps Bool
nt TuningFunction a
fP Double -> AuxiliaryTuningParameters -> Either String (PFunction a)
g) = Double
-> AuxiliaryTuningParameters
-> Bool
-> TuningFunction b
-> (Double
    -> AuxiliaryTuningParameters -> Either String (PFunction b))
-> Tuner b
forall a.
Double
-> AuxiliaryTuningParameters
-> Bool
-> TuningFunction a
-> (Double
    -> AuxiliaryTuningParameters -> Either String (PFunction a))
-> Tuner a
Tuner Double
p AuxiliaryTuningParameters
ps Bool
nt TuningFunction b
fP' Double -> AuxiliaryTuningParameters -> Either String (PFunction b)
g'
  where
    fP' :: TuningFunction b
fP' TuningType
b PDimension
d Double
r = TuningFunction a
fP TuningType
b PDimension
d Double
r (Maybe (Vector a)
 -> (Double, AuxiliaryTuningParameters)
 -> (Double, AuxiliaryTuningParameters))
-> (Maybe (Vector b) -> Maybe (Vector a))
-> Maybe (Vector b)
-> (Double, AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector b -> Vector a) -> Maybe (Vector b) -> Maybe (Vector a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((b -> a) -> Vector b -> Vector a
forall a b. (a -> b) -> Vector a -> Vector b
VB.map (b -> Getting a b a -> a
forall s a. s -> Getting a s a -> a
^. Getting a b a
Lens' b a
l))
    g' :: Double -> AuxiliaryTuningParameters -> Either String (PFunction b)
g' Double
x AuxiliaryTuningParameters
xs = JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b
forall b a.
JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b
liftPFunctionWith JacobianFunction b
jf Lens' b a
l (PFunction a -> PFunction b)
-> Either String (PFunction a) -> Either String (PFunction b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> AuxiliaryTuningParameters -> Either String (PFunction a)
g Double
x AuxiliaryTuningParameters
xs

-- | Tune 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)

-- | Tuning parameter.
--
--  The larger the tuning parameter, the larger the proposal and the lower the
-- expected acceptance rate; and vice versa.
type TuningParameter = Double

-- | The last tuning step may be special.
data TuningType = NormalTuningStep | LastTuningStep

-- | Compute new tuning parameters.
type TuningFunction a =
  TuningType ->
  PDimension ->
  -- | Acceptance rate of last tuning period.
  AcceptanceRate ->
  -- | Trace of last tuning period. Only available when requested by proposal.
  Maybe (VB.Vector a) ->
  (TuningParameter, AuxiliaryTuningParameters) ->
  (TuningParameter, AuxiliaryTuningParameters)

-- | Auxiliary tuning parameters.
--
-- Auxiliary tuning parameters are not shown in proposal summaries.
--
-- Vector may be empty.
type AuxiliaryTuningParameters = VU.Vector TuningParameter

tuningFunctionSimple :: PDimension -> AcceptanceRate -> TuningParameter -> TuningParameter
tuningFunctionSimple :: PDimension -> Double -> Double -> Double
tuningFunctionSimple PDimension
d Double
r Double
t = let rO :: Double
rO = PDimension -> Double
getOptimalRate PDimension
d in Double -> Double
forall a. Floating a => a -> a
exp (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rO)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t

-- | Default tuning function.
--
-- The default tuning function only uses the acceptance rate. In particular, it
-- does not handle auxiliary tuning parameters and ignores the actual samples
-- attained during the last tuning period.
tuningFunction :: TuningFunction a
tuningFunction :: TuningFunction a
tuningFunction TuningType
_ PDimension
d Double
r Maybe (Vector a)
_ = (Double -> Double)
-> (AuxiliaryTuningParameters -> AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap (PDimension -> Double -> Double -> Double
tuningFunctionSimple PDimension
d Double
r) AuxiliaryTuningParameters -> AuxiliaryTuningParameters
forall a. a -> a
id

-- | Also tune auxiliary tuning parameters.
tuningFunctionWithAux ::
  -- | Auxiliary tuning function.
  (TuningType -> VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters) ->
  TuningFunction a
tuningFunctionWithAux :: (TuningType
 -> Vector a
 -> AuxiliaryTuningParameters
 -> AuxiliaryTuningParameters)
-> TuningFunction a
tuningFunctionWithAux TuningType
-> Vector a
-> AuxiliaryTuningParameters
-> AuxiliaryTuningParameters
_ TuningType
_ PDimension
_ Double
_ Maybe (Vector a)
Nothing = String
-> (Double, AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
forall a. HasCallStack => String -> a
error String
"tuningFunctionWithAux: empty trace"
tuningFunctionWithAux TuningType
-> Vector a
-> AuxiliaryTuningParameters
-> AuxiliaryTuningParameters
f TuningType
b PDimension
d Double
r (Just Vector a
xs) = (Double -> Double)
-> (AuxiliaryTuningParameters -> AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap (PDimension -> Double -> Double -> Double
tuningFunctionSimple PDimension
d Double
r) (TuningType
-> Vector a
-> AuxiliaryTuningParameters
-> AuxiliaryTuningParameters
f TuningType
b Vector a
xs)

-- | Only tune auxiliary tuning parameters.
tuningFunctionOnlyAux ::
  -- | Auxiliary tuning function.
  (TuningType -> VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters) ->
  TuningFunction a
tuningFunctionOnlyAux :: (TuningType
 -> Vector a
 -> AuxiliaryTuningParameters
 -> AuxiliaryTuningParameters)
-> TuningFunction a
tuningFunctionOnlyAux TuningType
-> Vector a
-> AuxiliaryTuningParameters
-> AuxiliaryTuningParameters
_ TuningType
_ PDimension
_ Double
_ Maybe (Vector a)
Nothing = String
-> (Double, AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
forall a. HasCallStack => String -> a
error String
"tuningFunctionOnlyAux: empty trace"
tuningFunctionOnlyAux TuningType
-> Vector a
-> AuxiliaryTuningParameters
-> AuxiliaryTuningParameters
f TuningType
b PDimension
_ Double
_ (Just Vector a
xs) = (Double -> Double)
-> (AuxiliaryTuningParameters -> AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
-> (Double, AuxiliaryTuningParameters)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap Double -> Double
forall a. a -> a
id (TuningType
-> Vector a
-> AuxiliaryTuningParameters
-> AuxiliaryTuningParameters
f TuningType
b Vector a
xs)

-- IDEA: Per proposal type tuning parameter boundaries. For example, a sliding
-- proposal with a large tuning parameter is not a problem. But then, if the
-- tuning parameters are very different from one, a different base proposal
-- should be chosen.

-- | Minimal tuning parameter; subject to change.
tuningParameterMin :: TuningParameter
tuningParameterMin :: Double
tuningParameterMin = Double
1e-5

-- | Maximal tuning parameter; subject to change.
tuningParameterMax :: TuningParameter
tuningParameterMax :: Double
tuningParameterMax = Double
1e3

-- | Tune a 'Proposal'.
--
-- The size of the proposal is proportional to the tuning parameter which has
-- positive lower and upper boundaries of 'tuningParameterMin' and
-- 'tuningParameterMax', respectively.
--
-- Auxiliary tuning parameters may also be used by the 'Tuner' of the proposal.
--
-- Return 'Left' if:
--
-- - the 'Proposal' is not tuneable;
--
-- - the auxiliary tuning parameters are invalid.
--
-- Used by 'Mcmc.Chain.Save.fromSavedChain'.
tuneWithTuningParameters ::
  TuningParameter ->
  AuxiliaryTuningParameters ->
  Proposal a ->
  Either String (Proposal a)
tuneWithTuningParameters :: Double
-> AuxiliaryTuningParameters
-> Proposal a
-> Either String (Proposal a)
tuneWithTuningParameters Double
t AuxiliaryTuningParameters
ts Proposal a
p = case Proposal a -> Maybe (Tuner a)
forall a. Proposal a -> Maybe (Tuner a)
prTuner Proposal a
p of
  Maybe (Tuner a)
Nothing -> String -> Either String (Proposal a)
forall a b. a -> Either a b
Left String
"tuneWithTuningParameters: Proposal is not tunable."
  Just (Tuner Double
_ AuxiliaryTuningParameters
_ Bool
nt TuningFunction a
fT Double -> AuxiliaryTuningParameters -> Either String (PFunction a)
g) ->
    -- Ensure that the tuning parameter is strictly positive and well bounded.
    let t' :: Double
t' = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
tuningParameterMin Double
t
        t'' :: Double
t'' = Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
tuningParameterMax Double
t'
        psE :: Either String (PFunction a)
psE = Double -> AuxiliaryTuningParameters -> Either String (PFunction a)
g Double
t'' AuxiliaryTuningParameters
ts
     in case Either String (PFunction a)
psE of
          Left String
err -> String -> Either String (Proposal a)
forall a b. a -> Either a b
Left (String -> Either String (Proposal a))
-> String -> Either String (Proposal a)
forall a b. (a -> b) -> a -> b
$ String
"tune: " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
err
          Right PFunction a
ps -> Proposal a -> Either String (Proposal a)
forall a b. b -> Either a b
Right (Proposal a -> Either String (Proposal a))
-> Proposal a -> Either String (Proposal a)
forall a b. (a -> b) -> a -> b
$ Proposal a
p {prFunction :: PFunction a
prFunction = PFunction a
ps, prTuner :: Maybe (Tuner a)
prTuner = 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
-> AuxiliaryTuningParameters
-> Bool
-> TuningFunction a
-> (Double
    -> AuxiliaryTuningParameters -> Either String (PFunction a))
-> Tuner a
forall a.
Double
-> AuxiliaryTuningParameters
-> Bool
-> TuningFunction a
-> (Double
    -> AuxiliaryTuningParameters -> Either String (PFunction a))
-> Tuner a
Tuner Double
t'' AuxiliaryTuningParameters
ts Bool
nt TuningFunction a
fT Double -> AuxiliaryTuningParameters -> Either String (PFunction a)
g}

-- | See 'PDimension'.
getOptimalRate :: PDimension -> Double
getOptimalRate :: PDimension -> Double
getOptimalRate (PDimension Int
n)
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = String -> Double
forall a. HasCallStack => String -> a
error String
"getOptimalRate: Proposal dimension is zero or negative."
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = Double
0.44
  -- Use a linear interpolation with delta 0.0515.
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2 = Double
0.3885
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 = Double
0.337
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
4 = Double
0.2855
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
5 = Double
0.234
  | Bool
otherwise = String -> Double
forall a. HasCallStack => String -> a
error String
"getOptimalRate: Proposal dimension is not an integer?"
getOptimalRate PDimension
PDimensionUnknown = Double
0.234
getOptimalRate (PSpecial Int
_ Double
r) = Double
r

-- Warn if acceptance rate is lower.
rateMin :: Double
rateMin :: Double
rateMin = Double
0.1

-- Warn if acceptance rate is larger.
rateMax :: Double
rateMax :: Double
rateMax = Double
0.9

renderRow ::
  BL.ByteString ->
  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
-> ByteString
renderRow ByteString
name ByteString
ptype ByteString
weight ByteString
nAccept ByteString
nReject ByteString
acceptRate ByteString
optimalRate ByteString
tuneParam ByteString
manualAdjustment = 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
ro 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
14 ByteString
nAccept
    nr :: ByteString
nr = Int -> ByteString -> ByteString
alignRight Int
14 ByteString
nReject
    ra :: ByteString
ra = Int -> ByteString -> ByteString
alignRight Int
14 ByteString
acceptRate
    ro :: ByteString
ro = Int -> ByteString -> ByteString
alignRight Int
14 ByteString
optimalRate
    tp :: ByteString
tp = Int -> ByteString -> ByteString
alignRight Int
20 ByteString
tuneParam
    mt :: ByteString
mt = Int -> ByteString -> ByteString
alignRight Int
30 ByteString
manualAdjustment

-- | Header of proposal summaries.
proposalHeader :: BL.ByteString
proposalHeader :: ByteString
proposalHeader =
  ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
-> ByteString
renderRow
    ByteString
"Name"
    ByteString
"Description"
    ByteString
"Weight"
    ByteString
"Accepted"
    ByteString
"Rejected"
    ByteString
"Rate"
    ByteString
"Optimal rate"
    ByteString
"Tuning parameter"
    ByteString
"Consider manual adjustment"

-- | Proposal summary.
summarizeProposal ::
  PName ->
  PDescription ->
  PWeight ->
  Maybe TuningParameter ->
  PDimension ->
  Maybe (Int, Int, Double) ->
  BL.ByteString
summarizeProposal :: PName
-> PDescription
-> PWeight
-> Maybe Double
-> PDimension
-> Maybe (Int, Int, Double)
-> ByteString
summarizeProposal PName
name PDescription
description PWeight
weight Maybe Double
tuningParameter PDimension
dimension Maybe (Int, Int, Double)
ar =
  ByteString
-> 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
name)
    (String -> ByteString
BL.pack (String -> ByteString) -> String -> ByteString
forall a b. (a -> b) -> a -> b
$ PDescription -> String
fromPDescription PDescription
description)
    ByteString
weightStr
    ByteString
nAccept
    ByteString
nReject
    ByteString
acceptRate
    ByteString
optimalRate
    ByteString
tuneParamStr
    ByteString
manualAdjustmentStr
  where
    weightStr :: ByteString
weightStr = 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
weight
    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)
ar
    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)
ar
    acceptRate :: ByteString
acceptRate = 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
2 (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)
ar
    optimalRate :: ByteString
optimalRate = ByteString -> ByteString
BL.fromStrict (ByteString -> ByteString) -> ByteString -> ByteString
forall a b. (a -> b) -> a -> b
$ Int -> Double -> ByteString
BC.toFixed Int
2 (Double -> ByteString) -> Double -> ByteString
forall a b. (a -> b) -> a -> b
$ PDimension -> Double
getOptimalRate PDimension
dimension
    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
4) Maybe Double
tuningParameter
    checkRate :: Double -> Maybe a
checkRate Double
rate
      | Double
rate Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
rateMin = a -> Maybe a
forall a. a -> Maybe a
Just a
"rate too low"
      | Double
rate Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
rateMax = a -> Maybe a
forall a. a -> Maybe a
Just a
"rate too high"
      | Bool
otherwise = Maybe a
forall a. Maybe a
Nothing
    checkTuningParam :: Double -> Maybe a
checkTuningParam Double
tp
      | Double
tp Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= (Double
1.1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tuningParameterMin) = a -> Maybe a
forall a. a -> Maybe a
Just a
"tuning parameter too low"
      | Double
tp Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= (Double
0.9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tuningParameterMax) = a -> Maybe a
forall a. a -> Maybe a
Just a
"tuning parameter too high"
      | Bool
otherwise = Maybe a
forall a. Maybe a
Nothing
    tps :: Maybe ByteString
tps = Double -> Maybe ByteString
forall a. IsString a => Double -> Maybe a
checkTuningParam (Double -> Maybe ByteString) -> Maybe Double -> Maybe ByteString
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Maybe Double
tuningParameter
    ars :: Maybe ByteString
ars = (Double -> Maybe ByteString
forall a. IsString a => Double -> Maybe a
checkRate (Double -> Maybe ByteString)
-> ((Int, Int, Double) -> Double)
-> (Int, Int, Double)
-> Maybe 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)) ((Int, Int, Double) -> Maybe ByteString)
-> Maybe (Int, Int, Double) -> Maybe ByteString
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Maybe (Int, Int, Double)
ar
    manualAdjustmentStr :: ByteString
manualAdjustmentStr =
      let
       in case (Maybe ByteString
ars, Maybe ByteString
tps) of
            (Maybe ByteString
Nothing, Maybe ByteString
Nothing) -> ByteString
""
            (Just ByteString
s, Maybe ByteString
_) -> ByteString
s
            (Maybe ByteString
_, Just ByteString
s) -> ByteString
s