Ugo Dal Lago

Francesco Gavazzo

Ryo Tanaka

Effectful applicative similarity for call-by-name lambda calculi

Keywords: Call-by-name λ-calculus, applicative similarity, Howe's method, Algebraic effects

published or not. The documents may come

Effectful applicative similarity for call-by-name lambda calculi

Introduction

Program equivalence and refinement are crucial notions in the theory of programming languages, and their study is at the heart of programming language semantics. When higher-order functions are available, programs have a non-trivial interactive behaviour, and giving a satisfactory and at the same time handy definition of program equivalence becomes complicated. The problem has been approached, in many different ways, e.g. by denotational semantics or by contextual equivalence. These two approaches have their drawbacks, the first one relying on the existence of a denotational model, the latter quantifying over all contexts, thus making the task of proving programs equivalent quite hard. Handier methodologies for proving programs equivalent have been introduced along the years based on logical relations and applicative bisimilarity. Logical relations were originally devised for typed, normalising languages, but later generalised to more expressive formalisms, e.g., through step-indexing [START_REF] Appel | An indexed model of recursive types for foundational proof-carrying code[END_REF] and biorthogonality [START_REF] Benton | Relational semantics for effect-based program transformations: higher-order store[END_REF]. Starting from Abramsky's pioneering work on applicative bisimilarity [START_REF] Abramsky | The lazy lambda calculus[END_REF], coinduction has also been proved to be a useful methodology for program equivalence, and has been applied to a variety of calculi and language features [START_REF] Lassen | Relational reasoning about functions and nondeterminism[END_REF][START_REF] Dal Lago | On coinductive equivalences for higher-order probabilistic functional programs[END_REF][START_REF] Ong | Non-determinism in a functional setting[END_REF].

The just described scenario also holds when the underlying calculus is not pure, but effectful. There have been many attempts to study effectful λ-calculi [START_REF] Plotkin | Adequacy for algebraic effects[END_REF][START_REF] Moggi | Computational lambda-calculus and monads[END_REF]) by way of denotational semantics [START_REF] De'liguoro | Non deterministic extensions of untyped lambda-calculus[END_REF][START_REF] Jones | Probabilistic non-determinism[END_REF], logical relations [START_REF] Bizjak | Step-indexed logical relations for probability[END_REF], and applicative bisimilarity [START_REF] Lassen | Relational reasoning about functions and nondeterminism[END_REF][START_REF] Dal Lago | On coinductive equivalences for higher-order probabilistic functional programs[END_REF][START_REF] Crubillé | On probabilistic applicative bisimulation and call-by-value λ-calculi[END_REF]. But while the denotational and logical relation semantics of effectful calculi have been studied in the abstract (Goubault-Larrecq et al., 2008;[START_REF] Johann | A generic operational metatheory for algebraic effects[END_REF], the same cannot be said about applicative bisimilarity and related coinductive techniques. There is a growing body of literature on applicative bisimilarity for calculi with, e.g., nondeterministic [START_REF] Lassen | Relational reasoning about functions and nondeterminism[END_REF], and probabilistic effects [START_REF] Dal Lago | On coinductive equivalences for higher-order probabilistic functional programs[END_REF], but each notion of an effect has been studied independently, often getting different results. Distinct proofs of congruence for applicative bisimilarity, even if done through a common methodology, namely the so-called Howe's method [START_REF] Howe | Proving congruence of bisimulation in functional programming languages[END_REF], do not at all have the same difficulty in each of the cases cited above.

The observations above naturally lead to some questions. Is there any way to factor out the common part of the congruence proof for applicative bisimilarity in effectful calculi? Where do the limits on the correctness of applicative bisimilarity lie, in presence of effects?

This paper is part of a longstanding research effort directed to giving answers to the questions above. The first two authors, together with Paul Blain Levy, have recently introduced a general notion of applicative bisimilarity for lambda-calculi based on monads and relators [START_REF] Dal Lago | Effectful applicative bisimilarity: Monads, relators, and Howe's method (long version)[END_REF]. This provides, under mild conditions, a sound methodology for checking contextual equivalence of programs. The central idea is to take simulations as relations on terms and to use a relator to lift a (candidate) simulation to a relation on the monadic images of values. There is however little hope to prove a generic full-abstraction result in such a setting, although for certain notions of an effect, full abstraction is already known to hold [START_REF] Abramsky | The lazy lambda calculus[END_REF][START_REF] Crubillé | On probabilistic applicative bisimulation and call-by-value λ-calculi[END_REF].

In this paper, we study a different notion of simulation, which puts in relation not terms but their semantics. This way, interaction between terms and the environment can be modeled by an ordinary, deterministic, transition system and, with minimal side conditions, similarity can be proved to be not only a sound, but also to coincide with the contextual preorder. This, however, only holds when terms are call-by-name evaluated.

The paper is organised as follows. We first informally introduce the main ideas behind the present work by means of an extended example. We will then enter the technical developments of the work by first introducing a computational call-by-name λ-calculus Λ Σ enriched with arbitrary algebraic effects à la Plotkin and Power [START_REF] Plotkin | Adequacy for algebraic effects[END_REF]. The latter are defined by means of a signature of effect-triggering operation symbols which are interpreted as algebraic operations with respect to a monad.

In Section 4 we define an abstract notion of observation on top of which we define contextual and CIU approximation, as well as applicative similarity. Our first main result (Theorem 2), namely the equivalence between applicative similarity and CIU approximation, is proved in Section 5. Section 6, instead, is dedicated to our second main result (Corollary 1), namely full abstraction of applicative similarity (and thus CIU approximation) for contextual approximation. This result relies on a suitable CIU theorem (Theorem 3) whose proof is based on a variation of the so-called Howe's method [START_REF] Howe | Proving congruence of bisimulation in functional programming languages[END_REF].

Lifting Transition Systems

In this section we introduce the main ideas behind our approach by means of an example, namely a call-by-name nondeterministic λ-calculus [START_REF] Ong | Non-determinism in a functional setting[END_REF][START_REF] Lassen | Relational reasoning about functions and nondeterminism[END_REF], which we call Λ . The syntax of Λ extends the syntax of the ordinary λ-calculus [START_REF] Barendregt | The lambda calculus: its syntax and semantics[END_REF] with a binary nondeterministic choice operator , the intended semantics of an expression e f being to nondeterministically evaluate either to e or to f . We denote the collection of expressions and values (i.e. expressions of the form λx.e) by Λ and V, respectively, and let letters e, f, . . . and v, w, . . . range over elements of Λ and V, respectively. As usual, a program is a closed expression, i.e. an expression without free variables. We denote by Λ 0 the collection of all programs.

As we are in a nondeterministic setting we assume we can observe whether a program may converge 1 , meaning that the observable (operational) behaviour of a program e is whether there exists a value v such that e evaluates to v. For instance, the program (λx.x) Ω, where Ω is the purely divergent expression (λx.xx)(λx.xx), may converge to λx.x. We write e ⇓ v if e may converge to the value v, and e⇓ if there exists a value v such that e ⇓ v.

The notion of contextual approximation captures the idea of observational refinement in any possible environment. In our setting, as we have already remarked, observations are modelled by the predicate ⇓, whereas the role of the environment is played by contexts2 . Accordingly, we say that a program e contextually approximates a program f , written e ctx f , if for any context C, C[e]⇓ implies C[f ]⇓. We refer to the relation ctx as contextual approximation.

Although intuitively appealing, the relation ctx comes with a major drawback, namely its universal quantification over all possible contexts. Proving that a program e contextually approximates a program f requires to test the behaviour of e and f in any possible environment. For this reason, several proof techniques for ctx have been proposed in the last decades. Such techniques often come in the form of a relation R between programs such that R ⊆ ctx , meaning that e R f implies e ctx f . We say that R is sound for ctx if R ⊆ ctx , and that it is fully abstract

for ctx if R = ctx .
Starting with the seminal work by Abramsky [START_REF] Abramsky | The lazy lambda calculus[END_REF], applicative similarity has been proved to be a powerful proof technique for contextual approximation, both for the pure λ-calculus [START_REF] Abramsky | The lazy lambda calculus[END_REF] and for its nondeterministic and probabilistic extensions [START_REF] Ong | Non-determinism in a functional setting[END_REF][START_REF] Lassen | Relational reasoning about functions and nondeterminism[END_REF][START_REF] Dal Lago | On coinductive equivalences for higher-order probabilistic functional programs[END_REF]. The notion of applicative similarity for Λ is rooted in the observation that the operational semantics of Λ defines a (labelled) transition system T = Λ 0 , L, → over programs. Accordingly, the set of states3 is defined as the collection of programs, whereas the set of labels (or actions) L is defined as Λ 0 ∪ {τ }, where τ is a special symbol denoting evaluation. Finally, the transition relation → is defined as follows:

1. e τ - → v if e ⇓ v. 2. λx.f e - → f [e/x],
for any e ∈ Λ 0 .

Notice that due to the presence of nondeterminism, the above transition system is clearly nondeterministic. In general, a term e may converge to more than one value.

The notion of an applicative simulation is defined as ordinary simulation on T. Spelling out the definition, it is not hard to see that a relation R ⊆ Λ 0 × Λ 0 is an applicative simulation if and only if e R f implies that for any value λx.g such that e ⇓ λx.g, there exists a value λx.h such that f ⇓ λx.h and for any term k we have g

[k/x] R h[k/x].
Applicative similarity app is then coinductively defined as the largest applicative simulation. Notice that applicative similarity has better mathematical properties than ctx . In fact, app does not universally quantify over all contexts, but only on programs passed as arguments to lambda abstractions. Furthermore, app comes with an associated coinduction proof principle: in order to prove that a program f is applicatively similar to a program e, it is sufficient to exhibit an applicative simulation R relating e and f .

Remarkably, applicative similarity is included in contextual approximation, meaning that app is a sound proof technique for ctx . However, it is not hard to realise that full abstraction fails. In fact, we can easily convince ourselves that the program e = λx.(f g) contextually approximates h = λx.f λx.g, although e app h.

An informal argument supporting e ctx h is the following (a formal proof will be given later in full generality). As Λ is a call-by-name calculus, a context C can only test a program k in function position. That is, during the evaluation of a program C[k], k can only be evaluated in function position, i.e. in a situation of the form kc, where c may 'contain' other occurrences of k. Similarly, the result of the evaluation of kc can only be tested in function position. Iterating this observation, we see that testing a program k in a context C is equivalent to testing k against a finite sequence of inputs. This is the content of the so-called CIU Theorem [START_REF] Mason | Equivalence in functional languages with effects[END_REF]:

Theorem 1 (CIU, Λ -version). For all programs e, f , e ctx f holds if and only if for any finite sequence of programs k

1 , . . . , k n , ek 1 • • • k n ⇓ implies f k 1 • • • k n ⇓.
Notice that the CIU Theorem states that ctx is a form of trace refinement. Using Theorem 1 we can easily show e ctx h. This comes with no surprise: contextual approximation being a form of trace refinement, it is not sensitive to 'branching behaviours', and thus by no means can tell e and h apart. This is not the case for applicatively similarity which, instead, is sensitive to forms of branching. From an applicative simulation perspective, the difference between the expressions e and h as defined above is clear: e postpones the nondeterministic choice to the moment it receives an input, whereas h first makes the choice, and then waits to receive an input. Accordingly, it is sufficient to instantiate f as xΩ(λx.x) and g as x(λx.x)Ω to observe that e app h. For suppose λx.(f g) λx.f λx.g. In particular, there is a simulation R relating the two programs. By the very definition of applicative simulation, since λx.(f g) evaluates to itself and λx.f λx.g evaluates to either λx.f or λx.g, we must have λx.(f g) R λx.f or λx.(f g) R λx.g. In the first case, since R is a simulation, we must have

(f [k/x] g[k/x])Rf [k/x],
for any program k. This, however, cannot be the case, as taking k = λxy.x, we see that f [k/x] g[k/x] may converge to λx.x, whereas f [k/x] diverges. A similar argument shows that λx.(f g) R λx.g cannot hold as well (take k = λxy.y).

The reader might have recognised that e and h are nothing more than the encoding in Λ of the labelled transition systems below, which are the standard example in concurrency theory showing that bisimilarity is strictly finer than trace equivalence.

• a • b z z b $ $ • • • a z z a $ $ • b • b • •
Summing up, we have observed that in a nondeterministic setting, applicative similarity is sensitive to branching, whereas contextual approximation is not. This is because, in full generality, contextual approximation is a form of trace approximation, and thus intrinsically deterministic. This last observation suggests an easy way to make applicative similarity fully abstract for contextual approximation, namely to determinise the transition system T, and thus the very notion of applicative similarity. Such a construction is well-known, and dates back at least to the powerset construction in automata theory [START_REF] Rabin | Finite automata and their decision problems[END_REF]. Roughly speaking, the idea, which we are going to analyse in a more detailed fashion, is to make T deterministic by lifting the transition relation to sets of programs.

An Abstract View on Applicative Similarity

As we have already observed, T is a nondeterministic transition system. In fact, the transition relation τ -→⊆ Λ 0 × V 0 can be equivalently described as the evaluation function ev : Λ 0 → PV 0 mapping a program e to the set {v ∈ V 0 | e ⇓ v} (notice that ev(e) might be empty, finite, or infinite). Similarly, transitions of the form e -→ determine a function ap :

V 0 → Λ Λ 0 0 defined by ap(λx.f, e) = f [e/x].
As a consequence, we can equivalently describe T as the pair4 (ev : Λ 0 → PV 0 , ap : V 0 → Λ Λ 0 0 ). Applicative similarity can be then defined relying on a suitable notion of relation lifting, following the abstract approach of [START_REF] Dal Lago | Effectful applicative bisimilarity: Monads, relators, and Howe's method (long version)[END_REF]. A relation lifting is a map associating to any relation R ⊆ X × Y a relation ΓR ⊆ P(X) × P(Y ). Fixed such a map Γ, an applicative simulation can be defined as a relation R ⊆ Λ 0 × Λ 0 such that e R f implies ev(e) ΓR ev(f ), and λx.f R λx.g implies ap(λx.f, k) R ap(λx.g, k), for any program k. As a consequence, in order to recover the previous notion of applicative simulation for Λ , it is sufficient to instantiate Γ as follows:

V ΓR W ⇐⇒ ∀v ∈ V. ∃w ∈ W. v R w.
In this paper we consider a different approach to applicative similarity, which we could summarise in the slogan 'lift the transition system, not the relation'. In fact, instead of lifting relations between programs to relations between sets of programs, we can lift the whole transition system T to a transition system T * such that (i) states of T * are sets of programs, and (ii) the transition relation (and thus the associated notion of simulation) of T * is deterministic. We can build T * by lifting the functions ev and ap to:

ev * : P(Λ 0 ) → P(V 0 ) ap * : P(V 0 ) → P(V 0 ) Λ 0 .
The latter functions are defined in the obvious way, relying on the (strong) monad (MacLane, 1971) structure of P. Furthermore, as we are interested in observing whether a program may convergence, we equip T * with an observation function ob : P(V 0 ) → {⊥, }, where {⊥, } denotes the two-element Boolean algebra, defined by ob(V ) = if and only if V = ∅. Notice that we have ob(ev(e)) = if and only if e ⇓.

Finally, we define a notion of applicative simulation between elements in P(Λ 0 ) as a relation R ⊆ P(Λ 0 ) × P(Λ 0 ) such that:

1. E R F implies ob(ev * (E)) ≤ ob(ev * (F )); 2. V R W and V, W ∈ P(V 0 ) imply ap * (V, e) R ap * (W, e), for any program e.
As usual, applicative similarity ( app ) * is defined as the largest applicative simulation.

We immediately see that although coinductively defined, this new notion of applicative similarity is a trace-like refinement. For instance, it is straightforward to see that5 λx.(f g) ( app ) * λx.f λx.g. Additionally, contrary to app , ( app ) * is not only sound, but also fully abstract for ctx , as we will prove in full generality in Section 6.

At this point it is natural to ask whether the above construction can be given for other effectful languages, such as calculi with primitives for probabilistic nondeterminism or input/output. In fact, the reader may have noticed that the central ingredients used to build T * (at least from a mathematical perspective) were the (strong) monad structure of the powerset functor P and a suitable observation function ob. The rest of this paper answers the above question in the affirmative by generalising the above construction to a call-by-name λ-calculus enriched with algebraic effects à la Plotkin and Power [START_REF] Plotkin | Adequacy for algebraic effects[END_REF], this way providing a sound and complete characterisation of contextual approximation for a large class of effectful callby-name calculi.

Mathematical Preliminaries

We assume the reader to be familiar with basic domain theory and category theory. In particular, we give for granted the notions of an ω-complete partial order (ω-cpo, hereafter) and of a pointed ω-cpo (ω-cppo, hereafter), as well as the notions of a monotone, continuous, and strict function. We also assume the reader to be familiar with the basic vocabulary of category theory, viz. the notions of a category, functor, and monad. Since we will work with the category Set of sets and functions only, we simply refer to functors and monads in place of functors and monads on Set, respectively. The reader can consult [START_REF] Abramsky | Domain theory[END_REF][START_REF] Davey | Introduction to lattices and order[END_REF][START_REF] Maclane | Categories for the Working Mathematician[END_REF] for further details.

Following (Plotkin andPower, 2001, 2003), we consider (algebraic) operations as sources of side effects. Syntactically, such operations are described by means of a signature of operation symbols, i.e. a pair Σ = (F, α) consisting of a set F of operation symbols and a map α : F → N, assigning to each operation symbol a (finite) arity. Semantically, operation symbols are interpreted as algebraic operations on a monad/Kleisli triple [START_REF] Maclane | Categories for the Working Mathematician[END_REF][START_REF] Moggi | Computational lambda-calculus and monads[END_REF] 6 . Definition 1. A monad is a triple T = T, η, >>= where T is a map -called the carrier of T -associating to each set X a set T X, η -called the unit of T -is a family of functions η X : X → T X, and >>= -called the bind of T -is an operation mapping each function f : X → T Y to the function >>=f : T X → T Y subject to the following identities:

(>>=f ) • η = f >>= η = id >>=(>>=g • f ) = >>=g • >>=f,
for f and g functions with appropriate domain and codomain.

As it is customary, we write t >>= f in place of >>=f (t), for f : X → T Y and t ∈ T X. Moreover, when clear from the context we will omit subscripts, thus e.g. writing η(x) in place of η X (x).

Definition 2. Let Σ be a signature and T = T, η, >>= be a monad. We say that Σ is interpreted on T if associated to any n-ary operation symbol σ in Σ is a set-indexed family of functions σ X : (T X) n → T X. We say that Σ is algebraic for T if the following identity holds, for any n-ary operation symbol σ ∈ Σ, map f : X → T Y , and all elements t 1 , . . . , t n ∈ T X:

σ X (t 1 , . . . , t n ) >>= f = σ Y (t 1 >>= f, . . . , t n >>= f ).
Notice that algebraicity of Σ for T is essentially requiring each T X to carry a Σ-algebra structure such that >>=f is a Σ-algebra homomorphism, for any set X and function f : X → T Y .

Example 1. The following are examples of monads and algebraic operations on them. Due to space constraints, we are forced to omit many interesting examples, such as the exception and global state monad, for which we refer to [START_REF] Dal Lago | Effectful applicative bisimilarity: Monads, relators, and Howe's method (long version)[END_REF][START_REF] Plotkin | Algebraic operations and generic effects[END_REF].

1. The maybe monad M has carrier X ⊥ = {just x | x ∈ X}∪{⊥}, whereas unit and bind are thus defined:

η(x) = just x (just x) >>= f = f (x) ⊥ >>= f = ⊥.
The maybe monad is used to model partial computations: the symbol ⊥ denotes divergence, whereas just x denotes the result of a computation giving result x. Since the possibility of divergence is an intrinsic feature of any (Turing complete) programming language, no operation is required to produce it.

The output monad

O has carrier O(X) = A ∞ × X ⊥ ,
where A is a given alphabet and A ∞ denotes the set of finite and infinite strings over A.

Unit and bind of O are thus defined:

η(x) = (ε, just x) (u, ⊥) >>= f = (u, ⊥) (u, just x) >>= f = (uw, a),
where (w, a) = f (x), ε denotes the empty string, and uw denotes the concatenation of u and w (as usual, if u is infinite we define uw as u). The output monad models partial computations with output. An element (u, a) represents the result of a computation printing the string u during its processing. Due to the possibility of divergence, the string u may be infinite (and in such a case we expect a to be ⊥). We consider a signature containing a unary A-indexed operation symbol print c , where c ∈ A, whose intended semantics is to print the character c and then continue with its argument. Formally, we interpret print c ((u, x)) as (cu, x).

3. The non-empty powerset monad P has carrier P ne (X) = {U ⊆ X | U = ∅}, whereas unit and bind are thus defined:

η(x) = {x} U >>= f = x∈U f (x).
The non-empty powerset monad is used to model total nondeterministic computations. A set U ⊆ X denotes the possible results of a computation, meaning that if x ∈ U then x is one of the possible, nondeterministic outcomes of the computation. Nondeterministic partial computations are modelled by post-composing P with M. It is a routine exercise to verify that the latter is indeed a monad. We consider a signature containing a binary operation symbol for pure nondeterministic choice. We obtain the desired behaviour by interpreting as set-theoretic union.

4. The distribution monad D has carrier

D = {µ ∈ [0, 1] X | x µ(x) = 1},
where the map µ is assumed to have countable support7 . Unit and bind of D are thus defined:

η(x)(x ) = 1 if x = x 0 otherwise (µ >>= f )(y) = x∈X µ(x) • f (x)(y).
Since the unit of D is given by the so-called Dirac distribution, we will denote it by δ. The distribution monad is used to model total probabilistic computations, with the intended meaning that a distribution µ gives the probability of convergence to values. Probabilistic partial computations are modelled by post-composing D with M. It is a routine exercise to verify that the latter is indeed a monad, and actually isomorphic to the subdistribution monad. We can consider a signature containing a binary operation symbol or for fair probabilistic choice.

We obtain the desired behaviour by interpreting or as the operation ⊕ defined by (µ ⊕ ν)(x) 1 2 µ(x) + 1 2 ν(x). Since we will define the evaluation of a program e as the limit of its finite evaluations, we require monads to come with a suitable domain-theoretic structure.

Definition 3. Given a monad T = T, η, >>= and a signature Σ interpreted on T, we say that T is Σ-continuous if for any set X, the set T X carries an ω-cppo structure (with order X ) such that >>= is strict in its first argument (i.e. ⊥ >>= f = ⊥), and:

1. For any k-ary operation symbol σ in Σ σ X ( n<ω t 1 n , . . . , n<ω t k n ) = n<ω σ X (t 1 n , . . . , t k n ).
2. For any ω-chain (t n ) n<ω in T X and any ω-chain (f n ) n<ω in X → T Y (notice that the latter inherits a ω-cppo structure from T Y pointwise)

( n<ω t n ) >>= ( n<ω f n ) = n<ω (t n >>= f n ).
Example 2. All monads in Example 1 are Σ-continuous. In particular, X ⊥ is an ω-cppo with respect to the flat order x y ⇐⇒ x = ⊥ or x = y. To see that PM is Σ-continuous consider the ordering U V ⇐⇒ just -1 (U ) ⊆ just -1 (V ). Notice that we have {⊥} V for any set V . Finally, for DM define the order µ ν ⇐⇒ ∀x. µ(just x) ≤ ν(just x). Notice that we have η(⊥) µ, for any µ ∈ D(X ⊥ ). Checking continuity of >>= and of operations is a routine exercise. Ordering the output monad O requires a bit more attention, due to non-monotonicity of string concatenation with respect to the prefix order ⊆ on A ∞ (see [START_REF] Dal Lago | Effectful applicative bisimilarity: Monads, relators, and Howe's method (long version)[END_REF]). We define (u, a) (w, b) if and only if a = ⊥ and u ⊆ w, or a = just x, b = just y, u = w, and x = y.

We can now introduce the effectful calculus Λ Σ .

Syntax and Operational Semantics

In this section we introduce the computational language Λ Σ , and give it call-by-name monadic operational semantics. In the rest of this section let T = T, η, >>= be a fixed Σ-continuous monad and Σ be a fixed algebraic signature on it.

Definition 4. The collections Λ • and V • of expressions and values of Λ Σ are defined by the following grammar(s), where σ ranges over n-ary operation symbols in Σ (and x ranges over a fixed set of variables). We adopt standard conventions as in [START_REF] Barendregt | The lambda calculus: its syntax and semantics[END_REF]. In particular, we denote by F V (e) the collection of free variables of an expression e and identify expressions up to renaming of bound variables. We say that an expression e is closed (resp. open) if F V (e) = ∅ (resp. F V (e) = ∅) and refer to closed expressions as programs. We denote 8 by Λ(x) the set of expressions with free variables from x and write Λ in place of Λ(∅) (similar conventions apply to the set of values). Finally, we write f [e/x] for the capture-free substitution of the expression e for all free occurrences of the variable x in the expression f . Call-by-name evaluation contexts are expressions with a single hole [-] defined by the grammar E ::= [-] | Ee. We say that E is closed if it has no free variables. We write E[e] for the expression obtained by substituting the expression e for the hole [-] in E. Redexes are programs of the form (λx.f )e or σ(e 1 , . . . , e n ), the former producing a computation step, the latter producing the effect described by the operation σ. We notice that any program is either a value or a program of the form E[r], for a redex r.

Operational semantics is defined by means of an evaluation function -: Λ → T V mapping each program e ∈ Λ to a monadic value e ∈ T V. For instance, a probabilistic program is evaluated to a subdistribution of values (modelled via DM).

Definition 5. Define the N-indexed family of mapsn : Λ → T V as follows:

e 0 = ⊥ v n+1 = η(v) E[(λx.f )e] n+1 = E[f [e/x]] n E[σ(e 1 , . . . , e k ) n+1 = σ V ( E[e 1 ] n , . . . , E[e k ] n )
The sequence ( e n ) n<ω forms an ω-chain in T V, so that we can define e as n<ω e n .

Example 3. Let T be DM, and let e = Y (λx.I ⊕ x), where I is the identity combinator and Y is Curry's fixed point combinator. Writing a distribution as a formal sum 9 , for any n ≥ 1 we have

e n = n i=1 1 2 n •I +(1-n i=1 1 2 n )•⊥,
and thus e = sup n e n = 1 • I + 0 • ⊥.

Before defining behavioural equivalences and refinements for Λ Σ , it is useful to spell out some properties of the evaluation map -. 9 Any discrete distribution µ ∈ D(X) can be written as an expression of the form

i<ω p i • x i , where µ(x) = xi=x p i .
Lemma 1. The following identities hold.

v = η(v) E[(λx.e)f ] = E[e[f /x]] E[σ(e 1 , . . . , e k ) = σ V ( E[e 1 ] , . . . , E[e k ] ).
Proof. Simply observe that since σ V and >>= are continuous then so is -.

Lemma 2. For any program e and evaluation context E, we have:

E[e] n e n >>= (v → E[v] n ) E[e] = e >>= (v → E[v] ).
Proof sketch. The first inclusion can be easily proved by induction on n.

For the second identity one shows by induction on n that for any n ≥ 0

e n >>= (v → E[v] ) E[e] holds. As a consequence, e >>= (v → E[v] ) E[e]
. Since from the first inclusion (and Lemma 1) follows

E[e] e >>= (v → E[v] ) we conclude E[e] = e >>= (v → E[v] ).

Program Refinement

In this section we define a general notion of observation on top of which we define several notions of program refinement, the latter being defined as relations between programs. In order to qualify as a notion of a program refinement, a relation R needs to satisfy some minimal desiderata. Clearly, R must be a preorder. Furthermore, if e R f , then the observable behaviour of f , which is given through our general notion of observation, must refine the observable behaviour of e. This property is known as preadequacy. Finally, R must respect the way programs are (syntactically) constructed, a property known as compatibility. Notice that relying on compatibility one can show the refinement C[e] R C[f ] compositionally, by proving e R f . Summing up, we require program refinement relations to be preadequate precongruence relations.

Compatibility, however, requires to extend the notion of a program refinement to arbitrary expressions of the language. Indeed, dealing with a λ-abstraction λx.e compositionally, requires to work with the open expression e. As a consequence, we work with relations on arbitrary expressionswhich we call open relations -rather than on programs only -which we refer to as closed relations (see [START_REF] Pitts | Howe's method for higher-order languages[END_REF][START_REF] Gordon | A tutorial on co-induction and functional programming[END_REF][START_REF] Lassen | Relational reasoning about functions and nondeterminism[END_REF]). Definition 6. An open relation is a set R of triples (x, e, f ), where e, f ∈ Λ(x), which is closed under weakening, meaning that (x, e, f ) ∈ R implies (x ∪ {x}, e, f ) ∈ R.

We use the infix notation and write x e R f in place of (x, e, f ) ∈ R and abbreviate ∅ e R f as e R f . As a convention, when referring to relations we tacitly mean open relations.

x e R • f ⇐⇒ ∀ k ∈ Λ. e[ k/x] R f [ k/x].
The notion of reflexivity, symmetry, and transitivity straightforwardly extends to open relations (see e.g. [START_REF] Pitts | Howe's method for higher-order languages[END_REF]).

Definition 8. Let R be an open relation. We say that R is:

1. Compatible if it satisfies the clauses in Figure 1.

Substitutive if:

x

∪ {x} e R f ∧ x g R h =⇒ x e[g/x] R f [h/x].
3. Closed under substitution if:

x ∪ {x} e R f ∧ k ∈ Λ(x) =⇒ x e[k/x] R f [k/x].
Notice that the open extension of a closed relation is always closed under substitution. If R is a compatible preorder, we say that R is a precongruence.

Having defined the notion of a precongruence relation, we now move to preadequacy. As already remarked, the latter builds on a suitable notion of observation. Intuitively, an observation is a function that takes an element in T V representing the result of an effectful computation, and returns what we can observe of such a computation. Following standard practice, we assume that the observable part of a computation consists of the side effects happened during such a computation. As a consequence, the observable part of a computation is uniquely determined by the semantics of the (algebraic) operations involved.

C-var

x x R x x ∪ {x} e R f C-abs x λx.e R λx.f x e R f x g R h C-app x eg R f h x e 1 R f 1 . . . x e n R f n C-op x σ(e 1 , . . . , e n ) R σ(f 1 , . . . , f n ) Figure 1: Compatibility clauses.
Definition 9. Let 1 = { * } be the one element set and ! X : X → 1 be the unique function mapping each element in X to * . Define the observation function ob : T V → T 1 as T (! V ). The map ob trivially extends to programs as ob(e) = ob( e ).

Proposition 1. For all monadic values ϕ 1 , . . . , ϕ n , ω-chain (ϕ n ) n<ω in T V, and value v, we have:

ob(⊥) = ⊥ T (1) ob(η V (v)) = η 1 ( * ) ob(σ V (ϕ 1 , . . . , ϕ n )) = σ 1 (ob(ϕ 1 ), . . . , ob(ϕ n )) ob( n<ω ϕ n ) = n<ω ob(ϕ n ).
Proof. We first observe that by the very definition of monad we have ob = >>=(η 1 • !). As a consequence, we can rely on strictness and continuity of >>= to prove the first and last identities above. The second identity follows by Definition 1, whereas the third one is a direct consequence of algebraicity of operation symbols.

We can now define an abstract notion of preadequacy.

Definition 10. An open relation R is preadequate (with respect to ob) if ∅ e R f implies ob(e) ob(f ).
Having at our disposal the notion of a precongruence preadequate relation, we can introduce a canonical precongruence preadequate program refinement, namely Morris' contextual approximation [START_REF] Morris | Lambda calculus models of programming languages[END_REF].

Definition 11. Recall that contexts are syntax trees defined by the grammar:

C ::= [-] | eC | Ce | λx.C.
As usual, we write C[e] for the result of replacing the hole [-] in C with the expression e. Given an expression e, we say that a context C closes e, if C[e] is a program. We denote the collection of such contexts as Cl(e). The open relation ctx , called contextual approximation, is thus defined:

x e ctx f ⇐⇒ ∀C ∈ Cl(e) ∩ Cl(f ). ob(C[e]) ob(C[f ]).
Remark 1. Notice that in Definition 11 contexts are defined as syntax trees (i.e. expressions not considered modulo renaming of bound variables), and not as expressions. That means that when substituting an expression e for the hole [-] in a context C, we cannot consider ordinary capture-avoiding substitution. For instance, for C = λx.[-] we require C[x] to be λx.x (i.e. to capture the free variable x), whereas our definition of substitution gives λy.x. In order to be formally precise, we should then give an explicit definition of C[e], as it is done in [START_REF] Pitts | Howe's method for higher-order languages[END_REF]. Another, elegant, possibility is to characterise ctx as the largest compatible preadequate relation. In fact, under mild conditions on (algebraic) operations, it is easy to prove that the union of all compatible preadequate relations is itself preadequate and compatible, and thus the largest such relation (see [START_REF] Dal Lago | Effectful applicative bisimilarity: Monads, relators, and Howe's method (long version)[END_REF]). We see that ob(e) gives the probability of convergence of e, meaning that e ctx f holds if and only if for any context C, the probability of converge of C[f ] is an upper bound of the probability of convergence of C[e].

As stated in Remark 1, contextual approximation enjoys several nice properties, notably it is the largest preadequate compatible relation. For instance, since the composition

C 1 [C 2 [-]] of two contexts C 1 , C 2 is itself a context, ctx
is obviously compatible. Furthermore, since is a preorder, then so is ctx .

Before introducing applicative similarity, we define CIU approximation [START_REF] Mason | Equivalence in functional languages with effects[END_REF], i.e. the restriction of contextual approximation to evaluation contexts.

Definition 12. The closed relation ciu , called CIU approximation is thus defined:

e ciu f ⇐⇒ ∀E. ob(E[e]) ob(E[f ]),
where E ranges over closed evaluation contexts. As usual, we can extend ciu to arbitrary expressions by taking its open extension.

It is immediate to see that since any evaluation context is of the form [-]e 1 • • • e n , for some programs e 1 , . . . , e n , the relation ciu is a trace-like refinement which tests the applicative behaviour of programs against finite traces of the form e 1 , . . . , e n .

Applicative Similarity

In this section we define effectful applicative similarity and prove its equivalence with ciu . As already discussed in the introduction, at the very heart of our notion of effectful applicative similarity is the shift from the ordinary transition system of programs to the lifted transition system of 'monadic programs'. In order to qualify such a lifted system as an applicative system, according to the terminology of [START_REF] Abramsky | The lazy lambda calculus[END_REF], we first need to lift the notion of an application from ordinary programs to monadic programs.

Definition 13. Let e be a program and ϕ ∈ T V be a monadic value. Define monadic application ϕ • e ∈ T V as ϕ >>= (v → ve ).

We immediately notice that ef = e • f . Moreover, straightforward calculations give the following identities:

⊥ • e = ⊥, η(λx.f ) • e = f [e/x] , σ V (ϕ 1 , . . . , ϕ n ) • e = σ V (ϕ 1 • e, . . . , ϕ n • e).
We can now define effectful applicative similarity.

Definition 14. A relation R ⊆ T V × T V is an effectful applicative simula- tion (hereafter applicative simulation) if ϕ R ψ implies: 1. ob(ϕ) ob(ψ). 2. ∀e ∈ Λ. ϕ • e R ψ • e.
Applicative similarity app is defined as the largest applicative simulation. We extend applicative similarity to programs as e app f ⇐⇒ e app f .

Notice that applicative similarity is well-defined since the defining clauses of applicative simulation in Definition 14 induce a monotone operator Φ on the complete lattice 2 T (V)×T (V) of binary relations on monadic values. As a consequence, we can define similarity as the greatest fixed point of Φ. Applicative similarity being defined coinductively, it comes with an associated coinduction principle:

ϕ R ψ ∧ (R applicative simulation ) =⇒ ϕ app ψ.
Remark 2. Definition 14 gives the notion of an applicative simulation directly on monadic values, without explicitly defining the lifted transition system of monadic values (the transition system T * of the introduction). However, it is easy to see that we can lift the transition system of programs (ev : Λ → T V, ap : V → Λ Λ ), where ev(e) = e and ap(λx.f, e) = f [e/x], to the system (ev * : T Λ → T V, ap * : T V → (T Λ) Λ ), where10 :

ev * (ξ)
ξ >>= ev ap * (ϕ, e) = ϕ >>= ap(_, e).

Proposition 2. Applicative similarity is reflexive and transitive.

Proof. The proof is by coinduction. We immediately notice that {(ϕ, ϕ) | ϕ ∈ T V} is an applicative simulation. Furthermore, we observe that if R and S are applicative simulations, then so is R; S. In fact, ϕ R ψ and ψ S ξ imply ob(ϕ) ob(ξ), since is transitive. Additionally, R and S being applicative simulations, we have ϕ • e R ψ • e and ψ • e S ξ • e, and thus ϕ • e (R; S) ξ • e, for any e ∈ Λ.

We now prove the first of our main results, namely that app is a trace-like refinement, which actually coincides with ciu . Such a result gives a coinductive characterisation of ciu , which allows one to prove program refinement relying on both the induction and coinduction proof principle.

Theorem 2. app = ciu and ( app ) • = ( ciu )

) • = ( ciu ) • .
Example 5. Recall that according to standard applicative similarity a program of the form (λx.f ) (λx.f ) does not simulate λx.f g. We now prove that this is not the case for app as defined in Definition 14. Actually, we can even strengthen such a result showing that for any operation symbol such that η( * ) σ 1 (η( * ), . . . , η( * )) (lambda) abstraction distributes over σ, i.e. λx.σ(e 1 , . . . , e n ) app σ(λx.e 1 , . . . , λx.e n ). To see that it is sufficient to observe that the relation Finally, a similar argument shows that if σ 1 (η( * ), . . . , η( * )) η( * ), then σ(λx.e 1 , . . . , λx.e n ) app λx.σ(e 1 , . . . , e n ). Such a condition holds both for the powerset and the distribution monads with the operations in Example 1, but fails for the output monad with the print operation.

Up to this point we have introduced three notions of behavioural refinement: contextual approximation ctx , CIU approximation ciu , and applicative similarity app . The former has been shown to be a precongruence, i.e. a compatible preorder, whereas the latter two to coincide and to be preorders. In the next section we show that (the open extensions of) ciu and app are precongruences too, and that ctx , ( ciu ) • , and ( app ) • coincide.

A CIU Theorem

In this section we prove a CIU theorem (Theorem 3) stating that (the open extension of) ciu coincides with ctx . In virtue of Theorem 2 we will also obtain the equivalence between (the open extension of) app and ctx , meaning that the former is fully abstract with respect to the latter. Our proof of Theorem 3 follows [START_REF] Pitts | Howe's method for higher-order languages[END_REF] and it is based on a variation of Howe's technique [START_REF] Howe | Proving congruence of bisimulation in functional programming languages[END_REF], the standard relational construction to prove compatibility of applicative similarity. For readability, in the rest of this section we write in place of ciu .

The main idea behind Howe's technique is to extend to an open relation H that is compatible and substitutive by construction. As a consequence, in order to prove that • is compatible, it is sufficient to prove • = H . Proving the inclusion • ⊆ H is rather straightforward, whereas the other inclusion requires more effort. When unambiguous, we will write R H in place of (R • ) H for a closed relation R. The following result states some useful properties of Howe extension. The proof is standard and can be found in e.g. [START_REF] Dal Lago | On coinductive equivalences for higher-order probabilistic functional programs[END_REF]. E of R to (closed) evaluation contexts is defined in Figure 3.

(R • ) H . x x R e H-var x x R H e x ∪ {x} e R H g x λx.g R f x ∈ x H-abs x λx.e R H f x e R H k x g R H h x kh R f H-app x eg R H f x e i R H f i x σ(f 1 , . . . , f n ) R g H-op x σ(e 1 , . . . , e n ) R H g
E-nil [-] R H E [-] ∅ e R H f E R H E F E-app Ee R H E F f
E R H E F x e R H f x E[e] R H F [f ] Proof. A straightforward induction on the derivation of E R H E F . Lemma 5. If E[e] H f
∀n < ω. ob( E[e] n ) ob( F [f ] ).
We proceed by induction on n. The case for n = 0 is trivial. Let n > 0 and suppose the thesis holds for any m < n. We proceed by case analysis on e.

Case 1. Suppose e is a value λx.k. Since e H f , the latter must be the conclusion of an inference of the form

{x} k H h λx.h f H-abs λx.k H f
We now proceed by case analysis on E. f . In particular, (λx.k)g H h must be the conclusion of an inference of the form:

If E ≡ [-], then F ≡ [-] too, since E H E F .
{x} k H c λx.c d H-abs λx.k H d g H b db h H-app (λx.k)g H h
As a consequence, by substitutivity we obtain k ]] we can apply the induction hypothesis, obtaining:

[g/x] H c[b/x]. Since E H E F and G H E H imply E[G[-]] H E F [H[-
ob( E[G[(λx.k)g]] n ) = ob( E[G[k[g/x]]] n-1 ) ob( F [H[c[b/x]]] ) = ob( F [H[(λx.c)b]] ) ob( F [H[db]] ) ob( F [H[h]] ) ob( F [f ] ).
In the last inequality we have have used the following implications:

λx.c d =⇒ (λx.c)b db =⇒ F [H[(λx.c)b]] F [H[db]]
and similarity for db h and H[h] f . Theorem 3 (CIU).

) • = ctx = ( app ) • .
We conclude this section observing that applicative similarity, contextual approximation, and CIU approximation being trace-like refinement, we can straightforwardly extend Corollary 1 to their symmetric counterparts.

Corollary 2. Define applicative bisimilarity app , contextual equivalence ctx , and CIU equivalence ciu as app ∩ ( app ) -1 , ctx ∩ ( ctx ) -1 , and ciu ∩ ( ciu ) -1 , respectively, where R -1 denotes the converse of R. Then, ( ciu ) • = ctx = ( app ) • . 7. Call-by-name, Call-by-value, and Call-by-need Corollaries 1 and 2 give coinductive characterisations of contextual approximation and equivalence, respectively. Such characterisations are proved through a CIU theorem, which takes advantage of the shape of call-by-name evaluation contexts. It is then natural to ask whether our results can be extended to encompass call-by-value calculi, the latter being more natural in presence of effects.

An effectful notion of applicative bisimilarity for a call-by-value λ-calculus with algebraic effects has been given in [START_REF] Dal Lago | Effectful applicative bisimilarity: Monads, relators, and Howe's method (long version)[END_REF], where an abstract soundness theorem for it is proved. Such a notion, however, is not fully abstract, as already observed in [START_REF] Lassen | Relational reasoning about functions and nondeterminism[END_REF] in the context of nondeterministic calculi. It is straightforward to see that the notion of applicative (bi)similarity proposed in the present paper can be easily given for a call-byvalue calculus. Unfortunately, the resulting notion would not only invalidate Corollary 1 and Corollary 2, but would also be unsound, as witnessed by the probabilistic programs λx.(x or Ω) and (λx.x) or (λx.Ω), and the context (λx.x(x(λy.y))) [-].

Another issue that we have not addressed in this work is whether our techniques can be extended to call-by-need calculi. Concerning program equivalence and refinement, call-by-name and call-by-need calculi are equivalent in absence of computational effects (different from nontermination). Such an equivalence breaks in presence of effects such as pure and probabilistic nondeterminism, global states, and input-output. For instance, the program (λx.x or x)((λy.y) or Ω) converges with probability 0.5 if evaluated according to a call-by-need strategy, and with probability 0.25 if evaluated according to a call-by-name strategy.

Although our techniques seem to scale to specific effectful call-by-need languages (such as languages with primitives for raising exceptions), the study of coinductive techniques for call-by-need languages with algebraic languages is nontrivial, and we leave it for future work.

Conclusion and Future Work

We have shown how a notion of applicative similarity on call-by-name effectful lambda-calculi can be defined and proved fully-abstract. This is very much in line with logical relations as introduced in [START_REF] Johann | A generic operational metatheory for algebraic effects[END_REF], but simpler and applicable to untyped λ-calculi. The underlying transition system has monadic values as states, and is essentially deterministic. This is indeed the reason the framework is only applicable to call-by-name (and not to call-by-value) calculi, in contrast with [START_REF] Dal Lago | Effectful applicative bisimilarity: Monads, relators, and Howe's method (long version)[END_REF].

As a future work, the authors would like to investigate other transition system-based techniques to obtain sound (and complete) notions of (bi)simulations for effectful calculi. A notable result in such a direction is [START_REF] Gianantonio | Rpo, second-order contexts, and lambda-calculus[END_REF].

v

  ::= x | λx.e e ::= v | ee | σ(e, . . . , e).

  There is a canonical way to extend a closed relation to an open one. Definition 7. The open extension of a closed relation R is the open relation R • defined as follows:

  Example 4. 1. For the maybe monad we have ob : V ⊥ → 1 ⊥ satisfying ob(just v) = just * (and thus ob(v) = just * ) and ob(⊥) = ⊥. Therefore, ob(e) = just * if and only e converges (notation e ⇓). As a consequence, e ctx f if and only if ∀C. C[e] ⇓ =⇒ C[f ] ⇓. 2. For the partial nondeterministic monad PM we have P ne (1 ⊥ ) = {{⊥}, {just * }, {⊥, just * }} so that ob(⊥) = {⊥}, ob({just v}) = {just * } (and thus ob(v) = {just * }) and ob(e f ) = ob(e)∪ob(f ). We thus have ob(e) = {⊥} if and only if e may converge, and ⊥ ∈ ob(e) if and only if e must converge. As a consequence, according to the order defined in Example 2, we have ob(e) ob(f ) if and only e ⇓ =⇒ f ⇓. 3. For the partial distribution monad DM of Example 1 we first observe that D(1 ⊥ ) ∼ = [0, 1]. It is then easy to see that we obtain ob(⊥) = 0, ob(δ(just v)) = 1 (and thus ob(v) = 1) and ob(e or f ) = ob(e) ⊕ ob(f ).

R

  = {(η(λx.σ(e 1 , . . . , e n )), η(σ(λx.e 1 , . . . , λx.e n )))} ∪ {(ϕ, ϕ) | ϕ ∈ T V} is an applicative simulation. The clause on ob follows by hypothesis on σ, whereas for the clause on monadic application we simply notice that λx.σ(e 1 , . . . , e n ) • k = σ V ( e 1 [k/x] , . . . , e n [k/x] ) = σ(λx.e 1 , . . . , λx.e n ) • k.

Definition 15 .

 15 Given an open relation R, the Howe extension R H of R is inductively defined by the rules in Figure2. The Howe extension of a closed relation R is defined as

Figure 2 :

 2 Figure 2: Howe extension.

Lemma 3 .

 3 Let R be a preorder closed under substitution. Then R H is a compatible, reflexive, and substitutive open relation such that R H ; R ⊆ R H and R ⊆ R H . Since is a preorder, and thus • is a preorder closed under substitution, by Lemma 3 ( ) H is a compatible and substitutive open relation containing . Before proving Theorem 3, we need to extend Howe's construction to evaluation contexts. Definition 16. Given an open relation R, the Howe extension R H

Figure 3 :

 3 Figure 3: Howe extension: evaluation contexts.

Case 3 .

 3 Suppose e is of the form G[σ(c 1 , . . . , c k )]. We proceed exactly as in case 2. From G[σ(c 1 , . . . , c k )] H f by Lemma 5 we obtain an evaluation context H and a term d such that H[d] f , G H E H, and σ(c 1 , . . . , c k ) H d. The latter must be the conclusion of an inference of the form:c i H b i σ(b 1 , . . . , b k ) d H-op σ(c 1 , . . . , c k ) H d As in previous case we observe that E H E F and G H E H imply E[G[-]] H E F [H[-]].Reasoning as above we can now apply the induction hypothesis, obtaining:ob( E[G[σ(c 1 , . . . , c k )]] n ) = σ 1 (ob( E[G[c 1 ]] n-1 ), . . . , ob( E[G[c k ]] n-1 )) σ 1 (ob( F [H[b 1 ]] ), . . . , ob( F [H[b k ]] )) = ob( F [H[σ(b 1 , . . . , b k )]] ) ob( F [H[d]] ) ob( F [f ] ).

  • . Proof. To see that app ⊆ ciu notice that, by the very definition of applicative simulation, e app f implies ∀k ∈ Λ. ek app f k. As a consequence, for any evaluation context E = [-] k we have e app f implies E[e] app E[f ], and thus ob(E[e]) ob(E[f ]). Conversely, we prove ciu ⊆ app by coinduction, showing that R {( e • k, f • k) | e ciu f, k ∈ Λ} is an applicative simulation. Clearly R is closed under the lifting of application. Furthermore, we see that ob( e • k) ob( f • k), since e • k = e k , f • k = f k , and e ciu f . Since app = ciu we obviously have ( app

  The thesis follows since is reflexive. If E is of the form Gk, then G[e]k H f must be the conclusion of an instance of rule H-app. As a consequence, we have G[e] H c, k H h, and ch f . By induction hypothesis, from G[e] H c we obtain an evaluation context H and a term g such that G HThe main technical part of the Howe's method is the so-called Key Lemma.Lemma 6 (Key Lemma). For all (closed) evaluation contexts E, F and programs e, f , if E H E F and e H f , then ob(E[e]) ob(F [f ]).Proof. Since ob is continuous and E[e] = n<ω E[e] n to prove the thesis it is sufficient to prove:

		E H, e H g, and H[g]	c. From the latter we infer
	H[g]h	ch, and thus H[g]h	f , since ch	f and	is transitive. From
	k H h and G H E H we infer Gk H E Hh. We are done by taking F ≡ Hh.

, then there exist an evaluation context F and a term g such that E H E F , e H g, and F [g] f .

Proof. We proceed by induction on E. If E = [-], then we take F = [-] and g = f .

  We can now conclude the thesis since λx.h f implies (λx.h)cc f cc.

	Case 2. Suppose e is of the form G[(λx.k)g]. By Lemma 5, from G[(λx.k)g] H
	f we obtain the existence of an evaluation context H and a term h
	such that G H E H, (λx.k)g H h, and H[h]

By Proposition 1 we have:

ob( e n ) = ob(η(λx.k)) = η 1 ( * ) = ob(η(λx.h)) ob( f ),

since λx.h f . If E is of the form [-]gḡ, then F must be of the form [-]cc with g H c and ḡ H c. Moreover, we have:

E[e] n = (λx.k)gḡ n = k[g/x]ḡ n-1 .

From {x} k H h and g H c we infer, by substitutivity, k

[g/x] H h[c/x]. The latter, together with ḡ H c, gives k[g/x]ḡ H h[c/x]c.

We can now apply the induction hypothesis obtaining ob( k[g/x]ḡ n-1 ) ob( h[c/x]c ) = ob( (λx.h)cc ).

  • = ctx . Proof. To prove compatibility of it is sufficient to show that H restricted to programs coincide with . The latter amounts to show that ∅ e H f implies ob(E[e]) ob(E[f ]), for any evaluation context E. This directly follows from Key Lemma, since E H E E always holds. Since • is obviously preadequate, we see that • ⊆ ctx . Since ctx ⊆ • trivially holds, we conclude • = ctx . Corollary 1 (Full abstraction). ( ciu

Although other choices, such as must or may and must convergence, are possible.

Intuitively, a context is an expression with a hole[-]. We write C[e] for the expression obtained by filling in the hole in C with e.

From a formal perspective, we should work with a bi-labelled transition system distinguishing between programs and values. However, since the goal of this section is to give an informal exposition of the main ideas behind this work, such a level of mathematical accuracy is unnecessary.

Using the terminology of[START_REF] Abramsky | The lazy lambda calculus[END_REF] we refer to transition systems of this form as nondeterministic applicative transition systems.

Formally, we should write {λx.(f g)} ( app ) * {λx.f λx.g}.

We use the expression monad and Kleisli triple interchangeably.

Recall that the support of µ ∈ [0, 1] X is the set supp(µ) = {x ∈ X | µ(x) = 0}. We write x µ(x) in place of x∈supp(µ) µ(x).

Oftentimes we will write s to denote a finite sequence s 1 , . . . , s n of symbols s i .

From a categorical perspective, our definition of lifting relies on the fact that any monad on the category Set is strong. That means that for every monad T = T, η, >>= on Set we can regard >>= as lifting any function f : X × Y → T Z to >>=f : X × T Y → T Z.


Mederic Fourmy

Thomas Flayols

Nicolas Mansard

Joan Solà

Contact forces pre-integration for the whole body estimation of legged robots

State estimation, in particular estimation of the base position, orientation and velocity, plays a big role in the efficiency of legged robot stabilization. The estimation of the base state is particularly important because of its strong correlation with the underactuated dynamics, i.e. the evolution of center of mass and angular momentum. Yet this estimation is typically done in two phases, first estimating the base state, then reconstructing the center of mass from the robot model. The underactuated dynamics is indeed not properly observed, and any bias in the model would not be corrected from the sensors. While it has already been observed that force measurements make such a bias observable, these are often only used for a binary estimation of the contact state. In this paper, we propose to simultaneously estimate the base and the underactuation state by using all measurements simultaneously. To this end, we propose several contributions to implement a complete state estimator using factor graphs. Contact forces altering the underactuated dynamics are pre-integrated using a novel adaptation of the IMU pre-integration method, which constitutes the principal contribution. IMU pre-integration is also used to measure the positional motion of the base. Encoder measurements are then participating to the estimation in two ways: by providing leg odometry displacements, contributing to the observability of IMU biases; and by relating the positional and centroidal states, thus connecting the whole graph and producing a tightly-coupled whole-body estimator. The validity of the approach is demonstrated on real data captured by the Solo12 quadruped robot.

I. INTRODUCTION

Legged robot estimation has from the start been focused on the integration of odometry derived from the kinematic chain [START_REF] Roston | Dead reckoning navigation for walking robots[END_REF], later fused with Inertial Measurement Units (IMU) to capture higher movement dynamics [START_REF] Lin | Sensor data fusion for body state estimation in a hexapod robot with dynamical gaits[END_REF]. This approach lead to the use of Bayesian filters introduced for quadruped robots in [START_REF] Bloesch | State estimation for legged robotsconsistent fusion of leg kinematics and imu[END_REF] and adapted for humanoid robots in [START_REF] Rotella | State estimation for a humanoid robot[END_REF]. Focusing on proprioceptive odometry, they naturally produced drifting estimation of the robot position and yaw and have been complemented by loosely coupled corrections from exteroceptive sensors [START_REF] Nobili | Heterogeneous sensor fusion for accurate state estimation of dynamic legged robots[END_REF].

Another line of research is the estimation of centroidal quantities such as the center of mass, centroidal momentum and their derivative which is especially important for the balance and predictive control of humanoid robots. They all have in common a two-step approach where a first estimator outputs the base state and this is used by a second estimator to obtain the centroidal quantities. While [START_REF] Stephens | State estimation for force-controlled humanoid balance using simple models in the presence of modeling error[END_REF] uses simplified models such as the LIPM, [START_REF] Rotella | Humanoid momentum estimation using sensed contact wrenches[END_REF] applies the full dynamical model to develop Kalman Filters using kinematic measurements fused with sensor wrenches while permitting the estimation of an external force or kinematic model bias. [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-a spectral approach[END_REF], [START_REF] Bailly | Recursive estimation of the human body's center of mass and angular momentum derivative[END_REF] propose a fusion based on a careful spectral analysis of the measurement sources from which a cascade of complementary filters is derived.

Coming from the drone community in the context of Visual Inertial Odometry [START_REF] Leutenegger | Keyframe-based visual-inertial odometry using nonlinear optimization[END_REF], more recent works start to introduce nonlinear optimization to solve real time estimation problems using the factor graph paradigm [START_REF] Dellaert | Factor graphs for robot perception[END_REF]. The seminal work [START_REF] Lupton | Efficient integration of inertial observations into visual slam without initialization[END_REF] showed that IMU measurements could be efficiently integrated in such a framework using the so called pre-integration, which was improved later with proper estimation on the SO(3) manifold [START_REF] Forster | On-manifold preintegration for real-time visual-inertial odometry[END_REF]. Ref. [START_REF] Hartley | Hybrid contact preintegration for visual-inertialcontact state estimation using factor graphs[END_REF] was the first work applying this framework to a biped robot (Cassie). It derived factors to tightly fuse IMU measurements with preintegrated leg odometry and off-the-shelf visual odometry. Shortly after in [START_REF] Wisth | Robust legged robot state estimation using factor graph optimization[END_REF], another team proposed similar methods on a quadruped robot (ANYmal) replacing leg odometry by velocity measurements coming from a lower level IMU kinematic filter. They then built on these results by introducing a new pre-integrated factor taking explicitly into account the bias from low level filter velocity, made observable by vision [START_REF]Preintegrated velocity bias estimation to overcome contact nonlinearities in legged robot odometry[END_REF]. The leg odometry factors used in these works rely implicitly on gyroscope measurement integration for point feet robots.

However, none of these factor graph based works included the robot dynamics into the estimation yet. The most closely related work is that proposed by [START_REF] Nisar | Vimo: Simultaneous visual inertial model-based odometry and force estimation[END_REF] in the context of aerial vehicles. They incorporated the linear dynamics of the rotors as a pre-integrated factor within a visual inertial framework which permits them to observe an external force being applied to the drone.

In this paper we describe a method to perform whole-body state estimation of legged robots. The estimated state contains positional states of the base (position, velocity and orientation), and centroidal states (CoM position and velocity, and angular momentum). The method fuses information from IMU, joint encoders and contact forces. A unique factor graph is employed to represent the estimation problem, yielding a tightly coupled fusion which is solved by nonlinear least squares optimization. One key contribution is the pre-integration of the force measurements into a unique measurement, in what constitutes an adaptation of the IMU pre-integration technique widely used in visual-inertial odometry and SLAM. We present results on a real dataset from Solo12 quadruped with 12 DoF, see Fig. 1.

II. PROBLEM STATEMENT

A. Centroidal dynamics

The robot dynamics is described by the well-known:

M(q) vq + h(q, v q ) = τ q + k J k f k (1)
where q, v q , vq , τ q are the position, velocity, acceleration and torques of the robot in configuration space, f k are the contact forces (written as 3D forces in this paper), M is the generalized inertia matrix, h the sum of gravity, Coriolis and centrifugal forces, and J k the jacobians of the contact points.

Because of the underactuated nature of legged robots, the configuration is often separated into q = (p b , R b , q a ) where p is the position in world frame of the robot basis (typically, the torso), R b the rotation of the basis body with respect to the world and q a are the joint configuration of the actuated joints.

While (1) represents the whole dynamics, a subpart of it is of particular importance for legged robots. The centroidal dynamics is written by the two equations:

mc = mg + k f k , L = k (p k -c) × f k (2)
where c, L are the position of the center of mass (CoM) and angular momentum (AM) around the CoM (both expressed in world frame), m is the robot total mass, and the p k are the positions of the contact points in world frame. The centroidal dynamics is an exact part of (1), but more deeply reveals the underactuation: the robot can move only if applying the proper forces to the environment, as the joint torques alone cannot lead to any modification of the Center of Mass (CoM) or angular momentum (AM). The classical approach in estimation of legged-robot state is to first estimate the basis state, and then reconstruct the centroidal state (c, ċ, L) using the joint position and velocity measurements, and the robot model. This assumes that there is no direct measurement of the centroidal state. Consequently, we would not be able to recover the exact centroidal state if there is any bias in the robot model.

Yet, we can see from the centroidal dynamics that the force measurements are connected to the variation of the centroidal state. As observed in [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-a spectral approach[END_REF], a proper fusion of the force measurements with an estimation of the state of the basis makes the centroidal state observable. As the basis state is also observable from inertial and kinematic measurements [START_REF] Bloesch | State estimation for legged robotsconsistent fusion of leg kinematics and imu[END_REF], a direct consequence is that the fusion of force, inertial and kinematic measurements makes the basis and centroidal states observable.

In the following, we propose a practical method to implement a basis and centroidal estimator based on IMU, joint coders and force sensors. Based on factor graph, this method can also be added into a SLAM estimator to also estimate the absolute positions of the basis and the CoM.

B. Estimation problem MAP formulation

Our 12 degrees of freedom (DoF) quadruped robot Solo12 is equipped with an IMU, joint encoders monitoring the joint angles q a ∈ R 12 and angle velocities qa ∈ R 12 . We will often drop the dependency of certain measurements on q a and qa for readability once stated. Joint torques can also be estimated from motor current. These quantities are available at 1kHz. The goal is to reconstruct the full state of the robot, namely the base position p ∈ R 3 , base velocity v ∈ R 3 and base orientation R ∈ SO(3), as well as the robot's CoM position c ∈ R 3 , the CoM velocity ċ ∈ R 3 and the centroidal AM L ∈ R 3 at this same frequency. All quantities are expressed with respect to the world frame w: We pose the estimation as an optimization problem where variables χ are composed of a trajectory of states x i at times t i , also called Key Frames (KF), to which the biases b i are associated. This problem can be represented as a factor graph (Fig. 2) which is a bipartite graph where nodes correspond to the variables and factors to nonlinear constraints coming from the noisy sensor measurements. Then, under the assumption that χ has a multivariate Gaussian distribution and that measurement noises are Gaussian and independent, the maximum a-posteriori (MAP) problem is equivalent to solving the weighted least squares problem,

x = [p, v, R, c, ċ, L] [ w p b , w v b , w R b , w p c , w v c , w L]
χ * = arg min χ k r k (χ) 2 Σ k (3)
Each factor k corresponds to a residual function r() of any number of variables from χ and a covariance matrix Σ associated to sensor noises. In the remaining of the paper, we will refer to measurement quantities with a tilde, e.g. ã, and Gaussian noises affecting them as n a . The following section describes how these residuals are obtained from measurements.

III. FACTOR GRAPH FOR BASE-CENTROIDAL ESTIMATION

In the context of graph optimization estimation, factors can apply to a single KF (unary factors) or distant KFs. KFs are created at relatively low rate to avoid saturating the estimation problem with too many variables. We incorporate several kinds of factor to (3): state prior (0), leg odometry (LO), pre-integrated IMU deltas (I), pre-integrated force measurements (F ), centroidal kinematics (CK), and bias drifts (B), yielding the optimization problem (Fig. 2),

χ * = arg min χ r 0 i 2 Σ I i + i,j r I i,j 2
Σ I i,j + i,j r LO i,j 2
Σ LO i + i,j r F i,j 2
Σ F i,j + i r CK i 2 Σ CK i + i,j r B i,j 2
Σ B i,j (4)
These factors are described in the following paragraphs.

A. State priors r 0

In the posed problem, the base position and yaw angle are unobservable since we do not have access to absolute measurements (coming from GPS for instance). In order to ground the estimates to a reference and not have a floating problem, we impose a prior on the initial position and orientation of the base.

r 0 (x 0 ) = p 0 -p 0 Log(R 0 R 0 ) (5)
where Log denotes the logarithm application mapping elements of SO(3) to the R 3 representation of its Lie algebra so(3).

B. Leg odometry factor r LO

Assuming that we have access to an accurate contact detection, for each foot in contact between KFi and KFj, we can write that w p i wl = w p k wl + n LO where n LO denotes a Gaussian noise accounting for potential slip of the foot and kinematic inaccuracies. This pseudo measurement can be used to derive a 3D factor for each foot l in contact:

r LO l (p i , R i , p k , R k ) = p i + R i b p i bl -(p k + R k b p k bl ) (6
) Where b p i bl and b p k bl are the contact position in base frame at times i and j acquired from q a via forward kinematics.

Assuming good contact detection, leg odometry fused with the IMU is enough to properly observe 6DoF motion with as few as 2 contact feet at a time. Note that contrary to previous works [START_REF] Hartley | Hybrid contact preintegration for visual-inertialcontact state estimation using factor graphs[END_REF], [START_REF]Preintegrated velocity bias estimation to overcome contact nonlinearities in legged robot odometry[END_REF], this leg odometry factor does not require gyroscope integration C. Generalized pre-integration theory Pre-integration refers to the integration of high rate proprioceptive sensory data in an efficient way in the context of factor graphs. It was initially developed to deal with IMU measurements [START_REF] Lupton | Efficient integration of inertial observations into visual slam without initialization[END_REF], [START_REF] Forster | On-manifold preintegration for real-time visual-inertial odometry[END_REF]. In fact, if a standard integration is conducted naively to derive a factor measurement between two KFs, then the IMU data needs to be reintegrated at each solver iteration because the integral depends on the states we are to estimate. The IMU pre-integration theory solves this issue and can be generalized to any proprioceptive sensor as shown in [START_REF] Atchuthan | Towards new sensing capabilities for legged locomotion using real-time state estimation with low-cost IMUs[END_REF]- [START_REF] Fourmy | Absolute humanoid localization and mapping based on imu lie group and fiducial markers[END_REF]. Delta quantities ∆ ik are defined between KFs i and k so that x k = x i ∆ ik . A suitable composition is chosen so that ∆ ik does not depend on the initial states x i . Dependence on other, small-varying parameters b such as sensor bias is linearized as

∆ ik (b) = ∆ ik (b) ⊕ J ∆ b (b -b).
Then we can pre-integrate ∆ ik ∆ ik (b) once during data gathering, and use it to define residuals that are later evaluated many times by the optimizer.

1) Delta pre-integration: We perform pre-integration incrementally as follows. First, ∆ ii is initialized to the null motion. Its covariance Σ ∆ ik and the Jacobian J δ k bi are set to zero. At each reception of sensor data zk at t k , we integrate during δt to obtain the delta corresponding to a single data sample,

δ k = f (z k , b i , δt) , (7)
using the bias b i available in KF i. This single delta is integrated onto the delta pre-integrated so far using the delta composition law,

∆ ik = ∆ ij • δ k . (8)
The Jacobian of the pre-integrated delta with respect to biases, as well as the delta covariance, are also pre-integrated,

J ∆ ik bi = J ∆ ik ∆ij J ∆ij bi + J ∆ ik δ k J δ k bi (9)
Σ ∆ ik = J ∆ ik ∆ij Σ ∆ ij J ∆ ik ∆ij + J ∆ ik δ k J δ k z k Σ z J δ k z k J ∆ ik δ k (10)
where

J δ k bi , J δ k z k , J ∆ ik ∆ij and J ∆ ik δ k
are the Jacobians J y x = ∂y/∂x of (7,8), computed according to Lie theory [START_REF] Solà | A micro Lie theory for state estimation in robotics[END_REF].

2) Residual definition: The pre-integrated ∆ ik is used at the end of the pre-integration to define the residual,

r ik = (∆ ik ⊕ J ∆ ik bi (b i -b i )) ∆ik (11)
where b i is the current bias value, ∆ik = x k x i is the expected delta between keyframes, and {⊕, } are the composite plus and minus operators described in [START_REF] Solà | A micro Lie theory for state estimation in robotics[END_REF]. That is, {⊕, } are {+, -} for vectors, and for rotations we have

R ⊕ θ R Exp(θ) and R 2 R 1 Log(R 1 R 2 )
. The residual clearly depends on the keyframe states x i , x k and the bias b i . It has an associated covariance Σ ∆ ik .

D. IMU pre-integration factor r I

The states involved in this integration are the base states The pre-integration method in [START_REF] Forster | On-manifold preintegration for real-time visual-inertial odometry[END_REF] can be put in the formalism above by realizing that defining δ = f (z, b, δt) in [START_REF] Rotella | Humanoid momentum estimation using sensed contact wrenches[END_REF] as

x b = [p, v,
δ k =   δp δv δR   k =   1 2 (ã -b a -n a )δt 2 (ã -b a -n a )δt Exp(( ω -b ω -n ω )δt)   k (12)
and the composition law [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-a spectral approach[END_REF] as [13, eq. 32] and ∆ ik = x k x i is [13, eq. 33]. Full details can be found in [START_REF] Atchuthan | Towards new sensing capabilities for legged locomotion using real-time state estimation with low-cost IMUs[END_REF]Section 3.4].

∆ ik = ∆ ij • δ k in
∆p ik = ∆p ij + ∆v ij δt + ∆R ij δp k ∆v ik = ∆v ij + ∆R ij δv k ∆R ik = ∆R ij δR k (13) then x k = x i ∆ ik is

E. Force pre-integration factor r F

The Newton-Euler equations ( 2) relate the evolution of a system's CoM and AM due to gravity, external forces and torques in a chosen inertial frame. In the case of a legged robot with punctual contact feet, only forces f l are applied at each limb contact l, with no torque. We assume that at each limb contact we have access to noisy local force measurements f l = l f l + n f . To transform them to the body frame b, we compute b R l (q a ) ∈ SO(3) and b p bl (q a ) ∈ R 3 from the articulations configuration q a ∈ R 12 . The lever arm (p l -c) in (2) uses a measurement of the CoM position in base frame b p bc (q a ) ∈ R 3 . This measure is biased due to inaccuracies in the robot model and therefore we add a bias variable to estimate, b bc ∈ R 3 so that b p bc = b p bc +b bc +n bc . Assuming constant forces during the interval δt the integral of (2) yields the discrete dynamic model for the centroidal states:

c k = c k-1 + ċk-1 δt+ 1 2 gδt 2 + 1 2m R k-1 l b R k l ( f k l -n k f )δt 2 ċk = ċk-1 + gδt + 1 m R k l B R k L ( f k l -n k f )δt L k = L k-1 + R k l ( b p k bl -b p k bc +b k bc +n k bc )× b R k l ( f k l -n k f )δt (14)
Analogously to the IMU case, it is possible to pre-integrate force measurements to derive a factor on the states

x c = [c, ċ, L, R] using deltas ∆ c = [∆c, ∆ ċ, ∆L, ∆R].
The rotation measured by the gyroscope has to be included too for the pre-integration to work. Then the bias vector is

b = [b bc , b ω ].
We define measurements z k to be:

z k = b p k bc , ω ω ω k , f k l , b p k bl , b R k l l=1...4 (15)
The following operators are enough to particularize the method in Section III-C to the contact forces case. Integrating

z k during δt yields δ k in (7) as, δ k (z k , b i , δt) =     1 2m l b R k l f k l δt 2 1 m l b R k l f k l δt ( l b p k bl -( b p k bc -b bc )×( b R k l f k l ))δt Exp(( ω ω ω k -bω)δt)     (16)
and the delta composition law in [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-a spectral approach[END_REF] as

∆ • δ =     ∆c + ∆ ċδt + ∆Rδc ∆ ċ + ∆Rδ ċ ∆L + ∆RδL ∆RδR     . (17)
The expected delta ∆ ik = x k x i between KFs i, k reads

∆ ik =     ∆c ik ∆ ċik ∆L ik ∆ R ik     =      R i (c k -c i -ċi ∆t ik ) R i ( ċk -ċi -g∆t ik ) R i (L k -L i ) R i R k      . (18)
Finally, the propagation of the state x i to x k using ∆ ik , which can be used to retrieve a state at any time k between keyframes in the trajectory, is

x k = x i ∆ ik =     c i + ċi ∆t ik + R i ∆c ik + 1 2 g∆t 2 ik ċi + R i ∆ ċik + g∆t ik L i + R i ∆L ik R i ∆R ik     (19)
F. Centroidal kinematics factor r CK

We need to relate base states to centroidal quantities to ground their estimate. For that, we can rely on the inertialkinematic model for the computation of the CoM position wrt. base frame b p bc (q a ) ∈ R 3 , the CoM velocity wrt. base frame b v bc (q a , qa ) ∈ R 3 , the inertial matrix I(q a ) ∈ R 3×3 , and the kinematic momentum due to gesticulation of the robot limbs only L a (q a , qa ) ∈ R 3 . In particular, the computed CoM position is here considered an independent measurement, affected by Gaussian noise and a slowly varying bias b p bc = b p bc + b bc + n bc . The angular velocity from the IMU is used and its bias has to be incorporated in the factor b ω ω ω b = b ω ω ω b + b ω + n ω . In the end the equations used to derive the factor are:

c = R( b p bc -bpc -npc) + p ċ = v + R(( b ω ω ω b -bω -nω) × ( b p bc -b bc -n bc ) + ( b v bc -nvc)) L = R(I( b ω ω ω b -bω -nω) + La) (20)
Then, the residual r CK ∈ R 9 is simply expressed as: We tested our solution on the open source robot Solo12 [START_REF] Grimminger | An open torque-controlled modular robot architecture for legged locomotion research[END_REF] built recently at LAAS. Solo12 is not equipped with contact nor force sensors. Yet, it is possible to invert the robot dynamic equation ( 1) in order to reconstruct the contact forces knowing the robot configuration and derivatives vq , v q , q and joint torques given by the motor current. Some of these quantities are hard to obtain directly since they depend on the state being estimated (base orientation) or on numerical differentiation ( qa ). For these reasons, we precalculated these forces by benefiting from an internal filter of the IMU for the base orientation, and centered window differentiation for the articulation acceleration qa . Fig. 3 shows an example of the force reconstruction of one leg.

r CK =   b p bc -(R T (c -p) + bp) ( b ω ω ω b -bω) × ( b p bc -bp) + b v bc -R T ( ċ -v) I( b ω ω ω b -bω) + La -R T L   (21)

V. RESULTS

In this section we present results of the application of the estimator on datasets taken at LAAS-Toulouse on our platform Solo12: sinXYZ corresponds to a trajectory where the robot moves by following a sine wave reference with the feet fixed in place. We use the state estimation framework WOLF, developed at IRI-Barcelona, which relies on Google Ceres as the graph solver. For all runs, Key Frames were created at 4Hz for fair comparison.

A. Base estimation through inertial kinematic fusion

First, to validate the use of our kinematic factor, we include uniquely the IMU and LO factors to obtain an Inertial Kinematics (IK) estimator which conceptually includes the same information as estimators such as [START_REF] Bloesch | State estimation for legged robotsconsistent fusion of leg kinematics and imu[END_REF]. In Fig. 4, we compare our state estimation at 1kHz with Motion Capture (Mo-Cap) up-sampled from 200Hz to 1kHz. Velocity in base frame is also shown in Fig. 5.

Artificially removing contact factors (considering only 1, 2 or 3 feet in contact) can help us gain confidence in the use of this kinematic factor in situations where we rarely have all feet in contact, like for example with trotting gaits. In Fig. 6, we can see that only considering 1 foot in contact during the whole trajectory results in a drifting position, but as soon as 2 or more feet are in contact, the system is constrained enough for the drift to remain below around 5mm on all axes.

B. Centroidal estimation

Now, on the same trajectory, we deploy the full estimator with all factors described in Section III to jointly estimate the base and centroidal quantities. A ground truth on the centroidal quantities is difficult to obtain since no direct sensor can provide us with this information contrary to the base state. We can however validate our method by comparing it a two step procedure: first, estimate the base state with a state Kalman Filter as implemented in [START_REF] Bledt | Mit cheetah 3: Design and control of a robust, dynamic quadruped robot[END_REF], then compute the centroidal quantities directly from the robot kinematic model. The full estimator should be able to infer a bias on the b p bc measure so we artificially add a constant disturbance in the robot dynamic model on the lever of the base link of [0.03, 0.06, 0.04] cm, which then corresponds to a CoM bias of [-0.0197, -0.0394, 0.0263]. Fig. 8 shows that the bias estimated with our method closely matches the introduced bias. Fig. 9 shows a comparison between the base and CoM reconstruction with our method and with the twostep base KF with geometric CoM Note that the base-CoM difference on the z axis reflects the fact that limbs of the robot lower naturally its CoM. The estimated CoM velocity follows closely the velocity of the base as shown in Fig. 7. In this paper, we introduced the first realization of the tightly coupled estimation of base and centroidal quantities of a legged robot using factor graph optimization. This was enabled by the introduction of a new pre-integrated factor on the contact forces of the robot, making the CoM geometric measurement bias observable. We also designed a simple factor for point leg sensors that does not rely on gyro measurement. We demonstrated the viability of the approach on a dataset taken on the quadruped robot Solo 12. Within this framework, we plan to integrate exteroceptive measurement in a SLAM context to obtain a drift free estimation of base and centroidal states useful for planning and feedback control.

Fig. 1 :

 1 Fig. 1: LAAS-Gepetto robot Solo12 from the Open Dynamic Robot Initiative

  Some measurements are affected by time varying biases b = [b a , b ω , b p ] that can be estimated, where b a /b ω refer respectively to the accelerometer/gyro biases and b p is the bias on the CoM centroidal measurement.

Fig. 2 :

 2 Fig. 2: Factor graph including factors of IMU (green), leg odometry (red), centroidal kinematics (purple) and pre-integrated force-torque (blue). Columns are keyframes. State block nodes are p: position; R: orientation; v: velocity; ba: acc. bias; bw: gyro bias; c: CoM position; ċ: CoM velocity; L: angular momentum; bp: CoM measurement bias. Single-ended factors are priors.

  R] with deltas ∆ = [∆p, ∆v, ∆R]. The IMU produces biased and noisy measurements z = [ã, ω] of the base proper acceleration and angular velocity, with bias b = [b a , b ω ] and noise n = [n a , n ω ].

G

  . Bias prior and drift factors Bias prior and drift factors have trivial residuals r prior = bb ref and r drift = b k -b i respectively, each with its own associated covariance Σ.

Fig. 3 :

 3 Fig. 3: Force estimation on X-Y-Z axis (r-g-b) of one Solo12 leg expressed in world frame using proprioceptive sensors during the sinXYZ trajectory

Fig. 4 :Fig. 5 :

 45 Fig. 4: sinXYZ trajectory base position from the IMU+kinematics estimator (blue) vs Mo-Cap (red)

Fig. 6 :Fig. 7 :

 67 Fig. 6: sinXYZ trajectory base position error with different numbers of feet used for the leg odometry factors: 1 (blue), 2 (orange), 3 (green), 4 (red) from the IMU+kinematics estimator

Fig. 8 :Fig. 9 :

 89 Fig.8: Estimation of bias on CoM measurement from the full estimator along x-y-z axis ( red-green-blue) in base frame

This work has been supported by the MEMMO European Union project within the H2020 Program under Grant Agreement No. 780684, the Spanish Ministry of Science, Innovation, and Universities project EB-SLAM (DPI2017-89564-P), by the EU H2020 project GAUSS (DPI2017-89564-P) and by the Spanish State Research Agency through the Maria de Maeztu Seal of Excellence to IRI MDM-2016-0656.