R Abgrall

R Crisovan

ROBUST MODEL REDUCTION OF HYPERBOLIC PROBLEMS BY L 1 -NORM MINIMIZATION AND DICTIONARY APPROXIMATION

Keywords: Reduced order models, nonlinear hyperbolic problems, L 1 minimisation AMS subject classifications. 65M08, 65K99, 76L99

We are interested in the model reduction techniques for hyperbolic problems, in particular in fluids. This paper, which is a continuation of [2], proposes a dictionary approach coupled with a L 1 minimisation approach. We develop the method and analyse it in simplified one dimensional cases. We show in this case that error bound with the full model can be obtained provided a suitable minimisation approach is chosen. The capability of the algorithm is then shown on non linear scalar problems, one dimensional unsteady fluid problems and two dimensional steady compressible problems. The paper ends with a discussion on the cost of the method.

1. Introduction. Model reduction is becoming an essential tool to enable applications requiring either real-time predictions or the evaluation of a large number of partial differential equations (PDE) based computational models. The first category encompasses optimal control [START_REF] Ito | A Reduced-Order Method for Simulation and Control of Fluid Flows[END_REF][START_REF] Ly | Modeling and control of physical processes using proper orthogonal decomposition[END_REF] and model predictive control [START_REF] Hovland | Explicit model predictive control for large-scale systems via model reduction[END_REF][START_REF] Amsallem | Model Predictive Control under Coupled Fluid-Structure Constraints Using a Database of Reduced-Order Models on a Tablet[END_REF]. Routine analysis and parametrized studies [START_REF] Amsallem | Toward real-time computational-fluid-dynamics-based aeroelastic computations using a database of reduced-order information[END_REF], design optimization [START_REF] Legresley | Airfoil design optimization using reduced order models based on proper orthogonal decomposition[END_REF][START_REF] Amsallem | Design Optimization Using Hyper-Reduced-Order Models[END_REF] and the quantification of uncertainty [START_REF] Bui-Thanh | Parametric reduced-order models for probabilistic analysis of unsteady aerodynamic applications[END_REF] are applications pertaining to the second category, to name just a few. In all of these applications, the large dimensionality associated with the discretized partial equations prevents their solution in real-time. Model reduction reduces that cost by restricting the solution to a subspace of the solution space. This subspace is usually described by a small number of reduced basis vectors. In turn, a projection step reduces the dimensionality of the system of discrete equations considered, enabling their fast solution.

While the model reduction of elliptic and parabolic PDEs has been the subject of numerous studies [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics[END_REF][START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF] and its theory is well understood [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF][START_REF] Grepl | A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations[END_REF][START_REF] Rozza | Reduced Basis Approximation and a Posteriori Error Estimation for Affinely Parametrized Elliptic Coercive Partial Differential Equations[END_REF], reducing hyperbolic equations has proved to be much more challenging [START_REF] Barone | Stable Galerkin reduced order models for linearized compressible flow[END_REF]. More specifically, moving waves and discontinuities such as shocks require a large number of basis vectors to accurately approximate these features [START_REF] Dahmen | Double greedy algorithms: Reduced basis methods for transport dominated problems[END_REF]. This characterizes these problems as ones with large Kolmogorov n-widths [START_REF] Binev | Convergence Rates for Greedy Algorithms in Reduced Basis Methods[END_REF].

To circumvent this issue, approaches based on local bases [START_REF] Amsallem | Interpolation method for adapting reduced-order models and application to aeroelasticity[END_REF][START_REF] Dihlmann | Model reduction of parametrized evolution problems using the reduced basis method with adaptive time-partitioning[END_REF][START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF][START_REF] Maday | Locally adaptive greedy approximations for anisotropic parameter reduced basis spaces[END_REF] reduce the Kolmogorov n-width by considering local subspaces. The locality can be characterized in parameters [START_REF] Amsallem | Interpolation method for adapting reduced-order models and application to aeroelasticity[END_REF], time [START_REF] Dihlmann | Model reduction of parametrized evolution problems using the reduced basis method with adaptive time-partitioning[END_REF] or state-space [START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF]. In the present work, an approach based on dictionaries is considered [START_REF] Kaulmann | Online Greedy Reduced Basis Construction Using Dictionaries[END_REF][START_REF] S L Brunton | Compressive sensing and low-rank libraries for classification of bifurcation regimes in nonlinear dynamical systems[END_REF]. More specifically, solutions corresponding to various time and parameter instances are collected and stored in such a dictionary. Each solution will then be considered as a reduced basis vector. In turn, localization in time and space can be easily enforced by only considering basis vectors corresponding to restricted sub-domains of the time and parameters spaces. In addition to the reduction in number of basis vectors, this paper will demonstrate that a key advantage of a dictionary approach is a better approximation of states having sharp gradients and discontinuities. In particular, it will be demonstrated that avoiding basis truncation such as the one occurring in Proper Orthogonal Decomposition (POD) [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. Part I: coherent structures[END_REF] or Non-Negative Matrix Factorization [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF] avoid Gibbs phenomenon.

In addition to the choice of reduced basis, a key ingredient in projection-based model reduction is the definition of the reduced system of equations. For symmetric systems such as those arising in elliptic and parabolic PDEs, Galerkin projection is the method of choice. For non-symmetric systems, however, it has been shown that minimizing the L 2 -norm of the residual is preferable for stability considerations [START_REF] Carlberg | Efficient non-linear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations[END_REF][START_REF] Carlberg | The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows[END_REF]. In the present paper, model reduction based on the minimization of the L 1 -norm of the residual is introduced and its advantage is demonstrated in conjunction with a dictionary approach for reducing problems with sharp gradient and shocks. More specifically, the present work demonstrates that combining a dictionary and L 1 minimization promotes sparsity in the choice of basis functions participating in the reduced-order solution and results in more accurate and physical reduced-order solutions.

In the present paper, we consider as a future work the feasibility and the efficiency of the method and also the computational and the memory costs needed for the storage of the dictionary. This paper is organized as follows: we first discuss the problem of interest. In the following section, an approximation of the solution of nonlinear problems by reduced order models is presented. In the fourth section, we explain the role of L 1 minimization in this problem, then we present in detail the algorithm we have developed and provide an error estimate. Section 6 is a numerical section where we present efficient algorithms for the computation of the L 1 -norm minimizer, both in the cases of linear and nonlinear residuals. Then we provides several numerical examples that illustrate the behavior of our methods, on linear and nonlinear problems, both in 1D and 2D. The last section provides a discussion on the computational cost. A conclusion follows and we sketch some perspectives.

Problem of interest.

In this work, high-dimensional models (HDM) arising from the space discretization of hyperbolic PDEs are considered. PDEs of the following type are considered

     ∂W ∂t + L(W, µ) = f (t, µ), x ∈ Ω, t ∈ [0, T ], B(W, µ) = g(t, µ), x ∈ ∂Ω, t ∈ [0, T ], W (x, t = 0, µ) = W 0 (x, µ), x ∈ Ω (2.1) W ∈ R p is a scalar (p = 1) or vector (p > 1) field, Ω ⊂ R d is the domain of the equation 1 ≤ d ≤ 3
and ∂Ω the boundary of the domain. L is a differential operator such as the divergence of a flux and B a boundary operator, f and g are volume and surface forces, respectively and µ ∈ P ⊂ R m is a vector of m parameters defining the system of interest.

The HDMs result from a finite differences approximation or finite volume formulation of the PDE (2.1) under the following form

dw dt + f (w(t), t, µ) = g(t, µ), t ∈ [0, T ] w(t = 0) = w 0 , (2.2)
where w ∈ R N is the HDM state of large dimension N , t denotes time and f (•, •), g(•) are nonlinear functions of their arguments.

In the remainder of this paper, the time and parameter variables are grouped together, unless explicitly stated, as a variable τ = [t; µ]. Hence, the HDM state is parametrized as

w(τ ) = w(t, µ). (2.3)
In practice, the ODE (2.2) is discretized in time using a time discretization

t 0 = 0 < t 1 < • • • < t Nt = T .
Explicit and implicit time-discretization techniques are used in the present paper, resulting in a sequence of nonlinear systems of equations of large dimension N

r n (w) = 0, n = 1, • • • , N t , (2.4)
where

r n = [r n 1 , • • • , r n N ] T .
We give several examples later in the text. Note that the residual r n will depend on several time instances of the solution for unsteady problems, for example w n and w n-1 in the simplest case. Steady problems can also be written in the form r(w) = 0.

The goal of model reduction is to approximate the high-dimensional system (2.4) using a much smaller number of variables while retaining accuracy of the solution. For that purpose, projection-based model reduction techniques approximate the state w(τ ) in a subspace of R N using a reduced-order basis (ROB

) V = [v 1 , • • • , v k ] ∈ R N ×k .
The state is then approximated as

w(τ ) ≈ Vq(τ ) = k i=1 v i q i (τ ) (2.5)
where q(τ ) = [q 1 (τ ), • • • , q k (τ )] T ∈ R k denotes the vector of k reduced coordinates. Substituting the subspace approximation (2.5) into (2.4) usually results in a non-zero residual of dimension N r n (Vq) ≈ 0.

(2.6)

Two common approaches result in the definition of a reduced system of equations:

• Galerkin projection enforces the orthogonality of the residual to the ROB V as

V T r n (Vq) = 0, n = 1, • • • , N t . (2.7)
This defines a set of k nonlinear equations in terms of k unknowns which can be solved by Newton-Raphson's method. • Residual minimization approaches [START_REF] Legresley | Airfoil design optimization using reduced order models based on proper orthogonal decomposition[END_REF][START_REF] Bui-Thanh | Parametric reduced-order models for probabilistic analysis of unsteady aerodynamic applications[END_REF][START_REF] Carlberg | Efficient non-linear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations[END_REF][START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF][START_REF] Carlberg | The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows[END_REF] minimize the residual in the L 2 -norm sense

min q∈R k r n (Vq) 2 2 = N i=1 (r n i (Vq)) 2 , n = 1, • • • , N t . (2.8)
In practice, this nonlinear least-squares problem can be solved using Gauss-Newton or Levenberg-Marquardt iterations [START_REF] Nocedal | Numerical optimization[END_REF]. In Section 6, alternative residual minimization approaches based on L 1 -norm minimization which are more appropriate for the reduction of hyperbolic problems will be proposed.

3. Dictionary approach. Projection-based model reduction techniques [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. Part I: coherent structures[END_REF][START_REF] Rozza | Reduced Basis Approximation and a Posteriori Error Estimation for Affinely Parametrized Elliptic Coercive Partial Differential Equations[END_REF][START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF] based on snapshots pre-compute solutions of the HDM for specific values of the vector τ = [t; µ]. These snapshots are gathered in a snapshot matrix

S = [w(τ 1 ), • • • , w(τ Ns )]. (3.1)
Three approaches for compressing the snapshot matrix are described as follows:

• Proper Orthogonal Decomposition [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. Part I: coherent structures[END_REF] computes an optimal reduced-order basis of dimension k that minimizes the projection error of the snapshots onto the basis. • Balanced POD [START_REF] Willcox | Balanced model reduction via the proper orthogonal decomposition[END_REF], applicable to linear systems only, also takes into account snapshots of the dual system to construct the reduced basis for the primal and dual systems. • Non-negative matrix factorization [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] was recently applied to construct a non-negative reduced-order basis based on snapshots with positive entries in the context of contact problems [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF]. The reduced basis minimizes the positive reconstruction of the snapshots. All three approaches perform a compression of the information contained in the snapshot matrix S. More specifically, the N s vectors contained in S are compressed, leading to a reduced-order basis of dimension k ≤ N s .

In the present paper, an approach based on a dictionary of solutions is preferred as it does not incur any loss of information by compression. As such, the vectors {v i } k i=1 in the reduced basis are solutions of the HDM:

v i = w(τ i ), i = 1 • • • , k. (3.2)
So, the number of samplings corresponds to the number of elements in the reduced basis. In this case, the error can be controlled using the error estimate given in Section 5.2. The solution of the HDM will then be approximated as

w(τ ) ≈ k i=1 w(τ i )q i (τ ). (3.3)
In the present case, since the HDM is of very large dimension, over-complete dictionaries, as used in compressed sensing [START_REF] Candes | Robust Signal Recovery from Incomplete Observations[END_REF][START_REF] Donoho | Compressed sensing[END_REF] and for which k ≥ N will not be considered.

4. L1 -norm residual minimization. In the present paper, model reduction based on L 1norm residual minimization is introduced to reduce the dimensionality of hyperbolic equations as an alternative to Galerkin projection and L 2 -norm minimization. Motivations for the use of the L 1 -norm are provided in this section. Model reduction based on L 1 -norm minimization is introduced in Section 6 together with practical numerical procedure for their computation in Section 6.2.

Minimizing the L 1 -norm of the residual is known to lead to regressions that are much more robust to outliers [START_REF] Boyd | Convex optimization[END_REF]. In the context of hyperbolic systems, the work of Guermond et al. on Hamilton Jacobi equations and transport problems [START_REF] Guermond | A fast algorithm for solving first-order PDEs by L1minimization[END_REF][START_REF] Guermond | L 1 -Approximation of Stationary Hamilton-Jacobi Equations[END_REF] has shown, at least experimentally, that the numerical solution can retain an excellent non-oscillatory behavior by minimizing the L 1norm of the PDE residual. In [START_REF] Guermond | L 1 -Approximation of Stationary Hamilton-Jacobi Equations[END_REF], the schemes are designed by minimizing quantities that mimic the total variation of a functional, as in here. In the following sections, the idea is exploited in the case of model reduction, for which the set of grid points differs to the set of the full grid points, as an approximation subspace is introduced. For completeness, the motivation for L 1 -norm minimization is justified as follows for the problem

∂W ∂t + div F (W ) = 0 (4.1) defined on Ω ⊂ R d × R + .
The solution W belongs here to R p , so that F = (F 1 , . . . , F p ) T . The weak form of the equation is: for any ϕ ∈ C 1 0 (

• Ω)

p with compact support in the interior of

• Ω of Ω: 1 Ω ϕ(x, t) ∂W ∂t + div F (W ) dtdx = 0. (4.2)
Integrating by parts yields, we get

Ω ∂ϕ ∂t W dtdx + Ω ∇ϕ • F (W )dtdx = 0. (4.3)
Restricting to the set of test functions

T = ϕ ∈ C 1 0 ( • Ω) p , ||ϕ|| ∞ ≤ 1 , W is a solution if: sup ϕ∈T Ω ∂ϕ ∂t W dtdx + Ω ∇ϕ • F (W )dtdx = 0. (4.4)
Remember that for any function g ∈ L 1 (R d ), the total variation is defined as

T V (g) = sup ϕ∈C 1 0 (R d )∩L ∞ (R d ),||ϕ||∞≤1 R d ∇ϕ(x) • g(x)dx , (4.5)
and if, in addition,

g ∈ C 1 (R d ), then T V (g) = R d ||∇g||dx = ||∇g|| L 1 (R d )
. This shows that, defining the space-time flux F = (W, F ), W is a weak solution if and only if the (space-time) total variation of F vanishes, that is

T V F(W )) = 0. (4.6)
In other words, one can look for W as a function of

L 1 ∩ L ∞ such that W minimizes T V F(V )) over V ∈ L 1 ∩ L ∞ , i.e. W = argmin{T V F(V )), V ∈ L 1 ∩ L ∞ }. (4.7)
This does not guaranty uniqueness (and thus there is some abuse of language in this setting), since the entropy conditions are not encoded into this formulation. However, (4.7) indicates that a natural setting is to minimize the L 1 norm of the space-time divergence of the space-time flux F.

How does it translates in the discrete setting? For simplicity, we only mention the case of explicit schemes. We discuss later the solution procedure for the case of implicit schemes. The following useful classical result is mentioned. Consider {x j } j∈Z a strictly increasing sequence in R and

x j+1/2 = xj +xj+1 2 . Assuming that R = ∪ j∈Z [x j-1/2 , x j+1/2
[ and considering g defined by, for any j ∈ Z,

g(x) = g j if x ∈ [x j-1/2 , x j+1/2 [, (4.8)
then

T V (g) = j∈Z |g j+1 -g j |. (4.9)
Now, instead of having the exact solution, consider an approximation procedure that enables, from w n ≈ W ( . , t n ), to compute w n+1 ≈ W ( . , t n+1 ). For instance, assume that we have a finite volume method, d = 1 and for any grid point j ∈ {1, • • • , N }, we define the mesh x j+1/2 = j∆x, t n = n∆t and the control volumes

c j = (x j-1/2 , x j+1/2 ). Considering ϕ = 1 [x j-1/2 ,x j+1/2 ]×[t n ,t n+1 ] in (4.
2), we obtain the following:

0 = x j+1/2 x j-1/2 t n+1 t n ∂W ∂t + div F (W ) dtdx = x j+1/2 x j-1/2 W (x, t n+1 )dx - x j+1/2 x j-1/2 W (x, t n )dx + ∆t t n+1 t n 1 ∆t F (W (x j+1/2 , t))dt -∆t t n+1 t n 1 ∆t F (W (x j-1/2 , t))dt.
Using the approximations,

w n j ≈ 1 ∆x x j+1/2 x j-1/2 W (x, t n )dx and f j+1/2 (w n ) ≈ 1 ∆t t n+1 t n F (W (x j+1/2 , t))dt, we obtain: ∆x(w n+1 j -w n j ) + ∆t f j+1/2 (w n ) -f j-1/2 (w n ) = 0. (4.10)
In this case, the residual can be written as:

r(w n , w n+1 ) j = w n+1 j -w n j + ∆t ∆x f j+1/2 (w n ) -f j-1/2 (w n ) . (4.11)
Here f j+1/2 is any numerical flux at the cell interface x j+1/2 ; see [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] for the classical examples. One way to evaluate w n+1 is to minimize the total variation, i.e.

T V (r) = j∈I r(w n , w n+1 ) j , (4.12)
leading to

w n+1 = argmin v piecewise constant j∈I ∆x(v j -w n j ) + ∆t(f j+1/2 (w n ) -f j-1/2 (w n )) . (4.13)
In this case I differs from the full set of grid points as a subspace of R N is introduced. Clearly, if I is equal to the set of grid points, the solution is given by

w n+1 j = w n j - ∆t ∆x f j+1/2 (w n ) -f j-1/2 (w n ) . (4.14)
Substituting (2.5) in (4.11), for all time steps, find the coefficients q n i that minimizes the following residual:

r n j (V q) = k i=1 q n+1 i w n+1 j (µ i ) - k i=1 q n i w n j (µ j ) + ∆t ∆x f i+1/2 (w n ) -f i-1/2 (w n ) . (4.15)
For the case of implicit schemes, the solution is obtained similarly as for the explicit case, obtaining the following residual that has to be minimized:

r n j (V q) = k i=1 q n+1 i w n+1 j (µ i ) - k i=1 q n i w n j (µ j ) + ∆t ∆x f i+1/2 (w n+1 ) -f i-1/2 (w n+1 ) . (4.16)
remark 1. The L 1 norm is convex but not strictly convex, and hence the minimisation problem (4.12) may not have a unique solution. For this reason, in practice, we perturb the functional (4.12) to make it strictly convex. Let us denote it by J, and thus we will look for solutions that minimize J(r(w n , w n+1 )).

(4.17)

Examples are: 1. for ν > 0,

J(r(w n , w n+1 )) = j∈I r(w n , w n+1 ) j + ν j∈I (w n+1 j ) 2 , (4.18a) 2. More generaly, if U is a convex entropy, J(r(w n , w n+1 )) = j∈I r(w n , w n+1 ) j + ν j∈I U (w n+1 j ) (4.18b)
The functional (4.18b) can be used for systems.

Error estimation.

In this section, we provide an error estimate (in the scalar case) between the solution obtained by projection on the dictionary and the solution of the original scheme. This error estimate is another way to justify the method. These estimates are provided in a simple setting: we consider monotone scheme. In this section, we first precise the setting, then give a natural condition on the dictionary for obtaining this error estimate and then state and prove it.

5.1. Scheme setting. Consider the scalar conservation law equations with the initial condition:

∂u ∂t + ∂f (u) ∂x = 0, x ∈ R, t > 0 u(x, 0) = u 0 (x), x ∈ R.
(5.1)

After discretizing, we assume that the scheme writes, for u := (u j ) j∈Z ,

u n+1 = S(u n , λ) (5.2)
with λ = ∆t/∆x and the initial condition

u 0 j = given. (5.3)
We assume that the operator S is monotone for λ ∈ [0, b[, b > 0, i.e. if for any sequence u and v bounded for the L 1 or L ∞ norms with j ∈ Z, u j ≤ v j , then S(u j , λ) ≤ S(v j , λ). Let L 1 and L ∞ norms generically denoted by || . ||.

An example is given by the scheme

S(u j ) = u j -λ f (u j+1 , u j ) -f (u j , u j-1 ) (5.4)
where we assume that the numerical flux f (a, b) is monotone, i.e. increasing with respect to the first variable and decreasing with respect to the second one. S is monotone under a CFL like condition. Another example is given by the implicit scheme, where v j = S(u j ) is defined as the solution of

v j = u j -λ f (v j+1 , v j ) -f (v j , v j-1 ) (5.5)
which is unconditionally monotone. Thanks to Crandall-Tartar lemma (for example, see [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF]), we know that for any u and v,

||S(u, λ) -S(v, λ)|| ≤ ||u -v||
in the L 1 norm. The same is true in the L ∞ and L 2 norms.

Error estimate.

We collect and store in a dictionary the solutions {u n (µ i )} i of the problem (5.1) which correspond to various time and parameter instances and where the initial conditions are defined for the parameters {µ i } i=1,...,m ∈ P. Since the minimisation procedure admits a unique solution, this enables to define a projection operator p n for any time t n , by solving the minimization problem: knowing that u n µ ∈ span µi∈P

({u n (µ i )}), find u n+1 µ ∈ span µi∈P ({u n+1 (µ i )})
such that for ∀µ ∈ P:

u n+1 µ = argmin µi∈P J v µ -S(u n (µ), λ) : v µ ∈ span µi∈P ({u n+1 (µ i )}) = p n S(u n µ , λ) .
with J strictly convex.

We have immediately the following estimate:

J p n (S(u n µ , λ)) -S(u n (µ), λ) = min vµ∈ span µ i ∈P ({u n+1 (µi)}) J v µ -S(u n (µ), λ) ≤ min µi∈P J u n+1 (µ i ) -S(u n (µ), λ) = min µi∈P J S(u n (µ i ), λ)) -S(u n (µ), λ) ≤ min µi∈P J u n (µ i ) -u n (µ) (5.6)
provided λ enables to fulfill the monotonicity property for all the elements of the dictionary. The last inequality is possible due to the monotonicity of the scheme for the functionals defined in Remark 1 because a monotone scheme is L 1 stable, L 2 stable and T V stable.

Next, we consider the case associated with the minimization problem (6.2). If α i ≥ 0 and |P| i=1 α i = 1, we obtain a sharper error estimate of type:

J p n (S(u n µ , λ)) -S(u n (µ), λ) = J |P| i=1 α i u n+1 (µ i ) -u n+1 (µ) ≤ |P| i=1 |α i | J u n+1 (µ i ) -u n+1 (µ) ≤ |P| i=1 |α i | max µi∈P J u n+1 (µ i ) -u n+1 (µ) = max µi∈P J u n+1 (µ i ) -u n+1 (µ) ≤ max µi∈P J u 0 (µ i ) -u 0 (µ) .
(5.7)

Again, the last inequality is possible due to the monotonicity of the scheme for the functionals defined in Remark 1 because a monotone scheme is L 1 stable, L 2 stable.

We have shown the following result: Proposition 5.1. Consider J defined as in Remark 1. If S( . , λ) is monotone for λ ∈ [0, b[ then 1. At time t n+1 , the minimization is done on span {u n (µ i ), µ i ∈ P}, the reduced solution u n+1 µ = p n (S(u n µ , λ)), ∀µ ∈ P satisfies:

J u n+1 µ -S(u n (µ), λ) ≤ min µi∈P J u n (µ i ) -u n (µ) .
2. At time t n+1 , the minimization is done on the convex hull of the dictionary, the reduced solution u n+1 µ = p n (S(u n µ , λ)), ∀µ ∈ P satisfies:

J u n+1 µ -S(u n (µ), λ) ≤ max µi∈P J u 0 (µ i ) -u 0 (µ) .
For J defined by

J(r) = j∈I |r j | + ν j∈I |r j | 2 ,
we get the following estimate:

J (p n (S(u n µ , λ)) -S(u n (µ), λ)) = min vµ∈ span µ i ∈P ({u n+1 (µi)}) J (v µ -S(u n (µ), λ)) ≤ min µi∈P J (u n+1 (µ i ) -S(u n (µ), λ)) = min µi∈P J (S(u n (µ i ), λ)) -S(u n (µ), λ)) ≤ min µi∈P J (u n (µ i ) -u n (µ)) ≤ min µi∈P ||u n (µ i ) -u n (µ)|| 1 + ν||u n (µ i ) -u n (µ)|| 2 2 ≤ min µi∈P 1 + ν||u n (µ i ) -u n (µ)|| ∞ ||u n (µ i ) -u n (µ)|| 1 ≤ min µi∈P 1 + ν||u 0 (µ i ) -u 0 (µ)|| ∞ ||u 0 (µ i ) -u 0 (µ)|| 1
where the penultimate inequality holds due to the log-convexity of L p norm. remark 2.

1. In practice, the equation (5.7) can be seen as an algorithm of selecting the parameters {µ i } i=1,...,m by fixing a desired precision such that max µi∈P J u 0 (µ i ) -u 0 (µ) ≤ .

In this way, we can control from the beginning the parameters which will constitute the training set P from which we will compute the reduced basis. This step it is not expensive because it is needed to compute just the initial condition for each element {µ i } i=1,...,m ∈ P. 2. Using the functional (4.18b), the error estimators presented in this chapter are holding also for the system case due to the work of Bressan on L 1 stability estimates for hyperbolic systems of conservation laws [START_REF] Bressan | Hyperbolic systems of conservation laws; The one-dimensional Cauchy problem[END_REF] at least for initial conditions with small enough bounded variation. In that case, the L 1 norm of the solution can be estimated from the L 1 norm of the initial condition, and the entropy is diminishing.

6. Model Reduction by L 1 -norm minimization and dictionaries.

6.1. Procedure. In this section, model reduction based on minimizing the residual in the L 1 -norm is combined with the dictionary approach presented in Section 3.

As an alternative to Galerkin projection and residual minimization in the least-squares sense, a reduced system of equation is here obtained by minimizing, at each time step n = 1, • • • , N t , the L 1 -norm of the residual vector as

min q∈R k r n (Vq) 1 = n i=1 |r n i (Vq)| , n = 1, • • • , N t (6.1)
or as

min q∈R k r n (Vq) 1 subject to 1 T q = 1, q 0, n = 1, • • • , N t , (6.2)
where we project on the convex envelop of the dictionary. There are at least three difficulties associated with minimizing the L 1 -norm. A first one is associated to L 1 its non-differentiability at zero. To circumvent this issue, the Huber function [START_REF] Peter | Robust Statistics[END_REF], defined as follows can be introduced:

φ M (x) = x 2 if |x| ≤ M M (2|x| -M ) otherwise, (6.3)
Then, the sequence of reduced systems of equations based on the Huber function is

min q∈R k n i=1 φ M (r n i (Vq)) , n = 1, • • • , N t . (6.4)
The Huber function φ M behaves as a parabola close to x = 0 and as the L 1 -norm for large values of x. It is continuously differentiable on R (φ M ∈ C 1 (R)). It is also used in regressions as a loss function due to its non-sensitivity to outliers. In the present work, it will be used as a continuously differentiable alternative to the L 1 -norm. Figure 6.1 compares, in the scalar case, the L 2 and L 1 -norms to the norm based on the Huber function for the particular case M = 1. Practical algorithm for solving the systems of equations (6.1) and (6.4), both in the case of linear and nonlinear residual functions are presented in the following section. A second difficulty is that the L 1 norm is not strictly convex, so that the uniqueness is not guaranteed. This difficulty is taken into account in the solution procedure by adding a strictly convex penalization term, for example a L 2 constraint.

A third potential issue with using a dictionary, as opposed to a reduced basis, is the fact that the dictionary may be rank-deficient. One option to address this issue is to perform a Gramm-Schmidt orthogonalization or a rank-revealing QR factorization. A drawback of that approach is that dictionary members are then linearly combined. Alternatively, a regularization term and a random perturbation is here added to the minimization functionals to ensure a system with full rank and a unique solution as follows:

• For L 1 -norm minimization, the functional becomes

min q∈R k r n (Vq) 1 + ν q 2 2 = min q N i=1 |r n i (Vq)| + ν k j=1 q 2 j , n = 1, • • • , N t . (6.5)
• For Huber function minimization, the functional becomes

min q∈R k N i=1 φ M (r n i (Vq)) + ν q 2 2 , n = 1, • • • , N t . (6.6)
6.2. Algorithms. A classical solution to minimizing a linear residual vector in the L 1 -norm is by recasting the problem as a linear program (LP). More specifically, assuming that the residual is linear r n (Vq) = A n Vq + b n with A n ∈ R N ×N and b n ∈ R N , a solution to (6.1) is given by the solution q ∈ R k of the LP min q,s,t

1 T (s + t) s.t. A n Vq + b n -s + t = 0 s ≥ 0 t ≥ 0. (6.7)
Unfortunately, this LP involves k + 2N variables and 3N constraints, including N equality constraints, rendering this approach intractable in the case of model reduction.

Alternatively, the L 1 -norm minimization problem can be solved by Iteratively Reweighted Least Squares (IRLS) [START_REF] Daubechies | Iteratively re-weighted least squares minimization for sparse recovery[END_REF]. This approach proceeds iteratively by solving a sequence of weighted least-squares problem. An advantage of this approach is that its implementation can rely entirely on existing least-squares solvers. Furthermore, its complexity is similar to that of the L 2 -norm minimization problem. The procedure is presented in Algorithm 1 in the case of a nonlinear residual vector. At each iteration l, a weighted least-squares problem is solved, where the weight depend on the current value of the residual vector r l as follows:

W l = diag |r l i | -1 2
Algorithm 1 L 1 -norm minimization by Iteratively Reweighted Least-Squares Input: Residual function r(•) and associated Jacobian J(•), reduced basis V, initial guess q 0 , tolerance for convergence Output: Solution q 1: l = 0 2: Compute r 0 = r(Vq 0 ) and

Z 0 = J(Vq 0 )V 3: while l = 0 OR ∆q l-1 | 1 > (1 + q l-1 1 ) do 4:
Compute the weights

W l = diag |r l i | -1 2 5:
Solve the weighted least-squares problem

∆q l = argmin y W l Z l y + W l r l 2 2 6:
q l+1 = q l + ∆q l 7:

Compute r l+1 = r(Vq l+1 ) and Z l+1 = J(Vq l+1 )V 8:

l = l + 1 9: end while 10: q = q l Similarly, minimization of the Huber function can also be done by an IRLS procedure, as described in Algorithm 2. The procedure only differs from its L 1 -norm counterpart by the choice of weights. In the present work, the following choice of weights is proposed for a given residual vector r l

W l = diag δ(|r l i | < M ) + M |r l i | -1 2 δ(|r l i | ≥ M ) . (6.8)
Furthermore, the parameter M is computed as

M = 2 max(1, max(|r l i |)) (6.9)
with 2 = 10 -6 , as it has been found to be a robust choice across different applications.

Algorithm 2 Huber function minimization by Iteratively Reweighted Least-Squares Input: Residual function r(•) and associated Jacobian J(•), reduced basis V, initial guess q 0 , tolerance for convergence Output: Solution z

1: l = 0 2: Compute r 0 = r(Vq 0 ) and Z 0 = J(Vq 0 )V 3: while l = 0 OR ∆q l-1 | 1 > (1 + q l-1 1 ) do 4:
Compute the weights

W l = diag δ(|r l i | < M ) + M |r l i | -1 2 δ(|r l i | ≥ M ) 5: Let M = 2 max(1, max(|r l i |)) 6:
Solve the weighted least-squares problem

∆q l = argmin y W l Z l y + W l r l 2 2 7:
q l+1 = q l + ∆q l 8:

Compute r l+1 = r(Vq l+1 ) and Z l+1 = J(Vq l+1 )V 9:

l = l + 1 10: end while 11: q = q l 7. Numerical applications. 7.1. Nonlinear unsteady problems. 7.1.1. Unsteady Burgers' equation: Application of the error estimator of proposition 5.1. In [START_REF] Abgrall | Robust model reduction by l?1-norm minimization and approximation via dictionnaries: application to non linear hyperbolic problems[END_REF], we have applied the methodology presented above to the Burgers' equation which is the simplest non linear PDE. The initial condition is given by

u 0 (x, µ) = µ sin(2x) + 0.1, µ ∈ [0, 1]
This initial condition is chose such that a shock is generated in a finite time and this shock is unsteady for µ = 0 since its velocity us σ µ = 0.6µ. The dictionary is produced for the parameters {0, 0.2, 0.4, 0.6, 1.0}(k = 5) with a classical second order finite volume with MUSCL extrapolation. The CFL condition is always 0.5 and the grid has 100 mesh points. In these example, the target solution is for parameter µ = 0.5. A typical result is displayed in figure 7.1 The minimisation procedure is the L 1 minimisation (6.5) with ν = 10 -5 . As one can see, the results are good until the shock is formed and then more or less severe artefacts are produced. However, the solution is non oscillatory as expected. The quality of the results is independant of the L 1 minimization procedure. Nevertheless, proposition 5.1 suggest that some "greedy" selection of the parameter can be used in order to improve and control the error that is being done. Let us consider the same problem with the same target parameter µ = 0.5 and we set a precision = 0.05 for the estimate (5.7). This value of is taken only for illustration purpose and is not a 'magic' number. Doing so, one possibility for the dictionary is P = {0.47, 0.48, 0.52, 0.53}(k = 4) which satisfies the condition max µi∈P J u 0 (µ i ) -u 0 (µ ) ≤ where the functional J is defined as in the example of (4.17). We display the solutions obtained by L 1 -norm by LP minimization procedure for t = π 4 < 1, t = π 2 and t = π in Figures 7.2. We can observe that with the help of the estimate (5.7) we obtain a local dictionary. With large dictionaries, it may be computationally expensive to consider all dictionary members as potential basis vectors for a given value τ = (t , µ ) of time and parameters. Instead, a local dictionary approach can be considered by restricting the dictionary members considered for τ in a neighborhood of the time and/or parameter domains. The local dictionary is then defined as

V(τ ) = [w(τ 1 ), • • • , w(τ r )] , τ i ∈ T (t ) × P(µ ) ⊂ [0, T ] × P. (7.1)
In such a dictionary approach, restricting dictionary members in the time and/or parameter domains is straightforward, unlike the case of pre-computed reduced bases, for which an a priori partitioning of the parameter domains are necessary [START_REF] Dihlmann | Model reduction of parametrized evolution problems using the reduced basis method with adaptive time-partitioning[END_REF][START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF].

∂ ∂t   ρ ρu E   + ∂ ∂x   ρu ρu 2 + p u(E + p)   = 0, (7.2)
for which U = (ρ, ρu, E) T and the pressure is given by for the density, two for the velocities and one for the total energy on which the reduced solutions are (independantely) expanded. The minimisation method is the straight L 1 minimisation. In order to ilustrate the technique we have chosen a NACA012 profile with subsonic and transonic conditions. In the first case, the inflow mach number is M ∞ = 0.65 and the angle of attack may change. The dictionary P is constructed by sampling the parameters -3.0 • , -2.0 • , -1.0 • , 1.0 • , 2.0 • , 3.0 • } and the solution is sought for the predictive case µ = 0.0 • . We illustrate the elements in the dictionary (Figure 7.4) and then we obtain a ROM solution which is comparable with the exact numerical solution (Figure 7.5).

p = (γ -1) E - 1 2 ρu 2 (7.3)
In the second case, we consider that a shock exist (M ∞ = 0.85) in order to illustrate that our method can also deal with this kind of problems. The dictionary P is constructed by sampling the parameters {-3.0 • , -2.0 • , -1.2 • , 1.0 • , 2.0 • , 3.0 • }. Note that the symmetry is intentionally broken. The solution is sought for the predictive case µ = 0.0 • . As in the first case, we illustrate the elements in the dictionary, in order to show that they are different (Figure 7.6).

We obtain the following ROM solution using the L 1 -norm minimization onto the convex hull of the dictionary (Figure 7.7). This one proves to be more robust. There is some discrepancy between the ROM solution and the exact numerical solution because the problem is very sensitive to the angle (this can be seen from the dictionary elements (Figure 7.6)). If we increase the number of elements in the dictionary, the ROM solution will be more similar to the exact numerical one, but in this case there is no simple strategy to obtain error bounds as in the scalar case. 8. Computational cost. The minimisation procedure consists in looking for the minimum of functionals of the type (4.18) where all the degrees of freedom describing the solution appears. This is a challenge since in general the algorithms to solve this kind of minimisation procedure are much more expensive that those of the least square type. In all the calculations presented here, we have used a slightly different approach. Instead of using all the degrees of freedoms, we have used a small subset of them. For example, instead of the functional (4.18b), we minimize J(r(w n , w n+1 )) = where I is a small subset of the set of degrees of freedom. Of course the question is how to choose this set.

In the dictionary approach a lot of information is encoded. In the case of a subsonic solution, the problem to solve is essentially elliptic hence two different degrees of freedoms are linked together: this means that if one consider the grid points M 1 and M 2 , the flow field at M 1 has some knowledge of what occur at M 2 . In the case of a transonic flow, the subsonic and supersonic pockets are somehow disconnected.

Using this rational, we have selected randomly points in the mesh, taking into account an a priori knowledge of the location of the supersonic pockets when they exist. We have used only about 100 points in each case and made the minimisation. Several stencils have been chosen and the results seem to be un-sensitive to the choices. In addition, in the 1D case where the full grid can be used, we have not experienced any change in the solution. Of course this heuristic needs to be further studied and quantified, this will be the topic of a future paper. 9. Conclusions, perspectives. We have presented a general framework to approximate the solution of steady and unsteady hyperbolic problems. The solution can be smooth or discontinuous, and in the unsteady case, moving wave may exist. Starting from any standard scheme (explicit or implicit), the reduced order solution is obtained at each time step (or each iteration in the implicit case) from a minimization problem in the L 1 norm. We give a sufficient condition to be able to solve the problem, and discuss the practical aspects of the method. It is illustrated by several examples dealing with linear and non linear problems, scalar and systems in one and two space dimensions. A rough error estimate is given, relating the successive projections and the initial solution.

As a future work, the method efficiency can be improved by using hyper-reduction. Hyperreduction will hopefully solve the cost and the memory issues.

Fig. 6 . 1 .

 61 Fig. 6.1. Comparison of the L 2 , L 1 and Huber function (M = 1) norms

Fig. 7 . 1 .

 71 Fig. 7.1. Unsteady Burgers' equation: predicted solutions at target parameter µ = 0.5 at t = π 4 , π 2 and π.

Fig. 7 . 2 .

 72 Fig. 7.2. Unsteady Burgers' equation: predicted solutions at target parameter µ = 0.5 using the dictionary obtained from the estimate (5.7)

7. 1 . 2 .

 12 Euler equations. The one-dimensional Euler equations are considered on Ω = [0, 1]

Fig. 7 . 3 .

 73 Fig. 7.3. Results of the shock tube problem (7.6).

Fig. 7 . 4 .

 74 Fig. 7.4. The dictionary elements constructed from the parameters {-3.0 • , -2.0 • , -1.0 • , 1.0 • , 2.0 • , 3.0 • } for Mach=0.65

  j∈I r(w n , w n+1 ) j + ν j∈I

Fig. 7 . 5 .

 75 Fig. 7.5. The ROM solution and the exact numerical solution for Mach= 0.65 and AoA=0.0 •

Fig. 7 . 6 .

 76 Fig. 7.6. The dictionary elements constructed from the parameters {-3.0 • , -2.0 • , -1.2 • , 1.0 • , 2.0 • , 3.0 • } for Mach=0.85

so that in particular ∀(x, 0) ∈ R d × R + , ϕ(x, 0) = 0.

Acknowledgements. The first author was funded in part by the MECASIF project (2013-2017) funded by the French "Fonds Unique Interministériel" and SNF grant # 200021_153604 of the Swiss National Foundation. The second author has been supported in part by the SNF grant # 200021_153604 of the Swiss National Foundation RA would also like for several very insightful conversations with Y. Maday, LJLL, Université Pierre et Marie Curie. D. Amsallem is also warmly acknowledged.

with γ = 1.4. For this problem, we have taken a second order finite volume scheme with MUSCL extrapolation on the characteristic variables. The limiter is minmod on all waves. This problem is parametrized by the initial conditions U 0 (x; µ). In [START_REF] Abgrall | Robust model reduction by l?1-norm minimization and approximation via dictionnaries: application to non linear hyperbolic problems[END_REF], we had chosen to initialise one particular combination of the Sod shock tube problem and the Lax one. Here we choose another combination leading to solution with a richer structure in order to illustrate the strengh of the method.

The state U Sod (x) is defined by the primal physical quantities:

125 otherwise, u = 0.0 p = 1.0 if x ≤ 0.5, 0.1 otherwise, (7.4) and U Lax (x) defined by

The Sod condition presents a fan, followed by a contact and a shock. For the density and the pressure, the solution behaves monotonically, and the contact is moderate. The Lax solution has a very different behavior and the contact is much stronger. In this example, the initial condition are parametrized for µ ∈ [0, 1] as

• else if x < 0.5 + 0.2µ then

As suggested in [START_REF] Abgrall | Robust model reduction by l?1-norm minimization and approximation via dictionnaries: application to non linear hyperbolic problems[END_REF], we have reconstructed the density, momentum and total energy independently, and here we have used the L 1 minimization. The target solution if for µ = 0.6, the dictionary is obtained for µ ∈ {0, 0.2, 0.4, 0.5, 0.8, 1}. The final time of the simulation is T f in = 0.16. It is obtained again with 100 grid points and a time step of 0.001. The results are given on figure 7.3, where we have superimposed the reduced solution and the solution obtained by the full order method without any approximation. As it can be seen, the results are very close to each other and this confirms what was written in [START_REF] Abgrall | Robust model reduction by l?1-norm minimization and approximation via dictionnaries: application to non linear hyperbolic problems[END_REF].

7.2. Nonlinear steady problems: a two dimensional example. Te extension to multidimension is straightforward. We have started from a a code using a second order oscillation free method for solving the Euler equations. A description of this method can be found, for example, in [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF]. Note that this specific choice has no impact on the Reduced Order Model, since this algorithm is coded in Python on top of the CFD code which is called as a black box: any other CFD method would do the job.

The minimization is done on the mass, momentum components and total energy: we introduce 4 sets of independent parameters which are the expansion coefficients for the dictionaries. In other words, we first run the CFD code to get a finite number of CFD solutions. Each solution is described as vector of state variables (ρ, m x , m y , E). From this we form foru disctionaries, one


Emilio Jesús

Gallego Arias

SerAPI: Machine-Friendly, Data-Centric Serialization for COQ

Keywords: serialization, interactive protocols, user interfaces, program verification, theorem prover implementation

We present SerAPI, a library and protocol for machine-friendly communication with the COQ proof assistant. SerAPI is implemented using Ocaml's PPX pre-processing technology, and it is specifically targeted to reduce implementation burden for tools such as Integrated Development Environments or code analyzers.

SerAPI tries to address common problems that tools interacting with COQ face, providing a uniform and data-centric way to access term representations, proof state, and an extended protocol for document building.

SerAPI is work in progress but fully functional. It has been adopted by the jsCoq and PeaCoq Integrated Development Environments, and supports running inside a web browser instance. For the near future, we are focused on extending the document protocol and providing advanced display abilities to clients.

Introduction

Historically, the interface to the outside world of the COQ proof assistant has been its toplevel, a Read Print Eval Loop. COQ's REPL is targeted at human users, which issue proof building commandsto progress towards proving the "current goal" -intermixed with control commands and information requests. The limitations that a command-based model suppose for proof development quickly become apparent, and proof development environments such as Proof General [START_REF] Aspinall | Proof General: A Generic Tool for Proof Development[END_REF] appeared to provide a more convenient document-based model on top of the REPL.

Nobody would deny the success of this approach, however, interacting with the REPL interface brings some challenges and limitations. Interpreting prover output is complicated, especially in the presence of COQ's extensible parsing and printing mechanism. Similarly, information internal to the theorem was not exposed in a systematic way. Queries for such data -often needed by the tools -were performed by custom-purpose commands, often with semi-adhoc syntax and output format. The need to modify the core command language of COQ implied that new functionality was costly to implement, requiring a new release of system.

Until recently, the only alternative to the REPL was the use of plugins, written in Ocaml. Plugins provided full access to COQ's internals, however they were hard to distribute and restricted to Ocaml-friendly environments.

Fortunately, the situation has recently changed with the introduction of the State Transactional Machine [START_REF] Barras | Asynchronous Processing of Coq Documents: From the Kernel up to the User Interface[END_REF] in COQ 8.5 and a new, XML-based, interactive protocol. The new protocol was inspired by related proof IDE efforts [START_REF] Wenzel | Asynchronous User Interaction and Tool Integration in Isabelle/PIDE[END_REF], and several new developments are based on it, including an experimental port of Proof General [START_REF] Steckler | Proof General with XML Protocol Support[END_REF][START_REF] Faithfull | Coqoon: An IDE for interactive proof development in Coq[END_REF][START_REF] Wenzel | PIDE as front-end technology for Coq[END_REF].

Our particular first contact with COQ's STM came when we started the development of jsCoq [START_REF] Gallego Arias | jsCoq: towards hybrid theorem proving interfaces[END_REF], a port of COQ to the browser platform. jsCoq started as a Ocaml application, then transpiled to JavaScript by the js_of_ocaml [START_REF] Vouillon | From bytecode to JavaScript: the Js_of_ocaml compiler[END_REF] tool. The event-based nature of the STM API proved to be well suited to the browser environment, allowing the completion of a first prototype. However, we soon realized that the next step for jsCoq -running the prover process inside a web worker thread -would require the use of a protocol. The obvious choice would have been to use the existing XML protocol, however after careful consideration we chose to build our own extension. The main reasons were:

• jsCoq was developed using the STM Ocaml's API. However, there is significant impedance with the way the state machine is exposed through the XML protocol. Adapting would have required a complete rework of jsCoq. • Some of the ambitions jsCoq goals -advanced term display, asynchronous communication -were not well supported by the XML protocol. • The amount of serialization boilerplate was high. COQ exposes hundredths of datatypes and the XML protocol was designed prior to the introduction of the Ocaml's PPX meta programming facilities. The opportunity of designing a new system based on meta programming seemed worth the effort to us. • The final factor was the realization that the existence of such a tool could be potentially useful to other developers and users. This way, SerAPI was born.

SerAPI overview

The SerAPI philosophy is to use ML-style datatypes as the schema specification, using PPX to automatically generate a machinefriendly serialization. In our opinion, the cost of maintaining a separate data representation is too high. This has the effect to link SerAPI data representation to Ocaml's API stability, however care is taken in order to ensure that the more advanced representation are opt-in. For example, a client may request to receive all objects pretty printed as strings.

Once we have encoded the required datatypes, the API we want to expose are reified to -usually small -domain specific languages, with their corresponding interpreter serving as the link between the reified version and the actual COQ procedures. The DSL representation is serialized in the same way.

We briefly survey the components of the current distribution.

The Serialization Library The base component of SerAPI is SerLib, the serialization library. SerLib provides an overlay of COQ's modules, extending them with the corresponding serialization functions.

Serialization is performed by the ppx _ sexp _ conv package, with help from ppx _ import. For example, we declare the vernac _ expr datatype to be serializable to sexps with: type vernac _ expr = [%import: Vernacexpr.vernac _ expr] [@@deriving sexp] ( * will generate: * ) val vernac _ expr _ of _ sexp : Sexp.t → vernac _ expr val sexp _ of _ vernac _ expr : vernac _ expr → Sexp.t

The current distribution serializes 328 COQ datatypes by including this boilerplate in the corresponding modules. For the few private datatypes of interest -such as constr -we provide our own view and translation functions.

The SerAPI Protocol The second component is the SerAPI control protocol, an extension of COQ's STM API which is used to build and execute proof scripts. Given the IDE origins of SerAPI, it was natural to provide IDE support first. The current control protocol is specified by a DSL and its full description goes beyond the scope of this abstract. A peek on the definitions is: Several changes over the standard STM model have been performed, mainly to facilitate the work of the clients and to be more robust w.r.t. race conditions. Events or feedback from COQ are forwarded to the client unchanged.

type control _ cmd = |
The SerComp Coq Processor The last component is SerComp, an experimental batch processor for COQ files, supporting data aggregation and exporting of COQ files to a machine-friendly format.

A Data-Centric View

A central goal of SerAPI is to expose information known to COQ using a principled, "database-inspired" API. The SerAPI "Query" protocol is completely separated from the control protocol, and still in heavy development. An important feature is that changes in the protocol aggregation or query functions don't require changes in COQ at all. Thus, we enjoy freedom to experiment without interfering with core parts of the system.

The base datatype is the coq _ object datatype that simply adds a tag to COQ types. Then, each tag can be queried for, with common options including output format, paging, and filtering: A key difference with the current XML protocol is that SerAPI won't produce any non-control information unless queried. This implies a slightly different control flow that we prefer.

Current Status and Goals

SerAPI is open source and available at [START_REF] Arias | Coq Serializable Protocol (SerAPI)[END_REF]. The current prototype is fully functional, and we provide a functional web version1 . Experience with jsCoq gives us moderate confidence on the robustness of the approach; our main goal at this stage is to gather more interest and get feedback from users, helping us to converge on the design. The two most important tasks for the near future are:

• Finalize the control protocol: Experience implementing IDEs on top of SerAPI has led us to extend the original STM protocol in a few ways. However, some tricky details remain to be finalized, such as the definitive error handling procedure. We hope to complete a first draft of the extended STM protocol soon, and hopefully some of the ideas would be suitable for incorporation in COQ upstream. • Principled data API: COQ internally stores many useful data, however there is no principled interface for access. We believe however that for the most part, the ML API is not so far from it.

We hope we can contribute to improve and uniformize COQ's upstream API too, covering some important use cases such as efficient search and completion.

A medium term goal is to improve the COQ's printing system so tools can offer a richer interpretation of notations (see for example [START_REF] Pit-Claudel | Company-Coq: Taking Proof General one step closer to a real IDE[END_REF]). Some preliminary work is already in place but more research is needed before completing a fully-functioning prototype.

  type coq _ object = | CoqConstr of Constr.constr | CoqExpr of Constrexpr.constr _ expr | CoqTactic of KerName.t * Tacenv.ltac _ entry type query _ opt = { preds : query _ pred sexp _ list; limit : int sexp _ option; sid : Stateid.t [@default current _ state()]; pp : print _ opt [@default {pp _ fmt = PpSer}]; } type query _ cmd = Option | Tactics | Goals | ...

https://x80.org/rhino-hawk/

Acknowledgements: We would like to thank Clément Pit-Claudel, Valentin Robert, and Enrico Tassi for very valuable feedback and testing.