Guillaume Baudart

Louis Mandel

Eric Atkinson

Benjamin Sherman

Marc Pouzet

Michael Carbin

Reactive Probabilistic Programming

Keywords: CCS Concepts:, Theory of computation → Streaming models;, Software and its engineering → Data flow languages Probabilistic programming, Reactive programming, Streaming inference, Semantics, Compilation

Synchronous modeling is at the heart of programming languages like Lustre, Esterel, or SCADE used routinely for implementing safety critical control software, e.g., fly-bywire and engine control in planes. However, to date these languages have had limited modern support for modeling uncertainty -probabilistic aspects of the software's environment or behavior -even though modeling uncertainty is a primary activity when designing a control system.

In this paper we present ProbZelus the first synchronous probabilistic programming language. ProbZelus conservatively provides the facilities of a synchronous language to write control software, with probabilistic constructs to model uncertainties and perform inference-in-the-loop.

We present the design and implementation of the language. We propose a measure-theoretic semantics of probabilistic stream functions and a simple type discipline to separate deterministic and probabilistic expressions. We demonstrate a semantics-preserving compilation into a first-order functional language that lends itself to a simple presentation of inference algorithms for streaming models. We also redesign the delayed sampling inference algorithm to provide efficient streaming inference. Together with an evaluation on several reactive applications, our results demonstrate that ProbZelus enables the design of reactive probabilistic applications and efficient, bounded memory inference.

Introduction

Synchronous languages [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] were introduced thirty years ago for designing and implementing real-time control software. They are founded on the synchronous abstraction [START_REF] Berry | Real Time Programming: Special Purpose or General Purpose Languages[END_REF] where a system is modeled ideally, as if communications and computations were instantaneous and paced on a global clock. This abstraction is simple but powerful: input, output and local signals are streams that advance synchronously and a system is a stream function. It is at the heart of the data-flow languages Lustre [START_REF] Halbwachs | The Synchronous Dataflow Programming Language Lustre[END_REF] and SCADE [START_REF] Colaço | SCADE 6: A formal language for embedded critical software development (invited paper)[END_REF]; it is also the underlying model behind the discrete-time subset of Simulink.

The data-flow programming style is very well adapted to the direct expression of the classic control blocks of control engineering (e.g., relays, filters, PID controllers, control logic), and a discrete time model of the environment, with the feedback between the two. For example, consider a backward Euler integration method defined by the following stream equations and its corresponding implementation in Zelus [START_REF] Bourke | Zélus: a synchronous language with ODEs[END_REF], a language reminiscent of Lustre:

x 0 = xo 0 x n = x n-1 + x ′ n × h ∀n ∈ N, n > 0
let node integr (xo, x') = x where rec x = xo -> (pre x + x' * h)

The node integr is a function from input streams xo and x' to output stream x. The initialization operator -> returns its left-hand side value at the first time step and its right-hand side expression on every time step thereafter. The unit-delay operator pre returns the value of its expression at the previous time step. The following table presents a sample timeline showing the sequences of values taken by the streams defined in the program (where h is set to 0.1).

xo 0 0 0 0 0 0 0 . . . The node integr can be used to define other stream functions, e.g., a PID controller, which can be called in control structures like hierarchical automata, e.g., to express a system that switches between automatic and manual control. Compared to a general purpose functional language (or an embedded DSL), the expressiveness of a synchronous language is purposely constrained to modularily ensure safety properties that are critical for the targeted applications: determinacy, deadlock freedom (reactivity), the generation of statically scheduled code that runs in bounded time and space.

However, to date, these languages have had limited support for modeling uncertainty (e.g., a noisy sensor or channel, a variable delay), to simulate the interaction of a software controller and a partially unknown environment, or to infer parameters from noisy observations. Moreover, uncertainty is a first-order design concern for a controller that operates under the assumption of a probabilistic model of their environment (e.g., object tracking). Using a probabilistic environment model and data gathered from observing the environment, a controller can infer a distribution of likely environments given the observations. Existing approaches consist in hand-coding stochastic controllers that have a known solution (e.g., Kalman filters) which can be tedious and error-prone, or to simply perform off-line statistical testing on the generated code of a controller. Alternatively, in recent years, probabilistic programming has developed as an approach to endow general purpose languages with the ability to automate inference. Probabilistic Programming. Probabilistic programming languages are used to describe probabilistic models and automatically infer distributions of latent (i.e., unobserved) parameters from observations.

A popular approach [START_REF] Bingham | Pyro: Deep Universal Probabilistic Programming[END_REF][START_REF] Goodman | The Design and Implementation of Probabilistic Programming Languages[END_REF][START_REF] Murray | Automated learning with a probabilistic programming language: Birch[END_REF][START_REF] Tolpin | Design and Implementation of Probabilistic Programming Language Anglican[END_REF][START_REF] Tran | Deep Probabilistic Programming[END_REF][START_REF] Wu | Swift: Compiled Inference for Probabilistic Programming Languages[END_REF] consists in extending a general-purpose programming language with three constructs: [START_REF] Baudart | Reactive Probabilistic Programming[END_REF] x = sample(d) introduces a random variable x of distribution d, [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] observe(d, y) conditions on the fact that a sample from distribution d produced the observation y, and (3) infer m obs computes the distribution of the output values of a program or model m w.r.t. the observation of the input data obs.

Probabilistic programming languages offer a variety of automatic inference techniques ranging from exact inference [START_REF] Gehr | PSI: Exact Inference for Probabilistic Programs[END_REF] to approximate inference [START_REF] Narayanan | Probabilistic Inference by Program Transformation in Hakaru (System Description)[END_REF] and include hybrid approaches that combine exact and approximate techniques when part of the program is analytically tractable [START_REF] Murray | Delayed Sampling and Automatic Rao-Blackwellization of Probabilistic Programs[END_REF]. However, a standing challenge for these programming languages is that none of them meet the design goals of synchronous reactive languages by being immediately amenable to techniques to ensure that for example a program with an indefinite execution time runs in bounded memory.

Inference in the Loop. In this paper we extend Zelus1 to provide a synchronous probabilistic programming language, ProbZelus. ProbZelus enables one to combine deterministic reactive data-flow programs, such as integr (above), with probabilistic programming constructs to produce reactive probabilistic programs.

Compared to other probabilistic languages (e.g. WebPPL, Church, Stan) where inference is executed on terminating pure functions, our probabilistic models are stateful stream processors. Inference on probabilistic models runs in parallel with deterministic processes that interact with the environment. The distributions computed by infer at each step can thus be used by deterministic processes to compute new inputs for the next inference step. We term this capability inference-in-the-loop.

Streaming Inference. ProbZelus provides multiple inference algorithms, most notably the delayed sampling inference algorithm [START_REF] Murray | Delayed Sampling and Automatic Rao-Blackwellization of Probabilistic Programs[END_REF]. This hybrid strategy combines the approximate inference technique of particle filtering [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF] with exact inference when it is possible to symbolically determine the exact distribution for some or all of the latent variables of the program [START_REF] Doucet | Rao-Blackwellised Particle Filtering for Dynamic Bayesian Networks[END_REF].

However, the memory consumption of delayed sampling strictly increases with the number of random variables which is not practical for reactive applications that operate on infinite streams. We propose a novel streaming implementation of delayed sampling that can operate over infinite streams in constant memory for a large class of models. ProbZelus therefore provides an expressive language for reactive probabilistic programming with appropriate memory consumption properties.

Contributions. We present the following contributions: Design, Semantics, Compilation. We present ProbZelus, the first synchronous probabilistic programming language, combining language constructs for streams (reactivity) with those for probabilistic programming thus enabling inferencein-the-loop. We give a measure-based co-iterative semantics for ProbZelus that forms the basis of a compiler and demonstrate a semantics-preserving compilation strategy to a first-order functional language µF . Streaming Inference. We define the semantics of multiple inference algorithms on µF including particle filtering and delayed sampling. We then present a novel streaming delayed sampling implementation which enables partial exact inference over infinite streams in bounded memory for a large class of models.

let proba kalman (xo, u, acc, gps) = x where rec mu = xo -> (a *@ pre x) +@ (b *@ u) and x = sample (mv_gaussian (mu, noise)) and () = observe (gaussian (vec_get (x, 2), 1.0), acc) and () = present gps(pos) -> observe (gaussian (vec_get (x, 0), 0.01), pos) else () let node robot (xo, uo, acc, gps) = u where rec x_dist = infer 1000 kalman (xo, u, acc, gps) and u = uo -> lqr a b (mean (pre x_dist))

Figure 1. Robot controller with inference-in-the-loop. +@ and *@ are matrix operations, vec_get x i is the ith projection.

Evaluation. We evaluate the performance of ProbZelus on a set of benchmarks that illustrate multiple aspects of the language. We demonstrate that streaming delayed sampling drastically reduces the number of particles required to achieve better accuracy compared to a particle filter.

The result is ProbZelus, a synchronous probabilistic language that enables us to write, in the very same source, a deterministic model for the control software and a probabilistic model with complex interactions between the two. On one hand, a deterministic model of a controller can rely on predictions computed by a probabilistic model. On the other hand, a probabilistic model can be programmed in an expressive reactive language. ProbZelus is open source (https://github.com/IBM/probzelus). An extended version with appendices of the paper is also available [START_REF] Baudart | Reactive Probabilistic Programming[END_REF].

Example

In this section, we demonstrate how ProbZelus provides probabilistic modeling, inference-in-the-loop, and boundedmemory inference for a robot navigation system. The results of the inference are continuously used by a controller to correct the robot trajectory.

Inference in the Loop.

We consider a robot equipped with an accelerometer and a GPS. We assume that the motion of the robot can be described as: x t +1 = Ax t + Bu t where x t denotes the state of the robot (position, velocity, and acceleration) at a given time step t, and u t denotes the command sent to the robot. A and B are constant matrices. In addition, the robot receives at each step noisy observations from the accelerometer a t , and sporadically an estimation of the position from a GPS p t .

Figure 1 presents a controller, robot, that given an initial state xo, an initial command uo, and inputs from the accelerometer acc and the GPS gps computes a stream u of commands that drives the robot to a given target. The body of robot is the parallel composition of (1) the inference of a probabilistic process kalman that estimates x_dist the stream of distributions over the robot's state, and (2) a deterministic The command u is set to the initial command uo at the first time step, and is then computed by a Linear-Quadratic Regulator (LQR) [START_REF] Sontag | Mathematical control theory: deterministic finite dimensional systems[END_REF] -a stable and optimal controller for such dynamic systems -given the estimation of the state at the previous step. Because LQR controllers depend only on mean posterior state, the example in Figure 1 uses the mean function to compute the mean of x_dist before invoking the LQR controller.

Inference. The stream x_dist of distributions of state is inferred from the model defined by the probabilistic node kalman given the initial state xo, the command u, and the observations acc and gps. The keyword proba indicates a probabilistic model.

In this example, the model is a Kalman Filter illustrated in Figure 2. A Kalman filter is a time-dependent probabilistic model used to describe inference problems such as tracking, in which a tracker estimates the true position of an object given noisy, sensed observations. The robot's state x t is a latent random variable in that the tracker is not able to directly observe it. Each arrow connecting two random variables denotes a dependence of the variable at the head of the arrow on the variable at the tail. In this case, the observations at each time step depends on the current state, and the robot's state at a given time step depends only on its states at the previous time step.

Sampling. Inside the kalman node, the sample operator samples a value from a probability distribution. In this case, the expression samples the current state x from a multivariate Gaussian with mean obtained by applying the motion model to the previous state and the command. This code models the trajectory of a robot where at each time step, the state is Gaussian-distributed around an estimation computed from the motion model.

Observations. The expression observe conditions the execution on observed data. Its first parameter denotes a distribution that models the observation and its second parameter denotes the observed value itself. In this case, at each step, the first observe statement models a Gaussian-distributed observation of the current acceleration vec_get x 2 given by acc. The input gps is a signal that is only emitted when the GPS computes a new position. When a value pos is emitted on gps, the present construct executes its left branch, further conditioning the model by adding a Gaussian-distributed observation of the current position vec_get x 0 given by pos.

Streaming Inference

A classic operational interpretation of a probabilistic model is an importance sampler that generates random samples from the model together with an importance weight measuring the quality of the sample. In this model, each execution of a sample operator samples a value from the operator's corresponding distribution. Each execution of an observe evaluates the likelihood of the provided observation and multiplies the current importance weight by this value. Then, each execution step of infer yields a distribution represented as a set of pairs (output, weight) or particles. The particles can be re-sampled at each step to build a particle filter [START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF].

The integer parameter to infer determines how many particles to use: the more particles the user specifies, the more accurate the estimate of the distribution becomes. The PF points in Figure 3 present this improvement in accuracy as a function of increasing the number of particles for the robot example. However, as the latency results presents, the more particles the user specifies, the more computation is required for each step because each particle requires a full, independent execution of each time step of the model.

Streaming Delayed Sampling. Delayed sampling [START_REF] Murray | Delayed Sampling and Automatic Rao-Blackwellization of Probabilistic Programs[END_REF] can reduce the number of particles required to achieve a given desired quality of inference. Specifically, delayed sampling exploits the opportunity to symbolically reason about the relationships between random variables to compute closedform distributions whenever possible. To capture relationships between random variables, delayed sampling maintains a graph: a Bayesian network that can be used to compute closed-form distributions involving subsets of random variables. For instance, this inference scheme is able to compute the exact posterior distribution for our robot example. The SDS dots in Figure 3 show that the accuracy is independent of the number of particles since each particle computes the exact solution.

Figure 4 illustrates the evolution of the delayed sampling graph as it proceeds through the first four time steps of the robot example (for simplicity we assume that there is no GPS activation in these four steps). A notable challenge with the traditional delayed sampling algorithm is that the graph grows linearly in the number of samples. This property is not tractable in our reactive context because we would like to deploy our programs under the model that they run indefinitely, thus requiring that they execute with bounded resources. To address this problem, we propose a novel streaming delayed sampling (SDS) implementation of the delayed sampling algorithm. Specifically, in Figure 4 the node denoting the marginal posterior for x at step 1 can be eliminated from the graph at step 3 because the distributions for pre x and x have fully incorporated its effect on their values and, moreover, the program no longer maintains a reference to the node.

While the standard delayed sampling algorithm will keep this node alive through the edge pointers it maintains, SDS builds a pointer-minimal graph representation with a minimal number of edges that 1) ensure that the graph has sufficient connectivity to support operations in the traditional delayed sampling algorithm and 2) only maintain the reachability of nodes that can effect the distribution of future nodes in the graph. The result is that the memory consumption of SDS is constant across the number of steps while the memory consumption of the original delayed sampling implementation DS increases linearly in the number of steps (Figure 5).

Language: Syntax, Typing, Semantics

ProbZelus is a reactive probabilistic language with inferencein-the-loop which enables interaction between probabilistic models and deterministic processes. This capability introduces two design requirements. First, a probabilistic model must be able to receive inputs from an evolving environment. Second, instead of awaiting the final result of the inference, deterministic processes running in parallel need access to intermediate results. The resulting inference-in-the-loop enables feedback loops between inferred distributions from probabilistic models and deterministic processes, which our design controls by enforcing a separation between the semantics of probabilistic and deterministic execution.

In this section, we formalize the syntax of ProbZelus, introduce a type system that imposes a clear separation between deterministic and probabilistic expressions, and define the semantics of the language in a co-iteration framework where the semantics of probabilistic processes is adapted from Staton's measure-theoretic semantics for probabilistic programs [START_REF] Staton | Commutative Semantics for Probabilistic Programming[END_REF]. The co-iterative semantics forms the basis of a compiler that is described in Section 4.

Syntax

We focus on the following kernel of ProbZelus. The missing constructs (e.g., pre and ->) can be compiled into this kernel via a source-to-source transformation. A program is a sequence of declarations d of stream functions (node) and probabilistic stream functions (proba). An expression e is either a constant (c), a variable (x), a pair, an external operator application (op), a function application (f (e)), a delay (last x) that returns a value (x) from the previous step, or a set of locally recursive equations (e where rec E).

A set of equations E is either an equation x = e that define x with the stream e, the initialization of a variable with a constant init x = c, or parallel composition of sets of equations.

Operators (op) include boolean and arithmetic operators. In addition, ProbZelus offers a library of dedicated operators to analyze distributions, such as mean and variance, that can be used in any context (probabilistic or deterministic), e.g., on the result of the inference.

The control structure present e -> e 1 else e 2 is an activation condition that executes the expression e 1 only when the value of e is true and executes e 2 otherwise. It differs from if e then e 1 else e 2 , where both e 1 and e 2 are computed at each step (making their internal states evolve) and the returned value is chosen based on the value of e. 2 In the following example, o1 and o2 are different streams: ...

The reset e 1 every e construct re-initializes the values of the init equations and the corresponding last expressions in e 1 each time e is true.

The language is extended with the classic probabilistic expressions: sample to draw from a distribution, observe to condition on observations, and infer to compute the distribution described by a model.

Scheduling. In the expression e where rec E, E is a set of mutually recursive equations. In practice, a scheduler reorders the equations according to their dependencies. Initializations init x j = c j are grouped at the beginning, and an equation x j = e j must be scheduled after the equation x i = e i if the expression e j uses x i outside a last construct. A program satisfying this partial order is said to be scheduled. The compiler can also introduce additional equations to relax the scheduling constraints and rejects programs that cannot be statically scheduled [START_REF] Biernacki | Clock-directed modular code generation for synchronous data-flow languages[END_REF]. After scheduling, the expression e where rec E has the following form. For simplicity, we also assume that every initialized variable is defined in a subsequent equation, i.e., {x i } 1..k ∩ {y j } 1..n = {x i } 1..k . If it is not the case, in this kernel we can always add additional equations of the form x i = last x i .

Kernel. All ProbZelus programs can be encoded in this kernel language. For instance, the program integr of Section 1 can be rewritten in the kernel as follows:

let node integr (xo, x') = x where rec init first = true and init x = 0. and first = false and x = if last first then xo else last x + (x' * h)

A stream first is defined such that last first is only true at the first step. The -> operator is then compiled into an if statement. The pre operator is compiled into a last operator. The initialization value is arbitrary, the compiler's initialization analysis guaranties that this value is never used [START_REF] Colaço | Type-based initialization analysis of a synchronous dataflow language[END_REF].

Similarly, other constructs like hierarchical automata can be re-written using present and reset [START_REF] Colaço | Mixing signals and modes in synchronous data-flow systems[END_REF].

3.2 Typing: Deterministic vs. Probabilistic Deterministic and probabilistic expressions have distinct interpretations. A dedicated type system discriminates between the two kinds of expressions, assigning one of two kinds to each expression: D for deterministic, or P for probabilistic. The typing judgment G ⊢ k e : T states that in the environment G which maps variable names to their type, the expression e has kind k and type T . Function types T → k T ′ are extended with the kind k of the body and we introduce a new datatype T dist for the probability distribution over values of type T . The expressions sample, and observe are probabilistic but their arguments must be deterministic. Any deterministic expression can be lifted to a probabilistic expression using a sub-typing rule. The transition from probabilistic to deterministic is realized via infer (the complete set of typing rules is presented in Figure 12

i γ = () ⟦x⟧ s γ = λs. (γ (x), s) ⟦present e -> e 1 else e 2 ⟧ i γ = (⟦e⟧ i γ , ⟦e 1 ⟧ i γ , ⟦e 2 ⟧ i γ ) ⟦present e -> e 1 else e 2 ⟧ s γ = λ(s, s 1 , s 2 ). let v, s ′ = ⟦e⟧ s γ (s) in if v then let v 1 , s ′ 1 = ⟦e 1 ⟧ s γ (s 1 ) in (v 1 , (s ′ , s ′ 1 , s 2 )) else let v 2 , s ′ 2 = ⟦e 2 ⟧ s γ (s 2 ) in (v 2 , (s ′ , s 1 , s ′ 2
. , s n ), s). let γ 1 = γ [m 1 /x 1 _last] in . . . let γ k = γ k -1 [m k /x k _last] in let v 1 , s ′ 1 = ⟦e 1 ⟧ s γ k (s 1 ) in let γ ′ 1 = γ k [v 1 /y 1 ] in . . . let v n , s ′ n = ⟦e n ⟧ s γ ′ n-1 (s n ) in let γ ′ n = γ ′ n-1 [v n /y n ] in let v, s ′ = ⟦e⟧ s γ ′ n (s) in v, ((γ ′ n [x 1 ], . . . ,γ ′ n [x k ]), (s ′ 1 , . . . , s ′ n ), s ′ ) Figure 6.
Semantics of deterministic expressions.

The type T dist * is a subtype of T dist that represents distributions known to have a density, i.e., discrete distributions (w.r.t. the counting measure) and a subset of continuous distributions (w.r.t. the Lebesgue measure). In ProbZelus, to simplify the semantics of the language, the observe statement requires a value of type T dist * . In Appendix A, we extend the language with the factor statement for arbitrary conditioning.

Co-Iterative Semantics

We now give the semantics of ProbZelus in a co-iteration framework [START_REF] Caspi | A Co-iterative Characterization of Synchronous Stream Functions[END_REF]. In this framework, a deterministic stream of type T is defined by an initial state of type S and a transition function of type S → T × S.

CoStream(T , S) = S × (S → T × S)
Repeatedly executing the transition function from the initial state yields a stream of values of type T . The semantics of a deterministic expression (kindOf (e) = D) is defined using two auxiliary functions. If γ is an environment mapping variable names to values, ⟦e⟧ i γ denotes the initial state, and ⟦e⟧ s γ denotes the transition function:

⟦e⟧ γ : CoStream(T , S) = ⟦e⟧ i γ , ⟦e⟧ s γ
The deterministic semantics of ProbZelus presented in Figure 6 is an extension of [START_REF] Caspi | A Co-iterative Characterization of Synchronous Stream Functions[END_REF] with the control structures present and reset (see also Figure 13 of Appendix A.2).

The transition function of a variable always returns the corresponding value stored in the environment γ .

The present e -> e 1 else e 2 construct introduced in Section 2 returns the value of e 1 when e is true and the value of e 2 otherwise. The state (s, s 1 , s 2 ) stores the state of the three sub-expressions. The transition function lazily executes the expression e 1 or e 2 depending on the value of e and returns the updated state.

The state of a set of scheduled locally recursive definitions e where rec E comprises three parts: the value of the local variables at the previous step which can be accessed via the last operator, the state of the defining expressions, and the state of expression e. The initialization stores the initial values introduced by init and the initial states of all subexpressions. The transition function incrementally builds the local environment defined by E. First the environment is populated with a set of fresh variables x i _last initialized with the values stored in the state that can then be accessed via the last operator. Then the environment is extended with the definition of all the variables y i by executing all the defining expressions (where

{x i } 1..k ∩ {y j } 1..n = {x i } 1..k ).
Finally, the expression e is executed in the final environment. The updated state contains the value of the initialized variables defined in E that will the be used to start the next step, and the updated state of the sub-expressions.

Probabilistic extension. The semantics of a probabilistic expression (kindOf (e) = P) follows the same scheme, but the transition function returns a measure over the set of possible pairs (result, state):

CoPStream(T , S) = S × (S → (Σ T ×S → [0, ∞]))
A measure µ associates a positive number to each measurable set U ∈ Σ T ×S where Σ T ×S denotes the Σ-algebra of T × S, i.e., the set of measurable sets over pairs (result, state). We use the following notation for the semantics of a probabilistic expression e:

{[e]} γ : CoPStream(T , S) = {[e]} i γ , {[e]} s γ
The semantics of probabilistic expressions is presented in Figure 7 (the complete semantics is in Figure 14 of Appendix A.2). This measure-based semantics is adapted from [START_REF] Staton | Commutative Semantics for Probabilistic Programming[END_REF] to explicitly handle the state of the transition functions.

First, any deterministic expression can be lifted as a probabilistic expression. The transition function returns the Dirac delta measure (δ x (U ) = 1 if x ∈ U , 0 otherwise) on the pair returned by the deterministic transition function applied on the current state: ⟦e⟧ s γ (s) : T × S. The probabilistic operator sample(e) evaluates e which returns a distribution µ : T dist and a new state s ′ : S, and returns a measure over the pair (result, state) where the state is fixed to the value s ′ . observe(e 1 , e 2 ) evaluates e 1 and e 2 into a distribution with density µ : T dist * and a value v : T , and weights execution paths using the likelihood of v w.r.t. µ (µ pdf denotes the density function of the distribution µ).

{[e]} i γ = ⟦e⟧ i γ if kindOf (e) = D {[e]} s γ = λs. λU . δ ⟦e⟧ s γ (s) (U ) if kindOf (e) = D {[sample(e)]} i γ = ⟦e⟧ i γ {[sample(e)]} s γ = λs. λU . let µ, s ′ = ⟦e⟧ s γ (s) in ∫ T µ(dv) δ v ,s ′ (U ) {[observe(e 1 , e 2 )]} i γ = (⟦e 1 ⟧ i γ ,⟦e 2 ⟧ i γ ) {[observe(e 1 , e 2 )]} s γ = λ(s 1 , s 2 ). λU . let µ, s ′ 1 = ⟦e 1 ⟧ s γ (s 1 ) in let v, s ′ 2 = ⟦e 2 ⟧ s γ (s 2 ) in µ pdf (v) * δ (),(s ′ 1 ,s ′ 2 ) (U )
e where rec init x 1 = c 1 ... and init x k = c k and y 1 = e 1 ... and

y n = e n s γ = λ((m 1 , . . . ,m k ), (s 1 , . . . , s n ), s). λU . let γ 1 = γ [m 1 /x 1 _last] in . . . let γ k = γ k -1 [m k /x k _last] in let µ 1 = {[e 1 ]} s γ k (s 1 ) in ∫ µ 1 (dv 1 , ds ′ 1 )let γ ′ 1 = γ k [v 1 /y 1 ] in . . . let µ n = ⟦e n ⟧ s γ ′ n-1 (s n ) in ∫ µ n (dv n , ds ′ n )let γ ′ n = γ ′ n-1 [v n /y n ] in let µ = {[e]} s γ ′ n (s) in ∫ µ(dv, ds ′ ) δ v ,((γ ′ n [x 1 ],...,γ ′ n [x k ]),(s ′ 1 ,...,s ′ n ),s ′ ) (U ) Figure 7. Semantics of probabilistic expressions.
The state of a set of locally recursive definitions is the same as in Figure 6 and contains the previous value of the initialized variables and the states of the sub-expressions. The transition starts by adding the variables x i _last to the environment. We note ∫ µ(dv, ds)f (v, s) the integral of f w.r.t. the measure µ where variables v and s are the integration variables. The integration measure appears on the right of the integral to maintain the expression order of the source code and we allow local definitions (e.g., let x = v in . . . ) inside the integral to simplify the presentation. Local definitions are interpreted by successively integrating the measure on pairs (value, state) returned by the defining expressions. In other words, we integrate over all possible executions. Integrals need to be nested to capture the eventual dependencies in the successive expressions. The returned value is a measure on pairs (value, state) where the state captures the value of the initialized variables and the state of the sub-expressions.

Inference in the loop. The infer operator is the boundary between the deterministic and the probabilistic expressions. Given a probabilistic model defined by an expression e, at each step the inference computes a distribution of results and a distribution of possible next states. Expression e can contain free variables thus capturing inputs from deterministic processes. The distribution of results can be used by deterministic processes to produce new inputs for the next inference step.

⟦infer(e)⟧ i γ = λU . δ ⟦e⟧ i γ (U ) ⟦infer(e)⟧ s γ = λσ . let µ = λU . ∫ S σ (ds){[e]} s γ (s)(U ) in let ν = λU . µ(U )/µ(⊤) in (π 1 * (ν), π 2 * (ν))
The state of infer(e) is a distribution over the possible states for e. The initial state is the Dirac delta measure on the initial state of e. The transition function integrates the measure defined by e over all the possible states and normalize the result µ to produce a distribution ν : T × S dist (⊤ denotes the entire space). This distribution is then split into a pair of marginal distributions using the pushforward of µ across the projections π 1 and π 2 .

Remark. In ProbZelus deterministic processes that interact with the environment cannot rollback on actions based on past estimations (e.g., the command of the robot controller). However, the inferred distribution of a random variable captured in the state may evolve at each time step based on subsequent observations (e.g., the initial position of the robot). These two properties follow from the separation of the distribution of results and the distribution of states in the semantics of infer.

Alternatively, we could define a fix-point semantics of streams based on a Scott order, as simple lazy streams in Haskell, where the value of a stream can depend on future computations (as illustrated Appendix A.3). However, this approach is not practical in a reactive context where processes interact with the environment [START_REF] Caspi | Clocks in Dataflow Languages[END_REF][START_REF] Kahn | The Semantics of a Simple Language for Parallel Programming[END_REF].

Compilation

Following the semantics described in Section 3.3 each expression is compiled into a transition function that can be written in a simple functional first-order language extended with probabilistic operators we call µF . Importantly, the compilation process is the same for deterministic and probabilistic expressions. We can then give a classic interpretation to deterministic terms, and a measure-based semantics to probabilistic terms following [START_REF] Staton | Commutative Semantics for Probabilistic Programming[END_REF]. We then show that the semantics of the compiled code coincides with the co-iterative semantics described in Section 3.3.

A First-Order Functional Probabilistic Language

The syntax of µF is the following: A program is a set of definitions. An expression is either a constant, a variable, a pair, an operator, a function call, a conditional, a local definition, an anonymous function, or one of the probabilistic operators sample, observe, or infer. The infer operator is tailored for ProbZelus and always takes two arguments: a transition function of the form fun x -> e, and a distribution of states. This operator computes the distribution of results and the distribution of possible next steps. A type system similar to the one of ProbZelus is used to distinguish deterministic from probabilistic expressions (see Appendix B.1).

Compilation to µF

The compilation function C generates a function that closely follows the transition function defined by the co-iterative semantics presented in Section 3.3. Each expression is compiled into a function of type S → T × S which given a state returns a value and an updated state (see Appendix C for the complete definition). The compilation of present is thus:

C(present e -> e 1 else e 2 ) = fun (s,s1,s2) -> let v,zs' = C(e)(s) in if v then let v1,s1' = C(e 1 )(s1) in (v1,(s',s1',s2)) else let v2,s2' = C(e 2 )(s2) in (v2,(s',s1,s2'))
The probabilistic operators sample, and observe are treated as external operators. The compilation generates code that simply calls the µF version of these operators. The compilation of infer passes the distribution over states to the µF version of infer. The inference is thus aware of the distribution over states at the previous step. The compilation of a node declaration generates two definitions: the transition function f _step and the initial state f _init. The transition function is the result of compiling the body of the node with an additional argument to capture the input. The initial step is generated by the allocation function A which follows the definition of the initial state in the semantics of Section 3.3 (see the Appendix C for the complete definition).

We can now prove that the semantics of the generated code corresponds to the semantics of the source language described in Section 3.3.

Theorem. For all ProbZelus expression e, for all state s and environment γ :

• if kindOf (e) = D then ⟦e⟧ s γ (s) = ⟦C(e)⟧ γ (s), and, • if kindOf (e) = P then {[e]} s γ (s) = {[C(e)]} γ (s). Proof.
The proof is done by induction on the structure of e. As an example consider the expression sample(e). If this expression is well-typed and since typing is preserved by compilation (see Lemma C.1) kindOf (C(e)) = D. Using the induction hypothesis on ⟦C(e)⟧ γ (s) = ⟦e⟧ γ (s) we have:

{[C(sample(e))]} γ (s) = fun s -> let mu, s' = C(e)(s) in let v = sample(mu) in (v, s') γ (s) = λU . ∫ δ ⟦C(e)⟧ γ (s) (dµ, ds ′ ) {[let v = sample(mu) in (v, s')]} γ [µ/mu, s ′ /s'] (U ) = λU . let (µ,s) = ⟦C(e)⟧ γ (s) in {[let v = sample(mu) in (v, s')]} γ [µ/mu, s ′ /s'] (U ) = λU . let (µ, s ′ ) = ⟦e⟧ s γ (s) in ∫ µ(dv)δ v ,s ′ (U ) = {[sample(e)]} s γ (s)
Remark. The probabilistic semantics of µF is commutative (see [START_REF] Staton | Commutative Semantics for Probabilistic Programming[END_REF]Theorem 4]). We can thus show that the semantics of a ProbZelus program does not depend on the schedule used by the compiler to order the equations of local definitions.

Inference

The measure-based semantics of infer presented in Section 3.3 and Section 4.3 includes often intractable integrals. An additional challenge is to design inference techniques that can operate in bounded memory to be practical in a reactive context where the inference is a non-terminating process.

In this section, we show how to adapt particle filtering [START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF] to explicitly handle the state of the transition functions. We then present a novel implementation of delayed sampling, a recently proposed semi-symbolic inference, which enables partial exact inference over infinite streams in bounded memory for a large class of models including state-space models like the robot example of Figure 1 in Section 2.

Particle Filtering

In conventional probabilistic programming, the operational interpretation of a model is an importance sampler that randomly generates a sample of the model together with an importance weight measuring the quality of the sample.

Following the conventions of Section 3.3 we write ⟦e⟧ γ for the semantics of a deterministic expression, and {[e]} γ ,w for the semantics of a probabilistic expression. The additional argument w captures the weight. The probabilistic operator sample draws a sample from a distribution without changing the score. observe increments the score by the likelihood of the observation. A deterministic expression can be lifted in a probabilistic context: the corresponding sample is the return value of the expression and the score is unchanged. The let construct illustrates that the score is accumulated along the execution path (the complete semantics is presented in Figure 19 of Appendix D).

{[sample(e)]} γ ,w = (draw(⟦e⟧ γ ), w) {[observe(e 1 ,e 2 )]} γ ,w = let µ = ⟦e 1 ⟧ γ in ((), w * µ pdf (⟦e 2 ⟧ γ )) {[e]} γ ,w = (⟦e⟧ γ , w) if kindOf (e) = D {[let p = e 1 in e 2 ]} γ ,w = let v 1 , w 1 = {[e 1 ]} γ ,w in {[e 2 ]} γ [v 1 /p],w 1
Such a sampler is the basis of a particle filter or a bootstrap filter [START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF]. infer independently launches N particles. At each step, each particle samples the distribution of states σ obtained at the previous step and executes the sampler to compute a pair (result, state) along with its weight. The resulting pairs are then normalized according to their weights to form a categorical distribution µ over pairs of values and states (we write w i = w i / N i=1 w i for the normalized weights). This distribution is then split into the distribution of returned values and the distribution of next states.

⟦infer(fun s -> e, σ )⟧ γ = let µ = λU . N i=1 let s i = draw(⟦σ ⟧ γ ) in let (v i , s ′ i ), w i = {[fun s -> e]} γ ,1 (s i ) in w i * δ v i ,s ′ i (U ) in (π 1 * (µ), π 2 * (µ))
Remark. The resampling step requires the ability to clone particles in the middle of the execution. A classic technique is to compile the model in continuation passing style (CPS) [START_REF] Ritchie | C3: Lightweight Incrementalized MCMC for Probabilistic Programs using Continuations and Callsite Caching[END_REF] and use the probabilistic constructs sample and observe as checkpoints for resampling. In our context, the compilation presented in Section 4 externalizes the state of the transition functions. Duplicating the state effectively clones a particle during its execution. The code does not need to be compiled in CPS form and we avoid the alignment problem [START_REF] Daniel Lundén | Automatic Alignment of Sequential Monte Carlo Inference in Higher-Order Probabilistic Programs[END_REF].

Delayed Sampling

The basis of our streaming inference algorithm is delayed sampling, which we first review to explain its conceptual approach. Delayed sampling is an inference technique combining partial exact inference with approximate particle filtering to reduce estimation errors [START_REF] Lundén | Delayed sampling in the probabilistic programming language Anglican[END_REF][START_REF] Murray | Delayed Sampling and Automatic Rao-Blackwellization of Probabilistic Programs[END_REF]. In addition to the importance weight, each particle exploits conjugacy relationships between pairs of random variables to maintains a graph: a Bayesian network representing closedform distributions involving subsets of random variables. Observations are incorporated by analytically conditioning the network. Particles are thus only required to draw sample when forced to, i.e., when exact computation is not possible, or when a concrete value is required.

{[op(e)]} γ ,д,w = let (e ′ , д e , w e ) = {[e]} γ ,д,w in (app(op, e ′ ), д e , w e ) {[if e then e 1 else e 2 ]} γ ,д,w = let e ′ , д e , w e = {[e]} γ ,д,w in let v, д v = value(e ′ , д e ) in if v then {[e 1 ]} γ ,д v ,w e else {[e 2 ]} γ ,д v ,w e {[sample(e)]} γ ,д,w = let µ, д e , w ′ = {[e]} γ ,д,w in let X, д ′ = assume(µ, д e ) in (X, д ′ , w ′ ) {[observe(e 1 ,e 2 )]} γ ,д,w = let µ, д 1 , w 1 = {[e 1 ]} γ ,д,w in let X, д x = assume(µ, д 1 ) in let e ′ 2 , д 2 , w 2 = {[e 2 ]} γ ,д x ,w 1 in let v, д v = value(e ′ 2 , д 2 ) in let д ′ = observe(X, v, д v ) in ((), д ′ , w 2 * µ pdf (v))
To perform analytic computations, delayed sampling manipulates symbolic terms where random variables are referenced in the graph. The semantics of an expression {[e]} γ ,д,w takes an additional argument д for the graph and returns a symbolic term, an updated weight, and an updated graph. Given a graph, a symbolic term can be evaluated into a concrete value by sampling the random variables that appear in the term. The graph can be accessed and modified using the three following functions defined in [START_REF] Murray | Delayed Sampling and Automatic Rao-Blackwellization of Probabilistic Programs[END_REF].

v, д ′ = value(e, д) evaluate a symbolic term and return a concrete value.

X, д ′ = assume(µ, д) add a random variable X ∼ µ to the graph and return the variable.

д ′ = observe(X, v, д) condition the graph by observing the value v for the variable X .

Compared to the particle filter, any expression, probabilistic or deterministic, can contribute to a symbolic term. The evaluation function {[e]} γ ,д,w partially presented in Figure 8 must thus be defined on the entire language and not only on probabilistic constructs. For instance, the application of an operator op(e) returns a symbolic term app(op, e ′ ) that represents the application of op on the evaluation of e. Some terms are partially evaluated when symbolic computation is not possible. For instance, in the general case, to compute the importance weight of if e then e 1 else e 2 , each particle must compute a concrete value for the condition e.

The probabilistic sample(e) adds a new random variable to the graph without drawing a sample. observe(e 1 , e 2 ) adds a new random variable X ∼ µ where µ is defined by e 1 , then computes a concrete value v for e 2 and conditions the graph by observing the value v for X . As for the particle filter, the score is incremented by the likelihood of the observation. Symbolic Computations. The functions value, assume, and observe used in Figure 8 rely on the following mutually recursive lower level operations (Y is the parent of X ) [START_REF] Murray | Delayed Sampling and Automatic Rao-Blackwellization of Probabilistic Programs[END_REF]:

X, д ′ = initialize(µ, Y , д) add a new node X with a distribu- tion p X |Y = µ as a child of Y in д.
д ′ = marginalize(X, д) compute and store p X in д ′ from p Y and p X |Y where p Y and p X |Y are in д.

д ′ = realize(X, v, д) assign in д ′ a concrete value to a random variable X .

д ′ = condition(Y , д) compute p Y |X given p X , p X |Y
, and a concrete value X = v where v is in д.

In the class of Bayesian networks maintained by the delayed sampler, marginalization w.r.t. a parent node, and conditioning a parent on the value of a child are tractable operations.

To reflect these operations, nodes are characterized by a state (see Figures 4 and9). Initialized nodes are random variables with a conditional distribution p X |Y where the parent Y has no concrete value yet. Marginalized nodes are random variables with a marginal distribution p X that incorporate the distributions of the ancestors. Realized nodes are random variables that have been assigned a concrete value via sampling or observation.

The evaluation function value(e, д) forces the realization by sampling of all the random variables referenced in e to produce a concrete value. Similarly, the function observe(X, v, д) realizes a variable X with a given observation v. The realization of a random variable comprises three steps: (1) compute the distribution p(x) by recursively marginalizing the parents from a root node, (2) sample a value, or use the observation, and (3) use the concrete value to update the children and condition the parent which removes the dependencies.

The function assume(µ, д) adds a new node to the graph and is defined case by case on the shape of the symbolic term µ. If there is a conjugacy relationship between µ and a random variable Y present in the graph, e.g., µ = Bernoulli(Y ) with Y ∼ Beta(α, β), a new initialized node X ∼ µ is added as a child of Y . Otherwise, since symbolic computation is not possible, dependencies are broken by realizing the random variables that appear in µ, e.g., µ ′ = Bernoulli(value(Y, д)), and X ∼ µ ′ is added as a new root node.

Inference. The inference scheme is similar to the particle filter. At each step, the inference draws N states from σ to execute the transition function. For each particle, execution starts with the graph computed at the previous step and returns a pair of symbolic terms (result, state), the particle . . . . . . . . .

⟦infer(fun s -> e, σ )⟧ γ = let µ = λU . N i=1 let s i , д i = draw(⟦σ ⟧ γ ) in let (e i , s ′ i ), w i , д ′ i = {[fun s -> e]} γ ,1,д i (s i ) in let d i = distribution(e i , д ′ i ) in w i * d i (π 1 (U )) * δ s ′ i ,д ′ i (π 2 (U )) in (π 1 * (µ), π 2 * (µ))

Streaming Delayed Sampling

As illustrated in Section 2.2, a notable challenge with the traditional delayed sampling algorithm is that graph grows linearly in the number of samples. In the original formulation of delayed sampling [START_REF] Murray | Delayed Sampling and Automatic Rao-Blackwellization of Probabilistic Programs[END_REF], graph edges are only removed when a node is realized. All nodes that have been neither sampled nor observed are thus kept in the graph even if they are no longer referenced by the program. In a reactive programming context, such an implementation can consume unbounded memory.

Bounded Delayed Sampling. A simple mitigation is to limit the scope of symbolic computations to one time step and discard the graph at the end of each time step. We call this inference technique bounded delayed sampling (BDS).

BDS performs symbolic computations during the execution of a time step, and whenever possible, delays the sampling until the end of the instant. Like the particle filter, BDS guarantees a bounded-memory execution. For each particle, the size of the graph is bounded by the number of variables introduced during a time step, which by construction, is bounded for any valid ProbZelus program.

Streaming Delayed Sampling. Compared to the original delayed sampling algorithm, BDS loses the ability to perform symbolic computations using variables defined at different time steps. This can result in a significant loss of precision for models with inter-steps dependencies such as the robot example of Figure 1. To adapt delayed sampling to streaming settings while keeping its maximum accuracy, we designed a delayed sampler that is pointer-minimal where nodes that are no longer referenced by the program can be eventually removed. We call this inference streaming delayed sampling (SDS). SDS enables partial exact inference in bounded memory for a large class of models.

In the original implementation of delayed sampling, graph nodes need to access their parents and children. Marginalization requires access to the parent to incorporate the ancestor distribution. Realization requires access to both the parent and the children of a node to update their respective distributions with the concrete value assigned to the node.

In the pointer minimal implementation, initialized nodes only keep a pointer to their parent to follow the ancestor chain during marginalization and marginalized nodes only keep a pointer to a marginalized child (Delayed Sampling imposes that a node always has at most one marginalized child). Compared to the original implementation, marginal nodes only keep track of one child, and marginalization turns backward pointers to the parent node into forward pointers to the marginalized child. Note that this implementation prevents updating the children when the parent is realized, and prevents conditioning a parent when a child is realized. Instead, when marginalizing a node, the sampler first checks if the parent is realized to apply the update. Symmetrically, to realize a node, the sampler first checks if the children are realized and, if necessary, conditions the distribution before assigning the concrete value.

Figure 9 shows the evolution of the graph during one step for the robot example of Figure 1. At the end of the step, the value of pre x is updated. The previous value is not referenced anymore by the program and the node can be removed from the graph. In the original implementation, backward pointers between marginalized nodes prevent the collection (see Figure 4).

Limitations. With SDS, models like the robot example that only maintain bounded chains of dependencies between variables are guaranteed to be executed in bounded memory. The class of models that can be executed in bounded memory with our pointer-minimal implementation thus already comprises state-space models like Kalman filters, and models for learning unknown constant parameters from a series of observations (e.g., computing the bias of a coin from a succession of flips) where variables introduced at each step are immediately realized.

However, unbounded chains can still be formed if the program keeps a reference to a constant variable that is never realized. In the following example, at each step, a new variable x is added as a child of pre x and then marginalized for the observation. But p1 keeps a reference to the initial variable i which is never realized and thus forms an unbounded chain between i and x. let proba p1 (xo, obs) = (i, o) where rec init i = sample (gaussian (xo, 1.)) and x = sample (gaussian (i -> pre x, 1.)) and () = observe (gaussian (x, 1.), obs)

In addition, in ProbZelus, at each step the inference returns a snapshot of the current distribution without forcing the realization of any node in the graph. Compared to the original delayed sampling implementation, initialized nodes can be inspected without being realized. It is thus possible to form unbounded chains of initialized nodes which cannot be pruned even when nodes are no longer referenced in the program due to the backward pointers to the parent in initialized nodes. In the following example, at each step, the variable x is added as a child of pre x, but without observation these variables remain initialized for ever. To mitigate these issues, we can force the realization of trailing nodes at each step as in bounded delayed sampling or use a sliding window. Alternatively, the value (eval) function is available to the programmer and can be used to implement any strategy to force the evaluation of the nodes. For instance, the previous example can be adapted to execute in bounded memory:

Evaluation

We next evaluate the performance of ProbZelus on a set of benchmarks that illustrate multiple aspects of the language: inferring fixed parameters from observations, online trajectory estimation, inference-in-the-loop. For these examples, we compare the accuracy and the latency cost of the three inference techniques: PF, BDS, and SDS. Appendix E details our implementation as an extension of the Zelus compiler. Benchmarks. The models used in the experiments are summarized in Table 1 (a detailed description along with the code of the benchmarks is given in Appendix F.1). Two models infer fixed parameters from a series of observations. Beta-Bernoulli estimates the parameter of a Bernoulli distribution from a series of binary observations (e.g., the bias of a coin).

Gaussian-Gaussian estimates the mean and variance of a Gaussian distribution from a series of observations. The accuracy metric is the Mean Squared Error (MSE) of the inferred parameters compared to their exact values. Two models infer the state of a moving agent from noisy observations. Kalman-1D is a one-dimensional Kalman filter that models an agent that estimates its trajectory from noisy observations. Outlier adapted from [START_REF] Thomas | Expectation Propagation for approximate Bayesian inference[END_REF] models the same situation as Kalman-1D, but the sensor occasionally produces invalid readings. This models infer both the trajectory of the agent, and the bias of the sensor. The accuracy metric is the Mean Squared Error (MSE) of the inferred trajectory compared to the exact positions.

Two models use inference-in-the-loop (IITL). Robot is the robot example of Section 2 and the accuracy metric is the LQR loss. SLAM (Simultaneous Localization and Mapping) adapted from [START_REF] Doucet | Rao-Blackwellised Particle Filtering for Dynamic Bayesian Networks[END_REF] models an agent that estimates its position and a map of its environment. In this simplified version, the agent moves in a one-dimensional grid where each cell is either black or white. The robot's wheels may slip causing the robot to unknowingly stay in place (noisy motion), and the sensor is not perfect and may accidentally report the wrong color (noisy observations). At each step the robot uses the inferred position to decide its next move. The accuracy metric is the MSE of both the position and the map.

MTT (Multi-Target Tracker) adapted from [START_REF] Murray | Automated learning with a probabilistic programming language: Birch[END_REF] is a model where there are a variable number of targets with linear-Gaussian motion models with a state space of 2D position and velocity, producing linear-Gaussian measurements of the position at each time step. Targets randomly appear according to a Poisson process and each disappear with fixed probability at each step. Measurements do not identify which target they came from, and "clutter" measurements that come not from targets but from some underlying distribution add Experimental Setup. All the experiments were run on a server with 32 CPUs (2.60 GHz) and 128 GB memory. We ran all the benchmarks for 500 steps. In all cases, the inference runs in bounded memory (see Appendix F.3).

For each algorithm, we evaluated how much time it required to achieve 90% of runs close to a loss target (out of 1000 runs total):

| log(P 90% (loss)) -log(loss target )| < 0.5.

For each benchmark, the baseline is the median loss of SDS at 1000 particles as loss target for that benchmark. We measured the number of particles required to achieve this loss, and then measured the total execution time at this particle count for 500 steps (in Appendix F.2 we also evaluate loss and step latency across a fixed range of particle counts).

Results. Figure 10 shows the results. The height of each bar is the median total execution time, and the error bars are 90% and 10% quantiles, aggregated over 1000 runs. Each bar is labeled with the minimum number of particles required to achieve the accuracy threshold, accurate to 1.5 significant digits (100, 150, 200, 250, . . . ). We observe that SDS is able to compute an exact solution for Beta-Bernoulli, Kalman-1D, and Robot. In all these examples 1 particle is already enough to reach the target accuracy. Overall, the results show that the number of particles required to reach the desired accuracy with PF implies a significant slowdown compared to SDS. Moreover, the SLAM and MTT benchmarks show that, in some cases, PF is not an option: the target accuracy was not reached with 15, 000 and 2, 500 particles, respectively, at which point PF was already 10 times slower than SDS and we stopped the experiments.

As expected, BDS performance numbers are between those of PF and SDS. At worst, when there is no possible intra-step symbolic computations (e.g., Beta-Bernoulli), BDS behaves like a particle filter and requires as many particles as PF. At best, BDS performs as well as SDS (e.g., Outlier).

Additionally, Figure 10 also shows that for a given number of particles, the overhead induced by managing the delayed sampling graph is significant. Compared to BDS and SDS, depending on the benchmark, it is possible to use 2 to 4 times as many particles for PF with the same execution time. However, this is not enough to match the gain in accuracy.

Alternative Baselines. The results presented in Figure 10 do not quantify the speedup of SDS on the SLAM and MTT benchmarks because the other inference algorithms time out. To evaluate speedups on these two benchmarks, we used PF as an alternative baseline instead of SDS. Figure 11 presents the execution time of PF, BDS, and SDS to reach a loss close to the median of PF with 2000 and 4000 particles.

We observe that SDS requires a much smaller number of particles to reach similar accuracy which translates into speedups ranging from 10 1 (MTT-2000) to 10 4 (SLAM-4000). BDS requires either a similar or smaller numbers of particles. But the overhead introduced by the graph manipulations mostly translates in slowdowns compared to PF.

Related Work

Probabilistic Programming. Over the last few years there has been a growing interest on probabilistic programming languages. Some languages like BUGS [START_REF] Lunn | The BUGS project: Evolution, critique and future directions[END_REF], Stan [START_REF] Carpenter | Stan: A probabilistic programming language[END_REF], or Augur [START_REF] Huang | Compiling Markov chain Monte Carlo algorithms for probabilistic modeling[END_REF] offer optimized inference technique for a constrained subset of models. Other languages like WebPPL [START_REF] Goodman | The Design and Implementation of Probabilistic Programming Languages[END_REF], Edward [START_REF] Tran | Deep Probabilistic Programming[END_REF], Pyro [START_REF] Bingham | Pyro: Deep Universal Probabilistic Programming[END_REF], or Birch [START_REF] Murray | Automated learning with a probabilistic programming language: Birch[END_REF] focus on expressivity allowing the specification of arbitrary complex models. Compared to these languages, ProbZelus can be used to program reactive models that typically do not terminate, and inference can be run in parallel with deterministic components that interact with an environment. Reactive Languages with Uncertainty. Lutin is a language for describing non-deterministic reactive systems for testing and simulation [START_REF] Raymond | Lutin: A Language for Specifying and Executing Reactive Scenarios[END_REF], but while Lutin supports weighted sampling to describe constrained random scenarios, it does not support inference. ProPL [START_REF] Pfeffer | Functional Specification of Probabilistic Process Models[END_REF] is a language to describe probabilistic models for process that evolve over a period of time. This language also extends a probabilistic language with a notion of processes that can be composed in parallel, but compared to ProbZelus, ProPL focuses on a constrained class of models that can be interpreted as Dynamic Bayesian Networks (DBN), and relies on standard DBN inference techniques. In the same vein, CTPPL [START_REF] Pfeffer | CTPPL: A Continuous Time Probabilistic Programming Language[END_REF] is a language to describe continuous-time processes where the amount of time taken by a sub-process can be specified by a probabilistic model. These models cannot be expressed in ProbZelus which relies on the synchronous model of computation. It would be interesting to investigate how to extend ProbZelus to continuoustime models based on Zelus' support for ordinary differential equations (ODE) [START_REF] Bourke | Zélus: a synchronous language with ODEs[END_REF].

Inference. Researchers have proposed streaming inference algorithms, including variational [START_REF] Broderick | Streaming Variational Bayes[END_REF], or sampling-based [START_REF] Doucet | Rao-Blackwellised Particle Filtering for Dynamic Bayesian Networks[END_REF][START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF] approaches. Popular languages like Stan, Edward, or Pyro, offer support to stream data through the model during inference to handle large datasets. However, compared to ProbZelus, the model must be defined a priori and does not evolve during the inference. The Anglican and Birch probabilistic programming languages support delayed sampling [START_REF] Murray | Delayed Sampling and Automatic Rao-Blackwellization of Probabilistic Programs[END_REF]. These languages do not support streaming inference or reactive programming. Again, their interfaces only support inference on a complete probabilistic model.

Conclusion

Modeling uncertainty is a primary element of control systems for tasks that operate under the assumption of a probabilistic model of their environment (e.g., object tracking).

While synchronous languages have developed as a prominent way to develop control applications, to date there has been limited work in these languages on programming language support for modeling uncertainty.

In this paper we present ProbZelus, the first synchronous probabilistic programming language that lifts emerging abstractions for probabilistic programming into the reactive setting thus enabling inference-in-the-loop. Moreover, our streaming delayed sampling algorithm provides efficient semi-symbolic inference while still satisfying a key requirement of control applications in that they must execute with bounded resources.

Our results demonstrate that ProbZelus enables us to write, in the very same source, a deterministic model for the control software and a probabilistic model for its behavior and environment with complex interactions between the two.

A ProbZelus

In this section, we provide the complete denitions of the ProbZelus type system and semantics for the kernel language introduced Section 3 extended with the probabilistic operator factor(e) which is equivalent to observe(exp(1), e). Intuitively, factor can directly update the weight of the execution path with the value of an expression e.

A.1 Typing

The type system that discriminates deterministic from probabilistic expressions is dened Figure 12. To simplify the presentation, we ignored datatypes polymorphism.

The sub-typing rule indicates that any deterministic expression can be lifted into a probabilistic one. Expressions like constants, variables, and last are deterministic. The kind of classic Zelus expressions (pairs, op, local denitions, present, and reset) is the kind of their body. Similarly, the kind of equations is the kind of their dening expression, and parallel composition imposes the same kind for all the equations. Note that it is always possible to compose deterministic and probabilistic computations. For rules where all sub-expressions share the same kind k we enforce the use of the sub-typing rule to lift deterministic expressions.

The expressions sample, factor, and observe are probabilistic. The transition from probabilistic to deterministic is realized via infer: a deterministic expression whose body is always probabilistic. Probabilistic expressions can thus only occur under an infer.

Other Static Analyses. The Zelus compiler statically checks initialization, and causality of the program [START_REF] Bourke | Zélus: a synchronous language with ODEs[END_REF]. These two analyses guarantee that there exists a schedule of parallel equations that makes the streams productive. Extending these analyses to the probabilistic operators is straightforward: probabilistic operators can be treated as external operators.

A.2 Co-iterative Semantics

The co-iterative semantics of ProbZelus's deterministic processes is inspired by [START_REF] Colaço | Mixing signals and modes in synchronous data-flow systems[END_REF] and dened Figure 13.

A node is a stream function of type T ! D T 0 . In addition to the state, the transition function thus takes an additional input of type T and returns a pair (result, next state)

CoNode(T ,T 0 , S) = S ⇥ (S ! T ! T 0 ⇥ S).
The transition function of a variable always returns the corresponding value stored in the environment . The semantics of last x is a simple access to a special variable x_last. present e -> e 1 else e 2 introduced in Section 2 returns the value of e 1 when e is true and the value of e 2 otherwise. The state (s, s 1 , s 2 ) stores the state of the three sub-expressions. The transition function lazily executes e 1 or e 2 depending on the value of e and returns the updated state.

The state of a set of scheduled locally recursive denitions e where rec E comprises three parts: the value of the local variables at the previous step which can be accessed via the last operator, the state of the dening expressions, and the state of expression e. The initialization stores the initial values introduced by init and the initial states of all subexpressions. The transition function incrementally builds the local environment dened by E. First the environment is populated with a set of fresh variables x i _last initialized with the values stored in the state that can then be accessed via the last operator. Then the environment is extended with the denition of all the variables i by executing all the dening expressions (where

{x i } 1..k \ { j } 1..n = {x i } 1..k ).
Finally, the expression e is executed in the nal environment. The updated state contains the value of the initialized variables dened in E that will the be used to start the next step, and the updated state of the sub-expressions.

Probabilistic Extensions. The semantics of the probabilistic part of ProbZelus, dened Figure 14, follows the same structure as the deterministic semantics but denes measures over all possible executions as in [39]. In particular a succession of computation is interpreted as sequentially integrating over the results of the preceding computations.

As for deterministic nodes, the transition function of a probabilistic node of type T ! P T 0 takes an additional argument and returns a measure over pairs (result, next state).

CoPNode(T ,T 0 , S) = S ⇥ (S ! T ! ( T 0 ⇥S ! [0, 1]))

A.3 Alternative semantics

We could give dierent semantics to ProbZelus. For example, consider the following probabilistic node. With the semantics dened Section 3, this program produces the stream of distribution: Beta(2, 1), Beta(3, 1), . . . Note that, even though, p is dened as a constant, its distribution evolves at each steps.

Since the observe statement uses the constant true, we know that p is necessarily 1. An alternative semantics could thus returns the constant stream of distributions: 1 .

B The µF language

Similarly to ProbZelus, we extend µF with the probabilistic operator factor. We now present the complete type system and semantics for µF .

B.1 Typing

The type system dened Figure 15 is similar to the one Figure 12 to distinguish deterministic from probabilistic expressions, but with additional restrictions since the compiled code is in a more constrained form. Whenever possible we require sub-expressions to be deterministic, that is, in pairs, operator applications (including sample, factor, and observe), function calls, and the condition of a if/then/else. These restrictions simplify the presentation of the semantics but do not reduce the expressiveness of the language since it is always possible to introduce additional local denitions to name intermediate probabilistic expressions. For example if sample(bernoulli(0.5)) then ... can be rewritten let b = sample(bernoulli(0.5)) in if b then ...

G `D e : t G `P e : t typeOf (c) = t G `D c : t G(x) = t G `D x : t G `k e 1 : t 1 G `k e 2 : t 2 G `k (e 1 , e 2 ) : t 1 ⇥ t 2 typeOf (op) = t 1 ! D t 2 G `k e : t 1 G `k op(e) : t 2 G(f ) = t 1 ! k t 2 G `D e : t 1 G `k f (e) : t 2 G(x) = t G `D last x : t G `k E : G 0 G + G 0 `k e : t G `k
G `D e : T dist ⇤ G `D e : T dist G `k e : t G `k x = e : [t/x] G `k e : t G `k init x = e : [t/x] G + G 1 + G 2 `k E 1 : G 1 G + G 1 + G 2 `k E 2 : G 2 G `k E 1 and E 2 : G 1 + G 2 G + [t 1 /x] `D e : t 2 G `D let node f x = e : G + [t 1 ! D t 2 /f ] G + [t 1 /x] `P e : t 2 G `D let proba f x = e : G + [t 1 ! P t 2 /f ] G `D d 1 : G 1 G 1 `D d 2 : G 2 G `D d 1 d 2 : G 2

B.2 Semantics of µF

The semantics of µF follows [39]. In a deterministic context kindOf (e) = D, the semantics neo of an expression is the classic interpretation of a strict functional language. In a probabilistic context (kindOf (e) = P), we dene a the measure-based semantics {[e]} .

The probabilistic semantics of µF is presented in Figure 16. A deterministic expression is lifted to a probabilistic expression using the the Dirac delta measure applied to the value of the expression computed by the deterministic semantics. As in Section 3.3, a local denition let x = e 1 in e 2 is interpreted as integrating e 2 over the measure dened by e 1 . The semantics of the probabilistic operators is the following: sample(e) returns the distribution neo . factor(e) returns a measure dened on the singleton space () whose value is exp(neo ). observe(e 1 , e 2 ) is similar but the score is the density function of the distribution ne 1 o applied to ne 2 o .

Inference. infer handles the transition function generated by the compilation of Section 4. The rst argument of infer is a transition function, and the second argument a distribution over state . The inference rst integrates over the distribution and then normalizes the result µ to produce a distribution of pairs (result, next state). The special value > denotes the entire space (value, state). This distribution is then decomposed into a pair of distributions using the pushforward of µ.

C Compilation

Figure 18 presents the entire compilation function from ProbZelus to µF introduced Section 4. Figure 17 presents the allocation function.

Lemma C.1. The compilation preserves the kind (deterministic D, or probabilistic P) of the expressions. For any expression e, if G `k e : t, there exists G 0 and t 0 such that G 0 `k C(e) : t 0 . Proof. By induction on the structure of e. ⇤

Remark. The compilation presented in Figure 18 generates a function for each sub-expression. However, in most cases it is possible to simplify the code using static reduction. For instance, a constant can directly be compiled into a constant.

nco i = () nco s = s. (c, s) nxo i = () nxo s = s. ( (x), s) nlast xo i = () nlast xo s = s. ( (x_last), s) n(e 1 ,e 2 )o i = (ne 1 o i , ne 2 o i ) n(e 1 ,e 2 )o s = (s 1 , s 2 ). let 1 , s 0 1 = ne 1 o s (s 1 ) in let 2 , s 0 2 = ne 2 o s (s 2 ) in (( 1 , 2 ), (s 0 1 , s 0 2 )) nop(e)o i = neo i nop(e)o s = s. let , s 0 = neo s (s) in (op( ), s 0 ) nf (e)o i = (neo i , (f _init)) nf (e)o s = (s 1 , s 2 ). let 1 , s 0 1 = neo s (s 1 ) in let 2 , s 0 2 = (f _step)( 1 )(s 2 ) in ( 2 , (s 0 1 , s 0 2 )) npresent e -> e 1 else e 2 o i = (neo i , ne 1 o i , ne 2 o i ) npresent e -> e 1 else e 2 o s = (s, s 1 , s 2 ). let , s 0 = neo s (s) in if then let 1 , s 0 1 = ne 1 o s (s 1 ) in ( 1 , (s 0 , s 0 1 , s 2 )) else let 2 , s 0 2 = ne 2 o s (s 2 ) in ( 2 , (s 0 , s 1 , s 0 2 )) nreset e 1 every e 2 o i = (ne 1 o i , ne 1 o i , ne 2 o i ) nreset e 1 every e 2 o s = (s 0 , s 1 , s 2 ). let 2 , s 0 2 = ne 2 o s (s 2 ) in let 1 , s 0 1 = ne 1 o s (if 2 then s 0 else s 1 ) in ( 1 , (s 0 , s 0 1 , s 0 2 )
let 1 = [m 1 /x 1 _last] in . . . let k = k 1 [m k /x k _last] in let 1 , s 0 1 = ne 1 o s k (s 1 ) in let 0 1 = k [ 1 / 1 ] in . . . let n , s 0 n = ne n o s 0 n 1 (s n ) in let 0 n = 0 n 1 [ n / n ] in let , s 0 = neo s 0 n (s) in , (( 0 n [x 1 ], . . . , 0 n [x k ]), (s 0 1 , . . . , s 0 n ), s 0 ) nlet node f x = eo = [neo i /f _init, . s. neo s [ /x] /f _step] nlet proba f x = eo = [{[e]} i /f _init, . s. {[e]} s [ /x] /f _step] nd 1 d 2 o = let 1 = nd 1 o in nd 2 o 1 Figure 13.
Co-iterative semantics of deterministic ProbZelus programs. For local denitions each initialized variable is dened in a subsequent equation, i.e.,

{x i } 1..k \ { j } 1..n = {x i } 1..k . i = neo i if kindOf (e) = D {[e]} s = s. U . neo s (s) (U ) if kindOf (e) = D = s. U . ⇢ 1 if neo s (s) 2 U 0 otherwise {[(e 1 ,e 2 ) ]} i = ({[e 1 ]} i , {[e 2 ]} i ) {[(e 1 ,e 2 ) ]} s = (s 1 , s 2 ). U . let µ 1 = {[e 1 ]} s (s 1 ) in Ø µ 1 (d 1 , ds 0 1 ) let µ 2 = {[e 2 ]} s (s 2 ) in Ø µ 2 (d 2 , ds 0 2 ) ( 1 , 2 ),(s 0 1 ,s 0 2 ) (U ) {[op(e)]} i = {[e]} i {[op(e)]} s = s. U . let µ = {[e]} s (s) in Ø µ(d , ds 0 ) op( ),s 0 (U ) {[ f (e)]} i = (neo i , (f _init)) {[ f (e)]} s = (s 1 , s 2 ). U . let 1 , s 0 1 = neo s (s 1 ) in let µ 2 = (f _step)( 1 )(s 1 ) in Ø µ 2 (d 2 , ds 0 2 ) 2 ,(s 0 1 ,s 0 2 ) (U ) {[present e -> e 1 else e 2 ]} i = ({[e]} i , {[e 1 ]} i , {[e 2 ]} i ) {[present e -> e 1 else e 2 ]} s = (s, s 1 , s 2 ). U . let µ = {[e]} s (s) in Ø µ(d , ds 0 ) if then let µ 1 = {[e 1 ]} s (s 1 ) in Ø µ 1 (d 1 , ds 0 1 ) 1 ,(s 0 ,s 0 1 ,s 2 ) (U ) else let µ 2 = ne 2 o s (s 2 ) in Ø µ 2 (d 2 , ds 0 2 ) 2 ,(s 0 ,s 1 ,s 0 2 ) (U ) {[reset e 1 every e 2 ]} i = ({[e 1 ]} i , {[e 1 ]} i , {[e 2 ]} i ) {[reset e 1 every e 2 ]} s = (s 0 , s 1 , s 2 ). U . let µ 2 = ne 2 o s (s 2 ) in Ø µ(d 2 , ds 0 2 ) let µ 1 = ne 1 o s (if 2 then s 0 else s 1 ) in Ø µ(d 1 , ds 0 1 ) 1 ,(s 0 ,s 0 1 ,s 0 2 ) (U ) 8 > > > > > > > < > > > > > > > :
9 > > > > > > > = > > > > > > > ; i = © ≠ ´(c 1 , . . . , c k ), ({[e 1 ]} i , . . . , {[e n ]} i ), {[e]} i ™ AE 8 > > > > > > > < > > > > >
9 > > > > > > > = > > > > > > > ; s = ((m 1 , . . . ,m k ), (s 1 , . . . , s n ), s). U . let 1 = [m 1 /x 1 _last] in . . . let k = k 1 [m k /x k _last] in let µ 1 = {[e 1 ]} s k (s 1 ) in Ø µ 1 (d 1 , ds 0 1 ) let 0 1 = k [ 1 / 1 ] in Ø . . . let µ n = ne n o s 0 n 1 (s n ) in Ø µ n (d n , ds 0 n ) let 0 n = 0 n 1 [ n / n ] in let µ = {[e]} s 0 n (s) in Ø µ(d , ds 0 ) ,(( 0 n [x 1 ],..., 0 n [x k ]),(s 0 1 ,...,s 0 n ),s 0 ) (U ) {[sample(e)]} i = neo i {[sample(e)]} s = s. U . let µ, s 0 = neo s (s) in Ø T µ(d ) ,s 0 (U ) {[factor(e)]} i = neo i {[factor(e)]} s = s. U . let , s 0 = neo s (s) in exp( ) (),s 0 (U ) {[observe(e 1 , e 2 )]} i = (ne 1 o i ,ne 2 o i ) {[observe(e 1 , e 2 )]} s = (s 1 , s 2 ). U . let µ, s 0 1 = ne 1 o s (s 1 ) in let , s 0 2 = ne 2 o s (s 2 ) in µ pdf ( ) ⇤ (),(s 0 1 ,s 2 ) (U )
Figure 14. Co-iterative semantics of probabilistic ProbZelus expressions (i.e., kindOf (e) = P). For local denitions each initialized variable is dened in a subsequent equation, i.e.,

{x i } 1..k \ { j } 1..n = {x i } 1..k . G `D e : t G `P e : t typeOf (c) = t G `D c : t G(x) = t G `D x : t G `D e 1 : t 1 G `D e 2 : t 2 G `D (e 1 , e 2 ) : t 1 ⇥ t 2 typeOf (op) = t 1 ! D t 2 G `D e : t 1 G `D op(e) : t 2 G(f ) = t 1 ! k t 2 G `D e : t 1 G `k f (e) : t 2 G + [t 1 /x] `k e 1 : t 2 G `D e 2 : t 1 G `k (fun x -> e 1 )(e 2 ) : t 2 G `D e : bool G `k e 1 : t G `k e 2 : t G `k if e then e 1 else e 2 : t G `k e 1 : t 1 G + [t 1 /x] `k e 2 : t 2 G `k let x = e 1 in e 1 : t 2 G + [t 1 /x] `k e : t 2 G `D fun x -> e : t 1 ! k t 2 G `D e : t dist
G `P e 1 : t ⇥ t state G `D e 2 : t state dist G `D infer((fun x -> e 1 ),e 2 ) : t dist G `D e : T dist ⇤ G `D e : T dist G `D e : t G `D let f = e : G + [t/f ] G `D d 1 : G 1 G 1 `D d 2 : G 2 G `D d 1 d 2 : G 2 Figure 15.
Typing of µF with deterministic and probabilistic kinds.

{[let f = e]} = [{[e]} /f ] {[d 1 d 2 ]} = let 1 = {[d 1 ]} in {[d 2 ]} 1 {[e]} = U . neo (U ) if kindOf (e) = D {[e 1 (e 2 )]} = U . (ne 1 o (ne 2 o ))(U ) {[let p = e 1 in e 2 ]} = U . Ø T {[e 1 ]} (du){[e 2 ]} +[u/p] (U ) {[if e then e 1 else e 2 ]} = U . if neo then {[e 1 ]} (U ) else {[e 2 ]} (

D Inference D.1 Importance Sampling

Importance sampling. The most simple inference independently launches N particles. Each particle executes the importance sampler to compute a pair (result, weight). Results are then normalized in a categorical distribution, i.e., a discrete distribution over the results.

Õ i=1 let s i , w i = nSo [i] in let ( i , s 0 i ), w 0 i = {[fun s -> e]} ,w i (s i ) in w 0 i ⇤ i (U ) in (µ, [(s 0 i , w 0 i )] 1i N )
At each step, the inference executes one step of all the particles and normalizes the scores to return the distribution µ of possible results and an updated array of pairs (state, weight) for the next step. The weights of the particles are multiplied at each step and never reset. In other words, the inference reports at each step how likely is the execution path since the beginning of the program for each particle w.r.t. the model. Obviously the probability of each individual path quickly collapses to 0 after a few steps which makes this inference technique not practical in a reactive context where the inference process never terminates. The particle lter mitigates this issue by periodically re-sampling the set of particles.

E Implementation

ProbZelus is open source (hps://github.com/IBM/probzelus). It is implemented on top of Zelus (hp://zelus.di.ens.fr/). The new constructs sample, observe, and factor are Zelus nodes implemented directly in OCaml. The infer construct is a node that take as argument the Zelus node that represents the probabilistic model. The infer node thus takes as argument the allocation and step functions of the model as argument which corresponds to the compilation described in Section 4.

Relationship with the paper The code corresponding to the paper is available as a release hps://github.com/IBM/ probzelus/tree/pldi20. The example of Figure 1 is in examples /tracker/tracker_ds.zls.

The compiler implements the compilation scheme presented in Section 4 with a few optimizations: (1) intermediate step functions are statically reduced (2) useless state is removed when possible, and (3) state is updated imperatively.

Moreover, the compilation of proba nodes introduces an extra argument to the functions in order to pass the extra information w or (w, ) needed by the inference algorithms.

The code of the inferences algorithms is in the inference directory. The particle lter presented in Section 5.1 is in infer_pf.ml. The entry point of the Delayed Sampling algorithm presented in Section 5.2 is infer_ds_naive.ml and the core of the algorithm is ds_naive_graph.ml. The entry point and the core of the algorithm for the Streaming Delayed Sampling algorithm presented in Section 5.3 are respectively in infer_ds_streaming.ml and ds_streaming_graph.ml.

The Bounded Delayed Sampling algorithm presented in Section 5.2 can be implemented on top of both classical and streaming delayed sampling. The code is in the functor dened in ds_high_level.ml.

Finally, the code for the benchmarks presented Section 6 and Appendix F is available in examples/benchmarks. Artifact There is an artifact associated to the paper which is available with [START_REF] Baudart | Reactive Probabilistic Programming[END_REF]. It is distributed as a Linux image in the Open Virtualization Format that can be launch using a virtualization player like VirtualBox (hps://www.virtualbox. org). The credential to log into the virtual machine are: user: probzelus passord: probzelus

F Performance Evaluation

This section presents the experimental results. We ran each inference algorithm on a series of benchmarks and measured properties of the execution: accuracy, execution time, memory consumption. All the experiments were run on a server with 32 CPUs (2.60 GHz) and 128 GB memory.

F.1 Benchmarks

Beta-Bernoulli. The Beta-Bernoulli benchmark models an agent that estimates the bias of a coin. The model samples zt from a Beta(1, 1) distribution, and thereafter evaluates the observations with a Bernoulli distribution of parameter xt. Running SDS on this model is equivalent to exact inference in a Beta-Bernoulli conjugate model [START_REF] Goodman | The Design and Implementation of Probabilistic Programming Languages[END_REF] where each particle returns the exact solution. The benchmark's error metric is the mean squared error over time between the true coin probability and the expected probability conditioned on the stream of observations.

Gaussian-Gaussian. The Gaussian-Gaussian benchmark models an agent that estimates the mean and the standard deviation of a Gaussian. The initial values for the distribution of the mean mu follows a distribution N(0, 10) and the distribution of p sigma is N(0, 1). The distributions of mu and sigma are conditioned by the observations that follow a distribution N(mu, sigma). In the current implementations of delayed sampling we are doing exact inference only on the mean and not on the standard deviation (even if it would be possible). The benchmark's error metric is the mean squared error over time between the true mean and standard deviation and the expected probability conditioned on the stream of observations. Kalman. The Kalman benchmark models an agent that estimates its position based on noisy observations. The model chooses an initial position from N(0, 2500), and chooses subsequent positions from N(pre x, 1) where pre x denote the previous position. The model draws the observation at each time step from N(x, 1) where x is the true position. Running SDS on this model is equivalent to a Kalman lter [START_REF] Lunn | The BUGS project: Evolution, critique and future directions[END_REF] where each particle returns the exact solution. The benchmark's error metric is the mean squared error over time between the true position and the expected position conditioned on all previous observations. Outlier. The Outlier benchmark, adapted from Section 2 of [START_REF] Narayanan | Probabilistic Inference by Program Transformation in Hakaru (System Description)[END_REF], models the same situation as the Kalman benchmark, but with a sensor that occasionally produces invalid readings. The model chooses the probability of an invalid reading from a Beta(100, 1000) distribution, so that invalid readings occur approximately 10% of the time. At each time step, with the previously chosen probability, the model either chooses the observation from the invalid distribution N(0, 10000), or it chooses the observation from the Kalman model. Running SDS on this model is equivalent to a Rao-Blackwellized particle lter [START_REF] Gehr | PSI: Exact Inference for Probabilistic Programs[END_REF] that combines exact inference with approximate particle ltering. The benchmark's error metric is the mean squared error over time between the true position and the expected position conditioned on all previous observations. Robot. The Robot benchmark is detailed Section 2.

SLAM. Simultaneous Location And Mapping (SLAM) [START_REF] Pfeffer | Functional Specification of Probabilistic Process Models[END_REF]. Consider the simple case where a robot evolves in a discrete one-dimensional world and each position corresponds to a black or white cell. The robot can move from left to right and can observe the color of the cell on which it stands with a sensor. There are two sources of uncertainty: (1) the robot's wheels are slippery, so the robot can sometimes stay on the spot thinking about moving, (2) the sensor is making read errors, and can reverse the colors. The controller tries to infer the map (color of the cells) and the current position of the robot (Figure 20).

The robot maintains a map where each box is a random variable that represents the probability of being black or white (gray level in the Figure 20). The a priori distribution of these random variables is a Bernoulli(0.5) distribution: The sensor has a constant probability of reading error of sensor_noise. At each instant, the robot computes its current position x. The observation of the sensor follows a Bernoulli distribution parameterized by 1 -sensor_noise if the position is white and sensor_noise if the position is black. The benchmark's error metrics is the mean squared error over time between the exact map and position and the expected map and position.

Multi-Target Tracker. MTT (Multi-Target Tracker) adapted from [START_REF] Raymond | Lutin: A Language for Specifying and Executing Reactive Scenarios[END_REF] is a model where there are a variable number of targets with linear-Gaussian motion models producing linear-Gaussian measurements of the position at each time step. Targets randomly appear according to a Poisson process and each disappear with some xed probability at each time step. Measurements do not identify which target they came from, and "clutter" measurements that come not from targets but from some underlying distribution add to observations, complicating inference of which measurements are associated to which targets.

We model this with a ProbZelus program that has a state consisting of a list of position-velocity pairs that encode the track of each target. In this example we consider twodimensional targets, giving us a 4-dimensional vector representing position and velocity together.

The rst step is to dene helper functions that will be mapped over the list of tracks. The rst function tells us how frequently tracks die. They do so with probability p_dead which we set to e .02 . let proba death_fn _ = sample (bernoulli (p_dead))

We now dene how tracks are initialized when they are rst created. They are sampled from a multivariate Gaussian distribution with mean mu_new set to zero and covariance sigma_new set to a diagonal with variance 1 on the positions and variance 0.001 on the velocities. Next, we dene the motion model and update model. Each track tr is multiplied with the motion matrix a_u which encodes discrete time integration of position and velocity with time constant 1. We then sample a Gaussian distribution around the new position and velocity with covariance sigma_update which is a diagonal matrix with 0.01 variance for the position and 0.1 for the variance of the velocity. For observation model, we project out the position with the projection matrix proj_pos and observe it with covariance matrix sigma_obs which we set to a diagonal of 0.1. let proba state_update_fn (tr_num, tr) = (tr_num, sample (mv_gaussian(a_u *@ tr, sigma_update))) let observe_fn (_, tr) = (mv_gaussian (proj_pos *@ tr, sigma_obs))

We next dene the model for clutter data. We assume that each clutter point is drawn from a multivariate Gaussian with mean mu_clutter which is zero an covariance sigma_clutter which we set to 10. The model proceeds as follows. For every track, we use the filter list operator to remove all the tracks that died in this time step. We then sample the number of new tracksn_new from a Poisson distribution with parameter lambda_new which we set to 0.1. After forcing a sample of this value, we use the list constructer ini to build a list of new tracks. We append the survived and new tracks together, and then use the map list operator to apply the motion and observation models to each track. Next, we determine the amount of clutter by subtracting the number of observations from the number of surviving tracks. We then observe that this comes from a Poisson distribution with parameter lambda_clutter which we set to 1. Note that this sets the particle weight to 1 if the particle yields a negative amount of clutter. Next, we shue the track observations and the clutter together by forcing a sample of the shuffle random primitive, and nally observe that the resulting list yields the observed values. The accuracy metric is based on the Multiple Object Tracking Accuracy [START_REF] Bernardin | Evaluating Multiple Object Tracking Performance: The CLEAR MOT Metrics[END_REF]. This evaluates whether a track estimate contains the right targets across all time steps within a sufcient tolerance (we set the tolerance to 5 in our example). Conventional MOTA is in [0, 1] with 1 being the best; we have modied it to be in [0, 1] with 0 being the best by transforming it to MOTA⇤ = 1/MOTA 1.

Because we estimate a distribution of track estimates, we draw a sample from the track distribution to estimate the expected MOTA⇤.

Data. For each benchmark except Robot and SLAM, we obtained observation data by sampling from the benchmark's model. In these benchmarks, every run of each benchmark across all experiments uses the same data as input. For SLAM, we pre-sampled the map from the model, but sampled position data on the y as this data depends on the controller. For the Robot benchmark, we sampled all observations on the y because they all depend on the command from the controller. This means that for SLAM and Robot, each run uses dierent position observations.

F.2 DS vs. PF

We compare both the accuracy and runtime performance of BDS, SDS, and PF to investigate whether the delayed samplers can achieve better accuracy than the particle lter with the same amount of computational resources.

Accuracy Methodology. For a range of selected particle counts, we execute each benchmark multiple times and record the resulting accuracy. To measure accuracy we use the endto-end error metrics for each benchmark as described in Section F.1. We record the median and the 90% and 10% quantiles after 1000 runs.

Accuracy Results. Figures 21 and23 show the results of the accuracy experiment for the dierent benchmarks. The error bars show 90% and 10% quantiles, and the center is the median. The vertical lines corresponds to the number reported in Figure 10 where there is enough particles to achieve similar accuracy to delayed sampling with 1000 particles. In all cases, SDS is able to achieve equal or better accuracy than BDS which is itself equal or better than PF, but the results vary widely by benchmark. Note that SDS returns the exact posterior distribution for the Coin and Kalman benchmarks therefore its accuracy is independent of the number of particles. On the other-hand, BDS is not exact since the symbolic distributions are sampled at the end of each the step.

Performance Methodology. For a range of selected particle counts, we execute each benchmark multiple times (the same number as for the accuracy experiments described above) after a warm-up of 1 run and record the resulting performance: the latency of one step of computation. In the following graphs we report the median latency as well as the 90% and 10% quantiles of the collected data.

Performance Results. Figures 22 and24 shows how the latency for a single step varies with the number of particles each benchmark. The error bars show 90% and 10% quantiles, and the center is the median. With the three algorithms, the execution time increases linearly with the number of particles. In all cases, PF has lower latency than BDS which has lower latency than SDS.

Conclusions. These experiments show that the delayed samplers achieve better accuracy than the particle lter with the same computational resources. For some models SDS is able to compute the exact solution with only one particle (Kalman, Coin). BDS achieves better accuracy when relationships between variables dened in the same step can be exploited (Kalman). At worst the delayed samplers performs as a well as the particle lter (BDS on the Coin, SDS and BDS on the Outlier).

F.3 SDS vs. DS

We next evaluate the performance of SDS and BDS relative to our own OCaml implementation of the original delayed sampler (DS). We compare both the performance and memory consumption of the three algorithms at each time step to investigate whether, as the size of the input stream grows large, they can retain constant performance.

Performance Methodology. We execute each benchmark 1000 times after a warm-up of 1 run and record the latency. We execute each benchmark with 100 particles (even if only one particle is necessary for DS and SDS on the Coin and Kalman benchmarks to compute the exact distribution) and plot latency as a function of the time step. We report the median latency as well as the 90% and 10% quantiles of the collected data.

Performance Results. Figures 25 and27 shows the latency at each step of a run, aggregated over 1000 runs. PF, BDS, and SDS show nearly constant performance in time but DS gets linearly worse performance for the Kalman and Outlier benchmarks. For the Coin benchmark, the graph of DS remains of constant size because there is only one sample at the rst step and then only observe statements.

Memory Methodology. We next evaluate the memory consumption of the algorithms. For all benchmarks except the multi-target tracker, memory consumption is deterministic even in the presence of random choices. Therefore, we measure the ideal memory consumption of the execution of each benchmark after each step. The ideal memory consumption is the total amount of live words in the program's heap. In our implementation, we measure these numbers by forcing a garbage collection after each step. We use OCaml's standard facilities for forcing garbage collection as well as for measuring the amount of live words. We ran each algorithms 10 times with 100 particles.

For the multi-target tracker, the memory is not deterministic because it is determined by the number of hypothetical tracks, which is random. We report median and 10% and 90% values for memory consumption for this benchmark.

Memory Results. Figures 26 and28 shows the results of the memory consumption experiment. For all benchmarks, PF, BDS, SDS use constant memory over time, including for the multi-target tracker where their memory consumption is random at each time step. However, DS has increasing memory consumption over time for the Kalman, Outlier, and Robot benchmarks. The memory consumption of DS is constant for the Coin benchmark because the graph remains of constant size.

For the mutli-target tracker, the memory consumption of DS is based both on the number of hypothesized tracks and the length of the hypothesized tracks. We can see that the memory consumption of DS increases as the rst generation of tracks becomes longer, but eventually curtails its memory consuption when these tracks die. MTT's memory consumption thereafter increases again as the second generation of tracks starts to increase in length.

Conclusions. The original DS implementation consumes an increasing amount of memory over time for models that introduce new variables at each step (Kalman, Outlier, and Robot) in contrast to BDS and SDS whose memory consumptions are constant over time. For the multi-target tracker, the DS memory consumption is based on the length of the track which is in principle probabilistically bounded. However, DS still consumes much more memory than PF, BDS, and SDS because the tracks are long-lived.

Furthermore, DS step latency increases without bound as the number of steps becomes large on benchmarks where the memory increases. These observations conrm that the original DS implementation is not practical in a reactive settings.

⌅ PF

• BDS H SDS

Figure 3 .

 3 Figure 3. Particle filter (PF) and streaming delayed sampling (SDS) performances for the robot example of Figure 1. Accuracy is measured using the loss function of the LQR. Speed corresponds to the execution time of 500 steps.

step 4 Figure 4 .

 44 Figure 4. Evolution of the delayed sampling graph for the model of Figure 1. Each node denotes either a value (dark gray) or a distribution (light gray). Plain arrows represent dependencies in the underlying Bayesian network. The dotted arrow represents the pointers in the original data-structure implementing the graph. Labels indicate the program variables.

Figure 5 .

 5 Figure 5. Delayed sampling (DS) and streaming delayed sampling (SDS) memory consumption in thousands of live words in the heap per steps for the robot example.

d

  ::= let node f x = e | let proba f x = e | d d e ::= c | x | (e,e) | op(e) | f (e) | last x | e where rec E | present e -> e else e | reset e every e | sample(e) | observe(e, e) | infer(e) E ::= x = e | init x = c | E and E

  let node cpt () = o where rec o = 0 -> pre o + 1 let node present_vs_if (b) = (o1, o2) where rec o1 = present (b) -> cpt () else 0 and o2 = if b then cpt () else 0 b true true false true false false true ...

e

  where rec init x 1 = c 1 ... and init x k = c k and y 1 = e 1 ... and y n = e n

  of Appendix A.1). G ⊢ D e : T dist G ⊢ P sample(e) : T G ⊢ D e 1 : T dist * G ⊢ D e 2 : T G ⊢ P observe(e 1 , e 2 ) : unit G ⊢ D e : T G ⊢ P e : T G ⊢ P e : T G ⊢ D infer(e) : T dist G ⊢ D e : T dist * G ⊢ D e : T dist The type T dist represents distributions over values of type T . Distributions can be sampled (sample statement) and analyzed with external operators such as mean and variance.

  ⟦x⟧

)) e where rec init x 1 ==((c 1

 11 c 1 ... and init x k = c k and y 1 = e 1 ... and y n = e n i γ , . . . , c k ), (⟦e 1 ⟧ i γ , . . . , ⟦e n ⟧ i γ ), ⟦e⟧ i γ ) e where rec init x 1 = c 1 ... and init x k = c k and y 1 = e 1 ... and y n = e n s γ = λ((m 1 , . . . ,m k ), (s 1 , . .

C

  (infer(e)) = fun sigma -> let mu, sigma' = infer(C(e), sigma) in (mu, sigma')

Figure 8 .

 8 Figure 8. Delayed sampling sampler. Expressions return a pair (symbolic expression, weight).

  init(x, pre x)

Figure 9 .

 9 Figure 9. One step of the robot example of Figure1with SDS. Plain arrows represent dependencies and dotted arrows represent pointers at runtime. The sample statement adds the initialized nodes x (b). The observe statement adds the initialized node a (c), triggers the marginalizations of x (d) and a (e), and assigns to a its observed value (f). When the state is updated, the value x becomes pre x. The previous values are not referenced anymore and can be removed (g).

  let proba p2 (xo) = x where rec x = sample (gaussian (xo -> pre x, 1.))

  let proba p2' (xo) = x where rec x = sample (gaussian (xo -> pre x, 1.)) and _ = eval (xo -> pre x)

Figure 10 .

 10 Figure 10. Execution time comparison when 90% of 1000 runs reach an accuracy similar to the baseline (median accuracy of SDS with 1000 particles) after 500 steps. The number of particles required to reach this accuracy is shown on top of the bars. The error bars show the 10th and 90th percentiles.

Figure 11 .

 11 Figure 11. Execution time comparison with two different baselines: median accuracy of PF with 2000 and 4000 particles, respectively.

  let proba kahn_vs_scott () = p where rec init p = sample(beta(1, 1)) and () = observe(bernoulli(p), true)

  e where rec E : t G `k e : bool G `k e 1 : t G `k e 2 : t G `k present e -> e 1 else e 2 : t G `k e 1 : t G `k e 2 : bool G `k reset e 1 every e 2 : t G `D e : t dist G `P sample(e) : t G `D e 1 : t dist ⇤ G `D e 2 : t G `P observe(e 1 , e 2 ) : unit G `D e : float G `P factor(e) : unit G `P e : t G `D infer(e) : t dist

Figure 12 .

 12 Figure 12. Typing with deterministic and probabilistic kinds.

) ⇢ ⌧ e where rec init x 1 =

 1 c 1 and ... and init x k = c k and 1 = e 1 and ... and n = e n i = © ≠ ´(c 1 , . . . , c k ), (ne 1 o i , . . . , ne n o i ), neo i init x 1 = c 1 and ... and init x k = c k and 1 = e 1 and ... and n = e n s = ((m 1 , . . . ,m k ), (s 1 , . . . , s n ), s).

rec init x 1 =

 1 c 1 and ... and init x k = c k and 1 = e 1 and ... and n = e n

G

  `P sample(e) : t G `D e 1 : t dist ⇤ G `D e 2 : t G `P observe(e 1 , e 2 ) : unit G `D e : float G `P factor(e) : unit

Figure 16 .

 16 Figure 16. Probabilistic semantics of µF . The semantics is dened only for probabilistic expressions (kindOf (e) = P).

AFigure 17 .

 17 Figure 17. Memory allocation, i.e., initialization for the µF step functions.

N

Figure 18 .Figure 19 .

 1819 Figure 18. Compilation of ProbZelus to µF .

  let proba coin (yobs) = xt where rec init xt = sample (beta (1., 1.)) and () = observe (bernoulli xt, yobs)

  let proba gaussian_model (o) = (mu, sigma) where rec init mu = sample (gaussian (0., 10.)) and init sqrt_sigma = sample (gaussian (0., 1.)) and sigma = sqrt_sigma *. sqrt_sigma and () = observe (gaussian (mu, sigma), o)

  let proba delay_kalman (yobs) = xt where rec xt = sample (gaussian ((0., 2500.) -> (pre xt, 1.))) and () = observe (gaussian (xt, 1.), yobs)

  let proba outlier (yobs) = (is_outlier, xt) where rec xt = sample (gaussian ((0., 2500.) -> (pre xt, 1.))) and init outlier_prob = sample (beta (100., 1000.)) and is_outlier = sample (bernoulli outlier_prob) and () = present is_outlier -> observe (gaussian (0., 10000.), yobs) else observe (gaussian (xt, 1.), yobs)

Figure 20 .

 20 Figure 20. Screenshots of the execution of the SLAM with the PF and SDS inferences. For each screenshot, the top line shows the map, and the circle the exact position of the robot. The lower line represents the inferred map where the gray level indicates the probability for the cell to be black and the red dots the probability of presence of the robot on the cell.

  let proba bernoulli_priors i = sample (bernoulli 0.5) The robot starts from the position x0 and receives at each step a command Right or Left. It then moves to the left or right following the command with a probability of 10% of remaining in place (modeled by a Bernoulli distribution of parameter 0.1). let proba move (x0, cmd) = x where rec slip = sample (bernoulli 0.1) and xp = x0 -> pre x and x = match cmd with | Right -> min max_pos (if slip then xp else xp + 1) | Left -> max min_pos (if slip then xp else xp -1) end

  let proba slam (obs, cmd) = (map, x) where rec init map = Array_misc.ini (max_pos + 1) bernoulli_priors () and x = move (0, cmd) and o = Array_misc.get map x and p = if o then (1. -. sensor_noise) else sensor_noise and () = observe (bernoulli p, obs)

  let proba new_track_init_fn _ = (new_track_num (), sample (mv_gaussian (mu_new, sigma_new)))

  let proba clutter_init_fn _ = (mv_gaussian (mu_clutter, sigma_clutter))

  let proba obsfn (var, value) = observe (var, value) let proba model inp = t where rec init t = [] and t_survived = filter death_fn (last t) and n_new = sample (poisson lambda_new) and t_new = ini new_track_init_fn n_new and t_tot = append t_survived t_new and t = map state_update_fn t_tot and obs = map observe_fn t and n_clutter = (length inp) -(length obs) and () = observe(poisson lambda_clutter, n_clutter) and clutter = ini clutter_init_fn n_clutter and obs_shuffled = sample (shuffle (append obs clutter)) and present (not (n_clutter < 0)) -> do () = (iter2 obsfn (obs_shuffled, inp)) done

Figure 21 .

 21 Figure 21. Accuracy as a function of the number of particles.

Figure 22 .Figure 23 .

 2223 Figure 22. Runtime performance as a function of particles.

Figure 24 .Figure 25 .

 2425 Figure 24. Runtime performance as a function of particles.

Figure 26 .

 26 Figure 26. Memory consumption at each step of a run.

Figure 27 .

 27 Figure 27. Runtime performance at each step of a run.

Figure 28 .

 28 Figure 28. Memory consumption at each step of a run.

  The position p is only sporadically observed. process that computes the stream u of commands. It is written as two mutually recursive equations that define x_dist from u and u from the previous value of x_dist.

	. . .	x t -1	x t	x t +1		x t +2	. . .
		a t -1	a t	a t +1	p t +1	a t +2

Figure 2. Kalman filter for the robot example. Variables are either latent (white, e.g., state x) or observed (gray, e.g., acceleration a).

Table 1 .

 1 Benchmarks with: inference of fixed parameters from observations, estimation of a moving state (state-space model), and inference-in-the-loop (IITL).

		Fixed Moving IITL Metric
	Beta-Bernoulli	MSE
	Gaussian-Gaussian	MSE
	Kalman-1D	MSE
	Outlier	MSE
	Robot	LQR
	SLAM	MSE
	MTT	MOTA*

Language distribution and manual available at http://zelus.di.ens.fr.

The if construct can thus be considered as an external operator.

d ::= let f = e | d d e ::= c | x | (e, e) | op(e) | e(e) | if e then e else e | let p = e in e | fun p -> e | sample(e) | observe(e, e) | infer((fun x -> e), e) p ::= x | (p, p)

C(let node f x = e) = let f _step = fun (s,x) -> C(e)(s) let f _init = A(e)4.3 Semantics EquivalenceWe showed how to compile ProbZelus to µF a simple functional language with no loops, no recursion, and no higherorder functions, extended with the probabilistic operators. This language is similar to the kernel presented in[START_REF] Staton | Commutative Semantics for Probabilistic Programming[END_REF] for which a measure-based probabilistic semantics is defined (see also Figure16of Appendix B.2).

Acknowledgments

This work was supported in part by the MIT-IBM Watson AI Lab and the Office of Naval Research (ONR N00014-17-1-2699).