Frédéric Boussinot
email: frederic.boussinot@inria.fr

Frédéric Dabrowski
email: frederic.dabrowski@irisa.fr

Safe Reactive Programming: The FunLoft Proposal

We propose a multicore-ready programming language based on a two-level shared memory model. Concurrency units are schedulers and threads which are dispatched on available cores in a preemptive way. Each scheduler is in charge of its own portion of the memory. At runtime, several threads may link to a common scheduler. In this case, they enter a cooperative mode, evolve in synchronous rounds, and are granted access to the scheduler memory. At the opposite, an autonomous thread runs at its own pace but has access only to a local memory. The language ensures that programs are free of memory leaks, that code between two cooperation points is atomic, and that rounds are fair and always terminate (no run-time error nor divergence).

Introduction

Back in 1960, [START_REF] Hansen | The Origin of Concurrent Programming[END_REF] presents concurrent programming as a major revolution in computer programming. It already points out unexpected concurrent accesses to shared memory as one of the most difficult problems to face. At the language level, appropriate synchronization constructs may help to build reliable software. Such one construct is the now well-known monitor of Brinch Hansen and Hoare. Quoting [START_REF] Hoare | Monitors: an Operating System Structuring Concept[END_REF]: A primary task [...] is to construct resource allocation (or scheduling) algorithms for resources of various kind [...] Each scheduler will consist of a certain amount of local administrative data, together with some procedures and functions which are called by programs wishing to acquire and release resources. Such a collection of associated data and procedures is known as a monitor [...] Multi-threading concurrent programming maps in a very natural way onto multicore architectures. Execution threads are mapped onto available cores and access a common global memory. It is now widely recognized that programming with low-level synchronization mechanisms such as locks is a hard and error prone task. A possible alternative is the use of processes which are concurrent units with separate memory and run-time protection. Processes have been introduced in many concurrent languages as denoting concurrent programs with private memory and structured communication means (for example, in Erlang, Occam, or CML). As opposed to processes, threads do not have built-in memory protection and have direct access to a shared memory. A crucial point is that threads are more efficient than processes, quoting [START_REF] Birrell | An Introduction to Programming with Threads[END_REF] [...] threads are an attractive tool for allowing a program to take advantage of the available hardware. The alternative, with most conventional operating systems, is to configure your program as multiple separate processes, running in separate address spaces. This tends to be expensive to set up, and the costs of communicating between address spaces are often high, even in the presence of shared segments. By using a lightweight multi-threading facility, the programmer can utilize the processors cheaply.

The language we propose, namely FunLoft, relies on a two-level shared memory model. Concurrency units are schedulers and threads dispatched on available cores, in a preemptive way. Informally, at a first glance, one can see schedulers and threads as processes, in the sense that each of them is in charge of its own portion of the memory. At runtime, several threads may link to a common scheduler. In this case, they enter a cooperative mode, evolve in synchronous rounds and are granted access to the scheduler memory. In particular, in this second level, one can see threads as standard threads running in the memory space of a process, except that they execute cooperatively and have access to specific synchronization constructs. When autonomous, a thread runs at its own pace but only has access to its local memory. To guarantee that accesses to memory locations always respect the above-mentioned discipline, we rely on a simple type and effect system mapping memory locations to schedulers and threads and checking access rights. Moreover, we note that, ideally, a language for multi-threaded concurrent programming should meet at least the following requirements:

1. it should have a formal semantics, as intuitive as possible. 2. it should propose high-level constructs for communication and synchronization and these constructs should be efficient, even when there is a large number of concurrent behaviors. 3. it should propose means to define sections of code that should be atomic. 4. it should be safe, in the sense that no run-time error should occur. 5. it should be able to benefit from multi-core architectures.

We propose a language aiming at fulfilling the above-mentioned requirements. In particular, the design of the language ensures that : programs are free of memory leaks, code between two cooperation points is atomic and rounds are fair and always terminate (no run-time error nor divergence).

FunLoft [1], is presented in Section 2. The language is informally described in Section 2.1 and static analyses are considered in Section 2.2. Section 3 describes an example of colliding particles that illustrates the use of multi-core machines. Related work is considered in Section 4, and Section 5 concludes the paper.

FunLoft is an experimental programming language concerned with safety, and qualitative and quantitative resource control in a concurrent setting. Here, safety basically means absence of run-time errors. Resources that are intended to be controlled are the memory and the processing units. Qualitative resource control means control of accesses to the shared memory (see below). Quantitative resource control means absence of memory leaks.

The objective of FunLoft is however larger than to combine concurrency, resource control, and safety: first, the language proposes a particular form of concurrency, as an alternative to standard concurrency frameworks based on preemptive threads (Java threads or Pthreads [START_REF] Nichols | Pthreads Programming[END_REF]). Second, FunLoft allows users to benefit from real parallelism (in particular, the one offered by multi-core machines) when programming a single application. Third, efficiency is an important concern, which was present all along the language design.

Overview of FunLoft

A program in FunLoft is a list of definitions of variables, types, functions, and modules (see below). In order to be executed, a program must contain the definition of a module named main, which is the program entry-point. A variable is a name to which a value is associated. Values are first-order only and range over standard values (booleans and integers, for example), values of user-defined types, and thread, event, and scheduler names. Type definitions introduce polymorphic structured datatypes, possibly recursive (e.g. lists), built from union and concatenation of other types. Functions can be recursively defined; however, recursion is checked for termination, in a sense explained later, and thus function calls always terminate.

Schedulers A scheduler controls the threads which are linked to it and provides them with a portion of shared memory. A special scheduler (the implicit scheduler) is automatically launched by each executable program and a thread created as an instance of the main module is linked to it. The basic task of a scheduler is to schedule the threads which are linked to it. The scheduling is cooperative: linked threads have to return the control to the scheduler, in order to let the other threads execute. Leaving the control can be either explicit, via the instruction cooperate, or implicit, by waiting for an event (await statement) which is not present. All linked threads are cyclically considered by the scheduler until all of them have reached a suspension point (cooperate or await). Then, and only then, a new cycle can start. Cycles are called instants. A scheduler thus defines an automatic synchronization mechanism which forces all the threads it controls to run at the same pace: all the threads must have finished their execution for the current instant before the next instant can start. Note that the same thread can receive the control several times during the same instant; this is for example the case when the thread waits for an event which is generated (generate statement) by another thread later in the same instant. In this case, the thread receives the control a first time and then suspends, waiting for the event. The control then goes to the other threads, and returns back to the first thread after the generation of the event. At creation, each thread is linked to one scheduler (by default, the implicit scheduler). Several schedulers can be defined and simultaneously run in the same program. Schedulers thus define synchronous areas in which threads execute in cooperation. Basically, schedulers run autonomously, in a preemptive way, under the supervision of the OS.

During their execution, threads can unlink from the scheduler to which they are currently linked (unlink statement), and become free from any scheduler synchronization. Such free unlinked threads are, like schedulers, run by kernel threads under the supervision of the OS. Threads can also dynamically move from a scheduler to another scheduler (link statement).

Modules Modules are templates from which threads, called instances, are created. A module can have parameters which define corresponding parameters of its instances. A module also defines variables local to each created thread. As opposed to functions, modules cannot be recursively defined. The body of a module is basically a sequence of instructions with usual control statements such as loop, repeat statements, and conditional statements. There are two types of instructions: atomic instructions and non-atomic ones. Atomic instructions are logically run in a single instant. Function calls belong to this kind of instruction. Execution of non-atomic instructions may need several instants to complete. This is typically the case of the instruction await.

Communication and Synchronisation

The simplest way for threads to communicate is of course to use shared variables. For example, a thread can set a boolean variable to indicate that a condition is true, and other threads can test the variable to know the status of the condition. This basic pattern works well when all the threads accessing the variable are linked to the same scheduler. Indeed, in this case atomicity of the accesses to the variable is guaranteed by the cooperativeness of the scheduler. A general way to protect a data from concurrent accesses is thus to associate it with a scheduler to which threads willing to access the data should first link to 3 . Events are basically used by threads to avoid busy-waiting on conditions. An event is always associated with a scheduler (by default, the implicit scheduler) which is in charge of it during all its lifetime. An event is either present or absent during each instant of the scheduler which manages it. It is present if it is generated by the scheduler (e.g. as a consequence of interactions with the environment) at the beginning of the instant, or if it is generated by one of the threads executed by the scheduler during the instant; it is absent otherwise. The presence or the absence of an event can change from an instant to another, but all threads always "see" the presence or absence of events in the same way, independently of the order in which the threads are scheduled. For this reason, we say that events are broadcast. Values can be associated with event generations; they are collected during each instant (get all values statement) and their collection becomes available as a list at the next instant.

It is possible to define synchronised schedulers that are run autonomously but that share the same instants. Events are shared by synchronised schedulers, which is possible because synchronised schedulers themselves share instants. Thus, if s1 and s2 are two synchronised schedulers, it is for example possible to wait for an event in s1 and to generate the same event in s2.

Memory Model The global memory is divided in several parts:

-A private memory for each thread. This memory is initialised when the thread is created, from the parameters and the local variables of the module from which the thread is created. The system checks that no other thread can have access to this private memory. -A private memory for each scheduler. This memory can only be accessed by the threads linked to the scheduler. The system checks that it is not accessible from unlinked threads, nor from threads linked to another scheduler.

Thus, there exists no global variable shared by distinct schedulers, and a variable that is shared by several threads belongs to the unique scheduler that controls these threads and cannot be accessed by other threads. This is the way data-races are ruled-out in FunLoft.

Formalisation

Despite its two concurrency levels, FunLoft has a formal semantics which is very concise. A semantics restricted to a unique scheduler is presented in [START_REF] Dabrowski | Cooperative Threads and Preemptive Computations[END_REF]; a complete semantics, when there are several schedulers, is described in [START_REF] Dabrowski | Programmation Réactive Synchrone[END_REF] for a language very close to FunLoft, and the full formalisation of Funloft is described in a technical paper [START_REF] Boussinot | Formalisation of FunLoft[END_REF]. Among the properties checked at compiletime in FunLoft programs, are the following:

1. the program is well-typed under standard definitions 2. functions always terminate and linked threads always cooperate (despite functions are recursive and threads may never terminate). 3. The amount of memory used by the program and the number of (simultaneously) living threads are bound by a function of the program inputs (despite the use of complex datatypes and the dynamic creation of threads). 4. the memory model described above is actually respected; more precisely the compiler checks that there exists a mapping from memory locations to schedulers and threads such that any execution respects the corresponding access rights.

In particular, (1), ( 2) and (3) ensure that a thread can neither produce a runtime error, nor prevent the termination of instants, and thus delay endlessly the handling of further inputs. Thus, a thread can never prevent other threads from progressing.

Property (4) ensures that, when linked to a scheduler, a thread always executes atomically (at least from a logical point of view) the portions of code between two cooperation points.

Details on the formalisation of FunLoft are out of the scope of this paper. We refer the interested reader to [START_REF] Dabrowski | Programmation Réactive Synchrone[END_REF][START_REF] Boussinot | Formalisation of FunLoft[END_REF]. Ongoing work on computing explicit bounds on the amount of memory used by the program can also be found in [START_REF] Amadio | Feasible Reactivity in a Synchronous pi-Calculus[END_REF] (in FunLoft, we only check for the existence of bounds, and we do not compute them explicitly). However, here are some hints about the static analyses performed by the compiler.

Basically, the compiler relies on a type and effect system mapping disjoint subsets of memory locations to owners (schedulers and threads). It enforces a logical separation of the memory as follows:

1. when autonomous, a thread can only access its own memory, 2. when linked to a scheduler, a thread can only access its own memory and the memory of the scheduler it is linked to.

As linked threads behave cooperatively, it is immediate that the execution of code between two cooperation points is atomic. Termination of recursive functions is controled. In the present version of FunLoft, recursivity can only concern parameters of inductive types. The size of a parameter p is smaller than the size of a parameter q if p is a sub-term of q. For lists of parameters, one extends "lexicographically" the notion of size. In a complete sequence of calls (starting and ending on calls of the same function), one checks that at each step parameters are strictly decreasing. Thus, any function call is forced to terminate after a finite number of recursive calls.

To control the size of the global memory used by a program, memory locations are stratified over a finite domain. We rely on a simple data flow analysis to ensure that no cyclic dependence is possible for memory locations carrying complex datatypes (simple datatypes such as integers do not have this restriction). A value put into a memory location is always the result of the composition of a finite number of functions applied to values read from upper memory locations. As functions always terminate, the size of values held at a memory location is bound by a function of the size of the program inputs.

Finally, we rely on a simple static analysis to guarantee that the number of living threads remains bounded. Actually, we avoid (asynchronous) thread creation in loops. Then, at any point of execution, the number of living threads only depends on the size of the program inputs. Thanks to the design of the language, at the beginning of each instant, the size of the stack associated with a thread is bound by a constant. This entails that (given the existence of a garbage collection mechanism) the size of the global memory used by a program is a function of the program inputs.

Several other static controls are also made. For example, an event cannot be generated by a thread linked to a scheduler, and waited for by a thread linked to another scheduler, unless the two schedulers are synchronised (otherwise, a buffer to store generations would be necessary, and, as the schedulers would be asynchronous, the buffer size could not be bounded).

Implementation

We are currently developing an implementation prototype, available at [1]. Autonomous unlinked threads and schedulers co-existing in the same application are run by distinct kernel processing units (kernel threads) on distinct processors (multi-processors or multi-core architectures).

This prototype is used to implement several examples of graphical simulations, including cellular automata made of several thousands of cells, each cell being a thread. The simulation presented in Section 3 is one of these examples. These examples show, in particular, that systems with large numbers of concurrent programs can be efficiently implemented in FunLoft.

Colliding Particles

One considers the simulation of colliding particles. A simple brute-force algorithm is used in which each particle considers all other particles (leading thus to a complexity square in the number of particles). The particles are partitioned in two sets of equal size, and each set is linked to a distinct scheduler. In order to make the collisions realistic, the two schedulers are synchronised.

Collision Processing

A particle is made of four references on floats and of a constant color. The data-type particle t is defined by: The process all collisions function processes collisions between a particle given as first parameter, and a list of particles given as second parameter (collision between two particles is performed by the function collide not described here for simplicity). The parametric definition of the inductive type of lists is:

type 'a list = Nil_list | Cons_list of 'a * 'a list
The recursive definition of process all collisions is: let process_all_collisions (me,list) = match list with Nil_list -> () | Cons_list (other,tail) -> begin VIII collide (me,other); process_all_collisions (me,tail); end end

The system verifies that the function always terminates. This results from the fact that the second parameter (tail) of the recursive call is smaller than (is a subterm of) the second parameter (list) of the initial call (the first parameter remaining unchanged).

Particle Behaviour

At each instant, a particle generates the collide event with an associated value containing its coordinates, and then collects all the values generated on this very event. Thus, at the next instant, the coordinates of all the particles have been collected. They are then processed, and the particle is finally moved according to inertia. The code of module collide behavior is: The system checks that the body of the loop statement, which defines a cyclic behavior, cannot terminate instantly. Otherwise, the execution of an instance of the module would not cooperate and would prevent all other threads from execution. The absence of immediate termination comes from the fact that the collection of all the values takes the whole instant; thus, the processing of the collected values is delayed to the next instant.

let module collide_behavior (me,collide_event) = let r = ref
The particle behavior is actually made of three threads: one for bouncing on the borders of the simulation (bounce behavior), one for colliding with other particles (collide behavior), and one for drawing the particle on the screen (draw behavior): let module particle_behavior (collide_event,color) = let s = new_particle (color) in begin thread bounce_behavior (s); thread collide_behavior (s,collide_event); thread draw_behavior (s); end Note that the particle s (created by the function new particle in charge of randomly placing the particle on the screen) is shared between the three threads.

Simulation

First, two schedulers s1 and s2 are declared at top-level (keyword scheduler). The presence of the keyword and which links the two declarations makes the two schedulers synchronised. Then, two events draw event and collide event are declared (keyword event). The system will infer that these events are associated to the two schedulers s1 and s2. Then, after linking to scheduler s1, a thread instance of the module graphics and one of the module draw processor are created and linked to s1. In the same way, half of the particles are created and linked to s1. Finally, after linking to s2, the other half of particles are created and linked to s2. The main module is: The compiler verifies that no data-race can appear due to the sharing of particles: a particle is only shared by threads linked to the same scheduler. It also verifies that the two events draw event and collide events are always used by threads linked to s1 or to s2.

let s1 =

Execution Results

The machine characteristics are: a Mac running OS X 10.4.10, processor Intel Core 2 Duo, 2.33 GHz, 2GB of memory. Graphics is based on SDL [START_REF]Simple Directmedia Layer[END_REF].

The simulation contains 500 particles. The CPU usage is shown in the right part of Figure 1 (the graphical window is masked during the measure). This is to compare with the usage obtained with the single threaded version of the simulation, shown in the left part.

The time for simulating 100 instants (one instant corresponds to the execution of all the particles) is shown on the left part of Figure 2 (representing the output of the unix command time), with a memory footprint of about 70MB. Note that the total number of performed interactions is about 100 * 500 2 = 25 6 . The case of 1000 particles is shown on the right of Figure 2. These results show the benefit from using a multi-core computer on a simulation in with an heavy load of communication between threads.

Related work

FunLoft is strongly connected to synchronous languages [START_REF] Halbwachs | Synchronous Programming of Reactive Systems[END_REF]. In these languages, as in FunLoft, concurrent behaviors share a logical clock, the units of which are called "instants". Concurrent behaviors automatically synchronise at each instant and communicate using broadcast signals to which data can be associated. At each instant, a signal is consistently seen as present or absent by all the parallel processes, which also read the same associated data. In Esterel [START_REF] Berry | The Esterel Synchronous Programming Language: Design, Semantics, Implementation[END_REF], one of the first synchronous languages, reactivity (termination of the instants) and determinism are guaranteed for the kernel of the language. However, this kernel has deep constraints that restrict the language to applications such as the production of circuits.

Strongly related to synchronous languages is a family of formalisms grouped under the name "Reactive Programming" [START_REF]Reactive Programming[END_REF]. Reactive programming is implemented in various languages (SML, Java, ML, C, and Scheme). The formalism of FairThreads [START_REF] Boussinot | FairThreads: Mixing Cooperative and Preemptive Threads in C[END_REF], built over C, deserves a special status amongst the reactive programming formalisms for two reasons: first, in FairThreads, concurrent programs are basically threads, while the other formalisms use parallel composition, adapted for instants; second, FairThreads considers a mix of cooperative and preemptive threads, which gives a model in which synchrony and asynchrony are both present. FairThreads aims at being a more general purpose programming language than its ancestors, but loses reactivity and determinism. The FunLoft language is very close to FairThreads as regards the concurrency model. However, although the FairThreads documentation suggests that schedulers might be used to access shared memory in a cooperative way, no static checks ensure that threads not linked to a given scheduler (possibly linked to other schedulers) do not interfere with the threads linked to this scheduler.

There exists a line of experimental work aiming at obtaining both language expressivity and strong properties such as reactivity and determinism. For reactivity, we refer to [START_REF] Amadio | Feasible Reactivity in a Synchronous pi-Calculus[END_REF] for results in a synchronous π-calculus [START_REF] Amadio | A Synchronous pi-Calculus[END_REF]. In this paper, a static analysis is defined that guarantees that instants always terminate in time polynomial in the size of the program inputs. It extends preliminary results that can be found in [START_REF] Amadio | Resource Control for Synchronous Cooperative Threads[END_REF][START_REF] Amadio | Feasible Reactivity for Synchronous Cooperative Threads[END_REF]. For determinism, we refer to [START_REF] Amadio | Determinacy in a Synchronous pi-Calculus[END_REF] for results applying to the same calculus. However, the synchronous π-calculus does not consider shared memory: in it, valued signals are the only means of communication between concurrent behaviors.

The approach followed in FunLoft is very close to the work on static analysis for data races detection [START_REF] Flanagan | Types for Safe Locking[END_REF][START_REF] Flanagan | Type-based Race Detection for Java[END_REF][START_REF] Flanagan | Type Inference Against Races[END_REF][START_REF] Boyapati | A Parameterized Type System for Race-free Java Programs[END_REF][START_REF] Naik | Effective Static Race Detection for Java[END_REF][START_REF] Naik | Conditional Must not Aliasing for Static Race Detection[END_REF] and to the work on checking for atomicity in a lock-based concurrent paradigm [START_REF] Flanagan | Types for Atomicity[END_REF][START_REF] Flanagan | A Type and Effect System for Atomicity[END_REF]. Our point of view is that locks are a too low-level tool for general purpose concurrent programming. We hope that the programming model of FunLoft will provide a higher-level model of concurrency, thus helping in developing and debugging applications.

Conclusion

At the implementation level, threads allow users to benefit from multi-core architectures and from shared memory. At the logical level, one needs to ensure atomic execution of code blocks. We propose a framework based on FairThreads in which static checks insure three fundamental properties: that programs are free of memory leaks, that code between two cooperation points is atomic, and that rounds are fair and always terminate (no run-time error nor divergence). In this framework, fair threads can be as efficient as standard preemptive threads, but as secure as processes. The multi-scheduler aspect of the proposal fits well with multi-core architectures: schedulers and unlinked threads are executed by kernel threads that can be securely run in parallel by distinct cores without needing dynamic memory protection mechanisms.

Fig. 1 .Fig. 2 .

 12 Fig. 1. CPU usage. Left: 1 scheduler. Right: 2 schedulers

Note the analogy with monitors, described indeed by Hoare as schedulers, in the citation of Section 1.


Guy Perrier
email: perrier@loria.fr

A French Interaction Grammar

Keywords: Syntax, grammatical formalism, tree description, polarity, categorial grammar, unification grammar, interaction grammar

We present a relatively large coverage French grammar written with the formalism of Interaction Grammars. This formalism combines two key ideas: the grammar is viewed as a constraint system, which is expressed through the notion of tree description, and the resource sensitivity of natural languages is used as a syntactic composition principle by means of a system of polarities. We give an outline of the expressivity of the formalism by modelling significant linguistic phenomena and we show that the grammar architecture provides for re-usability and tractability, which is crucial for building large coverage resources: a modular source grammar is distinguished from the object grammar which results from the compilation of the first one, and the lexicon is independent of the grammar. Finally, we present the results of an evaluation of the grammar achieved with the LEOPAR parser with a test suite of sentences.

Introduction

The goal of our work is to model natural languages starting from linguistic knowledge and giving a central role to experimentation. For this, we need to express the linguistic knowledge by means of grammars and lexicons with the largest possible coverage: grammars have to represent all common linguistic phenomena and lexicons have to include the most frequent words with their most frequent use. As everyone knows, building such resources is a very hard task.

Firstly, we have to choose the formalism to represent the grammar. Currently, there is no leader among the formalisms used in the scientific community. Each of the most popular formalisms has its own advantages and drawbacks. We have designed a new formalism, Interaction Grammars (IG), the goal of which is to synthesize two key ideas, expressed in two kinds of formalisms up to now: using the resource sensitivity of natural languages as a principle of syntactic composition, which is a characteristic feature of Categorial Grammars (CG) [START_REF] Retoré | The Logic of Categorial Grammars[END_REF], and viewing grammars as constraint systems, which is a feature of unification grammars such as LFG [START_REF] Bresnan | Lexical-Functional Syntax[END_REF] or HPSG [START_REF] Sag | Syntactic Theory: a Formal Introduction[END_REF].

Although we use an original formalism, we are concerned with re-usability, which is expressed in two ways. Like with for programming languages, we distinguish two levels in the grammar. The source grammar aims at representing linguistic generalisations and it is written by a human, while the object grammar is directly usable by a NLP system and results from the compilation of the first one. In our case, we used XMG [START_REF] Duchier | XMG : Un compilateur de méta-grammaires extensible[END_REF], a tool devoted to this goal. XMG provides a high level language for writing a source grammar and a compiler which translates this grammar into an operational object grammar. The grammar is also designed in such a way that it can be linked with a lexicon independent of the formalism, where entries appear as feature structures.

The goal of the article is to show that it is possible to build realistic grammatical resources, which integrate a refined linguistic knowledge with a large coverage, and for this, we have chosen an experimental approach with the construction of a French grammar.

Interaction Grammars

IG [START_REF] Perrier | Interaction grammars[END_REF][START_REF] Perrier | La sémantique dans les grammaires d'interaction[END_REF] is a grammatical formalism which is devoted to the syntax and semantics of natural languages and which uses two notions: tree description and polarity.

Tree Descriptions

In a derivational view of the syntax of natural languages, the basic objects are trees and they are composed together in a more or less sophisticated way: by substitution in Context Free Grammars, by adjunction in Tree Adjoining Grammars, by application and abstraction in Categorial Grammars . . . Taking our view from the Model Theory [START_REF] Pullum | On the Distinction between Model-Theoretic and Generative-Enumerative Syntactic Frameworks[END_REF], we do not directly manipulate trees but properties which are used to describe them, in other words tree descriptions [START_REF] Rogers | Obtaining trees from their descriptions: an application to tree-adjoining grammars[END_REF]. This approach is very flexible as it allows the expression of elementary properties in a totally independent way, as they can be freely combined.

A tree description can be viewed either as an underspecified tree, or as the specification of a tree family, each tree being a model of the specification. Figure 1 gives an example of a tree description, which is associated with the relative pronoun qui (who), used inside a prepositional complement. This use gives rise to the phenomenon of pied piping as the following example illustrates: Jean [à la femme de qui] Pierre sait qu'on a présenté Marie , est ingénieur (Jean [to whose wife] Pierre knows someone presented Marie , is an engineer). This example is covered by the description of figure 1. 1A tree description is a finite set of nodes structured by two kinds of relations: dominance and precedence. Dominance relations can be immediate or large (respectively solid and dashed down arrows in figure 1). Constraints can be put on intermediate nodes for large dominance relations. Precedence relations (horizontal arrows in figure 1) can also be immediate or large.

Nodes, which represent constituents, are labelled with features describing their morpho-syntactic properties. Feature values are atoms or atom disjunctions and they can be shared with the help of a co-indexation mechanism. 2 Nodes can be Empty (the white box in figure 1) or Full, according to whether they have an empty phonological form or not. Full nodes can be Anchors (the dark box in figure 1), if they anchor a word of the language.

Polarities

Polarities are used to express the saturation state of syntactic trees.They are attached to features that label description nodes with the following meaning:

• a positive feature t → v expresses an available resource, which must be consumed;

• a negative feature t ← v expresses an expected resource, which must be provided; it is the dual of a positive feature;

• a neutral feature t = v expresses a linguistic property that is not a consumable resource.

• a virtual feature t ∼ v expresses a linguistic property that needs to be realised by combining with an actual feature (an actual feature is a positive, negative or neutral feature).

In figure 1, the empty node representing the trace of the prepositional phrase extracted from the relative clause carries a positive feature cat → pp and a negative feature f unct ← 1 ?, which means that this node provides a prepositional phrase that needs to receive a syntactic function. The tree root carries a virtual feature cat ∼ np which means that the node represents a virtual noun phrase which has to combine with an actual noun phrase.

The descriptions labelled with polarised feature structures are called polarised tree descriptions (PTDs) in the rest of the article.

Grammars as constraint systems

A particular interaction grammar is defined by a finite set of elementary PTDs, which generates a tree language. A tree belongs to the language if it is a model of a finite set of elementary PTDs with two properties:

• It is saturated : every positive feature t → v is matched with its dual feature t ← v in the model and vice versa. Moreover, every virtual feature has to find an actual corresponding feature in the model.

• It is minimal : the model has to add a minimum of information to the initial descriptions (it cannot add immediate dominance relations or features that do not exist in the initial descriptions).

Then, parsing reduces to the resolution of a constraint system. It consists of building all saturated and minimal models of a finite set of elementary PTDs. In practice, our grammar is totally lexicalized: each elementary PTD has a unique anchor, which is used for linking the description with a word of the language. In this way, in the parsing of a sentence, it is possible to select the only PTDs that are anchored by words of the sentence. The set of PTDs being selected, the building of a saturated and minimal model is performed step by step by means of a merging operation between nodes, which is guided by one of the following constraints:

• neutralise a positive feature with a negative feature having the same name and carrying a value unifiable with the value of the first feature;

• realise a virtual feature by combining it with an actual feature (a positive, negative or neutral feature) having the same name and carrying a value unifying with the value of the first feature.

The constraints of the description interact with node merging to entail a partial superposition of their contexts represented by the tree fragments in which they are situated. To summarise, IG combine the strong points of two families of formalisms: the flexibility of Unification Grammars and the saturation control of Categorial Grammars.

[2] Jean cat = np funct = ? [START_REF] Bresnan | Lexical-Functional Syntax[END_REF] cat -> np funct <-? [START_REF] Nasr | A formalism and a parser for lexicalised dependency grammars[END_REF] cat ~v mood = ind

[7] cat ~aux | v mood = ind [6] la cat = clit [5] cat -> np funct <-attr | obj [3] cat ~s mood = ind [11] cat <-np funct -> obj [8] cat = s mood = ind [9] cat <-np funct -> subj [10] cat = v mood = ind tense = pres [12] voit cat = v mood = ind tense = pres [[2]] Jean cat = np funct = subj [[1-9]] cat = np funct = subj [[4-10]] cat = v mood = ind tense = pres [[7-12]] voit cat = v mood = ind tense = pres [[6]] la cat = clit [[5-11]] cat = np funct = obj [[3-8]] cat = s mood = ind
Fig. 2: PTD associated with the sentence Jean la voit and its minimal saturated model Figure 2 presents an example of parsing for the sentence Jean la voit (Jean sees her). 3 The left side shows the set of initial PTDs associated with the sentence by the grammar. The grammar being lexicalized, each PTD is anchored by a word of the sentence and it has been extracted from a lexicon. These PTDs have been gathered in a unique PTD and precedence relations between anchors have been added to express word order in the sentence. These relations do not appear in figure 2.

The computation of the model shown on the right side of figure 2 from the initial description shown on the left side is performed by a sequence of 3 node mergings. 4 The interaction of tree constraints with these mergings entails two other mergings and a partial tree superposition.

The expressivity of Interaction Grammars

In the limits of this article, we have chosen to illustrate three aspects which are especially significant.

Unbounded dependencies and underspecified dominance relations

Underspecified dominance relations are used to represent unbounded dependencies and the feature structures that can be associated with these relations allow the expression of constraints on these dependencies: barriers to extraction for instance.

Relative pronouns, such as qui or lequel, give rise to pied piping as the following sentence shows: Jean [dans l'entreprise de qui] Marie sait que l'ingénieur travaille , est malade (Jean [in whose firm] Marie knows that the engineer works , is ill):

• There is a first unbounded dependency between the verb travaille and its extracted complement dans l'entreprise de qui. The trace of the extracted complement is denoted by the symbol. The dependency is modelled in the PTD associated with the qui relative pronoun represented in figure 1 by means of an underspecified dominance relation. The constraint linked to this dominance relation expresses that the dependency of the prepositional phrase on the verb of which it is the complement can only cross an unspecified sequence of embedded object clauses.

• Inside the prepositional phrase, there is a second unbounded dependency between the head of the constituent and the qui relative pronoun, which can be embedded arbitrarily deeply. This dependency is also represented in figure 1 with an underspecified dominance relation and the linked constraint expresses that all embedded constituents from the prepositional phrase to the qui relative pronoun are common nouns, noun phrases or prepositional phrases.

Polarities used for modelling negation

In French, negation can be expressed with the help of the particle ne paired with a specific determiner, pronoun or adverb. The position of the particle ne is fixed before an inflected verb but the second component of the pair, if it is a determiner like aucun or a pronoun like personne, can have a relatively free position in the sentence, as illustrated by the following examples: (a) Jean ne parle à aucun collègue (Jean speaks to no colleague). (b) Jean ne parle à la femme d'aucun collègue (Jean speaks to the wife of no colleague). (c) Aucun collègue de Jean ne parle à sa femme (No colleague of John's speaks to his wife). As figure 3 shows, the pairing of ne with aucun is expressed with a neg polarised feature attached to the node representing the maximal projection of the verbal kernel: aucun is waiting for such a feature, which will be provided by ne. The relatively free position of aucun is expressed by an underspecified dominance relation of the node representing the clause on the noun phrase that it introduces. The constraint linked to this dominance relation expresses the fact that aucun can only introduce arguments of the verbal head of the sentence or complements of these arguments.

The adjunction of modifiers by means of virtual polarities

In French, the position of adjuncts in the sentence is relatively free, as illustrated by the following example.

In the sentence Jean va rendre visite à Marie (Jean is going to visit Marie), the sentence modifier le soir (tonight) can appear at any position marked with a symbol, according to different communicative goals.

The virtual polarity f ∼ v did not exist in the previous version of IG [START_REF] Perrier | La sémantique dans les grammaires d'interaction[END_REF]. Modifier adjunction was performed by addition of a new level in the syntactic tree of the constituent being modified. Sometimes, introducing an additional level is justified linguistically, but in most cases it introduces artificial complexity and ambiguity. Taking again an idea of [START_REF] Nasr | A formalism and a parser for lexicalised dependency grammars[END_REF], with his system of black and white polarities, we have introduced virtual polarities. This allows a modifier to be added as a new daughter of the node that it modifies without changing the rest of the syntactic tree, in which the modified node is situated. This operation is called sister adjunction and it is used in some formalisms: dependency grammars, description substitution grammars [START_REF] Rambow | D-tree substitution grammars[END_REF]. This way of modelling modifiers is more flexible and it allows the previous examples to be treated without difficulty, including parenthetical clauses. [START_REF] Nasr | A formalism and a parser for lexicalised dependency grammars[END_REF] The architecture of the grammar

The modular organisation of the grammar

The grammar has been built with the XMG tool [START_REF] Duchier | XMG : Un compilateur de méta-grammaires extensible[END_REF], which allows grammars to be written with a high level of abstraction in a modular setting and to be compiled into low level grammars, usable by NLP systems. A grammar is organised as a class hierarchy by means of two composition operations: conjunction and disjunction. It is also structured according to several dimensions, which are present in all classes. Our grammar uses only two dimensions: the first one is the syntactic dimension, where objects are PTDs, and the second one is the dimension of the interface with the lexicon, where objects are feature structures.

To define the conjunction of two classes one needs to specify the way of combining the components of each dimension: for the syntactic dimension, PTD union is performed; for the lexicon interface dimension, it is realised as unification between feature structures. The current grammar is composed of 448 classes, including 121 terminal classes, which are compiled into 2059 PTDs. These classes are ranked by family. Some classes from a family can be used in the definition of classes belonging to another family. This is the case for instance for the Complement family, which include classes related to complements of predicative structures. It is used by three other families: Adjective, Noun and VerbDiathese, which respectively refer to adjectives, nouns and various verbal diatheses.

The link with a lexicon independent of the formalism

The grammar, in its current setting, is totally lexicalised: each elementary PTD of the grammar has a unique anchor node intended to be linked with a word of the language. Each PTD is associated to a feature structure, which describes a syntactic frame corresponding to words able to anchor the PTD, the description being independent of the formalism. This feature structure constitutes the PTD interface with the lexicon. The set of features used in the interfaces differs from that used in PTDs because they do not play the same role: they do not aim at describing syntactic structures but they are used for describing the morpho-syntactic properties of the words of the language in a way independent of the formalism.

The left side of figure 4 shows a non anchored PTD describing the syntactic behaviour of a transitive verb in the active voice. The PTD is accompanied by its interface, which is a two level feature structure.

The lexicon associates words of the language to syntactic frames in a form identical to the PTD interfaces. For instance, the central part of figure 4 shows a lexical entry for the verb voit in its transitive use.

The PTD anchoring is then performed by unification of the PTD interfaces with the compatible entries of the lexicon. There are not many corpora of this type for French. We have chosen the TSNLP [START_REF] Lehmann | tsnlp -Test Suites for Natural Language Processing[END_REF], which includes 1690 positive sentences and 1935 negative sentences. It is far from covering all of French grammar; in particular, it includes very few complex sentences but it stresses some phenomena such as coordination or the postion in the sentence of the adverbial complements. On the other hand, our grammar covers phenomena that are ignored by the TSNLP: the passive and middle voice of verbs, the subcategorisation of predicative nouns and adjectives, the control of the subject of infinitive complements, the relative and interrogative clauses. . .

For the parsing, we used LEOPAR 5 , which is a parser devoted to IG. With the current grammar, the parser accepts 88% of the 1690 positive TSNLP sentences and rejects 85% of the 1935 negative sentences. The 15% of accepted negative sentences are due to the fact that the grammar ignores phonological rules and semantics. The 12% of unanalysed positive sentences are due to various reasons: speech sentences, frozen or semi-frozen expressions, phenomena that are not yet taken into account (causatives, superlatives. . .).

Prospects

The next step is to use our French grammar to parse raw corpora. It is already possible to use LEOPAR 5 http://www.loria.fr/equipes/calligramme/leopar with a large lexicon for such a task. It is necessary to enrich the grammar because some common linguistic phenomena are not yet taken into account. We also need to improve the efficiency of the parser to contain the possible explosion resulting from the increase of the grammar size in combination with the increased sentence length.

Fig. 1 :

 1 Fig. 1: Tree description associated with the relative pronoun qui used inside a prepositional complement

Fig. 3 :

 3 Fig. 3: PTDs respectively associated with the particle ne and the determiner aucun

  Figure 4 on its right side shows a PTD anchored by the transitive verb voit. This PTD comes from the unification between the lexical entry for voit presented in the center of the figure and the interface of the non anchored PTD on the left side of the figure.

Fig. 4 :

 4 Fig. 4: From left to right, a non anchored PTD describing the syntactic behaviour of a transitive verb in the active voice, a lexical entry for the transitive verb voit and the PTD after anchoring with the verb voit

The extracted prepositional phrase is put between square brackets and its trace in the relative clause is represented by the symbol.

When two features share the same value, a common index n is put before their values. When a feature value is the disjunction of all elements of a domain, this value is denoted with "?".

We have simplified the figure by ignoring agreement features.

The head of each node includes the numbers of the nodes from the initial PTD which have been merged.


Léo Allemand-Giorgis

Georges-Pierre Bonneau

Stefanie Hahmann

G.-P Bonneau

Piecewise Polynomial Reconstruction of Scalar Fields from Simplified Morse-Smale Complexes

Morse-Smale complexes have been proposed to visualize topological features of scalar fields defined on manifold domains. Herein, three main problems have been addressed in the past: (a) efficient computation of the initial combinatorial structure connecting the critical points; (b) simplification of these combinatorial structures; (c) reconstruction of a scalar field in accordance to the simplified Morse-Smale complex. The present paper faces the third problem by proposing a novel approach for computing a scalar field coherent with a given simplified MS complex that privileges the use of piecewise polynomial functions. Based on techniques borrowed from shape preserving design in Computer Aided Geometric Design, our method constructs the surface cell by cell using piecewise polynomial curves and surfaces. The benefit and limitations of using polynomials for reconstruction surfaces from topological data are presented in this paper.

Introduction

Morse-Smale (MS) complexes are combinatorial structures encoding topological features of smooth scalar fields by connecting their critical points. The discrete Morse theory [START_REF] Forman | A user's guide to discrete morse theory[END_REF] has been used to define MS complexes for scalar data defined on meshes.

Because of their ability to characterize the topology of level sets, MS complexes have found many applications in visualization [START_REF] Gyulassy | Direct feature visualization using Morse-Smale complexes[END_REF]. Each cell in a MS complex is a set of integral lines of a discrete gradient flow (the equivalent of streamlines of a smooth gradient field) connecting a critical point with some value to another critical point with a higher value. Along an integral line, the data is monotonic, strictly increasing. Inside a 2d-cell, there are no critical points, and the level sets have the topology of an open interval. Previous works have their focus on the efficient computation of MS complexes [START_REF] Shivashankar | Parallel computation of 3D Morse-Smale complexes[END_REF][START_REF] Shivashankar | Parallel computation of 2D Morse-Smale complexes[END_REF], on the simplification and hierarchical representation of MS complexes [START_REF] Bremer | A topological hierarchy for functions on triangulated surfaces[END_REF][START_REF] Edelsbrunner | Hierarchical Morse complexes for piecewise linear 2-manifolds[END_REF] or on the extraction of features using MS complexes [START_REF] Gyulassy | Direct feature visualization using Morse-Smale complexes[END_REF]. The simplification of a MS complex is a combinatorial operation consisting in deleting two adjacent 0d-cells in the complex. All adjacent cells to the deleted nodes affected by the simplification are either deleted or unified. Since the simplification is a purely combinatorial operation, for visualization purposes it is thus often required to compute new scalar data that correspond to the simplified complex. The scalar field reconstruction problem is the topic of the present work. The main requirements on the scalar field reconstructed from a simplified MS complex are the following • interpolation of the scalar field values at the critical points • bijection between the nodes of the MS complex and the critical points • bijection between arcs of the MS complex and 1-cells. This problem has been solved in the past by applying to each 1d-and 2d-cell a discrete laplacian smoothing in the original mesh with the addition of constraints in order to ensure the monotonicity of data inside the cells [START_REF] Bremer | A topological hierarchy for functions on triangulated surfaces[END_REF][START_REF] Weinkauf | Topology-based smoothing of 2D scalar fields with C 1 -continuity[END_REF]. In contrary to previous mesh-based approaches, the present work investigates the use of smooth surface representations.

In the fields of Computer Aided Geometric Design and Approximation Theory shape preserving methods seek to approximate a scalar or geometric function while preserving some "shape properties" such as the sign, the monotonicity or the convexity of the given data. Most of these methods produce piecewise polynomial Bézier and B-spline functions, parametric curves and surfaces.

In the present work we propose to reconstruct scalar fields from simplified MS complexes in 2D, based on shape-preserving methods. Instead of using discrete laplacian smoothing on the original mesh, our approach is to use piecewise polynomial curves and surfaces. Specifically, we first approximate data along 1dcells by computing monotonic B-splines [START_REF] He | Monotone B-spline smoothing[END_REF], and we then show how to fill-in 2d-cells with monotonic triangular Bézier patches based on our previous work on monotonic interpolation of gridded data [START_REF] Allemand-Giorgis | Piecewise polynomial monotonic interpolation of 2D gridded data[END_REF]. The final resulting scalar field is a piecewise polynomial approximation of the original data that is coherent with the simplified MS complex. We present different approaches for parameterizing and computing these piecewise polynomial surfaces.

Related Work

Topology-based / Computational topology methods have become very effective for analysis, visualization and simplification of features in scientific data sets and geometric shape models. Most of these methods make use of results from Morse theory which dates back to the nineteenth century. It became useful in scientific visualization and computer graphics with a first efficient algorithm by Edelsbrunner et. al [START_REF] Edelsbrunner | Hierarchical Morse complexes for piecewise linear 2-manifolds[END_REF] for computing the MS complex for piecewise linear scalar fields. The MS complex is a topological representation of a scalar function which decomposes the domain of a function into regions having uniform gradient flow. It consists of the critical points (minima, maxima, saddles) and the separatrices connecting them.

A u t h o r p r e p r i n t

Efficient algorithms to compute the MS complex for 2D and 3D scalar fields can be found in [START_REF] Gyulassy | A topological approach to simplification of three-dimensional scalar functions[END_REF][START_REF] Gyulassy | Efficient computation of Morse-Smale complexes for three-dimensional scalar functions[END_REF][START_REF] Gyulassy | A practical approach to Morse-Smale complex computation: scalability and generality[END_REF][START_REF] Gyulassy | Computing Morse-Smale complexes with accurate geometry[END_REF][START_REF] Shivashankar | Parallel computation of 2D Morse-Smale complexes[END_REF].

Topological simplification is used for denoising, feature discrimination and smoothing of scalar fields. Algorithms apply, inter alia, the notion of topological persistence [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] for iteratively cancelling pairs of critical points [START_REF] Bremer | A topological hierarchy for functions on triangulated surfaces[END_REF][START_REF] Edelsbrunner | Hierarchical Morse-Smale complexes for piecewise linear 2-manifolds[END_REF] in MS complexes.

Closely related to topology-driven simplification algorithms are methods that reconstruct a function obeying the simplified topological data. A setup similar to ours are [START_REF] Bremer | A topological hierarchy for functions on triangulated surfaces[END_REF][START_REF] Weinkauf | Topology-based smoothing of 2D scalar fields with C 1 -continuity[END_REF] in the sense that the MS complex is simplified and a corresponding scalar function is computed by fully constraining the simplified complex. In [START_REF] Bremer | A topological hierarchy for functions on triangulated surfaces[END_REF] a MS complex is iteratively simplified by cancelling pairs of critical points ordered by persistence resulting in a hierarchical representation. The discrete embedding of boundary and the interior of the cells is smoothly reconstructed cell by cell using constrained Laplacian optimization. In [START_REF] Weinkauf | Topology-based smoothing of 2D scalar fields with C 1 -continuity[END_REF] a constrained bi-Laplacian optimization is applied globally on the entire domain. Bi-Laplacians increase the order of continuity from C 0 to C 1 but loose the maximum principle, i.e. loose the important monotonicity property of the resulting function inside the MS cells. Inequality constraints are thus added to enforce the scalar field to be monotonic inside each cell. In [START_REF] Jacobson | Smooth shape-aware functions with controlled extrema[END_REF] only the position and value of extremas are constrained to compute smooth shapes by minimizing a non-linear functional similar to [START_REF] Weinkauf | Topology-based smoothing of 2D scalar fields with C 1 -continuity[END_REF]. They further provide an explicit control over the topology of a scalar field. As an extension of [START_REF] Jacobson | Smooth shape-aware functions with controlled extrema[END_REF][START_REF] Weinkauf | Topology-based smoothing of 2D scalar fields with C 1 -continuity[END_REF] Günther et al. [START_REF] Günther | Fast and memory-efficient topological denoising of 2D and 3D scalar fields[END_REF] present a method for reconstructing iteratively a smooth scalar field from a set of selected extrema by combining non-linear optimization based on a monotonicity graph with topological simplification. Both methods [START_REF] Günther | Fast and memory-efficient topological denoising of 2D and 3D scalar fields[END_REF][START_REF] Jacobson | Smooth shape-aware functions with controlled extrema[END_REF] however do not constrain the entire MS complex, but use it to satisfy the monotonicity constrains in the optimization.

However, a big challenge of these numerical reconstruction approaches is the stability in the optimization process. This may create additional critical points in the output preventing it from strictly conforming to the input constraints. Additionally, the overall optimization process might be computationally expensive resulting in extensive running times. Tierny et al. [START_REF] Tierny | Optimal general simplification of scalar fields on surfaces[END_REF][START_REF] Tierny | Generalized topological simplification of scalar fields on surfaces[END_REF] propose an interactive combinatorial algorithm to edit a scalar field on a surface with a user-prescribed extrema that does not build on the computation of a topological structure such as MS complex or contour tree.

The major differences to our work are twofold: first, we provide an explicit polynomial representation of the surface patches and the cell boundaries.

This explicit representation has two advantages: it allows for an arbitrarily dense sampling of the resulting function. It further allows for exact evaluation and computation of differential properties of the scalar function for post-processing applications.

Second, our resulting surface exactly conforms to the input MS complex, meaning that no undesired additional critical points occur as it may be the case with numerical optimization methods.

Our method does not depend on the density/size of the discretization of the input function.

A u t h o r p r e p r i n t 3 Background and Overview

In this section, the required theoretical background on MS complexes is given, for scalar fields defined in a planar domain. First, Sect. 3.1 presents a summary of Morse theory, then Sect. 3.2 explains how MS complexes can be simplified. For further background about Morse theory, the reader is refereed to [START_REF] Morse | The Calculus of Variations in the Large[END_REF][START_REF] Smale | On gradient dynamical systems[END_REF].

Morse Theory

Let f W D ! R be a differentiable function, where D is a bounded open subset of R 2 . A point p 2 D is called a non-degenerated critical point of f if rf .p/ D 0 and the determinant of the Hessian of f in p is non zero. f is called a Morse function if its critical points are non degenerated. There are three kinds of critical points depending on the sign of the two eigenvalues of the Hessian matrix. The number of negative eigenvalues is called the index of the critical point. If the index is zero then p is a minimum, if the index is one, then p is a saddle point and if all eigenvalues are negative, then p is a maximum.

A

C 1 curve l W R ! D is called an integral line of f if .8t 2 R/ @ @t l.t/ D rf .l.t//.
It is a path which follows the gradient of the function, so the values of f along this line parameterized by t are strictly increasing. The limit of l when t ! 1 (resp. t ! C1) is called the origin (resp. destination) of l. The origin and the destination of an integral line are critical points of f , and index.destination/ > index.origin/. The union of a critical point fpg with the set of all points of integral lines which share the origin p is called the ascending manifold of p. Analogously, the union of a critical point fqg with the set of all points of integral lines which have the same destination q is the descending manifold of q. The set of all ascending manifolds and the set of all descending manifolds form two distinct partitions of D.

The MS complex is the partition of the domain D using the sets of integral lines sharing the same origin and destination. In other words, each set in this partition is the intersection of an ascending and a descending manifold. It has the structure of a CW-complex, with cells of three types, as illustrated in Fig. 1a: • 0-cells : these cells are critical points of the function, • 1-cells : these cells are specific integral lines called separatrices, • 2-cells : these cells are always bounded by a cycle of four 1-cells. Two of these 1-cells link a minimum with two saddles, and the two other 1-cells link the two saddles with a maximum. There is no critical point in the interior of a 2-cell.

In the present paper, we do not work with differentiable functions, but with discrete scalar fields sampling a differentiable function at the vertices of a triangulation D.

We therefore rely on Forman's discrete Morse Theory [START_REF] Forman | A user's guide to discrete morse theory[END_REF] which generalizes MS complexes using a combinatorial gradient vector field. In this setting, integral lines and cells of the MS complex are union of domain simplices (e.g. vertices, edges and triangles).

A u t h o r p r e p r i n t

Simplification

To explain the simplification of MS complexes, we use the canonical mapping from a MS complex to an undirected graph which maps 0-cells to nodes and 1-cells to edges. The 2-cells then correspond to cycles in the graph. The simplification of the MS complex is obtained by sequentially removing pairs of adjacent nodes in the graph [START_REF] Edelsbrunner | Hierarchical Morse complexes for piecewise linear 2-manifolds[END_REF][START_REF] Gyulassy | Combinatorial construction of Morse-Smale complexes for data analysis and visualization[END_REF]. The edges of least persistence are removed first [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF][START_REF] Robins | Computational topology at multiple resolutions: foundations and applications to fractals and dynamics[END_REF]. Let fp; qg be an edge in the graph such that index.q/ D index.p/ C 1. Removing the two nodes in the graph implies the removal of all edges incident to them, thereby creating a cycle of length strictly larger than four. New edges are then inserted between each node previously adjacent to p with the same index value than q except q and each node previously adjacent to q with the same index value than p except p. Each of these new edges replace three removed edges. These new edge insertions ensure that all new cycles are of length four and connect one minimum, one saddle, one maximum and a second saddle. This process is illustrated in Fig. 1a, b with index.p/ D 1, i.e. when p is a saddle point and q is a maximum. In this case a new edge is inserted between each saddle point that was adjacent to q and distinct from p, and the maximum that was adjacent to p and distinct from q. The previous process is purely combinatorial. It builds a new simplified topology from the topology of the MS complex, but does not produce a new embedding. A trivial embedding is obtained by associating to each new edge the union of the three 1-cells it replaces, and to each new cycle of union of previous 2-cells. But this trivial embedding does not correspond to a MS complex since there are critical points along the new 1-cells, as illustrated in Fig. 1b. Therefore a new function has to be reconstructed from the original function and from the simplified topology such that the topology of the MS complex of this reconstructed function matches the simplified topology, as illustrated in Fig. 1c. The function reconstruction is the subject of the present paper, as explained in the next section.

A u t h o r p r e p r i n t

Overview

Our method takes as input an initial discrete scalar field f W D R 2 ! R, where D is a planar triangulation. The scalar field values f , also called height values, are given at the vertices of D. We also take as input the persistence threshold, so that in the simplified MS complex, the value of least persistence should be larger than the input threshold. In practice the persistence threshold is specified as a percentage of the maximum persistence. 0% mean no simplification and 100% means maximum simplification.

We first use the recent parallel algorithm by Shivashankar et al. [START_REF] Shivashankar | Parallel computation of 2D Morse-Smale complexes[END_REF] to compute the discrete MS complex of f and to simplify its topology, as explained in Sect. 3.2. Our method is also working with output from other algorithms [START_REF] Bremer | A topological hierarchy for functions on triangulated surfaces[END_REF][START_REF] Edelsbrunner | Hierarchical Morse-Smale complexes for piecewise linear 2-manifolds[END_REF].

As explained in Sect. 3.2, the simplification of the MS Complex is a purely combinatorial process which removes critical points and produces a simplified topology connecting the remaining critical points. And the trivial embedding of this simplified topology consists in 1-cells and 2-cells inside which the original scalar field has critical points, as illustrated in Fig. 1b. Therefore we need to reconstruct a new scalar field whose MS complex matches exactly the simplified topology, i.e. it should be strictly monotonic along the simplified 1-cells, and should contain no critical points inside the simplified 2-cells.

In order to reconstruct such a scalar field we proceed in two steps. In the first step we traverse each 1-cell and compute a smooth uniform quadratic spline which is C 1 -continuous, strictly monotonic and approximates the original scalar field values. This is described in Sect. 4. In the second step we reconstruct a scalar field inside each 2-cell, which interpolates the four monotonic uniform quadratic splines along its boundaries and contains no critical points, as explained in Sect. 5. To compute this scalar field, we first parameterize the domain of the 2-cell into the unit square, and then compute a monotonic piecewise cubic triangular Bézier surface defined in the unit square, which interpolates the height values of the uniform quadratic splines along the edges of the square.

As a consequence our reconstructed scalar field is defined by piecewise polynomial functions along the 1-cells and inside the 2-cells.

Monotonic Smoothing of the 1-Cells

The first part of our method consists in reconstructing smooth and strictly monotone curves along the new 1-cells. Before explaining the details we must point out two facts. First, the 1-cells of the simplified MS complex in the domain are known to have a jagged shape prohibiting a smooth surface reconstruction [START_REF] Bremer | A topological hierarchy for functions on triangulated surfaces[END_REF][START_REF] Weinkauf | Topology-based smoothing of 2D scalar fields with C 1 -continuity[END_REF]. Second, it happens often that two 1-cells, having the same maximum as destination but distinct saddles as origin, merged before the maximum. In the other direction, two 1-cells, having the same minimum as origin and distinct saddles as destination, A u t h o r p r e p r i n t may share the same geometry before getting disconnected. Such merging or disconnection occur at special vertices called junction points. Figure 3-right shows two junction points between a minimum and three distinct saddles. These junction points segment the corresponding 1-cells into several pieces. We first compute new scalar field values at the junction points, so that subsequent processing may be done independently on each piece of 1-cells. Then for each piece we compute a piecewise quadratic uniform B-spline that is smooth, strictly monotonic and approximates the input scalar field values.

1. Scalar field value at junction points. If we follow a simplified 1-cell with junction points, the sequence of scalar field values at these junction points is not strictly monotone. Assigning iteratively a new value by linear interpolation of the critical end points as proposed in [START_REF] Weinkauf | Topology-based smoothing of 2D scalar fields with C 1 -continuity[END_REF] may run in conflicts in case a junction point is shared by more than two 1-cells. We therefore propose to compute the junction point values globally by solving the following linear programming problem:

P j Q f .p j / f .p j /j ! min subject to max_lower.p j / Ä Q f .p j / Q f .p j / Ä min_larger.p j /
where f is the initial scalar field, Q f the unknown value and p j the junction points in D. max_lower.p j / (resp. min_larger. p j /) denotes the maximum (resp. minimum) height value of all junctions and critical points adjacent to p j and are required to be lower (resp. larger) than Q f . p j /. The idea is to compute a new function value Q f .p j / as close as possible to f .p j / for all junction points, while satisfying the monotonicity requirements along all 1-cells.

Monotone reconstruction of the scalar field along 1-cells. We now compute

a monotone piecewise polynomial function along each part of the 1-cells by applying a monotone B-spline smoothing [START_REF] He | Monotone B-spline smoothing[END_REF] between two adjacent points (junction or critical) of the MS complex. Let .x i / iD0:::n be the n C 1 vertices of a part of a 1-cell in the domain D between two junction or critical points x 0 and x n . Let .t i / iD0:::n be the associated normalized arc length parameters, so that t 0 D 0 < t 1 < < t n D 1. New scalar field values have already been computed in x 0 and x n , and we have Q f .x 0 / < Q f .x n /. We now compute a quadratic uniform B-spline function g ˛W OE0; 1 ! R that is strictly increasing, and approximates the input scalar field values. Let g ˛.t/ D P kC1 iD0 ˛iN 2 i .t/, where t 2 OE0; 1, ˛i 2 R are the control points, 0 D u 2 D u 1 D u 0 < : : : < u k D u kC1 D u kC2 D 1 the uniform knots and N 2 i the B-spline basis functions of degree 2 defined on the knot sequence fu j g jD 2;:::;kC2 . We set ˛0 D Q f .x 0 / and ˛kC1 D Q f .x n / to ensure endpoint interpolation and compute the remaining coefficients ˛i, i D 1; : : : ; k as follows:

n 1 X iD1 j f .x i / g ˛.t i /j ! min subject to (1)
g 0 ˛.u 0 / D 0 (2)
g 0 ˛.u k / D 0 (3)
g 0 ˛.u i / 0; i D 1; : : : ; k 1 (4)
Condition (1) ensures approximation. ( 2) and (3) ensure zero gradients at the extremities and finally (4) ensures monotonicity. Figure 2 shows two results with the same input data but different numbers k of control-points.

As the reconstructed surface, see next section, will be composed of cubic triangular Bézier patches, we first convert the quadratic B-spline function g įnto quadratic Bézier curves by repeated knot insertion and then apply a degree elevation to obtain k cubic Bézier curves. See [START_REF] Farin | Curves and Surfaces for Computer-Aided Geometric Design: A Practical Guide[END_REF] for more details about Bézier and B-spline curves. Note that the parameter k represents the number of curve segments along the spline curve. The higher k the more flexible is the curve and the better it is able to approximate a point set. In contrary, the smaller k is the smoother is the resulting curve. In our experiments we choose either k D 4 or k D 7, see discussion in Sect. 6.

Smoothing of 1-cells in the domain.

As we aim to reconstruct the simplified MS complex with piecewise polynomial functions, we further smooth the jagged 1cells in the domain by approximating in a least square sense each piece of 1-cell by a quadratic uniform parametric B-spline curve c.t/ D P kC1 iD0 c i N 2 i .t/ R 2 , t 2 OE0; 1 with the same number of segments k and knot vector as for the monotone height function g ˛. Endpoints interpolation (critical points or junction

Monotone Interpolation Inside the 2-Cells

In the previous section we have explained how to compute smooth monotonic piecewise polynomial functions that approximate the input scalar field values on each 1-cell of the simplified MS complex.

In this section we go further and described our method for computing for each 2-cell a piecewise polynomial surface that interpolates its boundary 1-cells and has no interior critical points. The construction is done independently for each 2-cell. We begin by parameterizing the 2-cell onto the unit square OE0; 1 2 (Sect. 5.1). Then a monotonic surface defined in the unit square is computed (Sect. 5.2). Combining the monotone surface and the parameterization, a monotone reconstruction of the 2-cell is computed, which interpolates the 1-cells, has the 0-cells as critical points and no other critical points.

Parametrization of a 2-Cell

A 2-cell is a closed connected subdomain ˝ D R 2 bounded by four pieces of 1cells. We need to parameterize ˝in the unit square before computing the monotonic A u t h o r p r e p r i n t interpolant in next section. In other words, a diffeomorphism ˚W ˝! OE0; 1 2 needs to be defined. Such a diffeomorphism maps functions without critical points on the unit square to functions without critical points on ˝. If F is a function without critical points on the unit square, then Q f D F ı ˚W ˝! R is a function without critical points on ˝. We compute the parameterization ˚using the following procedure.

First, we sample points on @˝by evaluating all Bézier curves bounding ˝at a fixed number c of parameter values. Usually, we set c D 5. Then, a conforming Delaunay triangulation is computed which has all boundary points as constraints. The number of inserted points in the interior of ˝depends on the expected triangle quality. We set the triangle quality parameter so that the size of inner triangles is similar to the average sample point distance on the boundary. We use the CGAL [START_REF]CGAL: Computational Geometry Algorithms Library[END_REF] library to compute the triangulation.

A parameter value .u; v/ 2 OE0; 1 2 is computed for each vertex of the triangulation as follows. First, the boundary curves are mapped to the border of the unit square by associating the critical points (resp. junction points) with the corners of the square: the minimum is mapped to the parameter .0; 0/, the maximum to .1; 1/ and the two saddle points to the remaining corners. In some cases @˝does not contain the minimum or the maximum of the 2-cell because the boundary 1cells merge, as illustrated in the right image of Fig. 3. In these cases, we map to the corner .0; 0/ (resp. .1; 1/) the junction point that is closest to the missing minimum (resp. maximum) along @˝. Between each corner along @˝, we assign the parameter that was used to sample the boundary Bézier curves. Then Floater's Mean Value Coordinates [START_REF] Floater | Mean value coordinates[END_REF] algorithm is applied to compute a parameterization of the triangulation of ˝in the unit-square, constrained by the parameter values on @˝. The parameterization process is illustrated in Fig. 4.

Monotonic Scalar Field on the Unit Square

In Sect. 4 we have computed smooth 1-cells and new scalar field values strictly increasing along these 1-cells. In the previous section we have parameterized each 2-cell in the unit square. This parameterization maps the boundary 1-cells to the edges of the unit square. We now explain how to compute a scalar field on the unit square that interpolates the given increasing scalar field values on the edges of the domain and that has no critical point in the interior of the domain.

We apply techniques from [START_REF] Allemand-Giorgis | Piecewise polynomial monotonic interpolation of 2D gridded data[END_REF], where an algorithm is developed for interpolating a grid of scalar values by a monotonic C1 -continuous piecewise cubic triangular Bézier surface. The interpolant of [START_REF] Allemand-Giorgis | Piecewise polynomial monotonic interpolation of 2D gridded data[END_REF] is uniquely defined by interpolated scalar values and by two partial derivatives at the vertices of a uniform grid. Sufficient conditions on the interpolated values and the partial derivatives are derived to ensure that the interpolant has no interior critical points. More precisely, the interpolated values are supposed to be strictly increasing along the diagonals of the grid. Moreover if no partial derivative is known, then [START_REF] Allemand-Giorgis | Piecewise polynomial monotonic interpolation of 2D gridded data[END_REF] provides an algorithm that computes partial derivative values leading to a monotonic interpolant.

The monotone interpolant [START_REF] Allemand-Giorgis | Piecewise polynomial monotonic interpolation of 2D gridded data[END_REF] is applied as follows. If there is one or more junction points along one of the boundary 1-cells, the other 1-cells are subdivided, 1 so that the same number of cubic Bézier curves is used in the four 1-cells. We initialize a uniform grid that has in each direction a number of edges equals to the number of cubic Bézier curves. Along the boundaries of the grid, we prescribe the end-points of the cubic Bézier curves computed in Sect. 4 and the first derivative of these cubic Bézier curves at their end-points. To compute the scalar values to be interpolated at the interior points of the grid, a Laplace equation with the boundary values as constraints is solved. In almost all cases the resulting values are strictly increasing along the grid diagonals. If not, we interpolate linearly the scalar values along the diagonal. Then, the interior partial derivatives are computed by using the results of [START_REF] Allemand-Giorgis | Piecewise polynomial monotonic interpolation of 2D gridded data[END_REF]. In the end, we get a piecewise cubic triangular Bézier function that has no interior critical points, and that interpolates the strictly increasing 1-cells computed in Sect. 4. Such an interpolant is shown in Fig. 5a. Its reparameterization in the 2-cell domain, as explained in Sect. 5.1, is shown in Fig. 5b.

There is a particular case, called pouch cells which occurs when the two saddle points are merged. To handle this type of cells as well, we adapt our reconstruction process by only changing the parameterization step.

Results

We have implemented our algorithm and tested it on three types of functions. All the experiments were performed on a machine with an Intel Xeon W3520 CPU which provides 8 GB RAM. The MS complexes of the input function and its simplified versions were computed with [START_REF] Shivashankar | Parallel computation of 2D Morse-Smale complexes[END_REF] and preprocessed as explained in Sect. 3.3. We visualize the initial and reconstructed functions with color contour plots (linear color scale) augmented with a regularly sampled set of isolines. The MS complexes are visualized with the ascending (red curves) and descending 1-cells (green curves) together with the critical points, i.e. saddles (green dots), maxima (red dots) and minima (blue dots).

The first example is a hand-made scalar data set that has been created using Shepard's scattered data interpolation method [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF]. Given a set of N randomly chosen position and function values fx i I f i g, the method computes a function F .x/ D P N kD0 f k kx x i k = P N jD0 kx x j k such that F.x i / D f i with zero-gradients at the input points when > 1. Our so-called Shepard data set has then been built by taking randomly 50 data points and values and by evaluating F 2 on a 250 250 regular grid. The resulting function has at least these 50 extrema. A simplification of 25% is applied and the resulting MS complex reconstructed. Figure 6 shows the A u t h o r p r e p r i n t Fig. 6 Applying our algorithm a synthetic example that has been generated using Shepard's scattered data interpolation. We tested our algorithm with two different values of k (number of Bézier curves per 1-cell, grid size of Bézier patches). In both cases our method reconstructs a piecewise smooth scalar field that exactly obeys the topological structure of the MS complex. initial function (a) and its MS complex (b). A 25% simplification is applied (c). The middle and bottom line of Fig. 6 show two different reconstructions of the same simplified MS complex. They differ in the number k of curve segments composing each piece of the 1-cells (see Sect. 4). The role of k is twofold. On one hand it governs the smoothness and flexibility of the curves. The smaller k is, the smoother are the curves; but a higher value gives more flexibility to the curves and thus better approximates the height values. On the other hand, the number of Sibson patches is equal (or sometimes proportional) to k 2 per 2-cell. The middle row shows a reconstruction with k D 4, the bottom row with k D 7. The first step of our algorithm is shown in (d) and (g), where the jagged 1-cells of the simplified MS complex are smoothed by parametric B-spline curves in the plane, and monotone B-spline height functions that best approximate the initial scalar field along the 1-cells. The final reconstructions are shown in (e,f) and (h,i). Even though there is no much difference visible in the final reconstruction contour plots, it can be observed in the 1-cells reconstructions that the curves in the plane have less undulations for k D 4 in (d), but that the boundary 1-cells are less well approximated compared to k D 7 in (g), see the green curves at the bottom left and top right corners.

A u t h o r p r e p r i n t

The model statistics and computation times are listed in Table 1. As expected, reparameterization is the most costfull step. Triangulating the domain of each 2-cell followed by a parameterization takes 95% of the total computation time. It follows further that the influence of k on the computation time is neglectable, even though the number of Sibson patches is proportional to k 2 .

In the second example is a real terrain data set of the Mont Rainier. 2 It consists of 140148 vertices defining a 458 306 grid. This data set is more complex in the sense that the initial function is not smooth and has critical points that are not well distributed over the domain. Figure 7 shows the initial function (a), its MS complex (b) with 1931 critical points. After 2% simplification, the MS complex (c) has 69 critical points. Figure 7d,e,f show the reconstruction with the smoothed 1-cells and final function.

2 http://data.geocomm.com/catalog/. The last data set is sampled on a 100 100 grid from the mathematical function a sin.x=a/ sin.y=a/ with a D 10. The function has 31 critical points in the domain OE 49; 50 2 . Here we did not compute any simplification (0%), but we applied our algorithm to reconstruct a function from the initial MS complex, see Fig. 8.

A u t h o r p r e p r i n t

A u t h o r p r e p r i n t

Conclusion

This paper experienced for the first time the use of geometric modeling techniques to MS complex reconstruction. The main contribution is the application of a shape preserving spline surface. In particular we approximate the jagged 1-cells on the complex with smooth B-splines. We introduce a C 1 monotone Sibson spline interpolant and combine it with a reparameterization to reconstruct the 2-cells individually. The resulting surface obeys completely the given MS complex. The fact to provide an explicit closed form representation of the reconstructed function may be an advantage when follow-up evaluations of the function and any higher order derivatives are required. The independence of the present method of the mesh or grid on which the initial function is defined on may also be advantageous in case where the grid size is very large.

Beside these achievements we noticed several limitations of our method. First, the lack of existing C 1 reparameterisation techniques makes our resulting surface being only piecewise C 0 continuous at the end. Recent techniques in isogeometric analysis, where the problem of parameterizing a domain defined by boundary curves is also a central concern, enable to parameterize an arbitrary domain up to some restriction using C 1 splines [START_REF] Farin | Discrete Coons patches[END_REF][START_REF] Jaxon | Isogeometric analysis on triangulations[END_REF][START_REF] Xu | Parameterization of computational domain in isogeometric analysis: methods and comparison[END_REF]. However, none of these techniques can actually deal with arbitrary domains as they occur in MS complexes. A C 1 parameterization method for arbitrary domains would improve the general smoothness of our results. We could then develop the Sibson interpolant further to reach globally C 1 reconstructions.

A u t h o r p r e p r i n t

Second, regarding computation times, parameterization step seems to be the bottleneck of the presented method. As the reconstruction of the 2-cells is done independently for each 2-cell, the algorithm could be parallelized to improve this part of the method. We believe however, that it would be more efficient to overcome the parameterization step completely and to work directly on the domain. The triangular spline approximation approach in [START_REF] Jaxon | Isogeometric analysis on triangulations[END_REF] could be a further avenue to explore in the future.

Fig. 1

 1 Fig. 1 (a) A scalar function on a rectangular domain with its MS complex. It has six 2-cells and seventeen 1-cells. The lines of steepest ascent (resp. descent) from the saddle points are shown in orange (resp. green) (b) a pair of adjacent critical points in gray is removed. As a consequence a simplified topology is produced with four new 2-cells and three new 1-cells. (c) a function is reconstructed which MS complex matches the simplified topology (see different isoline pattern)

Fig. 2 Fig. 3

 23 Fig. 2 Height values (red dots) are approximated by a monotonic quadratic B-spline function (blue curve). Left: Approximation with k D 4 curve segments. Right: Approximation with k D 12 curve segments

  Figure 3 illustrates this smoothing for a real MS complex. The final monotone smooth height function can be written as .c.t/; g ˛.t//.

Fig. 4

 4 Fig. 4 (a) Boundary of a 2-cell, the red point is a maximum, the blue point is a minimum and the two green points are saddles. (b) triangulation of the 2-cell. (c) parameterization of the 2-cell

Fig. 5

 5 Fig. 5 Colormap of the piecewise cubic scalar field in the unit square. It has no critical points as can be seen from the isolines. (b) the same scalar field after reparameterization in the 2-cell domain

  Fig.[START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] Applying our algorithm a synthetic example that has been generated using Shepard's scattered data interpolation. We tested our algorithm with two different values of k (number of Bézier curves per 1-cell, grid size of Bézier patches). In both cases our method reconstructs a piecewise smooth scalar field that exactly obeys the topological structure of the MS complex. (a) Initial function, (b) initial MS complex, (c) simplified MS complex, (d) smoothed MS complex with k D 4 Bézier curves per 1-cell, (e) reconstruction, (f) reconstruction, (g) smoothed MS complex with k D 7 Bézier curves per 1-cell, (h) reconstruction (k D 4), (i) reconstruction (k D 7)

Fig. 7

 7 Fig. 7 The terrain set of Mt Rainier (a) has 1931 critical points (b). The simplified MS complex with 69 critical points is reconstructed. (a) Initial function, (b) initial MS complex, (c) simplified MS complex, (d) smoothed MS complex, (e) reconstruction, (f) reconstruction

Fig. 8

 8 Fig. 8 Reconstruction of data from sin.x/ sin.y/. The original data set (a) has been reconstructed (d, e) without topological simplification. (a) Initial function, (b) MS complex, (c) smooth MS complex, (d) reconstruction, (e) reconstruction

Table 1

 1 Measured performances model statistics

	Model		MSC		1-cells	2-cells
						time II		#
			simplif.			Delaunay time III Bézier
		#faces #V	level	#CP #2-cells k time I a +param. geometry patches
	Shepard 124,002 62,500 25%	25 22	4 0,13	8,43	1,15	6672
			25%	25 22	7 0,13	8,80	0,39	19,024
	MtRainier 280,296 140,913 2%	69 53	7 0,31	9,64	2,16	43,300
	Trigo	19,602 10,000 0%	31 24	7 0,03	2,11	0,11	6240
	a Time I is the computation time for smoothing and reconstruction of 1-cells

Subdivision of Bézier curves cuts them into several Bézier curves defining exactly the same curve[START_REF] Farin | Curves and Surfaces for Computer-Aided Geometric Design: A Practical Guide[END_REF].

Acknowledgements This research was partially funded by the ERC advanced grant EXPRES-SIVE (no. 291184).