Fabien Picarougne
email: fabien.picarougne@univ-nantes.fr

Lina

Aurélien Milliat
email: aurelien.milliat@univ-nantes.fr

Real Time 3D Geometry Correction for Folding Screens or Projectors with Distorting Lenses

Keywords: 3D, geometry correction, curved surface, folding screen, real time Concepts:, Computing methodologies → Rendering, 3D imaging, •Hardware → Displays and imagers

In this article, we describe a new method for real-time display, in a geometrically correct way, of a 3D scene on curved surfaces. This method differs from existing solutions in the literature by allowing its application with folding screens or projectors with distorting lenses. Our algorithm is not limited to a particular shape of the display surface and takes the position of the user into account to display an image that is geometrically correctly perceived from the user's viewpoint. The projection process of 3D objects is divided in three phases. The first two steps are independent from the 3D scene and act as a buffer that stores particular values and accelerate calculations of the third step. The execution of the later may then be performed in linear time in the number of vertices of the geometry to be displayed.

Projection on curved surface

Virtual reality equipment (VR) are increasingly used for industrial, recreational and educational. After several years of evolution, 3D and highly immersive image projection systems for multi-walled stereoscopic environments have reached maturity. However, developments are still active, and users want now to view immersive content on complex surfaces. For example, large flight simulators traditionally use cylindrical projection screens or dome projection screens, planetariums and Omnimax theaters use hemispherical screens. However, the projection of images on these non-planar surfaces requires correcting the geometry of the image to be visually correct from the perspective of the user.

Since the 90s, several solutions have been proposed to solve the problem of distortion with real-time correction. In the literature, we distinguish four main approaches for displaying content on non planar surfaces. 1) Ray tracing. But this method is extremely costly in calculation time and yet cannot be real-time for high resolutions. 2) Image rendering in two passes. It consists in generating and distort an image by mapping it on a grid that takes the form of the displayed surface. But this method requires a good knowledge of the parameters of the projection system in use in order to parametrized the mapping grid. 3) Projection mapping technique [START_REF] Raskar | The office of the future: A unified approach to image-based modeling and spatially immersive displays[END_REF]]. This rendering technique also acts in two passes, but instead of parameterizing a grid to model the displayed surface, it renders into a texture the 3D scene to display as seen from the viewer's position and projected it onto a 3D object that model the surface to be displayed, as if a projector was located at the viewer position. But this method makes the hypothesis that the image is displayed by a classical projector with a standard central projection: not applicable with distorting lenses or foldable screens. 4) Transfer on quadric surface. Here, we directly modify the coordinates of the 3D objects to display, according to the curvature of the display surface and the position of the user. This technique requires only one rendering pass but current implementations therefore limit the kind of display surface used to quadrics in order to simplify the final calculation.

We propose in this article to improve and generalize the use a curved display surfaces for viewing 3D objects by using classical GPU pipeline. Our method takes into consideration a display surface of any shape and any type (foldable screen, projectors potentially having distorting lenses. . . ) together with the position of the user and performs a real-time geometric correction. These data can also be updated dynamically, keeping the calculation of the geometric distortion of the image in real-time.

Description of the projection method

The originality of our method is to not taking into account any hypothesis on the way pixels are generated. Our projection process is mainly divided in three phases. The first processes the data describing the geometry of the surface and is executed each time the geometry of the projection surface changes. The second phase takes in consideration the position of the user in space and has to be executed each time the viewing position of the user change. The third phase focuses on the calculation of the final image generation and has to be executed on each image rendering.

In more details, in the first phase we acquire the geometry of the display surface. A huge literature exists to solve this problem and we do not claim any contribution on this point. We then associate a 3D spatial coordinate {x, y, z} to each displayed pixel {u, v} and build a 3D model M d of the display surface from these data. This mesh is simply built with a 2D Delaunay triangulation using the {u, v} coordinates. Then, each vertex Vi of the mesh is defined with its spatial coordinates {xi, yi, zi} in world coordinates and with a 2D vector {ui, vi} associated to the coordinates of the corresponding pixel of the display device in the device coordinates. It is important to note that this parametrization allows a complete independence in terms of the organization of pixels on the display surface. Our method allows the use of a wide variety of image display devices, ranging from curved LCD display to projectors fitted (or not) with distorting lenses.

The second phase of the algorithm consists in rendering in a texture buffer the mesh M d with a classical pinhole camera from the viewpoint of the user. To do that, we need to calculate projection parameters (PL) that optimize the definition of the stored data. This is achieved with a simple two-pass algorithm that we will not detail here. The 2D obtained image stores for each fragment of the mesh the interpolation of the {u, v} parameters associated with each ver-tex Vi of the mesh. In some words, it stores for each fragment of this projection, the coordinates of the corresponding pixel that illuminates the display surface. This image, called pixel lookup table in the following, is the result of this step and will be used in the last phase of the algorithm. It can be considered as the main distinctive contribution of our work. Using very simple 3D rendering techniques such as Z-buffer, it allows us to avoid making important calculations in the last phase of our method.

The last phase consists in calculating the projection of the geometry to be displayed on the display surface. We first project these data points on a 2D flat plane, keeping projection parameters PL. And for each vertex V G i of the geometry, we get the corresponding parameter {ui, vi} in the pixel lookup table. The computational cost of the search operation into the lookup table is constant. Therefore the global process of projection is linear as a function of the number of vertex of the geometry. Finally, we modify the coordinates

{x G i , y G i , z G i } of each vertex V G i
of the geometry to display by setting {ui, vi} to {xi, yi} and distance between the user's eye and V G i to the z coordinate. These new coordinates allow obtaining the final image to be displayed by performing a simple orthographic projection.

An additional step is required because, generally, a 3D scene is described by considering the context of planar projection. Therefore, it is common to model planar surfaces with very few points. However, in the context of a projection on a curved surface, wherein only the positions of the points defining the geometry are "distorted", visual artifacts become significant. In order to correct this distortion, we subdivide (tessel) automatically the geometric primitives forming a 3D object according to the length of each side of the primitive, once projected on the display surface. This optimizes the computational time and the quality of the deformation of the geometrical surfaces that are not initially intended to be curved.

Results

In order to evaluate our method we have setup a rig composed by a dome of two meters of diameter, with a standard projector (1440x1080 resolution) mounted with a deformation lens. All our tests have been done on a HP Z800 workstation working on Windows 7 64 bits with an Intel Xeon E5620 @2.40Ghz CPU with 8.0Go RAM and an nVidia Quadro K5000 GPU (driver version 355.97).

We have measured the computational cost of the last two phases of our method. The second phase is needed only if the user position change, and the third phase has to be executed on each rendering frame. In order to see the linearity of the computational cost of our method in function of the number of vertices, we have compared its performances with respect to standard rendering algorithm on different scene complexities: Scene #1 to #4 with respectively 3.088. 728, 3.960.142, 5.961.130 and 14.230.254 triangles primitives.

Figure 1 shows the framerates achieved by activating or deactivating the various phases of our method (our algorithm is completely inactive when the two tested phases are disabled) and the evolution of framerates with respect to scene complexities. These results prove that the third phase of our algorithm is computed in a linear time of the number of vertices of the 3D scene and is not scene dependent. Its cost is about 15% of the global computational process in the constraints of our evaluation. The second phase depends only on the number of vertexes of the mesh describing the rendering surface. This number in constant regardless of the complexity of the scene assessed. Then, its cost decreases with increasing the complexity of the scene. In order to evaluate the cost of the tessellation step, we have compared the results of our method (that uses a geometry and a tesselation shader) with a very simple graphical pipeline that just executes a vertex and a fragment shader and an other one that adds a geometry shader that does nothing except transferring data. This evaluation was performed only on the most complex scene (eg Scene #4) and the results are given in table 1. With our implementation, our method reduces the framerate of 55%, considering a rendering pipeline that needs a geometry shader stage. On a side note, the geometry shader stage, with a simple primitives forward, cost about 60% of performance in a rendering pipeline.

Conclusion

We presented in this paper a new real-time method to deform a 3D model displayed on a curved surface so that it is geometrically correctly seen from a user's viewpoint. This method differs from existing solutions by allowing its application with folding screens or projectors with distorting lenses.

The algorithm proceeds in three stages and the final distortion of the image is parametrized only by the position of the user and by the description of the display surface. The system can therefore be used in combination with projection optics (fisheye for example) or coarse display surfaces, thus limiting the cost of manufacturing immersive systems. These input parameters can also be easily obtained and injected in real-time into our system.

Figure 1 :

 1 Figure1: Graph comparing the framerates obtained by the using the different phases of our method with respect to scenes of different complexity. Results are given on an average of 1.000 frames. SP and TP respectively correspond to the second and third phase of our method.

Table 1 :

 1 Comparison of framerates obtained by very simple and standart rendering pipelines (one using a simple vertex shader and a second using a simple geometry shader) and by our algorithm, on Scene #4 on an average of 1.000 frames. SP and TP respectively correspond to the second and third phase of our method.

		std.	std. +geo	w/o SP w/o TP	w SP w/o TP	w/o SP w TP	w SP w TP
	fps	53.39	20.97	11.33	11.12	9.548	9.507
	%	471%	185%	100%	98.1%	84.3%	83.9%


A D Agaltsov
email: alexey.agaltsov@polytechnique.edu

G M Henkin
email: henkin@math.jussieu.fr

Explicit reconstruction of Riemann surface with given boundary in complex projective space 1

Keywords: Riemann surface, reconstruction algorithm, Burgers equation, Cauchy-type formulas

In this paper we propose a numerically realizable method for reconstruction of a complex curve with known boundary and without compact components in complex projective space.

Introduction

Let CP 2 be the complex projective plane with homogeneous coordinates (w 0 : w 1 : w 2 ). Let X ⊂ CP 2 be a complex curve with rectifiable boundary γ = ∂X. Without loss of generality, we suppose that the following conditions of general position hold:

(0 : 1 : 0) ∈ X, w 0 | γ = 0.
Put C 2 = {w ∈ CP 2 : w 0 = 0} with coordinates z 1 = w 1 w 0 , z 2 = w 2 w 0 . For almost all ξ = (ξ 0 , ξ 1 ) ∈ C 2 the points of intersection of X with complex line C 1 ξ = {z ∈ C 2 : ξ 0 + ξ 1 z 1 + z 2 = 0} form a finite set of points z (j) 1 (ξ), z (j) 2 (ξ) = h j (ξ), -ξ 0 -ξ 1 h j (ξ) , ξ = (ξ 0 , ξ 1 ), j = 1, . . . , N + (ξ).

By Darboux lemma [START_REF] Darboux | Théorie des surfaces[END_REF][START_REF]Chaînes holomorphes de bord donné dans CP n[END_REF] functions {h j } satisfy the following equations ∂h j (ξ) ∂ξ 1 = h j (ξ) ∂h j (ξ) ∂ξ 0 , ξ = (ξ 0 , ξ 1 ), j = 1, . . . , N + (ξ),

which are often called the shock-wave equations, the inviscid Burgers equations or the Riemann-Burgers equations. In this interpretation ξ 1 is the time variable and ξ 0 is the space variable.

The following Cauchy-type formula from [START_REF]Chaînes holomorphes de bord donné dans CP n[END_REF] plays the essential role in the reconstruction of X from γ:

G m (ξ) def == 1 2πi γ z m 1 d(ξ 0 + ξ 1 z 1 + z 2 ) ξ 0 + ξ 1 z 1 + z 2 = N + (ξ) j=1 h m j (ξ) + P m (ξ), m ≥ 1, (2)
where N + (ξ) = |X ∩ C 1 ξ | and P m (ξ 0 , ξ 1 ) is a polynomial of degree at most m with respect to ξ 0 at fixed ξ 1 . In addition, P 0 (ξ 0 , ξ 1 ) = -N -, where

N -= |X ∩ CP 1 ∞ |, CP 1
∞ = {w ∈ CP 2 : w 0 = 0}, and

P 1 (ξ 0 , ξ 1 ) = - N - k=1 a k ξ 0 + b k 1 -a k ξ 1 , if h k (ξ 0 , 0) ∼ a k ξ 0 + b k + O(ξ -1 0 ) at ∞. (3)
In particular, it follows from (2) that:

G 0 (ξ) = 1 2πi γ d(ξ 0 + ξ 1 z 1 + z 2 ) ξ 0 + ξ 1 z 1 + z 2 = N + (ξ) -N -. (4)
Let π 2 : C 2 → C be the projection onto the second coordinate: π 2 (z 1 , z 2 ) = -z 2 . We have that C \ π 2 γ = ∪ L l=0 Ω l , where Ω l≥0 are connected and Ω 0 is unbounded. From the definition of N ± it follows that N + (ξ 0 , 0) = N -, ξ 0 ∈ Ω 0 .

(5)

Assume that complex curve X does not contain compact components without boundary, or equivalently, satisfies the following condition of minimality: for any complex curve X ⊂ CP 2 such that ∂ X = ∂X = γ and for almost all ξ ∈ C 2 we have

| X ∩ C 1 ξ | ≥ |X ∩ C 1 ξ |. ( * )
Condition of minimality ( * ) is a condition of general position and is fulfilled for X if, for example, every irreducible component of X is a transcendental complex curve. Note that from theorems of Chow [START_REF] Chow | On compact complex analytic varieties[END_REF] and Harvey-Shiffman [START_REF] Harwey | A characterization of holomorphic chaines[END_REF] it follows that an arbitrary complex curve X ⊂ CP 2 satisfying ∂ X = ∂X admits the unique representation X = X ∪ V , where X is a curve satisfying ( * ), and V is a compact algebraic curve, possibly, with multiple components.

The main result of [START_REF] Dolbeault | Surfaces de Riemann de bord donné dans CP n[END_REF] gives a solution to the important problem of J. King [START_REF] King | Open problems in geomtric function theory[END_REF], when a real curve γ ⊂ CP 2 is the boundary of a complex curve X ⊂ CP 2 . Let γ ⊂ C 2 . Then γ = ∂X for some open connected complex curve X in CP 2 if and only if in a neighborhood W ξ * of some point ξ * ∈ C 2 one can find mutually distinct holomorphic functions h 1 , . . . , h N (ξ * ) satisfying shock-wave equations (1) and also the equation

∂ 2 ∂ξ 2 0 G 1 (ξ 0 , ξ 1 ) - p j=1 h j (ξ 0 , ξ 1 ) = 0, ξ = (ξ 0 , ξ 1 ) ∈ W ξ * .
In this work in development of [START_REF] Dolbeault | Surfaces de Riemann de bord donné dans CP n[END_REF][START_REF]Chaînes holomorphes de bord donné dans CP n[END_REF] we propose a numerically realizable algorithm for reconstruction of a complex curve X ⊂ CP 2 from the known boundary and satisfying the condition of minimality. This algorithm permits, in particular, to make applicable the result of [START_REF] Henkin | On the explicit reconstruction of a Riemann surface from its Dirichlet-to-Neumann operator[END_REF] about principal possibility to reconstruct the topology and the conformal structure of a two-dimensional bordered surface X in R 3 with constant scalar conductivity from measurements on ∂X of electric current densities, induced by three potentials in general position.

Our algorithm depends on the number of points at infinity N -of the complex curve X. It was tested on many examples and admits a simple and complete justification for N -= 0, 1, 2. Despite a cumbersome description for N -≥ 3, we show that, in principle, there are no obstacles for the justification and numerical realization for any N -≥ 0. Moreover, in Theorem 3.2 we propose a method for finding the parameter N -in terms of γ. This makes the algorithm much more applicable.

Cauchy-type formulas and Riemann-Burgers equations

We begin by giving a new proof of the Cauchy type formulas (2) from [START_REF]Chaînes holomorphes de bord donné dans CP n[END_REF], which allows to obtain explicit expressions for functions P m .

Theorem 2.1. Let X ⊂ CP 2 \ (0 : 1 : 0) be a complex curve with rectifiable boundary γ ⊂ C 2 and satisfying ( * ). Suppose that for almost all ξ ∈ C 2 all the points of intersection of X with C 1 ξ have multiplicity at most one. Then the following formulas hold for almost all ξ = (ξ 0 , ξ 1 ) ∈ C 2 :

G m (ξ) = N + (ξ) j=1 h m j (ξ) + P m (ξ), m ≥ 1, (6)
where P m (ξ) is holomorphic in a neighborhood of almost all (ξ 0 , 0), ξ 0 ∈ C \ π 2 (γ), and is polynomial in ξ 0 of degree at most m for any fixed ξ 1 . Furthermore, the following explicit formulas hold:

P m (ξ 0 , ξ 1 ) = N - s=1 m-1 k=0 i 1 +•••+im=k d i 1 w 1 dw i 1 0 (q s ) • • • d im w 1 dw im 0 (q s ) (m -k -1)! d m-k dw m-k 0 ln(ξ 0 w 0 + ξ 1 w 1 + w 2 )| qs - µ 0 s=1 i 1 +•••+im=m d i 1 w 1 dw i 1 0 (q s ) • • • d im w 1 dw im 0 (q s ),
where X ∩ CP 1 ∞ = {q 1 , . . . , q N -}. In particular, if N -= 0 then P m = 0, m ≥ 1.

Proof. Put g = ξ 0 w 0 + ξ 1 w 1 + w 2 and g = g w 0 = ξ 0 + ξ 1 z 1 + z 2 . Consider differential forms ω m def == z m 1 dg g = w m 1 w m 0 w 0 g d g w 0 = w m 1 w m 0 d g g - w m 1 w m+1 0 dw 0 , m = 0, 1, . . . Then G m (ξ) = 1 2πi γ ω m .
Let us compute this integral explicitly. Denote by p j , j = 1, . . . , N + (ξ) the points of intersection of X with CP 1 ξ , and by q s , s = 1, . . . , N -the points of intersection of X with infinity CP 1 ∞ . Denote by B ε j the intersection of X with the ball of radius ε in CP 2 centered at p j and by D ε s the intersection of X with the ball of radius ε centered at q s . The restriction of form ω m on X is meromorphic with poles at points p j and q s . Thus the following equality is valid:

G m (ξ) = 1 2πi γ ω m = N + (ξ) j=1 1 2πi bB ε j ω m + N - s=1 1 2πi bD ε s ω m .
If N -= 0, then the second group of terms is absent. The integral bB ε j ω m can be calculated as a residue at the first order pole:

bB ε j ω m = bB ε j z m 1 d g g - bB ε j w m 1 w m+1 0 dw 0 = bB ε j z m 1 d g g = 2πi h m j (ξ).
Let N -> 0. Computation of integral bD ε s ω 1 will be done in two steps. Let us calculate first bD ε s w m 1 w m+1 0 dw 0 . Consider the expansion of w 1 (w 0 ) into power series in w 0 in the neighborhood of point q s :

w 1 (w 0 ) = w 1 (q s ) + dw 1 dw 0 (q s )w 0 + d 2 w 1 dw 2 0 (q s )w 2 0 + • • • . 4 Note further that w m 1 (w 0 ) = ∞ k=0 i 1 +•••+im=k d i 1 w 1 dw i 1 0 (q s ) • • • d im w 1 dw im 0 (q s )w k 0 .
The coefficient near w m 0 can be presented in the form

bD ε s w m 1 w m+1 0 dw 0 = 2πi i 1 +•••+im=m d i 1 w 1 dw i 1 0 (q s ) • • • d im w 1 dw im 0 (q s ).

Now we can calculate the integral bD

ε s w m 1 w m 0 d g g . Using relation d g = d g
dw 0 dw 0 and expansion of w 1 (w 0 ) into power series in w 0 we obtain:

bD ε s w m 1 w m 0 d g g = bD ε s 1 w m 0 w 1 (q s ) + dw 1 dw 0 (q s )w 0 + d 2 w 1 dw 2 0 (q s )w 2 0 + • • • m d g dw 0 1 g dw 0 = ∞ k=0 bD ε s i 1 +•••+im=k d i 1 w 1 dw i 1 0 (q s ) • • • d im w 1 dw im 0 (q s )w k-m 0 d g dw 0 1 g dw 0 = m-1 k=0 bD ε s i 1 +•••+im=k d i 1 w 1 dw i 1 0 (q s ) • • • d im w 1 dw im 0 (q s )w k-m 0 d g dw 0 1 g dw 0 = m-1 k=0 2πi (m -k -1)! i 1 +•••+im=k d i 1 w 1 dw i 1 0 (q s ) • • • d im w 1 dw im 0 (q s ) lim w 0 →0 d m-k-1 dw m-k-1 0 d g dw 0 1 g .
From here, taking into account the relation d g dw 0 1 g = d ln g dw 0 , we obtain, finally

bD ε s w m 1 w m 0 d g g = 2πi m-1 k=0 i 1 +•••+im=k d i 1 w 1 dw i 1 0 (q s ) • • • d im w 1 dw im 0 (q s ) (m -k -1)! d m-k dw m-k 0 ln(ξ 0 w 0 + ξ 1 w 1 + w 2 )| qs .
It is a polynomial of degree at most m with respect to ξ 0 .

Remark 2.1. We excluded the case (0 : 1 : 0) ∈ X because the point (0 : 1 : 0) can lead to a non-polynomial contribution in ξ 0 in functions P m . Consider, for example, the projective curve X: w 2 0 = w 1 w 2 and the meromorphic 1-form Ω ξ 0 on it:

Ω ξ 0 = z 1 dz 2 z 2 +ξ 0 , ξ 0 ∈ C \ 0, where z 1 = w 1 w 0 , z 2 = w 2 w 0 .
The form Ω ξ 0 is holomorphic on X except the points q 1 = (0 : 0 : 1), q 2 = (0 : 1 : 0) and p ξ 0 = (-ξ 0 : 1 : ξ 2 0 ).

One can see that the residues of Ω ξ 0 at q 1 , q 2 are equal to 0 and 1/ξ 0 , respectively. It follows that the residue at p ξ 0 is equal to -1/ξ 0 . Now, consider the part X of X obtained by cutting off some small neighborhood of the point q 1 . Let γ denote the boundary of X. Choose any ξ = (ξ 0 , 0) such that

X ∩ C 1 ξ = ∅. Then X ∩ C 1 ξ = {(-1/ξ 0 , -ξ 0 )
} and, as far as the residue of Ω ξ 0 at q 1 is zero, we have that

0 = G 1 (ξ 0 , 0) = 1 2πi γ Ω ξ 0 = - 1 ξ 0 + P 1 (ξ 0 , 0), so that P 1 (ξ 0 , 0) = 1/ξ 0 .
We will also use the following result of [START_REF] Henkin | On the explicit reconstruction of a Riemann surface from its Dirichlet-to-Neumann operator[END_REF], which gives an effective characterization of functions h j≥1 of Theorem 2.1.

Theorem 2.2 (Remark 4 to Theorem 3a of [START_REF] Henkin | On the explicit reconstruction of a Riemann surface from its Dirichlet-to-Neumann operator[END_REF]). Let X ⊂ CP 2 \ (0 : 1 : 0) be a complex curve with rectifiable boundary γ ⊂ C 2 and satisfying ( * ). Fix any ξ * 0 ∈ Ω 0 and let W ξ * be a neighborhood of ξ * = (ξ * 0 , 0) in C 2 . Let G k≥1 , h j≥1 be defined as in (2) and (6). Suppose that there exist functions h 1 , . . . , h p , holomorphic in W ξ * and satisfying

∂ 2 ∂ξ 2 0 (G 1 -h 1 -• • • -h p ) = 0 in W ξ * , (7)
∂ h k ∂ξ 1 = h k ∂ h k ∂ξ 0 in W ξ * , k = 1, . . . , p. (8)
Then p ≥ N -. Furthermore, if p = N -:= |X ∩ CP 1 ∞ |, then h 1 , . . . , h N -coincide with h 1 , . . . , h N -in W ξ * (up to order).

Reconstruction algorithm

We now pass to the reconstruction algorithm for a complex curve X ⊆ CP 2 with given boundary ∂X and satisfying the condition of minimality ( * ). Let us consider the cases

N -= 0, 1, 2.
The reconstruction algorithm is based on formulas [START_REF] Harwey | Boundaries of complex analytic varieties[END_REF]. The next theorem permits to find the functions P m of ( 6). We will use the notation

ae kl = 1 2πi γ z k 1 z l 2 dz 2 , k, l ≥ 0. (9)
Theorem 3.1. Let X ⊂ CP 2 , (0 : 1 : 0) ∈ X, be a complex curve with rectifiable boundary γ ⊂ C 2 and satisfying ( * ). Let h j≥1 and P m≥1 be the functions defined in Theorem 2.1. Then the following statements are valid:

1. If N -= 0, then P m = 0 for all m ≥ 1. Besides, G 1 = 0 in a neighborhood of any ξ * = (ξ * 0 , 0) with ξ * 0 ∈ Ω 0 .

2. If N -= 1, then P 1 (ξ 0 , 0) = c 11 + c 12 ξ 0 , where constants c 11 and c 12 satisfy the following identity in ξ = (ξ 0 , 0) ∈ Ω 0 × 0:

c 11 ∂G 1 ∂ξ 0 (ξ) + c 12 ξ 0 ∂G 1 ∂ξ 0 (ξ) + G 1 (ξ) = G 1 (ξ) ∂G 1 ∂ξ 0 (ξ 0 , 0) -∂G 1 ∂ξ 1 (ξ). ( 10
) 3. If N -= 2, then P 1 (ξ 0 , 0) = c 11 + c 12 ξ 0 , P 2 (ξ 0 , 0) = c 21 + c 22 ξ 0 + c 23 ξ 2 0
, where constants c 11 , c 12 , c 21 , c 22 , c 23 satisfy the following identity in ξ = (ξ 0 , 0) ∈ Ω 0 ×0:

ae 10 (c 2 12 + c 23 ) = ∂G 2 ∂ξ 1 -2 ∂G 1 ∂ξ 1 (G 1 -c 11 -c 12 ξ 0 ) + G 1 (c 22 + 2c 23 ξ 0 ) + ∂G 1 ∂ξ 0 • (G 1 -c 11 -c 12 ξ 0 ) 2 -G 2 + c 21 + c 22 ξ 0 + c 23 ξ 2 0 + G 2 1 -2c 11 G 1 -2c 12 G 1 ξ 0 -G 2 • (-c 12 ), (11)
where all the functions are evaluated at point ξ = (ξ 0 , 0). We differentiate the equation ( 6) with respect to ξ 0 and restrict this equation and its differentiated version to Ω l × 0:

Proof. 1. By Theorem 2.1, if N -= 0 then P m = 0 for all m ≥ 1. It also follows from formula (6) that G 1 = 0 in a neighborhood of any ξ * = (ξ * 0 , 0) with ξ * 0 ∈ Ω 0 . 2. It follows from Theorem 2.1 that P 1 (ξ 0 , ξ 1 ) = C 11 (ξ 1 ) + C 12 (ξ 1 )ξ 0 in a neigh- borhood of almost all ξ * = (ξ * 0 , 0) ∈ Ω l × 0,
h 1 (ξ 0 , 0) = G 1 (ξ 0 , 0) -c 11 -c 12 ξ 0 , ∂h 1 ∂ξ 1 (ξ 0 , 0) = ∂G 1 ∂ξ 1 (ξ 0 , 0) -Ċ11 (0) -Ċ12 (0)ξ 0 , ∂h 1 ∂ξ 0 (ξ 0 , 0) = ∂G 1 ∂ξ 0 (ξ 0 , 0) -c 12 , (12)
where ξ 0 ∈ Ω 0 . By (1) there is the equality ∂h 1 ∂ξ 1 (ξ 0 , 0) = h 1 (ξ 0 , 0) ∂h 1 ∂ξ 0 (ξ 0 , 0). We substitute h 1 (ξ 0 , 0), ∂h 1 ∂ξ 1 (ξ 0 , 0), ∂h 1 ∂ξ 0 (ξ 0 , 0) in this equation by their expressions ( 12), and we obtain the equation

∂G 1 ∂ξ 1 (ξ 0 , 0) -Ċ11 (0) -Ċ12 (0)ξ 0 = = G 1 (ξ 0 , 0) -c 11 -c 12 ξ 0 ∂G 1 ∂ξ 0 (ξ 0 , 0) -c 12 . ( 13
)
This equation is valid, in particular, for ξ 0 ∈ Ω 0 . We divide it by ξ 0 and tend ξ 0 → ∞.

As a result, we obtain the equality Ċ12 (0) = -c 2 12 . Taking into account this equality, we can rewrite equation [START_REF] Wermer | The hull of a curve in C n[END_REF] in the form

∂G 1 ∂ξ 1 (ξ 0 , 0) -Ċ11 (0) = = G 1 (ξ 0 , 0) -c 11 -c 12 ξ 0 ∂G 1 ∂ξ 0 (ξ 0 , 0) -G 1 (ξ 0 , 0) -c 11 c 12 . (14)
Taking into account that ξ 0 ∂G 1 ∂ξ 0 (ξ 0 , 0) → 0 as ξ 0 → ∞ and passing ξ 0 → ∞ in (14), we obtain the equality Ċ11 (0) = -c 11 c 12 . Substituting ths explicit expression for Ċ11 (0) into (14), we get [START_REF] Kälberer | Stripe parametrization of tubular surfaces[END_REF].

3. By (1) functions h 1 (ξ) and h 2 (ξ) satisfy the Riemann-Burgers equation in a neighborhood of any ξ * = (ξ * 0 , 0) ∈ Ω 0 × 0, so that the following equalities are valid:

∂(h 1 h 2 ) ∂ξ 1 = h 1 ∂h 2 ∂ξ 1 + ∂h 1 ∂ξ 1 h 2 = h 1 h 2 ∂(h 1 + h 2 ) ∂ξ 0 , (15)
∂(h 2 1 + h 2 2 ) ∂ξ 0 = 2h 1 ∂h 1 ∂ξ 0 + 2h 2 ∂h 2 ∂ξ 0 = 2 ∂(h 1 + h 2 ) ∂ξ 1 . ( 16
) Note that h 1 h 2 = 1 2 h 1 + h 2 2 -1 2 h 2 1 + h 2 2 .
Therefore the system (15)-( 16) is equivalent to the system

∂(h 1 + h 2 ) 2 ∂ξ 1 - ∂(h 2 1 + h 2 2 ) ∂ξ 1 = h 1 + h 2 2 -h 2 1 + h 2 2 ∂(h 1 + h 2 ) ∂ξ 0 , (17)
∂(h 2 1 + h 2 2 ) ∂ξ 0 = 2 ∂(h 1 + h 2 ) ∂ξ 1 . ( 18
)
We substitute the expressions of (6) for h 2 1 + h 2 2 and h 1 + h 2 into (17), (18), using the notations P 1 (ξ 0 , ξ 1 ) = C 11 (ξ 1 ) + C 12 (ξ 1 )ξ 0 , P 2 (ξ 0 , ξ 1 ) = C 21 (ξ 1 ) + C 22 (ξ 1 )ξ 0 + C 23 (ξ 1 )ξ 2 0 . Then, equation (18) restricted to Ω 0 × 0 takes the form

∂G 2 ∂ξ 0 (ξ 0 , 0) -c 22 -2c 23 ξ 0 = 2 ∂G 1 ∂ξ 1 (ξ 0 , 0) -Ċ11 (0) -Ċ12 (0)ξ 0 . (19)
We divide this equation by ξ 0 and tend ξ 0 → ∞. It leads to the equality Ċ12 (0) = c 23 .

Taking this equality into account and passing ξ 0 → ∞ in (19), we obtain the equality Ċ11 (0) = 1 2 c 22 . Next, we substitute the expressions of (6) for h 2 1 + h 2 2 and h 1 + h 2 into (17) and restrict the obtained formula to Ω 0 × 0. It leads to the equality

2 G 1 -c 11 -c 12 ξ 0 ∂G 1 ∂ξ 1 -Ċ11 (0)-Ċ12 (0)ξ 0 - ∂G 2 ∂ξ 1 + Ċ21 (0)+ Ċ22 (0)ξ 0 + Ċ23 (0)ξ 2 0 = = G 1 -c 11 -c 12 ξ 0 2 -G 2 + c 21 + c 22 ξ 0 + c 23 ξ 2 0 ∂G 1 ∂ξ 0 -c 12 . ( 20
)
We divide this equation by ξ 2 0 and tend ξ 0 → ∞. This leads to the equality

2c 12 Ċ12 (0) + Ċ23 (0) = -c 2 12 + c 23 c 12 .
Taking this equality into account, dividing (20) by ξ 0 and passing ξ 0 → ∞, we obtain the equality

2c 11 Ċ12 (0) + 2c 12 Ċ11 (0) + Ċ22 (0) = -2c 11 c 12 + c 22 c 12 .
Using the obtained equalities, one can rewrite (20) in the form

2 G 1 -c 11 -c 12 ξ 0 ∂G 1 ∂ξ 1 -2G 1 Ċ11 (0) + Ċ12 (0)ξ 0 + 2c 11 Ċ11 (0) - ∂G 2 ∂ξ 1 + Ċ21 (0) = = G 1 -c 11 -c 12 ξ 0 2 -G 2 + c 21 + c 22 ξ 0 + c 23 ξ 2 0 ∂G 1 ∂ξ 0 + + G 1 -c 11 2 -2G 1 c 12 ξ 0 -G 2 + c 21 (-c 12 ). ( 21
)
We pass ξ 0 → ∞ and note that the following relations are valid

lim ξ 0 →∞ ξ 0 ∂G 1 ∂ξ 1 = lim ξ 0 →∞ ξ 0 1 2πi γ z 1 dz 1 ξ 0 + z 2 = 1 2πi γ z 1 dz 1 = 0, lim ξ 0 →∞ ξ 0 G 1 = lim ξ 0 →∞ ξ 0 1 2πi γ z 1 dz 2 ξ 0 + z 2 = 1 2πi γ z 1 dz 2 = ae 10 , lim ξ 0 →∞ ξ 2 0 ∂G 1 ∂ξ 0 = -lim ξ 0 →∞ ξ 2 0 1 2πi γ z 1 dz 2 (ξ 0 + z 2 ) 2 = - 1 2πi γ z 1 dz 2 = -ae 10 .
As a result, we obtain Due to the obtained relations, we can express constants Ċij (0) as functions of c ij :

Ċ11 (0) = 1 2 c 22 , Ċ12 (
Substituting these constants to (21), we obtain the third statement of Theorem 3.1.

Complement 3.1. The statement of Theorem 3.1 admits a development for the case N -≥ 3. In this case

P k (ξ 0 , ξ 1 ) = C k1 (ξ 1 ) + C k2 (ξ 1 )ξ 0 + • • • + C k,k+1 (ξ 1 )ξ k 0 , k = 1, . . . , N -. Denote Ċij (0) = ∂C ij
∂ξ 1 (0) and c ij = C ij (0) for i = 1, . . . , N -and j = 1, . . . , i + 1. Let us indicate the following general procedure for finding of constants c ij . Due to the Riemann-Burgers equations (1) the following identities in ξ 0 ∈ Ω 0 hold for k = 1, . . . , N --1:

- ∂G k ∂ξ 1 (ξ 0 , 0) + Ċk1 (0) + Ċk2 (0)ξ 0 + • • • + Ċk,k+1 (0)ξ k 0 = k k + 1 - ∂G k+1 ∂ξ 0 (ξ 0 , 0) + c k+1,2 + 2c k+1,3 ξ 0 + • • • + (k + 1)c k+1,k+2 ξ k 0 .
Taking into account that ∂G k ∂ξ 1 (ξ 0 , 0) → 0 and ∂G k+1 ∂ξ 0 (ξ 0 , 0) → 0 as ξ 0 → +∞ we obtain the equalities

Ċk,m (0) = km k + 1 c k+1,m+1 , k = 1, . . . , µ 0 -1, m = 1, . . . , k + 1.
Due to the Riemann-Burgers equations (1) the following identity in ξ 0 ∈ Ω 0 holds:

∂e µ 0 ∂ξ 1 (ξ 0 , 0) = e µ 0 (ξ 0 , 0) ∂p 1 ∂ξ 0 (ξ 0 , 0), (23)
where functions e k are given by the following formulas:

ke k (ξ 0 , ξ 1 ) = k-1 i=1
(-1) i+1 e k-i (ξ 0 , ξ 1 )p i (ξ 0 , ξ 1 ) + (-1) k+1 p k (ξ 0 , ξ 1 ),

p k (ξ 0 , ξ 1 ) = G k (ξ 0 , ξ 1 ) -C k1 (ξ 1 ) -C k2 (ξ 1 )ξ 0 -• • • -C k,k+1 (ξ 1 )ξ k 0 , (24)
where k = 1, . . . , µ 0 . Equality (23) allows to represent constants { Ċµ 0 ,j (0)} as functions of constants {c ij }. Finally, substituting the obtained expressions for constants { Ċij (0)} via constants {c ij } into equation ( 23) we obtain the identity in ξ 0 ∈ Ω 0 for computation of constants {c ij }.

For example, in the case N -= 3 the identity (23) in ξ 0 ∈ Ω 0 for finding of constants c ij takes the form

Ċ31 (0) + Ċ32 (0)ξ 0 + Ċ33 (0)ξ 2 0 + Ċ34 (0)ξ 3 0 = ∂G 3 ∂ξ 1 + + 3 4 (p 2 1 -p 2 ) ∂p 2 ∂ξ 0 -p 1 ∂p 3 ∂ξ 0 - 1 2 p 3 1 -3p 1 p 2 + 2p 3 ∂p 1 ∂ξ 0 ,
where all functions are evaluated at point (ξ 0 , 0), the functions p k are defined in formula (24) and the constants Ċ31 (0), Ċ32 (0), Ċ33 (0), Ċ34 (0) are given by formulas where ae 10 and ae 11 are defined in formula [START_REF] King | Open problems in geomtric function theory[END_REF].

Ċ31 (0) = 1 2
The next theorem permits to find

N -= |X ∩ CP 1 ∞ | from γ = ∂X. Theorem 3.2.
Let X ⊂ CP 2 \ (0 : 1 : 0) be a complex curve with rectifiable boundary γ ⊂ C 2 and satisfying ( * ). Let G m≥1 be the functions defined in (2) and let

N -= |X ∩ CP 1 ∞ |.
Fix any ξ * 0 ∈ Ω 0 and let W ξ * be a neighborhood of ξ * = (ξ * 0 , 0) in C 2 . Then the following statements are valid:

1. If G 1 = 0 in W ξ * , then either N -= 0, or γ bounds a complex curve in C 2 ,
where γ denotes γ with the opposite orientation.

2. If there exist complex constants c 11 , c 12 such that

∂ ∂ξ 1 G 1 -P 1 = (G 1 -P 1 ) ∂ ∂ξ 0 (G 1 -P 1 ) in W ξ * , (25)
where P 1 (ξ 0 , ξ 1 ) = c 11 ξ 0 +c 12 1+c 11 ξ 1 , then N -≤ 1. Furthermore, c 11 , c 12 are the same constants as in Theorem 3.1.

If there exist complex constants a

1 , a 2 , b 1 , b 2 , c 1 , c 2 such that c 1 + c 2 = ae 10 , ∂ ∂ξ 1 ((G 1 -P 1 ) 2 -G 2 -P 2 ) = ((G 1 -P 1 ) 2 -G 2 -P 2 ) ∂ ∂ξ 0 (G 1 -P 1 ) in W ξ * , P 1 (ξ 0 , ξ 1 ) = -a 1 ξ 0 +b 1 1-a 1 ξ 1 -a 2 ξ 0 +b 2 1-a 2 ξ 0 , P 2 (ξ 0 , ξ 1 ) = - 2 j=1 a j ξ 0 +b j 1-a j ξ 1 2 + 2a j c j 1-a j ξ 1 , (26)
then N -≤ 2. Furthermore, these constants are related to the constants of Theorem 3.1 by the equations:

a 1 + a 2 = -c 12 , b 1 + b 2 = -c 11 , a 2 1 + a 2 2 = -c 23 , a 1 b 1 + a 2 b 2 = -1 2 c 22 , b 2 1 + b 2 2 + 2a 1 c 1 + 2a 2 c 2 = -c 21 , c 1 + c 2 = ae 10 . (27)
Proof. 1. The equality G 1 = 0 in W ξ * implies that G 1 = 0 for all ξ ∈ Ω 0 × C. In turn, it implies, according to [START_REF]Chaînes holomorphes de bord donné dans CP n[END_REF], the moment condition

γ z k 1 1 z k 2 2 dz 2 = 0, k 1 , k 2 ≥ 0.
Then, according to [START_REF] Wermer | The hull of a curve in C n[END_REF] and [START_REF] Harwey | Boundaries of complex analytic varieties[END_REF], for an appropriate choice of orientation, γ is the boundary of a complex curve in C 2 .

Set h

= G 1 -P 1 in a neighborhood W ξ * of ξ * . Then h satisfies ∂ 2 ∂ξ 2 0 (G 1 -h) = 0 in W ξ * , ∂h ∂ξ 1 = h ∂h ∂ξ 0 in W ξ * .
It follows from Theorem 2.2 that N -≤ 1. Note also that equation ( 10) is the restriction of (25) to W ξ * ∩ (Ω 0 × 0).

3. Consider the following quadratic equation in variable t:

t 2 -(G 1 -P 1 )t + 1 2 ((G 1 -P 1 ) 2 -G 2 -P 2 ) = 0. ( 28
)
Suppose that the discriminant is non-zero at ξ * . Then, without loss of generality, it is non-zero in W ξ * (we can always choose a smaller neighborhood). We denote two different roots of this equation as

h 1 = h 1 (ξ), h 2 = h 2 (ξ).
Clearly, h 1 and h 2 are holomorphic in W ξ * . Furthermore, by the Vi?te formulas we have

h 1 + h 2 = e 1 := G 1 -P 1 in W ξ * , h 1 h 2 = e 2 := 1 2 ((G 1 -P 1 ) 2 -G 2 -P 2 ) in W ξ * .
Note that by definition

∂P 1 ∂ξ 1 = 1 2 ∂P 2
∂ξ 0 . Note also that by Lemma 3.3.1 of [START_REF]Chaînes holomorphes de bord donné dans CP n[END_REF] we have

∂G 1 ∂ξ 1 = 1 2 ∂G 2 ∂ξ 0 . It leads to the equation ∂ e 1 ∂ξ 1 = 1 2 ∂ ∂ξ 0 ( e 2 1 -2 e 2 ) in W ξ * . (29)
Furthermore, equation (26) can be rewritten in the form

∂ e 2 ∂ξ 1 = e 2 ∂ e 1 ∂ξ 0 in W ξ * . (30)
Now denote by h 1 , h 2 the shock-wave extensions of h 1 (•, 0) and h 1 (•, 0) to W ξ * which exist and are unique by the Cauchy-Kowalevski theorem. Set e 1 = h 1 + h 2 , e 2 = h 1 h 2 . Due to the shock-wave equations for h 1 and h 2 , the functions e 1 and e 2 satisfy ∂ e

1 ∂ξ 1 = 1 2 ∂ ∂ξ 0 ( e 2 1 -2 e 2 ) in W ξ * , ∂ e 2 ∂ξ 1 = e 1 ∂ e 2 ∂ξ 0 in W ξ * .
Thus, e 1 , e 2 and e 1 , e 2 are holomorphic solutions to the same system with the same restrictions at ξ 1 = 0. By the Cauchy-Kowalevski theorem, e 1 = e 1 and e 2 = e 2 . It follows from the Vi?te formulas that h 1 , h 2 coincide with h 1 , h 2 (up to order). Hence, h 1 , h 2 satisfy the shock-wave equations. Applying Theorem 2.2, we obtain that N -≤ 2.

It remains to consider the case when the determinant of equation ( 28) vanishes in W ξ * . Otherwise, it vanishes on (at most) a dimension one analytic set and in any neighborhood of ξ * there are balls where it does not vanish.

The zero discriminant condition reads

(G 1 -P 1 ) 2 = 2(G 2 -P 1 ) in W ξ * .
We set h = 1 2 (G 1 -P 1 ). Then by definition of P 1 and from the discriminant condition we get

∂ ∂ξ 2 0 (G 1 -2 h) = 0 in W ξ * , ∂ h ∂ξ 1 = h ∂ h ∂ξ 0 in W ξ * .
By Theorem 2.2 it implies N -≤ 2. Note also that if N -= 2, then all intersections of X with C 1 ξ , ξ ∈ W ξ * , are double. Finally, note that equation [START_REF] Nieser | Automatic generation of Riemann surface meshes[END_REF] is the restriction of (26) to the set W ξ * ∩ (Ω 0 × 0). Complement 3.2. The statement of Theorem 3.2 can be generalized to the case N -≥ 3 in the spirit of cases N -≤ 2. Such a generalization will be developed in a separate paper together with a statement of Theorem 3.1 for N -≥ 3, indicated in Complement 3.1.

We pass to the description of the algorithm of reconstruction of a complex curve X ⊂ CP 2 satisfying the minimality condition ( * ) from the known boundary γ

= ∂X ⊂ C 2 . Let {ξ k 0 } N k=1 , ξ k 0 ∈ C be an arbitrary grid in C, ξ i 0 = ξ j 0 , i = j, and ξ k 0 / ∈ π 2 γ, k = 1, . . . , N . The complex curve X intersects complex line {z 2 = -ξ k 0 } at points (h s (ξ k 0 , 0), -ξ k 0 ), 1 ≤ s ≤ N + (ξ k 0 , 0).
We are going to present the formulas for finding these points. in contained in Ω 0 . Without loss of generality, one can suppose that |ξ k 0 | < R for all k = 1, . . . , N .

Consider an auxiliary complex curve X R = {(z 1 , z 2 ) ∈ X | |z 2 | R}. Its boundary γ R consists of two disjoint parts (possibly, multiconnected): the first part is γ and the second part γ R is obtained by lifting the circle S R = {z ∈ C | |z| = R} to X via the projection π 2 : X → C.

The complex curve X R does not intersect infinity. Moreover, points of the form a, -ξ k 0 , k = 1, . . . , N belong to X if and only if they belong to X R . Therefore, in order to reconstruct X it is sufficient to reconstruct γ R and then to reconstruct X R , using the algorithm for the case when N -= 0. The algorithm can be formulated as follows:

1. New boundary. Choose a sufficiently large R > 0, so that B c R ⊂ Ω 0 and all ξ k 0 belong to B R . In the case of N -= 1, by virtue of formulas ( 6), we have that

h 1 (ξ 0 , 0) = 1 2πi γ z 1 dz 2 z 2 + ξ 0 -P 1 (ξ 0 , 0), |ξ 0 | = R,
where P 1 can be found using Theorem 3.1. This formula allows to recover γ R and, as a corollary, γ R = γ γ R = ∂X R .

In the case of N -= 2 we have two equalities:

h 1 (ξ 0 , 0) + h 2 (ξ 0 , 0) = 1 2πi γ z 1 dz 2 z 2 + ξ 0 -P 1 (ξ 0 , 0), |ξ 0 | = R, h 2 1 (ξ 0 , 0) + h 2 2 (ξ 0 , 0) = 1 2πi γ z 2 1 dz 2 z 2 + ξ 0 -P 2 (ξ 0 , 0), |ξ 0 | = R,
where P 1 and P 2 can be found using Theorem 3.1. Applying Newton identities and Vi?te formulas, we find h 1 (ξ 0 , 0) and h 2 (ξ 0 , 0). Thus, we have recovered γ R and, as a corollary, γ R = γ γ R = ∂X R .

2. Reduction. In order to find the complex curve X R with boundary ∂X R = γ R we apply the algorithm of reconstruction for the case of N -= 0.

Visualization

To our knowledge, there are, at least, two known algorithms for automatic visualisation of complex curves. The first one was proposed by Trott [START_REF] Trott | Visualization of Riemann surfaces of algebraic functions[END_REF] and requires the knowledge of an analytic expression for the curve. The second algorithm was proposed by Nieser-Poelke-Polthier [START_REF] Nieser | Automatic generation of Riemann surface meshes[END_REF] and requires the knowledge of the branching points and their indices for some fixed projection to C. Our algorithm requires the knowledge of the unordered set of points of intersections of the curve with complex lines C 1 ξ . Let us describe in few words the algorithm of visualisation of complex curves that we use in our examples. Denote by π 1 : C 2 → C the projection into the first factor: π 1 (z 1 , z 2 ) = z 1 . Suppose that X is a complex curve in C 2 such that the covering π 1 : X \ {ramification points} → C has multiplicity L. Consider, for simplicity, a rectangular grid Λ in C:

Λ = z ij 1 : Re z ij 1 = i N , Im z ij 1 = j N , i, j = 0, . . . , N ,
where N is a natural number. Suppose now that we are given the set X Λ = π -1 1 (Λ)∩X and we need to visualize the part of X lying above the rectangle 0

≤ Re z 1 ≤ 1, 0 ≤ Im z 1 ≤ 1.
Let us introduce some terminology. We define a path in Λ as a map γ : {1

, . . . , M } → Λ such that |γ(k + 1) -γ(k)| = 1
N for all admissable k, where M is some natural number.

Let γ : {1, . . . , M } → Λ be a path in Λ and let i : {1, . . . , M } → [1, M ] be the inclusion map. Define the function i

* γ : [1, M ] → C such that i * γ(k) = γ(k) for integer k and i * γ| [k,k+1
] is linear for all admissable k. It is clear that i * γ is a continuous function and hence it can be lifted to X by the map π 1 .

We define a path in X Λ as a map Γ : {1, . . . , M } → X Λ such that γ = π 1 • Γ is a path in Λ and Γ = i * L(i * γ), where i * is the pullback map with respect to i and L(i * γ) is some lift of i * γ to X by π 1 , i. e. L(i * γ) is a continuous map from [1, M ] to X such that π 1 • L(i * γ) = i * γ. We also say that Γ is obtained by lifting of γ.

We will call subsets of Λ and X Λ path-connected if every two points of these sets can be connected by a path in Λ and X Λ , respectively.

Let us describe the practical way to lift paths in Λ to paths in X Λ . Suppose that N is sufficiently large. Let γ : {1, . . . , M } → Λ be a path in Λ and let Γ(1) ∈ π -1 (γ(1)) ∩ X be an arbitrary point. We select Γ(k) ∈ π -1 (γ(k)) ∩ X in such a way that

|Γ(k) -Γ(k -1)| = min |z -Γ(k -1)| : z ∈ π -1 (γ(k)) ∩ X , k = 2, . . . , M.
Then Γ is a path in X Λ obtained by lifting of γ. All possible lifts of γ may be obtained by varying Γ(1). Note that if γ is closed, i. e. γ(1) = γ(M ), Γ need not to be closed.

Finding of ramification points and making branch cuts. The first step in visualization procedure consists in finding of ramification points of X with respect to projection π 1 . Since we have only a finite number of points on X we can find ramification points only approximately. More precisely, we will localize them in small circles. Examples of application of this algorithm are given at Fig. 1 and in Section 5. The visualization algorithm can be easily generalized to the case of general grids. For instance, in our examples we have used a modification with periodic grid.

Examples

1. No points at infinity. We consider the complex curve given by the equation

X 1 = (z 1 , z 2 ) ∈ C 2 : z 2 + 1 = z 2 1 -exp( z 2 4 ) 2 , |z 2 | ≤ 1.9 . (31)
The boundary γ 1 = ∂X 1 of this curve is the disjoint union of two real curves γ r 1 and γ g 1 , see Fig. 2. We are going to recover X 1 from γ 1 . First, note that the real curve |z 2 | = 1.9 in C divides C into two parts Ω 0 = {|z 2 | > 1.9} and Ω 1 = {|z 2 | < 1.9}. We check numerically that G 1 (ξ) is zero in a neighborhood W ξ * of some ξ * = (ξ * 0 , 0), ξ 0 ∈ Ω 0 , and thus N -= 0 by Theorem 3.2. Next, we compute the value

µ 1 = σ 0 (ξ 0 ) = 1 2πi γ 1 dz 2 z 2 + ξ 0 ,
for some ξ 0 ∈ Ω 1 (it will be the same for any ξ 0 ∈ Ω 1 ). It turns out that µ 1 = 4 and thus π 2 : X 1 → Ω 1 is a (ramified) covering with 4 sheets. Denote X 1 ∩ C 1 ξ = {(h j (ξ), -ξ 0 )} 4 j=1 , ξ = (ξ 0 , 0), |ξ 0 | < 1.9. Using Newton identities, we determine the values of symmetric functions σ k : σ 1 = s 1 , σ 2 = 1 2 (s 2 1 -s 2 ), σ 3 = 1 6 (s 3 1 -3s 1 s 2 + 2s 3 ), σ 4 = 1 24 (s 4 1 -6s 2 1 s 2 + 3s 2 2 + 8s 1 s 3 -6s 4 ). Finally, taking into account Vi?te formulas, we find h 1 (ξ 0 , 0), . . . , h 4 (ξ 0 , 0) as the roots of the following polynomial in t: t N + (ξ 0 ) -σ 1 (ξ 0 )t N + (ξ 0 )-1 + • • • + (-1) N + (ξ 0 ) σ N + (ξ 0 ,0) (ξ 0 ) = 0.

Reconstructed surface is represented at Fig. 2.

2. One point at infinity. Consider the the complex curve X 2 which is the part of the curve X 2 = (w 0 : w 1 : w 2 ) ∈ CP 2 : w 2 1 = w 2 2 + w 2 0 , bounded by the real curve γ 2 = {(z 1 (t), z 2 (t)) ∈ C 2 : t ∈ [0, 2π]}, z 2 (t) = 1.9e it , z 1 (t) = z 2 2 (t) + 1, z 1 (0) > 0, Next, taking into account the Vi?te formulas, we find the values h 1 (ξ 0 , 0), h 2 (ξ 0 , 0), |ξ 0 | = R, as the roots of the following polynomial in t: t 2 -σ 1 (ξ 0 )t + σ 2 (ξ 0 ) = 0.

Then we put γ 3,R = γ 3 γ 3,R and recover X 3,R from γ 3,R using the algorithm for the case of no points at infinity. The result of reconstruction is depicted at Fig. 4.

  where C 11 , C 12 are holomorphic in a neighborhood of zero. We need to find constants c 11 = C 11 (0) and c 12 = C 12 (0).

-2ae 10

 10 Ċ12 (0) + 2c 11 Ċ11 (0) + Ċ21 (0) = -c 2 12 + c 23 ae 10 -c 2 11 -2c 12 ae 10 + c 21 c 12 .

  0) = c 23 , Ċ23 (0) = -c 3 12 -3c 12 c 23 , Ċ22 (0) = -2 c 11 c 2 12 + c 12 c 22 + c 11 c 23 , Ċ21 (0) = ae 10 (c 2 12 + c 23 ) -c 12 (c 2 11 + c 21 ) -c 11 c 22 .

Figure 1 :

 1 Figure 1: Rieman surfaces of functions f (z) = sin(z), |z| 2 (left) and f (z) = √ z 4 + 1, |z| 2 (right) obtained by the visualization algorithm.

Figure 2 :z 2

 22 Figure 2: The given boundary γ 1 = γ r 1 γ g 1 (red and green) and the reconstructed complex curve X 1 ⊂ C 2 . The color indicates Im(z 1 ).

Figure 4 :2πi γ z 1 dz 2 z 2 z 2

 422 Figure 4: The real curve γ 3 (red) and the reconstructed complex curve X 3 with two points at infinity. The color indicates Im(z 1 ).

  ae 10 3c 11 c 2 12 + 3c 12 c 22 + 3c 11 c 23 + 2c 33

	-	1 2	ae 11 c 3 12 + 3c 12 c 23 + 2c 34 -

The algorithm takes as input the points {ξ k 0 } N k=1 and the curve γ (for example, represented as a finite number of points belonging to γ). On the output of the algorithm we obtain the set of points (h s (ξ k 0 , 0), -ξ k 0 ), 1 ≤ k ≤ N ; 1 ≤ s ≤ N + (ξ k 0 , 0).

No points at infinity

1. Computation of N + . According to formula (4), for every domain Ω l≥0 , the number µ l = N + (ξ 0 , 0), ξ 0 ∈ Ω l , is equal to the winding number of the curve π 2 γ with respect to a point ξ 0 ∈ Ω l :

2. Computation of power sums. If N -= 0 then, according to Theorem 3.1, for every point ξ k 0 ∈ Ω l≥0 we have that P m (ξ k 0 , 0) = 0. Using formula ( 6) we obtain the following formulas for the power sums of the functions to be determined:

3. Computation of symmetric functions. For every point ξ k 0 ∈ Ω ≥1 , the Newton identities

allow to reconstruct the elementary symmetric functions:

Desymmetrisation.

For every point ξ k 0 ∈ Ω l , using Vi?te formulas, one can find the complex numbers h 1 (ξ k 0 , 0), . . . , h µ l (ξ k 0 , 0) (up to order). The points

are the required points of the complex curve X.

One or two points at infinity

These cases can be reduced to the case N -= 0 in the following way. Since π 2 γ ⊂ C is a compact real curve, there exists such R > 0, that the set

Without restriction of generality we suppose that all ramification points are projected by π 1 into interior points of Λ. Take any interior point z 1 ∈ Λ and select a small closed path γ : {1, . . . , M } → Λ around z 1 so that there is at most one ramification point inside the polygon γ(1) . . . γ(M ). For example, one can take as γ the following path:

where i is the imaginary unit. Now consider different lifts of γ to X Λ . If at least one lift is not closed, mark z 1 as a possible ramification point (meaning that it is situated near the projection of some ramification point of X). Now vary z 1 and mark all possible ramification points. The resulting set will consist of several path-connected components each of which localizes the position of one ramification point of X with respect to π 1 .

Next, we need to cut the grid Λ with branching points removed into simply connected domains. The natural choice is the shortest cut graph, which can be computed by the algorithm proposed in [START_REF] Kälberer | Stripe parametrization of tubular surfaces[END_REF], which, in turn, generalizes the algorithm for closed surfaces of [START_REF] Erickson | Greedy optimal homotopy and homology generators[END_REF].

Denote the cut graph by Λ c . An important observation is that every closed path in Λ \ Λ c always lifts to a closed path in X Λ since it doesn't contain π 1 -projections of ramification points inside.

Visualization. Now denote Λ \ Λ c = ∪ S s=1 Λ s , where Λ s are different path-connected components. Take any z s 1 ∈ Λ s and z s 2 ∈ π -1 1 (z s 1 ) ∩ X. Now take other z 1 ∈ Λ s and connect z s 1 with z 1 by some path γ. Then γ lifts to a path Γ with Γ(1) = (z s 1 , z s 2 ) and Γ(2) = (z 1 , z 2 ) for some z 2 ∈ π -1 1 (z 1 ) ∩ X and z 2 doesn't depend on γ. Varying z 1 we thus obtain the map Σ(z s 1 , z s 2 ) : Λ s → X Λ which allows us to visualize the part of X. Varying z s 2 ∈ π -1 1 (z s 1 ) ∩ X (the latter is the finite set, namely, it consists of L elements) we obtain the other maps Σ(z s 1 , z s 2 ) which allow us to visualize other parts of X. Clearly, the set of obtained maps doesn't depend on the choice of z s 1 ∈ Λ s . Hence we can denote the obtained maps by Σ l s , l = 1, . . . , L. It is clear that

The part π -1 (Λ c ) ∩ X Λ consists of cuts and preimages of possible ramification points.

The cuts can be visualized as the already visualized part of the surface. The only problem is the visualization of π 1 -preimages of possible ramification points. But the latter take a little part of the surface when N is large and one can just forget about their visualization. On the other hand, in our examples they were visualized using low-level graphics approach.

intersecting infinity at (0 : 1 : 1). We suppose that the number of points at infinity is apriori unknown. We are going to reconstruct X 2,R = X 2 ∩ {z 2 < R}, where R = 3, from γ 2 . Note that the curve

According to the algorithm, the first step is to determine the number of points at infinity. We choose some ξ * = (ξ * 0 , 0) with ξ * 0 ∈ Ω 0 and check that in some small neighborhood W ξ * of ξ * we have G 1 ≡ 0 and thus N -≥ 1 by Theorem 3.2.

Next, we fix two generic ξ 1 0 , ξ 2 0 ∈ Ω 0 and determine c 11 , c 12 from the linear system c 11

This system is uniquely solvable and gives c 11 ≈ 1, c 12 ≈ 0.

Next, we check numerically that in the neighborhood W ξ * of ξ * the following identity holds:

Theorem 3.2 together with the estimate N -≥ 1 imply that N -= 1.

Next, according to the algorithm, in order to reduce the problem to the case of no points at infinity, we need to recover the auxilary curve γ

∞ . Thus, it remains to recover X 2,R from the known boundary γ 2,R using the algorithm for the case of no points at infinity, see the previous example. The reconstructed curve is depicted at Fig. 3.

3. Two points at infinity. Consider the complex curve X 3 , which is the part of the algebraic curve

, bounded by the real curve

The curve X 3 has two points at infinity, namely, q 1 = (0 : 1 : e 2πi 3 ), q 2 = (0 : 1 : e -2πi

3 ). We suppose that the number of points at infinity is a priori unknown. We are going to reconstruct X 3 from γ 3 . The first step is to determine the number of points at infinity. Choose a generic ξ * = (ξ * 0 , 0), ξ * 0 ∈ Ω 0 , and some neighborhood W ξ * of ξ * .

1. We check that G 1 ≡ 0 in W ξ * so that N -≥ 1 by Theorem 3.2.

2. We choose two generic ξ 1 0 , ξ 2 0 ∈ Ω 0 and solve (32) for c 11 , c 12 . It turns out that system (32) is uniquely solvable. However, identity (33) does not hold in W ξ * . Using Theorem 3.2, we obtain that N -≥ 2.

3. We consider [START_REF] Nieser | Automatic generation of Riemann surface meshes[END_REF] as an identity in a neighborhood of ξ * 0 . We find numerically that c 11 ≈ 0.9974, c 12 ≈ -0.3335, c 21 ≈ 1.0133, c 22 ≈ -0.6699, c 23 ≈ 0.5567.

We check that (26) holds in a neighborhood of W ξ * with constants a 1 , a 2 , b 1 , b 2 , c 1 , c 2 given by equations ( 27). If follows from Theorem 3.2 and from the inequality N -≥ 2 that N -= 2.

Suppose that we want to recover X 3,R = X 3 ∩ {|z 2 | < R}, R = 3. According to the algorithm, we need to find the auxilary boundary γ 3,R = X 3 ∩ {|z 2 | = R}. Denote γ 3,R = {(h j (ξ 0 , 0), -ξ 0 ) : |ξ 0 | = R, j = 1, 2}. We compute the values of symmetric functions σ 1 (ξ 0 ) = h 1 (ξ 0 , 0) + h 2 (ξ 0 , 0), σ 2 (ξ 0 ) = h 1 (ξ 0 , 0)h 2 (ξ 0 , 0) using the explicit


Sami Evangelista
email: evangeli@cnam.fr

High Level Petri Nets Analysis with Helena

This paper presents the high level Petri nets analyzer Helena. Helena can be used for the on-the-fly verification of state properties, i.e., properties that must hold in all the reachable states of the system, and deadlock freeness. Some features of Helena make it particularly efficient in terms of memory management. Structural abstractions techniques, mainly transitions agglomerations, are used to tackle the state explosion problem. Benchmarks are presented which compare our tool to Maria.

Introduction

Model checking is an automatic method for the verification of finite state systems. It consists of enumerating all the possible configurations or executions of the system to track the ones which do not match the specification.

Helena (a High LEvel Net Analyzer) is a model checker developed at the CNAM university in Paris. The formalism supported by Helena makes it suitable for the verification of concurrent software. Helena is part of the Quasar project (a tool for the analysis of concurrent Ada programs [START_REF] Evangelista | Quasar : a new tool for analyzing concurrent programs[END_REF]), but it is also a fully autonomous tool which can be used independently. In its current version, Helena can be used for the verification of state properties, i.e., properties that must hold in all the reachable states of the system, and deadlock freeness. Helena is a command-line oriented tool without any graphical user interface though we consider to include a graphical interactive simulator.

This paper is structured as follows. Section 2 explains the reasons which motivated us to design and implement Helena. The main features of Helena are presented in Section 3. A set of benchmarks which compare our tool to Maria [START_REF] Mäkelä | modular reachability analyser for algebraic system nets[END_REF] are presented in Section 4. Finally we conclude in Section 5.

Motivations

The model checking of concurrent programs involves a translation task from the original programming language to a formalism suitable for the expression of concurrency, e.g., Promela, Petri nets. In order to limit the state explosion problem of a model checking procedure, the produced model has to remain as small as possible, but still must be an exact translation of the program to guarantee a correct verification process. Quasar is a tool which aims at the automatic verification of concurrent Ada programs based on colored Petri nets. It can not currently handle the whole language, though a large part of the language which is related to concurrency is supported. The main task of Quasar is to perform automatic abstraction (or slicing) of a program and to translate it into a colored Petri net. In its current state, Quasar does not perform reachability analysis but interfaces with another tool which realizes this task. The high level formalism supported by Maria [START_REF] Mäkelä | modular reachability analyser for algebraic system nets[END_REF] pushed us into choosing this tool, though Quasar can also interface with Prod [START_REF] Varpaaniemi | PROD 3.4.00 -an advanced tool for efficient reachability analysis[END_REF].

After an intensive use of Quasar on various examples, we identified two main factors which limit the efficiency of our tool. Firstly, the translation from the source code to the colored Petri net is not straightforward. As the constructions of high level programming language do not have their exact counterpart in colored Petri nets, even in the Maria formalism, the translation step introduces additional transitions which generate intermediate states that are not relevant for the verification purpose. Secondly, the state vectors exhibited by models obtained from the translation of a program are usually quite large as software make heavy use of structured data types. Thus, though optimized to represent multiplicity in a compact way [START_REF] Mäkelä | Condensed storage of multi-set sequences[END_REF], the encoding scheme of Maria fails to represent efficiently state vectors with large color domains, and many tokens in places.

Helena has been designed from this previous experience to meet these two requirements : enable a straightforward and automatic translation of concurrent programs without resort to the introduction of many useless transitions and intermediate states, and handle state spaces with large state vectors.

Main Features of Helena

High Level Description Language. The class of high level Petri nets used in Helena was primarily designed to enable the simulation and verification of concurrent Ada programs. To achieve this, we naturally decided to include in Helena the possibility to define high level data types. There are four categories of data types : numerical type, enumerate type, structured type, and vector, i.e., array type. Another feature of Helena is the possibility for the user to define complex functions written in a pseudo-C syntax. These functions may then be used in arc expressions. Such a possibility allows to automatically translate sub programs and sequences of statements which do not include synchronizations into a single transition, provided that the input programming language can easily be mapped into Helena functions.

Compilation and Execution of the Model.

A well known and efficient approach to reduce the execution time of a model checking procedure is to compile the model into a source code which will correspond to the actual reachability analyzer. Compiling and executing the model has numerous advantages over an interpretation of the model. Mainly, the evaluation of the expressions in arc mappings can be drastically fasten. Tools that use this technique include Prod, Spin [START_REF] Holzmann | The model checker spin[END_REF], and Maria. It has been shown in [START_REF] Mäkelä | Applying compiler techniques to reachability analysis of high-level models[END_REF] that this technique greatly reduces the execution time even for small models for which we may think that the compilation of the generated code is a too severe overhead (which does not exist if the model is interpreted).

Helena models are also translated into executable code. For performance and portability issues we chose the language C. In order to ease the readability of the generated code, and the debugging of the compilation process, this code is commented, nicely formatted, and divided into several libraries.

Enabling Test Algorithm. Verification and simulation tools based on high level Petri nets spend a significant amount of time in determining under which assignments (or bindings) a high level transition is firable. This non trivial problem is known under the term of enabling test. The enabling test algorithm implemented in Helena has been described in [START_REF] Evangelista | An efficient algorithm for the enabling test of colored petri nets[END_REF]. It basically includes two main components.

1. We exploit the locality principle of Petri nets which states that the firing of a transition only affects the status (enabled/disabled) of its neighbor transitions. Thus, a depth first search algorithm can benefit from this locality principle by (a) maintaining a set of enabled transitions (b) updating this set when a transition instance is fired by only inspecting the neighbor transitions of the fired instance. Our implementation of this locality principle is based on the definitions of structural conflict and causality [START_REF] Dutheillet | Conflict sets in colored petri nets[END_REF]. The basic idea is to translate these relations into equivalent constraints systems [START_REF] Evangelista | Syntactical rules for colored petri nets manipulation[END_REF] before the search algorithm. During the search, the systems built are solved in order to identify disabled and enabled bindings. We illustrate our purpose with an example. Let us consider the net of figure 1. From the structure of the net and the arc expressions we can deduce that the firing of an instance (t, X t ) disables the firing of an instance (u, X y , Y u ) if the following constraint is verified : [X t = X u ∧ Y u = 0]. During the search, at each firing of an instance of transition t we instantiate this system with the firing binding to identify the instances of u which are disabled. We also observe that the firing of an instance (t, X t ) enables the firing of an

q p u t v X t X t , 0 X v X u , Y u Fig. 1. Translating arc expressions into constraints system instance (v, X v ) if the following holds : [X v = X t ].
Consequently, the firing of an instance (t, X t ) causes the insertion of the instance (v, X t ) into the set of enabled transitions without any additional check. 2. We showed in [START_REF] Evangelista | An efficient algorithm for the enabling test of colored petri nets[END_REF] that this approach is unfortunately not always sufficient to determine valid transitions assignments at the new marking. In most cases, this must be followed by an unification algorithm. The algorithm we proposed is an improvement of Mäkelä's unification algorithm [START_REF] Mäkelä | Optimising enabling tests and unfoldings of algebraic system nets[END_REF].

State Space Storage. Model checking tools usually represent states in two different forms :

an expanded form which is convenient for the implementation of the transition relation a compressed form in which expanded states are encoded before their insertion into the state space in order to save memory An interesting feature of Helena is the way it represents compressed states. The method, called ∆-markings method, has been described in [START_REF] Evangelista | Efficient state space storage in explicit model checking[END_REF]. The underlying idea is to store a large set of states in a non explicit way by only storing references on others states. Figure 2 illustrates the principle of the method. Markings met at a depth d such that d mod k δ = 0 (with k δ a user defined parameter) are stored explicitely whereas other markings are represented in the state space by a couple (pred, (t, c t )), where pred is the "address" in the hash table of one of the predecessors of the marking and (t, c t ) is the transition instance which firing leads from the predecessor to the marking. Retrieving the actual value of a marking from its "∆ encoding" can simply be done by following the links which point on the predecessors until a marking stored explicitely is reached. The marking is then obtained by firing the sequence of instances (t, c t ) which label the links followed to reach the explicit marking. For instance, the actual value of marking m can be retrieved by first backtracking to m and then to m 0 and to apply on it the firing of sequence (t, c t ).(t , c t ). For models exhibiting large state vectors, the compression ratios observed are quite impressive, whereas the method becomes less interesting for small models. The price to pay is an acceptable increase of the run time. Both the compression ratio and the run time increase can be influenced by parameter k δ . The state collapsing [START_REF] Visser | Memory efficient state storage in spin[END_REF] method of Spin is also implemented in Helena. This method is based on the observation that even if the set of syntactical possible values for a token in a place is huge (let us denote it size by m), the set of values really met during the search (which size is denoted by n) is in practice usually much smaller. This is so because the types of the places and transitions of the net are usually over-approximations made by the user of the possible values really met during the search.

With the collapse method a token is represented with log 2 (n) bits. This "collapsed token" is in fact an index of a table of size n which stores all the "true tokens" already met during the search on log 2 (m) bits. This table is initially empty, and filled during the search by the token values met. Thus, we make use of the two following functions to query the tokens table: index of (token t) → int looks for token t into the tokens table and returns the index at which the token has been found. It may also add an entry in the table if the token is not already in it. -token at(int i) → token returns the token at index i in the tokens table.

Since n can not be known a priori, Helena can detect overflows in this table, i.e., cases where the number of different tokens met is greater than n, and report it to the user who, in turn, can increase it and rerun the search. This restriction may seem quite bothering, but after intensive experiments, we observed that overflows are usually quickly detected during the search, causing only a small run time overhead.

Structural Abstraction Techniques.

The efficiency of the explicit model checking approach is seriously limited by the state explosion problem inherent to the concurrent execution of several components. An efficient way to tackle this problem is to perform structural abstractions on the initial model in order to obtain a simpler model which is equivalent to the initial one for to the specified property. For Petri nets and Colored Petri nets, transitions agglomerations are surely the most efficient abstraction techniques. They have been defined by Berthelot in his doctoral thesis [START_REF] Berthelot | Transformations and decompositions of nets[END_REF], and generalized by Haddad to Colored Petri nets in [START_REF] Haddad | A reduction theory for colored nets[END_REF]. Transition agglomerations merge consecutive transitions into a virtual atomic one which effect is the composition of the effects of these transitions. It results in a drastic reduction of the combinatory explosion due to the elimination of some intermediate states. In addition, the complexity of these transformations is linear with respect to the size of the net, and their application is totally automatic.

These agglomerations are implemented in Helena. This implementation is based on the symbolic calculus of some structural relations which ensure correct agglomerations in the unfolded net. Though this symbolic computation is not always possible, it usually works fine. To our best knowledge, Helena is currently the only tool that implement structural agglomerations for high level Petri nets.

Stack Representation. Model checking algorithms based on a depth first search, e.g., LTL model checking, store the set of non fully processed states in a stack. For some models, this stack can grow very large, and include almost all the states of the state space. For instance, models for which the state space forms a single strongly connected component usually follow this behavior. Thus, an awkward choice for the representation of this stack can lead to important memory wastes.

In Helena, the stack is represented as a vector of bits, and it uses a nice property of Petri nets : the transition relation is a reversible mechanism, that is, given a marking m and a transition t, there is a single marking m such that m [t m. This marking m is defined by

∀p ∈ P, m (p) = m(p) + W -(p, t) -W + (p, t).
Instead of representing the search stack as a stack of states we therefore chose to represent it as the sequence of transition bindings which leads from the initial marking to the current marking. When a state has been fully expanded, the binding on top of the stack is "unfired", and the next enabled binding is processed. This representation can naturally be combined with the collapse method to store transition bindings of the stack more compactly.

Expressing State Properties. The query language of Helena is rich enough to express a wide range of state properties. It is based on four basic operations :

-The count operation allows to count the number of distinct items in a place which fulfill a given condition. -The mult operation allows to count the cumulated multiplicities of distinct items in a place which fulfill a given condition. -The exists token operation allows to check that at least one token in a place fulfills a given condition. -The forall token operation allows to check that all the tokens present in a place fulfills a given condition.

For instance, the property forall_token(p in P : exists_token(q in Q : q->1 = p->2))

could be stated in an informal way as : for each token p held in place P there is a token q in place Q which is such that the first component of q and the second component of p are equal.

Benchmarks

We have compared our tool with Maria to study the performances of both tools concerning time and memory consumption. We considered these six examples :

-The distributed database system -The slotted ring protocol -The dining philosophers -The sieves of Eratosthene to find prime numbers -The leader election protocol of Chang and Roberts -The Peterson algorithm for the mutual exclusion problem

All these examples can be found in the Helena distribution (directory samples).

Experiments were done on a Pentium 4 with 2.5Ghz and a main memory of 1Gb. To allow a fair comparison, Maria was invoked with option --compile which allows to compile arc expressions into a dynamically linkable library. In addition, no capacity constraint was indicated in the Maria model and structural agglomerations were disabled in Helena in order to obtain the same numbers of markings and arcs. The results of our experiments are reported in table 1. For each model, row T reports the times observed for the exploration of the state space in seconds, row S reports the size of the state space in kilobytes, and row V reports the average size of the state vector in bytes. Compilation times of the models are included. Helena was invoked with several values of parameter k δ within the set {1, 5, 10, 20, 50}. Let us recall that k δ is the parameter of the compression method implemented in Helena. The higher this parameter is, the more compressed the state space will be, and the more slow the search will be. If k δ = 1 then no compression is performed. In addition, the state collapsing method was not helpful for these simple models and therefore disabled. We observe that the results obtained by Helena without any compression technique are comparable to the ones obtained by Maria. The encoding scheme of Maria gives a slight advantage to Maria in terms of memory usage, while Helena usually explores the state space in less time. When using the ∆-markings method, Helena outperforms its competitor concerning the size of the state space, especially for models with large state vectors (the distributed database system, and the sieves of Eratosthene). The best compression ratio is obtained for the distributed database system, with k δ = 50. For this case, the size of the state space is divided by almost 19. The drawback is an increase of the execution time. This one is acceptable for low values of k δ but tends to grow with it. For models with small vectors, e.g., the slotted ring protocol, our storage method is clearly less interesting though a reduction factor of 2, which is the average reduction factor observed for these models, can still be helpful for systems with large numbers of states. Finally, we observe that with our storage method, the average size of the state vector could be reduced to approximatively 10 bytes, whatever the model is.

Conclusions and Perspectives

The validation of software is a difficult problem. Few tools based on high level Petri nets have been designed to face this challenge as most of them focus on the representation of control, e.g. tools based on Well Formed Petri nets such as GreatSPN [START_REF] Chiola | Greatspn 1.7: graphical editor and analyzer for timed and stochastic petri nets[END_REF], over data. Thanks to multiple features, Helena can handle the validation of software systems.

The possibility to define high level data types and functions enables to model concurrent software written in high level programming languages such as Ada in a succinct way.

The ∆-markings method [START_REF] Evangelista | Efficient state space storage in explicit model checking[END_REF] is implemented in Helena. This one is particularly efficient when dealing with state spaces with large state vectors. This is why we believe that it should be adapted within the scope of the verification of software as these make heavy use of structured data types and usually exhibit large state vectors.

Lastly, structural abstractions help to tackle the state explosion problem. We plan to extend Helena in the following ways :

-At the current implementation stage, Helena supports the verification of state properties and the deadlock freeness. We plan to extend the field of properties that Helena can verify by including a module for the verification of linear time temporal logic properties (LTL). An interface with the extensible library SPOT [START_REF] Duret | Spot: an extensible model checking library using transition-based generalized Büchi automata[END_REF] is currently under study. The main interest of using SPOT is that it relies on transition-based generalized Büchi automata and allows translation of LTL formula to smaller automata and thus, smaller synchronized products.

-The current version of Helena tackles the state explosion problem by the use of structural abstraction techniques. We envisage to combine these with the stubborn set method of Valmari [START_REF] Valmari | On-the-fly verification with stubborn sets[END_REF][START_REF] Valmari | The state explosion problem[END_REF][START_REF] Michael | Finding stubborn sets of coloured petri nets without unfolding[END_REF]. However, defining an algorithm which computes good stubborn sets for colored Petri nets is a difficult task.

To our best knowledge, the only tool that compute stubborn sets for colored Petri nets (without unfolding) is CPN tools [START_REF] Michel Beaudouin-Lafon | CPN/tools: A tool for editing and simulating coloured petri nets[END_REF]. The algorithm has been described in [START_REF] Michael | Finding stubborn sets of coloured petri nets without unfolding[END_REF]. -New transitions agglomeration rules have been recently defined for ordinary and colored Petri nets in [START_REF] Haddad | New powerfull Petri nets reductions[END_REF][START_REF] Evangelista | New coloured reductions for software validation[END_REF]. An implementation of these agglomerations in Helena is scheduled.

G.

  Ciardo and P. Darondeau (Eds.): ICATPN 2005, LNCS 3536, pp. 455-464, 2005. c Springer-Verlag Berlin Heidelberg 2005

Fig. 2 .

 2 Fig. 2. Illustration of the ∆-markings method

Table 1 .

 1 Results of the Comparison of Helena with Maria

		Maria			Helena
			k δ = 1	k δ = 5	k δ = 10	k δ = 20	k δ = 50
		The distributed database system, N=12, 2 125 765 states
	T	900	452	492	533	778	782
	S	303 409	336 330	80 111	52 496	16 493	15 997
	V	146.15	162.01	38.59	25.29	7.95	7.71
		The dining philosophers, N=12, 4 126 351 states
	T	360	306	425	474	549	647
	S	57 992	67 947	35 679	31 587	29 527	28 290
	V	14.39	16.86	8.85	7.84	7.33	7.02
		The sieves of Eratosthene, N=40, 2 028 969 states
	T	116	90	121	151	206	368
	S	82 056	92 868	35 274	28 046	24 432	22 294
	V	41.41	46.87	17.80	14.15	12.33	11.25
		The leader election protocol, N=14, 1 518 111 states
	T	155	103	141	181	264	354
	S	29 367	33 258	15 930	13 667	12 786	11 616
	V	19.81	22.43	10.75	9.22	8.62	7.84
		The Peterson algorithm, N=4, 3 407 946 states
	T	136	91	133	154	188	315
	S	35 899	43 366	28 287	26 511	25 630	25 094
	V	10.78	13.03	8.50	7.97	7.70	7.54
		The slotted ring protocol, N=8, 3 294 720 states
	T	214	164	245	276	325	401
	S	37 622	42 760	26 370	23 456	23 342	22 730
	V	11.39	13.29	8.20	7.57	7.25	7.06

Availability. Helena is a free software available under the conditions of the GNU General Public License. It can be downloaded at the following URL : http://helena.cnam.fr.


Clara Lévy

Scarlett Gendrey

Séverine Bernardie

Marie-Aurélie Chanut

Aurélien Vallet

Laurent Dubois

Jean-Paul Duranthon

Prediction of displacement rates at an active landslide using joint inversion of multiple time series

Keywords: Landslide, displacement pattern, prediction, monitoring, joint inversion

This work focuses on the development of FLAME (Forecasting Landslides induced by Acceleration Meteorological Events) that analyze of the relationship between displacements and precipitations using a statistical approach in order to predict the surface displacement at active landslide. FLAME is an Impulse Response model (IR) that simulates the changes in landslide velocity by computing a transfer function between the input signal (e.g. rainfall or recharge) and the output signal (e.g. displacement). This model has been applied to forecast the displacement rates at Séchilienne (French Alps). The FLAME model is enhanced by achieving the calibration using joint inversion of multiple time series data. We consider that the displacements at two different sensors are explained by the same long-term response of the system to ground water level variations. The parameters describing the long-term response of the system are therefore identical for all sensors. The joint inversion process allows decreasing the ratio between the number of parameters to be inverted and the volume of data and is thus more statically steady. The results indicate that the models are able to reproduce the displacement pattern in general to moderate kinetic regime but not extreme kinetic regime. Our results do not give clear evidence of an improvement of the models performance with joint inversion of multiple time series of data. The reasons which could explain these inconclusive results are discussed in the paper.

Introduction

Recently, statistical models using impulse responses functions had been used to successfully predict landslides velocities [START_REF] Belle | The application of an innovative inverse model for understanding and predicting landslide movements (salazie cirque landslides, reunion island)[END_REF]Abellan et al. 2013). Not only these models facilitate displacement prediction based on a limited data set [START_REF] Belle | The application of an innovative inverse model for understanding and predicting landslide movements (salazie cirque landslides, reunion island)[END_REF]), but they also proved appropriate to take into account possible temporal variations of the landslide properties [START_REF] Bernardie | Prediction of changes in landslide rates induced by rainfall[END_REF]Abellan et al. 2015). Such models predict a behavior of the landslide similar to the period used to calibrate the model. Consequently, when the landslide deformation rate increases and thereby causes a rapid change of the landslide rheological and hydrological properties, the discrepancy between predicted and actual measurements increases. [START_REF] Bernardie | Prediction of changes in landslide rates induced by rainfall[END_REF] showed that this indication of a change in the landslide behavior could be used to predict fluidization events at the Super-Sauze landslide, several days in advance. At Super-sauze, a threshold criterion of the error of prediction was estimated by performing back analysis on several fluidization events.

The definition of thresholds values for real-time early warning systems becomes tricky for study-sites such as Séchilienne that never experienced a major crisis. One way to identify periods when the observed velocities significantly differ from the predicted velocities would be to assess the quality of the velocity predictions and to estimate an error envelope of these predictions. Thus, prior to apply such methods, we must ensure the quality of model predictions. Unfortunately, the studies aforementioned also showed that this latter type of statistical models usually presents some defaults. In particular, the inversion process can be rather instable when choosing to calibrate a model accounting for a complex hydrogeological response (i.e., with numerous parameters to be determined by inversion) using short time series of data.

Therefore, the improvement of the model predictions is an important step to help defining critical values of velocities. This study discusses the enhancement of the method by achieving the calibration of the model using joint inversion of multiple time series of data. The method is applied at the Séchilienne landslide where the geometry of the instrumentation meets the requirement of the method. In detail, it means that the displacement rate sensors are located on top of the same hydrogeological unit, thus allowing considering partially identical Impulse Response models for all the inverted time series. By doing so, the ratio between the number of parameters to be inverted and the volume of data decreases, thus a priori improving the robustness of the method. The first part of this paper details the current understanding of the hydrogeological behavior of the Séchilienne landslide. The second part of this paper presents the statistical model used in this study to predict surface velocities at the Séchilienne landslide. Finally we discuss the benefits of using joint inversion of multiple time series of data for the prediction of the displacement rates at active landslides using Impulse Response models.

Description of the Séchilienne landslide

Location, geometry

The Séchilienne landslide is located in the South-east of Grenoble in the French Alps (Fig1). The whole unstable mass is estimated between 48 and 63 million m 3 [START_REF] Le Roux | Deep geophysical investigation of the large séchilienne landslide (western alps, france) and calibration with geological data[END_REF]. Most of the unstable mass moves slowly with velocities ranging from 2 to 15 cm/year on average, except for a very active moving zone with velocities ranging from 150 to 300 cm/year on average (Vallet et al. 2015b). Rock falls occurs on a regular basis within this active zone which volume is estimated at 3 million m 3 [START_REF] Le Roux | Deep geophysical investigation of the large séchilienne landslide (western alps, france) and calibration with geological data[END_REF].

A particularity of the Séchilienne landslide is the absence of a well-define basal sliding surface.

Kinetics and hydrogeology

The influence of recharge on the kinematics of the Séchilienne landslide has been evidenced by numerous studies: [START_REF] Rochet | La déformation du versant sud du mont-sec dans le secteur des ruines de séchilienne (isère)[END_REF][START_REF] Alfonsi | Relation entre les paramètres hydrologiques et la vitesse dans les glissements de terrains. exemples de la clapière et de séchilienne (france)[END_REF][START_REF] Durville | The séchilienne landslide: monitoring and kinematics[END_REF]Chanut et al. 2013. Recently, Vallet et al. (2015b) showed that the massif at Séchilienne hosts a relatively shallow perched aquifer in the unstable zone and a deep aquifer in the intact massif. The existence of the shallow perched aquifer is related to the permeability contrast between the unstable mass at shallow depth and at deeper depth. Authors described the perched aquifer, with an extent and connectivity dependent on the seasonal recharge. For both aquifers, conductive fractures play a major role in the drainage (Vallet et al. 2015b). [START_REF] Cappa | Deep fluids can facilitate rupture of slow-moving giant landslides as a result of stress transfer and frictional weakening[END_REF], and more recently Vallet et al. (2015b) investigated the link between the Séchilienne landslide kinematics and the local hydrogeology. Both showed that low frequency displacement of the landslide might be linked with variations of the deep aquifer level, and whereas the high frequency displacements of the landslide might by linked with variations of the perched aquifer level in the landslide.

The existence of this deep aquifer has never been evidence by drilling experiments (since its depth is lower than -150 m). Its depth and geometry are not well known and this deep aquifer might correspond to water circulation inside a fracture network.

Monitoring network

The Séchilienne landslide is monitored since 1985 by the French public national body Cerema Centre Est [START_REF] Duranthon | Apport des méthodes topographiques et topométriques au suivi du versant rocheux instable des ruines de séchilienne[END_REF]Dubois et al. 2014). Precipitations are recorded at the Mont-Sec weather station, located a few hundred meters above the disturbed zone. This station is equipped with both a rain and a snow gauge, thereby allowing the estimation of snow cover and snow melt in water equivalent.

The landslide displacements are monitored using a variety of techniques (extensometers, radar, infra-red, inclinometers and GPS). In this study, we analyze displacements at extensometers A13 an A16, both located within the active moving zone (see Figure 1b). The measurements show the same general trends (see Fig. 5a). The extensometers are characterized by measurement errors of ± 0.5 mm.

Description of the forecasting model FLAME

We use the statistical model FLAME [START_REF] Bernardie | Prediction of changes in landslide rates induced by rainfall[END_REF] to predict surface velocities at the Séchilienne landslide.

FLAME is an Impulse Response model, which principle is to reproduce an output signal S (e.g. displacement) using the convolution product of an input signal E (e.g. rainfall or recharge) by a transfer function Γ as described in Eq.1 and Fig. 2:

( . ) = * ( ) = ∑ ( . ). ( -+ 1). [1]
where n is the discretized interval time, and k is the order (length) of the impulse response.

Fig. 2 Schematic framework of the statistical model FLAME

During the calibration phase, an optimization of the parameters describing the transfer function Γ is performed in order to minimize the Root Mean Square Error (RMSE) of the difference between the estimated displacement rates and the observed displacement rates. The calibration model is then used to forecast the displacement rates coming just after the calibration period (Fig. 3). The calibration of the models is repeated daily, using a moving calibration period of 200-days. This "sliding calibration window" enables to account for temporal changes in the rheological and hydrological response of the slope. The length of the calibration period was set at 200-days so that: (1) it is long enough to invert a transfer function Γ that matches seasonal variations of the landslide behavior and, (2) this calibration period is short enough to consider that the medium is a time-invariant system and can be characterized by a unique transfer function Γ. The transfer functions Γ i determined thank to the inversion process at date d i are used to predict the displacement rate at dates d i+1 to d i+20 .

The transfer functions Γ i of the FLAME model (Eq.2) are the convolution of a Gaussian function (e.g. representing the recharge of the water reservoir) by an exponential function (e.g. representing the discharge of the water reservoir):

( ) = -ln(2) * - ( ) [2]
For each Γ i , the three degrees of freedom to optimize are: -T≥0: the position of Gaussian (e.g. time delay between the recharge and the rainfall; in day); -D>0: the width at middle height of the Gaussian (e.g. time duration of the phenomenon; in day); -L: the half time duration of the drainage of the water reservoir (in day)

The transfer functions are normalized, considering a coefficient defined by the conservation law between inputs and outputs (Eq.3):

= ∑ ( ) ∑ ( ) [3]

Joint inversion

The displacement rates at the Séchilienne landslide are influenced by the deep aquifer as well as by the perched aquifer. Thus, the transfer function of the model should account for both aquifers and Eq.1 would become, for each sensor:

= . ( * ( ) + ) + . * ( ) + [4]
where β is a normalization constant, and γ and Γ, respectively, are the contribution coefficient and the transfer function of the input E (e.g., the water input), and cst, a constant corresponding to constant contribution. Indices da and pa stand for deep aquifer and perched aquifer, respectively.

This transfer function has 6 degrees of freedom, which have to be optimized during the calibration phase. On a first attempt (case 1), the transfer function Γ is simplified so as to account only for one aquifer:

= . ( * ( ) + ) [5]
where only 3 degrees of freedom are to be optimized during the calibration phase.

On a second attempt (case 2), the complete transfer function Γ is considered:

= . ( * ( ) + ) + . * ( ) + [6]
where 6 degrees of freedom are to be optimized during the calibration phase.

The calibration of these transfer functions can be performed independently for the time series of A13 and A16.

In this study, we also propose to perform the joint inversion of the A13 and A16 time series. The principle of the joint inversion is to assume a similar influence of the deep aquifer on the displacement rates at A13 and A16 (both sensors are located on top of the same hydrogeological unit), whereas such conclusion cannot be derived for the influence of the perched aquifer, described as temporary, with an extent and connectivity dependent on the seasonal recharge. Thus, during the joint inversion, the parameters corresponding to the influence of the deep aquifer will be considered identical for all sensors, whereas the parameters corresponding to the influence of the perched aquifer will be considered different. By doing so, the ratio between the number of parameters to be inverted and the volume of data decreases. The constraint on the inversion process is enhanced.

The joint inversion is performed for the period during which all sensors are operative, from the 24/02/1994, until the 30/06/2013.

Water inputs

The input signal E (see Eq. 4) representing water inputs will be considered as the recharge for the deep aquifer (i.e., the amount of water that is available for infiltration after evapotranspiration), and as the precipitations for the perched aquifer (e.g. the sum of rain and snow melt). The recharge is estimated using the model of Vallet et al. (2015a) (see Fig. 4).

Parameter space

During the inversion process, a large parameter space is investigated in order to optimize the value of the transfer function parameters. This parameter space is chosen so as to cover a large range of realistic values for the Séchilienne landslide (see Tab. 1).

Parameter

Lowest possible value (days)

Highest possible value (days)

T da T pa 0 0 80 90 L da L pa 1 1 80 300 D da D pa 1 1 80 90

Assessment of the quality of the results

The quality of the models will be assessed by analyzing the modelling errors. This error can be characterized by the RMSE of the difference between the observed displacement rate between date d i+1 and date d i+j , and the predicted displacement rate at the same dates; the lower the RMSE, the more accurate the model will be. Furthermore, we assessed the quality of the models computing the coefficient of determination (R²) between the observed and the predicted velocities. The closer to one the R² is, the more accurate the model will be. Finally, a visual inspection of the time-series shape was also carried out in order to investigate whether or not the predictions reproduced the characteristic behavior of the landslide, such as long term trends and seasonal variations.

Simulation of displacement rates at the Séchilienne Landslide

Models using one transfer function Γ da

Fig. 5 presents the results of models using one transfer function (see Eq. 5). Fig. 5a to 5c presents the results of independent inversions for sensors A13 and A16, when Fig. 5d to 5f presents the results of joint inversion for sensors A13 and A16. For both inversions, the predicted /day for mode models using erence, the qu lacements is rsion.

The exp formance with othesis to expl transfer funct to the short cipitations. Th ains higher eafter, we prop ctions, eventho onstraint due t ed displacemen ndent inversion nction of the ob value of R². c) predicted disp d A16.

ween the obser predicted disp 6c and 6f). Despite trained during lacement pred ig. 6, the aver mm/day after mm/day after r of predictio /day for mod sors A13 and A t inversion. A the long-term The quality of the results is clearly improved when using models with two transfer functions rather than models with one transfer function. The average RMSE of predictions gets closer to the measurement error of the sensors (± 0.5 mm/day). This improvement does not allow asserting that the optimized parameters (T,L,D, see Eq. 2) are sufficiently constrained during the inversion and are representative of the averaged behavior of the landslide.

The evolution of the optimized parameters T da ,L da and D da with time is studied for two distinct period P1 and P2 (see During period P2, the displacement rates are more irregular and with more amplitude, together with a general increase of the average velocity. During both periods, the evolutions of the optimized parameters T da , L da and D da are gradual and never reach the limits of the defined parameter space (see Tab. 1). Fig. 7a, 7b and7c shows the distribution of the value of parameters T da , L da and D da for these two periods, as well as their average value. The corresponding shapes of the average transfer function Γ da for periods P1 and P2 are presented on Fig. 7d. Both transfer functions Γ da have realistic shapes, and represent distinct average behaviors that are consistent with the general increase of the average landslide velocity between periods P1 and P2 (the time response decreases between periods P1 and P2).

Discussion

This study discusses the use of statistical models using impulse responses functions to predict landslides velocities. The inversion processes of these statistical models can be poorly constrained when choosing to calibrate a model accounting for a complex hydrogeological response (i.e., with numerous parameters to be determined by inversion). The aim of this study was to test if an enhancement of the method could be achieved by calibrating the models using joint inversion of multiple time series of data.

The method was applied at the Séchilienne landslide where the geometry of the instrumentation meets the requirement of the method. E.g., the displacement rate sensors are located on top of the same hydrogeological unit, thus allowing considering partially identical Impulse Response models for all the inverted time series. By doing so, the ratio between the number of parameters to be inverted and the volume of data decreases, thus a priori increasing the constraint for the search of an optimized solution during inversion.

Our results do not give clear evidence of an improvement of the models performance with joint inversion of multiple time series of data. The reasons which could explain these inconclusive results are numerous. One unfavorable factor is certainly the fact that the data set could not allow a better characterization of the landslide short-term behavior. Data time step is only 24 hours, when the landslide short-term time response is of the order of a few hours (Vallet et al., 2015b). Thus, the convergence of the models to an adequate estimate of the short-term response of the landslide is difficult. Another unfavorable factor is that we assume that the medium is a time-invariant system for the duration of the calibration. If this approximation has a good chance of being correct for the seasonal variations of the landslide, it might not hold true for the short-term response of the landslide. There might be a conflict between the need of a sufficiently long calibration period to constrain the inversion of the long-term response of the landslide, and the necessary assumption to consider the response of the short-term of the landslide as invariant during this calibration period.

Fig. 1

 1 Fig. 1 a) Localization of the Séchilienne landslide on the French Alps map. b) Photograph of the Séchilienne slope in July 2014, with the limits of "Les Ruines" high motion zone in yellow, and the upper limit of the whole unstable mass in white. The positions of extensometers A13 and A16 are indicated in red.

Fig. 3

 3 Fig. 3 Calibration of the model using observed displacement rate befor the date d i ; and prediction of the displacement rate with the calibrated model for dates d i+1 and d i+2 .

Fig. 4

 4 Fig. 4 Daily precipitations at the Mont-sec weather station, and the estimated recharge using the model of Vallet et al. (2015a) at the Séchilienne landslide.

  Fig. 5 Mo comp predic perf (observa

  Fig. 6d, 6g and 6h) for the model using two transfer functions and joint inversion. Period P1 goes from the 06/05/2002 to the 14/08/2003; period P2 goes from the 22/12/2009 to the 19/11/2010. During period P1, the displacement rates have seasonal fluctuations with an increasing trend.

Fig. 6

 6 Fig. 6 Models using two transfer functions (see Eq. 6): a) Observed displacement rates at extensometers A13 and A16 compared to the predicted displacement rates using independent inversions for sensors A13 and A16. b) The predicted displacement rates at A13 and A16 of Fig. 6a as a of the observed displacement rates. Model performance assessment through linear regression and the value of R². c) distribution of the model errors (observations vs. model outcomes). d) to f) same as a) to c) for predicted displacement rates using joint inversions for sensors A13 and A16. g) and h) zooms of Fig. 6d for the periods P1 (from the 06/05/2002 to the 14/08/2003) and P2 (from the 22/12/2009 to the 19/11/2010).

Fig. 7

 7 Fig. 7 a) Distribution and mean value of the T parameter of the tranfer function Γ da during periods P1 and P2 (see Fig. 6d, 6g and 6h). b) and c) same as a) for

Acknowledgments

The monitoring dataset at Séchilienne is part of the Cerema database (Website of Séchilienne rock slope instability:

http://www.versantsechilienne.developpement-durable.gouv.fr/) and the OMIV Observatory (Observatoire Multidisciplinaire des Instabilités de Versants: http://omiv.unistra.fr)