Gilles Chabert
email: gilles.chabert@emn.fr

Luc Jaulin
email: luc.jaulin@ensieta.fr

Contractor Programming

This paper describes a solver programming method, called contractor programming, that copes with two issues related to constraint processing over the reals. First, continuous constraints involve an inevitable step of solver design. Existing softwares provide an insufficient answer by restricting users to choose among a list of fixed strategies. Our first contribution is to give more freedom in solver design by introducing programming concepts where only configuration parameters were previously available. Programming consists in applying operators (intersection, composition, etc.) on algorithms called contractors that are somehow similar to propagators. Second, many problems with real variables cannot be cast as the search for vectors simultaneously satisfying the set of constraints, but a large variety of different outputs may be demanded from a set of constraints (e.g., a paving with boxes inside and outside of the solution set). These outputs can actually be viewed as the result of different contractors working concurrently on the same search space, with a bisection procedure intervening in case of deadlock. Such algorithms (which are not strictly speaking solvers) will be made easy to build thanks to a new branch & prune system, called paver. Thus, this paper gives a way to deal harmoniously with a larger set of problems while giving a fine control on the solving mechanisms. The contractor formalism and the paver system are the two contributions. The approach is motivated and justified through different cases of study. An implementation of this framework named Quimper is also presented.

Introduction

Constraint programming is a simple and efficient paradigm to handle a large class of combinatorial problems [START_REF] Rossi | Handbook of Constraint Programming (Foundations of Artificial Intelligence)[END_REF][START_REF] Dechter | Constraint Processing[END_REF][START_REF] Van Hentenryck | The OPL optimization programming language[END_REF]. In the presence of real-valued variables, constraint propagation algorithms combined with interval analysis [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF][START_REF] Kearfott | Rigorous Global Search: Continuous Problems[END_REF][START_REF] Hyvönen | Constraint Reasoning Based on Interval Arithmetic: The Tolerance Propagation Approach[END_REF][START_REF] Davis | Constraint Propagation with Interval Labels[END_REF] are also particularly well-suited, included for real-world applications (see, e.g., [START_REF] Merlet | Solving the Forward Kinematics of a Gough-type Parallel Manipulator with Interval Analysis[END_REF][START_REF] Jaulin | Localization of an Underwater Robot using Interval Constraint Propagation[END_REF]). We shall refer to this interval variant of constraint programming as interval programming. Even then, interval programming has had only moderate success. In our opinion, the reason is a lack of clear and unified formalism describing how solvers and derived programs are built. This paper is an attempt to fill this gap.

We propose a framework that allows one to build a continuous solver with a few lines, in a high-level syntax. More than just another tuning language, a programming framework is proposed.

Motivation

Three reasons justify the introduction of a solver programming framework in presence of real variables.

1. Constraint programming is a declarative paradigm which means that a programmer should spend most of the effort in modeling conveniently the problem. This effort may involve breaking symmetries, introducing global or soft constraints, etc.. All these concepts are related to modeling and represent, by the way, very active fields of research. In theory, the solver is a black box of which a programmer could ignore the details. Despite of this, there is always a need at some point to control the solver, as it is with all declarative languages (consider for instance Prolog cuts that allow ruling out choices in the search). We can say that the overall efficiency one gets is the combined result of efforts made on both aspects: modeling and solver control. With real variables, modeling languages are limited 1 . Constraints can usually be nothing but equations which means that mathematics is de facto the ultimate modeling language 2 . The consequence is that solver control becomes an inevitable step if one is to improve efficiency. 2. Continuous solvers have a two-layered structure, namely interval analysis and constraint programming, The lower layer includes interval arithmetics and interval numerical algorithms (e.g., an interval variant of the Newton iteration) with round-off considerations. The upper layer includes branch & prune algorithms for describing sets of reals defined by constraints (or optimizing a criterion under constraints). Since these layers correspond to quite different scientific communities, there is a need more than ever for an interface between them. In concrete words, it would be useful to give a constraint programmer the ability to develop a continuous solver without digging into details of interval analysis. 3. The last and possibly main point is related to the output of constraint solvers.

With discrete domains, the output is always the set of solutions (or a subset optimizing some criteria). But in continuous domains, there may be a large variety of different outputs. First, one may look for a sub-paving (a set of boxes) encompassing the solutions and this is precisely what most of the existing solvers provide. They act as root finders. Next, in case of a solution set with a non-null volume, several sub-pavings are expected, each satisfying a different property, basically: "may contain a solution", "does contain a solution" or "contains only solutions". Such solvers rather act as set describers. Actually, we will see through the examples of Section 2 that the semantics behind the sub-pavings may completely change from one problem to the other 3 . As we will show, neither a root finder nor a set describer is adapted for solving these problems. Of course, ad-hoc solutions always exist, but the purpose is precisely to avoid a multiplicity of programs where a single one would be enough. In practice, when people are facing a specific constraint problem that requires a specific algorithm, they have to reverse-engineer the code of an existing solver. Often, they redevelop it from scratch.

Contribution

We propose a formalism and an algorithm, called paver.

In the formalism, the different interval routines (evaluations, projections, existence tests, etc.) are all wrapped in the very same object called contractor (see Section 3). Of course, the concept of contractor is not a novelty on its own. Our (first) contribution is to redefine various constraint programming techniques (propagation, shaving, parameter splitting, etc.) as operations over contractors that yield new contractors (Section 4). Syntactically, the contractor is then the unique atom, whence a certain simplicity. A solver can then be programed, rather than configured, by combining different contractors (examples are given in §4.2, §4.3, §5.1 and §5.5).

The paver algorithm is a generic solver. It takes a list of contractors, an initial box and follows a classical recursion: the contractors are successively called on the current box until either it gets empty or no more contraction could be done.

In the latter case, the box is bisected and contractors are called back again.

The fact that different contractors work concurrently allows solving problems of quite different nature (see next section). This is our second contribution.

Hence, contractor programming consists in two distinct steps: contractor design and paver design. The former refers to the design of the most possible efficient contractor for a given set (constraint) and shall be discussed in Section 4. The latter refers to the selection of contractors that yield the desired output, regardless of efficiency.

This framework is already supported by a real system named Quimper that will be introduced in Section 6.

Thinking of contractor programming as an extension of constraint programming is valid to the extent that contractors help in modeling the output of a problem. But, fundamentally, there is not such an extension since constraints basically tell the "what" whereas contractors tell the "how".

Note that branching will not be covered in this paper but one has to keep in mind that this part of the paver should be customizable as well. Only plain bisection will be used in the examples (variable are selected with a round-robin heuristic).

Power of Contractor programming

As for every programming language, the power of contractor programming, i.e., the class of problems that can be solved using this paradigm would require the setting of computability theory to be described formally. Pavings obtained in our framework are the results of algorithms that recursively contract boxes to smaller boxes and this general definition does not even involve the concept of constraint. It can also fit well, e.g., in the context of computer graphics to draw fractals. This means that some pavings generated by our system cannot be characterized by constraints.

However, one may wonder which problems contractor programming is the most aimed at. The contractor paradigm is mostly aimed at building set computation tools over the reals. These tools typically calculate disjunction of equidimensional sets and sets are usually described by constraints. There are however very few restriction on their form: constraints can be first-order logic formulas involving numerical equations or inequalities and quantified parameters. They can also be handled implicitly by their associated contractor.

The main point is that the representation of the sets and the way it is calculated are both entirely controllable.

Related works

In discrete constraint programming, existing tools already give control mechanisms with sometimes an associated formalism. In constrained-based local search, the COMET system [START_REF] Van Hentenryck | Constraint-Based Local Search[END_REF] offers a language allowing one to give his own definition of moves and neighborhoods, the incremental computation of the impact after a move being done automatically.

In branch & prune systems, choice points and propagation can usually be controlled. The Choco [START_REF] Jussien | The CHOCO Constraint Programming Solver[END_REF], ILOG Solver [START_REF] Ilog | ILOG Solver[END_REF] and Gecode [START_REF] Schulte | Gecode[END_REF] libraries allow to specify the instantiation order for both variables and values. It can also be decided how to split domains (either instantiate or bisect domains, etc.).

Propagation has given rise to more sophisticated concepts. The aforementioned libraries allow the programmer to associate his own propagator, or filtering algorithm, to a constraint. In Gecode, propagators can be automatically generated from a formula. It is possible with Choco to control when propagators are launched thanks to a system of events (e.g., a bound is reduced, the variable is instantiated, etc.). This allows, as a side effect, to schedule constraints in the propagation loop (some costly propagators may ignore events of minor importance). Events handling have to be implemented by the propagators in Gecode, albeit recent works about advisors [START_REF] Lagerkvist | Advisors for Incremental Propagation[END_REF] facilitate this task. These mechanisms are all related to modeling (a propagator is associated to a constraint), preserving declarativity. As far as we know, there is not a way to build new propagators, say, by composing different propagators, at the modeling phase. Only programmers with a good knowledge of the target library are able to do so. This is a difference with contractor programming that aims at building solvers rather than operational constraints.

In contrast, efficiency has always taken precedence over flexibility in continuous constraint programming. Perhaps because of the prevailing numerical culture, solvers are indeed usually compared on a performance basis. As a consequence, none of the different existing implementations (Numerica [START_REF] Van Hentenryck | Numerica: A Modeling Language for Global Optimization[END_REF], RealPaver [START_REF] Granvilliers | Algorithm 852: RealPaver: An Interval Solver using Constraint Satisfaction Techniques[END_REF], Alias [START_REF] Merlet | Alias[END_REF], Rsolver [START_REF] Ratschan | RSolver[END_REF], GlobSol [START_REF] Kearfott | GlobSol[END_REF], Baron [START_REF] Sahinidis | Branch-And-Reduce Optimization Navigator[END_REF], etc.) supply a good level of openness or extensibility. These libraries contain solvers or global optimization algorithms with a very high-level interface, basically restricting users to enter a system of equations, set up some parameters and press enter. Of course, there may be many parameters. For instance, it can be decided whether Alias has to resort to the computation of the Hessian matrix or to ask RealPaver for the global consistency instead of a solution set. But there is no way to venture off the beaten tracks.

A side contribution of this paper is a new representation for solution spaces with non-empty interiors. More suitable representations have been proposed in [START_REF] Sam-Haroud | Consistency Techniques for Continuous Constraints[END_REF] using quadtrees/octrees or in [START_REF] Vu | Rigorous Solution Techniques for Numerical Constraint Satisfaction Problems[END_REF] using extreme vertices. However, quadtrees consist in discretizing constraints and reasoning using this auxiliary representation. This otherwise appealing method is therefore quite distant from interval programming. Results of this paper are unlikely to be applicable in the context of quadtrees.

Computing extreme vertices can be viewed as an independent operation to obtain a compact representation of a set of boxes. The dissertation [START_REF] Vu | Rigorous Solution Techniques for Numerical Constraint Satisfaction Problems[END_REF] contains also smart branching heuristics and a concept of active variables to preserve the alignment of small boxes. All these concepts are fully compatible with our framework and somehow orthogonal to contractors since they are related to branching issues. They can be viewed as complementary to ours.

Notations and vocabulary

The set of all intervals is denoted by IR.

A A sub-paving of [x] is a set of non-overlapping boxes included in [x]. A subpaving will either be considered as a collection of boxes {[x] (1) , . . . , } or as a union [x] (1) ∪ . . . depending on the context. Hence, a sub-paving can either be viewed as a discrete subset of IR n or as a compact subset of R n .

A paving of [x] is a collection of sub-pavings K 1 , . . . , K N such that [x] = 1≤k≤N [b]∈K k [b].
Finally, we will temporarily adopt the following definition. Let c be a constraint on the reals, i.e. c(x) is either true or false. A contractor for c is a function C : IR n → IR n , usually given by a polynomial-time algorithm, such that:

∀[b] ∈ IR n , C([b]) ⊆ [b] and ∀x ∈ [b]\C([b]) c(x) is not satisfied.
Note that the concept of contractor will be generalized in our formalism (see Section 3). Thanks to interval arithmetics, many contractors can be built with respect to an equation f (x) = 0, where f : R n → R m .

Example 1. Consider f : R 2 → R with f (x) = x 1 -exp(x 2 ) and define C :

IR 2 → IR 2 as follows: C([x]) := ([x 1 ] ∩ exp([x 2 ]), [x 2 ] ∩ log[x 1 ]). C is a contractor with respect to f (x) = 0.
Contractor versus propagator. In spite of appearance, contractor is not just a new name for propagator. A propagator is associated to a constraint and destined to be called in a propagation loop. As we will see, a contractor is not necessarily associated to a constraint and may have nothing to do with propagation. Furthermore, we will add a continuity condition in the definition of a contractor (see §3.2) that does not exist for propagators. Finally, the term filtering was not chosen because of a conflict with the everyday meaning of filtering in many areas dealing with continuous systems (e.g., Kalman filtering).

The ABC of interval programming

Let us end this introduction with a short description of what interval programming is all about.

The basic algorithm of interval programming is the root finder. Consider a system of equations f (x) = 0, where x is a variable in R n and f a mapping from R n to R m . A root finder wraps solutions into boxes of any desired precision (leaving aside floating-point considerations) by using a branch & prune process, where the prune operation is performed by a contractor.

Here is the sketch of the algorithm. Note that a stack is used for the depth-first search and that the precision criterion chosen here is width[x] < ε.

algorithm root-finder(function f, box [x] (0) ) output: a sub-paving push [x] (0) on the stack while the stack is not empty (1) and [x] (2) on the stack Of course, the existing implementations are much more sophisticated. But this simplified description is enough from a language perspective. In the ideal case of a square system (m = n) without singularity, solutions are punctual (and certification can be demanded). In the general case, a nm dimensional surface is described. This is illustrated in Figure 1.

| pop [x] from the stack | [x] ′ ← C([x]) where C is a contractor w.r.t. f (x) = 0 | if [x] ′ is not empty | | if width [x] ′ < ε | | | insert [x] ′ into the output sub-paving | | else | | | bisect [x] ′ into two subboxes [x] (1) and [x] (2) | | | push [x]
Consider now f : R n → R m , and [x]× [y] ∈ R n+m . We may look for the following set inverse:

f -1 ([y]) ∩ [x] := {x ∈ [x] | ∃y ∈ [y] y = f (x)}.
Note that set inversion is a particular subclass of quantified constraints (see, e.g., [START_REF] Goldsztejn | A branch and prune algorithm for the approximation of non-linear AE-solution sets[END_REF][START_REF] Ratschan | Quantified constraints under perturbation[END_REF]). The root-finder algorithm will split the whole solution set into small boxes, as depicted in Figure 2.a.

In this case (i.e., when the solution set has a non-null n-dimensional volume) we might expect the solver to dissociate boxes that may contain solutions (boundary boxes) from boxes that only contain solutions (inner boxes). Clearly, precision is only required for boundary boxes, splitting an inner box being useless and even counter-productive.

A more appropriate output is given in Figure 2.b. Hence, an enhanced algorithm set-describer can be derived from root-finder by inserting an inner test before the call to the contractor. The design of inner tests has been a matter of study in several publications (see, e.g., [START_REF] Goldsztejn | Inner and Outer Approximations of Existentially Quantified Equality Constraints[END_REF][START_REF] Herrero | Quantified Set Inversion Algorithm with Applications to Control[END_REF][START_REF] Grandón | Generalized Interval Projection: A New Technique for Consistent Domain Extension[END_REF]).

Motivating examples

The purpose of this section is to present five different problems tractable by interval programming and to show that neither root-finder nor set-describer are appropriate. Once again, we do not claim that no algorithm exists for these problems. We just claim that no algorithm exist that solves them all. Note that root-finder and set-describer will not be two different algorithms in our formalism.

A decomposable problem

In several applications, systems are sufficiently sparse to be decomposed by equational or geometric techniques [START_REF] Jermann | Decomposition of Geometric Constraint Systems: a Survey[END_REF][START_REF] Neveu | When Interval Analysis Helps Inter-Block Backtracking[END_REF][START_REF] Bliek | Using Graph Decomposition for Solving Continuous CSPs[END_REF]. Consider for example the following system of distance constraints in two dimensions. There are 5 points connected by 6 distance equations. Points A and B are fixed and the three others are the unknowns. Instead of directly solving the 6x6 overall system, we can process progressively in three steps. We can solve the first triangle {A, B, C}, which is a 2 × 2 system. Because of rigidity, there is only a finite number of solutions (two flips). Then, for each solution, the second triangle {B, C, D} can be solved because C is fixed. Finally, the triangle {C, D, E} can be solved in turn once C and D are fixed. In general, solving the decomposed problem (here, seven 2 × 2 systems including flips) is faster than solving the global problem (here, one 6 × 6 system). No standard root finder supports the implementation of such strategies. The state of the art in decomposition-based solving methods relies on dedicated algorithms.

Enhanced description of the ring

Let us consider again the set inversion problem introduced in §1.2. We mentioned that an inner test was necessary (whence the set-describer algorithm) but an inner contractor could also be introduced in this case and would actually be much more adapted. Indeed, the set inversion problem can be easily formulated as a conjunction of two inequalities. It has been shown (see [START_REF] Collavizza | Extending Consistent Domains of Numeric CSP[END_REF][START_REF] Benhamou | Universally Quantified Interval Constraints[END_REF]) in the case of inequalities that, by a negation trick, the problem can be reversed so that any contracted region of the reversed problem lies inside the feasible region of the original one.

Remember that [y] ∈ IR. We have (1) and [x] (2) are the two resulting boxes, then [x]\([x (1) ]

c(x) ⇐⇒ f (x) ∈ [y] ⇐⇒ y -≤ f (x) and f (x) ≤ y + Hence, a point is unfeasible if y -> f (x) or f (x) > y + . B C D E A Fig. 3. A decomposable system. Now, we can contract a box [x] with respect to y -> f (x) and f (x) > y + . If [x]
∪ [x (2) ]) is inside the solution set {x ∈ [x] | c(x)}.
In such situation, it would be fruitful to replace the current search space [x] by [x (1) ]∪[x (2) ] and to memorize that the complementary space is inner. The set-describer algorithm is therefore not really appropriate.

Level surfaces

Level surfaces is a generalization of the previous problem, and can be stated as follows. Consider f : R n → R and N decreasing values y 1 > . . . > y N . The goal is to classify all the points of an initial domain [x] according to the constraints:

f (x) ≥ y 1 or f (x) ≥ y 2 . . . or f (x) ≥ y N , (1)
with a priority corresponding to the intuitive idea that we are more interested in the highest values: the surface level y 1 , i.e., the number of points satisfying f (x) ≥ y 1 must be maximized first. Then, the surface level y 2 must be maximized and so on. Let us now focus on input/output. Of course, using directly (1) as an input constraint makes no sense, since ( 1) is equivalent to f (x) ≥ y N . One could rather solve N times the set inversion problem f (x) ∈ [y k , y k+1 ] and superimpose the different outputs (see Figure 4.b). However, this is not satisfactory for at least two reasons.

First, the boundaries around each level surface (the little white boxes in Figure 4.b) should not appear, except for the lowest one. This undesirable boundaries will be cleared off in Figure 11. One could circumvent this effect by slightly enlarging the intervals [y k , y k+1 ] before computing the inversions, making level surfaces overlap. Another trick would simply be to move the boundary subpaving of the k th surface level into the inner sub-paving of the (k -1) th level surface. However, in both cases, we "cheat" and lose reliability of the result.

Second, the decomposition of the input problem into N sub-problems and the aggregation of the intermediate results requires extra and undesirable manipulations from the user.

Set inclusion

Given two sets defined by constraints c 1 and c 2 , the set inclusion problem consists in proving or refuting that one set is included in the other. This is illustrated in Figure 5.a. We shall denote by set(c) the set associated to a constraint:

set(c) = {x ∈ R n | c(x) is true}.
Using the conjunction of c 1 and c 2 as an input in algorithm set-describer may address the problem in the following situations: if no solution is found, then no set can be included in the other since these sets do not intersect at all. In the other way around, if the whole initial box is included in the inner sub-paving then the two sets coincide. Except in these extreme cases, nothing can be said.

As before, one may solve separately c 1 and c 2 , and check for the inclusion of both the inner & boundary sub-pavings of c 1 into the inner sub-paving of c 2 (see Figure 5.b). Again, this requires undesirable manipulations by the user; but, worse, this Of course, one could legitimately object that in the example given in Figure 5, c 1 can be reversed as in §2.2 to c1 so that running a standard solver with c1 and c 2 will do the job as expected. Indeed, if no solution is found, no point satisfy c 2 and c1 , i.e., set(c 2 ) ⊆ set(c 1 ). However, in general, a constraint cannot be reversed while an inner test is still possible. As an example, a constraint c 1 (x) of the following form:

c 1 (x) ⇐⇒ ∃y ∈ [0, 1] f(x,y)=0 g(x,y)=0
cannot be reversed if y has multiple occurrences while testing for an inner box is possible (an algorithm is given in [START_REF] Goldsztejn | A branch and prune algorithm for the approximation of non-linear AE-solution sets[END_REF]).

Bounded-error Parameter estimation

This last example is rather devoted to solver design: the problem is not new in terms of input/output but requires a specific strategy.

Consider the following model of a real-time system [START_REF] Jaulin | Guaranteed Bounded-Error Parameter Estimation for Nonlinear Models with Uncertain Experimental Factors[END_REF]:

f (p, t) = 20 exp(-p 1 t) -8 exp(-p 2 t)
where p 1 and p 2 are two unknown parameters. Assume that a vector of 10 measurements y 1 , . . . , y 10 is available, these measurements corresponding to times t 1 , . . . , t 10 respectively.

The question is to compute a rigorous enclosure of all the feasible p 1 and p 2 from the experimental data, taking into account uncertainties on both y i and t i . Hence, given intervals [t i ] and [y i ], we have to describe the following set: (b) The result using a standard contraction and inner test without splitting the domains of ti (bisection precision is 0.01). First, a very few of inner boxes could be found. Second, the outer enclosure is not sharp at all.

{(p 1 , p 2 ) | ∀i ∈ [1..10], ∃y i ∈ [y i ], ∃t i ∈ [t i ] y i = f (p, t i )}.
If we restrict ourselves to a single measurement (instead of 10), the problem boils down to a simple set inversion, with a thick constant [t i ]. Hence, the set to be computed seems to be nothing but the intersection of 10 set inversions. The set-describer algorithm seems to suit well.

In fact, the multi-incidence of t in the expression of f makes this strategy inefficient. Indeed, whatever is the underlying method used for the contractor or the inner test, the domain of t has to be small for the resulting operator to be sharp. Hence, the domain of t has to be split by the solver. This is illustrated in Figure 6.

Second, introducing the vector t as variables makes the search space dimension jumps from 2 to 12 unless one wants to solve the 10 problems separately. This leads to a combinatorial explosion.

In conclusion, there is a need for splitting the component t i (and only this one) inside the contractor and inner test associated to the i th measure. To our knowledge, no available solver is generic enough to be tuned in this way.

A new formalism

We show now that the different examples of Section 2 can be solved as suggested by the very same algorithm, a generic solver called paver, when expressed in the contractor formalism.

We describe the paver and the formalism in this section. We will provide concrete examples of contractors in Section 4 and subsequently revisit the problems of Section 2.

The "contract and classify" paver

To help intuition, let us consider the problem in §2.3. It turns out from this example that several sets have to be described simultaneously, each set corresponding to a different constraint.

Let us assume that sets are disjoint. We will see below how to circumvent the case of overlapping sets.

Each point of the initial space has to be classified with respect to these sets. For a given point x, this means that the membership of x to the different sets can be tested successively. If one test succeeds, then x can be marked as treated and associated to the corresponding set.

A generic interval solver can directly be derived from this "test and classify" principle. Remember first that a contractor C can be easily built from any constraint c.

Only [x]

′ is left to be treated, and the second contractor is called. When all the contractors become ineffective on a box [x], we can say that the common fixpoint of all contractors is reached and [x] is bisected. The whole process is then repeated until the list of boxes becomes empty. We will see in the next section how the termination of the algorithm can be guaranteed.

algorithm paver(contractors C 1 , . . . , C k , box [x] (0) ) output: a paving of [x] (0) create k empty subpavings K 1 , . . . , K k push [x] (0)
on the stack while the stack is not empty | pop [x] from the stack Note that paver is not an AC3-like algorithm, for essentially two reasons. First, contractors given to the paver algorithm usually involve all the variables merely because they all correspond to sets (or more generally, to sub-pavings) which have the same dimension as the initial box. This is not a rule but a "user advice". One can put different constraints at this level in order to trace their contraction power. But, the paver level is intended to deal with sets while all the constraints related to a given set should be encapsulated in the same outer or inner contractor, as explained further. Second, because the semantics of contractors goes beyond that of constraints, the order in which they are called may have an influence on the quality of the result (the number of boxes in a subpaving). For instance, a contractor designed to detect a very weak condition should preferably be called on last resort (we might lose boxes satisfying a stronger condition).

Contractors

Definition 1 (Contractor). A contractor is a mapping C from IR n to IR n such that

(i) ∀[x] ∈ IR n , C([x]) ⊆ [x] (contraction) (ii) (x ∈ [x], C({x}) = {x}) ⇒ x ∈ C([x]) (consistency) (iii) C({x}) = ∅ ⇔ (∃ε > 0, ∀[x] ⊆ B (x, ε) , C([x]) = ∅) (continuity)
where B (x, ε) is the ball centered on x with radius ε.

A box [x] is said to be insensitive to C if C([x]) = [x]
and sensitive otherwise. By extension, a point x is said to be sensitive or insensitive whether {x} is sensitive or not. Property (i) states that a box can only be reduced by a contractor. Property (ii) states that no insensitive points can be removed. Finally, property (iii) is required for properties on the paver as we will explain further.

The set associated to a contractor C is the union of all of its insensitive points:

set(C) = {x ∈ R n , C(x) = x}.
The continuity property of C implies that set(C) is closed.

If C is a contractor for a constraint c, in the classical sense of the word (see §1.1), then set(C) ⊇ set(c). An important novelty in our formalism is to consider sets associated to contractors rather than to constraints. This allows a rigorous description of the output of our paver (see Proposition 2). A contractor is not only an algorithm, it can also be interpreted as a subset of R n , and all the standard operations on sets can be extended to contractors. We define:

(C 1 ∩ C 2 )([x]) := C 1 ([x]) ∩ C 2 ([x]) (intersection) (C 1 ∪ C 2 )([x]) := C 1 ([x]) ∪ C 2 ([x]) (union) (C 1 • C 2 )([x]) := C 1 (C 2 ([x])) (composition) C ∞ 1 := C 1 • C 1 • C 1 • . . . (iterated composition) C 1 ⊓ C 2 := (C 1 ∩ C 2 ) ∞ (iterated intersection) C 1 ⊔ C 2 := (C 1 ∪ C 2 ) ∞
(iterated union).

(

) 2
All these operations are stable, i.e., they only yield contractors.

We also introduce the following definition:

Definition 2. Let C 1 , . . . , C k be a collection of contractors. C 1 , . . . , C k are complementary if set(C 1 ) ∩ . . . ∩ set(C k ) = ∅, C 1 , . . . , C k are independent if ∀i = j, set(C i ) C ∩ set(C j ) C = ∅.
Additional properties of contractors play a significant role: This result being rather intuitive, the proof is not given (see Figures 7 and8). In practice, these operators are implemented as follows. For the proj-union, [y] is split into a sequence of small subintervals [y] i . Then, the contractor C is enforced on each sub-domain [x] × [y i ] and the hull of the results is returned. For the proj-intersection, a sequence of points y i are sampled from [y], C is enforce on [x] × {y i } and the intersection of the results is returned.

C is monotonous if [x] ⊆ [y] ⇒ C([x]) ⊆ C([y]), C is minimal if ∀[x] ∈ IR n , C([x]) = ([x] ∩ set(C)), C is idempotent if ∀[x] ∈ IR n , C(C([x])) = C([x]).
Although examples of contractors will be the topic of Section 4, we shall introduce right here the precision contractor C ε that has a special status. This contractor is aimed at controlling the precision of the paver or scheduling contractors (see §5.1 and §5.3) and will be implicitly referred to in the subsequent propositions. ) is obtained.

Definition 3 (Precision contractor). Given ǫ > 0, the precision contractor C ε is defined as follows:

C ε ([x]) := [x] if width([x]) > ǫ ∅ otherwise.
Many properties could be stated on contractors. As an example, the set of idempotent and monotonous contractors is a complete lattice for the inclusion. Listing these properties would be out of the scope of this paper, which is rather aimed at introducing the concepts. However, the fundamental properties of the paver given in the previous section are based on the properties of contractors and we are now in position for proving them.

Paver properties

Proposition 1 (Termination). Let {C 1 , . . . , C k } be a list of complementary contractors. The paver terminates for any initial box [x] (0) .

Proof. Assume by contradiction that paver does not terminate. Hence, there is a sequence of non-empty boxes [x] (i) and x ∈ [x] (0) such that for all i > 0:

(i) [x] (i) ⊂ [x] (i-1) (strict inclusion, i.e., [x] (i) = [x] (i+1) ) (ii) C 1 ([x] (i) ) = . . . = C k ([x] (i) ) = [x] (i) (iii) x ∈ [x] (i)
Indeed, (i) is a direct consequence of the properties of contractors and bisections and (ii) is due to the fixpoint postcondition of the inner while loop inside the paver algorithm. Finally, for all i, pick a point x(i) ∈ [x] (i) (which is nonempty). The sequence is bounded so it admits at least one accumulation point x. Since for all i, [x] (i) is a closed set that contains the whole subsequence (x (j) ) j≥i , it also contains this accumulation point x. Hence, (iii) holds. Now, each box [x] (i) can be identified to a 2n-tuple u (i) of its vertices coordinates. This identification can be made isometric with appropriate distances (the Hausdorff distance on IR n and its counterpart on R 2n ). For all component j, 1 ≤ j ≤ 2n, the sequence u (i) j is either increasing or decreasing, and bounded by a component of x. Hence, it converges to a point ūj . By the inverse isometry, it follows that the sequence of boxes [x] (i) converges then to a box [ x]. Furthermore, the width of this box is necessarily null (i.e., [ x] is a degenerated box {x})

since ∀i ≥ 0, rad([x] (i+n) ) ≤ 1 2 rad([x] (i) ).
Next, we have set(C 1 )∩. . .∩set(C k ) = ∅ which means that at least one contractor C l satisfies x ∈ set(C l ). By applying the continuity property of C l (see Definition 1) we obtain that there exists ǫ and i such that [x] (i) ⊆ B(x, ǫ) and C l ([x] (i) ) = ∅, i.e., a contradiction.

To ensure the complementarity between contractors, we usually resort to the precision contractor because set(C ε ) = ∅.

Informally, the following proposition states that, up to a given precision ǫ, the sub-paving associated to a contractor matches the intersection of the other contractors sets. Since the set of a contractor usually approximates the set related to an initial constraint, we have now a clear semantics for the paver.

Proposition 2 (Sub-pavings characterization). Let {C 1 , . . . , C k } be a list of independent contractors and C ε such that set(C ε ) = ∅. Let us denote by K 1 , . . . , K k and K ǫ the sub-paving returned by paver for C 1 , . . . , C k and C ε respectively. We have:

∀i, 1 ≤ i ≤ k, j =i set(C j )\ K i ⊆ K ǫ .
Proof. First of all, paver terminates thanks to the previous proposition. Let i be an index, 1 ≤ i ≤ k. For all j = i we have by hypothesis set

(C i ) C ∩ set(C j ) C = ∅ hence set(C i ) C ⊆ set(C j ) which implies set(C i ) C ⊆ ∩ j =i set(C j ). Since ∀[x] ∈ K i we have [x] ⊆ set(C i ) c (by applying the consistency of C i ), then K i ⊂ set(C i ) C ⊂ ∩ j =i set(C j ).
C l with 1 ≤ l ≤ k, i.e. [x] ⊆ set(C l ) C . Since [x] ∈ j =i set(C j ) then l = i which leads to a contradiction since [x] ∈ K i .

Contractor design

In this section, we present some contractors and operations that do not help in modeling problems but in improving the efficiency of a paver. We basically explain how to manage a constraint c, i.e., the different ways to build a contractor C satisfying set(C) = set(c). Likewise, all the operations op will preserve associated sets. E.g., an unary operator op will satisfy for all contractor C, set(op(C)) = set(C).

Hence, contractors and operations in this section will not have an influence on the result, in terms of sub-paving characterization. On the other hand, they may decrease the time complexity or the space complexity (i.e., the number of boxes in a sub-paving).

We will overlook technical details on purpose here, sticking to generalities. The interested reader may refer to Section 6 where important implementation choices we made for our own system are disclosed.

We shall distinguish numerical contractors, those which are directly built from a numerical constraint, and those which are the result of an operation (such as composition) between other contractors.

Numerical contractors can be related to equations, inequalities or systems of equations. The point is that a large part of interval analysis routines can be wrapped into these contractors. As announced above, a good separation can then be obtained between numerical skills (the design of numerical contractors) and constraint programming skills (the design of compound contractors).

In this layered framework, constraint propagation will be generalized to contractor propagation. A possible implementation of the corresponding operator (propag) will be described in detail in §6.3.

Numerical contractors

Several contractors can be associated to an equation f (x) = 0.

The simplest one consists in evaluating with interval arithmetics f ([x]) and checking whether 0 belongs to the image range or not. If 0 ∈ f ([x]) then [x] can be contracted to the empty set. Otherwise, it is left unchanged. This contractor can therefore be qualified as binary (in the sense that it keeps all or nothing). This test can be easily extended to inequalities. This contractor can be given many variants, based on various symbolic or numerical processing as the two following examples show.

Example 3. Assume that there exist several equivalent expressions for f , say f 1 , . . . , f k , each of them minimizing the multi-incidence of a different variable.

Since the overestimation of interval evaluation grows with the multi-incidence of variables (see [START_REF] Moore | Interval Analysis[END_REF]), we get to sharper results by computing f

1 ([x]) ∩ . . . ∩ f k ([x]) instead of f ([x]
), whence a more accurate (but slower) test.

Example 4. On the numerical side, we can resort to the centered form, i.e., compute f

(x 0 ) + f ′ ([x])([x] -x 0 ) instead of f ([x]
) where x 0 is an arbitrary point inside [x]. Higher-order Taylor formulas may also be payfull. One may also introduce intermediate monotonicity tests. Consider a two-dimensional vector x.

If ∂f ∂x1 ([x]) turns to be a positive interval then f is increasing with respect to ator of box consistency [START_REF] Van Hentenryck | Numerica: A Modeling Language for Global Optimization[END_REF]). 3. Linear relaxations (see [START_REF] Lebbah | Efficient Pruning Technique Based on Linear Relaxations[END_REF] and references therein). In such techniques, a nonlinear constraint is cast into a linear program which feasible region encompasses the original solution set.

x 1 so that f ([x]) can be replaced by [f ({x - 1 } × [x 2 ]), f ({x + 1 } × [x 2 ])],
Given an equation f (x) = 0, we will denote by C f in the sequel an arbitrary contractor among those above.

Of course, numerical techniques can also handle several equations simultaneously, i.e., a system of equations. Many of them are derived from a multivariate interval Newton, one famous variant being the Hansen-Sengupta algorithm [START_REF] Hansen | Bounding Solutions of Systems of Equations Using Interval Analysis[END_REF].

In case of linear equalities, many dedicated algorithms also exist (see [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF] or [START_REF] Ning | A Comparison of Some Methods for Solving Linear Interval Equations[END_REF]). One can be easily convinced that all these algorithms act as contractors.

Propagation

The propagation operator allows the implementation of interval variants of the classical AC3 algorithm such as hull consistency (see, e.g., [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF][START_REF] Delobel | Comparing Partial Consistencies[END_REF]) but also many more algorithms. This operator illustrates by itself the potential of our framework, in terms of contractor design. The key idea is to propagate contractors instead of constraints.

Given a list of contractors, the principle is to obtain the fixpoint of their composition, i.e., (

C 1 • • • • • C m )
∞ at a lower computation price. We shall denote by propag our operator. We have:

set(propag(C 1 , . . . , C m )) = set((C 1 • • • • • C m ) ∞ )
which means that propag is a pure efficiency contractor (it has no impact on the output).

An implementation of this operator will be proposed in §6.3.

Let us now recall that the hull consistency enforces bound consistency for each equation. In practice, only a relaxation of bound-consistency is calculated via the HC4Revise algorithm. The resulting propagation algorithm is called HC4. In our framework, given a list of constraints c 1 , . . . , c m , the HC4 algorithm is simply reimplemented as follows HC4(c 1 , . . . , c m ) ← propag(HC4Revise(c 1 ), . . . , HC4Revise(c m ))

A more interesting direct application is a multi-level propagation. Sometimes, we would like to group constraints we know to have strong "dependencies". For instance, one may find relevant to perform an intermediate fixpoint with constraints that strictly share the same variables. Consider three sets of constraints {c 1 , c 2 , c 3 }, {c 4 } and {c 5 , c 6 }. The following contractor will ensure that a constraint in a subset (say, c 4 ) will not be wakened between constraints of other subsets (say, c 1 and c 2 ).

propag(propag(c 1 , c 2 , c 3 ),c 4 ,propag(c 5 , c 6 ))

Shaving

The purpose of this section is to show that other operators exist besides propag.

Shaving is an operation that allows implementing refutation techniques, similar to SAC [START_REF] Bessière | Optimal and Suboptimal Singleton Arc Consistency Algorithms[END_REF] with discrete domains. With continuous domains, refutation is used to shrink endpoints only (instead of any value inside the domains) so that the structure of interval is always maintained (whence the fancy name shaving).

Detailing the algorithm here would take too much room. We shall only give a rough description.

A shaving operator shave takes a contractor C. Given a box [x], the resulting contractor shave(C) contracts "slices" with C, i.e., subboxes obtained from [x] by restricting the domain

[x] k of one component to a small subinterval [x - k , x - k +ǫ] or [x + k -ǫ, x + k ].
When the result of the subcontraction is an empty set, the slice is removed. Otherwise, contraction is tried on a smaller slice. This recursion leads to consistent endpoints. More formally, the resulting box

[y] := (shave(C))([x]) satisfies ∀k = 1, . . . , n C([y] 1 × • • • × [y] k-1 × {y - k } × [y] k+1 × • • • × [y] n ) = ∅ and C([y] 1 × • • • × [y] k-1 × {y + k } × [y] k+1 × • • • × [y] n ) = ∅.
The shaving operator can also be given many variants mainly because of the slicing which can either be optimistic ("try large slices first") or pessimistic ("check consistency of endpoints first").

Given m constraints c 1 , . . . , c m , a box consistency operator (see [START_REF] Benhamou | CLP(intervals) revisited[END_REF][START_REF] Van Hentenryck | Solving Polynomial Systems Using a Branch and Prune Approach[END_REF][START_REF] Van Hentenryck | Numerica: A Modeling Language for Global Optimization[END_REF]) can be reimplemented in our framework as follows:

propag(shave(univNewton(c 1 )),. . . ,shave(UnivNewton(c m )))

where UnivNewton is an univariate interval Newton contractor (see §4.1). The 3B consistency defined in [START_REF] Lhomme | Consistency Techniques for Numeric CSPs[END_REF] can be reimplemented as follows:

shave(HC4(c 1 , . . . , c m )).

Revisiting examples

We now revisit all the examples in our framework. By convention, the sub-paving of the precision contractor will always be represented in white in the pictures. As a preliminary example, consider the system of equations in Figure 1 page 8 and remember that the goal is to enclose the solution set into a sub-paving. By a direct application of our formalism, we then need two contractors:

1. A contractor C f , to remove unfeasible points.

2. The precision contractor C ε .

The desired output is given in Figure 9. This time, the result of the outer contraction appears as a sub-paving with only unfeasible points. The roots are all in the sub-paving of the precision contractor.

A decomposable problem

Assume that each subsystem (called block) has to be solved by a combination of propagation and interval Newton. Of course, a multivariate interval Newton (denoted by newton henceforth) requires a square system. When one variable is fixed (e.g., x C and y C when the second triangle is to be solved), it must not be considered as a variable by newton but rather as an interval constant. Such information typically enters as a "static" parameter (see §6.1). For that reason, the newton contractor takes two arguments: the set of constraints and the set of variables. Similarly, let us refine the behaviour of the precision contractor C ε (see §4). This contractor will also take here a subset S of variables in parameter.

A box [x] is emptied by C ε (S) if and only if the domains of the components in S all have a diameter lower than ε.

In the sequel, d P Q denotes the distance constraint between two points P and Q.

For each block (here: a triangle), the solving strategy is implemented thanks to the following operator block contractor:

block contractor({c 1 , . . . , c k }, {x 1 , . . . , x k }) ← HC4(c 1 ,. . . ,c k ) ∩ newton({c 1 , . . . , c k }, {x 1 , . . . , x k })
For instance, a paver run with block contractor({d BD , d CD }, {x D , y D }) (and C ε ) with x C and y C fixed (i.e., their domains have a diameter less than ε) will solve the second triangle in the required manner, i.e., using propagation and Newton.

The step-by-step solving method proposed in §2.1 can then be directly implemented with the following contractor: provided that variables are bisected block after block. In more details, the variables of the first block {x C , y C } must be bisected (say, in a round-robin fashion) until precision is achieved for both of them. The bisection procedure must then proceed to the second block, and so on. Branching procedures are not covered in this paper but this mechanism is orthogonal to contractors and one should be easily convinced that a block-wise bisection is easy to set up.

Enhanced description of the ring

In the set inversion problem, we need to characterize {x | f (x) ∈ [y]}. This means that two sets have to be described by contractors: {x | f (x) ∈ [y]} (the ring) To deal with unfeasible points, we need an outer contractor. This contractor, denoted by C f ∈[y] , can be simply defined as follows:

C f ∈[y] := [x] → π x (C g ([x], [y])).
The points inside the ring have to be classified with an inner contractor. As explained in §2.2, the latter comes naturally with the negation of the problem. We have

c(x) ⇐⇒ f (x) ∈ [y] ⇐⇒ y -≤ f (x) and f (x) ≤ y + ⇐⇒ f (x) ∈ [y -, +∞) and f (x) ∈ (-∞, y + ]. Hence, c(x) ⇐⇒ f (x) ∈ (-∞, y -) or f (x) ∈ (y + , +∞) =⇒ f (x) ∈ (-∞, y -] or f (x) ∈ [y + , +∞)
Therefore, the contractor

C f ∈[y] ([x]) := [x] → π x (C g ([x], (-∞, y -])) ∪ π x (C g ([x], [y + , +∞)))
only removes unfeasible points for c, i.e., feasible points for c.

Now, set(C f ∈[y] ) ∩ set(C f ∈[y]
) is the boundary of the ring, which is non empty. Hence, we have to add a precision contractor to ensure the termination of the algorithm. The result is depicted in Figure 10. Since C f ∈[y] and C f ∈[y] are independent, the ring represented on the figure can be interpreted as an approximation of the set to be described (see proposition 2).

Level surfaces

First of all, the contractor C f ∈[yi,+∞) (as defined in the previous paragraph) could be associated to the i th level surface {f ≥ y i }. Of course, in this way, the contractor of the lowest surface level would preempt most of the boxes, leading to a bad result (independence does not hold).

Introducing a priority between the different contractors (in addition to independence) does not require any extra concept. We just need to adapt the contractors to their actual semantics: a (sufficiently small) box [x] must be classified with the contractor of the i th surface level only if

[x] is inside the i th surface level and [x] does not intersect the (i-1) th surface level or [x] has a small width (lower than ε)

The desired contractor C i is then obtained by simply rewriting these conditions in terms of sub-contractors:

C i := C f ∈[yi,+∞) ∪ C f ∈[yi-1,+∞) ∩ C ε
or, using a sub-distributivity rule:

C i := C f ∈[yi-1,+∞] ∩ C f ∈[yi,+∞) ∪ C ε .
Finally, we need as before an outer contractor for all the surface levels, namely,

C n+1 := C f ∈(-∞,yn] . Moreover, the intersection of C 1 ∩ . . . ∩ C n+1 is {x | f (x) =
y n } and, as usual, this set can be treated by a precision contractor C ε . Figure 11 shows the desired output. The set inclusion problem can be solved efficiently in our formalism by using the following simple reasoning.

Set inclusion

A box [x] can be discarded from the search either if it does not belong to set(c 1 ) or if it is included in set(c 2 ). Indeed, in both cases, [x] cannot compromise the assertion to be proven since ∀x ∈ [x] we have c 1 (x) =⇒ c 2 (x).

Hence, an inner contractor C 1 for c 1 and an outer contractor C 2 for c 2 can be used jointly. Only the box that both (possibly) contain points inside set(c 1 ) and outside set(c 2 ) are bisected. The result is depicted in Figure 12.

No precision need to be introduced if int(set(c 1 )) ⊆ set(c 2 ) (where int stands for interior) since the algorithm will end in this case (we have set(C 1 )∩set(C 2 ) = ∅).

But we still need a precision contractor C ε in case of non inclusion.

Bounded-Error Parameter estimation

As we evoked above, modeling this problem causes no difficulty. The overall inner contractor is the union of the inner contractors related to each measure, say,

C gi ∈[yi] with g i (p) := f (p, t i ) and t i ∈ [t i ].
Here, the interval extension of g i is a thick function since an interval constant [t i ] substitutes for the variable t i .

Likewise, the intersection of the outer contractors C gi∈[yi] defines an outer contractor for the whole problem. For efficiency, we have to split t i inside the i th outer and inner contractor. Consider first outer contraction:

(p 1 , p 2 ) is feasible if at least one value in [t i ] satisfies f (p, t i ) ∈ [y i ]. Hence, if a point (p, t i ) ∈ R 3 is insensitive to C f ∈[yi]
then p must be insensitive to the outer contractor. The proj-union defined in §4 exactly fits this requirement. We define:

C outer := C ∪[t1] f (p,t1)∈[y1] ∩ . . . ∩ C ∪[t10] f (p,t10)∈[y10]
Similarly, the proj-intersection fits the requirement for the inner contraction. Indeed, a point (p, t i ) ∈ R 3 sensitive to C f ∈[yi] must be discarded by the inner contractor. We define:

C inner := C ∩[t1] f (p,t1) ∈[y1] ∪ . . . ∪ C ∩[t10] f (p,t10)) ∈[y10]
A paving resulting from the combination of these contractors is shown in Figure 13.

The Quimper System

Besides theoretical investigations, contractor programming has given rise to a real system named Quimper (QUick Interval Modeling and Programming in a bounded-ERror context). This system includes today three different programs: qPave (a graphical tool for paving sets in 2D), qSolve (a tool for listing numerical results, typically in high dimension) and qTraj (a graphical tool tailored to constraint problems derived from differential equations).

These are light-weight programs that only manage input/ouput, i.e., the language for writing contractors (the "Quimper language") and the interface for configuring, running the paver and handling generated pavings. They are linked to a C++ library called Ibex (Interval-Based Explorer) that implements a contractor programming framework.

Ibex is based itself on the Profil/Bias library [START_REF] Knüppel | Profil/Bias[END_REF] for the low-level interval arithmetics. However, part of this library has been wrapped into functions that manage all the borderline cases (infinite bounds, values out of definition domains, empty sets, etc.) so that arithmetic operations are always exception-free. This explains why a modified release of Profil/Bias is included in the Ibex package.

All these software components are under GPL licence and can be downloaded online [START_REF] Chabert | IBEX, an Interval-Based EXplorer[END_REF]. A user guide for Quimper (including the grammar of the Quimper language), as well as a complete documentation of Ibex classes are available on the same site. An archive containing all the examples of this paper in the Quimper syntax is also provided.

Ibex/Quimper has been developed by the first author but this software should be considered as a prototype since a real collaborative open-source project continuing this work is about to be launched. This new project will take advantage of the existing code.

The purpose of this section is to give some insight into the Ibex system. For convenience, we shall adopt some object-oriented coding notations. A contractor C is a class with the main function being contract, i.e., for a box

[x] C.contract([x]) contracts [x] with C.
Generalizing constraint propagation to contractor propagation required letting the interface of contractors inherit from constraints. First, we can ask a contractor whether the domain of a given variable can impact the result of the contraction or not. This notion simply generalizes the incidence graph of a constraint network.

A contractor C, as a class, therefore implements a Boolean function involves that takes in argument the index j of a variable, i.e.,

C.involves(j)

returns false or true. Note that, by default, a contractor always returns true.

In the parlance of object-oriented programming, this function needs not necessarily be overridden.

Second, we also added a few parameters to the contract function above, besides the box to be contracted, as explained in the next section.

Indicators

Every contractor has its own set of specific parameters. For instance, the precision contractor C ε (maxdiamGT in the Quimper syntax) takes a parameter ε. The interval Newton also needs a parameter for controlling the termination of the iteration: when a step does not reduce any interval by more than a given ratio, the procedure stops. This ratio can be set externally. These parameters are related to the semantics specific to contractors (e.g., the precision contractor has by definition a threshold) or their implementations. They are usually set once for all (as constructor parameters). We shall not consider this type of parameters any further.

There is another type of parameters required by constraint programming algorithms for efficiency reasons, called indicators. The purpose of these parameters is to notify contractors during the search about the context in which they are called. For instance, we may inform a contractor that only a contraction on a given variable is actually expected (any other contraction being superfluous).

If the time complexity of the contractor depends on the number of variables, some work is spared. We may also inform the contractor that only the domain of a given variable has been modified since the last call. Again, if the contractor works incrementally, this will speed up contraction.

Two indicators have been integrated into our system, corresponding precisely to the examples just given. The first is named scope and contains the subset of variables to be treated. The second is named impacted and contains the subset of variables whose domain has been impacted.

The semantics of indicators is constrained by one single fundamental property: they can be ignored by a called contractor without spoiling soundness (i.e., losing solutions), whence their name. If one develops a contractor that ignores the scope indicator, this contractor can still be passed as argument to an operator, should the latter be based on the communicability of this information (e.g., propagate below). The consequence is only a loss of efficiency.

Hence, only two indicators are proposed today but the nice point with indicators is that new ones can be invented at any time and integrated progressively in the subsequent implementations of contractors. In other words, backward compatibility is complete. We just had to put all the indicators in a dictionary structure. The real signature of the contract function is:

contract(box [x]
, dictionary indicators).

As an example, one can ask a contractor C to focus on the two first variables only in the following way:

C.contract([x], {scope = {0, 1}}).

Numerical Contractors

There is nothing particular in the implementation of numerical contractors in Quimper, except that indicators have to be managed.

-The binary test 0 ∈ f ([x]) is actually not proposed in Quimper, but such a contractor would necessarily ignore all the indicators. -HC4Revise is implemented. The complexity of this contractor is linear in the length of the constraint expression. Hence, contracting for a single variable amounts to contracting for all the variables (up to a small constant factor).

In Quimper, HC4Revise ignores the scope indicator. -The univariate interval Newton is implemented. As contrary to HC4Revise, this contractor deals with one particular variable at a time. Setting the scope indicator to a single variable divides the contraction time by the number of variables. -Linear relaxations are not yet implemented. Now, since a linear programming solver is called iteratively to reduce the bounds of each initial variable, the scope indicator would have to be taken into account as with univariate Newton.

The multivariate interval Newton is also implemented and ignores the scope and impacted indicators.

Propagation

Here is how the propagation operator is implemented in Quimper.

Our convention for indices is to use i and j for contractors and variables respectively. When a couple (i, j) is revised, the i th contractor has to work on the j th variable. If the revision succeeds (a significant part of the i th domain is removed) then the agenda is updated with the following classical procedure:

update agenda(integer i, integer j, subset scope) | forall i ′ = i such that C i ′ .involves(j) | | forall j ′ such that C i ′ .involves(j ′ ) | | | if j ′ ∈ scope\{j} then add (i ′ , j ′ ) in the agenda Now, if a revision (i, j) fails, there still may be a residual contraction. Furthermore, some subcontractors may not take into account the scope indicator so that unsolicited contractions can appear as the propagation loop goes along. At some point, the accumulation of these small contractions on a variable j can add up to a significant contraction. The agenda must be updated in consequence. In this case, the j th variable is the source but there is not a particular contractor. The agenda must then be updated with all the couples (i, j), i describing the set of contractors. This is what the next procedure does.

in this case, mainly because of the parameter splitting process which occurs systematically. The paving includes 18249 boxes.

Conclusion

We have presented a new framework for interval programming. The benefits of this framework are twofold. First, a large class of constraint-based problems can now be addressed with a unique simple algorithm, called paver. Second, a full control on the solving process is at hand, including the core of the propagation loop.

This framework is entirely built on the concept of contractor: inputs are contractors instead of constraints. Similarly, outputs are related to contractor sets instead of constraints sets. As just said, this -apparently small-change in modeling can lead to significant improvements in the design of new solvers but also in a declarative way. In a sense, the imperative aspect of solvers is now subsumed in the declarative one: the end user write constraints while the constraint programmer write contractors.

This approach is supported by a real system called Quimper which solves quickly all the different problems mentioned in this paper.

There are plenty of possible extensions for this work. The first one is perhaps to deal with global optimization. This extension would probably resort to dynamic contractors, i.e., contractors parameterized by a value that can be updated during the solving process. Orthogonally to contractors that potentially remove all unfeasible points, local finders that find peculiar feasible points seem to be a key feature as well.

Fig. 1 .

 1 Fig. 1. An enclosure of {x ∈ [-0.5, 2.5]×[-2.5, 0.5], f (x) = 0} with f (x) = exp(x1x2)sin(x1 -x2).

Fig. 2 .

 2 Fig. 2. A set inversion problem f (x) ∈ [y] with f (x) = p x 2 1 + x 2 2 and [y] = [3, 6] in [x] = [-7, 7] × [-7, 7]. (a) The output of the set inversion problem as defined in §2.1 (i.e., small boxes enclosing all the solutions). (b) A more adapted output for the same set inversion problem, including inner test.

Figure 4 .

 4 a shows an example of 4 surface levels. Note that this problem could be easily generalized to a MAX-SAT problem.

Fig. 4 .

 4 Fig. 4. Surfaces levels f (x) ≥ yi. (a) An example with y1 = 0.8, y2 = 0.6, y3 = 0.4, y4 = 0.2 and f (x) = sin(x 2 1 + x 2 2 )/(exp(x1) + x 2 2 ), in [x] = [-2, 2] × [-2, 2]. Darkness increases with f . (b) The same example decomposed into set inversion problems.

Fig. 5 .

 5 Fig. 5. Set inclusion. (a). Two subsets of [x] = [-3, 1]×[-1, 3] defined by constraints: c1(x) ⇐⇒ exp(x1x2) -sin(x2 -x1) ∈ [-0.1, 0.1] (in dark gray) and c2(x) ⇐⇒ exp(x1 -x2) × sin(x1x2) ∈ [-0.1, 0.1] (in light gray). We have set(c1) ⊆ set(c2). (b) The same set inclusion problem decomposed into set inversion problems.

Fig. 6 .

 6 Fig. 6. The parameter estimation problem. (a) The set of feasible points (p1, p2).(b) The result using a standard contraction and inner test without splitting the domains of ti (bisection precision is 0.01). First, a very few of inner boxes could be found. Second, the outer enclosure is not sharp at all.

Example 2 .

 2 The contractor C has the following property: if C([x]) = [y] (i.e.,[y] is the contraction of [x]) then every eliminated point in [y]\[x] does not satisfy c. Thus, with a simple negation trick, a contractor C such that all eliminated points satisfy c is as easy to build. Consider the set inversion problem, in §2.2. A point inside the ring satisfiesf (x) ∈ [y]. Consider a contractor C built from f (x) ∈ [y]. If C([x]) = [y] then [x]\[y] is inside the ring.Given a box [x] and a constraint c, if the complementary contractor C reduces [x] to [y], then [x]\[y] can be classified, i.e., mark as "inside c". Then, the paver proceeds as follows (the algorithm is detailed below). The solver takes as input a list of contractors and an initial box [x] (0) . The first contractor is enforced on [x] := [x] (0) . If the contraction is effective, a smaller box [x] ′ is returned and the difference4 [x]\[x] ′ is stored in a sub-paving associated to the first contractor.

( 3 )

 3 Given a contractor C on R n × R m , we also define the two following contractors on R n :C ∪[y] ([x]) := y∈[y] π x (C([x], y) C ∩[y] ([x]) := y∈[y] π x (C([x],y) with π x ([x], [y]) := [x]. The first operation is called proj-union and the second proj-intersection. Once again, these operations are stable in the set of contractors. We have: set(C ∪[y] ) = {x, ∃y ∈ [y], (x, y) ∈ set(C)}, set(C ∩[y] ) = {x, ∀y ∈ [y], (x, y) ∈ set(C)}.

Fig. 7 .

 7 Fig. 7. Proj-Union. We consider a contractor C on IR 2 with S := set(C). (a) Representation of C ∪ [y] when [y] is degenerated (a real y). The contraction on [x] results in an enclosure of {x ∈ [x] | (x, y) ∈ S}. Since C is continuous, when both [x] and [y] are degenerated, the enclosure is minimal, i.e., C ∪y (x) = {x} iff (x, y) ∈ S. (b) Principle of an implementation resorting to parameter splitting and union of sub-contractions (the interval [y] is decomposed as finely as possible). (c) Result of the proj-intersection in practice: only an outer approximation of set(C ∪[y] ) is obtained, after numerous repetitions of the splitting operation presented in (b).

Fig. 8 .

 8 Fig. 8. Proj-intersection. (a) In a similar way as in Figure 8, given a fixed interval [y], C ∩[y] allows to constrain x ∈ R with ∀y ∈ [y], (x, y) ∈ S. (b) Principle of an implementation resorting to parameter splitting and intersection. (c) Result of the proj-intersection in practice: only an outer approximation of set(C ∩[y] ) is obtained.

  a smaller interval. More sophisticated contractors associated to a constraint c (equation or inequality) compute an outer estimation of the feasible set, i.e., C([x]) is a nontrivial superset of {x ∈ [x] | c(x)}. Three important techniques can be found in the literature: 1. Forward-backward traversals of the syntax tree with intermediate interval computations, such as HC4Revise [3]. 2. Univariate interval Newton iterations (as performed by the narrowing oper-

Fig. 9 .

 9 Fig. 9. System of equations, with an outer and precision contractor. The sub-paving of the outer contractor is in light gray.

  block contractor({d AC , d BC },{x C , y C }) ∩ ( C ε ({x C , y C }) ∪ ( block contractor({d BD , d CD }, {x D , y D }) ∩ ( C ε ({x C , y C , x D , y D }) ∪ block contractor({d DE , d CE }, {x E , y E })))

Fig. 10 .

 10 Fig. 10. A description of the ring involving outer and inner contractors.

  and {x | f (x) ∈ int([y])} (the complementary of the ring). Put [y] = [y -, y + ] and g(x) := f (x)y.

Fig. 11 .

 11 Fig. 11. An appropriate paving, with smooth boundaries between surface levels.

Fig. 12 .

 12 Fig.12. The set inclusion problem solved by a more adapted algorithm. Boxes in dark gray does not belong to set(c1). Boxes in light gray are inside set(c2). We resort to splitting only in regions where boundaries of set(c1) and set(c2) are close.

Fig. 13 .

 13 Fig. 13. Parameter estimation solved with inner and outer contractors based on projunion and proj-intersection. The precision of the paver is still 0.01.

  The two last inclusions are strict because set(C i ) C is open while the other sets are closed. It follows that j =i set(C j )\ K i is nonempty.

Take [x] in this set and assume [x] ∈ K ǫ . The box [x] has been removed from the search by

This limitation holds for numerical systems involving analytic expressions, which cover most of the mathematical models of physics problems. But modeling languages as such are not limited with real variables any more than with discrete ones (one may introduce table of constraints, piecewise constraints, etc.).

There is still a notable exception: geometrical constraints. Geometry represents a semantic level above algebra; as an example, the "intersection of three spheres" can be introduced as a global constraints instead of three equivalent distance equations[START_REF] Batnini | Décomposition Sémantique pour la Résolution de Systèmes d'Equations de Distance[END_REF]. But, except in such cases, no improvement has to be expected from the modeling side.

We can say that no modeling language dedicated to the output exists so far, and this is another important distinction with discrete problems where this modeling aspect is not as ubiquitous (solvers using explanations are counter-examples).

Note that resorting to a simple tree structure avoids to explicitly describe the difference between [x] and [x] ′ : the latter simply becomes a subnode of the former.

update agenda(integer j, subset scope)

| forall i such that C i .involves(j) | | forall j ′ such that C i .involves(j ′ ) | | | if j ′ ∈ scope then add (i, j ′ ) in the agenda Now, let us give the main procedure. Two local boxes are defined, the first is used to measure the result of each revision (the contraction of C i over x j ) while the latter measure the contractions collected since the last time x j triggered an update. // so that all may have been impacted

contract(box

Some performance results

All the examples in this paper have been executed almost instantaneously (less than 0.1 second) when the precision was not high (this corresponds to the figures with apparent boxes).


J Exinger

C Hartard

F Lafferrière

C Fenninger

L J Charbonnière

author. Dr Helene Jeulin
email: h.jeulin@chru-nancy.fr

Addresses

Evaluation of a Lateral Flow Immunoassay COVIDTECH® SARS-CoV-2 IgM/IgG Antibody rapid test

Keywords: SARS-CoV-2, COVID-19, serological diagnosis, Flow lateral immunoassay Highlights

 Lateral flow immunoassay permits rapid detection of anti-SARS-CoV-2 IgM and IgG  The specificity reached 100% for IgM and IgG detection  Sensitivity reached 100% after day+14 post-symptoms, for IgM/IgG combined interpretation Abstract (150 words) The new epidemic coronavirus SARS-CoV-2 is responsible for severe respiratory illness (i.e. COVID-19). RT-PCR is the gold standard in COVID-19 diagnosis and is performed on respiratory samples, while serological tests may contribute to diagnosis and permit seroprevalence studies. The lateral flow immunoassay (LFIA) PolyDTech simultaneously detects anti-SARS-CoV-2 IgM and IgG using a Sprotein recombinant antigen and has been independently evaluated in two laboratories. The specificity evaluated on 53 pre-pandemic samples reached 100% for IgM and IgG. Analyzing samples from patients with RT-PCR-confirmed infection, IgM/IgG antibodies were detected in 18/26 (69%)

and 58/58 (100%) samples before day 13 th and from the 14 th day post-symptom onset respectively. These results were concordant with those of another LFIA method at 89% success rates (BioSynex BSS IgM/IgG) and at 94% with a chemiluminescent Immunoassay (Maglumi Snibe). Overall, this LFIA method detecting anti-SARS-CoV-2 IgM and IgG is suitable for SARS-CoV-2 serological diagnosis, when the patient is, as recommended, > day+14 after onset of symptoms.

Introduction

In December 2019, a new coronavirus epidemic, caused by SARS-CoV-2 (Severe Acute Respiratory Syndrome Coronavirus-2) started in Wuhan, China. This new pathogen responsible for COVID-19 has rapidly spread across the globe and resulted in a pandemic affecting more than 53 million people (confirmed cases) and leading to 1,300,576 confirmed deaths on the 13 November 2020 (A ACTUALISER au dernier moment (bulletin hebdomadaire) et penser à ajouter la REF). France has been impacted significantly by the « first wave », with a mean of 4,500 daily confirmed cases during the last week of March 2020 and up to 950 daily confirmed deaths during the first fortnight of April 2020 (https://geodes.santepubliquefrance.fr/). The analysis of the dynamics of mortality from COVID-19 in French hospitals showed that the Grand-Est region was one of the most affected in France (beginning week 9, peak during the week 13), thus is respresents an ideal locale for the investigation of the effects of the aforementioned virus Coronaviruses belong to the Coronaviridae family and are classified into four genera including α-, β-, γ-, and δ-coronaviruses [START_REF] Kirtipal | From SARS to SARS-CoV-2, insights on structure, pathogenicity and immunity aspects of pandemic human coronaviruses[END_REF]. They consist of 60-220 nm diameter enveloped particles presenting the largest known RNA viral genome (26-32 kb in length).

Coronaviruses primarily affects the respiratory system; among the 7 pathogenic coronaviruses currently described, many of them (i.e. CoV-NL63, CoV-229E, CoV-OC43, and CoV-HKU1) cause mild clinical symptoms leading to flu-like illness, while SARS-CoV, MERS-CoV, and the newly described SARS-CoV-2 cause severe respiratory illness and result in death in patients, primarily in cormorbid patients with underlying health issues. In addition to sever acute respiratory syndrome, COVID-19 is also associated to systemic inflammation, leading to sepsis, acute cardiac injury, heart failure and multi-organ dysfunction in patients at high risk [START_REF] Wang | Review of the 2019 novel coronavirus (SARS-CoV-2) based on current evidence[END_REF].

The SARS-CoV-2 genome shares 82% sequence identity with SARS-CoV and MERS-CoV and > 90% sequence identity for essential enzymes and structural proteins [START_REF] Naqvi | Insights into SARS-CoV-2 genome, structure, evolution, pathogenesis and therapies: Structural genomics approach[END_REF]. SARS-CoV-2 comprises four major structural proteins, including the surface spiculated glycoprotein (S) which presents a cell receptor binding domain and the capsid protein (N) [START_REF] Liu | Evaluation of nucleocapsid and spike protein-based enzyme-linked immunosorbent assays for detecting antibodies against SARS-CoV-2[END_REF]. Protein N is highly antigenic and widely expressed during infection [START_REF] Liu | Evaluation of nucleocapsid and spike protein-based enzyme-linked immunosorbent assays for detecting antibodies against SARS-CoV-2[END_REF]. The antibodies directed against protein N do not show cross-reactivity with antigens of seasonal coronaviruses i.e. NL63, 229E, OC43 and HKU1 [START_REF] Guo | Profiling early humoral response to diagnose novel coronavirus disease (COVID-19)[END_REF]. The Spike protein is the molecular determinant of viral attachment to the primary cell host receptor Angiotensin converting enzyme-2 (ACE-2), permitting fusion and entry into cells [START_REF] Zost | Potently neutralizing and protective human antibodies against SARS-CoV-2[END_REF]. Thus antibodies directed against the S protein are thought to neutralize virus entry into the cell [START_REF] Dogan | Novel SARS-CoV-2 specific antibody and neutralization assays reveal wide range of humoral immune response during COVID-19[END_REF]. Hospitalized patients with severe infection have higher antibody levels and higher neutralization capacity than outpatients [START_REF] Ma | Serum IgA, IgM, and IgG responses in COVID-19[END_REF].

The gold standard of COVID-19 diagnosis is RT-PCR on nasopharyngeal swab, or lower respiratory track samples in severely ill patients. Serological tests help determine whether a person with negative RT-PCR or without RT-PCR diagnosis has been previously or currently infected with SARS-CoV-2. The detection of IgM and / or IgG is optimal from day +15 of the infection. The seroconversion rate is higher for total Ab, than for IgM or IgG (93.1%, 82.7% and 64.7% respectively, [START_REF] Zhao | Antibody Responses to SARS-CoV-2 in Patients With Novel Coronavirus Disease[END_REF]. The median seroconversion is also shorter for total Ig (11 days) than for IgM (12 days) and IgG (14 days) [START_REF] Zhao | Antibody Responses to SARS-CoV-2 in Patients With Novel Coronavirus Disease[END_REF]. The level of IgM and IgA is stable from after infection, whereas the rate IgG continues to increase after D15 [START_REF] Guo | Profiling early humoral response to diagnose novel coronavirus disease (COVID-19)[END_REF], and the increase in IgG levels seems to correlate with the decrease in viral load in the respiratory tract [START_REF] Jin | Diagnostic value and dynamic variance of serum antibody in coronavirus disease 2019[END_REF].

The lateral flow immunoassay (LFIA) PolyDTech test is based on an immunochromatographic method that detects simultaneously anti-SARS-CoV-2 IgM and IgG using a SARS-CoV-2 recombinant antigen; the full principle for IgM antibodies is depicted in Figure 1. Typically, when IgM antibodies are present in the blood/serum sample, they react with gold nanoparticles (GNPs) labeled at their surfaces with a SARS-CoV-2 recombinant antigen. After the buffer is added into the well, the IgM/GNP complex elutes to combine with a mouse anti-human IgM monoclonal antibody set at the surface of the test (T) line. The adduct formed develops a brown-red color at the T line, in which corresponds to a positive test. As a control, in the absence of IgM SARS-CoV-2 antibodies the complex is not formed, resulting in no color at the T line.

The principle of the SARS-CoV-2 IgG detection kit is similar, except that the GNPs for the detection are labeled with a mouse anti-human IgG monoclonal antibody which forms a complex in the presence of SARS-CoV-2 IgG antibodies. The complex migrates upon buffer addition and is set at the surface of the T line labeled with the SARS-CoV-2 recombinant antigen, resulting in a red-brown color at the T line.

As a control (C) experiment, GNPs labeled with a rabbit IgG antibody (for IgM control), and GNPs labeled with a mouse anti-human IgG monoclonal antibody (for IgG control) are also present in the kit. Upon elution with the buffer, these control GNPs migrate and set at the surface of the C line which are coated with a goat anti-rabbit IgG antibody and a goat anti-mouse IgG antibody, revealing a brown-red color at the C line. The absence of a red C line indicates an invalid test.

The method was evaluated on clinical samples in two different hospital centers (namely H and N).

The specificity was evaluated on blood samples collected before the SARS-CoV-2 pandemic arrived in France. The sensitivity was evaluated using blood samples from patients with confirmed COVID-19, diagnosed using molecular methods. The PolyDTech method was compared to another commercial LFIA (BioSynex BSS IgM/IgG) as well as a commercial chemiluminescent Immunoassay (CLIA) method (Maglumi Snibe).

Material and Methods

Negative samples

The negative samples consisted of serum sample collected in two hospital Centers before the SARS-CoV-2 pandemic was reported in France. All serum samples were stored at -20°C before use. Additionally, to complete the evaluation of specificity, eight sera from patient previously tested Influenza positive and four with rheumatoid factor have been tested in Center H.

Samples from COVID-19 patients

For sensitivity evaluation, first, sequential samples of two SARS-CoV-2 positive patients have been tested. Sera were sampled at day (D+)17, D+24, D+38, D+66 and D+112 post-infection for one patient (EL) and at D+10, D+16, D+31 and D+108 for the other one (CH).

In center H, the positive control group consists of 63 laboratory-confirmed COVID-19 patients or hospital staff (method?) that were sampled at least seven days after their first clinical signs. For each participant, information on clinical situation and the disease's kinetics were collected from a survey. The PolyDTech LFIA method permits the qualitative detection of IgM and IgG anti-SARS-CoV2.

According to the manufacturer's instructions, the serum sample (20 µL) was deposited in the well of the device. In a second step, 80 to 100 µL of the buffer were added into the well. The reading of the test is made after 15 to 20 minutes.

The presence of a red-brown color on C line is required to validate the test. The presence of a brownred line on T line revealed the detection of IgM and/or IgM SARS-CoV-2 antibodies in the sample.

Comparison of methods

The comparison was done on 27 samples ranging from D+5 and D+112 post-symptoms. The PolyDTech LFIA Method was compared to another LFIA method, the BIOSYNEX COVID-19 BSS IgG/IgM, and a CLIA IgM+IgG method, Maglumi Snibe.

The antigen used in the Biosynex COVID-19 BSS method corresponds to the RBP domain of S protein.

Analyses were realized according to the manufacturer recommendations. Briefly 10 µL of sera was deposited in the S well and 2 droplets of diluent were deposited in the B well of the test cassette.

The result was read after 10 minutes and before 20 minutes after the deposit. The control line must be present to validate the test. The detection of IgM and/or IgG is associated with the presence a colored line in the IgM and/or IgG zone.

The Maglumi SARS-CoV-2 antibodies IgM and IgG CLIA test used two antigens (S and N protein) coated on magnetic beads. The quantity of antibodies in a sample (i.e. 10 µL) is correlated with the relative light units (RLU) measured by chemiluminescence. A Cut off value was 1.00 UA/mL as recommended by the manufacturer.

Results

Specificity of the assay

The specificity was evaluated on 53 serum specimens collected before the emergence of SARS-CoV-2 in France. All samples were negative for IgG and IgM (Table 1). Additionally, no cross-reaction was observed when tested eight sera from influenza-positive patients, or four sera positives for rheumatoid factor.

Sensitivity of the assay

Overall, 93 serum samples were collected from laboratory-confirmed COVID-19 patients.

For the 84 individual serum samples, the period from onset to serum samples extends from 5-116 days (Table 2). The global sensitivity for all the samples was 53.6% (45/84) for IgM, 85.7% (72/84) for IgG, and 90.5% (76/84) for IgM/IgG. Serum samples were subdivided in subgroups according to sample collection time. The samples taken before D+7 were all negative. For samples taken between D+7 and D+10, and between D+11 to D+13, the IgM/IgG sensitivity was 75.0% and 81.8% respectively. Before between D+7 to D+10, IgM were detected in 7/12 samples, while IgG were detected in 6/12 samples. Between D+11 to D+13, IgM were detected in 6/11 samples, while IgG were detected in 8/11 samples. Form D+14, the sensitivity of IgG detection reached 100.0%, and consequently the overall sensitivity of IgM/IgG reached 100.0% In these samples, IgM were detected in 11/18 samples before D+21, and 21/41 from D+21.

Overall, the Negative Predictive Value (NPV) was calculated to be 86.9% and the positive predictive value (PPV) reached 100.0% for IgM/IgG detection (Table 3). The NPV reached 93.0% for samples from D+7 and 100% for samples from D+14.

Iterative samples were obtained for two patients (n=5 and n=4). One was negative for IgM for all samples from D+17 to D+112 but positive for IgG for all the 5 samples, while the other one was still negative for IgM and IgG from D+10 to D+108 after symptoms (Table 4).

Comparison with other LFIA and CLIA method

The FLIA PolyDtech was compared to another FLIA method, Biosynex BSS G/M and to the Maglumi Snibe CLIA method. Eighteen samples were analyzed using the three methods. All samples from D+13 (n=8) were consistently positive for IgM and/or IgG. For earlier samples, 2/9 were discordant between both FLIA methods (sampled at D+8 and D+10 post symptoms onset), and 1/9 was discordant between PolyDtech FLIA and CLIA method, sampled at D+10.

Overall, the concordance of IgM/IgG detection between the two LFIA methods reached (16/18) 89.0% and the concordance of IgM/IgG detection between the PolyDtech LFIA method and the Maglumi Snibe CLIA method reached (16/17) 94.0% (Table 4).

Elsewhere, two patients with iterative samples were analysed (n=5 and n=4). For one patient, all five samples ranging from D+17 to D+112 were positive for IgM and/or IgG whatever the method used.

For the other patient, usual serological kinetic was only observed using the Biosynex FLIA.

Discussion

The recent outbreak of a novel coronavirus leads to the need for new diagnosis methods.

LFIA has the advantages to be easy to implant, quick to realize and doesn't require device. While RT-PCR is the only diagnostic method currently validated for SAR-CoV-2 infection, serological tests are needed to evaluate the seroprevalence of the infection and to improve the knowledge of COVID-19 spread, thus improving combat of the disease. Serologic tests can also allow confirmation of a SARS-CoV-2 infection if RT-PCR was not able to be realized or was negative despite characteristic symptoms. The viral antigens recommended in the context of serological tests are protein S "spike", especially the RBD domain, or protein N [START_REF] Liu | Evaluation of nucleocapsid and spike protein-based enzyme-linked immunosorbent assays for detecting antibodies against SARS-CoV-2[END_REF].

In this study we evaluated a new LFIA, for IgM and IgG anti-SARS-CoV-2 detection based on anti-spike protein immunoglobulins detection. This test was independently evaluated on two hospital laboratories in North-Eastern France, where the epidemic massively developed during the first three months of 2020.

In SARS-CoV-2 serological tests development, Spike protein and the nucleoprotein have been used as antigens because they are highly antigenic. The specificity of anti-Spike antibodies has been shown to be higher than the specificity of N protein [START_REF] Ayouba | Multiplex detection and dynamics of IgG antibodies to SARS-CoV2 and the highly pathogenic human coronaviruses SARS-CoV and MERS-CoV[END_REF]. The specificity of the method was evaluated on pre-pandemic serum from 2019 November and August. All pre-pandemic samples were negative for anti-SARS-CoV-2 IgM and IgG. In the context of new viral emergence, specificity could be impacted by cross-reactivity with other pathogens, circulating in the same place and time, meaning the seasonal coronavirus (CoV-NL63, CoV-229E, CoV-OC43, and CoV-HKU1) or previously emerging coronaviruses like SARS-CoV or MERS-CoV. Previous studies have demonstrated that cross-reaction with seasonal coronavirus are weak. They are more frequent with SARS-CoV and at a lower rate with MERS-COV in accordance with phylogeny of theses coronaviruses [START_REF] Ayouba | Multiplex detection and dynamics of IgG antibodies to SARS-CoV2 and the highly pathogenic human coronaviruses SARS-CoV and MERS-CoV[END_REF]. In countries where MERS-CoV and SARS-Cov have not experienced intensively circulation, this cross-reactivity may not have an impact on the specificity of serological methods. However, the specificity of other ELISA and CLIA methods have been previously evaluated and reached 99% and 92.2% for IgM, and 99% and 100% for IgG, respectively [START_REF] Infantino | Diagnostic accuracy of an automated chemiluminescent immunoassay for anti-SARS-CoV-2 IgM and IgG antibodies: an Italian experience[END_REF]. Elsewhere specificity of another immunochromatographic assay was evaluated to be 98% [START_REF] Imai | Clinical evaluation of an immunochromatographic IgM/IgG antibody assay and chest computed tomography for the diagnosis of COVID-19[END_REF]. The specificity of the LFIA evaluated here was 100%, in conformity with French recommendations.

According to Liu et al., the sensitivity of protein S-based ELISA for IgM was higher than that of ELISA based on the N protein [START_REF] Liu | Evaluation of nucleocapsid and spike protein-based enzyme-linked immunosorbent assays for detecting antibodies against SARS-CoV-2[END_REF]. The sensitivity of ELISA S-based assay was found to reach 60.5% between D+6 and D+10 after disease onset, 96% between D+16 and D+20, and 100% after D+31 [START_REF] Ong | Comparison of diagnostic accuracies of rapid serological tests and ELISA to molecular diagnostics in patients with suspected coronavirus disease 2019 presenting to the hospital[END_REF]. For the method evaluated here, IgM was the most sensitive before D+11 (58% vs 50% for IgG), while from D+11, IgG was the most sensitive (88%). The global sensitivity of IgM/IgG reached 97% from D+10 and 100% from D+14. The French sanitary authority (Haute Autorité de Santé, 2020) recommends the use of automated ELISA tests, detecting IgM and IgG or total Ig, and displaying a minimum sensitivity of 90 to 95% and a minimum specificity of 98%. As serological diagnosis is recommended from D+14 post symptoms onset, the method is meeting the requirements of HAS recommendations. Otherwise, the anti-S antibodies seem to persist longer in serum and thus more accurately detect the presence of past infection [START_REF] Ripperger | Detection, prevalence, and duration of humoral responses to SARS-CoV-2 under conditions of limited population exposure[END_REF], even if it was not evaluated here.

The kinetics of anti-SARS-CoV-2 seroconversion is still not understood completely. Multiple situations have been described, including early IgG seroconversion before D+7 and IgG seroconversion earlier than IgM [START_REF] Long | Antibody responses to SARS-CoV-2 in patients with COVID-19[END_REF]. These different situations are also observed in the present study. Surprisingly, one patient was negative at D+10-108 for IgM and IgG; crucially this serum was also negative with numerous other dependable methods, including Maglumi Snibe IgG (Eurobio) and Platelia SARS-CoV-2 Total Ab (bioRad) (data not show), highlighting the advances still to be made in the field. Indeed, it has been described that non-hospitalized patients may have lower antibodies rate [START_REF] Grzelak | SARS-CoV-2 serological analysis of COVID-19 hospitalized patients, pauci-symptomatic individuals and blood donors[END_REF] which can adversely affect the probability of detection Overall, this LFIA method detecting anti-SARS-CoV-2 IgM and IgG displays appreciable specificity and sensitivity. It meets the requirements of recommendations for serological diagnosis of COVID-19, when RT-PCR is negative or cannot been realized, provided that patient has developed symptoms for more than 14 days.

  These pre-pandemic serum samples have been used to evaluate the specificity of the PolyDTech LFIA method. The panel of negative serum sample comprised 23 samples from April 2019 for Center H and 30 samples from December 2019 for Centre N.

  In center N, positive control group consists of 21 laboratory-confirmed COVID-19 patients (Viral RNA was extracted using the NucliSens® technology (bioMerieux France) and amplified by RT-PCR protocols targeting RdRp gene developed by Institut Pasteur Paris). For each, the day of the initial appearance of clinical indicators was collected from medical records. Thus, 93 samples from PCR-confirmed positive patients, ranging from D+5 to D+112 post-symptoms onset, were tested. 2.3. Lateral Flow Immunoassay for detection of IgM and IgG antibodies for SARS-CoV-2