François Fages
email: francois.fages@inria.fr

A Constraint-based Mathematical Modeling Library in Prolog with Answer Constraint Semantics

Keywords: constraint logic programming, algebraic modeling languages, answer constraints, MiniZinc, meta-predicates, constraint solving, constraint simplification, attributed variables, ISO-Prolog

Constraint logic programming emerged in the late 80's as a highly declarative class of programming languages based on first-order logic and theories with decidable constraint languages, thereby subsuming Prolog restricted to equality constraints over the Herbrand's term domain. This approach has proven extremely successfull in solving combinatorial problems in the industry which quickly led to the development of a variety of constraint solving libraries in standard programming languages. Later came the design of a purely declarative front-end constraint-based modeling language, MiniZinc, independent of the constraint solvers, in order to compare their performances and create model benchmarks. Beyond that purpose, the use of a high-level modeling language such as MiniZinc to develop complete applications, or to teach constraint programming, is limited by the impossibility to program search strategies, or new constraint solvers, in a modeling language, as well as by the absence of an integrated development environment for both levels of constraint-based modeling and constraint solving. In this paper, we propose to solve those issues by taking Prolog with its constraint solving libraries, as a unified relation-based modeling and programming language. We present a Prolog library for high-level constraint-based mathematical modeling, inspired by MiniZinc, using subscripted variables (arrays) in addition to lists and terms, quantifiers and iterators in addition to recursion, together with a patch of constraint libraries in order to allow array functional notations in constraints. We show that this approach does not come with a significant computation time overhead, and presents several advantages in terms of the possibility of focussing on mathematical modeling, getting answer constraints in addition to ground solutions, programming search or constraint solvers if needed, and debugging models within a unique modeling and programming environment.

Introduction

Jean-Baptiste Fourier is very well-known for his numerous contributions to mathematics and physics, but much less known for being the father of what is called today Constraint Programming. In two lectures at the French Academy of Sciences given in 1823 and 1824, he considered the example of determining the region where a given weight can be placed on a triangular table with constraints on the maximum forces that can be exerted on each leg, in order to promote a general problem solving method based on, first, modeling the problem at hand with inequalities over real numbers, and second, solving them by applying general purpose simplification rules a.k.a. Fourier Motzkin's elimination rules today. Quoting him, "the advantage of this method consists in that it is sufficient, in all cases, to express the conditions of the question, which is easy, and to then combine these expressions, by means of general rules which are always the same; and we thus form the solution which could only have been reached by a series of very complicated reasonings". Furthermore, he mentions that these simplification rules presented for linear inequalities apply as well to non-linear inequalities a.k.a. interval arithmetic: "if the conditions are expressed by non-linear inequalities, the question does not change its nature, and can still be treated by the same principles.". In his second lecture, he considers the optimisation problem and describes a method for moving from vertices to vertices of the feasible region polyhedron, in order to improve the cost up to its optimal value, a.k.a the geometrical interpretation of the Simplex algorithm [START_REF] Lassez | On Fourier's algorithm for linear arithmetic constraints[END_REF].

In the realm of Linear Programming today, several Algebraic Modeling Languages such as OPL, AMPL, Mosel, etc. have been defined as input modeling languages with a syntax near to mathematical notations for optimization problems, using indices, sets, algebraic expressions and data handling variables. Then external solvers interfaced with these modeling languages can be called to actually solve the problem instances defined by the model and a dataset.

In the realm of Constraint Programming, MiniZinc is a similarly high-level constraint-based modeling language to model constraint satisfaction and optimization problems in a solver-independent way [START_REF] Nethercote | MiniZinc: Towards a standard CP modelling language[END_REF][START_REF]The Zinc team: MiniZinc web page[END_REF]. A MiniZinc model is usually transformed in a FlatZinc model in which the high-level constructs have been eliminated and replaced by a flat constraint satisfaction problem that can be solved by a variety of constraint solvers parsing FlatZinc syntax. This is the way for instance SICStus-Prolog is interfaced with FlatZinc as a general purpose constraint solver to solve problems modeled in MiniZinc [START_REF] Carlsson | Sicstus 4.2.3[END_REF]. Because of the importance of global constraints, several of them are defined in MiniZinc, kept in FlatZinc, and handled by the constraint solver, either directly as global constraints if they are implemented as such, or by decomposition into basic constraints supported by the constraint solver.

Because of the importance of search in constraint programming, MiniZinc contains some predefined common search heuristic options that are kept in FlatZinc for their possible interpretation by the constraint solvers, in addition to solver specific annotations in the MiniZinc model that are similarly kept in FlatZinc [START_REF] Nethercote | MiniZinc: Towards a standard CP modelling language[END_REF]. Nevertheless, the absence of programming constructs in a modeling language like MiniZinc limits the search strategies that can be tried to solve a complex problem, which may finally lead to abandon the high-level modeling approach in favor of a lower level constraint programming language.

Furthermore, when it comes to teaching constraint-based modeling and algorithms for decision making, there is a need to teach constraint solving algorithms and program them as well. The recourse to a modeling language like MiniZinc makes it possible to focus on the declarative modeling aspects but cannot be used to show the implementation of constraint solvers, nor implement new constraints. Because of its roots in first-order logic, Prolog should be a natural choice to address these needs as both a programming language in its own right, and a modeling language based on relations a.k.a. constraints. Nevertheless the historical focus on list data structure, and the definition of meta-predicates before the advent of constraint logic programming, make of Prolog a poor modeling language for complex constraint satisfaction problems, moreover exhibiting incompatibilities of use between some standard meta-predicates and constraint solving libraries.

In this paper we propose solutions to those issues by presenting a general purpose mathematical modeling library1 in Prolog, based on bounded quantifiers and arrays for subscripted variables, which is essentially compatible with MiniZinc models and obviously Prolog programming facilities for programming search, constraint solvers, and using external interfaces.

The next section contains a motivating example for modeling the N-queens problem and breaking all the symmetries of that placement problem on a square (i.e. with respect to the 8 elements of the dihedral group of isometries of order 4). The recourse to mathematical notations made possible with our modeling library using subscripted variables and bounded quantifiers rather than lists and recursion, marks a striking gain in declarativity for engineers and students in engineering schools more acquainted to mathematical notations than programming structures. We provide some performance figures which show that this is achieved with no computation time overhead.

In Sec. 3 we compare our modeling library in Prolog to the solver-independent constraint-based modeling language MiniZinc [START_REF] Nethercote | MiniZinc: Towards a standard CP modelling language[END_REF][START_REF]The Zinc team: MiniZinc web page[END_REF]. We show the large compatibility between both modeling languages, and describe the advantages of the Prolog library approach in terms of ability to compute answer constraints, not just ground solutions, a definite advantage of constraint logic programming [START_REF] Jaffar | Constraint logic programming[END_REF] illustrated there with Fourier's example, as well as to program search and debug models in a unique modeling and programming environment.

In Sec. 4 we describe the main predicates of our modeling library for arrays and iterator meta-predicates for answer constraint semantics. The development of such a library raises however a number of subtle issues that cannot be solved with standard Prolog meta-predicates that have been designed before the advent of constraints. This leads us to make a proposal for extending the specification of the behaviour of ISO-Prolog meta-predicates to attributed variables, in the realm of Prolog answer constraint semantics. Finally, we conclude on the importance of keeping the development of Prolog systems in both perspectives of a relation-based declarative programming language and a constraint-based mathematical modeling language.

Motivating Example

List Recursion versus Array Iteration in the N-Queens Problem

The N-queens puzzle consists in placing N queens on an NxN chessboard such that they do not attack each other, i.e. no pair of queens is on the same line or diagonal. This is a classical example used in constraint programming to illustrate the power of the active use of constraints to prune the search tree when they are posted in advance in the paradigm of "constrain and generate", by opposition to pure backtracking "generate and test" programs limited to small size problems. This makes it possible for instance to place by search 100 queens in 0.25s CPU time on a standard laptop, a performance out of reach of search methods without active use of constraints to filter the domains of all variables during search.

The standard definition of this problem by Prolog clauses with constraints over the integers uses a list of N variables representing the queens (say per column as in Fig. 1) with a finite domain of possible values in the interval [1,N] (representing their position in the rows) and posts the no-attack constraints by double recursion on that list of variables, before enumerating the values in the domain of the variables (predicate labeling/2 with ff first-fail smallest-domain variable and-choice heuristics), as follows: The double recursion on the lists of queens and of their successor queens on the right, is standard in Prolog where the list data structure is promoted, Fig. 1. Placement of 8 queens on a 8x8 chessboard with no two queens in the same column, row or diagonal, modeled as a constraint satisfaction problem, with 8 variables Q1, . . . , Q8, representing the queens say in each column, taking integer values in the interval [START_REF] Apt | Constraint Logic Programming using Eclipse[END_REF][START_REF] Jaffar | Constraint logic programming[END_REF] representing the rows where they are placed.

mainly for historical reasons including comparison to Lisp, and where most builtin predicates are defined on lists.

In algebraic modeling languages however, mathematical notations using subscripted variables, with iteration over subscripts instead of recursion, are usually preferred for better declarativity. This is what offers our mathematical modeling library in Prolog. The double list recursion for posting the constraints between each pair of queens in the program above can thus be replaced by iteration on the indices of an array of decision variables as follows:

:-use_module(modeling). Compared to the previous Prolog program using lists and recursion, the advantage in terms of declarativity of the constraint-based model using subscripted variables should be clear, at least for engineers or students more familiar with mathematical notations than programming structures.

Breaking Symmetries in the N-Queens Problem

The gain in declarativity becomes even more striking however, when it comes to breaking all symmetries of the chessboard square with respect of the 8 isometries of the dihedral group of order 4 for a square [A,B,C,D] as in Fig. 1, using lexicographic ordering constraints [START_REF] Walsh | General symmetry breaking constraints[END_REF], namely the reflection symmetry (AB)(CD) around the vertical axis, a variable symmetry in the model of Fig. 1 which can be eliminated here by constraining the first variable to be less than the last; the horizontal axis reflection (AD)(BC), i.e. a value symmetry that can be eliminated here by constraining the first variable to be less or equal to the mid row value; the diagonal reflection (BD), a variable-value symmetry that can be broken by a lexicographic ordering constraint between the variables of the primal model and the dual model defined by D i = j ⇔ Q j = i; the second diagonal reflection (AC), similarly broken using a dual model; the rotation by 90°(ABCD) broken in the same manner; the rotation by 180°(AC)(BD); the rotation by 270°(ADCB); -(the identity).

These symmetry breaking constraints can be added to the model using lexicographic ordering lex_leq/2 constraints on arrays as follows: Implementing arrays in ISO-Prolog is straightforward using term meta-predicates functor/3 and arg/3, and reading/writing cells by unification (backtrackable and non-backtrackable assignments using setarg/3 or nb_setarg/3 are also possible but not considered here).

On the other hand, implementing bounded quantifiers and iteration metapredicates like for_all/2 above raises some issues with the existing ISO-Prolog meta-predicates which have been defined for the provability semantics of Prolog, before the advent of the answer constraint semantics of Constraint Logic Programming [START_REF] Jaffar | Constraint logic programming[END_REF]. This justifies a new implementation of control meta-predicates compatible with constraint solving libraries in Prolog, and probably the definition of a second level of normalization for ISO-Prolog with attributed variables for constraint handling. These aspects are discussed in Section 4.3. We report here some performance figures obtained with our modeling library in SWI-Prolog version 9.0.4, including the clpfd library of constraints over integer variables, on a MacBook Pro 2,3 GHz Quad-Core Intel Core i7 32 GB 3733 MHz.

Performance Figures

Table 1 shows the absence of significant difference in computation time between the execution of our models based on arrays and universal quantifiers for iterations, and the Prolog programs based on lists and list recursions.

On the other hand, the speed-up due to logarithmic access in an array versus linear access in a list (which does not occur in the list based model of the Nqueens problem) begins to show up from size 1000 in pure cell access problems. Type declarations are necessary in MiniZinc to overload operators and constraint predicates over the four constraint domains considered in MiniZinc: the integers, the Booleans, the real numbers (floating points) and the domain of finite sets of integers. This is not implemented in our modeling library in Prolog since constraint predicate symbols are currently not overloaded. More specifically constraint predicates over integers and Booleans values 0/1 of library clpfd are prefixed and distinguished by the # symbol, constraints over real numbers are distinguished by enclosing between curly brackets in the clpr library used, constraints on finite sets are currently not implemented.

Comparison to MiniZinc Modeling Language

It is worth noting however that prescriptive typing and type inference are possible in Prolog, as demonstrated for instance in [START_REF] Fages | Typing constraint logic programs[END_REF] using a powerful system of subtyping constraints [START_REF] Coquery | TCLP: overloading, subtyping and parametric polymorphism made practical for constraint logic programming[END_REF]. This could be used to parse MiniZinc models and type check mixed MiniZinc-Prolog programs.

Programming Search

The search procedure is the second most important component of constraint programming techniques to solve hard combinatorial problems. For that reason, MiniZinc offers the possibility to provide a limited set of search annotations, which can be interpreted by the solvers, for specifying and-choice heuristics (e.g. first-fail principle), or-choice heuristics (e.g. best-first criterions) and also some search strategies such as dichotomic search.

Of course, such search options cannot cover all needs for controlling search. In [START_REF] Martinez | Search by constraint propagation[END_REF], it is however shown how some more elaborated search strategies such as Limited Discrepancy Search [START_REF] Harvey | Limited discrepancy search[END_REF] can be specified by constraints in the model. In [START_REF] Rendl | MiniSearch: A solver-independent meta-search language for MiniZinc[END_REF], a solver-independent language is proposed to control search in MiniZinc with a limited set of instructions to execute each time a solution is found. This is illustrated by the implementation of Large Neighborhood Search (LNS), an important strategy for optimization problems.

In [START_REF] Schrijvers | Search combinators[END_REF], search combinators are proposed to specify general search instructions to execute at each choice point in the model. This requires interaction with the solver at each node of the search tree which is hardly supported by constraint solvers implemented in a different programming language.

The capability to program search at that level of granularity can however make a decisive difference in hard combinatorial problems. For instance, in disjunctive scheduling, it can be worth spending time to determine the best mutual exclusion constraint to select for the next choice point, in order to prune the search tree by constraint propagation, instead of just duplicating the real search space with a bad selection. Sophisticated strategies, including lookingahead techniques for determining that most informative disjunctive constraint by trying them all may thus become necessary to implement.

This requires the recourse to a programming language which is possible with a programming library-based approach to modeling, like here in Prolog, a natural choice for relation-based modeling, and not possible with a pure modeling language approach like MiniZinc or other algebraic modeling languages, beyond a limited set of search annotations.

Answer Constraint Semantics in Fourier's Example

Thanks to the computed answer constraint semantics of the class of constraint logic programming languages [START_REF] Jaffar | Constraint logic programming[END_REF], and of Prolog with its constraint-solving libraries, the execution of MiniZinc models directly in our modeling library in Prolog makes it possible to obtain much more informative answer constraints, compared to the ground solutions obtained with MiniZinc solvers.

This can be illustrated by the example given by Fourier in his lecture on constraint-based modeling and solving for systems of linear inequalities [START_REF] Lassez | On Fourier's algorithm for linear arithmetic constraints[END_REF]. He took the problem of placing a given weight P at unknown coordinates (X, Y ) on an isocele triangular table with a maximum force F exerted on each of the 3 legs (see Fig. 2).

By assuming that the triangle table has leg 1 at coordinates (0, 0), leg 2 at coordinates (20, 0) and leg 3 at coordinates (0, 20), simple moment equilibrium equations give the following model: Fig. 2. Example given by Fourier to illustrate a general purpose constraint-based modeling method, applied to the placement of a weight on a triangular table with limit constraints on the forces exerted on its legs 1, 2, 3, together with a general purpose method for solving systems of linear inequalities over the real numbers. The inner triangle represents the validity domain of the placement of a weight of 2 units. This triangle is symbolically represented by the computed answer constraint for the query fourier(2, X, Y, 1).

For a weight of 3 units, a single point is returned as unique solution. Interestingly, for the placement of a weight of 2 units maximizing its coordinate along the X axis, the value 10 is found for X together with one constraint on Y and A which defines the vertical segment where the weight can be placed. Similarly, the computed answer constraints returned for the placement of a weight of 2 units without optimization criterion define symbolically the validity domain depicted by the inner triangle in Fig. 2. It is worth noting that the computed answer constraint set semantics are a quite unique advantage of the constraint logic programming paradigm, compared to other declarative modeling languages that are usually interfaced to solvers restricted to compute ground solutions and not constraints.

Model Debuging and Visualization in a Unique IDE

In our modeling library approach in Prolog, a MiniZinc model can be executed in two ways:

either by parsing the MiniZinc syntax and interpreting it directly by the predicates of our library, or by parsing the FlatZinc model generated by the MiniZinc compiler which is the standard way of evaluating backend MiniZinc constraint solvers.

A FlatZinc parser has been developed for SICStus Prolog which shows excellent performance on MiniZinc challenge competitions [START_REF]The Zinc team: MiniZinc web page[END_REF].

On the other hand, the advantage of parsing MiniZinc syntax in a Prolog system has already been demonstrated in a system like Eclipse [START_REF] Apt | Constraint Logic Programming using Eclipse[END_REF]. Making it available in a library for standard Prolog systems is however a significant contribution to factorize system development efforts, and use the single Integrated Development Environment (IDE) of Prolog to debug MiniZinc models.

For instance, the general purpose search tree visualization and interaction library CLPGUI, developed for constraint logic programming [START_REF] Fages | CLPGUI: a generic graphical user interface for constraint logic programming[END_REF], can be used here to visualize, and restore states in, the search tree developed by our modeling library in Prolog. This is illustrated in Fig. 3 in the model above for enumerating all 92 solutions of the 8 queens problem.

Fig. 3. Visualization of the search tree for enumerating all the 92 solutions of the 8queens problem in Prolog with constraints, using the interactive graphical visualization system CLPGUI of [START_REF] Fages | CLPGUI: a generic graphical user interface for constraint logic programming[END_REF]. This illustrates the good job of constraints for pruning the search tree in this example. Indeed, after around 5 choice points with a limited number of possible values, the branches of the search tree become linear and lead to one solution.

From the point of view of a modeler, or of a teacher, this is a very useful tool that is difficult to implement without a unique declarative programming/modeling IDE. For the same reasons, model debuging is greatly facilitated by the absence of a first transformation of the model to a FlatZinc model, followed by a second transformation to an external constraint solver which can create runtime errors difficult to recover.

Modeling Library

4.1 Subscripted Variables with Library arrays.pl Subscripted variables, i.e. arrays, are not normalized in ISO-Prolog but can be easily represented by terms with functor array. Multidimensional arrays are then represented by arrays of arrays. On top of that, library modeling.pl defines predicates for defining arrays of Booleans, integers or real number with domain constraints. The main predicates defined in library arrays.pl are: array(?Array, ?Dimensions) true if Array is an array of dimension Dimensions. Either creates an array of given dimensions, or returns the list of dimensions of a given array, or checks the relation.

cell(+Array, +Indices, ?Cell). true if Cell is the cell or subarray cell of a given array at given indices (integer or list of integers). Throws an error if the indices are out of range of the array. The cell is read and written by unification. The functional notation syntax Array[Indices] is also authorized in constraints to denote the cell of an array without having to introduce an existentially quantified variable for that, e.g. with the quantifier meta-predicates described below.

array_list(?Array, ?List). true if List is the flat list of all elements in Array. Either creates a one-dimensional array, or the list of array cells, or checks the consistency of both representation.

array_lists(?Array, ?List). true if List is the list of lists of elements in the array following its dimensions.

tensor(+A, +Op, +B, +Rel, ?C) equivalent to (A Op B) Rel C, where A, B, C are arrays or lists of same dimensions, Op is a binary operation executed element-wise, and Rel a binary predicate.

Furthermore, for the sake of generality outside the scope of our modeling library, backtrackable and non-backtrackable imperative array cell assignments are also defined by extra-logical predicates set_cell/3 and nb_set_cell/3 using ISO-Prolog predicate setarg/3.

Bounded Quantification Meta-Predicates in quantifiers.pl

Iteration with subscripted variables is more natural to implement using quantifiers on the indices rather than recursion. To this end, we introduce the following general-purpose meta-predicates to quantify variables either universally or existentially, and distinguish them from context variables in goals: for_all(+Args, +Goal) calls Goal for all the arguments specified in Args, i.e. calls the conjunction Goal(Arg1), ..., Goal(ArgN). The list of arguments Args is a list of universally quantified variables given with a finite domain (either finite interval of integers or list of terms), plus possibly a deterministic condition goal, using expressions of the form Args in Domain where Condition. exists(+Vars, +Goal) to make a list of variables renamed-apart and local to the goal. let(+Bindings, +Goal) to additionally apply a binding constraint =, =.., in, #=, #<, #>, #=<, #>= between an existentially quantified variable and a term, possibly using Array[Indices] functional notation for array cells. list_of(+Vars, ?List) to build the list of terms specified in Args.

The introduction of those meta-predicates in our modeling library is justified by the fact that ISO-Prolog meta-predicates bagof/3, setof/3, findall/3 have been introduced before the advent of constraint logic programming, and of the implementation of constraint solvers based on attributed variables, at a time focussing on the provability semantics of Prolog. Indeed, in his seminal paper [START_REF] Warren | Higher-order extensions to Prolog: Are they needed?[END_REF], David Warren introduced an extra logical meta-predicate to collect information across branches of the search tree. Paraphrasing him, the extension takes the form of a new built-in predicate: setof(X,P,S) to be read as: "The set of instances of X such that P is provable is S". This is achieved with a special mechanism introduced in the WAM [START_REF] Warren | An abstract Prolog instruction set[END_REF] for copying and memorizing terms in a list across backtracking. The non deterministic goal P can thus be used as a generator of the instances to be checked for satisfiability.

This mechanism has been used to introduce various meta-predicates for higher-order programming in Prolog [START_REF] Sagonas | Effcient execution of HiLog in WAM-based prolog implementations[END_REF]. Those meta-predicates thus refer to the provability semantics of Prolog, not the answer constraint semantics. Similarly, meta-predicates forall/3, foreach/3 found in several Prolog dialects have no normalized behaviour with respect to attributed variables and constraints.

In the success set semantics of a universal quantifier for goals, it is not clear whether the successfull bindings should be kept across backtracking. On the other hand, in the answer constraint semantics, it is clear that the computed answer constraints should be added conjunctively to the store of constraints. This is precisely what is achieved by our for_all/2 meta-predicate that generates instances of the goal by iteration, unlike for example the forall/2 predicate of SWI-Prolog that generates instances by backtracking:

?-L=[X, Y], forall(member(V, L), V=a). L = [X, Y]. ?-L=[X, Y], for_all(V in L, V=a). L = [a, a], X = Y, Y = a.

Interface to Constraint Solvers with Library clp.pl

We also found it useful to create a front-end interface to existing libraries for solving constraints over the real numbers, integers and booleans, for several reasons. First, this front-end library makes it possible to use evaluable expressions and Array[Indices] functional notation in constraints and in the definitions of variable domains. Second, global constraints on lists of variables can now undifferently accept arrays instead of lists.

Beyond that, some questions arise concerning the answer constraint semantics and the copying or not by Prolog meta-predicates of the constraints attached to the variables. For constraint libraries using attributed variables, this amounts to the question of whether those meta-predicates should copy the attributes or not, in particular the constraint propagator attributes.

In the answer constraint semantics, the answer is clearly yes, and this is what is done for instance in SWI-Prolog and Scryer-Prolog. There is however currently a side effect concerning the duplication of the constraint propagators, e.g.

?-L=[A, B], L ins 1..2, A #=< B, bagof(W, member(W, L), L2). L = L2, L2 = [A, B], A in 1..2, B#>=A, B#>=A, B#>=A, B in 1..2.
This is patched in our library clp.pl by making a set union instead of concatenating the constraint propagators attached to two variables when they are unified. In this respect, it is worth noting that the set union of constraint propagators should not be implemented by maintaining the constraint propagation attributes sorted, since that could severely affect the scheduling and performance of constraint propagators.

Furthermore, a more efficient and more general solution would be to consider the simplification of constraints on two variables that get unified. This is not usually done in constraint programming and global constraint solvers, e.g. Symbolic simplification of constraints is however instrumental in SMT solvers or in Constraint Handling Rules (CHR) [START_REF] Frühwirth | Constraint Handling Rules[END_REF]. This has been shown responsible for drastic performance improvement in some contexts, for instance in [START_REF] Coquery | TCLP: overloading, subtyping and parametric polymorphism made practical for constraint logic programming[END_REF][START_REF] Coquery | Subtyping constraints in quasi-lattices[END_REF] for solving subtyping constraints using CHR, with better performance than dedicated algorithms in CAML, thanks to the constraint simplifications performed by CHR immediately upon unification of two variables.

Proposal for a Second-level of Normalization of ISO-Prolog

Because of the importance of constraint-based methods in many applications of Prolog, and because of the numerous implementations of constraint solvers in Prolog libraries using attributed variables, we propose to specify the behaviour of ISO-Prolog predicates with respect to attributed variables in a new level 2 norm for ISO-Prolog.

More specifically, we mainly propose to open discussion to:

1. normalize attributed variables; 2. specify the copying or not of attributed variables in ISO-Prolog predicates findall/3, bagof/3, setof/3; 3. define meta-predicates compatible with constraint answer semantics; 4. normalize predicates for arrays.

Conclusion

Last year was the 50th year of Prolog, and this remarkable longevity for a programming language was clearly recognized as a mark of the fundamental importance of the relational programming paradigm based on first-order Horn clause logic and constraints in decidable theories. It also showed the need to transmit knowledge to new generations of Prolog system developpers, teachers and users, and unite the Prolog community with new momentum.

The development of our constraint-based mathematical modeling library in Prolog aims at going in that direction from the triple points of view of the users, by providing a higher-level modeling library with full programming features; of the teachers, by providing a unique environment for both constraint-based modeling and programming; and of the system developpers, by factorizing library development efforts, and providing use cases of standard predicates and libraries that require some corrections and probably a normalization effort.

On-going work concerns the addition of a MiniZinc parser to our modeling library in order to directly execute and debug MiniZinc models in a single modeling/programming environment, compare performances using the large database of MiniZinc models, compute answer constraints not just ground solutions for getting more declarative answers, and stop being blocked from programming search for hard decision making problems.

:

  -use_module(library(clpfd)). queens(N, Queens):length(Queens, N), Queens ins 1..N, safe(Queens), labeling([ff], Queens). safe([]). safe([QI | Tail]) :noattack(Tail, QI, 1), safe(Tail). noattack([], _, _). noattack([QJ | Tail], QI, D) :-QI #\= QJ, QI #\= QJ + D, QI #\= QJ -D, D1 #= D + 1, noattack(Tail, QI, D1).

  queens(N, Queens):-int_array(Queens, [N], 1..N), for_all([I in 1..N-1, D in 1..N-I], (Queens[I] #\= Queens[I+D], Queens[I] #\= Queens[I+D]+D, Queens[I] #\= Queens[I+D]-D)), satisfy(Queens).

  sym_elim(N, Queens) :-Queens[1] #< Queens[N], % vertical reflection symmetry Queens[1] #=< (N+1)//2, % horizontal reflection symmetry int_array(Dual, [N], 1..N), for_all([I, J] in 1..N, Queens[I] #= J #<==> Dual[J] #= I), lex_leq(Queens, Dual), % first diagonal reflection int_array(SecondDiagonal, [N], 1..N), for_all(I in 1..N, SecondDiagonal[I] #= N + 1 -Dual[N+1-I]), lex_leq(Queens, SecondDiagonal), int_array(R90, [N], 1..N), for_all(I in 1..N, R90[I] #= Dual[N+1-I]), lex_leq(Queens, R90), % rotation symmetry by 90°i nt_array(R180, [N], 1..N), for_all(I in 1..N, R180[I] #= N + 1 -Queens[N+1-I]), lex_leq(Queens, R180), int_array(R270, [N], 1..N), for_all(I in 1..N, R270[I] #= N + 1 -Dual[I]), lex_leq(Queens, R270).

3. 1 Types

 1 Compared to MiniZinc, the writing of the N-queens model is very similar:int: n = 8; array [1..n] of var 1..n: queens; constraint forall (i in 1..n-1, d in 1..n-i) (queens[i] != queens[i+d] /\ queens[i] != queens[i+d] + d /\ queens[i] != queens[i+d] -d); solve satisfy;

  10.0*_A-10.0*_B, X=10.0*_B, _=2.0-_A-_B, _B=<1.0, _A=<1.0}. ?-fourier(2, X, Y, 1), maximize(X). X = 10.0, {Y=10.0-10.0*_A, _=1.0-_A, _A=<1.0, _A>=0.0}. ?-fourier(2, X, Y, 1), maximize(X+Y). X = Y, Y = 10.0.

?

  -L=[X, Y], L ins 1..3, all_distinct(L), X=Y. L = [Y, Y], X = Y, Y in 1..3, all_distinct([Y, Y]).

Table 1 .

 1 CPU time in seconds and reported number of logical inferences in SWI-Prolog for solving various variants of the N-queens problem, compared between the recursive program on lists and the mathematical model using subscripted variables and quantifiers. The last benchmark simply shows the logarithmic versus linear access times in large datasets represented by arrays versus lists.

	Benchmark	CPU time in seconds Nb logical inferences
	100-queens first solution
	math model	0.373	4,552,518
	list program	0.298	4,472,406
	all-distinct math model	3.611	63,565,626
	all-distinct list program	3.470	63,502,308
	8-queens all 92 solutions
	math model	0.041	757,835
	list program	0.042	757,303
	all-distinct math model	0.297	4,025,643
	all-distinct list program	0.289	4,012,622
	8-queens all 12 non-symmetrical sol.
	math model	0.045	667,850
	list program	0.048	664,708
	all-distinct math model	0.125	1,471,542
	all-distinct list program	0.128	1,509,989
	10 2 accesses in an array of size 10 2	0.000	301
	in a list of size 10 2	0.000	544
	10 3 accesses in an array of size 10 3	0.001	3,001
	in a list of size 10 3	0.003	4,144
	10 4 accesses in an array of size 10 4	0.005	30,001
	in a list of size 10 4	0.161	40,001
	10 5 accesses in an array of size 10 5	0.027	300,001
	in a list of size 10 5	13.615	400,001
	10 6 accesses in an array of size 10 6	0.148	3,000,001
	in a list of size 10 6	1567.177	4,000,144

The Prolog libraries presented here form a pack named modeling, currently available for SWI-Prolog at https://lifeware.inria.fr/wiki/Main/Software#modeling.

Acknowledgments. I am grateful to my students at Ecole Polytechnique for their interest in my courses on Constraint Programming including practical work that evolved over the years from Prolog to MiniZinc and now back to Prolog; to Mathieu Hemery and Sylvain Soliman for their participation in the teaching and fruitful discussions; to Guy-Alain Narboni for his vision of the importance of the Prolog heritage and the organization of the 50th year of Prolog in Paris; to Markus Triska, Ulrich Neumerkel and Christian Jendreiko for their organization of the Scryer Meetup in Dusseldorf; and to the reviewers for their comments.


Stéphan Rouzière
email: stephan.rouziere@u-psud.fr

Pascale Launois

Ana M Benito

Wolfgang K Maser

Erwan Paineau

Unravelling the hydration mechanism in a multi-layered graphene oxide paper by in-situ X-ray scattering

Hydration of a multi-layered graphene oxide (GO) paper in water vapor atmosphere was studied by in-situ X-ray scattering over a wide range of relative humidity (RH). The intercalation of water molecules physically adsorbed between GO layers induces the expansion of the interlayer distance. Two regimes of adsorption are clearly evidenced according to the variation of the interlayer distance as a function of the relative humidity. The expansion of the interlayer distance is limited to ~1 Å for RH ranging from 6 to 53%, whereas it is ~3.5 Å for RH 75 to 100%. The hydration mechanism corresponds to the progressive adsorption of water molecules onto hydrophilic sites at the surface of GO layers in the low-RH regime, then an additional water monolayer is formed in the high-RH regime.

Introduction

Unexpected physical properties are observed when molecules are confined in nanoscale environments. One of the most blatant examples is probably the tremendously enhanced water transport in carbon nanostructures. Ultra-fast water transport has been revealed in hydrophobic carbon nanostructures, either inside carbon nanotubes or through nanopores created on a graphene sheet [START_REF] Majumder | Enhanced flow in carbon nanotubes[END_REF][START_REF] Holt | Fast mass transport through sub-2-nanometer carbon nanotubes[END_REF][START_REF] Surwade | Water desalination using nanoporous single-layer graphene[END_REF][START_REF] Guo | Nanofluidic transport through isolated carbon nanotube channels: Advances, controversies, and challenges[END_REF]. Contrariwise, graphene oxide (GO) is a hydrophilic nanoporous material in which water can easily adsorb between layers [START_REF] Lerf | Hydration behavior and dynamics of water molecules in graphite oxide[END_REF][START_REF] Barroso-Bujans | Sorption and desorption behavior of water and organic solvents from graphite oxide[END_REF]. Water molecules can permeate through interconnected nanochannels formed in-between GO layers, with a nearly frictionless motion over the non-oxidized surfaces [START_REF] Chen | Observation and analysis of water transport through graphene oxide interlamination[END_REF][START_REF] Devanathan | Molecular dynamics simulations reveal that water diffusion between graphene oxide layers is slow[END_REF][START_REF] Deng | Mechanism of water transport in graphene oxide laminates[END_REF]. Based on the permeation properties of water, free-standing GO membranes [START_REF] Dikin | Preparation and characterization of graphene oxide paper[END_REF] are now considered as possible new materials for outperforming current technologies for water purification [START_REF] Nair | Unimpeded permeation of water through helium leak-tight graphene based membranes[END_REF][START_REF] Mi | Graphene oxide membranes for ionic and molecular sieving[END_REF][START_REF] You | Graphene and graphene oxide for desalination[END_REF], molecular sieving [START_REF] Rk | Precise and ultrafast molecular sieving through graphene oxide membranes[END_REF][START_REF] Abraham | Tunable sieving of ions using graphene oxide membranes[END_REF] or humidity sensing [START_REF] Bi | Ultrahigh humidity sensitivity of graphene oxide[END_REF].

Although several models are proposed to describe the complex and defective GO structure [START_REF] Dreyer | The chemistry of graphene oxide[END_REF], the most accepted structural model (Lerf-Klinowski model [START_REF] Lerf | Structure of graphite oxide revisited[END_REF]) defines graphene oxide as a non-stoichiometric material. A GO sheet is a graphene sheet covered randomly with oxidized functional groups (OFGs), namely, epoxide and hydroxyl groups attached to the basal plane and carbonyl or carboxyl groups at the sheet edges, together with residual nonfunctionalized graphitic areas. For stacked GO sheets, the presence of OFGs leads to an expansion of the spacing between layers, from 3.35 Å in pristine graphite to an interlayer distance value d ~ 5.9 to 6.7 Å in GO assemblies [START_REF] Boehm | Untersuchungen am Graphitoxyd. IV. Vergleich der Darstellungsverfahren[END_REF]. It also imparts swelling properties to the material. By analogy with swelling clays [START_REF] Ferrage | Hydration properties and interlayer organization of water and ions in synthetic Na-smectite with tetrahedral layer charge. Part 1. Results from X-ray diffraction profile modeling[END_REF], the interlayer distance can be further increased by crystalline swelling, up to ~ 12.5 Å, according to the degree of hydration of the material [START_REF] Lerf | Hydration behavior and dynamics of water molecules in graphite oxide[END_REF][START_REF] Barroso-Bujans | Sorption and desorption behavior of water and organic solvents from graphite oxide[END_REF][START_REF] Buchsteiner | Water dynamics in graphite oxide investigated with neutron scattering[END_REF][START_REF] Rezania | Hydration of bilayered graphene oxide[END_REF][START_REF] Talyzin | The structure of graphene oxide membranes in liquid water, ethanol and water-ethanol mixtures[END_REF][START_REF] Rouzière | Intercalated water in multi-layered graphene oxide paper: an X-ray scattering study[END_REF][START_REF] Liu | Graphene oxide papers with high water adsorption capacity for air dehumidification[END_REF]. However, it was pointed out by several authors [START_REF] Lerf | Hydration behavior and dynamics of water molecules in graphite oxide[END_REF][START_REF] Talyzin | The structure of graphene oxide membranes in liquid water, ethanol and water-ethanol mixtures[END_REF][START_REF] Krishnamoorthy | The chemical and structural analysis of graphene oxide with different degrees of oxidation[END_REF] that the degree of hydration in GO is strongly dependent on the sample preparation method. The method of synthesis (Brodie, Staudenmaier or Hummers' methods [START_REF] Dreyer | The chemistry of graphene oxide[END_REF]), the sample shaping (powder, paper or membrane) or the hydration protocol (saturated water vapor atmosphere vs. immersion in liquid water) lead to variations in the resulting interlayer distance, as illustrated in table S1. In the present work, we focus on the crystalline swelling of a GO paper obtained by Hummers' method in water vapor atmosphere.

The expansion of the interlayer distance d during hydration is directly related to the adsorption of water molecules intercalated in-between the regularly stacked GO layers.

Several mechanisms of adsorption have been proposed in previous studies, taking into account the maximum value d max after complete hydration or how the distance d changes as a function of relative humidity. However, even for the same d max measured, description of the hydration mechanism can notably differ. For instance, Cerveny et al. [START_REF] Cerveny | Dynamics of water intercalated in graphite oxide[END_REF] attributed the expansion of the interlayer distance to the uptake of a water monolayer whereas Liu et al. [START_REF] Liu | Graphene oxide papers with high water adsorption capacity for air dehumidification[END_REF] recently suggested the presence of two layers of water. Furthermore, an interstratification model in which empty interlayer spaces are filled by a variable number of water monolayers was proposed [START_REF] Lerf | Hydration behavior and dynamics of water molecules in graphite oxide[END_REF][START_REF] Buchsteiner | Water dynamics in graphite oxide investigated with neutron scattering[END_REF][START_REF] Talyzin | Temperature dependent structural breathing of hydrated graphite oxide in H2O[END_REF], while hydration process in water vapor phase was explained by a continuous incorporation of water molecules which changes gradually the thickness of GO layers [START_REF] Rezania | Hydration of bilayered graphene oxide[END_REF]. Obviously, a complete view of the hydration mechanism involved is still missing.

The aim of the present work is to detail the hydration mechanism in a GO paper in water gas atmosphere. A particular attention has been paid to the experimental method used for the study. In-situ measurements were performed with a single sample in a dedicated X-ray humidity-controlled chamber. Variation of the interlayer distance was measured as a function of RH ranging from 0 to 100%, and the reversibility of the swelling process was checked at each RH step. Finally, the small number of layers (four layers) stacked in GO assemblies allowed controlling a homogeneous hydration. Two regimes of adsorption have been clearly evidenced, a low-RH regime corresponding to localized adsorption of water molecules onto the surface of GO layers, and a high-RH regime reflecting the additional formation of a water monolayer which is more likely to contribute to the water transport properties in graphene oxide.

Experimental

Sample preparation

A dispersion of graphene oxide was prepared from graphite flakes (flake size > 100 mesh, Sigma-Aldrich, ref. 332461) using a modified Hummers' method by oxidation with NaNO 3 , H 2 SO 4 and KMnO 4 as reported earlier [START_REF] Vallés | Flexible conductive graphene paper obtained by direct and gentle annealing of graphene oxide paper[END_REF]. The obtained graphite oxide was dispersed in distilled water, ultrasonicated and centrifuged, resulting in an aqueous dispersion of exfoliated graphene oxide. For preparation of freestanding flexible GO papers, 500 mL of the aqueous GO dispersion, with a concentration of 0.15 mg ml -1 , was filtered through an Isopore filter (ϕ 47 mm, pore size 1.2 µm) by a vacuum directed vacuum flow, during which GO sheets form assemblies of a few layers with periodic stacking. The obtained paper with thickness of 20 µm is made of GO assemblies. The GO paper was rinsed with water and dried overnight before being peeled off from the filter. The physico-chemical properties of the GO paper were fully characterized in ref. [START_REF] Núñez | Graphene oxide-carbon nanotube hybrid assemblies: cooperatively strengthened OH⋯O=C hydrogen bonds and the removal of chemisorbed water[END_REF], with a C:O ratio of 4.8 calculated from XPS analysis. In this work, we used a single piece of 20µm-thick GO paper with a size of 4x2 mm 2 approximatively.

X-ray controlled humidity chamber

In-situ X-ray scattering (XRS) measurements were carried out with a monochromatized X-ray beam delivered by a molybdenum rotating anode generator (λ = 0.711 Å, Rigaku HU3R, 40kV-40mA) coupled with a multilayer W/Si mirrors optics (Xenocs FOX 2D Mo 25-INF); X-ray beam size is 0.8 × 1 mm 2 at the sample position [START_REF] Taché | MOMAC: a SAXS/WAXS laboratory instrument dedicated to nanomaterials[END_REF]. The sample was fixed on both ends with ceramic paste on a brass sample holder and dried by heating at 80°C under vacuum for 12 h prior to XRS measurements. Note that a temperature treatment of 80ºC does not affect the chemical composition of our GO paper, as clearly evidenced in previous decomposition studies [START_REF] Núñez | Graphene oxide-carbon nanotube hybrid assemblies: cooperatively strengthened OH⋯O=C hydrogen bonds and the removal of chemisorbed water[END_REF]. Furthermore, the findings of the decomposition studies on our own material [START_REF] Núñez | Graphene oxide-carbon nanotube hybrid assemblies: cooperatively strengthened OH⋯O=C hydrogen bonds and the removal of chemisorbed water[END_REF] are in full agreement with the further literature work [START_REF] Zhu | Pseudonegative thermal expansion and the state of water in graphene oxide layered assemblies[END_REF], which clearly underline that changes in the chemical composition are not taking place in the temperature range from RT to 80 ºC. Then, the sample was placed in a controlled humidity chamber with entrance and exit windows made of 500µm-thick Mylar film for X-rays paths (Fig. 1). The relative humidity (RH) was generated by an aqueous saturated salt solution placed in a reservoir connected to the chamber which was beforehand evacuated with a primary vacuum pump (Fig. 1). To explore the whole range of RH, we used the following salts in saturated conditions [START_REF] Line | Handbook of Chemistry and Physics[END_REF]: LiBr (RH = 6%), LiCl (RH = 11%), CH 3 COOK (RH = 22%), K 2 CO 3 (RH = 43%), Mg(NO 3 ) 2 (RH = 53%), NaCl (RH = 75%), KCl (RH = 84%). The highest relative humidity (RH = 100%) was achieved by using ultrapure water (conductivity σ = 5.5 × 10 -3 mS.m -1 ). The chamber was systematically pumped (P ~ 1 mbar) between each salt solution change which allowed performing all XRS measurements on the same sample and controlling the reversibility of the GO swelling. XRS images were recorded by a large-area X-ray sensitive plate detector with 150 µm pixel size (MAR345, marXperts GmbH), placed behind the exit window of the chamber (Fig. 1b).

Measurements were performed in transmission geometry with GO paper oriented parallel to the incident X-ray beam. Extraction of the scattered intensity 𝐼 as a function of the scattering modulus 𝑄 (𝑄 = 4πsin𝜃/λ where 2𝜃 is the scattering angle) were performed with homedeveloped software.

Results and discussion

X-ray scattering of a multi-layered GO paper

Thanks to the specific geometry used for XRS measurements, the respective contributions from in-plane and out-of-plane structural features of GO paper are clearly observed (Fig. 2 inset). The preferential orientation of stacked GO layers leads to angular anisotropy of the 001 reflection located at the scattering wave-vector 𝑄 !!" . The wave-vector value is directly related

to the interlayer distance d by 𝑑 = !! ! !!"
. Peak position decreases (interlayer distance increases) during hydration, as it is illustrated in Fig. 2. A second anisotropic peak is also observed in the same direction at a larger value Q ~ 1.78 Å -1 , which is due to remaining impurities from the synthesis, arising from the periodic stacking of non-oxidized graphitic sheets. Accordingly, this latter peak is insensitive to hydration ( [START_REF] Rouzière | Intercalated water in multi-layered graphene oxide paper: an X-ray scattering study[END_REF] and Fig. 2). In the orthogonal direction, 10 peak reflects the in-plane hexagonal arrangement of carbon atoms in graphene sheets. Its position corresponds to the expected lattice parameter 𝑑 !" = 2.46 Å and the C-C bond length of 1.42 Å (𝑑 !"" / 3) of a graphene sheet, and does not change during hydration either. Finally, position and width of 001 reflection under vacuum (black dotted curve in Fig. 2a) are 𝑄 !!" ~0.98 Å -1 and FWHM = 0.25 Å -1 respectively. Determination of the average layer number N of GO assemblies was performed by fitting the profile of 001 peak according to the formalism developed in ref. [START_REF] Rouzière | Intercalated water in multi-layered graphene oxide paper: an X-ray scattering study[END_REF] using the expression of the scattered intensity 𝐼 !!" 𝑄 ∝ 𝐹 𝑄 ! 𝐺(𝑄), F(Q) being the structure factor of a GO layer, with the interference function G(Q) :

𝐺 𝑄 = 1 𝑁 𝑠𝑖𝑛 𝑁𝑄𝑑 2 𝑠𝑖𝑛 𝑄𝑑 2 !
The results correspond to a stacking of 4 GO layers with an interlayer distance of 6.4 Å, the latter value being in good agreement with previous studies [START_REF] Boehm | Untersuchungen am Graphitoxyd. IV. Vergleich der Darstellungsverfahren[END_REF][START_REF] Zhu | Pseudonegative thermal expansion and the state of water in graphene oxide layered assemblies[END_REF]. and high (84%) RH. The variation in position of 001 reflection is due to the intercalation of water in GO stacking while the diffracting coherent domain remains unchanged with the same number of stacked layers [START_REF] Rouzière | Intercalated water in multi-layered graphene oxide paper: an X-ray scattering study[END_REF]. No additional disorientation between GO layers is observed either. It is noteworthy that 001 reflection returns systematically to initial values of position and width of the dried sample, after re-pumping the chamber before changing the salt solution. Corresponding XRS image with acquisition time of 15 minutes is similar to the one of the dried sample, inferring desorption is achieved in a much shorter time. It indicates, firstly, that vacuum process removes effectively all the water molecules physically adsorbed between GO layers, and secondly, that hydration does not cause any damage to GO structure.

Variation of the interlayer distance as a function of time was measured at different RH ranging from 6 to 100% (Fig. 3a). At each RH, the interlayer distance increases strongly within the first two hours after hydration starts, then more gently to approach a saturation value after about 10 hours. Remarkably, the expansion occurs immediately, from the lowest value of RH 6%, which confirms the high water adsorption capacity of graphene oxide. The saturation value becomes larger when RH is higher, with a maximum value of about 10 Å at RH 100%. Two regimes in the variation of interlayer distance d according to RH are clearly evidenced. For RH from 6 to 53%, the saturation value is restricted within the value range 6.9-7.4 Å whereas it is larger than 9.4 Å for RH above 75%. The time dependence of the variation of the interlayer distance during hydration at each RH value can be fitted according to a simple empirical relation:

𝑑 = 𝑑 1 -𝑒 ! ! ! !.! + 𝑑 !"#$$% (1)
τ being a time constant, 𝑑 !"#$$% the value of the interlayer distance in vacuum (𝑑 !"#$$% = 6.4

Å), and 𝑑 the additional expansion of the interlayer distance at 𝑡 = ∞. The result of fits indicates that the time constant τ is roughly equal to 2 h (±1h) whatever the RH value, during which the interlayer distance d strongly increases. A more accurate determination of τ is out of reach, since the time resolution of the measure is limited by the image acquisition time during which the position of reflection 001 is averaged. A time constant of a few hours was also measured on GO papers by Liu et al. [START_REF] Liu | Graphene oxide papers with high water adsorption capacity for air dehumidification[END_REF]. Moreover, it must be underlined that the saturation value of d measured in the low-RH regime is rapidly overcome for RH higher than 75%, within the first 15 minutes for RH of 100% for instance, which corresponds to enhanced water adsorption in the high-RH regime.

Two-step hydration mechanism

Variation of the expansion Δd as a function of RH is plotted in Fig. 3b. Its value is constant, about 1 Å, for RH varying from 6 to ~55%, while it is around 3.5 Å for RH above 75%. Our results highlight two different regimes of water adsorption which were not clearly identified in earlier studies (see Fig. 4). Previous results show a large dispersion of d values and a rather linear dependence of d with RH up to 80%. The maximum value around 10 Å at RH 100% is consistent with the ones measured on other GO papers by Liu et al. [START_REF] Liu | Graphene oxide papers with high water adsorption capacity for air dehumidification[END_REF] whereas larger values measured by Lerf et al. [START_REF] Lerf | Hydration behavior and dynamics of water molecules in graphite oxide[END_REF] were obtained with samples immersed in liquid water. Various experimental conditions, sample preparation and shape or hydration method (see table S1)

between the different studies could explain the dispersion in the values while RH is varying.

In order to minimize such effects, a special methodology was applied in the present study, allowing us to evidence the existence of two different hydration states depending on RH. Insitu measurements were carried out on the same sample. Vacuum pumping before changing the salt solution ensured a proper measure of the water adsorption at each RH, starting systematically from the same initial state of dried GO paper, and also confirming the reversibility of the swelling process. Let us underline that variation of the interlayer distance is strictly related to the intercalation of water molecules in-between GO layers and not to water molecules adsorbed on the edges of the layers or out of the interlayer spaces. The intercalated water molecules are physically adsorbed and are easily removed by vacuum pumping. Finally, the hydration process was homogeneous. Indeed, one observes a single 001 peak with no anomalous widening during the whole hydration process (Fig. 2). It is probably due to the limited number of stacked GO layers (four layers) in GO assemblies. A simple model of hydration can be deduced from the two regimes of adsorption clearly observed. Considering that the van der Waals diameter of a water molecule is about 3 Å, the expansion of the interlayer distance measured in the low-RH regime which is limited to 1 Å, does not leave enough space for the formation of a layer of water, contrarily to the conclusion drawn by Liu et al. [START_REF] Liu | Graphene oxide papers with high water adsorption capacity for air dehumidification[END_REF]. Instead, water molecules can be adsorbed on specific sites on GO layer's surface such as holes or hydrophilic OFGs. In particular, water molecules can interact with OFGs through double hydrogen bonding [START_REF] Bi | Ultrahigh humidity sensitivity of graphene oxide[END_REF][START_REF] Tarasevitch | Interaction of water, methanol and benzene molecules with hydrophilic centres at a partially oxidised model graphite surface[END_REF][START_REF] Medhekar | Hydrogen bond networks in graphene oxide composite paper: structure and mechanical properties[END_REF], exhibiting localized motions as probed by neutron spectroscopy [START_REF] Lerf | Hydration behavior and dynamics of water molecules in graphite oxide[END_REF][START_REF] Buchsteiner | Water dynamics in graphite oxide investigated with neutron scattering[END_REF]. This type of intercalated water is thus referred as bound water (sketch in Fig. 3b). This first-step hydration mechanism spreads out over a wide RH range [6-53%], resulting in the limited expansion of the interlayer distance.

A second-step hydration regime occurs at high RH, where intercalated molecules are likely to form local clusters containing an increasing number of water-water bonds. Clusters eventually merge to form a water monolayer in between the layers formed of GO and physically bonded water molecules discussed above (sketch in Fig. 3b). The presence of the water monolayer induces the larger expansion of the interlayer distance to 3.5 Å. Due to the change in the hydrogen bonds network [START_REF] Medhekar | Hydrogen bond networks in graphene oxide composite paper: structure and mechanical properties[END_REF], motion of the water molecules within the monolayer is less constrained than the one of bound water, that is water molecules become more mobile and acquire a liquid-like behavior. It is worth mentioning here that high permeation rate of water through membranes designed for filtration and separation techniques [START_REF] Nair | Unimpeded permeation of water through helium leak-tight graphene based membranes[END_REF] is believed to be achieved thanks to the mobile water monolayer, which can only exist in the high-RH regime [75-100%] where the interlayer distance is large enough to accommodate a monolayer of water. The evidence for two different hydration regimes depending on RH is also relevant in the context of high humidity sensitivity of GO films [START_REF] Bi | Ultrahigh humidity sensitivity of graphene oxide[END_REF], where the specific roles of bound water and of mobile water layer are discussed.

In summary, hydration of a multilayered graphene oxide paper in vapor gas atmosphere was studied by in-situ X-ray scattering measurements. Intercalation of water molecules physically adsorbed between the layers induces the expansion of the interlayer distance according to two regimes of adsorption clearly observed for the first time, and corresponding to a two-step hydration mechanism. For RH below 55%, water molecules are bound to functional oxygen groups of GO layers whereas an additional water monolayer is formed at RH above 75%. Our results should stimulate further theoretical and numerical simulations and allow one to get a deeper understanding of some of the outstanding properties of GO for water filtration, molecular sieving or humidity sensing.

Fig. 1 :

 1 Fig.1: (a) Operating principle of the X-ray controlled humidity chamber: schematic front

Fig. 2 :

 2 Fig.2: XRS intensity along 001 direction (intensity along the horizontal red line on XRS

Fig. 3 :

 3 Fig.3: (a) Variation of the interlayer distance as a function of time at different RH. Eye-

Fig. 4 :

 4 Fig.4: Experimental values (red circles and dashed line as a guide for the eyes) of the

Acknowledgements

EP, SR and PL benefited of a grant from the French Research National Agency under the program Investissements d'Avenir (SEDECO project, ANR-10-LABX-0035: Labex NanoSaclay). AMB and WKM acknowledge funding from the Spanish MINEICO (project ENE2016-79282-C5-1-R), the Gobierno de Aragón (Grupo Reconocido TO3_17R) and associated EU Regional Development Funds.

SUPPLEMENTARY INFORMATION

Table S1. Literature data on interlayer distance of GO at maximum swelling (d max ) obtained by X-ray diffraction at 25°C under saturated relative humidity (RH = 100%) if not specified, or immersed in liquid water. Protocol of GO synthesis, sample shaping and hydration method are also reported.

Reference d max (Å) GO preparation Sample shaping

Hydration method

Lerf et al. [