API

@model(expr[, warn = false])

Macro to specify a probabilistic model.

If warn is true, a warning is displayed if internal variable names are used in the model definition.

Examples

Model definition:

@model function model(x, y = 42)
    ...
end

To generate a Model, call model(xvalue) or model(xvalue, yvalue).

struct Model{F,argnames,defaultnames,missings,Targs,Tdefaults,Ctx<:AbstractContext,Threaded}
    f::F
    args::NamedTuple{argnames,Targs}
    defaults::NamedTuple{defaultnames,Tdefaults}
    context::Ctx=DefaultContext()
end

A Model struct with model evaluation function of type F, arguments of names argnames types Targs, default arguments of names defaultnames with types Tdefaults, missing arguments missings, and evaluation context of type Ctx.

Here argnames, defaultargnames, and missings are tuples of symbols, e.g. (:a, :b). context is by default DefaultContext().

An argument with a type of Missing will be in missings by default. However, in non-traditional use-cases missings can be defined differently. All variables in missings are treated as random variables rather than observations.

The Threaded type parameter indicates whether the model requires threadsafe evaluation (i.e., whether the model contains statements which modify the internal VarInfo that are executed in parallel). By default, this is set to false.

The default arguments are used internally when constructing instances of the same model with different arguments.

Examples

julia> Model(f, (x = 1.0, y = 2.0))
Model{typeof(f),(:x, :y),(),(),Tuple{Float64,Float64},Tuple{}}(f, (x = 1.0, y = 2.0), NamedTuple())

julia> Model(f, (x = 1.0, y = 2.0), (x = 42,))
Model{typeof(f),(:x, :y),(:x,),(),Tuple{Float64,Float64},Tuple{Int64}}(f, (x = 1.0, y = 2.0), (x = 42,))

julia> Model{(:y,)}(f, (x = 1.0, y = 2.0), (x = 42,)) # with special definition of missings
Model{typeof(f),(:x, :y),(:x,),(:y,),Tuple{Float64,Float64},Tuple{Int64}}(f, (x = 1.0, y = 2.0), (x = 42,))

(model::Model)([rng, varinfo])

Sample from the prior of the model with random number generator rng.

Returns the model's return value.

Note that calling this with an existing varinfo object will mutate it.

nameof(model::Model)

Get the name of the model as Symbol.

getargnames(model::Model)

Get a tuple of the argument names of the model.

getmissings(model::Model)

Get a tuple of the names of the missing arguments of the model.

contextualize(model::Model, context::AbstractContext)

Return a new Model with the same evaluation function and other arguments, but with its underlying context set to context.

setthreadsafe(model::Model, threadsafe::Bool)

Returns a new Model with its threadsafe flag set to threadsafe.

Threadsafe evaluation ensures correctness when executing model statements that mutate the internal VarInfo object in parallel. For example, this is needed if tilde-statements are nested inside Threads.@threads or similar constructs.

It is not needed for generic multithreaded operations that don't involve VarInfo. For example, calculating a log-likelihood term in parallel and then calling @addlogprob! outside of the parallel region is safe without needing to set threadsafe=true.

It is also not needed for multithreaded sampling with AbstractMCMC's MCMCThreads().

Setting threadsafe to true increases the overhead in evaluating the model. Please see the Turing.jl docs for more details.

requires_threadsafe(model::Model)

Return whether model has been marked as needing threadsafe evaluation (using setthreadsafe).

rand([rng=Random.default_rng()], model::Model)

Sample a VarNamedTuple of raw values from the prior of model.

logprior(model::Model, params)
logprior(model::Model, varinfo::AbstractVarInfo)

Return the log prior probability of variables params for the probabilistic model, or the log prior of the data in varinfo if provided.

Note that this probability always refers to the parameters in unlinked space, i.e., the return value of logprior does not depend on whether VarInfo has been linked or not.

See also logjoint and loglikelihood.

Examples

julia> @model function demo(x)
           m ~ Normal()
           for i in eachindex(x)
               x[i] ~ Normal(m, 1.0)
           end
       end
demo (generic function with 2 methods)

julia> # Using a `NamedTuple`.
       logprior(demo([1.0]), (m = 100.0, ))
-5000.918938533205

julia> # Using a `OrderedDict`.
       logprior(demo([1.0]), OrderedDict(@varname(m) => 100.0))
-5000.918938533205

julia> # Truth.
       logpdf(Normal(), 100.0)
-5000.918938533205

logprior(model::DynamicPPL.Model, chain::MCMCChains.Chains)

Return an array of log prior probabilities evaluated at each sample in an MCMC chain.

Examples

julia> using MCMCChains, Distributions

julia> @model function demo_model(x)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           for i in eachindex(x)
               x[i] ~ Normal(m, sqrt(s))
           end
       end;

julia> # Construct a chain of samples using MCMCChains.
       # This sets s = 0.5 and m = 1.0 for all three samples.
       chain = Chains(repeat([0.5 1.0;;;], 3, 1, 1), [:s, :m]);

julia> logprior(demo_model([1., 2.]), chain)
3×1 Matrix{Float64}:
 -3.2956988239086447
 -3.2956988239086447
 -3.2956988239086447

loglikelihood(model::Model, params)
loglikelihood(model::Model, varinfo::AbstractVarInfo)

Return the log likelihood of variables params for the probabilistic model, or the log likelihood of the data in varinfo if provided.

See also logjoint and logprior.

Examples

```jldoctest; setup=:(using Distributions) julia> @model function demo(x) m ~ Normal() for i in eachindex(x) x[i] ~ Normal(m, 1.0) end end demo (generic function with 2 methods)

julia> # Using a NamedTuple. loglikelihood(demo([1.0]), (m = 100.0, )) -4901.418938533205

julia> # Using a OrderedDict. loglikelihood(demo([1.0]), OrderedDict(@varname(m) => 100.0)) -4901.418938533205

julia> # Truth. logpdf(Normal(100.0, 1.0), 1.0) -4901.418938533205

loglikelihood(model::DynamicPPL.Model, chain::MCMCChains.Chains)

Return an array of log likelihoods evaluated at each sample in an MCMC chain.

Examples

julia> using MCMCChains, Distributions

julia> @model function demo_model(x)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           for i in eachindex(x)
               x[i] ~ Normal(m, sqrt(s))
           end
       end;

julia> # Construct a chain of samples using MCMCChains.
       # This sets s = 0.5 and m = 1.0 for all three samples.
       chain = Chains(repeat([0.5 1.0;;;], 3, 1, 1), [:s, :m]);

julia> loglikelihood(demo_model([1., 2.]), chain)
3×1 Matrix{Float64}:
 -2.1447298858494
 -2.1447298858494
 -2.1447298858494

logjoint(model::Model, params)
logjoint(model::Model, varinfo::AbstractVarInfo)

Return the log joint probability of variables params for the probabilistic model, or the log joint of the data in varinfo if provided.

Note that this probability always refers to the parameters in unlinked space, i.e., the return value of logjoint does not depend on whether VarInfo has been linked or not.

See also logprior and loglikelihood.

Examples

julia> @model function demo(x)
           m ~ Normal()
           for i in eachindex(x)
               x[i] ~ Normal(m, 1.0)
           end
       end
demo (generic function with 2 methods)

julia> # Using a `NamedTuple`.
       logjoint(demo([1.0]), (m = 100.0, ))
-9902.33787706641

julia> # Using a `OrderedDict`.
       logjoint(demo([1.0]), OrderedDict(@varname(m) => 100.0))
-9902.33787706641

julia> # Truth.
       logpdf(Normal(100.0, 1.0), 1.0) + logpdf(Normal(), 100.0)
-9902.33787706641

logjoint(model::Model, chain::MCMCChains.Chains)

Return an array of log joint probabilities evaluated at each sample in an MCMC chain.

Examples

julia> using MCMCChains, Distributions

julia> @model function demo_model(x)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           for i in eachindex(x)
               x[i] ~ Normal(m, sqrt(s))
           end
       end;

julia> # Construct a chain of samples using MCMCChains.
       # This sets s = 0.5 and m = 1.0 for all three samples.
       chain = Chains(repeat([0.5 1.0;;;], 3, 1, 1), [:s, :m]);

julia> logjoint(demo_model([1., 2.]), chain)
3×1 Matrix{Float64}:
 -5.440428709758045
 -5.440428709758045
 -5.440428709758045

DynamicPPL.LogDensityFunction(
    model::Model,
    getlogdensity::Any=getlogjoint_internal,
    vnt_or_vi::Union{VarNamedTuple,VarInfo,OnlyAccsVarInfo}=VarInfo(model).values,
    accs::Union{NTuple{<:Any,AbstractAccumulator},AccumulatorTuple}=DynamicPPL.ldf_accs(getlogdensity);
    adtype::Union{ADTypes.AbstractADType,Nothing}=nothing,
)

A struct which contains a model, along with all the information necessary to:

• calculate its log density at a given point;
• and if adtype is provided, calculate the gradient of the log density at that point.

This information can be extracted using the LogDensityProblems.jl interface, specifically, using LogDensityProblems.logdensity and LogDensityProblems.logdensity_and_gradient. If adtype is nothing, then only logdensity is implemented. If adtype is a concrete AD backend type, then logdensity_and_gradient is also implemented.

getlogdensity should be a callable which takes a single argument: an OnlyAccsVarInfo, and returns a Real corresponding to the log density of interest. There are several functions in DynamicPPL that are 'supported' out of the box:

• getlogjoint_internal: calculate the log joint, including the log-Jacobian term for any variables that have been linked in the provided VarInfo.
• getlogprior_internal: calculate the log prior, including the log-Jacobian term for any variables that have been linked in the provided VarInfo.
• getlogjoint: calculate the log joint in the model space, ignoring any effects of linking
• getlogprior: calculate the log prior in the model space, ignoring any effects of linking
• getloglikelihood: calculate the log likelihood (this is unaffected by linking, since transforms are only applied to random variables)

vnt is a VarNamedTuple which contains vectorised representations of all the random variables in the model. You can create one manually, but for this argument, it is easier to pass either:

• an oavi::OnlyAccsVarInfo which contains a VectorValueAccumulator, in which case the values are taken from that accumulator (this approach is recommended going forward);
• a vi::VarInfo, in which case the vi.values field is used.

accs allows you to specify an AccumulatorTuple or a tuple of AbstractAccumulators which will be used when evaluating the log density. (Note that any accumulators from the previous argument are discarded.) By default, this uses an internal function,DynamicPPL.ldf_accs`, which attempts to choose an appropriate set of accumulators based on which kind of log-density is being calculated.

If the adtype keyword argument is provided, then this struct will also store the adtype along with other information for efficient calculation of the gradient of the log density. Note that preparing a LogDensityFunction with an AD type AutoBackend() requires the AD backend itself to have been loaded (e.g. with import Backend).

Fields

Note that it is undefined behaviour to access any of a LogDensityFunction's fields, apart from:

• ldf.model: The original model from which this LogDensityFunction was constructed.
• ldf.adtype: The AD type used for gradient calculations, or nothing if no AD type was provided.
• ldf.transform_strategy: The transform strategy that specifies which variables in the LogDensityFunction are linked or unlinked.

Extended help

Up until DynamicPPL v0.38, there have been two ways of evaluating a DynamicPPL model at a given set of parameters:

1. With unflatten!! + evaluate!! with DefaultContext: this stores a vector of parameters inside a VarInfo's metadata, then reads parameter values from the VarInfo during evaluation.

2. With InitFromParams: this reads parameter values from a NamedTuple or a Dict, and stores them inside a VarInfo's metadata.

In general, both of these approaches work fine, but the fact that they modify the VarInfo's metadata can often be quite wasteful. In particular, it is very common that the only outputs we care about from model evaluation are those which are stored in accumulators, such as log probability densities, or raw values.

To avoid this issue, we use OnlyAccsVarInfo, which is a VarInfo that only contains accumulators. It implements enough of the AbstractVarInfo interface to not error during model evaluation.

Because OnlyAccsVarInfo does not store any parameter values, when evaluating a model with it, it is mandatory that parameters are provided from outside the VarInfo, namely via InitContext.

The main problem that we face is that it is not possible to directly implement DynamicPPL.init(rng, vn, dist, strategy) for strategy::InitFromParams{<:AbstractVector}. In particular, it is not clear:

• which parts of the vector correspond to which random variables, and
• whether the variables are linked or unlinked.

Traditionally, this problem has been solved by unflatten!!, because that function would place values into the VarInfo's metadata alongside the information about ranges and linking. That way, when we evaluate with DefaultContext, we can read this information out again. However, we want to avoid using a metadata. Thus, here, we extract this information from the VarInfo a single time when constructing a LogDensityFunction object. Inside the LogDensityFunction, we store a mapping from VarNames to ranges in that vector, along with link status.

When evaluating the model, this allows us to combine the parameter vector together with those ranges to create an InitFromVector, which lets us very quickly read parameter values from the vector.

Note that this assumes that the ranges and link status are static throughout the lifetime of the LogDensityFunction object. Therefore, a LogDensityFunction object cannot handle models which have variable numbers of parameters, or models which may visit random variables in different orders depending on stochastic control flow. Indeed, silent errors may occur with such models. This is a general limitation of vectorised parameters: the original unflatten!! + evaluate!! approach also fails with such models.

OnlyAccsVarInfo(accs...)

OnlyAccsVarInfo is a wrapper around a tuple of accumulators.

Its name stems from the fact that it implements the minimal AbstractVarInfo interface to work with the tilde_assume!! and tilde_observe!! functions for InitContext.

Note that this does not implement almost every other AbstractVarInfo interface function, and so using this with a different leaf context such as DefaultContext will result in errors.

For more information about accumulators, please see the DynamicPPL documentation on accumulators.

to_vector_params(
    vector_values::VarNamedTuple,
    ldf::LogDensityFunction
)

Extract vectorised values from a VarNamedTuple that contains VectorValues and LinkedVectorValues, and concatenate them into a single vector that is consistent with the ranges specified in the LogDensityFunction.

This is useful when you want to regenerate new vectorised parameters but using a different initialisation strategy.

Note that the transform status of the variables in the VarNamedTuple must be consistent with the transform strategy stored in the LogDensityFunction. This function checks for that.

model | (x = 1.0, ...)

Return a Model which now treats variables on the right-hand side as observations.

See condition for more information and examples.

condition(model::Model; values...)
condition(model::Model, values::NamedTuple)

Return a Model which now treats the variables in values as observations.

See also: decondition, conditioned

Limitations

This does currently not work with variables that are provided to the model as arguments, e.g. @model function demo(x) ... end means that condition will not affect the variable x.

Therefore if one wants to make use of condition and decondition one should not be specifying any random variables as arguments.

This is done for the sake of backwards compatibility.

Examples

Simple univariate model

julia> using Distributions

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> m, x = model(); (m != 1.0 && x != 100.0)
true

julia> # Create a new instance which treats `x` as observed
       # with value `100.0`, and similarly for `m=1.0`.
       conditioned_model = condition(model, x=100.0, m=1.0);

julia> m, x = conditioned_model(); (m == 1.0 && x == 100.0)
true

julia> # Let's only condition on `x = 100.0`.
       conditioned_model = condition(model, x = 100.0);

julia> m, x = conditioned_model(); (m != 1.0 && x == 100.0)
true

julia> # We can also use the nicer `|` syntax.
       conditioned_model = model | (x = 100.0, );

julia> m, x = conditioned_model(); (m != 1.0 && x == 100.0)
true

In the above we have specified the conditioning variables via keyword arguments. You can also provide a NamedTuple, AbstractDict{<:VarName}, or a VarNamedTuple; internally these are all converted to a VarNamedTuple.

For example, here we use a Dict:

julia> conditioned_model_dict = condition(model, Dict(@varname(x) => 100.0));

julia> m, x = conditioned_model_dict(); (m != 1.0 && x == 100.0)
true

julia> # There's also an option using `|` by letting the right-hand side be a tuple
       # with elements of type `Pair{<:VarName}`, i.e. `vn => value` with `vn isa VarName`.
       conditioned_model_pairs = model | (@varname(x) => 100.0);

julia> m, x = conditioned_model_pairs(); (m != 1.0 && x == 100.0)
true

Condition only a part of a multivariate variable

When conditioning on multiple variables at a time, we can use missing to signal that a part of the variable should not be conditioned on.

However, note that in this case each element of the multivariate random variable must be on its own tilde-statement. In other words, if we write m ~ MvNormal(...), then we cannot condition on only m[1]. (In principle, for some distributions this can be possible, specifically when the distribution can be factorised into independent components, like an MvNormal with a diagonal covariance matrix. However, this is not currently implemented.)

julia> @model function demo_mv(::Type{TV}=Float64) where {TV}
           m = Vector{TV}(undef, 2)
           m[1] ~ Normal()
           m[2] ~ Normal()
           return m
       end
demo_mv (generic function with 4 methods)

julia> model = demo_mv();

julia> conditioned_model = condition(model, m = [missing, 1.0]);

julia> # (✓) `m[1]` sampled while `m[2]` is fixed
       m = conditioned_model(); (m[1] != 1.0 && m[2] == 1.0)
true

Intuitively one might also expect to be able to write model | (m[2] = 1.0, ). You cannot do this with a NamedTuple because the VarName m[2] cannot be represented as a Symbol (i.e., Symbol("m[2]") is not the same as @varname(m[2])).

julia> # (×) `m[2]` is not set to 1.0.
       m = condition(model, var"m[2]" = 1.0)(); m[2] == 1.0
false

But you can do this if you use a Dict or a VarNamedTuple as the underlying storage instead:

julia> vnt = @vnt begin
           @template m = zeros(2)
           m[2] := 1.0
       end
VarNamedTuple
└─ m => PartialArray size=(2,) data::Vector{Float64}
        └─ (2,) => 1.0

julia> m = condition(model, vnt)(); (m[1] != 1.0 && m[2] == 1.0)
true

Nested models

condition also supports the use of nested models through the use of to_submodel.

julia> @model demo_inner() = m ~ Normal()
demo_inner (generic function with 2 methods)

julia> @model function demo_outer()
           # By default, `to_submodel` prefixes the variables using the left-hand side of `~`.
           inner ~ to_submodel(demo_inner())
           return inner
       end
demo_outer (generic function with 2 methods)

julia> model = demo_outer();

julia> model() ≠ 1.0
true

julia> # To condition the variable inside `demo_inner` we need to refer to it as `inner.m`.
       conditioned_model = model | (@varname(inner.m) => 1.0, );

julia> conditioned_model()
1.0

julia> # If you attempt to condition on `inner` itself, it must refer to the prefixed
       # latent variables, not the return value. For example, this will work:
       conditioned_model2 = model | (inner = (m = 1.0,), );

julia> conditioned_model2()
1.0

julia> # However, if `inner` does not contain `m` inside it as a field, this will not
       # result in any conditioning:
       conditioned_model_fail = model | (inner = "something else", );

julia> conditioned_model_fail == 1.0
false

conditioned(context::AbstractContext)

Return VarNamedTuple of values that are conditioned on under context.

Note that this will recursively traverse the context stack and return a merged version of the condition values.

conditioned(model::Model)

Return the conditioned values in model.

Examples

julia> using Distributions

julia> using DynamicPPL: conditioned, contextualize, PrefixContext, CondFixContext, Condition

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
       end
demo (generic function with 2 methods)

julia> m = demo();

julia> # Returns all the variables we have conditioned on + their values.
       conditioned(condition(m, x=100.0, m=1.0))
VarNamedTuple
├─ x => 100.0
└─ m => 1.0

julia> # Nested ones also work. (Note that `PrefixContext` also prefixes
       # the variables of any `CondFixContext` that is _inside_ it.)
       new_context = PrefixContext(@varname(a), CondFixContext{Condition}(VarNamedTuple(m=1.0,)));
       cm = condition(contextualize(m, new_context), x=100.0);

julia> conditioned(cm)
VarNamedTuple
├─ a => VarNamedTuple
│       └─ m => 1.0
└─ x => 100.0

julia> # Since we conditioned on `a.m`, it is not treated as a random variable.
       # However, `a.x` will still be a random variable.
       keys(VarInfo(cm))
1-element Vector{VarName}:
 a.x

julia> # We can also condition on `a.m` _outside_ of the PrefixContext:
       cm = condition(contextualize(m, PrefixContext(@varname(a))), (@varname(a.m) => 1.0));

julia> conditioned(cm)
VarNamedTuple
└─ a => VarNamedTuple
        └─ m => 1.0

julia> # Now `a.x` will be sampled.
       keys(VarInfo(cm))
1-element Vector{VarName}:
 a.x

decondition(model::Model)
decondition(model::Model, variables...)

Return a Model for which variables... are not conditioned on. If no variables are provided, then all conditioned variables will be removed.

Note that this function cannot decondition variables that are provided as arguments to the model function itself; it can only decondition variables that were provided via condition or |.

This is essentially the inverse of condition.

Examples

julia> using Distributions

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> conditioned_model = condition(demo(), m = 1.0, x = 10.0);

julia> conditioned_model()
(m = 1.0, x = 10.0)

julia> # By specifying the `VarName` to `decondition`.
       model = decondition(conditioned_model, @varname(m));

julia> (m, x) = model(); (m ≠ 1.0 && x == 10.0)
true

julia> # `decondition` also accepts symbols, although VarNames are preferable for
       # type stability reasons.
       model = decondition(conditioned_model, :m);

julia> (m, x) = model(); (m ≠ 1.0 && x == 10.0)
true

julia> # `decondition` multiple at once:
       (m, x) = decondition(model, :m, :x)(); (m ≠ 1.0 && x ≠ 10.0)
true

julia> # `decondition` without any symbols will `decondition` all variables.
       (m, x) = decondition(model)(); (m ≠ 1.0 && x ≠ 10.0)
true

Note that decondition is only guaranteed to work when you decondition variables that were explicitly provided to condition earlier. In this example we condition on @varname(m) but decondition on @varname(m[1]), which fails because m[1] was not explicitly conditioned on:

julia> @model function demo_mv(::Type{TV}=Float64) where {TV}
           m = Vector{TV}(undef, 2)
           m[1] ~ Normal()
           m[2] ~ Normal()
           return m
       end
demo_mv (generic function with 4 methods)

julia> model = demo_mv();

julia> conditioned_model = condition(model, @varname(m) => [1.0, 2.0]);

julia> conditioned_model()
2-element Vector{Float64}:
 1.0
 2.0

julia> deconditioned_model = decondition(conditioned_model, @varname(m[1]));

julia> deconditioned_model()  # (×) `m[1]` is still conditioned
2-element Vector{Float64}:
 1.0
 2.0

fix(model::Model; values...)
fix(model::Model, values::NamedTuple)

Return a Model which now treats the variables in values as fixed.

See also: unfix, fixed

Examples

Simple univariate model

julia> using Distributions

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> m, x = model(); (m ≠ 1.0 && x ≠ 100.0)
true

julia> # Create a new instance which treats `x` as observed
       # with value `100.0`, and similarly for `m=1.0`.
       fixed_model = fix(model, x=100.0, m=1.0);

julia> m, x = fixed_model(); (m == 1.0 && x == 100.0)
true

julia> # Let's only fix on `x = 100.0`.
       fixed_model = fix(model, x = 100.0);

julia> m, x = fixed_model(); (m != 1.0 && x == 100.0)
true

Other ways of specifying fixed values

Specifying fixed values can be done exactly in the same way as for condition; please see its docstring for more examples.

Difference from condition

The only difference between fixing and conditioning is as follows:

• Conditioned variables are considered to be observations, and are thus included in the computation log-joint and log-likelihood, but not the log-prior.
• Fixed variables are considered to be constant, and are thus not included in any log-probability computations.

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> model_fixed = fix(model, m = 1.0);

julia> model_conditioned = condition(model, m = 1.0);

julia> logjoint(model_fixed, (x=1.0,))
-0.9189385332046728

julia> logjoint(model_conditioned, (x=1.0,))
-2.3378770664093453

julia> # The difference is the missing log-probability of `m`:
       logpdf(Normal(), 1.0)
-1.4189385332046727

fixed(context::AbstractContext)

Return a VarNamedTuple containing the values that are fixed under context.

Note that this will recursively traverse the context stack and return a merged version of the fix values.

fixed(model::Model)

Return the fixed values in model.

Examples

julia> using Distributions

julia> using DynamicPPL: fixed, contextualize, PrefixContext, CondFixContext, Fix

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
       end
demo (generic function with 2 methods)

julia> m = demo();

julia> # Returns all the variables we have fixed on + their values.
       fixed(fix(m, x=100.0, m=1.0))
VarNamedTuple
├─ x => 100.0
└─ m => 1.0

julia> # The rest of this is the same as the `condition` example above.
       fm = fix(contextualize(m, PrefixContext(@varname(a), CondFixContext{Fix}(VarNamedTuple(m=1.0)))), x=100.0);

julia> Set(keys(fixed(fm))) == Set([@varname(a.m), @varname(x)])
true

julia> keys(VarInfo(fm))
1-element Vector{VarName}:
 a.x

julia> # We can also fix `a.m` _outside_ of the PrefixContext:
       fm = fix(contextualize(m, PrefixContext(@varname(a))), (@varname(a.m) => 1.0));

julia> fixed(fm)
VarNamedTuple
└─ a => VarNamedTuple
        └─ m => 1.0

julia> # Now `a.x` will be sampled.
       keys(VarInfo(fm))
1-element Vector{VarName}:
 a.x

unfix(model::Model)
unfix(model::Model, variables...)

Return a Model for which variables... are not considered fixed. If no variables are provided, then all fixed variables will be removed.

This is essentially the inverse of fix.

Conceptually this is very similar to decondition and thus the same limitations apply; please see its docstring for more details.

Examples

julia> using Distributions

julia> @model function demo()
           m ~ Normal()
           x ~ Normal(m, 1)
           return (; m=m, x=x)
       end
demo (generic function with 2 methods)

julia> fixed_model = fix(demo(), m = 1.0, x = 10.0);

julia> fixed_model()
(m = 1.0, x = 10.0)

julia> # By specifying the `VarName` to `unfix`.
       model = unfix(fixed_model, @varname(m));

julia> (m, x) = model(); (m != 1.0 && x == 10.0)
true

julia> # When `NamedTuple` is used as the underlying, you can also provide
       # the symbol directly (though the `@varname` approach is preferable if
       # if the variable is known at compile-time).
       model = unfix(fixed_model, :m);

julia> (m, x) = model(); (m != 1.0 && x == 10.0)
true

julia> # `unfix` multiple at once:
       (m, x) = unfix(model, :m, :x)(); (m != 1.0 && x != 10.0)
true

julia> # `unfix` without any symbols will `unfix` all variables.
       (m, x) = unfix(model)(); (m != 1.0 && x != 10.0)
true

predict([rng::AbstractRNG,] model::Model, chain::MCMCChains.Chains; include_all=false)

Sample from the posterior predictive distribution by executing model with parameters fixed to each sample in chain, and return the resulting Chains.

The model passed to predict is often different from the one used to generate chain. Typically, the model from which chain originated treats certain variables as observed (i.e., data points), while the model you pass to predict may mark these same variables as missing or unobserved. Calling predict then leverages the previously inferred parameter values to simulate what new, unobserved data might look like, given your posterior beliefs.

For each parameter configuration in chain:

1. All random variables present in chain are fixed to their sampled values.
2. Any variables not included in chain are sampled from their prior distributions.

If include_all is false, the returned Chains will contain only those variables that were not fixed by the samples in chain. This is useful when you want to sample only new variables from the posterior predictive distribution.

Examples

using AbstractMCMC, Distributions, DynamicPPL, Random

@model function linear_reg(x, y, σ = 0.1)
    β ~ Normal(0, 1)
    for i in eachindex(y)
        y[i] ~ Normal(β * x[i], σ)
    end
end

# Generate synthetic chain using known ground truth parameter
ground_truth_β = 2.0

# Create chain of samples from a normal distribution centered on ground truth
β_chain = MCMCChains.Chains(
    rand(Normal(ground_truth_β, 0.002), 1000), [:β,]
)

# Generate predictions for two test points
xs_test = [10.1, 10.2]

m_train = linear_reg(xs_test, fill(missing, length(xs_test)))

predictions = DynamicPPL.AbstractPPL.predict(
    Random.default_rng(), m_train, β_chain
)

ys_pred = vec(mean(Array(predictions); dims=1))

# Check if predictions match expected values within tolerance
(
    isapprox(ys_pred[1], ground_truth_β * xs_test[1], atol = 0.01),
    isapprox(ys_pred[2], ground_truth_β * xs_test[2], atol = 0.01)
)

# output

(true, true)

marginalize(
    model::DynamicPPL.Model,
    marginalized_varnames::AbstractVector{<:VarName};
    varinfo::DynamicPPL.AbstractVarInfo=link(VarInfo(model), model),
    getlogprob=DynamicPPL.getlogjoint,
    method::MarginalLogDensities.AbstractMarginalizer=MarginalLogDensities.LaplaceApprox();
    kwargs...,
)

Construct a MarginalLogDensities.MarginalLogDensity object that represents the marginal log-density of the given model, after marginalizing out the variables specified in varnames.

The resulting object can be called with a vector of parameter values to compute the marginal log-density.

Keyword arguments

• varinfo: The varinfo to use for the model. By default we use a linked VarInfo, meaning that the resulting log-density function accepts parameters that have been transformed to unconstrained space.

• getlogprob: A function which specifies which kind of marginal log-density to compute. Its default value is DynamicPPL.getlogjoint which returns the marginal log-joint probability.

• method: The marginalization method; defaults to a Laplace approximation. Please see the MarginalLogDensities.jl package for other options.

• Other keyword arguments are passed to the MarginalLogDensities.MarginalLogDensity constructor.

Example

julia> using DynamicPPL, Distributions, MarginalLogDensities

julia> @model function demo()
           x ~ Normal(1.0)
           y ~ Normal(2.0)
       end
demo (generic function with 2 methods)

julia> marginalized = marginalize(demo(), [:x]);

julia> # The resulting callable computes the marginal log-density of `y`.
       marginalized([1.0])
-1.4189385332046727

julia> logpdf(Normal(2.0), 1.0)
-1.4189385332046727

@model function f()
    x ~ Normal()
    y ~ truncated(Normal(); lower=x)
end

VarInfo(
    mld::MarginalLogDensities.MarginalLogDensity{<:LogDensityFunctionWrapper},
    unmarginalized_params::Union{AbstractVector,Nothing}=nothing
)

Retrieve the VarInfo object used in the marginalisation process.

If a Laplace approximation was used for the marginalisation, the values of the marginalized parameters are also set to their mode (note that this only happens if the mld object has been used to compute the marginal log-density at least once, so that the mode has been computed).

If a vector of unmarginalized_params is specified, the values for the corresponding parameters will also be updated in the returned VarInfo. This vector may be obtained e.g. by performing an optimization of the marginal log-density.

All other aspects of the VarInfo, such as link status, are preserved from the original VarInfo used in the marginalisation.

init_strategy = DynamicPPL.InitFromParams(varinfo.values, nothing)
oavi = DynamicPPL.OnlyAccsVarInfo((
    DynamicPPL.LogPriorAccumulator(),
    DynamicPPL.LogLikelihoodAccumulator(),
    DynamicPPL.RawValueAccumulator(false),
    # ... whatever else you need
))
_, oavi = DynamicPPL.init!!(rng, model, oavi, init_strategy, DynamicPPL.UnlinkAll())

Example

julia> using DynamicPPL, Distributions, MarginalLogDensities

julia> @model function demo()
           x ~ Normal()
           y ~ Beta(2, 2)
       end
demo (generic function with 2 methods)

julia> # Note that by default `marginalize` uses a linked VarInfo.
       mld = marginalize(demo(), [@varname(x)]);

julia> using MarginalLogDensities: Optimization, OptimizationOptimJL

julia> # Find the mode of the marginal log-density of `y`, with an initial point of `y0`.
       y0 = 2.0; opt_problem = Optimization.OptimizationProblem(mld, [y0])
OptimizationProblem. In-place: true
u0: 1-element Vector{Float64}:
 2.0

julia> # This tells us the optimal (linked) value of `y` is around 0.
       opt_solution = Optimization.solve(opt_problem, OptimizationOptimJL.NelderMead())
retcode: Success
u: 1-element Vector{Float64}:
 4.88281250001733e-5

julia> # Get the VarInfo corresponding to the mode of `y`.
       vi = VarInfo(mld, opt_solution.u);

julia> # `x` is set to its mode (which for `Normal()` is zero).
       vi[@varname(x)]
0.0

julia> # `y` is set to the optimal value we found above.
       DynamicPPL.getindex_internal(vi, @varname(y))
1-element Vector{Float64}:
 4.88281250001733e-5

julia> # To obtain values in the original constrained space, we can either
       # use `getindex`:
       vi[@varname(y)]
0.5000122070312476

julia> # Or invlink the entire VarInfo object using the model:
       vi_unlinked = DynamicPPL.invlink(vi, demo()); vi_unlinked[:]
2-element Vector{Float64}:
 0.0
 0.5000122070312476

to_submodel(model::Model[, auto_prefix::Bool])

Return a model wrapper indicating that it is a sampleable model over the return-values.

This is mainly meant to be used on the right-hand side of a ~ operator to indicate that the model can be sampled from but not necessarily evaluated for its log density.

Arguments

• model::Model: the model to wrap.
• auto_prefix::Bool: whether to automatically prefix the variables in the model using the left-hand side of the ~ statement. Default: true.

Examples

Simple example

julia> @model function demo1(x)
           x ~ Normal()
           return 1 + abs(x)
       end;

julia> @model function demo2(x, y)
            a ~ to_submodel(demo1(x))
            return y ~ Uniform(0, a)
       end;

When we sample from the model demo2(missing, 0.4) random variable x will be sampled:

julia> vi = VarInfo(demo2(missing, 0.4));

julia> @varname(a.x) in keys(vi)
true

The variable a is not tracked. However, it will be assigned the return value of demo1, and can be used in subsequent lines of the model, as shown above.

julia> @varname(a) in keys(vi)
false

We can check that the log joint probability of the model accumulated in vi is correct:

julia> x = vi[@varname(a.x)];

julia> getlogjoint(vi) ≈ logpdf(Normal(), x) + logpdf(Uniform(0, 1 + abs(x)), 0.4)
true

Without automatic prefixing

As mentioned earlier, by default, the auto_prefix argument specifies whether to automatically prefix the variables in the submodel. If auto_prefix=false, then the variables in the submodel will not be prefixed.

julia> @model function demo1(x)
           x ~ Normal()
           return 1 + abs(x)
       end;

julia> @model function demo2_no_prefix(x, z)
            a ~ to_submodel(demo1(x), false)
            return z ~ Uniform(-a, 1)
       end;

julia> vi = VarInfo(demo2_no_prefix(missing, 0.4));

julia> @varname(x) in keys(vi)  # here we just use `x` instead of `a.x`
true

However, not using prefixing is generally not recommended as it can lead to variable name clashes unless one is careful. For example, if we're re-using the same model twice in a model, not using prefixing will lead to variable name clashes: However, one can manually prefix using the prefix(::Model, input):

julia> @model function demo2(x, y, z)
            a ~ to_submodel(prefix(demo1(x), :sub1), false)
            b ~ to_submodel(prefix(demo1(y), :sub2), false)
            return z ~ Uniform(-a, b)
       end;

julia> vi = VarInfo(demo2(missing, missing, 0.4));

julia> @varname(sub1.x) in keys(vi)
true

julia> @varname(sub2.x) in keys(vi)
true

Variables a and b are not tracked, but are assigned the return values of the respective calls to demo1:

julia> @varname(a) in keys(vi)
false

julia> @varname(b) in keys(vi)
false

We can check that the log joint probability of the model accumulated in vi is correct:

julia> sub1_x = vi[@varname(sub1.x)];

julia> sub2_x = vi[@varname(sub2.x)];

julia> logprior = logpdf(Normal(), sub1_x) + logpdf(Normal(), sub2_x);

julia> loglikelihood = logpdf(Uniform(-1 - abs(sub1_x), 1 + abs(sub2_x)), 0.4);

julia> getlogjoint(vi) ≈ logprior + loglikelihood
true

prefix(ctx::AbstractContext, vn::VarName)

Apply the prefixes in the context ctx to the variable name vn.

prefix(model::Model, x::VarName)
prefix(model::Model, x::Val{sym})
prefix(model::Model, x::Any)

Return model but with all random variables prefixed by x, where x is either:

• a VarName (e.g. @varname(a)),
• a Val{sym} (e.g. Val(:a)), or
• for any other type, x is converted to a Symbol and then to a VarName. Note that this will introduce runtime overheads so is not recommended unless absolutely necessary.

Examples

julia> using DynamicPPL: prefix

julia> @model demo() = x ~ Dirac(1)
demo (generic function with 2 methods)

julia> rand(prefix(demo(), @varname(my_prefix)))
VarNamedTuple
└─ my_prefix => VarNamedTuple
                └─ x => 1

julia> rand(prefix(demo(), Val(:my_prefix)))
VarNamedTuple
└─ my_prefix => VarNamedTuple
                └─ x => 1

typed_identity(x)

Identity function, but with an overload for with_logabsdet_jacobian to ensure that it returns a sensible zero logjac.

The problem with plain old identity is that the default definition of with_logabsdet_jacobian for identity returns zero(eltype(x)): https://github.com/JuliaMath/ChangesOfVariables.jl/blob/d6a8115fc9b9419decbdb48e2c56ec9675b4c6a4/src/with_ladj.jl#L154

This is fine for most samples x, but if eltype(x) doesn't return a sensible type (e.g. if it's Any), then using identity will error with zero(Any). This can happen with, for example, ProductNamedTupleDistribution:

julia> using Distributions; d = product_distribution((a = Normal(), b = LKJCholesky(3, 0.5)));

julia> eltype(rand(d))
Any

The same problem precludes us from eventually broadening the scope of DynamicPPL.jl to support distributions with non-numeric samples.

Furthermore, in principle, the type of the log-probability should be separate from the type of the sample. Thus, instead of using zero(LogProbType), we should use the eltype of the LogJacobianAccumulator. There's no easy way to thread that through here, but if a way to do this is discovered, then typed_identity is what will allow us to obtain that custom behaviour.

@addlogprob!(ex)

Add a term to the log joint.

If ex evaluates to a NamedTuple with keys :loglikelihood and/or :logprior, the values are added to the log likelihood and log prior respectively.

If ex evaluates to a number it is added to the log likelihood.

Examples

julia> mylogjoint(x, μ) = (; loglikelihood=loglikelihood(Normal(μ, 1), x), logprior=1.0);

julia> @model function demo(x)
           μ ~ Normal()
           @addlogprob! mylogjoint(x, μ)
       end;

julia> x = [1.3, -2.1];

julia> loglikelihood(demo(x), (μ=0.2,)) ≈ mylogjoint(x, 0.2).loglikelihood
true

julia> logprior(demo(x), (μ=0.2,)) ≈ logpdf(Normal(), 0.2) + mylogjoint(x, 0.2).logprior
true

and to reject samples:

julia> @model function demo(x)
           m ~ MvNormal(zero(x), I)
           if dot(m, x) < 0
               @addlogprob! (; loglikelihood=-Inf)
               # Exit the model evaluation early
               return
           end
           x ~ MvNormal(m, I)
           return
       end;

julia> logjoint(demo([-2.1]), (m=[0.2],)) == -Inf
true

returned(model::Model, chain::MCMCChains.Chains)

Execute model for each of the samples in chain and return an array of the values returned by the model for each sample.

Examples

General

Often you might have additional quantities computed inside the model that you want to inspect, e.g.

@model function demo(x)
    # sample and observe
    θ ~ Prior()
    x ~ Likelihood()
    return interesting_quantity(θ, x)
end
m = demo(data)
chain = sample(m, alg, n)
# To inspect the `interesting_quantity(θ, x)` where `θ` is replaced by samples
# from the posterior/`chain`:
returned(m, chain) # <= results in a `Vector` of returned values
                               #    from `interesting_quantity(θ, x)`

Concrete (and simple)

julia> using Turing

julia> @model function demo(xs)
           s ~ InverseGamma(2, 3)
           m_shifted ~ Normal(10, √s)
           m = m_shifted - 10

           for i in eachindex(xs)
               xs[i] ~ Normal(m, √s)
           end

           return (m, )
       end
demo (generic function with 1 method)

julia> model = demo(randn(10));

julia> chain = sample(model, MH(), 10);

julia> returned(model, chain)
10×1 Array{Tuple{Float64},2}:
 (2.1964758025119338,)
 (2.1964758025119338,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.09270081916291417,)
 (0.043088571494005024,)
 (-0.16489786710222099,)
 (-0.16489786710222099,)

returned(model::Model, parameters...)

Initialise a model using the given parameters and return the model's return value. The parameters must be provided in a format that can be wrapped in an InitFromParams, i.e., InitFromParams(parameters..., nothing) must be a valid AbstractInitStrategy (where nothing is the fallback strategy to use if parameters are not provided).

As far as DynamicPPL is concerned, parameters can be either a singular NamedTuple or an AbstractDict{<:VarName}; however this method is left flexible to allow for other packages that wish to extend InitFromParams.

Example

julia> using DynamicPPL, Distributions

julia> @model function demo()
           m ~ Normal()
           return (mp1 = m + 1,)
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> returned(model, (; m = 1.0))
(mp1 = 2.0,)

julia> returned(model, Dict{VarName,Float64}(@varname(m) => 2.0))
(mp1 = 3.0,)

DynamicPPL.pointwise_logdensities(
    model::DynamicPPL.Model,
    chain::MCMCChains.Chains,
    ::Type{Tout}=MCMCChains.Chains
    ::Val{whichlogprob}=Val(:both),
)

Runs model on each sample in chain, returning a new MCMCChains.Chains object where the log-density of each variable at each sample is stored (rather than its value).

whichlogprob specifies which log-probabilities to compute. It can be :both, :prior, or :likelihood.

You can pass Tout=OrderedDict to get the result as an OrderedDict{VarName, Matrix{Float64}} instead.

See also: DynamicPPL.pointwise_loglikelihoods, DynamicPPL.pointwise_prior_logdensities.

Examples

julia> using MCMCChains

julia> @model function demo(xs, y)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, √s)
           for i in eachindex(xs)
               xs[i] ~ Normal(m, √s)
           end
           y ~ Normal(m, √s)
       end
demo (generic function with 2 methods)

julia> # Example observations.
       model = demo([1.0, 2.0, 3.0], [4.0]);

julia> # A chain with 3 iterations.
       chain = Chains(
           reshape(1.:6., 3, 2),
           [:s, :m];
           info=(varname_to_symbol=Dict(
               @varname(s) => :s,
               @varname(m) => :m,
           ),),
       );

julia> plds = pointwise_logdensities(model, chain)
Chains MCMC chain (3×6×1 Array{Float64, 3}):

Iterations        = 1:1:3
Number of chains  = 1
Samples per chain = 3
parameters        = s, m, xs[1], xs[2], xs[3], y
[...]

julia> plds[:s]
2-dimensional AxisArray{Float64,2,...} with axes:
    :iter, 1:1:3
    :chain, 1:1
And data, a 3×1 Matrix{Float64}:
 -0.8027754226637804
 -1.3822169643436162
 -2.0986122886681096

julia> # The above is the same as:
       logpdf.(InverseGamma(2, 3), chain[:s])
3×1 Matrix{Float64}:
 -0.8027754226637804
 -1.3822169643436162
 -2.0986122886681096

julia> # Alternatively: pldsdict = pointwiselogdensities(model, chain, OrderedDict) OrderedDict{VarName, Matrix{Float64}} with 6 entries: s => [-0.802775; -1.38222; -2.09861;;] m => [-8.91894; -7.51551; -7.46824;;] xs[1] => [-5.41894; -5.26551; -5.63491;;] xs[2] => [-2.91894; -3.51551; -4.13491;;] xs[3] => [-1.41894; -2.26551; -2.96824;;] y => [-0.918939; -1.51551; -2.13491;;]

DynamicPPL.pointwise_loglikelihoods(
    model::DynamicPPL.Model,
    chain::MCMCChains.Chains,
    ::Type{Tout}=MCMCChains.Chains
)

Compute the pointwise log-likelihoods of the model given the chain. This is the same as pointwise_logdensities(model, chain), but only including the likelihood terms.

See also: DynamicPPL.pointwise_logdensities, DynamicPPL.pointwise_prior_logdensities.

DynamicPPL.pointwise_prior_logdensities(
    model::DynamicPPL.Model,
    chain::MCMCChains.Chains
)

Compute the pointwise log-prior-densities of the model given the chain. This is the same as pointwise_logdensities(model, chain), but only including the prior terms.

See also: DynamicPPL.pointwise_logdensities, DynamicPPL.pointwise_loglikelihoods.

extract_priors([rng::Random.AbstractRNG, ]model::Model)

Extract the priors from a model. This is done by sampling from the model and recording the distributions that are used to generate the samples.

Examples

Changing parameters

julia> using Distributions, StableRNGs

julia> rng = StableRNG(42);

julia> @model function model_dynamic_parameters()
           x ~ Normal(0, 1)
           y ~ Normal(x, 1)
       end;

julia> model = model_dynamic_parameters();

julia> extract_priors(rng, model)[@varname(y)]
Normal{Float64}(μ=-0.6702516921145671, σ=1.0)

julia> extract_priors(rng, model)[@varname(y)]
Normal{Float64}(μ=1.3736306979834252, σ=1.0)

Changing support

julia> using LinearAlgebra, Distributions, StableRNGs

julia> rng = StableRNG(42);

julia> @model function model_dynamic_support()
           n ~ Categorical(ones(10) ./ 10)
           x ~ MvNormal(zeros(n), I)
       end;

julia> model = model_dynamic_support();

julia> length(extract_priors(rng, model)[@varname(x)])
6

julia> length(extract_priors(rng, model)[@varname(x)])
9

extract_priors(model::Model, varinfo::AbstractVarInfo)

Extract the priors that were used to generate a VarInfo.

This is done by evaluating the model at the values present in varinfo and recording the distributions that are present at each tilde statement.

A named distribution that carries the name of the random variable with it.

run_ad(
    model::Model,
    adtype::ADTypes.AbstractADType;
    test::Union{AbstractADCorrectnessTestSetting,Bool}=WithBackend(),
    benchmark=false,
    atol::AbstractFloat=1e-8,
    rtol::AbstractFloat=sqrt(eps()),
    getlogdensity::Function=getlogjoint_internal,
    rng::Random.AbstractRNG=Random.default_rng(),
    varinfo::AbstractVarInfo=link(VarInfo(model), model),
    params::Union{Nothing,Vector{<:AbstractFloat}}=nothing,
    verbose=true,
)::ADResult

Description

Test the correctness and/or benchmark the AD backend adtype for the model model.

Whether to test and benchmark is controlled by the test and benchmark keyword arguments. By default, test is true and benchmark is false.

Note that to run AD successfully you will need to import the AD backend itself. For example, to test with AutoReverseDiff() you will need to run import ReverseDiff.

Arguments

There are two positional arguments, which absolutely must be provided:

1. model - The model being tested.
2. adtype - The AD backend being tested.

Everything else is optional, and can be categorised into several groups:

1. How to specify the VarInfo.

   DynamicPPL contains several different types of VarInfo objects which change the way model evaluation occurs. If you want to use a specific type of VarInfo, pass it as the varinfo argument. Otherwise, it will default to using a linked TypedVarInfo generated from the model. Here, linked means that the parameters in the VarInfo have been transformed to unconstrained Euclidean space if they aren't already in that space.

2. How to specify the parameters.

   For maximum control over this, generate a vector of parameters yourself and pass this as the params argument. If you don't specify this, it will be taken from the contents of the VarInfo.

   Note that if the VarInfo is not specified (and thus automatically generated) the parameters in it will have been sampled from the prior of the model. If you want to seed the parameter generation for the VarInfo, you can pass the rng keyword argument, which will then be used to create the VarInfo.

   Finally, note that these only reflect the parameters used for evaluating the gradient. If you also want to control the parameters used for preparing the gradient, then you need to manually set these parameters in the VarInfo object, for example using vi = DynamicPPL.unflatten!!(vi, prep_params). You could then evaluate the gradient at a different set of parameters using the params keyword argument.

3. Which type of logp is being calculated.

   By default, run_ad evaluates the 'internal log joint density' of the model, i.e., the log joint density in the unconstrained space. Thus, for example, in

   @model f() = x ~ LogNormal()

   the internal log joint density is logpdf(Normal(), log(x)). This is the relevant log density for e.g. Hamiltonian Monte Carlo samplers and is therefore the most useful to test.

   If you want the log joint density in the original model parameterisation, you can use getlogjoint. Likewise, if you want only the prior or likelihood, you can use getlogprior or getloglikelihood, respectively.

4. How to specify the results to compare against.

   Once logp and its gradient has been calculated with the specified adtype, it can optionally be tested for correctness. The exact way this is tested is specified in the test parameter.

   There are several options for this:

   • You can explicitly specify the correct value using WithExpectedResult().
   • You can compare against the result obtained with a different AD backend using WithBackend(adtype).
   • You can disable testing by passing NoTest().
   • The default is to compare against the result obtained with ForwardDiff, i.e. WithBackend(AutoForwardDiff()).
   • test=false and test=true are synonyms for NoTest() and WithBackend(AutoForwardDiff()), respectively.

5. How to specify the tolerances. (Only if testing is enabled.)

   Both absolute and relative tolerances can be specified using the atol and rtol keyword arguments respectively. The behaviour of these is similar to isapprox(), i.e. the value and gradient are considered correct if either atol or rtol is satisfied. The default values are 100*eps() for atol and sqrt(eps()) for rtol.

   For the most part, it is the rtol check that is more meaningful, because we cannot know the magnitude of logp and its gradient a priori. The atol value is supplied to handle the case where gradients are equal to zero.

6. Whether to benchmark.

   By default, benchmarking is disabled. To enable it, set benchmark=true. When enabled, the time taken to evaluate logp as well as its gradient is measured using Chairmarks.jl, and the ADResult object returned will contain grad_time and primal_time fields with the median times (in seconds).

7. Whether to output extra logging information.

   By default, this function prints messages when it runs. To silence it, set verbose=false.

Returns / Throws

Returns an ADResult object, which contains the results of the test and/or benchmark.

If test is true and the AD backend returns an incorrect value or gradient, an ADIncorrectException is thrown. If a different error occurs, it will be thrown as-is.

AbstractADCorrectnessTestSetting

Different ways of testing the correctness of an AD backend.

WithBackend(adtype::AbstractADType=AutoForwardDiff()) <: AbstractADCorrectnessTestSetting

Test correctness by comparing it against the result obtained with adtype.

adtype defaults to ForwardDiff.jl, since it's the default AD backend used in Turing.jl.

WithExpectedResult(
    value::T,
    grad::AbstractVector{T}
) where {T <: AbstractFloat}
<: AbstractADCorrectnessTestSetting

Test correctness by comparing it against a known result (e.g. one obtained analytically, or one obtained with a different backend previously). Both the value of the primal (i.e. the log-density) as well as its gradient must be supplied.

NoTest() <: AbstractADCorrectnessTestSetting

Disable correctness testing.

ADResult{Tparams<:AbstractFloat,Tresult<:AbstractFloat,Ttol<:AbstractFloat}

Data structure to store the results of the AD correctness test.

The type parameter Tparams is the numeric type of the parameters passed in; Tresult is the type of the value and the gradient; and Ttol is the type of the absolute and relative tolerances used for correctness testing.

Fields

• model::Model: The DynamicPPL model that was tested

• getlogdensity::Function: The function used to extract the log density from the model

• varinfo::AbstractVarInfo: The VarInfo that was used

• params::Vector{Tparams} where Tparams<:AbstractFloat: The values at which the model was evaluated

• adtype::ADTypes.AbstractADType: The AD backend that was tested

• atol::AbstractFloat: Absolute tolerance used for correctness test

• rtol::AbstractFloat: Relative tolerance used for correctness test

• value_expected::Union{Nothing, Tresult} where Tresult<:AbstractFloat: The expected value of logp

• grad_expected::Union{Nothing, Vector{Tresult}} where Tresult<:AbstractFloat: The expected gradient of logp

• value_actual::AbstractFloat: The value of logp (calculated using adtype)

• grad_actual::Vector{Tresult} where Tresult<:AbstractFloat: The gradient of logp (calculated using adtype)

• grad_time::Union{Nothing, Tresult} where Tresult<:AbstractFloat: If benchmarking was requested, the time taken by the AD backend to evaluate the gradient of logp (in seconds)

• primal_time::Union{Nothing, Tresult} where Tresult<:AbstractFloat: If benchmarking was requested, the time taken by the AD backend to evaluate logp (in seconds)

ADIncorrectException{T<:AbstractFloat}

Exception thrown when an AD backend returns an incorrect value or gradient.

The type parameter T is the numeric type of the value and gradient.

Fields

• value_expected::AbstractFloat

• value_actual::AbstractFloat

• grad_expected::Vector{T} where T<:AbstractFloat

• grad_actual::Vector{T} where T<:AbstractFloat

• atol::AbstractFloat

• rtol::AbstractFloat

A collection of models corresponding to the posterior distribution defined by the generative process

s ~ InverseGamma(2, 3)
m ~ Normal(0, √s)
1.5 ~ Normal(m, √s)
2.0 ~ Normal(m, √s)

or by

s[1] ~ InverseGamma(2, 3)
s[2] ~ InverseGamma(2, 3)
m[1] ~ Normal(0, √s)
m[2] ~ Normal(0, √s)
1.5 ~ Normal(m[1], √s[1])
2.0 ~ Normal(m[2], √s[2])

These are examples of a Normal-InverseGamma conjugate prior with Normal likelihood, for which the posterior is known in closed form.

In particular, for the univariate model (the former one):

mean(s) == 49 / 24
mean(m) == 7 / 6

And for the multivariate one (the latter one):

mean(s[1]) == 19 / 8
mean(m[1]) == 3 / 4
mean(s[2]) == 8 / 3
mean(m[2]) == 1

A tuple of all models defined in DynamicPPL.TestUtils.

logprior_true(model, args...)

Return the logprior of model for args.

This should generally be implemented by hand for every specific model.

See also: logjoint_true, loglikelihood_true.

loglikelihood_true(model, args...)

Return the loglikelihood of model for args.

This should generally be implemented by hand for every specific model.

See also: logjoint_true, logprior_true.

logjoint_true(model, args...)

Return the logjoint of model for args.

Defaults to logprior_true(model, args...) + loglikelihood_true(model, args..).

This should generally be implemented by hand for every specific model so that the returned value can be used as a ground-truth for testing things like:

1. Validity of evaluation of model using a particular implementation of AbstractVarInfo.
2. Validity of a sampler when combined with DynamicPPL by running the sampler twice: once targeting ground-truth functions, e.g. logjoint_true, and once targeting model.

And more.

See also: logprior_true, loglikelihood_true.

logprior_true_with_logabsdet_jacobian(model::Model, args...)

Return a tuple (args_unconstrained, logprior_unconstrained) of model for args....

Unlike logprior_true, the returned logprior computation includes the log-absdet-jacobian adjustment, thus computing logprior for the unconstrained variables.

Note that args are assumed be in the support of model, while args_unconstrained will be unconstrained.

See also: logprior_true.

logjoint_true_with_logabsdet_jacobian(model::Model, args...)

Return a tuple (args_unconstrained, logjoint) of model for args.

Unlike logjoint_true, the returned logjoint computation includes the log-absdet-jacobian adjustment, thus computing logjoint for the unconstrained variables.

Note that args are assumed be in the support of model, while args_unconstrained will be unconstrained.

This should generally not be implemented directly, instead one should implement logprior_true_with_logabsdet_jacobian for a given model.

See also: logjoint_true, logprior_true_with_logabsdet_jacobian.

varnames(model::Model)

Return the VarNames defined in model, as a Vector.

posterior_mean(model::Model)

Return a NamedTuple compatible with varnames(model) where the values represent the posterior mean under model.

"Compatible" means that a varname from varnames(model) can be used to extract the corresponding value using e.g. AbstractPPL.getvalue(posterior_mean(model), varname).

setup_varinfos(model::Model, example_values::NamedTuple; include_threadsafe::Bool=false)

Return a tuple of instances for different implementations of AbstractVarInfo with each vi, supposedly, satisfying vi[vn] == get(example_values, vn) for vn in varnames.

If include_threadsafe is true, then the returned tuple will also include thread-safe versions of the varinfo instances.

test_values(vi::AbstractVarInfo, vals::NamedTuple, vns)

Test that vi[vn] corresponds to the correct value in vals for every vn in vns.

DynamicPPL.DebugUtils.check_model(
    [rng::Random.AbstractRNG,]
    model::Model;
    error_on_failure=false,
    fail_if_discrete=false
)

Check model for potential issues. Returns true if the model check succeeded, false otherwise.

The model is only evaluated a single time, so if the model contains any indeterminism, results may differ across runs. The rng argument can be used to control reproducibility if needed.

Issues that this function checks for

• Repeated usage of the same or overlapping VarNames

• NaN on the left-hand side of observe statements

• (if fail_if_discrete is set) Usage of discrete distributions

• Empty models emit a warning, but do not fail (since they are not incorrect per se)

Keyword arguments

• error_on_failure::Bool: Whether to throw an error (instead of just returning false) if the model check fails.

• fail_if_discrete::Bool: Whether to fail (i.e., return false or throw an error, depending on error_on_failure) when the model contains discrete distributions. Discrete distributions do not have a differentiable log-density and are incompatible with gradient-based approaches such as HMC / NUTS or optimisation.

Examples

Correct model

julia> using DynamicPPL.DebugUtils: check_model; using Distributions

julia> @model demo_correct() = x ~ Normal()
demo_correct (generic function with 2 methods)

julia> model = demo_correct();

julia> check_model(model)
true

julia> cond_model = model | (x = 1.0,);

julia> # Empty models will issue a warning, but not a failure
       check_model(cond_model)
┌ Warning: The model does not contain any parameters.
└ @ DynamicPPL.DebugUtils DynamicPPL.jl/src/debug_utils.jl:215
true

Incorrect model

```jldoctest; setup=:(using Distributions) julia> using DynamicPPL.DebugUtils: check_model; using Distributions

julia> @model function demoincorrect() # Sampling x twice. x ~ Normal() x ~ Exponential() end demoincorrect (generic function with 2 methods)

julia> # Notice that VarInfo(modelincorrect) evaluates the model, but doesn't actually # alert us to the issue of x being sampled twice. model = demoincorrect(); varinfo = VarInfo(model);

julia> checkmodel(model; erroron_failure=true) ERROR: varname x used multiple times in model

has_static_constraints([rng, ]model::Model; num_evals=5)

Attempts to detect whether model has static constraints (i.e., the support of all variables is the same regardless of what their values are). Returns true if the model has static constraints, false otherwise.

Note that this is a heuristic check based on sampling from the model multiple times and checking if the model is consistent across runs.

Arguments

• rng::Random.AbstractRNG: The random number generator to use when evaluating the model.
• model::Model: The model to check.

Keyword Arguments

• num_evals::Int: The number of evaluations to perform. Default: 5.

model_warntype(model[, varinfo]; optimize=true)

Check the type stability of the model's evaluator, warning about any potential issues.

This simply calls @code_warntype on the model's evaluator, filling in internal arguments where needed.

Arguments

• model::Model: The model to check.
• varinfo::AbstractVarInfo: The varinfo to use when evaluating the model. Default: VarInfo(model).

Keyword Arguments

• optimize::Bool: Whether to generate optimized code. Default: false.

model_typed(model[, varinfo]; optimize=true)

Return the type inference for the model's evaluator.

This simply calls @code_typed on the model's evaluator, filling in internal arguments where needed.

Arguments

• model::Model: The model to check.
• varinfo::AbstractVarInfo: The varinfo to use when evaluating the model. Default: VarInfo(model).

Keyword Arguments

• optimize::Bool: Whether to generate optimized code. Default: true.

gen_evaluator_call_with_types(model[, varinfo])

Generate the evaluator call and the types of the arguments.

Arguments

• model::Model: The model whose evaluator is of interest.
• varinfo::AbstractVarInfo: The varinfo to use when evaluating the model. Default: VarInfo(model).

Returns

A 2-tuple with the following elements:

• f: This is either model.f or Core.kwcall, depending on whether the model has keyword arguments.
• argtypes::Type{<:Tuple}: The types of the arguments for the evaluator.

AbstractVarInfo

Abstract supertype for data structures that capture random variables when executing a probabilistic model and accumulate log densities such as the log likelihood or the log joint probability of the model.

See also: VarInfo

VarInfo{
    Tfm<:AbstractTransformStrategy,
    T<:VarNamedTuple,
    Accs<:AccumulatorTuple
} <: AbstractVarInfo

The default implementation of AbstractVarInfo, storing variable values and accumulators.

VarInfo is quite a thin wrapper around a VarNamedTuple storing the variable values, and a tuple of accumulators. The only really noteworthy thing about it is that it stores the values of variables vectorised as instances of AbstractTransformedValue. That is, it stores each value as a special vector with a flag indicating whether it is just a vectorised value (VectorValue), or whether it is also linked (LinkedVectorValue). It also stores the size of the actual post-transformation value. These are all accessible via AbstractTransformedValue.

VarInfo additionally stores a transform strategy, which reflects the linked status of variables inside the VarInfo. For example, a VarInfo{LinkAll} should contain only LinkedVectorValues in its values field.

Because the job of VarInfo is to store transformed values, there is no generic setindex!! implementation on VarInfo itself. Instead, all storage must go via setindex_with_dist!!, which takes care of storing the value in the correct transformed form. This in turn means that the distribution on the right-hand side of a tilde-statement must be available when modifying a VarInfo.

You can use getindex on VarInfo to obtain values in the support of the original distribution. To directly get access to the internal vectorised values, use getindex_internal, setindex_internal!!, and unflatten!!.

For more details on the internal storage, see documentation of AbstractTransformedValue and VarNamedTuple.

Fields

• transform_strategy::AbstractTransformStrategy

• values::VarNamedTuple

• accs::DynamicPPL.AccumulatorTuple

setindex_with_dist!!(
    vi::VarInfo,
    tval::Union{VectorValue,LinkedVectorValue},
    dist::Distribution,
    vn::VarName,
    template::Any,
)

Set the value of vn in vi to tval. Note that this will cause the linked status of vi to update according to what tval is. That means that whether or not a variable is considered to be 'linked' is determined by tval rather than the previous status of vi.

setindex_with_dist!!(
    vi::VarInfo,
    tval::UntransformedValue,
    dist::Distribution,
    vn::VarName,
    template::Any
)

Vectorise tval (into a VectorValue) and store it. (Note that if setindex_with_dist!! receives an UntransformedValue, the variable is always considered unlinked, since if it were to be linked, apply_transform_strategy will already have done so.)

is_transformed(vi::AbstractVarInfo[, vns::Union{VarName, AbstractVector{<:Varname}}])

Return true if vi is working in unconstrained space, and false if vi is assuming realizations to be in support of the corresponding distributions.

If vns is provided, then only check if this/these varname(s) are transformed.

set_transformed!!(vi::AbstractVarInfo, trans::Bool[, vn::VarName])

Return vi with is_transformed(vi, vn) evaluating to true.

If vn is not specified, then is_transformed(vi) evaluates to true for all variables.

VarNamedTuple{names,Values}

A NamedTuple-like structure with VarName keys.

VarNamedTuple is a data structure for storing arbitrary data, keyed by VarNames, in an efficient and type stable manner. It is mainly used through getindex, setindex!!, templated_setindex!!, and haskey, all of which only accept VarNames as keys. Other notable methods are merge and subset.

VarNamedTuple has an ordering to its elements, and two VarNamedTuples with the same keys and values but in different orders are considered different for equality and hashing. Iterations such as keys and values respect this ordering. The ordering is dependent on the order in which elements were inserted into the VarNamedTuple, though isn't always equal to it. More specifically

• Any new keys that have a joint parent VarName with an existing key are inserted after that key. For instance, if one first inserts, in order, @varname(a.x), @varname(b), and @varname(a.y), the resulting order will be (@varname(a.x), @varname(a.y), @varname(b)).
• Index keys, like@varname(a[3])or@varname(b[2,3,4:5]), are always iterated in the same order anArraywith the same indices would be iterated. For instance, if one first inserts, in order,@varname(a[2]),@varname(b), and@varname(a[1]), the resulting order will be(@varname(a[1]), @varname(a[2]), @varname(b))`.

Otherwise insertion order is respected.

setindex!! and getindex on VarNamedTuple are type stable as long as one does not store heterogeneous data under different indices of the same symbol. That is, if either

• one sets a[1] and a[2] to be of different types, or
• if a[1] and a[2] both exist, one sets a[1].b without setting a[2].b,

then getting values for a[1] or a[2] will not be type stable.

@vnt begin ... end

Construct a VarNamedTuple from a block of assignments. Each assignment should be of the form var := value, where var is a variable name. This is best illustrated by example:

julia> using DynamicPPL

julia> @vnt begin
           a := 1
           b := 2
       end
VarNamedTuple
├─ a => 1
└─ b => 2

You can set entirely arbitrary variables:

julia> @vnt begin
           a.b.c.d.e := "hello"
       end
VarNamedTuple
└─ a => VarNamedTuple
        └─ b => VarNamedTuple
                └─ c => VarNamedTuple
                        └─ d => VarNamedTuple
                                └─ e => "hello"

For variables that have indexing, it is often necessary to provide a template, so that the VNT 'knows' what kind of array is being used to store the values, and can set the values in the appropriate places (consider e.g. OffsetArrays where x[1] may not mean what it usually does).

This is done by inserting a @template macro call in the block. The @template macro accepts whitespace-separated arguments, which must either be

• a single symbol (e.g. @template x), in which case the template is the value of x (and x must already be defined in the current scope); or
• an assignment of the form y = expr, in which case the template for y is the value of expr. In this case y does not need to be defined in the current scope, but any symbols referenced in expr must be.

For example:

julia> x = zeros(5); outside_y = zeros(3, 3);

julia> @vnt begin
            @template x y=outside_y
            x[1] := 1.0
            y[1, 1] := 2.0
       end
VarNamedTuple
├─ x => PartialArray size=(5,) data::Vector{Float64}
│       └─ (1,) => 1.0
└─ y => PartialArray size=(3, 3) data::Matrix{Float64}
        └─ (1, 1) => 2.0

If no template is provided, the VNT will use a GrowableArray. This can produce correct results in simple cases, but is not recommended for general use. Please see the VarNamedTuple documentation for more details.

julia> @vnt begin
            # No template provided.
            x[1] := 1.0
            y[1, 1] := 2.0
       end
┌ Warning: Creating a growable `Base.Array` of dimension 1 to store values. This may not match the actual type or size of the actual `AbstractArray` that will be used inside the DynamicPPL model.
│
│  If this is not the type or size that you expect, please see: https://turinglang.org/docs/uri/growablearray
└ @ DynamicPPL.VarNamedTuples /path/to/DynamicPPL.jl/src/varnamedtuple/partial_array.jl:823
┌ Warning: Creating a growable `Base.Array` of dimension 2 to store values. This may not match the actual type or size of the actual `AbstractArray` that will be used inside the DynamicPPL model.
│
│  If this is not the type or size that you expect, please see: https://turinglang.org/docs/uri/growablearray
└ @ DynamicPPL.VarNamedTuples /path/to/DynamicPPL.jl/src/varnamedtuple/partial_array.jl:823
VarNamedTuple
├─ x => PartialArray size=(1,) data::DynamicPPL.VarNamedTuples.GrowableArray{Float64, 1}
│       └─ (1,) => 1.0
└─ y => PartialArray size=(1, 1) data::DynamicPPL.VarNamedTuples.GrowableArray{Float64, 2}
        └─ (1, 1) => 2.0

vnt_size(x)

Get the size of an object x for use in VarNamedTuple and PartialArray.

By default, this falls back onto Base.size, but can be overloaded for custom types. This notion of type is used to determine whether a value can be set into a PartialArray as a block, see the docstring of PartialArray and ArrayLikeBlock for details.

apply!!(func, vnt::VarNamedTuple, name::VarName)

Apply func to the subdata at name in vnt, and set the result back at name.

Like map_values!!, but only for a single VarName.

julia> using DynamicPPL: VarNamedTuple, setindex!!

julia> using DynamicPPL.VarNamedTuples: apply!!

julia> vnt = VarNamedTuple()
VarNamedTuple()

julia> vnt = setindex!!(vnt, [1, 2, 3], @varname(a))
VarNamedTuple
└─ a => [1, 2, 3]

julia> apply!!(x -> x .+ 1, vnt, @varname(a))
VarNamedTuple
└─ a => [2, 3, 4]

densify!!(vnt::VarNamedTuple)

Convert any PartialArrays in vnt where all elements are filled to just normal arrays.

map_pairs!!(func, vnt::VarNamedTuple)

Apply func to all key => value pairs of vnt, in place if possible.

func should accept a pair of VarName and value, and return the new value to be set.

map_values!!(func, vnt::VarNamedTuple)

Apply func to elements of vnt, in place if possible.

PartialArray{
    ElType,
    num_dims,
    D<:AbstractArray{ElType,num_dims},
    M<:AbstractArray{Bool,num_dims}
}

An array-like like structure that may only have some of its elements defined.

One can set values in a PartialArray either element-by-element, or with ranges like arr[1:3,2] = [5,10,15]. When setting values over a range of indices, the value being set must either be an AbstractArray or otherwise something for which vnt_size(value) or Base.size(value) (which vnt_size falls back onto) is defined, and the size matches the range. If the value is an AbstractArray, the elements are copied individually, but if it is not, the value is stored as a block, that takes up the whole range, e.g. [1:3,2], but is only a single object. Getting such a block-value must be done with the exact same range of indices, otherwise an error is thrown.

If the element type of a PartialArray is not concrete, any call to setindex!! will check if, after the new value has been set, the element type can be made more concrete. If so, a new PartialArray with a more concrete element type is returned. Thus the element type of any PartialArray should always be as concrete as is allowed by the elements in it.

The internal implementation of an PartialArray consists of two arrays: one holding the data and the other one being a boolean mask indicating which elements are defined. These internal arrays may need resizing when new elements are set that have index ranges larger than the current internal arrays. To avoid resizing too often, the internal arrays are resized in exponentially increasing steps. This means that most setindex!! calls are very fast, but some may incur substantial overhead due to resizing and copying data. It also means that the largest index set so far determines the memory usage of the PartialArray. PartialArrays are thus well-suited when most values in it will eventually be set. If only a few scattered values are set, a structure like SparseArray may be more appropriate.

DynamicPPL.VarNamedTuples.templated_setindex!!(vnt, value, vn, template; allow_new=Val(true))

Assign value to the location in vnt specified by vn.

The argument template must be provided in order to guide the creation of PartialArrays, as well as to concretise any dynamic indices in vn. It must be an object that has the shape of the top-level symbol in vn. For example:

vnt = VarNamedTuple()
templated_setindex!!(vnt, 10, @varname(x[1]), rand(2, 2))

Here, rand(2, 2) is the template for the top-level symbol x, which tells setindex!! that x should be a PartialArray that is backed by a matrix.

The actual data inside template is not needed, and template is never mutated by this call.

NoTemplate

A singleton struct representing the fact that there is no template for a top-level variable. When NoTemplate is used, several things happen:

• If you attempt to call make_leaf with an Index optic, this will error.
• When recursing into substructures, NoTemplate will be propagated.

Collectively this means that you can only set values for variables that only have Property optics, unless a template is provided.

It might seem more idiomatic to use nothing or missing for this. However, this causes a bug with BangBang.setindex!!: https://github.com/JuliaFolds2/BangBang.jl/issues/43 so we use a dedicated struct instead.

SkipTemplate{N}(value)

A struct representing the fact that value is the template for the variable N levels down from the top-level variable. In other words, SkipTemplate{0}(value) is equivalent to just value, and SkipTemplate{1}(value) means that value is the template for a when setting the variable @varname(x.a).

NamedTuple(vnt::VarNamedTuple)

Convert a VarNamedTuple to a standard NamedTuple, provided all keys in the VarNamedTuple are VarNames with top-level symbols. If any key is a VarName with a non-identity optic (e.g., @varname(x.a) or @varname(x[1])), this will throw an ArgumentError.

Examples

julia> using DynamicPPL, BangBang

julia> vnt = VarNamedTuple(); vnt = setindex!!(vnt, 10, @varname(x))
VarNamedTuple
└─ x => 10

julia> NamedTuple(vnt)
(x = 10,)

julia> vnt2 = setindex!!(vnt, 20, @varname(y.a))
VarNamedTuple
├─ x => 10
└─ y => VarNamedTuple
        └─ a => 20

julia> NamedTuple(vnt2)
ERROR: ArgumentError: Cannot convert VarNamedTuple containing non-identity VarNames to NamedTuple. To create a NamedTuple, all keys in the VarNamedTuple must be top-level symbols.
[...]

AbstractAccumulator

An abstract type for accumulators.

An accumulator is an object that may change its value at every tilde_assume!! or tilde_observe!! call based on the random variable in question. The obvious examples of accumulators are the log prior and log likelihood. Other examples might be a variable that counts the number of observations in a trace, or a list of the names of random variables seen so far.

An accumulator type T <: AbstractAccumulator must implement the following methods:

• accumulator_name(acc::T) or accumulator_name(::Type{T})
• accumulate_observe!!(acc::T, dist, val, vn)
• accumulate_assume!!(acc::T, val, tval, logjac, vn, dist, template)
• reset(acc::T)
• Base.copy(acc::T)

In these functions:

• val is the new value of the random variable sampled from a distribution (always in the original unlinked space), or the value on the left-hand side of an observe statement.
• tval is the original AbstractTransformedValue that was obtained from the initialisation strategy. This is passed through unchanged to accumulate_assume!! since it can be reused for some accumulators (e.g. when storing linked values, if the linked value was already provided, it is faster to reuse it than to re-link val).
• dist is the distribution on the RHS of the tilde statement.
• vn is the VarName that is on the left-hand side of the tilde-statement. If the tilde-statement is a literal observation like 0.0 ~ Normal(), then vn is nothing.
• logjac is the log determinant of the Jacobian of the link transformation, if the variable is stored as a linked value in the VarInfo. If the variable is stored in its original, unlinked form, then logjac is zero.
• template is a value that conveys the shape of the top-level symbol in vn, and is used specifically for accumulators that carry VarNamedTuples.

To be able to work with multi-threading, it should also implement:

• split(acc::T)
• combine(acc::T, acc2::T)

See the documentation for each of these functions for more details.

accumulate_assume!!(acc::AbstractAccumulator, val, tval, logjac, vn, right, template)

Update acc in a tilde_assume!! call. Returns the updated acc.

See AbstractAccumulator for the meaning of the arguments.

accumulate_assume!! may mutate acc, but not any of the other arguments.

See also: accumulate_observe!!

accumulate_observe!!(acc::AbstractAccumulator, right, left, vn)

Update acc in a tilde_observe!! call. Returns the updated acc.

See AbstractAccumulator for the meaning of the arguments.

accumulate_observe!! may mutate acc, but not any of the other arguments.

See also: accumulate_assume!!

accumulator_name(acc::AbstractAccumulator)

Return a Symbol which can be used as a name for acc.

The name has to be unique in the sense that a VarInfo can only have one accumulator for each name. The most typical case, and the default implementation, is that the name only depends on the type of acc, not on its value.

reset(acc::AbstractAccumulator)

Return a new accumulator like acc, but with its contents reset to the state that they should be at the beginning of model evaluation.

Note that this may in general have very similar behaviour to split, and may share the same implementation, but the difference is that split may in principle happen at any stage during model evaluation, whereas reset is only called at the beginning of model evaluation.

split(acc::AbstractAccumulator)

Return a new accumulator like acc suitable for use in a forked thread.

The returned value should be such that combine(acc, split(acc)) is equal to acc. This is used in the context of multi-threading where different threads may accumulate independently and the results are then combined.

Note that this may in general have very similar behaviour to reset, but is semantically different. See reset for more details.

See also: combine

combine(acc::AbstractAccumulator, acc2::AbstractAccumulator)

Combine two accumulators which have the same type (but may, in general, have different type parameters). Returns a new accumulator of the same type.

See also: split

VNTAccumulator{AccName}(f::F, values::VarNamedTuple=VarNamedTuple()) where {AccName,F}

A generic accumulator that applies a function f to values seen during model execution and stores the results in a VarNamedTuple.

AccName is the name of the accumulator, and is exposed to allow users to define and use multiple forms of VNTAccumulator within the same set of accumulators. In theory, each VNTAccumulator with the same function f should use the same accumulator name. This is not enforced.

The function f should have the signature:

f(val, tval, logjac, vn, dist) -> value_to_store

where val, tval, logjac, vn, and dist have their usual meanings in accumulate_assume!! (see its docstring for more details). If a value does not need to be accumulated, this can be signalled by returning DoNotAccumulate() from f.

DoNotAccumulate()

Sentinel value indicating that no accumulation should be performed for a given variable.

getacc(at::AccumulatorTuple, ::Val{accname})

Get the accumulator with name accname from at.

getacc(vi::AbstractVarInfo, ::Val{accname})

Return the AbstractAccumulator of vi with name accname.

setacc!!(at::AccumulatorTuple, acc::AbstractAccumulator)

Add acc to at. Returns a new AccumulatorTuple.

If an AbstractAccumulator with the same accumulator_name already exists in at it is replaced. at will never be mutated, but the name has the !! for consistency with the corresponding function for AbstractVarInfo.

setacc!!(vi::AbstractVarInfo, acc::AbstractAccumulator)

Add acc to the AccumulatorTuple of vi, mutating if it makes sense.

If an accumulator with the same accumulator_name already exists, it will be replaced.

setaccs!!(vi::AbstractVarInfo, accs::AccumulatorTuple)
setaccs!!(vi::AbstractVarInfo, accs::NTuple{N,AbstractAccumulator} where {N})

Update the AccumulatorTuple of vi to accs, mutating if it makes sense.

setaccs!!(vi:AbstractVarInfo, accs::AccumulatorTuple) should be implemented by each subtype ofAbstractVarInfo`.

deleteacc!!(at::AccumulatorTuple, ::Val{accname})

Delete the accumulator with name accname from at. Returns a new AccumulatorTuple.

getlogp(vi::AbstractVarInfo)

Return a NamedTuple of the log prior, log Jacobian, and log likelihood probabilities.

The keys are called logprior, logjac, and loglikelihood. If any of them are not present in vi an error will be thrown.

setlogp!!(vi::AbstractVarInfo, logp::NamedTuple)

Set both the log prior and the log likelihood probabilities in vi.

logp should have fields logprior and loglikelihood and no other fields.

See also: setlogprior!!, setloglikelihood!!, getlogp.

acclogp!!(vi::AbstractVarInfo, logp::NamedTuple; ignore_missing_accumulator::Bool=false)

Add to both the log prior and the log likelihood probabilities in vi.

logp should have fields logprior and/or loglikelihood, and no other fields.

By default if the necessary accumulators are not in vi an error is thrown. If ignore_missing_accumulator is set to true then this is silently ignored instead.

getlogjoint(vi::AbstractVarInfo)

Return the log of the joint probability of the observed data and parameters in vi.

See also: getlogprior, getloglikelihood.

getlogjoint_internal(vi::AbstractVarInfo)

Return the log of the joint probability of the observed data and parameters as they are stored internally in vi, including the log-Jacobian for any linked parameters.

In general, we have that:

getlogjoint_internal(vi) == getlogjoint(vi) - getlogjac(vi)

getlogjac(vi::AbstractVarInfo)

Return the accumulated log-Jacobian term for any linked parameters in vi. The Jacobian here is taken with respect to the forward (link) transform.

See also: setlogjac!!.

setlogjac!!(vi::AbstractVarInfo, logjac)

Set the accumulated log-Jacobian term for any linked parameters in vi. The Jacobian here is taken with respect to the forward (link) transform.

See also: getlogjac, acclogjac!!.

acclogjac!!(vi::AbstractVarInfo, logjac; ignore_missing_accumulator::Bool=false)

Add logjac to the value of the log Jacobian in vi.

Errors if vi does not have a LogJacobianAccumulator, unless ignore_missing_accumulator is set to true, in which case this is silently ignored instead.

See also: getlogjac, setlogjac!!.

getlogprior(vi::AbstractVarInfo)

Return the log of the prior probability of the parameters in vi.

See also: getlogjoint, getloglikelihood, setlogprior!!.

getlogprior_internal(vi::AbstractVarInfo)

Return the log of the prior probability of the parameters as stored internally in vi. This includes the log-Jacobian for any linked parameters.

In general, we have that:

getlogprior_internal(vi) == getlogprior(vi) - getlogjac(vi)

setlogprior!!(vi::AbstractVarInfo, logp)

Set the log of the prior probability of the parameters sampled in vi to logp.

See also: setloglikelihood!!, setlogp!!, getlogprior.

acclogprior!!(vi::AbstractVarInfo, logp)

Add logp to the value of the log of the prior probability in vi.

Errors if vi does not have a LogPriorAccumulator, unless ignore_missing_accumulator is set to true, in which case this is silently ignored instead.

See also: accloglikelihood!!, acclogp!!, getlogprior, setlogprior!!.

getloglikelihood(vi::AbstractVarInfo)

Return the log of the likelihood probability of the observed data in vi.

See also: getlogjoint, getlogprior, setloglikelihood!!.

setloglikelihood!!(vi::AbstractVarInfo, logp)

Set the log of the likelihood probability of the observed data sampled in vi to logp.

See also: setlogprior!!, setlogp!!, getloglikelihood.

accloglikelihood!!(vi::AbstractVarInfo, logp)

Add logp to the value of the log of the likelihood in vi.

Errors if vi does not have a LogLikelihoodAccumulator, unless ignore_missing_accumulator is set to true, in which case this is silently ignored instead.

See also: accloglikelihood!!, acclogp!!, getloglikelihood, setloglikelihood!!.

keys(vi::AbstractVarInfo)

Return an iterator over all vns in vi.

getindex(vi::AbstractVarInfo, vn::VarName[, dist::Distribution])
getindex(vi::AbstractVarInfo, vns::Vector{<:VarName}[, dist::Distribution])

Return the current value(s) of vn (vns) in vi in the support of its (their) distribution(s).

If dist is specified, the value(s) will be massaged into the representation expected by dist.

empty!!(vi::AbstractVarInfo)

Empty vi of variables and reset all accumulators.

This is useful when using a sampling algorithm that assumes an empty vi, e.g. SMC.

isempty(vi::AbstractVarInfo)

Return true if vi is empty and false otherwise.

getindex_internal(vi::AbstractVarInfo, vn::VarName)
getindex_internal(vi::AbstractVarInfo, vns::Vector{<:VarName})
getindex_internal(vi::AbstractVarInfo, ::Colon)

Return the internal value of the varname vn, varnames vns, or all varnames in vi respectively. The internal value is the value of the variables that is stored in the varinfo object; this may be the actual realisation of the random variable (i.e. the value sampled from the distribution), or it may have been transformed to Euclidean space, depending on whether the varinfo was linked.

See https://turinglang.org/docs/developers/transforms/dynamicppl/ for more information on how transformed variables are stored in DynamicPPL.

See also: getindex(vi::AbstractVarInfo, vn::VarName, dist::Distribution)

setindex_internal!!(vi::VarInfo, val, vn::VarName)

Set the internal (vectorised) value of variable vn in vi to val.

This does not change the transformation or linked status of the variable.

abstract type AbstractTransformation

Represents a transformation to be used in link!! and invlink!!, amongst others.

A concrete implementation of this should implement the following methods:

• link!!: transforms the AbstractVarInfo to the unconstrained space.
• invlink!!: transforms the AbstractVarInfo to the constrained space.

See also: link!!, invlink!!

struct NoTransformation <: DynamicPPL.AbstractTransformation

Transformation which applies the identity function.

struct DynamicTransformation <: DynamicPPL.AbstractTransformation

Transformation which transforms the variables on a per-need-basis in the execution of a given Model.

This is in constrast to StaticTransformation which transforms all variables before the execution of a given Model.

Different VarInfo types should implement their own methods for link!! and invlink!! for DynamicTransformation.

See also: StaticTransformation.

struct StaticTransformation{F} <: DynamicPPL.AbstractTransformation

Transformation which represents a fixed bijector to be applied to the variables, as opposed to deriving the bijector again at runtime.

See also: DynamicTransformation.

Fields

• bijector::Any: The function, assumed to implement the Bijectors interface, to be applied to the variables

link([t::AbstractTransformation, ]vi::AbstractVarInfo, model::Model)
link([t::AbstractTransformation, ]vi::AbstractVarInfo, vns::NTuple{N,VarName}, model::Model)

Transform variables in vi to their linked space without mutating vi.

Either transform all variables, or only ones specified in vns.

Use the transformation t, or default_transformation(model, vi) if one is not provided.

See also: default_transformation, invlink.

invlink([t::AbstractTransformation, ]vi::AbstractVarInfo, model::Model)
invlink([t::AbstractTransformation, ]vi::AbstractVarInfo, vns::NTuple{N,VarName}, model::Model)

Transform variables in vi to their constrained space without mutating vi.

Either transform all variables, or only ones specified in vns.

Use the (inverse of) transformation t, or default_transformation(model, vi) if one is not provided.

See also: default_transformation, link.

link!!([t::AbstractTransformation, ]vi::AbstractVarInfo, model::Model)
link!!([t::AbstractTransformation, ]vi::AbstractVarInfo, vns::NTuple{N,VarName}, model::Model)

Transform variables in vi to their linked space, mutating vi if possible.

Either transform all variables, or only ones specified in vns.

Use the transformation t, or default_transformation(model, vi) if one is not provided.

See also: default_transformation, invlink!!.

invlink!!([t::AbstractTransformation, ]vi::AbstractVarInfo, model::Model)
invlink!!([t::AbstractTransformation, ]vi::AbstractVarInfo, vns::NTuple{N,VarName}, model::Model)

Transform variables in vi to their constrained space, mutating vi if possible.

Either transform all variables, or only ones specified in vns.

Use the (inverse of) transformation t, or default_transformation(model, vi) if one is not provided.

See also: default_transformation, link!!.

DynamicPPL.update_link_status!!(
    orig_vi::VarInfo, transform_strategy::AbstractTransformStrategy, model::Model
)

Given an original VarInfo orig_vi, update the link status of its variables according to the new transform_strategy.

abstract type AbstractTransformStrategy end

An abstract type for strategies that determine how each variable should be transformed.

For subtypes of AbstractTransformStrategy, the only method that needs to be overloaded is DynamicPPL.target_transform(::AbstractTransformStrategy, vn::VarName), which returns an AbstractTransform that specifies how the variable with name vn should be transformed.

The transform strategy dictates how the log-Jacobian is accumulated during model evaluation. Regardless of what initialisation strategy is used (and what kind of transformed value init() returns, the log-Jacobian that is accumulated is always the log-Jacobian for the forward transform specified by target_transform(strategy, vn).

That is, even if init() returns an UntransformedValue, if the transform strategy is LinkAll() (which returns DynamicLink for all variables), then the log-Jacobian for linking will be accumulated during model evaluation.

Subtypes in DynamicPPL are LinkAll, UnlinkAll, LinkSome, and UnlinkSome.

LinkAll() <: AbstractTransformStrategy

Indicate that all variables should be linked.

UnlinkAll() <: AbstractTransformStrategy

Indicate that all variables should be unlinked.

LinkSome(vns::Set{<:VarName}, fallback) <: AbstractTransformStrategy

Indicate that the variables in vns must be linked. The link statuses of other variables are determined by the fallback strategy.

UnlinkSome(vns::Set{<:VarName}, fallback) <: AbstractTransformStrategy

Indicate that the variables in vns must not be linked. The link statuses of other variables are determined by the fallback strategy.

abstract type AbstractTransform end

An abstract type to represent the intended transformation for a variable.

DynamicLink <: AbstractTransform

A type indicating that a target transformation should be derived by recomputing the invlink transform from the distribution on the right-hand side of the tilde.

Unlink <: AbstractTransform

A type indicating that the target transformation should be nothing.

target_transform(linker::AbstractTransformStrategy, vn::VarName)

Determine whether a variable with name vn should be linked according to the linker strategy. Returns DynamicLink() if the variable should be linked, or Unlink() if it should not.

This function can in the future be extended to support fixed transformations.

DynamicPPL.apply_transform_strategy(
    strategy::AbstractTransformStrategy,
    tv::AbstractTransformedValue,
    vn::VarName,
    dist::Distribution,
)

Apply the given strategy to the transformed value tv for a tilde-statement vn ~ dist.

Specifically, this function does a number of things:

• Calculates the raw value associated with tv.

• Checks whether the strategy expects the VarName vn to be linked or unlinked. If the current link status of tv matches the expected link status, tv is returned unchanged. Otherwise, either linking or unlinking is applied as necessary. Note that this function does not perform vectorisation unless it is needed.

  A table summarising the possible transformations is as follows:

  tv isa ...	target_transform(...) isa DynamicLink	target_transform(...) isa Unlink
  LinkedVectorValue	-> LinkedVectorValue	-> UntransformedValue
  VectorValue	-> LinkedVectorValue	-> VectorValue
  UntransformedValue	-> LinkedVectorValue	-> UntransformedValue

• If vn is supposed to be linked, calculates the associated log-Jacobian adjustment for the forward linking transformation (i.e., from unlinked to linked).

This function returns a tuple of (raw_value, new_tv, logjac).

transformation(vi::AbstractVarInfo)

Return the AbstractTransformation related to vi.

default_transformation(model::Model[, vi::AbstractVarInfo])

Return the AbstractTransformation currently related to model and, potentially, vi.

link_transform(dist)

Return the constrained-to-unconstrained bijector for distribution dist.

By default, this is just Bijectors.bijector(dist).

invlink_transform(dist)

Return the unconstrained-to-constrained bijector for distribution dist.

By default, this is just inverse(link_transform(dist)).

merge(varinfo, other_varinfos...)

Merge varinfos into one, giving precedence to the right-most varinfo when sensible.

This is particularly useful when combined with subset(varinfo, vns).

See docstring of subset(varinfo, vns) for examples.

subset(varinfo::AbstractVarInfo, vns::AbstractVector{<:VarName})

Subset a varinfo to only contain the variables vns.

The ordering of variables in the return value will be the same as in varinfo.

Examples

julia> @model function demo()
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, sqrt(s))
           x = Vector{Float64}(undef, 2)
           x[1] ~ Normal(m, sqrt(s))
           x[2] ~ Normal(m, sqrt(s))
       end
demo (generic function with 2 methods)

julia> model = demo();

julia> params = @vnt begin
           s := 1.0
           m := 2.0
           x := [3.0, 4.0]
       end
VarNamedTuple
├─ s => 1.0
├─ m => 2.0
└─ x => [3.0, 4.0]

julia> vi = last(init!!(model, VarInfo(), InitFromParams(params), UnlinkAll()));

julia> keys(vi)
4-element Vector{VarName}:
 s
 m
 x[1]
 x[2]

julia> # Extract one with only `m`.
       vi_subset1 = subset(vi, [@varname(m),]);

julia> keys(vi_subset1)
1-element Vector{VarName}:
 m

julia> vi_subset1[@varname(m)]
2.0

julia> # Extract one with both `s` and `x[2]`.
       vi_subset2 = subset(vi, [@varname(s), @varname(x[2])]);

julia> keys(vi_subset2)
2-element Vector{VarName}:
 s
 x[2]

julia> vi_subset2[[@varname(s), @varname(x[2])]]
2-element Vector{Float64}:
 1.0
 4.0

subset is particularly useful when combined with merge(vi::AbstractVarInfo)

julia> # Merge the two.
       vi_subset_merged = merge(vi_subset1, vi_subset2);

julia> keys(vi_subset_merged)
3-element Vector{VarName}:
 m
 s
 x[2]

julia> vi_subset_merged[[@varname(s), @varname(m), @varname(x[2])]]
3-element Vector{Float64}:
 1.0
 2.0
 4.0

julia> # Merge the two with the original.
       vi_merged = merge(vi, vi_subset_merged);

julia> keys(vi_merged)
4-element Vector{VarName}:
 s
 m
 x[1]
 x[2]

julia> vi_merged[[@varname(s), @varname(m), @varname(x[1]), @varname(x[2])]]
4-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0

Notes

Type-stability

unflatten!!(vi::AbstractVarInfo, x::AbstractVector)

Return a new instance of vi where the internal values stored in vi.values have been overwritten with the values in x.

This is the inverse operation of internal_values_as_vector.

internal_values_as_vector(vi::AbstractVarInfo)

Return all internal values stored in vi.values as a flattened Vector.

evaluate!!(model::Model, varinfo)

Evaluate the model with the given varinfo, wrapping it in a ThreadSafeVarInfo if the model is marked as needing threadsafe evaluation.

The exact semantics depend on the model's context. Fundamentally, this method executes the model evaluation function (i.e., the function used to define the model) using the given varinfo as an argument. At each tilde-statement, tilde_assume!! or tilde_observe!! is called, whose behaviour depends on the model's context.

Broadly speaking, if the leaf context is an InitContext, then this function:

• uses the initialisation strategy inside the InitContext;
• uses the transform strategy inside the InitContext;
• uses the accumulators inside varinfo (resetting them before evaluation);
• overwrites the values in varinfo with the new values obtained from the initialisation strategy.

If the leaf context is a DefaultContext, then this function:

• uses the values inside the varinfo as the initialisation strategy;
• derives a transform strategy from the varinfo's stored variables (if a linked variable is stored, then the transform strategy will treat that variable as linked; likewise for unlinked)
• uses the accumulators inside varinfo (resetting them before evaluation);
• does not overwrite the values in the varinfo (that is unnecessary since the values used for evaluation are already stored in varinfo).

The long-term plan for this method is to:

• Replace DefaultContext with InitContext by splitting up the functionality of DefaultContext into its constituent components
• Remove the VarInfo argument, and instead use only an AccumulatorTuple
• Separate the initialisation and transform strategies into separate arguments, instead of storing them inside the model's context.

struct DefaultContext <: AbstractContext end

DefaultContext, as the name suggests, is the default context used when instantiating a model.

julia> @model f() = x ~ Normal();

julia> model = f(); model.context
DefaultContext()

As an evaluation context, the behaviour of DefaultContext is to require all variables to be present in the AbstractVarInfo used for evaluation. Thus, semantically, evaluating a model with DefaultContext means 'calculating the log-probability associated with the variables in the AbstractVarInfo'.

InitContext(
    [rng::Random.AbstractRNG=Random.default_rng()],
    strategy::AbstractInitStrategy,
    transform_strategy::AbstractTransformStrategy,
)

A leaf context that indicates that new values for random variables are currently being obtained through sampling. Used e.g. when initialising a fresh VarInfo.

The strategy argument specifies how new values are to be obtained (see AbstractInitStrategy for details), while the transform_strategy argument specifies whether values should be treated as being in linked or unlinked space. That also means that transform_strategy determines whether the log-Jacobian of the link transform is included when evaluating the model.

DynamicPPL.tilde_assume!!(
    context::AbstractContext,
    right::Distribution,
    vn::VarName,
    template::Any,
    vi::AbstractVarInfo
)::Tuple{Any,AbstractVarInfo}

Handle assumed variables, i.e. anything which is not observed (see tilde_observe!!). Accumulate the associated log probability, and return the sampled value and updated vi.

vn is the VarName on the left-hand side of the tilde statement.

template is the value of the top-level symbol in vn.

This function should return a tuple (x, vi), where x is the sampled value (which must be untransformed, i.e., insupport(right, x) must be true!) and vi is the updated VarInfo.

DynamicPPL.tilde_assume!!(
    ::DefaultContext,
    right::Distribution,
    vn::VarName,
    template::Any,
    vi::AbstractVarInfo
)

Handle assumed variables. For DefaultContext, this function extracts the value associated with vn from vi, If vi does not contain an appropriate value then this will error.

DynamicPPL.tilde_assume!!(
    context::AbstractContext,
    right::DynamicPPL.Submodel,
    vn::VarName,
    ::Any,
    vi::AbstractVarInfo
)

Evaluate the submodel with the given context.

DynamicPPL.tilde_observe!!(
    context::AbstractContext,
    right::Distribution,
    left,
    vn::Union{VarName, Nothing},
    vi::AbstractVarInfo
)::Tuple{Any,AbstractVarInfo}

This function handles observed variables, which may be:

• literals on the left-hand side, e.g., 3.0 ~ Normal()
• a model input, e.g. x ~ Normal() in a model @model f(x) ... end
• a conditioned or fixed variable, e.g. x ~ Normal() in a model model | (; x = 3.0).

The relevant log-probability associated with the observation is computed and accumulated in the VarInfo object vi (except for fixed variables, which do not contribute to the log-probability).

left is the actual value that the left-hand side evaluates to. vn is the VarName on the left-hand side, or nothing if the left-hand side is a literal value.

Observations of submodels are not yet supported in DynamicPPL.

This function should return a tuple (left, vi), where left is the same as the input, and vi is the updated VarInfo.

DynamicPPL.tilde_observe!!(
    ::DefaultContext,
    right::Distribution,
    left,
    vn::Union{VarName,Nothing},
    vi::AbstractVarInfo,
)

Handle observed variables. This just accumulates the log-likelihood for left.

AbstractParentContext

An abstract context that has a child context.

Subtypes of AbstractParentContext must implement the following interface:

• DynamicPPL.childcontext(context::AbstractParentContext): Return the child context.
• DynamicPPL.setchildcontext(parent::AbstractParentContext, child::AbstractContext): Reconstruct parent but now using child as its child context.

childcontext(context::AbstractParentContext)

Return the descendant context of context.

setchildcontext(parent::AbstractParentContext, child::AbstractContext)

Reconstruct parent but now using child is its childcontext, effectively updating the child context.

Examples

julia> using DynamicPPL: InitContext, CondFixContext, Condition

julia> ctx = CondFixContext{Condition}(VarNamedTuple(; a = 1));

julia> DynamicPPL.childcontext(ctx)
DefaultContext()

julia> ctx_prior = DynamicPPL.setchildcontext(ctx, InitContext(MersenneTwister(23), InitFromPrior(), UnlinkAll()));

julia> DynamicPPL.childcontext(ctx_prior)
InitContext{MersenneTwister, InitFromPrior, UnlinkAll}(MersenneTwister(23), InitFromPrior(), UnlinkAll())

leafcontext(context::AbstractContext)

Return the leaf of context, i.e. the first descendant context that is not an AbstractParentContext.

setleafcontext(left::AbstractContext, right::AbstractContext)

Return left but now with its leaf context replaced by right.

Note that this also works even if right is not a leaf context, in which case effectively append right to left, dropping the original leaf context of left.

Examples

julia> using DynamicPPL: leafcontext, setleafcontext, childcontext, setchildcontext, AbstractContext, InitContext

julia> struct ParentContext{C} <: AbstractParentContext
           context::C
       end

julia> DynamicPPL.childcontext(context::ParentContext) = context.context

julia> DynamicPPL.setchildcontext(::ParentContext, child) = ParentContext(child)

julia> Base.show(io::IO, c::ParentContext) = print(io, "ParentContext(", childcontext(c), ")")

julia> ctx = ParentContext(ParentContext(DefaultContext()))
ParentContext(ParentContext(DefaultContext()))

julia> # Replace the leaf context with another leaf.
       leafcontext(setleafcontext(ctx, InitContext(MersenneTwister(23), InitFromPrior(), UnlinkAll())))
InitContext{MersenneTwister, InitFromPrior, UnlinkAll}(MersenneTwister(23), InitFromPrior(), UnlinkAll())

julia> # Append another parent context.
       setleafcontext(ctx, ParentContext(DefaultContext()))
ParentContext(ParentContext(ParentContext(DefaultContext())))

setleafcontext(model::Model, context::AbstractContext)

Return a new Model with its leaf context set to context. This is a convenience shortcut for contextualize(model, setleafcontext(model.context, context)).

init!!(
    [rng::Random.AbstractRNG,]
    model::Model,
    varinfo::AbstractVarInfo,
    init_strategy::AbstractInitStrategy,
    [transform_strategy::AbstractTransformStrategy=get_transform_strategy(varinfo),]
)

Evaluate the model and replace the values of the model's random variables in the given varinfo with new values, using a specified initialisation strategy. If the values in varinfo are not set, they will be added using a specified initialisation strategy.

transform_strategy tells the model evaluation whether variables should be interpreted as linked or unlinked. Right now, it is slightly complicated because the default behaviour depends on the varinfo provided. If varinfo isa VarInfo, then the transform strategy is inferred from the VarInfo, i.e., linked variables in the VarInfo are treated as linked during evaluation. Conversely, if varinfo isa OnlyAccsVarInfo, then you must specify the transform strategy explicitly, since an OnlyAccsVarInfo does not contain any information about which variables are transformed.

Returns a tuple of the model's return value, plus the updated varinfo object.

AbstractInitStrategy

Abstract type representing the possible ways of initialising new values for the random variables in a model (e.g., when creating a new VarInfo).

Any subtype of AbstractInitStrategy must implement the DynamicPPL.init method, and in some cases, DynamicPPL.get_param_eltype (see its docstring for details).

init(rng::Random.AbstractRNG, vn::VarName, dist::Distribution, strategy::AbstractInitStrategy)

Generate a new value for a random variable with the given distribution.

This function must return an AbstractTransformedValue.

If strategy provides values that are already untransformed (e.g., a Float64 within (0, 1) for dist::Beta, then you should return an UntransformedValue.

Otherwise, often there are cases where this will return either a VectorValue or a LinkedVectorValue, for example, if the strategy is reading from an existing VarInfo.

DynamicPPL.get_param_eltype(strategy::AbstractInitStrategy)

Return the element type of the parameters generated from the given initialisation strategy.

The default implementation returns Any. However, for InitFromParams which provides known parameters for evaluating the model, methods are sometimes implemented in order to return more specific types.

In general, if you are implementing a custom AbstractInitStrategy, correct behaviour can only be guaranteed if you implement this method as well. However, quite often, the default return value of Any will actually suffice. The cases where this does not suffice, and where you do have to manually implement get_param_eltype, are explained in the extended help (see ??DynamicPPL.get_param_eltype in the REPL).

Extended help

There are a few edge cases in DynamicPPL where the element type is needed. These largely relate to determining the element type of accumulators ahead of time (before evaluation), as well as promoting type parameters in model arguments. The classic case is when evaluating a model with ForwardDiff: the log-probability accumulators must be promoted to contain Duals, and any Vector{Float64} arguments must be promoted to Vector{Dual}. Other tracer types, for example those in SparseConnectivityTracer.jl, also require similar treatment.

If the AbstractInitStrategy is never used in combination with tracer types, then it is perfectly safe to return Any. This does not lead to type instability downstream because the actual accumulators will still be created with concrete Float types (the Any is just used to determine whether the float type needs to be modified).

In case that wasn't enough: in fact, even the above is not always true. Firstly, the accumulator argument is only true when evaluating with ThreadSafeVarInfo. See the comments in DynamicPPL.unflatten!! for more details. For non-threadsafe evaluation, Julia is capable of automatically promoting the types on its own. Secondly, the promotion only matters if you are trying to directly assign into a Vector{Float64} with a ForwardDiff.Dual or similar tracer type, for example using xs[i] = MyDual. This doesn't actually apply to tilde-statements like xs[i] ~ ... because those use Accessors.set under the hood, which also does the promotion for you. For the gory details, see the following issues:

• https://github.com/TuringLang/DynamicPPL.jl/issues/906 for accumulator types
• https://github.com/TuringLang/DynamicPPL.jl/issues/823 for type argument promotion

get_param_eltype(varinfo::AbstractVarInfo, context::AbstractContext)

Get the element type of the parameters being used to evaluate a model, using a varinfo under the given context. For example, when evaluating a model with ForwardDiff AD, this should return ForwardDiff.Dual.

By default, this uses eltype(varinfo) which is slightly cursed. This relies on the fact that typically, before evaluation, the parameters will have been inserted into the VarInfo's metadata field.

For InitContext, it's quite different: because InitContext is responsible for supplying the parameters, we can avoid using eltype(varinfo) and instead query the parameters inside it. See the docstring of get_param_eltype(strategy::AbstractInitStrategy) for more explanation.

AbstractTransformedValue

An abstract type for values that enter the DynamicPPL tilde-pipeline.

These values are generated by an AbstractInitStrategy: the function DynamicPPL.init should return an AbstractTransformedValue.

Each AbstractTransformedValue contains some version of the actual variable's value, together with a transformation that can be used to convert the internal value back to the original space.

Current subtypes are VectorValue, LinkedVectorValue, and UntransformedValue. DynamicPPL's VarInfo type stores either VectorValues or LinkedVectorValues internally, depending on the link status of the VarInfo.

Subtypes of this should implement the following functions:

• DynamicPPL.get_transform(tv::AbstractTransformedValue): Get the transformation that converts the internal value back to the original space.

• DynamicPPL.get_internal_value(tv::AbstractTransformedValue): Get the internal value stored in tv.

• DynamicPPL.set_internal_value(tv::AbstractTransformedValue, new_val): Create a new AbstractTransformedValue with the same transformation as tv, but with internal value new_val.

• DynamicPPL.VarNamedTuples.vnt_size(tv::AbstractTransformedValue): Get the size of the original value before transformation.

VectorValue{V<:AbstractVector,T,S}

A transformed value that stores its internal value as a vectorised form. This is what VarInfo sees as an "unlinked value".

These values can be generated when using InitFromParams with a VarInfo's internal values.

LinkedVectorValue{V<:AbstractVector,T,S}

A transformed value that stores its internal value as a linked andvectorised form. This is what VarInfo sees as a "linked value".

These values can be generated when using InitFromParams with a VarInfo's internal values.

UntransformedValue{V}

A raw, untransformed, value.

These values can be generated from initialisation strategies such as InitFromPrior, InitFromUniform, and InitFromParams on a standard container type.

get_transform(tv::AbstractTransformedValue)

Get the transformation that converts the internal value back to the raw value.

get_internal_value(tv::AbstractTransformedValue)

Get the internal value stored in tv.

set_internal_value(tv::AbstractTransformedValue, new_val)

Create a new AbstractTransformedValue with the same transformation as tv, but with internal value new_val.

ParamsWithStats

A struct which contains parameter values extracted from a VarInfo, along with any statistics associated with the VarInfo. The statistics are provided as a NamedTuple and are optional.

AbstractMCMC.from_samples(
    ::Type{MCMCChains.Chains},
    params_and_stats::AbstractMatrix{<:DynamicPPL.ParamsWithStats}
)

Convert an array of DynamicPPL.ParamsWithStats to an MCMCChains.Chains object.

AbstractMCMC.to_samples(
    ::Type{DynamicPPL.ParamsWithStats},
    chain::MCMCChains.Chains,
    model::DynamicPPL.Model
)

Convert an MCMCChains.Chains object to an array of DynamicPPL.ParamsWithStats.

For this to work, chain must contain the varname_to_symbol mapping in its info field.