Manuel Serrano
email: manuel.serrano@inria.fr

The Computer Scientist Nightmare My Favorite Bug

A multitier bug

Mainly Firefox appeared to be affected by the disease. Other browsers, such as Chrome, Midori, or Opera, were doing well. At that early moment the logical

Hop has recently been used by a small French company, which we will call AIMM in the rest of this paper, to implement a new widget for a web multimedia application. This widget contained two parts: a music selector that lets end users browse a database of artists and music, and an HTML5 music player supporting on-demand conversion from one music format to another.

Figure 1 presents a screenshot of the widget. The music player is on the top, the database browser on the bottom. Clicking an artist name pushes a new panel presenting the songs that artist has produced. Clicking the small black arrow to the left of the name restores the previous panel. Graphical effects improve user experience. Depending on the speed of the platform and the web browser used, a new panel slides from left to right or blends into the previous one.

More than a real application, this was an endeavor, or an evaluation in order to understand how Hop fits in the context of realistic multimedia web application. The development was a priori easy because all the elements needed to implement the widget were provided by the Hop development kit off the shelf. I was in charge of the development that I estimated would take only a handful of days. With the AIMM engineers, we agreed on a common API that Hop could use to access the actual database; then I started to develop the widget. After a couple of days, everything went as we wished. The prototype was operational.

Early tests showed that it was reliable enough to enter the second stage: to be tested on the AIMM server.

The AIMM server was a classical Linux Debian hosted by an x86/64 processor. A usual setting, almost the same as the one I used for developing the widget. After having installed additional Linux packages required by Hop, I installed the development kit and the prototype. Then I started to test the widget. At first, it seemed to be doing perfectly well. Everything seemed to be working as it should, until I realized that, on some browsers, clicking the artist names on the selector produced no effect. The new panel never showed up. Plagued with this erroneous behavior, the application was utterly useless. So began the battle I fought to understand and eliminate that bug. suspicion was that something specific to Firefox prevented it from executing the application correctly. However, this first intuition was contradicted by an obvious observation: the same Firefox was doing well when the web page was served by a server running on my machine, that is, when the application was executed entirely locally. This was particularly shocking because the client-side code served by the two servers, the AIMM server and my local machine, were supposed to be identical. My first investigation drove me to inspect the implementation of the client-side part of the application in order to spot which Firefox specificity was breaking its execution.

Hop is a multitier language. A single formalism is used to program the serverside and the client-side of web applications. The web page visualized on a browser is first elaborated on a server as an abstract syntax tree that is eventually compiled on-the-fly into Html. Both ends of the application use the same Document Object Model to reify Html trees as first class values. In Hop, Html tags are standard functions. Hop extends the official set of Html tags with numerous new widgets that are implemented by composing elementary objects. The artist selector was implemented using one of these objects: the SPAGE widget.

A SPAGE, which stands for sliding page, is an Html container (i.e., a box) augmented with a title and a stack of inner elements. Only the top-most element is visible at a time. A visual effect gives the feeling that new elements slide over the previous ones. The source code using the SPAGE in the AIMM application was similar to:

(<SPAGE> :title "languages" :id "lang" (<DIV> (<SPAGE> :title "functional" :parent "lang" (<UL> (<LI> "Lisp") (<LI> "ML") (<LI> "Haskell")))) (<DIV> (<SPAGE> :title "sequential" :parent "lang" (<UL> (<LI> "Fortran") (<LI> "Pascal") (<LI> "C")))

))

The functions <DIV>, <UL>, and <LI> are the Hop functions for the eponymous Html tags1 . The function <SPAGE> is the constructor for the sliding page. In this example, the first page displays the words functional and sequential. When one is clicked, the corresponding page slides onto the screen and the associated enumeration is displayed. The default graphical configuration of SPAGE is defined by the following CSS rule set: The AIMM widget overrides these rules to use a gray background gradient and white texts, as shown in Figure 1. The effect field lets applications choose the graphical effect they desire when a new element is pushed or popped. The four possible values are "fade", "slide", "auto", or "none". It defaults to "auto", which tells Hop to select the most suitable effect according to the characteristic of the browser. The CSS declaration is used by the client-side function spage-effect which implements the graphical effect on the browser. It is defined as:

(define (spage-effect node width step) (let ((e (node-computed-style-get node :effect))) (cond ((eq? e 'fade) (spage-fade node (if (> step 0) -0.3 0.2))) ((eq? e 'slide) (spage-slide node width step)) ((eq? e 'none) (spage-none node (-step))) ((< (hop-config :browser-speed) 80) (spage-fade node (if (> step 0) -0.3 0.2))) (else (spage-slide node width step)))))

When the effect is "auto", it chooses the best effect according to the browser graphical animation speed. If the browser graphical animation is slow, it fades elements. Otherwise, it slides them. At the time that problem occurred, browsers based on Gecko (Firefox) were not able to animate smoothly the visual effect, while browsers based on Webkit (Google-Chrome, Midori, Safari) could do it. So, the configuration effect: auto of the SPAGE widget yield two different visual effects on the two browser families. This explained why Firefox and Google-Chrome behaved differently: they were simply not using the same code for the animation. A brief additional examination of the client-side code unveiled that the very call to spage-fade was the reason for the problem.

The central piece of the Hop development kit is the web broker. It is a full-fledged web server which embeds the Hop compilers. When a program runs on the broker, i.e., the server-side of the application, it generates the program which is installed on the web browser, i.e., the client-side of the application. This client-side part is compiled on-the-fly into JavaScript, the universal language of web browsers, by the Hop broker. As it turned out, the call to spage-fade was erroneously compiled into JavaScript. Instead of: SpageFade( node, step > 0 ? -0.3 : 0.2 ); it was compiled into:

SpageFade( node, step > 0 ? -0.0 : 0.0 );

These zeroes were of course the reason for the problem: the step of the fading being 0, the element stayed invisible forever. So, the problem was not due to Firefox executing incorrectly the code it received but to Hop delivering this incorrect code! Firefox being innocent, the question remained: what was so special with the compilation of the SpageFade call that drove Hop to generate these wrong zeroes?

2 Read, Compile, Print

Although Hop and JavaScript share many similarities (they are both fully polymorphic functional languages, using dynamic type checking, and automatic memory management), the compilation from Hop to JavaScript is not straightforward. First, JavaScript is not fully lexically scoped. Hence, a closure analysis is needed to compile Hop closures into JavaScript closures. Second, JavaScript is not properly tail-recursive, so static analyses are needed to efficiently compile Hop tail recursion into JavaScript loops. Third, because on the web the client-side code first traverses the network, its size matters. The client-side Hop compiler uses fancy features to generate code as compact as possible. All in all, this makes the implementation of the JavaScript compiler about 17KLOC of Hop code (the system is fully bootstrapped). The reason for the erroneous generation of "0.0" was likely to be located somewhere in these lines.

After thorough investigations I ended up with the conclusion that the Hop reader, that is the server-side function in charge of reading Hop files, was at some point broken. At the beginning of the application, the reader was correct and then, after a while, it was returning 0 for all floating point numbers. The compiler was then probably not broken per se, only the reader was wrong.

The Hop syntax is defined by a regular language. It can then be parsed efficiently by a finite state automaton such as those generated by the Hop form "regular-grammar". The actual implementation of the Hop parser looks like:

+ digit)) (in "eE") (? (in "+-")) (+ digit)))) (the-flonum)) ...)
The first part defines variables (float in the code snippet above) which bind regular expressions used in the rules. These rules bind regular expressions, defined in a infix syntax, to expressions that are executed when a pattern matches. Regular grammars are compiled on the fly into Hop procedures. The implementation of regular grammars is about 3.5KLOC of Hop code for the compiler, seconded by 1.1KLOC of native code for the IO layer.

Compiling the Hop reader produces a function such as:

(lambda (ip) (define (the-flonum::double) (rgc buffer flonum ip)) (define (the-fixnum::long) (rgc buffer fixnum ip)) (define (the-symbol::symbol) ...) ... (define (STATE0 ip) ...) (case (STATE0 ip) ... ((10) (the-flonum)) ...))

The function rgc buffer flonum extracts from the input buffer the corresponding floating point number. It is part of the native runtime system. The Hop native compiler (73KLOC of Hop code) can produce either native code, JVM bytecode, or .NET bytecode. The function rgc buffer flonum thus exists for the three backends. C was the backend used for the AIMM experiment. Its implementation is:

#define rgc buffer flonum( ip ) \ strtod( RGC INPUT PORT BUFFER( ip ), 0 );
As the reader at some point was no longer able to read numbers correctly, either RGC INPUT PORT BUFFER or the libc function strtod had to be wrong! This was very shocking. I could not believe that Hop IO buffers were wrong. The overall size of the Hop implementation is about 325KLOC of Hop code which are all read by the Hop reader. If the buffers were not correctly managed, the bootstrap which reads all the 12,371,080 characters composing the source code could not reasonably succeed! I could not either believe that strtod was wrong. The Gnu libc is too widely used and too many applications depend on it for such an obvious error to occur. However, the Linux version running on the AIMM server being rather old, I spent half an hour googling to check if that particular version of the library was known to have a wrong strtod implementation. It was not.

The usual suspects

At this point of the paper, it is probably useful to make a short point. I was chasing a bug that i) only appeared with Firefox but Firefox was innocent; ii) appeared in the client-side of the application but it was located in the server-side of the application; iii) was due to an error in the server-side code in charge of producing the client-side; iv) was due to server-side code that seemed to first behave correctly and then, at some point of the execution, started to behave badly; v) was not observable when executed on my machine but systematic when executed on the AIMM server. The obvious conclusion of all these observations was that something specific to the AIMM server was responsible for the wrong behavior.

When portability across a single operating system seems broken, the first idea that usually comes to mind is to check if there could be a potential confusion between pointers size. In our case we had an ideal suspect: my personal machine was using 32bit pointers while the AIMM server was using 64bit pointers! This suspicion was even strengthened by my previous investigation. As said before, I was more or less suspecting an IO buffer overflow and a confusion between 32bit pointers and 64bit pointers is very prone to buffer overflows. A quick test conducted on another x86/64 Linux machine of our own showed that Hop and the AIMM widget were running like a charm on that other 64bit pointers computer. The problem was elsewhere.

The Hop broker (75KLOC) is deeply multi-threaded. When started, it spawns several preemptive threads (usually 20) that wait for connections. Each client request is handled in a separate thread. So, several client-side compilations may occur simultaneously. Multi-threaded applications are known to be a nightmare to program and an even worse nightmare to debug, in particular because multithreading may break sequentially correct code. The C standard library is plagued with several non-thread-safe functions. In particular, all those that use static buffers, such as strtok. Would it be possible that strtod was also using a static buffer that made it non thread-safe? An excerpt of the Linux man page is given Figure 2. It tells nothing about multi-threading and it proposes no threadsafe alternative. So, at least with the Gnu libc, it is very likely that strtod is thread safe. Another observation mostly invalidates the theory of a potential multi-threading problem. The bug was far too reproducible! The charm of multi-threading errors is their non-deterministic behavior. Here, the bug always occurred, at the same moment, on the same machine. Not something that could have happened with an error in the implementation of the parallelism.

After all, since the bug was reproducible, it was an ideal candidate for a debugger. I recompiled everything in debug mode and tried to learn more. No success. Firstly, the debugger was not easy to use because the error did not show up until the function had been called several thousand times. Secondly, it merely validated some already known facts: yes, the problem also showed up when the application was running inside the debugger; yes, after a while the-flonum was returning 0. Nothing else to learn. In addition to the debugger, I tried the readeval-print loop Hop can spawn. Same result.

Since the debugger and the interactive loop were mostly useless, I instrumented the application to log all the calls to strtod and I wrote another program for replaying the logs. Surprisingly, this new program ran well, as all the calls to strtod completed correctly. So, if strtod was broken, it was broken in a weird context-dependent way. This experiment taught me another important fact: by observing the log player running correctly, I realized that many Hop instances were running well on the AIMM server. Something specific to the AIMM widget made it run incorrectly on the AIMM server! At a moment, for no particular reason, by luck, or because I was unconsciously examining all the possible directions, the small detail that I had in front of my eyes since the beginning and that I had neglected so far attracted my attention: if strtod was wrong by returning 0.0 how come it returned -0.0 and not just 0.0 for "-0.3"? Where did this "-" come from? After all, it seemed that something was correct in the result since the sign was correct! A further experiment showed that actually the result was only wrong for the decimal part of the numbers. That is "0.2" was read as 0.0 but "1.2" was read as 1.0. After this observation everything went fast and easy.

As stated in the Linux man page (see above), strtod is locale-dependent. In French, the radix character is ",", contrary to English that uses ".". Hence, parsing -0.3 in English must return -0.3 while parsing it in French must return -0.0, because the french parsing of the number stops after parsing the characters and 0, ignoring the rest of the string. My computer and the AIMM server were both located in France but my machine was using an English setting while the AIMM server was using a French setting. On the AIMM server, the Linux global environment variable LOCALE was set to "fr FR". This was the reason for the bug. The only problem yet to be solved was to understand why the behavior of strtod changed during the execution. Why did it start with an English setting and at some point switch to French?

The Hop server starts in a bare minimal core image. When an application is loaded it dynamically loads the libraries it depends on. Many multimedia Hop applications depend on the widely used Gstreamer library to manipulate multimedia material. In the AIMM application, it was used to encode MP3 into OGG on the fly (while Google-Chrome supports MP3, Firefox only supports OGG). So, when the AIMM application started it loaded Gstreamer.

The Gstreamer documentation specifies that an application should first call gst init to initialize the library but it says nothing about the locale. However inspecting the source code, and in particular, the file gst/gst.c, I spotted the following: This clearly meant that Gstreamer allowed itself to change the locale of the application, with the side effect of changing the behavior of strtod, and in that case, breaking the Hop client-side compiler. The fix was straightforward and it took exactly two lines of code. It merely required saving the current locale before calling the gst init function and to restore it afterward as in:

char *locale = setlocale( LC ALL, 0 ); gst init( &argc, &argv ); setlocale( LC ALL, locale );

Concluding remarks

This article is unacademic on purpose. It is written in a first-person narrative mode. it contains no citations, no mathematical formulae of any kind, and no algorithms, but it does contain code! Readers will eventually judge, but I think that in spite of containing no theorem, it demonstrates a flaw in the principles we use to design our languages and to write our programs because the bug that is described in this paper could probably occur in many different contexts and with most programming languages. How many readers of this paper can be sure that none of their programs are affected by a similar problem?

This paper tells a story of a program that apparently broke none of our commonly accepted rules but that still went wrong. The program suffered no memory leak. It was type safe, with respect to the common understanding of what a type should be. It suffered no data race. But still it went wrong! Who is to blame for the error then? This is a legitimate question. Someone has to be blamed! Obviously I have my share of responsibility. It might be argued that the Hop implementation of the-flonum is incorrect and this is the reason for the whole problem. It definitively is and the function should be re-written. My knowledge of the C function strtod was too imprecise and it had yielded to the problem described in the paper. I'm the main culprit. However, it is inconceivable to re-implement all the low level functions used by Hop. It is likely that some other functions rely on a hidden state that might be changed externally, in particular amongst the functions used to implement the network and DNS APIs. However, re-writing these functions would take an unreasonable amount of work and independently of the feasibility of the task it is also worth wondering if this approach makes any sense. The purpose of Hop is to provide an infrastructure to help re-using code written elsewhere by other people, maybe using other formalisms. Re-writing all the core APIs in Hop would then contradict the overall goal of the whole system.

If I'm not the only one to blame, who else? Clearly the C functions strtod and atof have a peculiar design that is partly responsible for the problem. First, relying on external information, the locale, they make programs that use them unsealed. Second, using a weak definition for the external representation of numbers, they are prone to incorrectly parse numbers. Even worse, they offer no means to isolate the code relying on them. The locale used by this function is always a kind of oracle that cannot be controlled by the program.

In my opinion, the Gstreamer library also shares a part of the responsibility. First, it might be argued that no API should ever change the global context in which programs that use them execute. Second, if such a modification cannot be avoided, then it should be stated very explicitly in the documentation. In that particular case, Gstreamer seemed to call setlocale for no particularly good reason (parsing ID3 tags potentially expressed in the locale of the host is a weak motivation). Second, as reported, the documentation of the gst init function does not even mention the change of the locale. This might be the biggest fault.

What tool could have helped to prevent the error or to detect it? Usual debuggers were of no help because, first, the execution was correct although with an unintended behavior (after all, strtod was configured to return -0.0), second, the lines of code that were incorrect were not implemented in the main language but in a hidden implementation language that the system tries hard to make invisible to the end programmer.

It is frequently argued that strong static type checking is a powerful tool that helps detecting errors soon. In this particular case it would have been of no help. With respect to the type system found in general purpose languages, everything in the execution was correct. Maybe this tells us that the types we are used to are inadequate? Maybe this tells that we should rely on types that denote concrete entities and not abstract mathematical things? Maybe the type of strtod should not be a function that takes a string of characters and that returns a real but a function that takes a floating point number represented as a string of characters and that returns a real? Maybe the type of gst init should not only be a function that accepts an integer and an array of strings and that returns an integer but also something that denotes the effect on the current locale?

I think the real reason for the error presented in this paper is that the languages we have designed so far, and maybe even the way we think a program should be written, do not allow programs to compose safely. The error presented in this article came from an apparently innocuous function of a widely used library that was permitted to modify the whole behavior of the application. As long as the systems we design will be liable to such behaviors this kind of error will keep occurring. I think this story shows the limit of the current direction we are following in designing our languages and programs. I also think that it shows that we should explore different paths where safe composition is the starting point of the design. In the era we are to enter where programs need to use heterogeneous sets of resources and data, this will be vital for the reliability of our yet to be invented brave new applications. The expected form of the (initial portion of the) string is optional leading white space as recognized by isspace(3), an optional plus ('+') or minus sign ('-') and then either (i) a decimal number, or (ii) a hexadecimal number, or (iii) an infinity, or (iv) a NAN (not-a-number).

A decimal number consists of a nonempty sequence of decimal digits possibly containing a radix character (decimal point, locale-dependent, usually '.'), optionally followed by a decimal exponent. A decimal exponent consists of an 'E' or 'e', followed by an optional plus or minus sign, followed by a nonempty sequence of decimal digits, and indicates multiplication by a power of 10.

...

CONFORMING TO

C89 describes strtod(), C99 describes the other two functions.

Fig. 2. The Linux man page for strtod

Fig. 1 .

 1 Fig. 1. The Hop web widget

  (regular-grammar ((float (or (: (* digit) "." (+ digit)) (: (+ digit) "." (* digit)))) ...) ... ((: (* digit) (or letter special) (* (or letter special digit (in )))) (the-symbol)) ((: (? (in "-+")) (or float (: (or float (

  , strtold -convert ASCII string to floating-point number SYNOPSIS #include <stdlib.h> double strtod(const char *nptr, char **endptr); ... DESCRIPTION The strtod(), strtof(), and strtold() functions convert the initial portion of the string pointed to by nptr to double, float, and long double representation, respectively.

Note that in Hop, identifiers may contain the special characters "<", ">", "!", or "+". Hence, "<DIV>" is a regular identifier as is "string->integer", "remq!", or even "+".


Julio Toss

João Comba

Bruno Raffin

Motivation

Parallel Voronoi Computation for Physics-based Simulations

Voronoi diagrams are fundamental data structures in computational geometry with applications on different areas like on physics-based simulations. For non-Euclidean distances the Voronoi diagram must be performed over a grid-graph, where the edges encode the required distance information. The major bottleneck in this case is a shortestpath algorithm that must be computed multiple times during the simulation. We present a GPU algorithm for solving the shortest-path problem from multiple sources using a generalized distance function. Our algorithm was designed to leverage the grid-based nature of the underlying graph that represents the deformable objects. Experimental results report speed-ups up to 65x over a current reference sequential method. We discuss the performance benefits of this GPU algorithm for two Voronoi-based interpolation techniques: Sibson and Distance Ratio. Such approaches can be used for soft object simulation algorithms for real time physics engines.

VORONOI DIAGRAMS

The Voronoi diagram is a classical partitioning of the space into closest-point regions. It has a vast application domain, being usually employed for answering proximity queries, for example in classical problems like Finding Nearest Site, Facility Location, Motion Planning and Coverage in sensor networks. The Voronoi diagram is also a key for the "Natural Neighbor Interpolation"(NNI) method proposed by Sibson [START_REF] Sibson | A Brief Description of Natural Neighbor Interpolation[END_REF]. The Sibson's interpolation is a well-known method used for interpolating irregular spaced data with applications in several different fields such as medical imaging, meteorological or geological modeling [START_REF] Beutel | Natural neighbor interpolation based grid DEM construction using a GPU[END_REF], flow map reconstruction [START_REF] Barakat | Adaptive refinement of the flow map using sparse samples[END_REF] and scattered data visualization [START_REF] Park | Discrete sibson interpolation[END_REF]. Natural Neighbor based interpolations have also been applied to the field of solid mechanics by Sukumar et Al. [START_REF] Sukumar | Voronoi cell finite difference method for the diffusion operator on arbitrary unstructured grids[END_REF] in the Natural Element Method (NEM) which uses Sibson and non-sibsonian (Laplace) interpolators to perform crack simulations.

Still in the field of physics simulation, Voronoi-based interpolations have also been applied to meshless simulation of complex deformable

• Julio Toss is a joint PhD student between Universidade Federal do Rio Grande do Sul (UFRGS / Brazil) and Univéristé Grenoble Alpes (UGA / France), e-mail:jtoss@inf.ufrgs.br • João Comba is with UFRGS e-mail: comba@inf.ufrgs.br • Bruno Raffin is with UGA/INRIA -France, e-mail: bruno.raffin@inria.fr bodies [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF]. In that work, the authors compute a Voronoi tessellation on a non-Euclidean space by using a discrete distance map which encodes "material-aware" distances that are biased according to the local rigidity inside the simulated body (see use case on Figure 2). Although the idea of a Voronoi partitioning of the space remains the same as in other classical applications, its computation here is fundamentally different. The distances are not defined on the Euclidean space, instead they rely on shortest paths computation over an implicit grid-graph. Computing this variant of Voronoi diagrams is therefore significantly more costly than the classical discrete Voronoi case. Using the Sibson's Interpolation method on a graph space actually requires to compute a shortest-path tree for each interpolated value queried.

A lot of effort has already been dedicated for improving the performance of Sibson's interpolation, particularly when interpolating on a discrete grid in the Euclidean space (Discrete Sibson Interpolation). A popular approach relies on the parallelization with Graphics Processing Units (GPUs) of the Discrete Voronoi diagram like on [START_REF] Beutel | Natural neighbor interpolation based grid DEM construction using a GPU[END_REF] for Digital Elevation Model (DEM) construction.

Early studies on the parallelization of the discrete Voronoi diagram already exploited the GPU parallel processing capabilities [START_REF] Hoff | Fast computation of generalized Voronoi diagrams using graphics hardware[END_REF]. With the popularization of these architectures and evolution of the programming tools, like CUDA and OpenCL, other algorithms for Voronoi computation were proposed allowing a better utilization of their computing power [START_REF] Weber | Parallel algorithms for approximation of distance maps on parametric surfaces[END_REF][START_REF] Rong | Jump flooding in GPU with applications to Voronoi diagram and distance transform[END_REF].

The graph variant of the Voronoi diagram, referred as Graph Voronoi Diagram [START_REF] Erwig | The graph Voronoi diagram with applications[END_REF], defines a vertex partitioning in a connected graph G(V, E). Given a subset S ⊂ V of source vertices, each vertex v ∈ V is assigned to the partition P i of the source vertex s i ∈ S that has the shortest path distance. Applications of this kind of Voronoi diagram arises in several network problems and social data analysis like, for instance, in community detection algorithms [START_REF] Deritei | Community detection by graph Voronoi diagrams[END_REF]. Sequential algorithms and complexity analysis for building Voronoi diagrams on graphs are discussed in [START_REF] Erwig | The graph Voronoi diagram with applications[END_REF].

Parallel solutions for computing the Graph Voronoi have to deal with the shortest path problem. Dijkstra's algorithm implemented with a priority queue provides an efficient solution to this problem but is inherently sequential with lots o synchronizations. For graphs with grid topology, a first parallel algorithm to compute the Voronoi diagram was presented on [START_REF] Toss | Parallel Shortest Path Algorithm for Voronoi Diagrams with Generalized Distance Functions[END_REF]. The algorithm creates a distance map over a discrete image by computing the shortest paths from each pixel to the closest Voronoi seed (Figure 2(d)). The implementation is done on GPUs and takes advantage of the grid connectivity to leverage parallelism. This implementation has a straightforward application in the Voronoi-based interpolation methods used on [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF] and other Natural Neighbor interpolations used in the SOFA real-time physics-based simulation framework [START_REF] Sofa | Simulaton Open Framework Architecture[END_REF].

In the following sections we will focus on the Discrete Voronoi diagram using graph distances (geodesic). We will show how this structure is used in the context of shape function computation for realtime physics-based simulations. We then present our parallel Voronoi implementation [START_REF] Toss | Parallel Shortest Path Algorithm for Voronoi Diagrams with Generalized Distance Functions[END_REF] applied to the Natural Neighbor Interpolation method.

DISCRETE VORONOI DIAGRAMS IN SOFA

The Discrete Voronoi diagram is used in a real-time simulation framework called SOFA [START_REF] Sofa | Simulaton Open Framework Architecture[END_REF]. SOFA is a modular and extendable architecture that allows researchers of different fields related to physics-based simulation to implement and compare their own algorithms.

Deformation using Shape Functions

In [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF], Faure et Al. proposed a method for simulation of complex objects which are composed of mixed types of materials with different stiffness. This method relies on meshless models using sparse samples to capture the displacements at the simulation nodes which are then interpolated within the object using a novel material-aware shape function. This shape function uses a special distance metric that is scaled according to the local material rigidity of the simulated object (Figures

2(b) 2(c))
. This technique allows to easily take into account the material heterogeneity during the simulation. The material properties of the object, like stiffness, is usually represented as a volumetric image of voxels containing the property values. To determine the materialaware distance, the shortest-path is computed over this 3D grid of voxels. Each voxel represents a vertex of a grid-graph with 26 neighbors, and the edges' weights are defined as a function of the stiffness of the adjacent voxels.

To understand the role of shape functions in numerical simulations consider the following steps:

1. The displacement of a deformable object is sampled at discrete locations called degrees of freedom (DoFs). Each DoF has a shape function associated that defines where and how it will influence other points in the object. The area of influence is often referred as the support of the shape function (Figure 3).

2. The goal is then to interpolate these sampled displacements within the rest of the object.

3. The displacement at a given point is interpolated as a weighted sum of the nodes displacements, where the weights are the values of the shape function for each DoFs influencing this point.

Finally, it remains the problem of defining how weights should be computed. This will tell how shape functions from different DoFs would be blended on the rest of the domain. Voronoi diagrams have been used in the so called Voronoi Shape Functions to compute these weights. In the following section we describe different schemes of interpolating nodal values using the Voronoi diagram generated from the simulation nodes.

Voronoi-based Interpolations

One of the most well-known interpolation methods based on the Voronoi diagram is the Natural Neighbor Interpolation (NNI) proposed by Sibson [START_REF] Sibson | A Brief Description of Natural Neighbor Interpolation[END_REF]. Consider the problem of finding neighbors in a set of non uniform distributed data points. By taking the Voronoi tessellation induced by these points, the natural neighbors are the data points whose Voronoi cells share a common frontier. The Sibson interpolation is defined as a ratio of areas in 2D (volumes in 3D). We insert the query point q in the initial Voronoi Tessellation of sample points. The interpolating weight of each data point is given by the ratio between the area stolen from the neighbor Voronoi cell and the area of the newly inserted Voronoi cell (Figure 2(e)).

The Laplace (non-Sibsonian) interpolation ( [START_REF] Belikov | The non-Sibsonian interpolation: A new method of interpolation of the values of a function on an arbitrary set of points[END_REF] as cited in [START_REF] Sukumar | Voronoi cell finite difference method for the diffusion operator on arbitrary unstructured grids[END_REF]) uses the same notion of natural neighbor but it computes the ratio between segments in 2D (areas in 3D). Instead of taking the area of the neighbor cells it uses the ratio between the length of the Voronoi frontier (line in 2D / facet in 3D) and the distance from q its natural neighbor nodes.

Both natural neighbors interpolation methods described above require the computation of a new Voronoi diagram for each query point q added to the input diagram of the data samples. This results in Q executions of the Voronoi diagram where Q is the interpolation resolution desired.

To reduce the complexity in terms of number of Voronoi diagrams computed per query point [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF] propose a alternative interpolation method where the number of Voronoi diagrams computed is a constant The Distance Ratio interpolation introduced by [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF], applies a particular scheme that computes the ratio between the distance from the point to the Voronoi border and the distance to the node (center of the Voronoi cell). Although less formal guaranties on the properties were presented for this interpolant, this algorithm is implemented on SOFA [START_REF] Sofa | Simulaton Open Framework Architecture[END_REF] and shows good practical results, with the advantage of being more computationally efficient.

PARALLEL VORONOI DIAGRAM COMPUTATION

As shown in the previous sections, the Voronoi diagram is an important building block for several other algorithms, in particular on the physics-based simulation domain. In real-time simulations, changes on the topology may modify the computed distances, for examples due to cut on the object. It is then necessary to be able to efficiently recompute the shape functions. Parallel Voronoi computation on euclidean distance have been extensively studied on previous work [START_REF] Hoff | Fast computation of generalized Voronoi diagrams using graphics hardware[END_REF][START_REF] Weber | Parallel algorithms for approximation of distance maps on parametric surfaces[END_REF][START_REF] Rong | Jump flooding in GPU with applications to Voronoi diagram and distance transform[END_REF]. However, such approaches cannot be directly applied on the physics simulation use cases described in this paper, where distance measures are actually shortest path computed on a graph. We show in this section that parallel Voronoi Diagrams on geodesic distances can be computed in parallel using GPUs.

Geodesic distance is closely related to the well studied singlesource shortest-path (SSSP) problem from the graph theory domain. Parallel algorithms for solving the SSSP problem on general graphs have been proposed by [START_REF] Harish | Large Graph Algorithms for Massively Multithreaded Architectures[END_REF][START_REF] Kumar | A modified parallel approach to Single Source Shortest Path Problem for massively dense graphs using CUDA[END_REF][START_REF] Madduri | Parallel shortest path algorithms for solving large-scale instances[END_REF][START_REF] Ortega-Arranz | A new GPU-based approach to the Shortest Path problem[END_REF] . This solutions are usually based either on Dijkstra's or Bellman-Ford algorithm. Using these general graph algorithms for resolving the Voronoi Shape Functions in physics based simulation would imply an over-killing management of a full graph data-structure (like an adjacency matrix or list). Targeting efficiency, we presented in [START_REF] Toss | Parallel Shortest Path Algorithm for Voronoi Diagrams with Generalized Distance Functions[END_REF] a parallel algorithm using GPUs for computing the graph-Voronoi diagram on voxelized 3D grid. This solution is well suited for the physics simulations used in SOFA as it exploits the grid topology that implicitly represent edges. The solution is based on parallel wavefront expansions starting at each Voronoi seed. By using the massive amount of parallel threads of the GPU we are able to compute each Voronoi cell concurrently.

Parallel Voronoi Performance

We present here some experimental results of speed-up of the parallel Voronoi Computation. The experiments where conducted on a NVIDIA GPU GTX480 with 1.5 GBytes of global memory and 15 Multiprocessors with 32 cores each, totaling 480 CUDA cores. The speed-up presented are related to the base sequential version from SOFA [START_REF] Sofa | Simulaton Open Framework Architecture[END_REF] and executed on an Intel CoreTM i7 CPU model 930 with 4 cores running at 2.89Ghz and 12 GB memory .

The benchmark consists in computing the Voronoi diagram on a 3D volume for a given set of randomly distributed data points. The dataset used varies in volume dimensions, distance map distribution and number of seeds. Figure 4 shows an example of two Voronoi diagrams computed with same dimensions and seeds but using different distance distribution. The right diagram was computed considering an image with a gradient of stiffness increasing from right to left whereas the left one has constant stiffness. These two distances distributions are referred as Gradient and Constant on the bar plots of Figure 5. Note that the distance between two neighbor voxels is given by a function of the stiffness between them. As seen in Figure 4 the shape of the Voronoi diagram greatly depends on the stiffness properties of material map of the object simulated .

In the parallel algorithm implemented each CUDA thread processes a single voxel. With larger input volumes more parallelism is exposed, therefore we note a increase of speed-up (Figure 5). At 256 3 the speedup for the constant topology slightly decreases. We attribute this to an overhead of scheduling an excessive number of idle threads. This happens because with a volume of 256 and only 10 seeds the Voronoi diagram becomes overly sparse.

The amount of parallel work available also increases with the number of seed in the Voronoi diagram (Figure 6). When many Voronoi cells are being computed concurrently more threads are active at the same time. In the classical Natural Neighbor Interpolation, each queried point is inserted, one at a time, to the Voronoi diagram of initial data samples.

For each seed added, the initial diagram is updated to generate the new Voronoi cell which will then be used to compute the interpolation.

As seen in the experimental results, the number of seeds computed in the Voronoi diagram has a big impact on the amount parallelism that will be exposed. Computing a single Voronoi cell in parallel reduces the possibilities of parallelization. This situation is even worse when we consider to update an existing Voronoi diagram with the addition of a new seed (we refer to this as the query seed). In this case, only a limited region (around the query seed) would be recomputed (e.g. Figure 2(e)). Performing NNI interpolation simply as sequence of (parallel) Voronoi computations on the GPU doesn't pay-off the overheads of memory transfers and thread scheduling inherent to this architecture.

The evident strategy to generate interpolated values over a grid would be to perform multiple queries concurrently in parallel. This approach was used in [START_REF] Beutel | Natural neighbor interpolation based grid DEM construction using a GPU[END_REF] to generate a interpolation of regular grid using the assumption that every sample has limited radius of influence. This allowed to decompose the domain in independent blocks where queries could be answered in parallel batches. In [START_REF] Park | Discrete sibson interpolation[END_REF], Park et Al. propose a more efficient implementation of Sibson's interpolation on raster images. Their method avoid the explicit construction of a new Voronoi diagram for each query point. Instead they use a Kd-tree structure to find the closest seed to the current query point and use this distance as a radius of influence to increment the interpolation weight. Again, these techniques make assumptions that are valid for euclidean space but not trivially generalized for geodesic (graph) distances.

Parallelization can still be implemented for Natural Neighbor Interpolation over graph spaces if we duplicate some data structures. More precisely, for each NNI query, we copy the input Voronoi diagram of data sample and update it with the addition of a new seed at the coordinates of this query.

Table 1 shows the computation time spent for the parallel Sibson algorithm on a NVidia Tesla K40 GPU. The interpolation is performed over a uniform grid of dimensions 100 × 40 × 20 (80000 voxels) and 20 data points. The GPU algorithm performs batches of parallel NNI queries in sequence iteratively until all the grid is computed. We show the average amount of time spend by each batch of parallel queries and the total time for interpolating the whole grid. Note that the parallel version manages to amortize the overhead when more than 10 NNI queries are done in parallel.

CONCLUSION AND FUTURE WORK

In this work we review a non-usual application of a classical geometrical structure: the Voronoi diagram applied to meshless simulations. We show that the parallelization with GPUs can be used to speed-up the computation of the Voronoi diagram variant required by this domain.

Considering the characteristics of this application -interpolation on a discrete grid and non-euclidean distances -we build upon our parallel Voronoi algorithm to show how it can be used to parallelize the wellknow Natural Neighbor Interpolation method.

The parallel Voronoi implementation presented has a valuable application in soft object simulation methods. It provides a performance solution to the physics-based simulation community that wants to employ Voronoi shape-functions in their meshless simulations.

Finally, recent work have proposed to also use Voronoi Shapefunctions on grids with "extended" connectivity (non-manifold grids [START_REF] Manteaux | Interactive Detailed Cutting of Thin Sheets[END_REF]) which would allow to represent objects with more complex topologies in meshless frameworks. Possible extensions of this work will consider the application of our parallel algorithm in these new domains.

Figure 1 .

 1 Figure 1. Voronoi diagram of a set of data points (in blue) in 2D Euclidean space. The space is partitioned into closest-point regions called Voronoi cells.

Figure 2 .

 2 Figure 2. Use case example: (a) The T-bone steak from the SOFA dataset [6, 16] contains a mixture of flexible meat, softer grease and a rigid bone. As input we take the voxelized material map of stiffness values (b) and the coordinates of the simulation nodes. Distances to the nodes inside the object (c) are biased according to the stiffness values and used to compute the Voronoi diagram (d) of the nodes. This diagram will be used to compute the natural neighbors interpolation (e) on the other voxels of the domain.

Figure 3 .

 3 Figure 3. Shape-functions' weights for the 3 simulation nodes (from Figure 2) computed with two different interpolation methods: Sibson (top) and Distance Ratio (bottom). Weights are normalized starting at 1 at the node locations (Voronoi cell center) and decreases until it vanishes outside of the support.

Figure 4 .

 4 Figure 4. Comparison of Voronoi diagrams generated with the same set of seeds on two different material maps. Left: with an uniform stiffness. Right: with a stiffness gradient

Figure 5 .

 5 Figure 5. Speed-up for different input sizes. Gradient and constant topologies are presented for synthetic benchmarks only. Steak's topology corresponds to the dataset shown in the use case of Figure 2.

Figure 6 .

 6 Figure 6. Average speed-up when increasing the number of seeds of the Voronoi diagram. The instance used has a volume size of 128 3 and a gradient material map.

Table 1 .

 1 Computation time for parallel NNI queries. Avg. Sibson Query corresponds to time spent for querying a whole batch.

	Batch Size	Avg. Sibson Query (ms) Total time ( ms)
	CPU 1 (seq)	1.003	81057.226
	GPU 1	9.745	779657.756
	GPU 10	12.413	99309.125
	GPU 100	29.537	23629.870
	GPU 150	34.807	18587.130

ACKNOWLEDGMENT

We would like to thank project Capes/Cofecub 764/13 for the financial support.