1 Context and motivation #

1.1 Goals and target audience #

Statistical modeling in the sciences requires tools that are both mathematically rigorous and practically durable. However, there is still a lack of common standards and infrastructure to express, share, combine and evaluate statistical models across modeling tools and languages. Models should be FAIR (Findable, Accessible, Interoperable, Reusable). We would like to evaluate the same models using different computational engines in multiple end-user languages — in physics, especially C++, Python, and Julia — on a wide variety of modern hardware platforms (including accelerator hardware).

This document proposes FlatPPL, a declarative statistical modeling language that aims to move us closer to these goals, both directly and in connection with existing statistical languages, standards and tools.

FlatPPL is inspired by modeling needs in physics. High Energy Physics (HEP) in particular has a decades-long tradition of rigorous statistical analysis, with code lifetimes measured in decades and a strong culture of reproducibility and model preservation. The particle physics community and related fields like astrophysics and nuclear physics are the primary initial target audience. FlatPPL itself, however, is carefully designed not to be physics-specific, but to be broadly usable for statistical scientific models in general.

New readers may want to read the first four sections (motivation, overview, value types, and language design), then consult the following reference-style chapters (measure algebra, functions, distributions) as needed. Later sections provide worked examples, interoperability guidance and more.

1.2 Probabilistic languages #

A probabilistic language is a formal language for declaring generative models — descriptions of how data could have been produced by a stochastic process (often called forward modeling). The literature partially distinguishes between probabilistic modeling languages and probabilistic programming languages, though the distinction is not always sharp. A probabilistic programming language is often understood to provide both model specification and automatic inference, though not all do. The term probabilistic modeling language is less common, but clearly expresses that inference is not part of the feature set.

FlatPPL is primarily declarative: it describes models, not inference procedures. The scientist writes a model that reads like a simulation recipe: start with a set of parameter values, compute derived quantities, and describe how observations arise from distributions that depend on those parameters. The source model is not an inference procedure or control-flow program. It denotes a static mathematical object that different algorithms can traverse or evaluate in different ways (see below).

FlatPPL does, however, also support likelihood object declarations and density evaluation. Density evaluation defines the semantics of likelihood objects and is also useful for density-based computations within deterministic parts of models. This goes beyond what most probabilistic modeling languages offer, which often have a purely Bayesian focus, but is important for a language that aims to mesh well with formats and frameworks like HS³ and RooFit, and to equally support both frequentist and Bayesian settings.

Algorithms can use a probabilistic model in two fundamental ways, commonly called generative mode and scoring mode:

• Generative mode (simulation): traverses the declared model graph forward and draws random values from probability distributions to produce synthetic data.
• Scoring mode (density evaluation): given parameters and observed values, calculate log-likelihood or log-posterior density values for frequentist and Bayesian inference methods.

In addition to density evaluation, FlatPPL also supports engine-deterministic random-value generation — making both generation and scoring first-class within the language. This enables test-data generation and scoring (likelihood and posterior density values) directly in FlatPPL, though in production these operations are typically driven externally via host-language engine APIs.

Together, generative and scoring mode form the basis for the full range of statistical workflows: maximum likelihood estimation, profile likelihood ratios, Bayesian posterior sampling, hypothesis testing, model comparison, goodness-of-fit checking, and simulation-based inference.

The key design requirements here are:

1. Language-independent. Not tied to a specific programming language. The design must allow for implementation of generative and scoring mode in a wide variety of host languages.
2. Inference-agnostic. Must serve both Bayesian and frequentist use cases.
3. Not tied to a specific engine. No coupling to particular inference algorithms or computational backends.
4. Long-lived. Code lifetimes in HEP have long been measured in decades and data preservation is becoming an increasing concern in many scientific fields. The design must be durable enough to outlast current software and hardware ecosystems.
5. Expressive. Must allow us to express a wide corpus of models across many scientific domains.

Accelerator compatibility. Models that are expressed as a static DAG of bindings — with value shapes that are statically known or resolved at module-load time, no loops, no dynamic control flow, no shapes that change during evaluation, but with explicit support for elementwise operations — map naturally to accelerator-oriented IRs such as MLIR/StableHLO/XLA. Engines targeting high-performance backends (e.g., via JAX in Python or Reactant.jl in Julia) can lower operations on a model, like sampling or density/likelihood evaluation, to these IRs — without fundamental impedance mismatches for the large class of common models with static topology and statically known shapes.

1.3 A starting point: RooFit and HS³ #

The current principal building blocks for statistical modeling in High Energy Physics are RooFit (a C++ modeling toolkit in ROOT) and the HEP Statistics Serialization Standard (HS³), a JSON-based interchange format. pyhf is a JSON specification and Python implementation of the HistFactory template-fitting model class that originated in RooFit. While pyhf comes with a Python engine, both HS³ and pyhf JSON are model descriptions for which engines can be, and currently are being, implemented in multiple languages.

RooFit provides a rich and mature framework for building probability models. Its architecture is based on directed acyclic graphs (DAGs) that express computational dependencies between named objects.

As such, HS³ and pyhf successfully demonstrate how statistical semantics can be disentangled from a DSL that is tied to a specific host language like RooFit/HistFactory. HS³ is young, compared to RooFit, but already in use by the ATLAS collaboration for publishing likelihoods on HEPData.

However, RooFit and HistFactory were built for a specific scientific domain and lack features required in many other domains. HS³ (and pyhf) largely inherited these limitations:

• No distribution/PDF distinction. RooFit conflates distributions with their PDFs, and PDFs do not separate parameters from observables — the distinction arises from usage context (which variables appear in the dataset at fit time). This allows operations such as normalizing a likelihood function over parameter space and treating it as a probability density — an operation that is statistically ill-defined in general, since the likelihood is not a probability measure on parameter space.
• Only scalar variables. All variables are scalar workspace-global objects — there are no vector-valued parameters or variates. Record-like structures (e.g. named components of a multivariate normal) must be flattened into individually named scalars.
• No support for linear algebra.
• No support for complex numbers.
• Stochastic dependencies — one distribution’s variate is another’s parameter — require explicit conditional product construction. Stan-like stochastic graphs are foreign to the RooFit/HS³ approach.

The HS³ standard is written in a more cross-disciplinary fashion, with a clearer separation of some statistical concepts. In particular, it makes the forward-modeling approach that underpins RooFit more explicit. Some of the limitations above may be addressed in future versions of HS³, but others are fundamental, in particular the inability to express stochastic graphs.

The JSON format of HS³ (and pyhf) makes for easy machine-readability, but is not conducive to human readability and human authoring, especially for larger models. Models are often built in RooFit or more specialized tools and then exported to HS³. As RooFit is a fairly narrow C++-embedded DSL that directly expresses computational graphs, exporting to HS³ (built to match RooFit) is straightforward.

This leaves us with an authoring problem: statistical DSLs that are more expressive or allow for broad use of host-language concepts are much harder to map to a language-independent representation, and even harder to map to one not specifically matched to them. On the other hand, models expressed in a standard like HS³ are themselves highly portable, but how do we generate them in scientific domains where models are written in languages or created by tools that may use different semantic models?

1.4 Other probabilistic languages #

There is a rich landscape of stochastic/probabilistic languages; in addition to RooFit, HS³ and pyhf, the following are particularly relevant in our context:

Stan (Carpenter et al., 2017) is a very strong candidate for longevity: it has a large and active user and developer community, bindings for multiple languages (R, Python, Julia and others), and solid funding. However:

• Stan is fundamentally Bayesian, though modern Stan does support some frequentist workflows. The language is designed around the model block, which defines a joint probability distribution over parameters and observations with no syntactic separation between prior and observation model.
• Stan is a full probabilistic programming language with rich syntax, tightly coupled to a specific compiler and runtime (stanc → C++). There is no independent implementation, so it is not suited to be a language-independent interchange format.

Pyro/NumPyro, Turing.jl, PyMC are powerful Bayesian-centric systems. They couple a DSL embedded in their host languages (Python or Julia) with a set of inference algorithms that are specific to each. Independent implementations of the modeling DSLs are not feasible, however, as they leverage very large subsets of the host language and so are tied to that language.

GraphPPL.jl (used by RxInfer) separates model specification from inference backend, which is architecturally what we want. It is Julia-specific and Bayesian-focused and trades generality for high-speed inference.

Hakaru (Narayanan et al., 2016) has elegant semantics built on the Giry monad, expressing programs as measure expressions with support for both frequentist and Bayesian reasoning. Hakaru, however, is based on and tied to Haskell, and does not appear to be actively maintained.

This list should not be seen as exhaustive, nor as an attempt to compare, let alone rate or critique, these languages. Any misrepresentations are unintentional and the authors will be happy to correct them.

All of these languages are powerful and well-designed. FlatPPL is not meant to replace any of them, but to help bridge gaps between expressive power, portability, longevity and DSL-independent model authoring. While FlatPPL is suitable for direct model authoring (see below), different communities use different tools for good reasons. FlatPPL is not an attempt to establish a universal authoring standard for statistical models. Nor is it an attempt to establish a single universal model exchange standard.

1.5 Use cases for FlatPPL #

Model representation and evaluation. Models can be directly authored in or converted to FlatPPL and then evaluated via FlatPPL implementations/engines. FlatPPL is designed for efficient implementation in multiple host languages and use of accelerator hardware. Once stable, FlatPPL will aim for long-term backward-compatibility.

Model conversion. FlatPPL is designed to be a suitable intermediate stage when converting models between different stochastic languages and formats. Tracing compilers and other relevant technologies have become increasingly popular and powerful in recent years, so the technological basis seems largely in place. FlatPPL intentionally supports explicit distribution composition (like RooFit) as well as stochastic graphs (like Stan, Pyro and others) so it can connect across the spectrum of stochastic languages.

Design and reasoning. FlatPPL explores a set of mathematical and statistical semantics that is broad, coherent and rigorous, and has the potential to inform design or extension of other probabilistic languages and standards. FlatPPL also comes with a simple and readable canonical syntax that makes it suitable as a reasoning aid in cases where mathematical notation would still be too informal or too dependent on context. The semantics of FlatPPL are, however, independent of this canonical syntax.

2 Language overview #

2.1 FlatPPL in a nutshell #

The name FlatPPL reflects the language’s most distinctive design choices. Probabilistic models are expressed as static graphs of named mathematical objects — variates, measures, functions, and likelihoods — in a single flat module-level namespace, with no function blocks or block scopes, no loops and no dynamic branching. A FlatPPL document is a sequence of named bindings in static single-assignment (SSA) form. The order of statements is semantically irrelevant; the graph structure is determined by name references, not by textual position. Data is represented by ordinary values (arrays, records, tables).

This simplicity makes FlatPPL amenable to serialization, static analysis, and compilation to accelerator backends, while still being expressive enough to cover a wide range of models across scientific domains. The resulting graph structure is similar to an HS³ JSON document or a RooFit workspace, though FlatPPL concepts like random draws and measure/function reification do not currently exist in HS³ and RooFit. See the profiles and interoperability section on how FlatPPL maps to them.

Creating FlatPPL models. FlatPPL is intended both for direct authoring and as a target representation for models defined in other stochastic languages. Writing FlatPPL documents directly is practical for smaller and medium-sized models and for didactic settings. Users may prefer to use host-language frameworks to author models though, especially larger and complex models, or to use other DSLs to run, exchange and preserve models. FlatPPL is designed to be broad enough to make it a practical target from a wide variety of modeling frontends.

Converting FlatPPL models. FlatPPL defines profiles (see Profiles), specific subsets of the language for easy mapping to target probability languages and formats.

Canonical syntax. FlatPPL comes with a canonical surface form: a small, easy-to-parse syntax. See section Syntax for the grammar and host-language embedding options. This document uses the canonical syntax as a notation to define FlatPPL semantics and to present language examples. Note though that the semantics of FlatPPL are independent from this canonical syntax. The design of FlatPPL allows for alternative surface forms with the same semantics.

Canonical intermediate representation. FlatPPL also comes with a canonical intermediate representation (IR) based on S-expressions, the Flat Probabilistic Intermediate Representation (FlatPIR) (see Intermediate Representation). FlatPIR supports metadata like type and phase annotations and is suitable for automated term-rewriting, enabling code optimization and conversion between FlatPPL profiles. FlatPPL maps directly to and from FlatPIR, enabling round-tripping.

2.2 Implementation targets #

The scientific communities where we expect FlatPPL to see most use primarily work in C++/ROOT, Python, and Julia. Each ecosystem brings its own strengths and infrastructure for statistical modeling:

C++ / RooFit. RooFit is the most mature and widely deployed statistical modeling toolkit in high energy physics, but currently lacks some features required in other fields where it could play a role. FlatPPL is likely to target RooFit via HS³ conversion initially, though direct support is also possible. In either case the implementation strategy is evolution, not replacement: non-breaking additions to widen the semantic scope of RooFit and bring it as close to the scope of FlatPPL as feasible. See the HS³/RooFit profile section for details. Stan, as mentioned before, is very powerful but does not cover all of our requirements. Many strictly Bayesian FlatPPL models could be converted to Stan model blocks though and run on the Stan engine. Accelerator support for RooFit seems less likely for now in general.

Python. pyhf covers the HistFactory subset of HS³, zfit has partial support, and pyhs3 provides a first Python HS³ implementation. In regard to FlatPPL there is more room for direct support in the Python ecosystem than in C++/RooFit. JAX offers a natural path to accelerator-oriented execution via MLIR/StableHLO.

Julia. There is only a prototype HS³ implementation in Julia (HS3.jl). Julia has a rich ecosystem of statistics packages like Distributions.jl and MeasureBase.jl that provide an excellent basis for an inference-agnostic implementation of FlatPPL, orthogonal to inference packages like ProfileLikelihood.jl, BAT.jl and others. FlatPPL and HS³ models could be supported in Julia via the same graph engine. The Julia equivalent to JAX is Reactant.jl, like JAX it targets accelerators via MLIR/StableHLO.

2.3 A first example #

Before delving into the language more formally, here is a small example to convey the flavor of the FlatPPL language. This high energy physics model describes a simple particle mass measurement where the observed spectrum is a superposition of signal and background events, with a systematic uncertainty on the signal resolution:

%%%
# Particle mass measurement

Unbinned Poisson-process model: a Gaussian signal peak at known mass on a
falling exponential background, with a Gaussian-shifted resolution
systematic. Returns a likelihood object `L` parameterized by the expected
signal and background counts.
%%%
flatppl_compat = "0.3"

# Inputs: expected signal and background event counts
% Expected number of signal events.
n_sig = elementof(reals)
% Expected number of background events.
n_bkg = elementof(reals)

% Standard-normal systematic shift applied to the detector resolution.
raw_syst ~ Normal(mu = 0.0, sigma = 1.0)
resolution = 2.5 + 0.3 * raw_syst

# Signal: Gaussian peak at known mass, uncertain resolution
signal_shape = Normal(mu = 125.0, sigma = resolution)

# Background: falling exponential
background_shape = Exponential(rate = 0.05)

# Unbinned data
observed_data = [120.1, 124.8, 125.3, 130.2, 135.7, 142.0]

# Combined intensity: unnormalized superposition (weights = expected event counts)
intensity = superpose(
    weighted(n_sig, signal_shape),
    weighted(n_bkg, background_shape)
)

# Unbinned model: Poisson process over scalar mass values
events ~ PoissonProcess(intensity = intensity)

% Likelihood as a function of the expected event counts.
L = likelihoodof(kernelof(events), observed_data)

The example uses both kinds of textual annotation: # for plain comments (parser-discarded) and % / %%% for doc-comments.

Reading top to bottom, this is a generative recipe: declare inputs, draw a systematic shift, compute the resolution, define signal and background shapes, combine them as an unnormalized superposition (where the weights encode expected event counts), and draw events from the resulting Poisson process. Since the event space is scalar (mass values), the PoissonProcess produces an array variate, and the observed data is a plain array of mass values. (This top-to-bottom reading is for intuition only; semantically the bindings form a dependency graph and may be resolved in any topological order.) The same specification supports generative mode (engines draw synthetic events) and scoring mode (engines compute the log-likelihood at given parameter values).

Note. From a Bayesian perspective, the same model can be read as having a prior Normal(mu = 0.0, sigma = 1.0) on raw_syst, with n_sig and n_bkg as externally supplied hyperparameters or model parameters. This illustrates FlatPPL’s inference-agnostic design.

2.4 Core concepts #

FlatPPL has five kinds of first-class objects.

Abstract values denote real numbers, integers, booleans, complex numbers, fixed-size arrays, and records. They may be deterministic (literal constants, results of ordinary functions, data, external inputs) or stochastic (variates introduced by draw(...)). In generative mode, each abstract value evaluates to a single concrete value; for stochastic abstract values, that concrete value is generated randomly once.

Kernels, measures and distributions. Transition kernels are mappings from an input space to measures. FlatPPL does not distinguish between a kernel with an empty interface and a measure: in FlatPPL, such a kernel is a measure (see variates and measures). Normalized measures (kernels) are probability measures (Markov kernels), also called probability distributions. Variates can only be drawn from probability measures. Otherwise, FlatPPL treats measures and kernels uniformly in measure algebra (see measure algebra). Variates can be reified as probability measures via lawof(...), and as Markov kernels via kernelof(...) (see reification to measures and kernels and kernelof).

Likelihood objects represent the density of a model evaluated at observed data, as a function of the model’s input parameters. The observed data is bound to the likelihood object when it is constructed. To prevent a mix-up of likelihood and log-likelihood values, FlatPPL does not treat a likelihood object as a function that returns the one or the other. Instead, (log-)likelihood values are computed via densityof(L, theta) and logdensityof(L, theta) to make the choice explicit. (See likelihoods and posteriors for the full treatment.)

Functions compute result values from input values in a deterministic fashion. See calling conventions and anonymous functions for details. Values can be reified as functions via functionof(...).

Tuples are ordered bundles of objects. Tuples may contain any objects except for other tuples (see Tuples).

Measures, likelihood objects, functions, and tuples are first-class in the sense that they can be bound to names, passed to operations and referenced by other bindings. However, they may not appear inside arrays, records, or tables.

Modules represent whole FlatPPL documents; each FlatPPL source file is a module. A FlatPPL module can load other modules (via load_module(module_filename)) and access objects in loaded modules via dot-syntax scoping (loaded_module.some_object). Module objects give access to another namespace, but are not themselves first-class objects in the computational graph: they may not be passed to functions or appear inside data structures. See multi-file models for details.

2.5 Language map #

The table below provides a compact overview of the language. Each family name links to the section where the constructs are documented. The lists are selected highlights, not exhaustive — see the linked sections for complete listings.

2.6 A tour of FlatPPL #

The following code blocks illustrate the main language features. Each construct is explained in detail in later sections. These snippets are independent and do not form a single model.

2.6.1 Values and collections #

Scalars, arrays, nested arrays, matrices, records, and basic operations:

# Scalars
x = 3.14
n = 42
b = true

# Collections
v = [1.0, 2.0, 3.0]
nested = [[1, 2], [3, 4]]
M = rowstack([[1, 2, 3], [4, 5, 6]])
r = record(mu=3.0, sigma=1.0)

# Indexing, field access, slicing
y = A[i]
z = A[i, j]
w = r.mu
col_j = M[:, j]

# Decomposition into named scalars
a, b, c ~ MvNormal(mu = mean, cov = cov_matrix)
p, q = some_record
l, m, n = some_tuple

# Arithmetic, comparisons, function calls
rate = efficiency * mu_sig + background
is_positive = x > 0
y = exp(x)
z = ifelse(is_positive, a, b)
combined = cat(record1, record2)
joined = cat(array1, array2)

2.6.2 Function calls #

Calling conventions:

Normal(mu = 0, sigma = 1)           # keyword
Normal(record(mu = 0, sigma = 1))   # record auto-splatting
exp(x)                              # positional

2.6.3 Complex arithmetic #

Constructing and calculating with complex values:

# Complex construction
z1 = complex(3.0, 2.0)
z2 = 3.0 + 2.0 * im               # equivalent
phase = cis(3 * pi / 4)           # unit-modulus complex from angle

# Complex arithmetic
A_total = A_sig * coupling + A_bkg

# Squared modulus (real-valued result)
intensity = abs2(A_total)

# Decomposition and conjugation
x = real(z1)
y = imag(z1)
z_bar = conj(z1)

2.6.4 Draws, measures, and the stochastic core #

draw, lawof, functionof, and kernelof bridge between values, measures, functions, and kernels:

# Inputs
mu = elementof(reals)
sigma = elementof(interval(0.0, inf))

# Random draw from a distribution
a ~ Normal(mu = mu, sigma = sigma)

# Extract the distribution governing a value
M = lawof(a)

b = 2 * a + 1

# Reify a deterministic sub-DAG as a function (boundary stops trace at `a`)
f = functionof(b, x = a)

# Stochastic sub-DAG reified as a kernel
K = kernelof(b, x = a)

2.6.5 Broadcasts #

broadcast applies functions or kernels elementwise over arrays and tables:

# Function over array (keyword binding)
C = broadcast(f, x = A)

# Same, positional
C = broadcast(f, A)

# Kernel over array
D ~ broadcast(K, x = A)

2.6.6 Value-level operations #

Accessing and renaming values, and transforming measures based on that:

# Element and subset access
field_a = get(some_record, "a")
sub = get(some_record, ["a", "c"])

# Array to record conversion
named = relabel(some_array, ["a", "b", "c"])

# Structural relabeling of a measure
mvmodel = relabel(MvNormal(mu = some_mean, cov = some_cov), ["a", "b", "c"])

# Variable transformation
log_normal = pushfwd(x -> exp(x), Normal(mu = 0, sigma = 1))

# Projection (marginalizes out b)
marginal_ac = pushfwd(fn(get(_, ["a", "c"])), mvmodel)

2.6.7 Measure algebra and composition #

Combining, reweighting, and transforming measures some more:

# IID draws
xs ~ iid(Normal(mu = 0, sigma = 1), 100)

# Additive rate superposition
sp = superpose(weighted(n_sig, sig), bkg)

# Normalized mixture
mix = normalize(superpose(
    weighted(0.7, M1), weighted(0.3, M2)))

# Joint of independent components
j = joint(M1, M2)

# Marginalizing composition
pp = kchain(prior, forward_kernel)

# Hierarchical joint (retains both variates)
hj = jointchain(
    pushfwd(fn(relabel(_, ["a"])), M1),
    pushfwd(fn(relabel(_, ["b"])), K_b))

# Truncated (unnormalized) measure
positive_normal = truncate(Normal(mu = 0, sigma = 1),
    interval(0, inf))

# Fundamental measures and density-defined distributions
leb = Lebesgue(support = reals)
bern = fn(bernstein(coefficients = [c0, c1, c2], x = _))
smooth_shape = normalize(weighted(bern, Lebesgue(support = interval(lo, hi))))

2.6.8 Anonymous functions #

The fn(...) form wraps a hole expression — an expression containing _ — to create an anonymous function with positional parameters:

# Single hole — one-argument function
poly = fn(polynomial(coefficients = [a0, a1, a2], x = _))
squared = fn(_ ^ 2)

# Multi-hole: two-argument anonymous function
ratio_sq = fn((_ / _) ^ 2)

2.6.9 Interpolation, binning, and systematic variations #

Constructors for binned models and HistFactory-style yield arithmetic:

edges = linspace(0.0, 10.0, 5)
counts = bincounts(edges, event_data)

# Binned observation model via pushforward
binned_model = pushfwd(fn(bincounts(edges, _)),
    PoissonProcess(intensity = M_intensity))

# Interpolation for systematic variations (particle-physics standard module)
hepphys = standard_module("particle-physics", "0.1")
kappa = hepphys.interp_poly6_exp(0.95, 1.0, 1.05, alpha)
morphed = hepphys.interp_poly6_lin(tmpl_dn, nominal, tmpl_up, alpha)

2.6.10 Data #

Data is represented by ordinary values — no special data type:

observed_counts = [5, 12, 8, 3]
data_table = table(a = [1.1, 1.2], b = [2.1, 2.2])

2.6.11 Presets #

Advisory parameter/input values for use with a compatible function, kernel, or likelihood:

# Parameter starting values (advisory, not part of model semantics)
starting_values = record(mu_sig = 1.0, raw_syst = fixed(0.0), n_bkg = 50.0)

2.6.12 Analysis: likelihoods and posteriors #

Likelihood construction, combination, and posterior construction:

L = likelihoodof(kernelof(obs), data)
R = interval(2.0, 8.0)
L_sub = likelihoodof(normalize(truncate(kernelof(obs), R)), filter(fn(_ in R), data))
L_total = joint_likelihood(L1, L2)

# Unnormalized posterior
posterior = bayesupdate(L, prior)

# Deterministic function composition
pipeline = fchain(calc_kinematics, apply_cuts)

2.6.13 Modules and interface adaptation #

Module loading and parameter renaming:

# Load a module and optionally bind some of its inputs
sig = load_module("signal_channel.flatppl", mu = signal_strength, theta = nuisance)

sig_model = sig.model
L_sig = likelihoodof(sig.model, sig.data)

3 Value types and data model #

FlatPPL has a small, fixed set of value types. This section defines what kinds of values exist in the language, their invariants, and how they interact. Operations on values are documented in built-in functions.

3.1 Scalar types #

Real. Floating-point numbers like 3.14, -0.5, 1e-3.

Integer. Integer numbers like 42, 0, -7.

Bool. true or false (lowercase). In arithmetic contexts, false is promoted to zero and true to one, permitting expressions such as true + true, 3 * false, and sum(mask) to count true entries. Conditional and logical constructs (ifelse, land, lor, lnot, lxor) strictly require boolean arguments; zero and one are not implicitly converted to booleans.

Complex. A complex number. Constructed via complex(re, im) or via arithmetic with the imaginary unit im:

z1 = complex(3.0, 2.0)
z2 = 3.0 + 2.0 * im           # equivalent
phase = cis(3 * pi / 4)       # unit-modulus complex from angle

When a real and a complex value meet in arithmetic, the real is promoted to complex with zero imaginary part.

Scalar value categories and sets. FlatPPL distinguishes boolean, integer, real, and complex scalar values operationally. In particular, conditionals and logical operators require boolean values. However, the predefined value sets satisfy the canonical inclusions booleans $\subset$ integers $\subset$ reals, and there is a canonical embedding of reals into complexes. Arithmetic may use these canonical embeddings implicitly where specified by the language.

3.2 Predefined constants #

Note on infinities. posreals, nonnegreals, and reals admit inf (and, for reals, -inf) as legal values. Strictly speaking, these are subsets of the extended reals $\overline{\mathbb{R}} = \mathbb{R} \cup \{-\infty, +\infty\}$, not of $\mathbb{R}$. This is a deliberate choice for compatibility with common numerical and statistical libraries. When FlatPPL refers to Lebesgue(support = reals), the reference measure is the ordinary Lebesgue measure on the finite-real part; the points $\pm\infty$ carry zero Lebesgue mass. Arithmetic on infinities follows IEEE 754 conventions.

The selector all and the hole token _ are syntactic elements, not value constants; they are documented in calling conventions and anonymous functions.

3.3 Arrays #

Arrays are fixed-size, ordered, n-dimensional collections of scalar values (real, integer, boolean and complex values) or arrays.

Literal one-dimensional arrays are denoted as [1.0, 2.0, 3.0] and may contain arbitrary valid FlatPPL expressions that evaluate to allowed element types (e.g. [a, b, 2 * c]).

One-dimensional arrays of scalars act as vectors for linear algebra (see built-in functions). In addition, transposed vectors are a distinct type in FlatPPL (see linear algebra). The term vector will represent both non-transposed vectors (one-dimensional arrays) and transposed vectors in the following, unless noted otherwise.

Vectors of vectors are not interpreted as matrices implicitly, but can be turned into matrices explicitly using rowstack or colstack (see array operations).

FlatPPL supports standard linear algebra operations (addition, multiplication) on scalars, vectors, and matrices.

3.4 Records #

Records comprise ordered named fields, written as record(name1=val1, name2=val2, ...). Field values may be scalars, arrays, or records. Field access uses dot syntax: r.name1 (lowers to get(r, "name1")); nested fields chain, e.g. r.a.b. Field order is part of the record’s identity: record(a=1, b=2) and record(b=2, a=1) are distinct values. This is significant for alignment with parameter spaces and for deterministic serialization. Fields are accessed by name, not by position — get(r, i) is not supported to avoid ambiguity with row indexing on tables.

3.5 Presets #

Presets are global fixed records and sets that hint suitable parameter/input values and value domains to FlatPPL tooling. Presets are advisory and not tied to a particular function, kernel, or likelihood. It is up to users and tooling to pair presets with compatible interfaces and decide how to use them, for example as reference points, starting values for optimizers, value ranges for optimization or plotting, etc.

Preset points. Any literal (or fixed, in general) global binding some_name = record(name1=val1, name2=val2, ...) can be interpreted as a possibly suitable input for functions, kernels and likelihoods that have inputs/parameters with these names and shapes or that take a record of this shape as an input. Tooling may offer such preset points to users to choose from.

Values in a preset record that are wrapped in fixed(...) indicate that these values should be held constant while others are varied, e.g. during optimization. fixed(x) is semantically identical to identity(x) during FlatPPL code evaluation, it is merely a hint to tooling. The hint is only meaningful when fixed(...) appears directly as a preset-record field value.

For example:

L_init = record(a = 2.0, b = [4, 5, 6], c = fixed(8.0))

Preset domains. Any literal/fixed global binding like some_name = cartprod(name1=some_set, name2=some_other_set, ...) can be interpreted as a possibly suitable domain for input/parameter-compatible functions, kernels and likelihoods. Like with preset points, tooling may offer such preset domains to users to choose from.

For example:

L_domain = cartprod(a = interval(0, 5), b = cartpow(interval(-10, 10), 3), c = interval(0, 20))

3.6 Tables #

Tables are datasets that consist of named columns. All columns must have the same length (row count). Each column is a vector or a table; a vector column’s elements may themselves be arrays (e.g. a 3-vector per entry). A column may not itself be a higher-dimensional array, as this would require a leading-axis convention for row iteration and broadcasting, which FlatPPL intentionally avoids. Each row of a table is a record; if some columns of the table are tables themselves, the corresponding entries of the row records are records themselves.

Tables are constructed from columns via table(col1 = [...], col2 = [...]):

events = table(mass = [1.1, 1.2, 1.3], pt = [45.2, 32.1, 67.8])

Implementations may choose whichever table realization and memory layout they prefer, also on a case-by-case basis.

Tables can also be constructed from records of equal-length vectors via table(r) and converted to such records via record(t).

Indexing. Tables support both column and row access:

• Column access by field name: t.colname, equivalent to get(t, "colname"), returns the column with that name as a vector.
• Row access by integer index: t[i], equivalent to get(t, i), returns the i-th row as a record (i starts at 1).

lengthof(t) returns the number of table rows.

Broadcasting. When a table is passed to broadcast, it is traversed row-wise and each row treated as a record passed to the function used in the broadcast.

Data carriers by model shape. FlatPPL uses ordinary values as data carriers:

• Single scalar datum → scalar value
• Single structured datum → record or array
• Unbinned scalar event sample → plain array
• Unbinned multivariate event sample → table
• Binned count data → plain count array

3.7 Sets #

FlatPPL has a limited notion of sets, used to specify input domains, supports, truncation regions, and analysis regions. The predefined sets are:

• reals — $\mathbb{R}$, the real numbers ($\pm\infty$ admitted).
• posreals — $(0, +\infty]$, the positive reals including $+\infty$.
• nonnegreals — $[0, +\infty]$, the non-negative reals including $+\infty$.
• unitinterval — $[0, 1]$, the unit interval.
• posintegers — $\{1, 2, 3, \ldots\}$, the positive integers.
• nonnegintegers — $\{0, 1, 2, \ldots\}$, the non-negative integers.
• integers — $\mathbb{Z}$, the set of all integers.
• booleans — $\{\mathrm{false}, \mathrm{true}\}$.
• complexes — $\mathbb{C}$, the set of all complex numbers.
• anything — a broad placeholder set for generic interfaces (e.g., anonymous functions via holes). Not formally the union of all other sets; it signals that no specific type constraint is imposed.
• rngstates — the set of RNG state values; members are algorithm-dependent opaque values (see random number generation).

Additional sets may be constructed using the following language constructs:

Interval. interval(lo, hi) denotes the closed interval $[lo, hi]$.

Cartesian product. cartprod(S1, S2, ...) produces a Cartesian product of sets S1, S2, etc., mirroring joint(M1, M2, ...) for measures. Each member is the cat of one element per component set (so vector components concatenate). The resulting set is a set of arrays, not a set of tuples. The element type is the common type of the component element types. For example, cartprod(reals, posreals) is the set of 2-element arrays [a, b] with a in $\mathbb{R}$ and b in $(0, +\infty]$, and cartprod(reals, integers) is the set of real 2-vectors [a, b] with a in $\mathbb{R}$ and b in $\mathbb{Z}$ — real-valued since integers $\subset$ reals, with the second element restricted to integers.

The keyword form cartprod(a = S1, b = S2, ...) produces a set of records with field a in S1, field b in S2, etc., mirroring joint(a = M1, b = M2, ...).

Cartesian power. cartpow(S, size) produces the Cartesian power of S with shape size, mirroring iid(M, size) for measures. size is a positive integer (1-D) or a vector of positive integers (multi-axis). For example, cartpow(reals, 3) represents $\mathbb{R}^3$ and cartpow(reals, [3, 3]) the set of $3 \times 3$ real matrices. When S is a record set, the power is the set of tables with those columns: cartpow(cartprod(a = reals, b = posreals), n) is the set of n-row tables.

Standard simplex. stdsimplex(n) denotes the standard $(n{-}1)$-dimensional probability simplex $\Delta_{n-1} = \{x \in \mathbb{R}^n : x_i \geq 0,\; \sum_i x_i = 1\}$. Lebesgue(support = stdsimplex(n)) is the $(n{-}1)$-dimensional coordinate Lebesgue measure on the simplex: the image of $dx_1 \cdots dx_{n-1}$ under the chart that appends $x_n = 1 - \sum_{i<n} x_i$ (dropping any other coordinate gives the same measure). It assigns zero mass to sets that do not intersect the simplex. It is not the surface (Hausdorff) measure of the embedded simplex, which is larger by the factor $\sqrt{n}$.

relabel applies to set products in the same way as to measures (see interface adaptation).

Sets that govern values. valueset(x) returns the canonical value set associated with node x:

• For x = elementof(S), valueset(x) is S.
• For x = draw(M), valueset(x) is the measurable set of M.
• For deterministically computed nodes, valueset(x) returns a conservative superset of the values that x can take (since the exact set is often not tractable).

Note: valueset is a low-level language construct used when lowering functionof or kernelof with boundary inputs. User-level code should typically use elementof(...) and specify sets explicitly.

3.8 Beyond values #

Measures, likelihood objects, functions, and tuples are also first-class objects in FlatPPL — they can be bound to names, passed to combinators, and referenced by other bindings. However, they are not value types: they may not appear inside arrays, records, or tables. Modules are namespace references rather than computational objects; they may be bound to names but cannot be passed as arguments or stored in containers. See core concepts for details.

4 Language design #

This section details the semantics of FlatPPL’s core constructs: modules, objects, and callables, and built upon them variates, measures, deterministic and stochastic graphs, and more.

4.1 Objects, expressions, names and modules #

The FlatPPL language consists of objects (measures, likelihood objects, tuples, or values like numbers, arrays, records and tables). Objects include callables (value functions, constructors, transition kernels and special operations) that operate on objects.

Ordinary callables have named inputs and a single output; that output may be a tuple bundling multiple components (see Tuples). Special operations are callables that handle inputs in a different way and typically provide higher-level semantics. The output of any callable depends deterministically on its inputs and calls may not have any side effects. No callables may have nullary inputs, as this would make them equivalent to known values.

Numerical precision (e.g., 32-bit vs. 64-bit floating point) is not specified by FlatPPL; the choice is left to implementations and their users.

A FlatPPL module is an unordered set of bindings of names to expressions. Expressions are single or nested calls that bind expressions (literal or by name reference) to inputs of callables.

FlatPPL is loop-free and has no block structure, so a module is implicitly a directed acyclic graph (DAG), which may not be fully connected. The nodes in that graph are the named and unnamed (results of) calls, the edges are the connections between outputs and inputs of calls.

The graphs of several modules can be combined, see Module composition below for details.

Note: Record field names and table column names are local to their object and not part of the global module namespace, nor are the argument names of functions and kernels.

4.2 Binding names #

Public bindings. Names that do not begin with an underscore are public: they form the interface of a FlatPPL module. They must match the regular expression ^[A-Za-z][A-Za-z0-9_]*$.

Private bindings. Binding names that begin with a single underscore and do not end with an underscore (regular expression ^_[A-Za-z]([A-Za-z0-9_]*[A-Za-z0-9])?$), e.g. _tmp, are private to a module. They are not part of the module’s public interface and may be eliminated, inlined, renamed, or otherwise not preserved by tooling such as term-rewriting or dead-code elimination.

The bare underscore _ is itself a valid binding name. Each occurrence of _ on the left-hand side (standalone or inside a decomposition, e.g. value, _ = rand(rstate, m)) lowers to a distinct auto-generated private name. So _ may be used to discard values.

Auto-generated names. Names starting with a double underscore (regular expression ^__[A-Za-z0-9][A-Za-z0-9_]*$) are reserved for automatically generated module-level binding names. Like other underscore-prefixed names they are private and elidable. Auto-generated binding names with a purely numerical ID are denoted in hexadecimal form with a __0x prefix in the canonical FlatPPL and FlatPIR syntax (regular expression ^__0x[0-9a-f]+$), but may be encoded differently in other representations of the language.

Placeholder names. Names starting and ending with a single underscore (regular expression ^_[A-Za-z]([A-Za-z0-9_]*[A-Za-z0-9])?_$) are reserved for placeholder variables inside functionof and kernelof (see placeholders and holes). They must not be used for module-level bindings.

Name resolution. FlatPPL has two predefined modules: self refers to the current module, and base is a predefined module containing FlatPPL’s built-ins. self.foo always refers to a current-module binding; base.foo always refers to a built-in. self and base themselves are reserved names and cannot be bound to something else.

An unqualified name (no self. or base. prefix) resolves as follows:

1. If the name is bound in the current module, it resolves to that binding.
2. Otherwise, it resolves to the FlatPPL built-in of that name.

Unresolvable names are static errors.

This makes built-in names shadowable: a module may bind any name except for self and base. Adding new built-ins to FlatPPL is therefore a non-breaking change.

4.3 Calling conventions #

Nullary calls (f()) are not allowed.

All ordinary callables — built-in or user defined value functions, constructors or transition kernels — accept arguments in two or three forms (denoted here in the canonical syntax):

• Keyword arguments (named arguments): f(a = x, b = y, ...). Arguments are bound to inputs by name, the order of the arguments is not relevant.

• Auto-splatting (of records and table columns): f(record(a = x, b = y, ...)) and f(table(a = x, b = y, ...)) are equivalent to f(a = x, b = y, ...). The order of fields or columns is not relevant. A call with field or column names that do not match the callable’s argument names is a static error. Auto-splatting is shallow and occurs only when a record or table is the call’s sole argument (whatever its field count, a single field included); a record given alongside other arguments, or bound to a parameter by keyword, is an ordinary value and is not splatted. A sole positional record or table therefore always splats: whether its field or column names match the callable’s argument names decides only whether the call is valid, never whether the splat occurs. A callable with exactly one input whose documented domain admits records or tables is exempt and receives a sole positional record or table whole, so sum(t) reduces over the table rather than splatting. User-defined callables are never exempt. Passing a record or table as one ordinary argument requires the keyword spelling, as in f(pars = record(...)). Auto-splatting is a rule of the ordinary calling convention; no special operation splats a sole record or table argument. table(r) and record(t) perform the record–table conversions of tables directly, as dedicated conversions rather than as an instance of auto-splatting.

• Positional arguments: f(x, y, ...). Positional arguments are accepted only if the callable has ordered inputs, so that the arguments can be mapped to the inputs in order.

All built-in ordinary callables have a defined input order and accept both positional and keyword arguments.

Special operations have zero to three distinguished inputs: unnamed, ordered inputs of fixed arity. They may have additional variadic named or unnamed inputs, the order of which may or may not be significant. The total number of inputs is never zero:

• elementof, external, draw: One distinguished input.
• vector: Unnamed variadic inputs with significant order.
• tuple: Unnamed variadic inputs with significant order (minimum two).
• record and table: Named variadic inputs with significant order.
• functionof and kernelof: One distinguished input, plus optional variadic named inputs with significant order.
• lawof, fixed: One distinguished input.
• broadcast: One distinguished input for the function to be broadcast, plus named or unnamed inputs that match the inputs of that function.
• broadcasted: One distinguished input.
• cat, fchain, kchain: Variadic unnamed inputs with significant order.
• cartprod, joint, jointchain: Variadic unnamed or named inputs with significant order.
• get: One distinguished input plus variadic unnamed input with significant order.
• superpose: Variadic unnamed inputs with no significant order.
• load_module: One distinguished input plus optional variadic named inputs with no significant order.
• standard_module: Two distinguished inputs.
• aggregate, metricsum, markovchain, kscan: Three distinguished inputs.
• ksuperpose: Two distinguished inputs (the kernel and the weight vector); the resulting kernel is applied separately to the parameter family.
• load_data: One distinguished input plus optional variadic named inputs with significant order.
• checked: Named parameters value and condition, per §07; the canonical calling form is keyword-based.

A distinguished input has no name and so cannot be passed by keyword. The measure combinators likewise take their inputs positionally: a keyword spelling such as normalize(M = mu) is a static error. Where this specification refers to a distinguished input by a name, as in aggregate(f_reduction, output_axes, expr), the name identifies the input in prose only. A call binds the input by position, never by keyword argument.

4.4 Tuples #

Some operations produce a single output that naturally groups several distinct components — e.g. a kernel and its base measure together, or an updated RNG state alongside a generated value. Tuples package such outputs as an ordered, fixed-length bundle of FlatPPL objects.

The surface form (a, b, c) lowers to tuple(a, b, c). Tuples must contain at least two components, so (), (x,) are not allowed. Tuple elements are accessed via t[i], lowering to get(t, i), with a positive integer literal index (starting at 1). Decomposition as in a, b, c = (...) is positional.

Tuples are objects, not values. They have no valueset, are not drawn from measures, and are not part of the measure algebra. Specifically:

• A tuple may not appear inside an array, record, or table, but may appear inside another tuple (tuples nest).
• elementof(...) and external(...) may not produce tuples.
• Measures, kernels, and likelihoods never use tuples as their domain.
• ==/equal does not compare tuples.
• Tuples do not auto-splat like records and tables do.

Tuples otherwise flow like other objects. They may be bound to names, passed to callables that accept them, returned from user-defined functions, and decomposed or projected.

4.5 Variates and measures #

FlatPPL distinguishes variates from measures and kernels. A variate represents a specific value — one realization in any given evaluation of the model. A measure or kernel, by contrast, represents the entire distribution over possible values. More formally, measures are monadic while variates are not.

Keeping variates and measures distinct matters because arithmetic means different things for each: In mathematics, $2 \cdot x$ transforms a variate (producing a new variate with twice the value), while $2 \cdot \mu$ rescales the measure (multiplying its mass on every set by 2). FlatPPL supports both via different syntax — arithmetic on variates, weighted(...) on measures.

A binding of the form c = f(a, b) introduces a deterministic node in the computational DAG. A binding of the form x ~ Normal(mu = c, sigma = s), equivalent to x = draw(Normal(mu = c, sigma = s)), introduces a stochastic node. In generative mode, a stochastic node yields a sampled value; in scoring mode, it contributes a density term that is either evaluated (if observed) or marginalized out (if latent).

FlatPPL intentionally supports two equivalent mechanisms to express stochastic computations:

1. Stochastic-node notation expresses models as a mix of deterministic computations and draw statements, reading like a generative recipe.
2. Measure-composition notation writes models as a mix of deterministic computations and measure algebra, using weighted, joint, jointchain, kchain, pushfwd, and related operations to combine and transform measures.

Both can be used together in a FlatPPL module, but they map to different types of probabilistic coding systems. Stochastic-node notation mirrors probabilistic programming languages like Stan and Pyro, while measure composition mirrors the HS³ and RooFit approach. By supporting both approaches, FlatPPL can be emitted from both types of systems. Term-rewriting via FlatPIR can raise and lower code to match either of them. This enables FlatPPL and FlatPIR to act as an interoperability platform.

4.6 Internal parameters and external inputs #

FlatPPL declares unresolved values via two special operations:

• elementof(S) — declares a module-internal parameter that takes a specific value during a single evaluation of a subgraph that contains it; the value may vary between evaluations (e.g., during parameter inference).
• external(S) — declares a module-external input. The value must be supplied by applications that use the module or by modules that load the containing module and bind this input to a fixed value in their own namespace. Module inputs can be thought of as hyperparameters and their values don’t change between subgraph evaluations.

n_dims = external(posintegers)
mu = elementof(reals)
sigma = elementof(interval(0.0, inf))
dist = iid(Normal(mu = mu, sigma = sigma), n_dims)
x ~ dist
y = 2 * x

The distinction between external and elementof determines phase classification (see below), the ancestors functionof resolves to values (see application and reification), and cross-module binding rules for load_module (see Multi-file models).

4.7 Phases #

FlatPPL classifies every binding into one of three phases by ancestor analysis:

• fixed — no elementof(...) ancestor and no draw(...) ancestor, but may have load_data(...) ancestors.
• parameterized — at least one elementof(...) ancestor, no draw(...) ancestor.
• stochastic — at least one draw(...) ancestor.

External inputs and loaded data (length and content) are fixed; elementof inputs are parameterized; draw nodes are stochastic. Phase propagates through the DAG: a binding’s phase is the dominant of its ancestors’ phases (stochastic > parameterized > fixed).

Records, arrays, tables, and tuples may bundle components of differing phases; under the ancestor rule the container carries the joined phase, and projections inherit it. This is conservative — a projection such as r.a technically depends on the entire record’s ancestors, even though its value is just one field. Engines may sharpen this by flattening projections with statically known selectors (r.field, t[i] with integer literal index, or the corresponding get(...) and decomposition forms) before phase analysis or the ancestor trace, which recovers the selected component’s phase directly.

Both fixed and parameterized bindings are deterministic, but their values have different life cycles: A FlatPPL module can be thought of as having an initialized state, where external inputs have been set and data and referenced modules have been loaded. Fixed values are the values that are given or deterministically computable at this point and so do not change after module initialization. Parameterized values differ between evaluations of the same subgraph (e.g. of a likelihood) of the initialized module, given different inputs. Note that this is a mental model, applications are not required to use an explicit initialization state to implement these semantics.

Phase governs which ancestors are resolved to values at reification (see application and reification) and load-time binding rules (see Multi-file models).

4.8 Application and reification #

FlatPPL provides operations that turn subgraphs into first-class objects and vice versa.

A function represents a reified deterministic DAG, either implicit (built-in) or explicitly constructed. Ordinary function application y = f(a, b, ...) introduces a deterministic node y into the graph. functionof(y) goes in the opposite direction: it reifies the ancestor subgraph of y as a first-class function — the backward program slice of y (Weiser, 1981).

Conversely, a probability measure represents a reified stochastic DAG, either implicit (built-in) or explicit. x ~ m (equivalent to x = draw(m)) introduces a stochastic node x by drawing a variate from a normalized measure (i.e. a probability measure) m. In the other direction, m = lawof(x) reifies the ancestor subgraph of x as a probability measure — the law of x as a random variable. The identity law relating the two directions is stated under reification to measures below.

draw differs fundamentally from rand: rand produces a concrete random value, while draw introduces a stochastic node that represents the existence of a random value. Applying rand to a measure reified from a DAG that contains draw nodes resolves those draws to concrete random values; evaluating a density on the same measure instead takes those stochastic nodes as inputs to the density (unless they are marginalized out in the stochastic graph).

Reifying measure-valued expressions to kernels. If functionof is applied to a measure node, it generates a transition kernel — a measure-generating callable in FlatPPL — instead of a function. If the measure is normalized, the resulting kernel is a Markov kernel. For stochastic-phase m, functionof(m, ...) reifies the conditional kernel and functionof(lawof(m), ...) the marginal (see reification to measures); for fixed or parameterized phase the two coincide.

So functionof reifies a value or measure node in the computational graph to a value function or transition kernel. While lawof reifies a value node in the computational graph to a probability measure, and a measure node to the law of a draw from it (see below).

kernelof (see below) combines lawof and functionof.

4.8.1 Reification to measures #

lawof(x) reifies the ancestor sub-DAG of x as the probability measure that is the total law of x — the probability measure that x, considered as a random variable, is distributed according to.

Identity law. lawof(draw(m)) is equivalent to m for m of fixed or parameterized phase; for stochastic-phase m it is the marginal law of a draw from m. Equal laws do not make values interchangeable as joint components: a joint of two reified laws of the same draw is the singular diagonal joint. Otherwise joint(m, m) contributes a fresh coordinate per occurrence, so the two draws are independent given m’s stochastic ancestors — which remain shared — and independent outright when m has none.

lawof also accepts a measure argument: lawof(m) is lawof(draw(m)), the law of a draw from m. Each draw from m is a fresh coordinate, while the draw nodes among m’s ancestors remain the same nodes of the trace. A probability measure of fixed or parameterized phase is its own law, so lawof(m) is equivalent to m and lawof is idempotent. For stochastic-phase m — a random measure — lawof(m) is the marginal law of a draw from it: the mixture $\nu(B) = \int \kappa(z, B)\, dP(z)$ of the kernel $\kappa$ carrying m’s draw ancestors to m, over their joint law $P$ — the same integral as kchain. lawof(m) requires m’s %mass to be %normalized (see total-mass classes); any other settled class — %null, %finite, %locallyfinite, or %unknown — is a static error, since an unnormalized measure is not its own law. %deferred is an inference state, not a total-mass class, and triggers no error; an engine that admits a %deferred-mass argument assumes normalization rather than proving it, and must leave the result’s %mass %deferred. lawof never normalizes its argument; normalize(m) states that intent. On a non-nullary kernel, lawof lifts pointwise, as the uniform kernel extension does for measure-algebra operations.

Trace of the reified law. A reified measure or kernel carries its traced sub-DAG as part of its value; a stochastic node shared between the traces of several joint components enters the composed trace once.

Phase of the reified law. Although the ancestor subgraph of the argument of lawof will typically include stochastic nodes, the resulting measure is itself deterministic (of parameterized or fixed phase): lawof absorbs stochasticity into the reified law rather than propagating it outward. Thus functionof can reify subgraphs that include stochastic nodes as long as they are reified to measures (see below).

4.8.2 Reification to functions and kernels #

When called with a single argument (without boundary specifications, see below), functionof traces the ancestor subgraph of its argument back to all leaves of parametric phase — that is, all elementof leaves. These leaf nodes become the inputs of the reified callable (a function or a kernel). Fixed ancestors (including external(...) and load_data(...) nodes) are resolved to their values and do not become inputs.

functionof can be called with additional keyword arguments to designate and label boundary nodes, stopping the graph trace there — so these nodes become, under their new names, the inputs of the resulting function or kernel. Boundary inputs themselves may be of parametric or stochastic phase, but not fixed phase. functionof effectively substitutes each boundary node a with an input node elementof(valueset(a)) under the given name.

FlatPPL has no closures: a reified callable captures no enclosing environment, and a fixed ancestor is resolved to its value rather than retained as a binding.

Referential transparency is a core property of FlatPPL. This requires that the sub-graph to be reified by functionof must not contain stochastic nodes that are not reified to measures. This means that the sub-graph must not contain draw nodes and nodes derived from them which are not ancestors of lawof nodes in that subgraph (since lawof absorbs stochastic phase). lawof nodes in the sub-graph of a functionof only operate within that subgraph, including marginalization.

Consider a simple deterministic computation:

c = a ^ 2
d = max(b, 1.5)
e = c * d

Here e is a specific value during any given evaluation of the code. But the computation that produces e from a and b is useful in its own right: we might want to apply it elementwise over arrays, or use it as a transformation in pushfwd. The name e refers to a value, not to the computation that produced it, so we need a way to extract the computation as a first-class function:

f = functionof(e)                   # f: {a, b: Real} → Real
C = broadcast(f, a = A, b = B)      # apply f elementwise over arrays A, B

functionof(e) captures the entire computation leading to e — the sub-DAG that contains e and all its ancestors — as a reusable function object. The sub-DAG must be fully deterministic and so must not contain any draw nodes.

The argument names of the resulting function are the names of the leaf nodes of the reified sub-DAG; the input nodes of the function are decoupled from these leaf nodes. Fixed ancestor nodes are resolved to their values and not exposed as inputs. As the graph nodes are not ordered, the function only supports keyword arguments, not positional arguments.

The output type of the reified function matches the type of the argument of functionof:

f = functionof(e)                                    # scalar output
f = functionof(record(x = something, y = other))     # record output
f = functionof([something, other])                    # array output
f = functionof((something, other))                   # tuple output

Boundary inputs may also be tuples (functionof(..., t = some_tuple_expr)), in which case the reified function takes a tuple argument.

Specifying reification boundaries. Sometimes only a selected part of the ancestor sub-DAG should be reified. In our example, e depends on c and d, which in turn depend on a and b. If we want the function represented by the subgraph that starts at a and d — ignoring how d was computed from a and b — we can specify boundary inputs that stop the ancestor backtrace early:

g = functionof(e, p = a, q = d)     # g: {p, q: Real} → Real
M2 = pushfwd(g, some_measure)       # transform a measure over (p, q)

The keyword arguments p = a, q = d declare that the trace stops at nodes a and d, which become the inputs of g under the new names p and q. The computation from a and b to d is excluded — g only contains the path from a and d to e.

Boundary input specification is all-or-none: either every reified input is specified explicitly, or none is. Boundary input names must be distinct — a repeated name is a static error, which likewise forbids a lambda or named function from repeating an argument name. With explicit boundary specification, the reified function supports positional arguments in addition to keyword arguments, with positional order determined by the order in which boundary inputs are specified. Without a boundary specification, inputs are traced back to the parameterized-phase ancestor leaves of the reified expression (i.e. elementof nodes). Fixed-phase ancestors (e.g. external and load_data) are resolved to their values instead. The reified function then only supports keyword arguments, as no argument order can be inferred. A specified boundary node a can be thought of as being substituted with a new node, generated via elementof(valueset(a)), in the reified graph. Substitution applies to all boundary nodes before the ancestor trace runs, so a boundary node whose only paths to the output pass through another boundary node becomes disconnected from the output in the substituted graph — the resulting callable is constant in that input. This is permitted, not an error, and is what enables hierarchical-model composition (see Kernels and kernelof).

The function argument names do not have to differ from the boundary node names:

h = functionof(e, a = a, d = d)     # h: {a, d: Real} → Real

The resulting function h now has arguments named a and d, but these are local to the function and decoupled from the original nodes a and d.

Identity law. functionof(f(a, b), ..., a = a, b = b, ...) is equivalent to f.

Lambda notation. A lambda function is a shorthand notation for functionof with placeholders. A single-argument lambda is written arg -> expr; two or more arguments are listed in parentheses: (arg1, arg2, ...) -> expr (the single-argument parenthesised form (arg) -> expr is not valid).

Either form resolves to functionof(expr', arg1 = _arg1_, ...), where expr' is expr with every free occurrence of each arg_i rewritten to the placeholder _arg_i_. Inside the body, arg_i refers to the lambda’s input, not to any module binding of the same name. There is no nullary lambda. A lambda body must not itself be a lambda: a curried form such as x -> y -> x + y is a static error under the placeholder scoping rule.

For example, x -> 2 * x + 1 is equivalent to functionof(2 * _x_ + 1, x = _x_), and (x, y) -> x * y + 1 is equivalent to functionof(_x_ * _y_ + 1, x = _x_, y = _y_).

4.8.3 Kernels and kernelof #

kernelof(x, kwargs...) reifies (typically stochastic) value nodes to Markov kernels. x must not be a measure. kernelof(x, kwargs...) is equivalent to functionof(lawof(x), kwargs...) interpreted within the reified subgraph delimited by kwargs — the boundary substitution applies before the inner lawof is interpreted, so an enclosing kernelof boundary scopes what the inner lawof marginalizes over.

Identity law. kernelof(draw(K(a, b, ...)), a = a, b = b, ...) is equivalent to K. Equal output laws do not make kernels interchangeable as joint components: a joint of two reifications of one draw is the singular diagonal at every input, while a joint of two constructor kernels contributes a fresh coordinate per occurrence.

Consider this Bayesian example:

theta1 ~ Normal(mu = 0.0, sigma = 1.0)
theta2 ~ Exponential(rate = 1.0)
a = 5.0 * theta1
b = abs(theta1) * theta2
obs ~ iid(Normal(mu = a, sigma = b), 10)

joint_model = lawof(record(theta1 = theta1, theta2 = theta2, obs = obs))
prior_predictive = lawof(record(obs = obs))
prior = lawof(record(theta1 = theta1, theta2 = theta2))
forward_kernel = kernelof(record(obs = obs), theta1 = theta1, theta2 = theta2)

Here we define

• joint_model: the joint probability distribution over parameters and observation. joint_model is equivalent to jointchain(prior, forward_kernel).

• prior_predictive: the probability distribution of the observation obtained by marginalizing over theta1 and theta2 — they are internal stochastic nodes in the traced sub-DAG, not boundary inputs, so lawof integrates them out. prior_predictive is equivalent to kchain(prior, forward_kernel).

• prior: the probability distribution of the parameters theta1 and theta2.

• forward_kernel: the Markov kernel of the forward model; it maps values for theta1 and theta2 to probability distributions of the observation.

The same four objects can be expressed equivalently in pure measure algebra, without draw or lawof:

theta1 = elementof(reals)
theta2 = elementof(posreals)
a = 5.0 * theta1
b = abs(theta1) * theta2
obs_dist = iid(Normal(mu = a, sigma = b), 10)

prior = joint(theta1 = Normal(mu = 0.0, sigma = 1.0),
              theta2 = Exponential(rate = 1.0))
forward_kernel = functionof(obs_dist)
joint_model = jointchain(prior, forward_kernel)
prior_predictive = kchain(prior, forward_kernel)

Here forward_kernel is built directly from the measure-valued expression obs_dist via functionof (see above), and joint_model and prior_predictive are assembled from prior and forward_kernel using measure combinators.

Reification with interdependent boundary nodes. Reification places no constraint on the DAG-dependency structure among the chosen boundary nodes: a boundary may be an ancestor or descendant of another. Substitution (see reification boundaries) replaces every designated node with a fresh independent input before the ancestor trace runs, so original dependencies between boundaries are erased and the reified callable takes all of them as independent inputs. A boundary whose only paths to the output went through another boundary then has no occurrence in the lowered body — the callable is constant in that input.

The eight-schools model (D. Rubin, 1981) is a typical instance:

mu ~ Normal(0, 5)
tau ~ normalize(truncate(Cauchy(0, 5), interval(0, inf)))
theta ~ iid(Normal(mu, tau), J)
y ~ Normal.(theta, std_errs_data)

prior = lawof(record(mu = mu, tau = tau, theta = theta))
forward_kernel = kernelof(record(y = y), mu = mu, tau = tau, theta = theta)
joint_model = jointchain(prior, forward_kernel)

After substitution, mu, tau, theta are independent inputs of forward_kernel. In this specific example y depends only on theta, so the output of forward_kernel does not depend on mu or tau.

Reification and module scope. functionof and kernelof reify within the current module only: a parameterized value reached through a loaded-module reference cannot become an input — neither by the automatic trace nor as an explicit boundary node — so such a reification is a static error. A loaded module’s callables and fixed values may be used in the reified DAG (applied, or referenced and resolved to their values); only taking a cross-module parameterized value as an input is disallowed.

Note that lawof reifies a measure, which has no input list, so it is unrestricted — a measure may reference cross-module values, keeping the identity law intact across module boundaries. A reified measure that has a parametric dependency on a node defined in another module cannot then be reified to a kernel, due to the restriction above.

4.9 Interface adaptation #

FlatPPL provides relabel for structural renaming of outputs. At the value level, relabel assigns or renames fields on scalars, arrays, records, and tables:

v = relabel([1.0, 2.0, 3.0], ["x", "y", "z"])
# equivalent to:
v = record(x = 1.0, y = 2.0, z = 3.0)

and renames record fields and table columns:

v = relabel(record(a = 1.0, b = 2.0, c = 3.0), ["x", "y", "z"])
# equivalent to:
v = record(x = 1.0, y = 2.0, z = 3.0)

and wraps a scalar into a single-field record:

v = relabel(1.0, ["x"])
# equivalent to:
v = record(x = 1.0)

The same output-side renaming lifts directly to sets, functions, measures, and kernels:

named_S = relabel(cartpow(reals, 3), ["x", "y", "z"])
named_f = relabel(f, ["x", "y", "z"])
named_M = relabel(M, ["x", "y", "z"])
named_K = relabel(K, ["x", "y", "z"])

For functions, relabel(f, names) is post-composition with relabel on the function result; for measures it is equivalent to pushfwd(fn(relabel(_, names)), M); for kernels it acts on the output measures.

See built-in functions for full reference documentation on relabel.

4.10 Function composition and annotation #

fchain(f1, f2, f3, ...) composes deterministic functions left-associatively: fchain(f1, f2, f3)(x) equals f3(f2(f1(x))).

fchain combines well with auto-splatting: if f1 returns a record and f2 accepts keyword arguments matching the record fields, the two functions compose directly. fchain is the deterministic analogue of kchain.

bijection(f, f_inv, logvolume) annotates a function f with its inverse f_inv and the log-volume-element logvolume of the forward map. The result is semantically identical to f, but engines can use the inverse and volume element when computing densities of pushforward measures. logvolume may be a function or a scalar (0 for volume-preserving maps). See pushfwd for examples.

4.11 Placeholders and holes #

4.11.1 Placeholder variables #

Creating functions and kernels with boundary inputs via functionof and kernelof requires the creation of unique global variable names. Placeholder variables are special variable names of the form _name_ (leading and trailing underscore) that are local to a functionof(...) or kernelof(...) and can be thought of as implicitly creating a unique global input via elementof(anything). All placeholders must appear both in the expression to be reified and the boundary input keyword arguments.

For example

f = functionof(_a_ * b + _c_, a = _a_, c = _c_)

is equivalent to:

_tmp1 = elementof(anything)
_tmp2 = elementof(anything)
f = functionof(_tmp1 * b + _tmp2, a = _tmp1, c = _tmp2)

Placeholders are not holes (see below). An expression with placeholders like _a_ * b + _c_ must not appear outside of a functionof(...) or kernelof(...).

Scoping rule. The scope of a placeholder is the nearest enclosing functionof or kernelof. The same placeholder name may appear in different scopes without conflict:

functionof(functionof(_a_ * b, a = _a_)(some_value) + _a_, a = _a_)

A placeholder in an inner functionof or kernelof must be bound there, so this code is invalid:

# DISALLOWED:
functionof(functionof(_a_ * b + _c_, a = _a_)(some_value) + _d_, c = _c_, d = _d_)

4.11.2 Holes and fn #

The reserved name _ denotes a hole — a position in a deterministic expression where an argument is not yet supplied. Holes are only valid inside the special operation fn(...), which delimits the scope of hole lowering. The form fn(expr) wraps a hole expression and produces an anonymous function whose parameters are the holes in expr, in strict left-to-right reading order. This is analogous to the $f(\cdot, b)$ notation used in mathematics to denote a function with a free argument.

Each _ introduces a distinct positional parameter, named arg1, arg2, … in left-to-right reading order. These names are normative and may be used as keyword arguments. Holes do not inherit keyword names from enclosing call positions.

Note: Holes work differently than placeholders (see above).

A single hole, resulting in a one-argument function:

neg = fn(0 - _)
poly = fn(polynomial(coefficients = cs, x = _))

The trivial case fn(_) is the identity function, equivalent to the built-in identity.

Multiple holes — left-to-right positional order:

g = fn(f(_, b, _))
h = fn((_ / _) ^ 2)

Each _ is distinct: fn(_ * _) multiplies two different inputs rather than squaring one. Use placeholders if arguments need to appear in the expression more than once, e.g. functionof(_x_ * _x_, x = _x_).

Lowering. fn(expr) lowers to a functionof with placeholder variables. For example

g = fn(f(_, b, _))

lowers to

g = functionof(f(_arg1_, b, _arg2_), arg1 = _arg1_, arg2 = _arg2_)

which in turn lowers to

_tmp1 = elementof(anything)
_tmp2 = elementof(anything)
g = functionof(f(_tmp1, b, _tmp2), arg1 = _tmp1, arg2 = _tmp2)

4.12 Broadcasting #

broadcast(f_or_K, name = array, ...) or broadcast(f_or_K, array, ...) maps a function or kernel elementwise over arrays (and row-wise over tables; see tables). Keyword arguments bind inputs by name. If the callable has a declared positional order, positional binding is also permitted.

Dot-syntax. As a concise shorthand for broadcast, FlatPPL provides dot-call notation f.(args) and dot-operator notation a .op b, following the elementwise dotted operators of MATLAB and the broadcasting dot syntax of Julia. f.(<args>) lowers to broadcast(f, <args>) and supports calls with positional and keyword arguments. So Normal.(means, sigmas) is syntactic sugar for broadcast(Normal, means, sigmas). A dotted binary operator a .op b lowers to broadcast(opfn, a, b) and a dotted unary operator .op x to broadcast(opfn, x), where opfn is the function the plain operator lowers to. So A .+ B lowers to broadcast(add, A, B) and .! X lowers to broadcast(lnot, X). The following examples show both explicit broadcast calls and the equivalent dot-notation.

Deterministic broadcast with a named function:

b = 2 * a + 1
f = functionof(b, a = a)
C = broadcast(f, a = A)

b = 2 * a + 1
f = functionof(b, a = a)
C = f.(a = A)

With positional argument binding:

C = broadcast(f, A)

C = f.(A)

Using an anonymous function:

C = broadcast(fn(2 * _ + 1), A)

C = fn(2 * _ + 1).(A)

Multi-input broadcast:

d = a * x + b_param
g = functionof(d)
E = broadcast(g, a = slopes, x = points, b_param = intercepts)

d = a * x + b_param
g = functionof(d)
E = g.(a = slopes, x = points, b_param = intercepts)

Stochastic broadcast — kernel over array, producing an array-valued measure:

K = fn(Normal(mu = _, sigma = 0.1))
D ~ broadcast(K, A)

K = fn(Normal(mu = _, sigma = 0.1))
D ~ K.(A)

Return type:

• broadcast(function, ...) returns an array value.
• broadcast(kernel, ...) returns an array-valued measure: the independent product measure of the kernel applications at each array position.

The stochastic case returns a single product measure, not an array of measures. This respects the rule that measures are not stored inside arrays or records while still enabling vectorized stochastic model building.

Independence is explicit: Kernel broadcast means independent elementwise lifting. It does not cover dependent sequential kernels, autoregressive chains, or coupled array structure. For those, use jointchain or kchain with explicit indexing.

Collection arguments: FlatPPL does not automatically insert leading or trailing dimensions for array arguments, unlike some other languages. It does, however, automatically expand singleton dimensions: All collection arguments (arrays and tables) must have the same number of axes. Tables count as having one axis (the table’s rows) here. Along each axis, all collections must have the same size or be singular (size one). Size-one array axes are implicitly expanded by repetition to match the size of the other collection arguments along these axes. A size-one axis expanded against a zero-length axis yields length 0. addaxes (see array operations) may be used to reshape all input arrays to the same number of axes.

For example, given a function f, a matrix A and a vector b

C = broadcast(f, A, addaxes(b, 1, 0))

behaves like NumPy-style broadcasting, while

C = broadcast(f, A, addaxes(b, 0, 1))

behaves like Julia-style broadcasting.

Non-collection inputs: Scalar values, functions, kernels, measures and likelihood objects are allowed as broadcasting inputs, they are simply not iterated over but held constant while collection arguments are iterated over. If there are no collection arguments, broadcast behaves like a single function or kernel call.

Disallowed inputs: Records and tuples are not allowed as inputs of broadcasts.

Tuple-returning callables: if f returns a tuple, broadcast(f, ...) returns a tuple of arrays (componentwise), not an array of tuples.

broadcasted(f) returns a callable that is equivalent to applying broadcast to f — that is, broadcasted(f)(args...) is equivalent to broadcast(f, args...).

4.13 Reductions #

reduce(f, xs) is a fold over xs using the binary associative function f, called positionally as f(acc, next). For a vector xs = [x1, x2, ..., xn], it computes f(...f(f(x1, x2), x3)..., xn). For a table, xs is traversed row-wise and f takes two records (the accumulator and the next row). The first element (or row) is used as the initial accumulator value, so xs must be non-empty and f must return a value of the same type as the elements (or rows) of xs. Since f is required to be associative, implementations may evaluate in parallel. Unlike broadcast, reduce accepts only a deterministic function f, not a kernel.

scan(f, init, xs) is a left scan over xs using the binary function f with explicit initial accumulator init. f is called positionally as f(acc, next): the running accumulator (of the type of init) and the next element of xs (or row record, for tables), returning the new accumulator value of the same type as init. Output entry i is f(...f(f(init, x1), x2)..., xi), i.e. the accumulator after consuming xs[i]; the result has the same length (or row count) as xs and does not include init itself. f must be a deterministic function, not a kernel.

4.14 Multi-axis aggregation #

aggregate(f_reduction, output_axes, expr) generalizes vector reductions to multi-axis tensor contraction, as a generalization of Einstein summation.

aggregate evaluates expr at every combination of values of its named axes and reduces the resulting scalars by f_reduction along the axes that do not appear in output_axes, yielding an array of the shape declared by output_axes.

• f_reduction: an order-invariant vector-to-scalar reduction — i.e. a function $f: S^n \to S$ where f_reduction(v) is invariant under permutations of v. The eligible built-ins are sum, prod, mean, var, std, maximum, minimum, median, lany and lall.
• output_axes: an axis list of distinct axis names [.name1, .name2, ...] listing the retained axes in output order. Repeated names are a static error. The empty axis list [] is legal and denotes full reduction to a scalar.
• expr: an expression in which array indexing may contain axis names, like A[.i, 1, .j] or get(A, .i, 1, .j). Every axis name in output_axes must occur at least once in expr; any further axis names occurring in expr are reduced over with f_reduction. All array dimensions indexed with the same axis name must have the same length.

Examples:

A = rowstack([[1, 3, 5], [9, 5, 1]])
B = rowstack([[1, 0], [0, 1], [1, 1]])

# Matrix multiplication
C = aggregate(sum, [.i, .k], A[.i, .j] * B[.j, .k])
# → C = [[6, 8], [10, 6]]

# Weighted sum of squared differences, reducing over .j
w = [1, 2, 1]
D = aggregate(sum, [.i, .k], (A[.i, .j] - B[.j, .k])^2 * w[.j])
# → D = [[34, 25], [114, 113]]

# Column-wise variance of a matrix
V = aggregate(var, [.j], A[.i, .j])
# → V = [32, 2, 8]

# Row-wise sum with one fixed column
S = aggregate(sum, [.i], A[.i, 1])
# → S = [1, 9]

# Product over .j of (A + B) entries: prod-reduction over two matrices
P = aggregate(prod, [.i, .k], A[.i, .j] + B[.j, .k])
# → P = [[36, 24], [100, 108]]

Axis names are lexically scoped to the enclosing aggregate(...) and are not values; see axis names for the surface rules.

:= notation. As a shorthand for sum-aggregate, FlatPPL provides result[.name1, .name2, ...] := expr, equivalent to result = aggregate(sum, [.name1, .name2, ...], expr). The bracketed axis list may be empty for full reduction to a scalar.

So

D[.i, .k] := (A[.i, .j] - B[.j, .k])^2 * W[.j]

lowers to

D = aggregate(sum, [.i, .k], (A[.i, .j] - B[.j, .k])^2 * W[.j])

and the scalar (full-reduction) case

s[] := A[.i] * B[.i]

lowers to

s = aggregate(sum, [], A[.i] * B[.i])

aggregate composes cleanly with functionof as the namespace of axis names is local to the enclosing aggregate and the namespace of placeholders is local to the enclosing functionof:

mymatmul = functionof(
    aggregate(sum, [.i, .k], _A_[.i, .j] * _B_[.j, .k]),
    A = _A_, B = _B_
)

Relationship to broadcasting. Aggregation overlaps with broadcasting when no reduction takes place. For example, given

v = some_vector
A = some_matrix  # with first axis of same length as v
B = addaxes(v, 0, 1)

an aggregation (using singleton-axis indexing)

aggregate(f_reduction, [.i, .j], A[.i, .j] * B[.i, !])

is equivalent to

broadcast((a, b) -> a * b, A, B)

for every eligible f_reduction that is the identity on a one-element input; var and std are undefined over a single element, and lany and lall require boolean input.

4.15 Metric-aware Einstein summation #

metricsum(metric, output_axes, expr) is a metric-aware variant of sum-aggregate for tensor expressions with upper (contravariant) and lower (covariant) indices.

Variance-marked axis names are required inside metricsum: .<name>^ denotes an upper-index and .<name>_ denotes a lower-index axis (see axis names). Axis-names with variance markers are lexically scoped to the enclosing metricsum.

All-contravariant canonical storage. Outside metricsum, arrays carry no tensor metadata and are interpreted to contain the elements of the all-contravariant tensor. So a vector v outside of metricsum represents $v^\mu$, and v[.mu^] inside of metricsum accesses the stored vector entries directly. Likewise the array metric stores the elements of the contravariant tensor $g^{\mu\nu}$. So within metricsum, the covariant vector v[.mu_] accesses v transformed via the metric: $v_\mu = g_{\mu\nu} v^\nu$, where $g_{\mu\nu}$ is equivalent to the contents of inv(metric).

The metric argument is interpreted as upper-upper $g^{ij}$, consistent with the rule above. It must be a square, symmetric, and invertible rank-2 array. Lower-index access of the metric (metric[.i_, .j_]) denotes the inverse $g_{ij}$; mixed access (metric[.i^, .j_]) denotes the Kronecker delta $\delta^i{}_j$.

Output variance specifies the variance of the components computed by expr, not the stored layout. The actual result returned by metricsum is automatically raised to all-upper canonical storage. Reading the result later inside another metricsum with the same metric recovers the originally defined components.

:= notation. As a shorthand, metric: result[output_indices...] := expr lowers to result = metricsum(metric, [output_indices...], expr).

Expression restrictions. metric itself and all arrays indexed with co-/contravariant axis names in expr must be arrays of scalars. expr must produce scalar values for all combinations of axis index values.

Static checks. Every repeated non-output index in expr must occur exactly twice — once upper and once lower; every output index must occur in expr with the same variance and may not also be contracted; bare neutral aggregate axes (.i without a variance marker) are not allowed inside metricsum. A non-output index occurring once is legal and is summed over by the lowering below. Its variance is therefore semantically significant, and the result is coordinate-dependent rather than tensorial.

Equivalence to aggregate under identity metric. metricsum(eye(n), ...) is equivalent to an aggregate(sum, ...) with co-/contravariant axis names replaced by aggregate axis names.

Lowering to aggregate. Each _ (lower-variance) axis name in expr becomes an inv(metric) contraction; each _ output axis becomes a metric contraction after the sum, raising the result to all-upper canonical storage. The whole lowering may be expressed as a single aggregate(sum, ...) with metric factors inlined, or as a chain that precomputes mixed-variance intermediates as common subexpressions; both forms are semantically equivalent.

Example. Composition of three Lorentz transformations $L^\mu{}_\rho = L_1{}^\mu{}_\nu \, L_2{}^\nu{}_\sigma \, L_3{}^\sigma{}_\rho$ as (1,1)-tensors, under metric g:

g: L[.mu^, .rho_] := L1[.mu^, .nu_] * L2[.nu^, .sigma_] * L3[.sigma^, .rho_]

Equivalent non-shorthand form:

L = metricsum(g, [.mu^, .rho_],
    L1[.mu^, .nu_] * L2[.nu^, .sigma_] * L3[.sigma^, .rho_])

Lowering to aggregate, inline metric insertion:

__g_down = inv(g)
L = aggregate(sum, [.mu, .rho_up],
    L1[.mu, .a] * __g_down[.a, .nu] *
    L2[.nu, .b] * __g_down[.b, .sigma] *
    L3[.sigma, .c] * __g_down[.c, .rho] *
    g[.rho, .rho_up]
)

Lowering to aggregate, precomputed mixed-variance intermediates:

__g_down = inv(g)
__L1_mixed = L1 * __g_down
__L2_mixed = L2 * __g_down
__L3_mixed = L3 * __g_down
__L_mixed = aggregate(sum, [.mu, .rho],
    __L1_mixed[.mu, .nu] * __L2_mixed[.nu, .sigma] * __L3_mixed[.sigma, .rho]
)
L = __L_mixed * g

4.16 Lowered linear form #

A FlatPPL module admits a stable lowering to a linear SSA-style core form in which every non-atomic subexpression is bound to an auto-generated unique name (see Placeholders and holes for the lowering stages). In the resulting form, every binding’s right-hand side is either a literal or a function call: name = c or name = f(name, ...). Operators, indexing, field access, and array literals all desugar to function calls (add, get, vector, etc.), giving the core form a uniform shape. This is a semantic property of FlatPPL modules; it is independent of the surface syntax.

4.17 Standard modules #

FlatPPL’s core built-ins aim to be domain-neutral. Functionality and vocabulary more specific to particular disciplines is provided via standard modules that are not part of base.

Standard modules are versioned independently from FlatPPL itself, and loaded via standard_module(name, compat). The names of standard modules are official and coordinated; compat is a compatibility version string with the same semantics as flatppl_compat (see below). standard_module only accepts positional arguments, not keyword arguments.

For example:

hepphys = standard_module("particle-physics", "0.1")
peak = hepphys.CrystalBall(m0 = 125.0, sigma = 2.0, alpha = 1.5, n = 3.0)

See section Standard modules for the current set of FlatPPL standard modules.

Engine support. FlatPPL engines are not required, but encouraged, to implement all standard modules if feasible on the given platform.

4.18 Module composition #

A FlatPPL module is a flat namespace of named bindings (see above). Modules can load other modules written in FlatPPL though, to make models composable:

load_module(source) loads a FlatPPL file and returns a module reference.

source may be a file path or a URL (see Remote file caching).

Each load_module call instantiates the loaded module independently: two calls share no nodes, even with identical source and substitutions, so reified callables from different calls have disjoint stochastic ancestors. To share one instance, bind the module reference once and reuse the name.

In the canonical syntax, bound names in the loaded module are accessed via dot syntax:

sig_module = load_module("signal.flatppl")
bkg_module = load_module("background.flatppl")

sig_model = sig_module.model
bkg_model = bkg_module.model

Access to loaded modules is not transitive: The loading module may access names in the loaded module, but not names in modules loaded by that module, so sig_module.model is valid but sig_module.some_other_module.some_name is not.

Load-time substitution. load_module may also be called with keyword arguments to substitute explicit input nodes of the loaded module:

sig_module = load_module("signal_channel.flatppl", mu = signal_strength, theta = nuisance)

The left-hand side of each keyword argument must refer to an input of the loaded module. The phase of this input determines what it can be bound to on the right-hand side:

• external inputs of the loaded module may only be bound to fixed values in the loading module.
• elementof inputs of the loaded module may only be bound to parameterized values in the loading module.
• No other kinds of nodes in the loaded module may be bound to nodes in the loading module.

Value sets must be compatible in both cases, so the computational structure of the loaded module is not modified.

Stochastic boundary. Only bindings of fixed or parameterized phase in the loaded module are accessible from the loading module (module.name). Bindings of stochastic phase — direct draws or values with draw ancestors that have not been reified via lawof/kernelof — are invisible to the loading module. This preserves referential transparency and avoids semantic ambiguity when two loaded modules load a common third module.

Path resolution. Relative file paths in load_module(...) are resolved relative to the directory of the FlatPPL file containing that load_module(...) call, not the host process’s working directory. For embedded FlatPPL code, relative paths are resolved relative to the directory of the source file containing the embedded FlatPPL code block. The forward slash / is the mandatory path separator on all platforms. Parent-directory traversal via .. is allowed. Absolute file paths are permitted but discouraged, as they prevent relocatable model repositories.

Aliasing is just assignment: sig_model = sig_module.model creates a local alias — a reference to the same underlying object in the loaded module’s DAG, not a clone.

Bundles. source may be a bundle holding a main FlatPPL module file and its module and data dependencies. If source is a directory, load_module loads its root main.flatppl (which will typically itself use load_module and load_data to load dependencies located under that directory). If source is a ZIP file (extension .zip; .flatppl.zip recommended where practical), it loads main.flatppl from the archive root. If there is no root main.flatppl in the archive, but the archive’s sole top-level entry is a directory containing a main.flatppl, it is loaded from there. A missing main.flatppl is an error.

Within a bundle (directory or ZIP), relative paths — in both load_module and load_data — resolve only inside the bundle and must not escape its root via ...

4.19 FlatPPL version compatibility #

A FlatPPL module may declare which versions of FlatPPL it is compatible with via the reserved binding flatppl_compat:

flatppl_compat = "0.1"

The value is a string following Julia-style semantic versioning conventions: for pre-1.0 versions, the minor version is breaking ("0.1" means $\geq$ 0.1.0, $<$ 0.2.0); for versions $\geq$ 1.0, the major version is breaking ("1" means $\geq$ 1.0.0, $<$ 2.0.0). Multiple ranges are comma-separated and combined with OR:

flatppl_compat = "0.8, 0.9.2, 1.0.0, 2"

declares compatibility with versions v0.8.x, v0.9.2 up to (excluding) v0.10, v1.x.y, and v2.x.y.

The declaration is optional. Short-lived models, didactic examples and the like may omit it. For embedded FlatPPL blocks, version compatibility may be managed at the host-language level (e.g. via Python or Julia package/environment dependency version bounds on FlatPPL packages). FlatPPL files of models intended for long-term use, publication or archival should definitely include a compatibility declaration.

The compatibility declaration of a loaded module is accessible via dot syntax (like any other bound value in the module): some_module.flatppl_compat.

4.20 Code documentation #

FlatPPL treats documentation as a first-class property of bindings, not as ambient commentary. Doc-comments in surface FlatPPL (%, %%%; see section Syntax) attach to bindings and are preserved when lowering to FlatPIR.

Default markup and markup tags. By default, documentation is written in Markdown. The markup language can be explicitly selected with a markup tag. FlatPPL tooling should know how to handle the following tags:

• %md/%%%md: GitHub-Flavored Markdown with math ($...$, $$...$$)
• %typ/%%%typ: Typst

Markup tags should reflect the canonical file extension of the markup language.

Attachment. A doc-comment attaches to at most one binding, in one of two positions:

• Leading: before a FlatPPL binding (only whitespace, newlines and ; statement separators between).
• Trailing: a single-line % ... after a binding’s right-hand side, before the next newline or ; statement separator. Block %%% forms must not be in trailing position.

Each binding may carry at most one doc-comment (leading or trailing, not both). Two leading % ... lines on the same binding are an error — use a %%% block for multi-line content. A doc-comment that doesn’t attach to a binding is invalid code.

Module-level documentation. A doc-comment attached to the flatppl_compat binding serves to document the module itself:

%%%
# Eight-schools model

Hierarchical Normal model after Rubin (1981).
%%%
flatppl_compat = "0.3"

Comments are not documentation. Plain comments (#, ###) are discarded at parse time and do not appear in FlatPIR. They are author-eyes-only notes on the surface source. Anything intended to outlive the surface file — for tools, for downstream readers via FlatPIR, or for export to external systems — must use a doc-comment. (FlatPIR has its own ; line comments in the canonical text syntax, but those are reserved for tooling annotations and do not carry user-written surface comments.)

4.21 Remote file caching #

A load_module(url) or load_data(url) source may be an http/https URL rather than a local path. FlatPPL is meant to be supported by multiple engines and tools in a variety of host languages, and the design leaves a lot of freedom to individual FlatPPL implementations. But caching of remote content to local files should be consistent across various tools and engines, so they should adhere to the following conventions (subject to change in future FlatPPL versions) and thus use a shared local cache:

Cache directory. The main cache directory (referred to below as <flatppl-cachedir>) is set by the environment variable FLATPPL_CACHEDIR. If not set, the following default is used:

• Linux and BSD: A directory flatppl directly under the XDG cache directory (defaults to $HOME/.cache/flatppl).
• macOS: $HOME/Library/Caches/flatppl
• Windows: %LOCALAPPDATA%\flatppl

Layout and keys. Under <flatppl-cachedir>/v1/:

• objects/<kk>/<key>.<ext> holds the fetched URL content. key is the lowercase-hexadecimal SHA-256 of the request URL with any #-fragment removed (no other normalization), and <kk> is the first two characters of key. <ext> is the requested file’s full trailing extension — everything after the first . in the URL’s final path segment — so a multi-part extension is kept whole; a final segment with no . uses just <key>.
• objects/<kk>/<key>_meta.json is a mandatory JSON metadata file with fields url (the original URL), resolved_url (the URL fetched after any redirects), retrieved (ISO 8601 UTC time), content_type, and the HTTP validators etag and last_modified (any may be null if the server omits it). Readers ignore unknown fields.
• trust/<kk>/<key> is a per-URL trust marker — its presence means the URL is trusted — keyed by the same <kk>/<key> as the object. Trust is based on the original URLs, not on redirect URLs.
• tmp/ holds temporary files during downloads, must be on the same filesystem as objects/ to achieve atomic renames.

The cache is keyed by URL hash; there is no separate index, and objects/ is its complete state.

Resolve and fetch. To resolve a URL, FlatPPL implementations use the cached object under the matching hash if present. Otherwise, implementations check if the URL is trusted (see below), fetch the URL via the temporary directory (see below for details), and store the URL content in the objects directory under the URL hash with the file extension added. URL redirects are followed for download but the hash of the original URL is used in the cache. A fetch that fails — a network error, a final response whose status is not 2xx, or an unresolvable redirect — is an error, and nothing is written to the cache (no partial object and no metadata). If FLATPPL_CACHE_OFFLINE is set, a cache miss is an error and no URL fetch is attempted.

Trust. Before fetching a URL that has no trust marker, interactive tooling must obtain the user’s approval and then create its trust/<kk>/<key> marker. Non-interactive tooling must error if a requested URL is not marked as trusted. If the environment variable FLATPPL_TRUST is set, all URLs are trusted implicitly by interactive and non-interactive tooling, but no trust markers are created.

Atomicity and concurrent tools. Content is initially downloaded to the tmp/ directory, then fsynced and then atomically renamed into the destination file under the objects directory. The _meta.json file is written the same way and renamed into place before its content object, so a present <key>.<ext> always has its metadata. Trust markers are created with an exclusive O_CREAT | O_EXCL open. The whole cache is lock-free and <flatppl-cachedir> should be located on a file system that supports atomic renaming.

Environment variables.

5 Canonical syntax #

This section specifies the canonical surface form of FlatPPL, used throughout this document as a notation for defining FlatPPL semantics and presenting examples. It also defines a mechanism for embedding FlatPPL in Python and Julia (see host-language embedding below).

The semantics of FlatPPL do not depend on this canonical syntax. Alternative syntactical representations may be advantageous for specific software ecosystems and use cases. Such alternative representations must map directly to and from canonical FlatPPL without change in semantics.

Canonical FlatPPL uses the filename extension .flatppl.

5.1 Statements #

A FlatPPL module is a sequence of statements separated by newlines or semicolons (equivalent). One statement per line is the recommended style; semicolons exist primarily as a fallback for channels that may not preserve line breaks.

5.2 Comments #

# starts a line comment; ### alone on a line opens a block comment closed by a matching ### alone on a line. Both forms are discarded by the parser. Line comments are terminated by the next newline or ;, whichever comes first. Block fences may have leading horizontal whitespace (so embedded FlatPPL inside an indented host block is fine).

x = 3.14  # Inline comment.

###
Block comment.
###

For documentation that attaches to bindings and survives into FlatPIR, see below.

5.3 Documentation #

Doc-comments are lexically symmetric to plain comments (% ↔︎ #, %%% ↔︎ ###) but attach to bindings and survive into FlatPIR. Semantics, attachment rules, default markup, and module-level documentation are specified in Code documentation; this section covers the surface lexical forms only.

• Single-line: %[<markup>]? <content> runs to end of line or to the next ;, whichever comes first.
• Block: %%%[<markup>]? alone on a line opens; a matching %%% alone on a line closes. Content is verbatim, line by line. Fence lines may have leading horizontal whitespace, which is ignored; content lines are taken as written.
• Markup tag (optional, no space after the leading % / %%%): md (default, GitHub-Flavored Markdown with $...$ / $$...$$ math) or typ (Typst). An unrecognized tag is a parse error. The naming convention is “typical file extension”; additional tags may be added in future spec versions.

% Prior mean.
mu = 0

sigma = 1  % Prior std. dev.

%%%
The observation model: independent Gaussians, shared unknown scale.
%%%
obs ~ iid(Normal(mu, sigma), 5)

5.4 Supported constructs #

FlatPPL has a very lean syntax:

• Bindings: name = expr and decomposition a, b, c = expr (see below).
• Tilde bindings: name ~ expr and decomposition a, b ~ expr, equivalent to name = draw(expr) and a, b = draw(expr) respectively (see variates and measures).
• Named functions: f(arg1, arg2, ...) = expr is shorthand for binding f to a lambda (see named functions).
• Literals: numbers (3.14, 42, 0xF7, 0x3e, 1_000_000, 1.45e7), strings ("foo"), booleans (true, false), arrays ([1, 2, 3]), records (record(a = 1, b = 2)), tuples ((a, b)).
• Infix arithmetic, exponentiation, and comparisons: +, -, *, /, ^, unary -, and the comparisons <, >, ==, !=, <=, >=, and in (set membership). Comparisons may be chained: a < b <= c lowers to land(a < b, b <= c). ^ is right-associative and binds tighter than unary -.
• Logical operators: &&, ||, ! (lowering to land, lor, lnot; see Logic and conditionals).
• Broadcasting: Dot-call f.(...) and dot-prefixed operator application a .+ b are syntactic sugar for broadcast (see Broadcasting syntax).
• Lambda: arg -> expr (single arg) or (arg1, arg2, ...) -> expr (multi-arg) is shorthand for functionof with placeholders (see Lambda syntax).
• Aggregation: axis-indexed binding C[.i, .k] := expr is shorthand for sum-aggregate; the metric-prefixed form g: C[.mu^, .nu_] := expr is the metricsum shorthand for metric-aware Einstein summation. (See Axis names and aggregation.)
• Function calls: f(x, y) (positional), f(a = x, b = y) (keyword) and f(object, a = x, b = y) (for some special operations).
• Indexing and field access: A[i], A[i, j], A[:, j], r.field.

Also see the formal grammar below.

See binding names for rules on binding names, name resolution, and the reserved modules self and base.

5.5 Excluded constructs #

The above are the only syntactical constructs allowed in FlatPPL. The following Python and Julia constructs, for example, are not allowed directly in canonical FlatPPL, but can easily be represented in other ways:

• No type annotations. Types are inferred from the semantic rules.
• No loops or conditionals. Use ifelse(cond, a, b) for piecewise definitions (see logic and conditionals).
• No function definition blocks. A function body is a single expression (a lambda or f(x) = expr, see named functions); for named intermediate steps use functionof (see language design).
• No implicit operator broadcasting. Infix +, -, *, /, ^ and unary - follow standard linear-algebra and scalar semantics: + and - require operands of identical shape (scalars, or arrays of matching shape), * supports the products enumerated in operator-equivalent functions, / requires a scalar divisor, and ^ is scalar-only.

5.6 Decomposition syntax #

The left side of an = or ~ assignment may decompose an array, record, or tuple into named components:

a, b, c ~ MvNormal(mu = mean_vector, cov = cov_matrix)
# equivalent to
# a, b, c = draw(MvNormal(mu = mean_vector, cov = cov_matrix))

x, y = some_record
l, m, n = some_tuple
value, _ = rand(rstate, m)

Decomposition is by position. For records, the field order determines which value each name receives; for arrays and tuples, positional index does. This is syntactic sugar: it lowers to an assignment followed by indexed or field-access bindings.

5.7 Indexing and slicing #

FlatPPL uses 1-based indexing.

A[:, j] selects all elements along the first axis at fixed index j. This lowers to get(A, all, j), where all is a predefined selector meaning “entire axis.”

A[:, j]          # → get(A, all, j)
A[i, :]          # → get(A, i, all)
T[:, :, k]       # → get(T, all, all, k)
T[i, :, k]       # → get(T, i, all, k)

! is the only selector: it extracts the unique element of a length-1 axis. This lowers to get(A, only, ...). The indexed axis must have length one.

A[!, j]          # → get(A, only, j)
v[!]             # → get(v, only)

5.8 Special operations #

elementof(S), valueset(x), draw(M), lawof(x), functionof(...), kernelof(...), and fn(...) are special operations with their own syntax rules — they are not ordinary function calls. Their semantics are defined in language design. load_module(...) is documented in multi-file models.

5.9 Broadcasting syntax #

FlatPPL provides dot-prefixed shorthand for broadcast:

• f.(<args>) — dot-call.
• a .op b — dotted binary operators.
• .op x — dotted unary operators.

Each dotted operator has the same precedence as its plain counterpart. in has no dotted form (its right operand is a set, not a broadcastable value). See Broadcasting for dot-notation lowering.

5.10 Lambda syntax #

Lambda functions, denoted as arg -> expr and (arg1, arg2, ...) -> expr, are syntactic sugar for functionof. (arg) -> expr is not legal syntax (no parentheses around the argument in single-argument lambdas). At least one argument is required (no nullary lambdas).

The body extends as far right as possible — lambdas have lower precedence than every other expression form, so (a, b) -> a^2 + b^2 parses with a^2 + b^2 as the body. Inside the body, the argument names refer to the lambda’s inputs and shadow any module-level binding of the same name. See Reification to functions and kernels for the desugaring.

5.11 Named functions #

f(arg1, arg2, ...) = expr is syntactic sugar for binding f to a lambda — it desugars to f = (arg1, arg2, ...) -> expr. As for lambdas, at least one argument is required (f() = expr is not legal; bind a plain value with f = expr) and the argument names are local to the body. The defined function accepts positional and keyword calls, and f is a first-class value usable wherever a functionof result is.

f(x, y) = x^2 * y^2                          # equivalent to f = (x, y) -> x^2 * y^2
g(x, y, z) = record(p = x + y, q = y * z)    # record/array/tuple body → multi-output
h(x, y) = [x / y, x * log(y)]

The construct is purely a surface rewrite: it adds no FlatPIR node and inherits every function property — scoping, duplicate-argument rules, phase, doc-comment attachment — from the lambda desugaring. There is no tilde form, since a measure is not a function (f(arg1, ...) ~ expr is not legal).

5.12 Axis names and aggregation #

Axis names are written .<name> and are symbolic index labels used by aggregate. They are lexically scoped to the enclosing aggregation and are not values: an axis name is legal only as an entry in aggregate’s output_axes axis list, as an index inside [...] within the body, or as a binder on the left-hand side of :=. Used anywhere else it is a static error.

The aggregation form C[.i, .j, ...] := expr is shorthand for sum-aggregate; see there for the desugaring. The bracketed axis list may be empty (x[] := expr) for full reduction to a scalar.

Axes may carry a variance marker — .<name>^ (upper / contravariant) or .<name>_ (lower / covariant) — inside metricsum for metric-aware Einstein summation. The marker form metric: C[...] := expr is the shorthand for metricsum. Axis names themselves may not end in _, since trailing _ is reserved as the lower-variance marker.

5.13 Host-language embedding #

FlatPPL defines a recommended mechanism for embedding FlatPPL code in Python and Julia. Python and Julia implementations of FlatPPL should use this approach so that independent tooling, like FlatPPL grammars and extensions for code editors, can support it consistently.

Embedding FlatPPL in Python should be realized via a call to a function flatppl on a raw string, flatppl(r"""<FlatPPL code>"""), e.g.

model = flatppl(r"""
mu = elementof(reals)
x ~ Normal(mu = mu, sigma = 1)
""")

FlatPPL is not indentation-sensitive, so leading horizontal whitespace from an indented Python block does not affect the model.

Embedding FlatPPL in Julia should be realized via a Julia string macro flatppl"""<FlatPPL code>""", e.g.

model = flatppl"""
mu = elementof(reals)
x ~ Normal(mu = mu, sigma = 1)
"""

Host data enters an embedded model only through the explicit input mechanisms (external, load_data, and parameterized load_module); embedded FlatPPL does not capture host-language variables, and host string interpolation should not be used to splice values into a model, as that would bake constants into the graph instead of creating input nodes.

The only sequence that cannot appear inside embedded FlatPPL source — including inside doc-comments — is a literal """, which would terminate the host string. Markdown fenced code inside doc-comments should use triple-backticks (```).

5.14 Formal grammar #

The canonical surface syntax is defined in EBNF below (ISO 14977-style, with ::= for production and | for alternation).

Grammar.

(* Top level *)
Module          ::= StmtSep* (Statement (StmtSep+ Statement)*)? StmtSep* EOF
Statement       ::= Binding | TildeBinding | Decomposition | TildeDecomposition
                  | FunctionDefinition
                  | AggregateBinding | MetricsumBinding

(* Bindings *)
Binding            ::= Name "=" Expression
Decomposition      ::= Name ("," Name)+ "=" Expression
TildeBinding       ::= Name "~" Expression
TildeDecomposition ::= Name ("," Name)+ "~" Expression
FunctionDefinition ::= Name "(" Name ("," Name)* ")" "=" Expression
AggregateBinding   ::= Name AxisList ":=" Expression
MetricsumBinding   ::= Name ":" Name AxisList ":=" Expression

(* Expressions — lambda at top, logical OR/AND above comparisons,
   exponentiation below multiplicative *)
Expression      ::= Lambda | Or
Lambda          ::= LambdaParams "->" Expression
LambdaParams    ::= Name | "(" Name "," Name ("," Name)* ")"
Or              ::= And (("||" | ".||") And)*
And             ::= Comparison (("&&" | ".&&") Comparison)*
Comparison      ::= Additive (CompOp Additive)*       (* chained *)
Additive        ::= Multiplicative (AddOp Multiplicative)*
Multiplicative  ::= Unary (MulOp Unary)*
Unary           ::= ("-" | ".-") Unary | ("!" | ".!") Unary | Exponential
Exponential     ::= Postfix (("^" | ".^") Unary)?
Postfix         ::= Primary (FieldAccess | DotCall | Indexing | Call)*
Primary         ::= Literal | Name | Axis | AxisList | "(" Expression ")"

FieldAccess     ::= "." Name
DotCall         ::= "." "(" CallArgs ")"
Indexing        ::= "[" IndexExpr ("," IndexExpr)* "]"
IndexExpr       ::= Expression | ":" | "!"
Axis            ::= "." AxisName VarianceMarker?
AxisName        ::= Letter
                  | Letter (Letter | Digit | "_")* (Letter | Digit)
                    (* Identifier that must not start with "_" or end with "_";
                       trailing "_" is reserved as the lower-variance marker. *)
VarianceMarker  ::= "^" | "_"
AxisList        ::= "[" (Axis ("," Axis)*)? "]"

CompOp          ::= "<" | ">" | "==" | "!=" | "<=" | ">=" | "in"
                  | ".<" | ".>" | ".==" | ".!=" | ".<=" | ".>="
AddOp           ::= "+" | "-" | ".+" | ".-"
MulOp           ::= "*" | "/" | ".*" | "./"
ContinuationOp  ::= AddOp | MulOp | CompOp | "^" | ".^"
                  | "&&" | "||" | ".&&" | ".||"
                  | "->" | "=" | "~" | ":=" | ":"
                    (* Every infix binary operator, the lambda arrow, and the
                       binding operators. Only a TRAILING occurrence continues a
                       line (see "Statement separation"); a line beginning with
                       an operator starts a new statement, and is a parse error
                       unless that operator is unary. *)

(* Calls *)
Call            ::= "(" CallArgs ")"
CallArgs        ::= PositionalArgs | KeywordArgs | MixedArgs
PositionalArgs  ::= Expression ("," Expression)*
KeywordArgs     ::= KeywordArg ("," KeywordArg)*
KeywordArg      ::= Name "=" Expression
MixedArgs       ::= Expression ("," Expression)* ("," KeywordArg)+

(* Literals *)
Literal         ::= Number | String | Boolean | ArrayLiteral | RecordLiteral | TupleLiteral
Boolean         ::= "true" | "false"
Number          ::= IntegerLit | RealLit
ArrayLiteral    ::= "[" Expression ("," Expression)* ","? "]"
RecordLiteral   ::= "record" "(" KeywordArg ("," KeywordArg)* ","? ")"
TupleLiteral    ::= "(" Expression "," Expression ("," Expression)* ","? ")"

(* Lexical *)
Name            ::= (Letter | "_") (Letter | Digit | "_")*
IntegerLit      ::= DecIntLit | HexIntLit
DecIntLit       ::= Digit ("_"? Digit)*
HexIntLit       ::= "0x" HexDigit ("_"? HexDigit)*
HexDigit        ::= Digit | "a" .. "f" | "A" .. "F"
RealLit         ::= DecIntLit "." DecFracPart? Exponent?
                  | DecIntLit Exponent
                  | "." DecFracPart Exponent?
DecFracPart     ::= Digit ("_"? Digit)*
Exponent        ::= ("e" | "E") ("+" | "-")? DecIntLit
String          ::= '"' StringChar* '"'
StringChar      ::= any character except '"' and '\' | '\' EscapeChar
EscapeChar      ::= '"' | '\' | "n" | "t" | "r" | "0"
Letter          ::= "a" .. "z" | "A" .. "Z"
Digit           ::= "0" .. "9"
Newline         ::= LF | CR | CRLF
StmtSep         ::= Newline | ";"

(* Plain comments — discarded by the parser *)
LineComment     ::= "#" { any character except newline or ";" }
BlockComment    ::= HWS* "###" HWS* Newline
                    { any line whose trimmed content is not "###" }
                    HWS* "###" HWS* Newline

(* Doc-comments — attached to bindings; see "Documentation" *)
DocLine         ::= "%" MarkupTag? { any character except newline or ";" }
DocBlock        ::= HWS* "%%%" MarkupTag? HWS* Newline
                    DocBlockLine*
                    HWS* "%%%" HWS* Newline
DocBlockLine    ::= { any character except newline } Newline
                  ; a line whose trimmed content equals "%%%" closes the block
MarkupTag       ::= "md" | "typ"
HWS             ::= " " | "\t"            (* horizontal whitespace *)

Statement separation. Statements are separated by one or more newlines or semicolons; the two are fully equivalent. Blank lines and comment-only files are permitted. Newlines inside an unclosed ( or [ (paren/bracket depth > 0) are treated as whitespace (implicit line continuation), letting expressions span multiple lines:

rate = superpose(
    weighted(mu_sig * efficiency, signal_template),
    bkg_template
)

At paren/bracket depth 0, a newline is likewise treated as whitespace when the line’s last token is a ContinuationOp, so the statement continues on the next line that carries a token. A trailing line comment, and any blank or comment-only lines in between, do not end the continuation. A ^ or _ immediately following an axis name is that axis’s variance marker, not a ContinuationOp.

rate = mu_sig * signal_yield +      # trailing operator continues the line
       mu_bkg * bkg_yield

Note on MixedArgs. Syntactically, any Call may use MixedArgs (one or more leading positional expressions followed by one or more keyword arguments). Semantically, only the special operations functionof, kernelof, broadcast, load_module, and load_data accept this shape; other callables must use PositionalArgs or KeywordArgs only. Among these, functionof, kernelof, load_module, and load_data take exactly one leading positional argument; only broadcast accepts multiple positional arguments before the keyword arguments.

Note on reserved words. The keywords in, true, false, all, and only are recognized before Name and cannot be used as bindings. The top-level binding names inputs and outputs are reserved for the determinization signature.

Note on holes and placeholders. The lexical rule for Name admits _ (the hole used inside fn(...)) and trailing-underscore identifiers _x_ (placeholders used inside functionof/kernelof). The grammar parses both as ordinary Name; the syntactic restrictions on where they may appear are documented in functions and reification.

Note on axis names. The grammar admits Axis (.<name>) as a Primary, but Axis is legal only inside an aggregation — as an entry of aggregate’s output_axes, as an [...] index in its body, or as a binder of an AggregateBinding. Anywhere else it is a static error. The grammar likewise admits AxisList as a Primary, but it is legal only as the output_axes argument of an aggregate or metricsum call and as the axis-list binder of an AggregateBinding or MetricsumBinding; anywhere else it is a static error. Unlike ArrayLiteral, AxisList may be empty: aggregate(sum, [], expr) denotes full reduction to a scalar.

Note on tuples. (x) is a parenthesised expression. (x, y) is a tuple. The single-element form (x,) is not in the grammar — single-element tuples are not supported (consistent with the design rule that tuples have at least two components).

Note on parser disambiguation. The grammar is intended to be parsed with bounded lookahead. CallArgs is technically ambiguous in pure EBNF (an Expression can begin with a Name, and so can a KeywordArg), but a one-token lookahead after the leading Name (checking for =) suffices to choose between PositionalArgs/MixedArgs and KeywordArgs. After a ., a one-token lookahead distinguishes DotCall (. followed by () from FieldAccess (. followed by a Name). Dot-prefixed operators (.+, .^, .==, …) are single tokens, recognized by maximal munch. Maximal munch also resolves the trailing-dot real literal against a dotted operator: in 1./x the 1. is the real literal 1.0 (so this is 1.0 / x), as in Julia. A dotted operator on an integer-literal operand needs whitespace or an explicit fractional part — 1 ./ x or 1.0 ./ x. After a closing ), a one-token lookahead for -> distinguishes a Lambda from a parenthesised expression or tuple literal; a parenthesised lambda is well-formed only if the parenthesised content was a list of two or more bare Names. A bare Name immediately followed by -> is a single-argument Lambda. At statement level, a Name followed by ( begins a FunctionDefinition (no other statement form starts Name "("), disambiguated by checking that the closing ) is followed by =. A .Name token is FieldAccess when it follows a Postfix-able expression, and Axis otherwise (at the start of a Primary). Inside [...], a ! token followed immediately by , or ] is the only axis keyword; otherwise it is the unary logical-not operator starting an Expression. In AxisList’s legal positions (as above), [...] parses as AxisList, not ArrayLiteral.

6 Measure algebra and analysis #

This section documents the measure-level operations that form the compositional core of FlatPPL.

6.1 Measure-theoretic foundations #

A measurable space is a pair $(X, \Sigma_X)$ of a set and a $\sigma$-algebra. All spaces arising in FlatPPL are standard Borel spaces ($\mathbb{R}$, $\mathbb{Z}$, and finite products thereof), where the $\sigma$-algebra is the standard Borel $\sigma$-algebra and can be left implicit. A measure on $X$ is a $\sigma$-additive function $\mu: \Sigma_X \to [0, \infty]$. A probability measure has $\mu(X) = 1$. All measures in FlatPPL are $\sigma$-finite (admitting a countable cover of finite-measure sets), which ensures that product and marginalization operations are well-defined and that the Radon-Nikodym theorem applies whenever a measure is absolutely continuous with respect to its reference measure (so densities exist). In the rest of this document, “measure” means “$\sigma$-finite measure.”

A transition kernel (or kernel) from $X$ to $Y$ is a measurable function $\kappa: X \to M(Y)$, where $M(Y)$ is the space of measures on $Y$. When each $\kappa(x, \cdot)$ is a probability measure, the kernel is called a Markov kernel. In FlatPPL, kernels are represented as functions that map value points to measures.

The classical Giry monad (Giry, 1982) operates on probability measures, which are normalized. FlatPPL extends this to $\sigma$-finite measures in general, e.g. to represent non-normalized posteriors and intensity measures. Staton et al. (2016) and Staton (2017) provide the formal basis for this extension using the more general class of s-finite measures; all $\sigma$-finite measures are s-finite, so FlatPPL’s algebraic operations are well-founded within that framework.

Density convention. All density formulas in this section are with respect to a reference measure implied by the constituent distribution types: Lebesgue for continuous variates, counting measure for discrete variates. When a kernel $\kappa(\theta)$ is parameterized by $\theta$, the family is assumed dominated by a single $\theta$-independent reference measure.

Reference measure for product measures. When joint(M1, M2, ...) (or iid(M, size), jointchain(M, K1, ...) etc.) combines components with individual reference measures $\rho_1, \rho_2, \ldots$ (each either Lebesgue or Counting on the corresponding component support), the reference measure of the product is the product $\rho_1 \otimes \rho_2 \otimes \cdots$ on the joint variate space. For components sharing no stochastic ancestor, the joint density w.r.t. this product reference is the product of the component densities; a shared-ancestor joint keeps the same product reference, with density given by its equivalent record law (see joint). Mixed continuous-discrete joints are handled uniformly under this rule: e.g. joint(c = Normal(mu = 0, sigma = 1), k = Poisson(rate = 3)) has reference measure $\mathrm{Lebesgue}(\mathbb{R}) \otimes \mathrm{Counting}(\mathbb{Z})$ on $\mathbb{R} \times \mathbb{Z}$, with joint density $\phi(c;\,0,1) \cdot \mathrm{Pois}(k;\,3)$ at $(c, k)$.

Normalization convention. Normalization is always explicit in FlatPPL. Built-in distribution/measure constructors do not normalize their inputs, and measure-algebra operations never rescale their inputs or outputs.

6.2 The measure monad #

The Giry-style measure monad is defined by two operations:

• Unit: $\eta_X(x) = \delta_x$ (Dirac measure at $x$). In FlatPPL: Dirac(value = v).
• Bind: $(\nu \mathbin{\texttt{>>=}} \kappa)(B) = \int_X \kappa(x)(B)\, d\nu(x)$. In FlatPPL: kchain(M, K).

6.3 Fundamental measures and measure algebra #

6.3.1 Fundamental measures #

FlatPPL provides three fundamental measures: the reference measures Lebesgue and Counting, and the point-mass measure Dirac.

• Lebesgue(support = S) — the canonical continuous reference measure on the support set S, restricted to S. For full-dimensional subsets of Euclidean or product spaces this is the ordinary Lebesgue measure on the ambient space. For lower-dimensional embedded affine sets such as stdsimplex(n), it is the coordinate Lebesgue measure of the set’s free coordinates (see standard simplex), not its surface area.

  S may be any FlatPPL set: one-dimensional (e.g. reals, interval(0, 1), posreals), a Cartesian power (e.g. cartpow(reals, n)), a record-structured product (e.g. cartprod(a = reals, b = posreals)), a lower-dimensional embedded set (e.g. stdsimplex(n)) and so on (see sets).

  iid(Lebesgue(reals), n) is equivalent to Lebesgue(cartpow(reals, n)).

• Counting(support = S) — the counting measure on $\mathbb{Z}$, restricted to support S. Mass 1 at every integer in S. Reference measure for all discrete distributions.

• Dirac(value = v) — point-mass probability measure at v for any variate type.

The predefined constants reals (equivalent to interval(-inf, inf)) and integers (the set of all integers) serve as the default supports for the Lebesgue and counting measures respectively. The support parameter specifies where the measure is nonzero; density is zero outside. Measure algebra operations require their operands to share the same variate space (same type and dimension).

Uniform kernel extension. Mathematically, a measure is equivalent to a transition kernel with an empty first argument. So in FlatPPL, we unify measures and kernels and identify measures with nullary kernels. Measure algebra operations accept both kernels in general and measures as a (very important) special case of kernels. On a kernel, the operation applies to the output measure at each input point:

• pushfwd(f, K) denotes $\theta \mapsto \mathrm{pushfwd}(f, \kappa(\theta))$
• weighted(w, K) denotes $\theta \mapsto \mathrm{weighted}(w, \kappa(\theta))$

This applies to all measure-to-measure operations except jointchain and kchain, which require non-nullary kernels in all but the first argument (see dependent composition).

Operations that map measures to values, like totalmass, densityof, and logdensityof, require closed measures (i.e. nullary kernels) as inputs. densityof(M, x) and logdensityof(M, x) evaluate the density of a measure at a point with respect to an implicit reference measure.

To evaluate a density at many points (e.g. a grid for numerical integration or plotting), broadcast the operation rather than calling it per point: broadcast(fn(logdensityof(M, _)), grid) (equivalently fn(logdensityof(M, _)).(grid)) returns one log-density per grid element. The point argument stays scalar.

6.3.2 Density reweighting #

• weighted(weight, base) — produces the measure $\nu(A) = \int_A f(x)\, dM(x)$, with $d\nu = f \cdot dM$, where $f$ is a non-negative weight (a constant or a function of the variate $x$ of $M$) and $M$ the base measure. normalize(weighted(f, Lebesgue(support = S))) produces a probability distribution whose density w.r.t. Lebesgue on $S$ is proportional to $f$.

  Weight arity. A one-parameter weight receives the variate whole. If the variate is a $k$-element array with $k \geq 2$, a weight of exactly $k$ scalar parameters instead receives one component per parameter, in order; any other arity is an error. Only weighted and logweighted admit the $k$-parameter form; every other function-taking construct, including pushfwd, passes the whole variate.

  # equivalent over Lebesgue(support = cartprod(interval(0, 1), interval(0, 1)))
  w_components(x, y) = x * y
  w_variate(v) = v[1] * v[2]

• logweighted(logweight, base) — like weighted, but the weight or weighting function is given in log-space: $d\nu = \exp(g) \cdot dM$.

• bayesupdate(L, prior) — reweights a prior measure by a likelihood object, producing the unnormalized posterior: $d\nu(\theta) = L(\theta) \cdot d\pi(\theta)$. Lowers to logweighted(fn(logdensityof(L, _)), prior). See posterior construction for details.

6.3.3 Normalization and mass #

• normalize(M) — given a measure $M$ with finite total mass $Z = \mathrm{totalmass}(M) > 0$, returns the probability measure $M / Z$. If $Z = 0$ or $Z = \infty$, the result is undefined. On a non-nullary kernel, normalizes the output measures.

• totalmass(M) — returns the total mass $Z = \int dM(x)$ as a scalar value. Requires a closed measure (not a non-nullary kernel).

6.3.4 Additive superposition #

• superpose(M1, M2, ...) — measure addition: $\nu(A) = M_1(A) + M_2(A) + \ldots$ All components must share the same variate space. The result is generally not normalized. For example:

  intensity = superpose(weighted(amplitude, signal_shape), bkg_shape)
  events ~ PoissonProcess(intensity = intensity)

  To build a normalized mixture distribution, use normalize(superpose(weighted(w1, M1), weighted(w2, M2))). For example:

  mix = normalize(superpose(weighted(a1, normal1), weighted(a2, normal2)))

• ksuperpose(kernel, weights) — lifts a kernel to a weighted mixture: applied to a parameter family it yields $\nu = \sum_i w_i\,\kappa(\theta_i)$, one component per weight ($N$ of them, not necessarily statically known), with $\theta_i$ the parameters of component $i$. Each collection argument — passed as to broadcast, or as one table whose columns are the parameters — stacks one value per component along one extra leading axis: a scalar parameter takes a length-$N$ vector, a vector parameter an $N \times d$ matrix, a matrix parameter an $N \times d \times d$ array, and any other axis structure is a static error. A non-collection argument, or a collection whose leading axis has size one, is shared by every component, and the weights must be non-negative but need not sum to one. For example:

  mix = normalize(ksuperpose(Normal, weights)(mu = means, sigma = sigmas))

6.3.5 Joint composition #

• joint(M1, M2, ...) — the joint law of its components. A component contributes a fresh coordinate; a stochastic node shared between component traces (through a reified component — lawof, kernelof — or a stochastic constructor parameter) remains a single node of the composed trace. Components that share no stochastic node are independent, and their joint is the product measure: $(M_1 \otimes M_2)(A \times B) = M_1(A) \cdot M_2(B)$.

  The output variate is formed by combining the component variates via cat (see array operations). All components must have the same shape class: all scalars (yielding a vector), all vectors (yielding a concatenated vector), or all records with distinct field names (yielding a merged record). Mixing shape classes is a static error.

  For example, the measure product of a normal and an exponential probability measure, defined over a space of vectors, would be

  M1 = Normal(mu = 0, sigma = 1)
  M2 = Exponential(rate = 1.0)
  vj = joint(M1, M2)

  Keyword form. joint(name1 = M1, name2 = M2, ...) names the component variates, producing a measure over a space of records:

  rj = joint(name1 = M1, name2 = M2)

  is equivalent to joint(relabel(M1, ["name1"]), relabel(M2, ["name2"])).

  In this keyword form, a record-valued component becomes a nested record under its name — the name adds a level, it does not merge the inner fields (unlike the positional cat form above).

  For kernels, joint(K1, K2, ...) results in a kernel that fans a single input out to all component kernels, so each of them receives the same input. The result’s inputs are the union of the component kernels’ inputs by name; a component receives the inputs it declares and is unaffected by the others, as in reification with interdependent boundary nodes. Components that share a stochastic node must agree on that node’s ancestry: every ancestor of the shared node that any component binds as a boundary input must be bound by every sharing component, under the same input name. A joint in which a sharing component binds such an ancestor under a different name, or does not bind it at all — in particular a measure component, which binds nothing — is a static error. Measure components are permitted and are the nullary case: they ignore the input. A measure component may be parameterized and may share stochastic nodes with kernel components; only a shared node with a boundary-bound ancestor is excluded, by the naming clause above. The keyword form applies unchanged, producing a kernel whose output variate is a record. At each input point the result is the joint of the component output measures, governed by the sharing rules above; the fanned input is a value, not a stochastic node, and so induces no dependence by itself. The result’s total-mass class is the product of the components’ classes, as in the measure case; when components sharing a stochastic node include more than one non-normalized member, no class stronger than unknown is statically justified. A fan-out of Markov kernels is a Markov kernel.

  Equivalent record law. joint(a = lawof(a), b = lawof(b)) is equivalent to lawof(record(a = a, b = b)); the positional form is the corresponding cat law (see reification to measures).

  z ~ Normal(mu = m, sigma = s)
  a ~ Normal(mu = z, sigma = s_a)
  b ~ Normal(mu = z, sigma = s_b)

  For these draws, joint(a = lawof(a), b = lawof(b)) has cross-covariance $\mathrm{Var}(z) = s^2$; a joint of two constructor measures with the same marginals has cross-covariance $0$ only when their parameters reach no shared stochastic node.

  Singular joints. When one component’s variate is determined by the others given the shared ancestors (the same draw referenced twice, a deterministic transform of another component), the joint law has no density w.r.t. the product reference measure. Sampling is well-defined; a density query is a static error where statically detectable, and is otherwise refused by the engine.

• iid(M, size) — the product measure $M^{\otimes N}$ over arrays of shape size, where N = prod(size). size is a positive integer (1-D length) or a vector of positive integers (multi-axis shape). When M is a reified law, each of the $N$ copies carries its own copy of the reified sub-DAG, stochastic ancestors included; iid never shares nodes between copies. A size derived from data rather than written in source may resolve to 0, giving the empty product measure, whose log-density is $0$: the empty sum in the density rule for composed measures.

  When M is record-valued and size is a scalar length, the variate is an $N$-row table — one row per draw of M’s record — mirroring cartpow on a record set. A size that resolves to zero is then an error rather than the empty product measure above, as a table has no zero-row form. A multi-axis size retains the array shape.

  For example, to represent the draw of 100 IID samples from a normal distribution, use

  obs ~ iid(Normal(mu = a, sigma = b), 100)

6.3.6 Dependent composition #

• kchain(M, K1, K2, ...) — left-associative Kleisli composition (monadic bind). Keeps only the last kernel’s variates, marginalizing out all intermediate variates. In contrast to standard Kleisli composition, the first argument may also be a measure (a nullary kernel). See jointchain below for the variant that retains all variates.

  Mathematically, we define the chain of a measure $\mu(A)$ and a transition kernel $\kappa$ as

  $$\nu(B) = \int \kappa(a, B)\, d\mu(a)$$

  This involves a marginalization integral, which is generally intractable. Left-associative.

  prior_predictive = kchain(prior, forward_kernel)

  Equivalence with stochastic nodes:

  model = kchain(M1, K2, K3)

  is equivalent to

  a ~ M1
  b ~ K2(a)
  c ~ K3([a, b])
  model = lawof(c)

• jointchain(M, K1, K2, ...) — kernel-conditioned joint measure. The first argument is a base measure or kernel; the remaining arguments are non-nullary kernels whose inputs bind to the variates of everything to their left.

  jointchain is left-associative. In contrast to kchain, the output variate is the cat of the variates of all the components, as with joint.

  Keyword form. jointchain(name1 = M, name2 = K1, ...) names the component variates, producing a measure over a space of records. It is equivalent to jointchain(relabel(M, ["name1"]), relabel(K1, ["name2"]), ...).

  Mathematically, we define the joint chain of a measure $\mu(A)$ and a transition kernel $\kappa$ as

  $$\nu(A \times B) = \int_A \kappa(a, B)\, d\mu(a)$$

  The density of the joint chain is the product of the constituent conditional densities — no marginalization integral is involved, unlike with kchain. So density is tractable if the densities of all the components are.

  Equivalence with stochastic nodes:

  model = jointchain(M1, K2, K3)

  is equivalent to

  a ~ M1
  b ~ K2(a)
  c ~ K3([a, b])
  model = lawof([a, b, c])

  Relationship to kchain:

  jointchain(M, K)

  is equivalent to

  kchain(M, a -> joint(Dirac(value = a), K(a)))

  Like fchain, kchain and jointchain combine well with auto-splatting: a record-shaped variate from step $i$ splats into step $i+1$’s keyword inputs by field name. A non-record variate — for example the cat’d variate of a positional joint — carries no field names, so it feeds a kernel only when the kernel has a single input, to which the whole value is bound; feeding one to a kernel with two or more inputs is a static error, as a single value cannot be split across inputs by name. Use the named form (joint(name1 = M1, ...)) or relabel to name the components, producing a record variate whose fields splat by name.

• markovchain(kernel, init, n) — measure over length-n trajectories of a time-homogeneous Markov chain.

  kernel is a Markov kernel (state) -> measure_over_state; init is a value in the state space; n is a positive integer. Step $i$ is $\text{traj}_i \sim \kappa(\text{traj}_{i-1})$ with $\text{traj}_0 = \text{init}$. The initial value is not part of the trajectory. The resulting measure is a measure over arrays [traj[1], ..., traj[n]], excluding the initial state. If init and traj[i] are records, then the trajectories are tables, not arrays.

  Example — Brownian motion (100 steps of $\Delta t = 0.01$, starting at zero):

  D = 4.1    % Diffusion constant
  dt = 0.01  % Time step
  f_step = x -> Normal(x, sqrt(2*D * dt))
  traj ~ markovchain(f_step, 0.0, 100)

• kscan(kernel, init, xs) — Kleisli scan that generalizes markovchain with exogenous inputs threaded through each step, also a stochastic version of scan.

  kernel is a Markov kernel (state, x) -> measure_over_state; step $i$ is $\text{traj}_i \sim \kappa(\text{traj}_{i-1}, \text{xs}_i)$ with $\text{traj}_0 = \text{init}$. Trajectories have length lengthof(xs). As with markovchain, init is a value in the state space and not part of the trajectory.

  Example — Brownian motion with variable timesteps:

  D = 4.1  % Diffusion constant
  dts = [0.01, 0.02, 0.015, 0.018, 0.012]  % Time steps
  f_step = (x, dt) -> Normal(x, sqrt(2*D * dt))
  traj ~ kscan(f_step, 0.0, dts)

6.3.7 Support restriction #

• truncate(M, S) — restricts the support of measure M to the set S: $\nu(A) = M(A \cap S)$. Does not normalize automatically.

  half_normal = normalize(truncate(Normal(mu = 0, sigma = 1), interval(0, inf)))

6.3.8 Transformation and projection #

• pushfwd(f, M) — pushforward of measure $M$ through function $f$:

  $$(f_* M)(Y) = M(f^{-1}(Y))$$

  For kernels, pushfwd acts on their output measures.

  For example, a log-normal probability measure can be constructed as

  mu = Normal(mu = 0, sigma = 1)
  nu = pushfwd(exp, mu)  # → LogNormal

  The equivalent in stochastic-node form is:

  mu = Normal(mu = 0, sigma = 1)
  x ~ mu
  y = exp(x)
  nu = lawof(y)

  A pushforward can also be used to project, respectively marginalize:

  mu = relabel(iid(Normal(mu = 0, sigma = 1), 3), ["a", "b", "c"])
  pushfwd(fn(get(_, ["a", "c"])), mu)   # marginalizes out b

• locscale(m, shift, scale) — affine (location-scale) pushforward, shorthand for pushfwd(x -> scale * x + shift, m). For example, locscale(Normal(0, 1), mu, sigma) is equivalent to Normal(mu, sigma), and locscale(StudentT(nu), mu, sigma) is the location-scale Student-t. When m is vector-valued, scale may be a matrix: locscale(MvNormal(zeros(n), eye(n)), mu, lower_cholesky(cov)) is equivalent to MvNormal(mu, cov) — the affine map x -> lower_cholesky(cov) * x + mu. shift and scale must be value-compatible with the variate of m; for general matrix-vector affine maps use pushfwd directly.

• bijection(f, f_inv, logvolume) annotates a function f with its inverse f_inv and the log-volume-element logvolume of the forward map. The result is a function that is semantically f.

  FlatPPL engines will often need the inverse of f and the volume element when computing densities of pushforward measures. Function inverses are hard to derive automatically and the computation of Jacobian determinant via automatic differentiation can be very inefficient, while the user or system that authors/generates FlatPPL may have access to both in closed form.

  logvolume is the generalized log-volume-element of the forward function — it generalizes the log-absolute-determinant of the Jacobian to mappings between spaces of different dimension. It may be a function or a scalar value (logvolume = 0 for volume-preserving bijections). The convention is that logvolume describes the forward map.

  The user asserts that f_inv is the inverse of f and that logvolume is correct with respect to how f is used in the FlatPPL module. FlatPPL implementations are not required to verify this.

  For standard cases like exp, FlatPPL engines can be expected to know the inverse and volume element, but it would be written in FlatPPL as

  exp_bijection = bijection(exp, log, identity)

  A more interesting example that includes an explicit definition of domain and codomain of the function is squaring on the positive reals:

  pos_x = elementof(interval(0, inf))
  sq = bijection(
      functionof(pos_x ^ 2, x = pos_x),
      functionof(sqrt(pos_x), x = pos_x),
      fn(log(2 * _))
  )

6.3.9 Engine contract for pushfwd density evaluation #

densityof(pushfwd(f, M), y) and logdensityof(pushfwd(f, M), y) require the engine to invert f and apply the volume element. For a bijection f with inverse f_inv and forward log-volume logvolume, the density is given by the change-of-variables formula

$$\log \mathrm{densityof}(\mathrm{pushfwd}(f, M), y) = \log \mathrm{densityof}(M, f^{-1}(y)) - \mathrm{logvolume}(f^{-1}(y))$$

equivalently $\mathrm{densityof}(\mathrm{pushfwd}(f, M), y) = \mathrm{densityof}(M, f^{-1}(y)) \cdot \exp\!\left(-\mathrm{logvolume}(f^{-1}(y))\right)$. The forward log-volume is evaluated at the preimage $f^{-1}(y)$ and subtracted (e.g. for exp_bijection, logvolume = identity, giving the log-normal density $\log \mathrm{densityof}(M, \log y) - \log y$). Engines must support density evaluation in the following three cases: