David Defour

Manuel Marin

FuzzyGPU: a fuzzy arithmetic library for GPU

Data are traditionally represented using native format such as integer or floating-point numbers in various flavor. However, some applications rely on more complex representation format. This is the case when uncertainty needs to be apprehended. Fuzzy arithmetic is one of the major tools to address this problem, but the execution time of basic operations such as addition or multiplication makes its usage prohibitive. In this article, thanks to a new representation format and modern GPU characteristics we show that it is possible to greatly reduce the execution time of those operations. These techniques have been implemented in fuzzyGPU, a freely distributed library of common operations over fuzzy number.

I. INTRODUCTION

Management of uncertainty has appeared as a necessity in the design of expert systems, where the information in the knowledge base is ambiguous and imprecise. The information is usually represented in the form of ambiguous sentences such as at least x, probably between 0.4 and 0.7, etc. This particular kind of uncertainty is best modelled using fuzzy numbers, as they allow to deal with such type of quantifiers [START_REF] Zadeh | The role of fuzzy logic in the management of uncertainty in expert systems[END_REF]. The fuzzy approach has been successfully applied in finance [START_REF] Yoshida | A new evaluation of mean value for fuzzy numbers and its application to american put option under uncertainty[END_REF], transportation [START_REF] Bonvicini | Risk analysis of hazardous materials transportation: evaluating uncertainty by means of fuzzy logic[END_REF], supply chain management [START_REF] Chen | A fuzzy approach for supplier evaluation and selection in supply chain management[END_REF], load flow [START_REF] Miranda | Fuzzy modelling of power system optimal load flow[END_REF], [START_REF] Saber | Resource scheduling under uncertainty in a smart grid with renewables and plug-in vehicles[END_REF]... for configuring expert systems to run under uncertainty scenarios. A fuzzy arithmetic based on the α-cut approach can be implemented using interval arithmetic [START_REF] Dubois | Operations on fuzzy numbers[END_REF], however this can be cumbersome in terms of computation intensity, as each fuzzy operation requires at least two interval operations to be computed. Depending on the number of uncertainty levels in the model, the computational burden can increase dramatically and so the calculation time.

GPUs have driven substantial acceleration to numerous regular applications thanks to several mechanism. The first and most powerful one is thread or data level parallelism (TLP) exploited by the numerous compute unit available. Instruction level parallelism (ILP) has also proven to be very useful to hide instruction latency, achieving better performance at lower occupancy [START_REF] Volkov | Better performance at lower occupancy[END_REF]. Another source of acceleration, but less frequently used, come from the hardware accelerated compute unit dedicated to the evaluation of interpolation, special functions or fused operations such as FMA. By going deeper into these subtleties, one can notice that latest CUDA GPUs allow to statically select rounding attribute for every floating-point operation. This characteristic simplifies interval arithmetic operations [START_REF] Collange | Chapter 9 -interval arithmetic in cuda[END_REF], as it suppresses the overhead associated with the changes of rounding mode.

Applications relying on fuzzy logic have already been successfully ported on GPU [START_REF] Anderson | Speedup of fuzzy clustering through stream processing on graphics processing units[END_REF], [START_REF] Martnez-Zarzuela | Fuzzy art neural network parallel computing on the gpu[END_REF]. However, to our knowledge there were no prior work considering fuzzy arithmetic and the underlying representation of data dedicated to a GPU execution. In this article we will investigate how GPU can improve performance of fuzzy arithmetic. To achieve this goal, we propose to replace the traditional lower-upper encoding, which requires to manage a different lower and upper bound for each α-cut, by a midpoint-radius encoding, whenever possible. This modification of data representation greatly reduces the required bandwidth as well as the number of instructions to perform basic operation. These algorithmic modifications combined with hardware specificities of Nvidia's GPU are greatly reducing the cost of fuzzy arithmetic, making fuzzy arithmetic more affordable in term of execution time.

The main contributions brought by this article are as follows:

• Fuzzy arithmetic library. We describe the implementation of a library of fuzzy arithmetic operations, callable from either CPU or GPU code optimized for each type of architecture. It includes two classes of fuzzy numbers, one based on the lower-upper encoding, the other one based on the midpoint-radius encoding. This library is based on C++ overload of basic operation making it easily extendible, portable and efficient. • Benefice of dedicated architecture. GPUs are known as powerful vector coprocessor. We detail in this article, the benefice brought by GPUs and explain how TLP, ILP, static rounding mode and memory usage impact performance at a fine granularity of basic operations on complex format. • Representation format. We propose a new representation format for fuzzy numbers with symmetric membership function which both reduces the number of operations and memory bandwidth compared to the traditional representation format. The rest of the article is divided as follows: in sections II and III we expose the theory underlying our fuzzy arithmetic library. In section IV we present the implementation issues and discuss about trade-off related to ILP, TLP, bandwidth and latency. Then, in section V we present the performance of our library of fuzzy arithmetic. Section VI concludes the paper and introduces some ideas for future work.

II. FUZZY NUMBERS AND THEIR IMPLEMENTATIONS

Fuzzy numbers and fuzzy arithmetic provide an answer where computations have to be done with imprecise parameters such as applications in fuzzy control, decision making, approximate reasoning, optimization statistics with imprecise probabilities [START_REF] Dubois | Operations on fuzzy numbers[END_REF], [START_REF] Zadeh | The role of fuzzy logic in the management of uncertainty in expert systems[END_REF]. Fuzzy arithmetic is often approached through the α-cut concept [START_REF] Dubois | Operations on fuzzy numbers[END_REF], which consists of representing the fuzzy number as a collection of intervals, one for each level of uncertainty in the model, called the α-cuts. Fuzzy operations are then performed by executing interval operations over the α-cuts of the same level. In this section we will recall basic concept of interval arithmetic, as well as fuzzy numbers and arithmetic.

A. Interval arithmetic

Interval arithmetic is a tool that accounts for one single level of uncertainty in numerical computation. It represents each quantity as a range of possibilities, and operates on these ranges, generating new ranges of possibilities for the results. Interval arithmetic is also used to provide certified results for operations subject to rounding errors such as floatingpoint operations. There can be represented using lower-upper encoding, one of its bounds and the diameter of the interval, or by the center of the interval and its radius. Most of the existing implementation such as Boost [START_REF] Brnnimann | The design of the boost interval arithmetic library[END_REF], MPFI [START_REF] Revol | The mpfi library[END_REF] and GAOL [START_REF] Goualard | Gaol 3.1. 1: Not just another interval arithmetic library[END_REF], are based on the lower-upper encoding. Our proposition is based on midpoint-radius representation format. Let describes this two alternatives.

1) Lower-upper encoding: It uses two floating-point numbers, l and u, holding the lower and upper bounds of the interval. The interval [l, u] is the set of elements comprised between l and u, i.e.

[l, u] = {x ∈ R : l ≤ x ≤ u}
The interval operations are defined as

[a, b] + [c, d] = [ (a + b), (c + d)], [a, b] -[c, d] = [ (a -d), (b -c)], [a, b] • [c, d] = [min( (a • c), (b • d), (a • d), (b • c)), max( (a • c), (b • d), (a • d), (b • c))]
where is the rounding attribute toward -∞ and , the rounding attribute toward +∞. The changes of rounding mode guarantee that the resulting interval includes all possible solutions. For many architecture, such as CPU, the rounding mode corresponds to a processor state, which implies flushing the entire pipeline before being able to start an operation with a new rounding mode.

2) Midpoint-radius encoding: Midpoint-radius format holds the midpoint m and the radius ρ [START_REF] Rump | Fast and parallel interval arithmetic[END_REF]. The interval m, ρ is the set of elements whose distance to m is not higher than ρ:

m, ρ = {x ∈ R : |x -m| ≤ ρ}
Using this encoding, the interval operations are defined as

a, α + b, β = (a + b), ( | (a + b)| + α + β) a, α -b, β = (a -b), ( | (a -b)| + α + β) a, α • b, β = (a • b), (η + | (a • b)|+ (|a| + α)β + α|b|)
where is the rounding attribute toward -∞, the rounding attribute toward +∞ and the rounding attribute toward the nearest. is the relative rounding error unit -= 1 2 -, and η is the smallest representable (denormalized) floating point number. In double precision, = 2 -52 and η = 2 -1074 .

B. Fuzzy arithmetic

While interval arithmetic operates on intervals, sets of one single level of uncertainty, fuzzy arithmetic operates on fuzzy numbers, which are sets of several levels of uncertainty. For each level of uncertainty we can obtain a corresponding interval or α-cut.

Every fuzzy number p has a membership function, µ p, which assigns to each value in the real line a degree of membership to the fuzzy number. It exists a unique value m ∈ R for which µ p(m) = 1, the highest degree of membership. The membership function is a monotony around this value, i.e. µ p(x) is increasing for x < m and decreasing for x ≥ m.

A fuzzy number is defined as

p = {(x, µ p(x)) | x ∈ R}
An α-cut of p is obtained by selecting all the elements with a degree of membership higher than a certain α i , i.e.

p αi = {x ∈ R : µ p(x) ≥ α i }
Fuzzy arithmetic operations are carried out by performing the operation on the α-cuts of the same level. That is, if p and q are fuzzy numbers with k levels of uncertainty; p αi and q αi , for i ∈ {0, . . . , k -1} are their k α-cuts; and ⊕ ∈ {+, -, ×}, then r = p ⊕ q is computed through k interval operations,

r αi = p αi ⊕ q αi for i ∈ {0, . . . , k -1}.

III. PROPOSED REPRESENTATION FORMAT OF FUZZY

NUMBER

An important characteristic of fuzzy numbers is the shape of the membership function. Depending on the application, fuzzy numbers may have a symmetric membership function and then all the α-cuts will be centred around the same midpoint [START_REF] Chang | The genetic algorithm based tuning method for symmetric membership functions of fuzzy logic control systems[END_REF], [START_REF] Chen | Rule-base self-generation and simplification for data-driven fuzzy models[END_REF]. When fuzzy numbers exhibit symmetrical shape, we will prove in this section some properties which we will exploit for the proposed implementation. Let's first give the definition of a symmetric fuzzy number: If symmetric fuzzy numbers do not exhibit any interesting property for a representation based on lower-upper interval, this is not the case for midpoint-radius encoding. The proposed Theorem 1 is stating that all the α-cuts defining a symmetrical fuzzy numbers will share the same center.

Definition (Symmetric fuzzy number). Let p be a fuzzy num- ber, µ p its membership function, m ∈ R such that µ p(m) = 1. If µ p is symmetric around m, i.e. µ p(m -x) = µ p(m + x), ∀x ∈ R
Theorem 1. Let p be a fuzzy number, µ p its membership function, m ∈ R such that µ p(m) = 1. Let p αi = m i , ρ i be any arbitrary α-cut, expressed in the midpoint-radius encoding. Then

p is symmetric ⇐⇒ m i = m, ∀p αi .
Proof. We will prove the reciprocal, i.e.

p is non-symmetric ⇐⇒ ∃ p αi , with m i = m. i) p is non-symmetric =⇒ ∃ p αi , with m i = m. If p is non-symmetric, it means that exists x 0 ∈ R such that µ p(m -x 0 ) = µ p(m + x 0 ) (1)
Let

p αi = [l i , u i ] = m i , ρ i be an α-cut, such that l i = m -x 0 . By definition, µ p(l i ) = µ p(u i ) = α i (2)
Now, assume m = m i . Since l i = m -x 0 , then x 0 = ρ i and

u i = m + x 0 .
Replacing values in equation 2, we obtain

µ p(m -x 0 ) = µ p(m + x 0 ) which contradicts equation 1. Hence, m = m i . ii) ∃ p αi , with m i = m =⇒ p is non-symmetric. Let p αi = [l i , u i ] = m i , ρ i be an α-cut with m i = m.
Again, by definition,

µ p(l i ) = µ p(u i ) = α i (3) Say x 0 = m -l i . Since m = m i , m + x 0 = u i . Also, µ p is monotonic, so µ p(m + x 0 ) = µ p(u i ).
Replacing values in equation 3, we obtain

µ p(m -x 0 ) = µ p(m + x 0 )
Hence, p is non-symmetric.

Figure 1 shows an example of symmetric fuzzy number. We see that all the α-cuts have the same common midpoint. Symmetry is an important property for fuzzy numbers as it is preserved by common arithmetic operations such as {+, -, ×} as stated by the proposed theorem 2.

Theorem 2. Let p and q be two fuzzy numbers and ⊕ ∈ {+, -, ×}. If p and q are both symmetric, then r = p ⊕ q is also symmetric.

Proof. Let's say, without loss of generality, that the system has k levels of uncertainty. We can express the k α-cuts of p, q and r in the midpoint-radius encoding as

p αi = p i , β i q αi = q i , γ i r αi = r i , δ i for i ∈ {0, . . . , k -1}.
We have to prove the theorem for the 3 basic operations. We will split the demonstration in two cases, the addition and subtraction first and the multiplication second.

i) ⊕ ∈ {+, -},

r i = (p i ⊕ q i ) δ i = ( |r i | + β i + γ i )
From Theorem 1, since p and q are symmetric, we have

p i = p q i = q Then r i = (p i ⊕ q i ) = (p ⊕ q) = r
for all i ∈ {0, . . . , k -1}. Hence, r is symmetric.

ii) ⊕ = ×,

r i = (p i × q i ) δ i = (η + |r i | + (|p i | + β i )γ i + β i |q i |)
Again, from Theorem 1, since p and q are symmetric, we have

p i = p q i = q Then r i = (p i × q i ) = (p × q) = r
for all i ∈ {0, . . . , k -1}. Hence, r is symmetric.

Implementations of fuzzy arithmetic are available in various programming languages [START_REF] Anile | Implementing fuzzy arithmetic[END_REF], [START_REF] Hanss | The transformation method for the simulation and analysis of systems with uncertain parameters[END_REF], [START_REF] Kolarović | Fuzzy numbers and basic fuzzy arithmetics (+[END_REF]. As for interval arithmetic, they are based on the lower-upper encoding for fuzzy numbers. The implementation we propose in this paper provides dedicated class to manipulate fuzzy numbers in both the lower-upper and the midpoint-radius encoding. The latter has been implemented making use of the result presented in Theorem 2. That is, for practical purposes, a symmetric fuzzy number can be represented as a unique midpoint, which is common to all the α-cuts, plus a set of radius. This allows, when performing fuzzy arithmetical operations, to save bandwidth and execution time as measured in section V.

IV. IMPLEMENTATION DETAILS

Our fuzzy arithmetic library is written in CUDA C and C++, and callable from either CPU and GPU code, transparently adapting itself to the underlying architecture. The source code is available at [START_REF] Defour | Fuzzygpu: a fuzzy arithmetic library for gpu[END_REF].

The implementation consists of successive wrappers. At the user-end, we have the two fuzzy classes, one for each fuzzy encoding. We decided to offer the two encoding, lower-upper representation format for general case and midpoint-radius when the user is dealing with symmetric fuzzy number. The lower-upper fuzzy class relies on a lower-upper interval class, which serves the purpose of holding the α-cuts. The lowerupper interval class, in turn, relies on a rounded arithmetic class for operating between the α-cuts, according to the rules of interval arithmetic. The midpoint-radius fuzzy class is linked directly to the rounded-arithmetic class, as the α-cuts are stored as a unique midpoint, valid for all the α-cuts, plus a set of radius, making use of the results from the previous section. Ultimately, the rounded arithmetic class is a wrapper of C and CUDA compiler intrinsics with dedicated machine instructions.

The fuzzy template classes are parametrized by the number of α-cuts N and the data type T. In the lower-upper encoding, N and T specify the size and type of an array of intervals which hold the lower and upper bounds for each α-cut. In the midpoint-radius encoding, N and T are the size and type of an array of scalar radius, and T is also the type of the scalar midpoint, common to all the α-cuts.

The basic arithmetic operators are overloaded to work on both fuzzy classes. The operations are carried out by a sequential for loop that iterates over the set of α-cuts. We decided not to spread the loop among different threads in GPU code, since most applications in real life do not involve more than three to four degrees of uncertainty in their data model [START_REF] Bonvicini | Risk analysis of hazardous materials transportation: evaluating uncertainty by means of fuzzy logic[END_REF], [START_REF] Chen | A fuzzy approach for supplier evaluation and selection in supply chain management[END_REF], [START_REF] Miranda | Fuzzy modelling of power system optimal load flow[END_REF], [START_REF] Saber | Resource scheduling under uncertainty in a smart grid with renewables and plug-in vehicles[END_REF], [START_REF] Yoshida | A new evaluation of mean value for fuzzy numbers and its application to american put option under uncertainty[END_REF]. This number is way too low to exploit thread level parallelism (TLP). However, TLP may be exploited in vector operations involving fuzzy numbers, through kernels that process each element of a fuzzy array in a different thread.

Performance of basic operations over complex datatype such as fuzzy numbers are impacted by numerous factors. Next, we will discuss and compare the two representation formats regarding the number of instructions, the memory requirement and the ILP.

A. Number of instructions

Table I shows the number of instructions, such as basic operation, minimum, maximum or absolute value, involved in the addition and multiplication for different data types, including fuzzy numbers. In the case of fuzzy data types, the number of α-cuts is represented by x. The lower-upper fuzzy requires only 2 operations less than midpoint-radius on the addition, but requires 2.8 times more operations when it comes to the multiplication. Therefore, by just comparing the number of instructions the midpoint-radius should bring a real speed-up over the lower-upper representation format.

B. Memory usage

Table II shows the memory space required to store different data types. The units have been normalized to the size of one scalar. In the case of fuzzy, again, x represents the number of α-cuts. The lower-upper fuzzy requires double the space than the midpoint-radius. As the bandwidth requirement of the midpoint-radius is half the one of the lower-upper representation format, an application that needs to access memory at a high rate will definitely benefit from that property.

C. Instruction level parallelism (ILP)

Ideal ILP as defined in [START_REF] Goossens | Limits of instruction-level parallelism capture[END_REF], is a good measure of how a given sequence of instructions could be handled on today's but also future architectures. The higher the ideal ILP is, TABLE II: Memory requirements for different data types, where x corresponds to the number of α-cuts.

Data type

Memory usage Scalar 1

Lower-upper interval 2

Midpoint-radius interval 2

Lower-upper fuzzy 2x

Midpoint-radius fuzzy 1 + x the higher the amount of instructions that could be pipelined during the execution of the application. However, a bigger ideal ILP requires a larger amount of hardware resources to exploit it.

Table III shows the ideal ILP of the addition and multiplication for different data types. This ILP is determined as the ratio of the total number of instructions to the number of levels in the dependency tree. The number of α-cuts in fuzzy data types does not affect the denominator, as each α-cut is independent from all others. Lower-upper fuzzy exhibits more ideal ILP than midpoint-radius on both the addition and the multiplication. However, ILP will be exploited differently depending on the GPU generation as well as the data type in use (single or double precision). For example, GPU with CUDA capability 3.0 can schedule up to 2 independent instructions for a given warp scheduler. The impact of ILP on example applications is also studied in the next section.

V. TESTS AND RESULTS

We have described our fuzzy arithmetic library and exposed the main differences between the two fuzzy encoding methods. Now we will present the tests we designed for the purpose of evaluating the performance of our fuzzy library on GPU under different circumstances. GPU performance is often driven by the balance between memory access and computational effort. In this sense, we consider two types of applications: computebound application, where arithmetic instructions are dominant over memory accesses; and memory-bound application.

The following architectures were used in this study:

• Intel Xeon E5645 CPU.

• NVIDIA GeForce GTX 480 GPU (compute capability 2.0).

• NVIDIA GeForce GTX 680 GPU (compute capability 3.0). CUDA version is 4.2 and g++ version is 4.6.3 in all cases.

A. Compute-bound application: AXPY series

Listing 1 shows an AXPY kernel that computes the i-th element of the series c ← ac+b, where all the values are fuzzy numbers, through an iterative loop. All the threads read the same starting value, perform the recursive computation, and then write the result to a different position on an output vector. As we increment the number of iterations n, the number of floating point instructions grows; however, the number of global memory accesses, one for loading the starting value and one for storing the result, remains constant. In this way we can arbitrarily increment the number of floating-point instructions per accesses to global memory . We measured the execution time of our kernel (as well as other performance indicators) under different parameter configurations, including:

• lower-upper and midpoint-radius fuzzy encoding;

• a number of α-cuts from 1 to 24;

• single and double precision. The kernel arguments were tuned for different runs in order to achieve maximum occupancy of the architecture. The results are plotted in figure 2. Specifically, figure 2a shows different performances achieved by different configurations on different CPU and GPU architectures. Performance is represented on a logarithmic scale as the number of iteration performed per second depending on the number of α-cut. The Xeon CPU is running a single threaded version of the kernel using (i) our library, and (ii) the java library in [START_REF] Kolarović | Fuzzy numbers and basic fuzzy arithmetics (+[END_REF] 1 . We observe differences and performance gains coming from heterogeneous sources. First, there is a pure architectural gain on passing from CPU to GPU to execute the code. Second, within the GPU, there is an algorithmic gain on choosing the midpointradius over the lower-upper encoding. Finally, there is a gain on passing from double to single precision in the fuzzy calculations.

The most interesting of these gains is the one driven by the algorithm, as it puts on relief the differences between lowerupper and midpoint-radius fuzzy implementations. Figure 2b shows the speed-up achieved by midpoint-radius over lowerupper using different precision on two GPU architectures. The shapes we observe are the results of differences in the number of operations, memory size and stronger data dependencies in the case of the midpoint-radius. The difference in the number of operations accounts for the main trend and the other factors account for the irregularities.

The main trend is observed between 1 and 8-9 α-cuts. In this range, the curve follows a growing pattern which is consistent with the theoretical analysis of number of operations per different type of fuzzy operations we made before. This analysis is summed up in table I. From this table, we can calculate the ratio of number of operations per cycle of the AXPY loop required by the midpoint-radius encoding, to the number of operations required by the lower-upper encoding, which is

16x 7 + 7x
This function perfectly describes the shape of the speedup curve between 1 and 8-9 α-cuts, according to figure 2b. Moreover, this behaviour does not depend on the architecture neither on the precision. This result validates our theoretical analysis and also serves to illustrate the performance of the library when no limits on the architecture are being disturbed.

After 9-10 α-cuts, we have a complete different situation. Here, the effect of register spilling and local memory accesses with high latency becomes preponderant. Remember CUDA devices of compute capability 2.0 and 3.0 may allocate up to precision values are stored in two consecutive registers. When this limit is attained, the kernel starts using local memory to allocate extra data. Figure 2c shows the number of registers used by the AXPY kernel per implementation. Figure 2d shows the percentage of replayed instructions due to local memory accesses.

We observe that when using the lower-upper encoding and single precision values, the 63 registers limit is reached at about 9 α-cuts. Local memory transactions start immediately after, at 10 α-cuts. At this same point, in figure 2b, we observe an increase of the speed-up in single precision, as the midpoint-radius implementation is not suffering from register spilling yet. The size of the midpoint-radius fuzzy being about the half of the lower-upper, spilling only starts at double the α-cuts, i.e. 18. From this point on, the speed-up curve in single precision moves back to the ideal shape that can be explained purely by the ratio of number of operations.

In double precision the scenario is a little more complex. Both lower-upper and midpoint-radius fuzzy start spilling registers relatively early, between 6 and 8 α-cuts. In figure 2b, we observe that there is a slight drop in the speed-up in double precision, at about 9 α-cuts. The midpoint-radius fuzzy implementation is loosing performance in this zone and we believe this is due to spilling the midpoint of one of the fuzzy numbers involved in the calculation. Note that midpointradius fuzzy arithmetic proceeds by computing the midpoint first and then uses it to calculate all the radius. If one midpoint is spilled to off-chip local memory, there might not be enough independent arithmetic instructions to hide the latency of this memory access. This effect cannot be appreciated in single precision as register spilling does not get too serious even for the higher number of α-cuts considered. This demonstrates that regarding memory, performance is mainly driven by the amount of memory required to store data and secondly by the temporal dependency among them.

As we stated in section IV, real life applications typically involve 3 to 4 α-cuts. In this range, register spilling and local memory accesses with high latency is not an issue. The above results were obtained using the default 16 KB L1 cache configuration. Nevertheless, if pressure over the register file is too high, the user of this library has the possibility to tune the size of L1 cache and shared memory [START_REF][END_REF]. On one hand, increasing the size of the L1 cache allows to use more αcuts before seeing the impact of register spilling; on the other hand, increasing shared memory does not affect the number of α-cuts for which the spilling occurs, but allows to reduce its overall impact on performance.

B. Memory-bound application: sort by keys

Listing 2 shows a thrust [START_REF] Hoberock | Thrust: A parallel template library[END_REF] program that sorts a vector of fuzzy numbers on the device. We observe in this example how easily the fuzzy type is integrated to other GPU libraries, as all the arithmetic operators are overloaded.

The vector is sorted by keys, which are integers of 32 bits. The sorting algorithm used by thrust in this case is radix sort. The thrust kernels will read the fuzzy array from global memory, sort it on the device, and then write the sorted array back to global memory. The sorting process performs one step per key bit, i.e. 32 steps in this case. At each step, each element in the fuzzy array will be read and copied into a new position. This is a typical case of memory-bound application. We measured the GPU sorting time for the lower-upper and midpoint-radius fuzzy encoding, with a number of α-cuts ranging between 1 and 24 and for single and double precision.

The time spent in transferring data between host and device is not included in our measure. Results for the GTX 480 are plotted in figure 3. Specifically, figure 3a shows performance for different fuzzy encoding and precision, in terms of millions of sorted elements per second. Figure 3b presents the same information in terms of speed-up of midpoint-radius over lower-upper.

We observe in figure 3 that for the same number of α-cuts, midpoint-radius encoding allows to sort twice the amount of fuzzy numbers than lower-upper encoding, whether we are in single or double precision. For this application, the size of the data array being sorted becomes the main factor driving performance. Midpoint-radius requires half the size of lowerupper representation, thanks to a shared midpoint between all the α-cuts. The memory ratio of midpoint-radius to lowerupper representation size is given by table II:

2x 1 + x
where x is the number of α-cuts. We observe that the speedup curves in figure 3b follow approximately this ratio except between 8 and 14 α-cuts in double precision and for more than 15 α-cuts in single precision, where sudden performance drops are causing the speed-up to fluctuate. Thrust sort uses shared memory to speed-up the sorting process. When the values being treated by threads within one SM do not longer fit in shared memory, it starts using global memory with a direct impact over performance. As we stated earlier, lowerupper fuzzy needs twice the amount of memory than midpointradius fuzzy. As a result of this, it saturates shared memory at half the number of α-cuts. When both architectures saturate shared memory, the speed-up curve goes back to the normal behaviour, explained by the ratio of memory sizes. Note that single precision midpoint-radius does not saturate shared memory in this experiment and keeps a high performance all along the considered range.

VI. CONCLUSION

In this article, we presented a fuzzy arithmetic library written in CUDA and C++, based on the α-cuts concept and relying on interval and rounded floating point arithmetic optimized for different hardware. Our implementation offers two fuzzy containers. The first one uses a traditional lowerupper encoding for intervals and allows to represent any fuzzy number, regardless of the shape of the membership function. The second one offers a new optimized midpoint-radius encoding for fuzzy numbers having a symmetric membership function, which notably improves computational efficiency. The optimization consists of using a single variable to store all the α-cuts' common midpoint. We also detailed the performance impact of TLP, ILP and memory usage for recent NVIDIA's GPU regarding operations over fuzzy number's for the proposed representation format. We have shown through a theoretical study and real measures for memory and compute bound kernels that a gain of 2 to 20 can be obtained by preferring the midpoint-radius encoding over the traditional lowerupper encoding. However, accuracy of the new representation format has not been evaluated and is part of the future works.

Fig. 1 :

 1 Fig. 1: Symmetric fuzzy number and α-cuts, all centered around the same midpoint.

#

  include "fuzzy_lib.h" template<class T, int N> __global__ void compute_test(int n, fuzzy_lu<T, N> * a, fuzzy_lu<T, N> * b, fuzzy_lu<T, N> * output) { fuzzy_lu<T, N> c = ( * a); for (int i = 0; i < n; i++) c = ( * a) * c + ( * b); output[blockIdx.x * blockDim.x + threadIdx.x] = c; } int main(){ fuzzy_lu<double, 4> * d_a, * d_b, * d_result; cudaMalloc((void ** )&d_result, ...); ... compute_test<double, 4><<<GRID_SIZE, BLOCK_SIZE>>>( n_iters, d_a, d_b, d_result); } Listing 1: AXPY loop kernel (compute-bound test).

  (a) Related performance of the AXPY kernel. (b) Algorithmic gain. (c) Registers used per implementation. (d) Instruction replays due to local memory accesses.

Fig. 2 :

 2 Fig. 2: Compute-bound benchmark.

#Listing 2 :

 2 include "fuzzy_lib.h" #include <thrust/sort.h> int main(){ thrust::device_vector<fuzzy_lu<double, 4> > d_a(N); thrust::device_vector<unsigned int> k(N); ... thrust::sort_by_key(k.begin(), k.end(), d_a.begin()); } Thrust sort by keys (memory-bound test).

  (a) Thrust's sort by keys performance. (b) Algorithmic gain.

Fig. 3 :

 3 Fig. 3: Memory-bound benchmark.

TABLE I :

 I Number of instructions per arithmetical operation, for different data types, where x corresponds to the number of α-cuts.

	Data type	Number of instructions
		Addition Multiplication
	Scalar	1	1
	Lower-upper interval	2	14
	Midpoint-radius interval	4	10
	Lower-upper fuzzy	2x	14x
	Midpoint-radius fuzzy	2 + 2x	5 + 5x

TABLE III :

 III ILP per arithmetical operation, for different data types, where x correspond to the number of α-cuts..

	Data type		ILP
		Addition	Multiplication
	Scalar	1	1
	Lower-upper interval	2	14/3
	Midpoint-radius interval	5/4	2
	Lower-upper fuzzy	2x	4.6x
	Midpoint-radius fuzzy	0.5 + x/2	1 + x

Note that the library in[START_REF] Kolarović | Fuzzy numbers and basic fuzzy arithmetics (+[END_REF] does not provide certified interval arithmetic, with correct rounding attributes, as our library does.


Marin Bertier

Matthieu Perrin

Cédric Tedeschi

On the Complexity of Concurrent Multiset Rewriting

Keywords: Chemical Programming Model, Multiset Rewriting, Complexity Modèle de programmation chimique, réécriture de multi-ensemble, complexité

While the complexity of the digital infrastructures surrounding us seems to increase continuously, autonomic computing appears as a promising paradigm to be explored. The quest for adequate high-level abstractions to program autonomic systems led to the reemergence of rule-based programming as a promising paradigm. Also, nature appears to be a great source of inspiration. The chemical analogy was at the origin of the so-called chemical programming model, a rule-based programming model enhanced with a chemically-inspired execution model, making it a promising candidate to develop autonomic systems. The metaphor envisions a computation as a set of concurrent reactions between molecules of data arising non-deterministically, until no more reactions can take place, in which case, the solution contains the final outcome of the computation. More formally, such models strongly rely on concurrent multiset rewriting: the data are a multiset of molecules, and reactions are the application of a set of conditioned rewrite rules. At run time, these rewritings are applied concurrently, until no rule can be applied anymore (the elements they need do not exist anymore in the multiset). One of the main barriers towards the actual adoption of such models come from their complexity at run time: each computation step may require a complexity in O(n k ) where n denotes the number of elements in the multiset, and k the size of the subset of elements needed to trigger one rule. In this paper, we explore the possibility of improving the complexity of searching elements (or reactants) through a static analysis of the reaction condition. In particular, through the modeling of the problem with graphs, we give a characterisation of this complexity, based on the results of the well-known subgraph isomorphism problem. This allows us then to give an algorithm for searching reactants, that can be easily parallelised, and whose complexity is lower than the basic case, most of the time.

Introduction

The amount of technologies surrounding us seems to grow continuously. The underlying infrastructure supporting such a computing power (e.g., the Internet) is highly large, dynamic and heterogeneous. When coming to build applications on top of such platforms, its characteristics make the practical leveraging of such a computing power, a quite difficult problem.

One answer comes from autonomic computing, which consists in building systems that can "manage themselves in accordance with high-level guidance from humans" [START_REF] Parashar | Autonomic Computing: An Overview[END_REF]. Put differently, humans should ideally only be required to write a set of high-level rules defining the behaviour of the system. Then the system should be able to run indefinitely, adhering to these high-level, technical details-free rules, whatever the conditions on the underlying platform are.

The implementation of an autonomic system offers several challenges that can not be tackled at once. A prerequisite is to distinguish the development of the low-level machinery from the definition of adequate programming abstractions, allowing this high-level human guidance. The quest for such adequate high-level abstractions led to the reemergence of rule-based programming as a promising paradigm, for instance to specify distributed systems in a declarative manner [START_REF] Grumbach | Netlog, a rule-based language for distributed programming[END_REF]. Also, nature appears to be a great source of inspiration. It was early suggested that biological and chemical system were analogies worth to be explored [START_REF] Babaoglu | Design patterns from biology for distributed computing[END_REF][START_REF] Hariri | Handbookf of Bioinspired Algorithms and Applications[END_REF][START_REF] Viroli | A Biochemical Approach to Adaptive Service Ecosystems[END_REF]. Note that the autonomic computing paradigm was initially proposed by analogy with the autonomic nervous system. The chemical analogy was at the origin of the so-called chemical programming model, a rule-based programming model enhanced with a chemically-inspired execution model, which shows good properties to develop autonomic systems [START_REF] Banâtre | Chemical specification of autonomic systems[END_REF][START_REF] Mostéfaoui | Towards a computing model for open distributed systems[END_REF].

Metaphorically speaking, a chemical program is envisioned as a chemical solution where molecules of data float and react according to some reaction rules specifying the program, to produce new data (products of reactions). Formally speaking, these artificial chemistries [START_REF] Dittrich | Artificial chemistries -a Review[END_REF] strongly rely on concurrent multiset rewriting : the solution is a multiset whose elements are the molecules, and reactions are conditioned rewrite rules to be applied on it. At run time, rules are applied concurrently on the multiset. Once no reactions can be applied anymore, i.e., when no subset of elements satisfying any of the reaction rules' condition can be found in the multiset, the program is said to be inert. In this state, the solution is stable and contains the final result of the program.

While such a model is today regarded as a promising way to specify autonomic systems, one of the main barrier towards its actual adoption is related to its execution complexity : each computation step (i.e. the application of one rewrite rule, assumes that some reactants satisfaying the rule's condition has been located in the multiset. Let us assume the number of objects in the multiset is denoted by n, and that the arity of the rule (i.e., the number of reactants needed to be found for its application, is denoted by k. Then, in the worst case, an exhaustive exploration of all possible combinations of k molecules among n is needed, and the complexity involved is in O(n k ), which can become very problematic, especially with rules with a high arity.

Contribution. In this paper, we explore the possibility of improving the complexity of searching reactants through a static analysis of the reaction condition. In particular, through the modeling of the problem with graphs, we give a characterisation of this complexity, based on the subgraph isomorphism problem. This allows us then to give an algorithm for searching reactants, whose complexity is lower than the basic case of O(n k ) in the sequential case, most of the time. Also, we show that this algorithm can be easily parallelised.

Organisation of the paper. The article is organised as follows. Section 2 introduces the problem and gives its modeling. Also, its NP-completeness is recalled. In Section 3, the characterisation of the complexity, by analogy with the subgraph isomorphism problem, is given. Then, the M2JA algorithm putting this result into practice is described. Its complexity, both in sequential and parallel settings is discussed. Section 5 draws some conclusion and gives some hints for future works.

Preliminaries

The problem to be solved is the search for elements in a multiset, that match a rule's condition. More specifically, the algorithm to be designed takes two input parameters, namely, 1. a chemical rule R = replace x 1 :: T 1 , . . . , x k :: T k by P (x 1 , . . . , x k ) if C(x 1 , . . . , x k ) ; 2. a multiset M composed of n molecules and returns:

• ⊥ if M ∪ R is inert, i.e. if no subset of molecules of M can satisfy R • a tuple (m 1 , . . . , m k ) of molecules in M , where ∀i, m i is of type T i and C(m 1 , . . . , m k ) is true, otherwise.

Study of the rules

If rules were not conditioned, the problem would be simple, since in this case, we would only need to compare the number of available molecules for each type to the number of molecules required. Thus, the remainder of the article will mostly focus on the study of the reaction condition. This reaction condition is a formula of the propositional logic, whose literals are the application of a boolean function to the variables x 1 to x k . The definition of these literals is left very open, as some implementations of the language allow Java expressions, which makes the reactants searching problem undecidable in theory. As evaluation of Java code is far from the subject of this paper, we assume the evaluation of reaction conditions finishes.

Inria

As any propositional formula, a reaction condition can be put in disjunctive normal form. Note that the type of a molecule can be seen as a condition on this molecule. Our formula is then of the following form, where the X i,j are subsets of variables:

R =   replace x 1 , . . . , x k by P (x 1 , . . . , x k ) if L i=1 li j=1 f i,j (X i,j )   , ∀i, j, X i,j ⊂ {x 1 , . . . , x k }. Molecules m 1 , . . . , m k verify C 1 (x 1 , . . . , x k ) ∨ C 2 (x 1 , . . . , x k ) if
and only if they verify either C 1 (x 1 , . . . , x k ) or C 2 (x 1 , . . . , x k ), so the various terms of the disjunction can be searched separately. The rule R can be applied if and only if one of the rules of the following set can be applied:

L i=1    replace x 1 , . . . , x k by P (x 1 , . . . , x k ) if li j=1 f i,j (X i,j )  
 Thus, from now on, we will restrict our analysis to the rules of the following form:

R = replace x 1 , . . . , x k by P (x 1 , . . . , x k ) if f 1 (X 1 ) ∧ • • • ∧ f l (X l ).
For such a rule, let us define V(R) = {x 1 , . . . , x k } the set of its variables and P(R) = l i=1 {(f i , X i )} the set of its predicates, each predicate being associated to a litteral in the reaction condition. For a predicate p = (f, X), we also define F(p) = f and arg(p) = X.

While the number L of such formulae/rules thus generated is potentially exponential over k, it only depends on the initial rule, and not on the size of the multiset. Moreover, as the inertia is detected if and only if no reactants can be found for all of the rules, the search process for different rules is independent, they can be fully parallelised.

In the following definition, we make possible to group predicates. We can then assume, modulo equivalence, that the arguments of the predicates of a rule are pairwise disjoint.

Definition 1 (equivalence of rules) Let us define ≡ as the smallest equivalence relation on the rules such that, for any rule R 1 and R 2 :

V(R 1 ) = V(R 2 ) ∧ P(R 1 ) = P ∪ {(f 1 , E 1 ), (f 2 , E 2 )} ∧ P(R 2 ) = P ∪ {(f 1 ∧ f 2 , E 1 ∪ E 2 )}    =⇒ R 1 ≡ R 2 .
Definition 2 (rank of a rule) We define the rank of a rule R as the greatest arity of its predicates:

rk(R) = max p∈P(R) | arg(p)|. (1)
A rule with a rank of 2 can be represented as a graph, in which the vertices are the variables and the edges are the predicates. Note that most of the problems encountered in literature, can be solved by rules with a rank of 1 or 2. For rules with a higher rank, some predicates connect more than two variables, making the rule impossible to be pictured with such a graph. Figure 1 illustrates this with two examples, of the following rules of arity 4 and of respective rank 2 and 3. While the second one can still be represented graphically, it cannot be expressed as a graph. replace x, y, z, t by P (x, y, z, t

) if f (z, y) ∧ g(x, z) ∧ h(x, t) ∧ i(y, t) ∧ j(z, t) (2)
replace x, y, z, t by P (x, y, z, t) if f (x, y, z) ∧ g(x, y, t) ∧ h(x, z, t) ∧ i(y, z, t)

Structuration of the multiset

The previous section was focused on defining the world of rules. We now devise the characterisation of the world of molecules, and the relationships between both worlds. The central definition in the following is the axiom. It represents a predicate whose each variable has been associated to an actual molecule.

Definition 3 (axiom) Let p = (f, X) be a predicate. An axiom is a pair (p, m) where m is a function that associates a molecule to each variable of p, such that f (m(x 1 ), . . . , m(x n )) is true. We extend to axioms, the notations previously defined for predicates. Let a = (p, m) be axiom. F(a) = F(p) and arg(a) = arg(p). Also, we note

P(a) = p, a[x] = m(x) for all x ∈ arg(p) and M(a) = {a[x] : x ∈ arg(p)}.
We extend this notations to set of axioms: let A be a set of axiom.

M(A) = a∈A M(a) is the set of the molecules contained in A.
Given a rule R and a set of molecules M , we can define the set of axioms that can be constructed. As one can expect, only some of these axioms are of interest for our problem. In the following, we give a characterisation of a set of axioms in regard to the rule it is based on.

Definition 4 (set of induced axioms) For a solution M and a rule R, we define the set of induced axioms by R in M as the set of the axioms composed of a predicate of R and of molecules of M :

A(R, M ) = {a : P(a) ∈ P(R) ∧ ∀x ∈ arg(a), a[x] ∈ M } (4)
Definition 5 (satisfaction of a variable) Let x ∈ V(R), A be a set of axioms and m ∈ M(A). The molecule m is said to satisfy x in A, denoted m |= A x, if for each predicate of the rule pertaining x, we can find at least one axiom in A in which x is associated to m:

∀p ∈ P(R), x ∈ arg(p) ⇒ (∃a ∈ A, P(a) = p ∧ a[x] = m). (5)
Informally speaking, a molecule m does not satisfy a variable x if there is no subsets of molecules containing m in the solution, that verify a predicate in the rule that contains x. Definition 6 (refined and exclusive set of axioms) A set of axioms A is said to be refined if each of its molecule satisfies at least one variable, i.e. if condition 6 is true. A is said to be exclusive if each of its molecule satisfies at most one variable, i.e. if condition 7 is true.

∀m ∈ M(A), |{x ∈ V(R) : m |= A x}| ≥ 1 (6) ∀m ∈ M(A), |{x ∈ V(R) : m |= A x}| ≤ 1 (7)

Inria

More importantly, we define the validity and maximality of a set of axioms:

Definition 7 (maximal and valid set of axioms) A set of axioms A is said to be valid for the rule R if at least one molecule satisfies each of its variables, i.e. if condition 8 is true. A is said to be maximal if at most one molecule satisfies each of its variables, i.e. if condition 9 is true.

∀x ∈ V(R), |{m ∈ M(A) : m |= A x}| ≥ 1 (8) ∀x ∈ V(R), |{m ∈ M(A) : m |= A x}| ≤ 1 (9)
A possible reaction for a rule R and a multiset M is a one-to-one relation between the variables of R and the molecules of M , such that each molecule satisfies the variable it is associated to. From the above definitions, such a relation corresponds to a subset of the set of axioms induced by R in M is all refined, exclusive, maximal and valid at the same time. This notion of affectation is devised in the following. Note that the exclusivity is preserved by all the transformations we do in the following. Thus, we assume that the set of induced axioms is already exclusive. Note that if it is not, it is possible to reduce to this case by copying the molecules for each of their variable and adding inequality constraints to the rule.

Definition 8 (affectation)

Let A be a set of axioms, x 1 , . . . , x p ∈ V(R) and m 1 , . . . , m p ∈ M(A). An affectation of m 1 , . . . , m p to x 1 , . . . , x p , denoted A[x 1 := m 1 , . . . , x p := m p ], is the union of all the refined subsets of the set {a ∈ A : ∀i, 1. The property that ∀i, ∀m, m |= A ′ x i ⇒ m = m i . This is true since, as A ′ ⊂ B, the molecules in A ′ are also in B, so m

x i ∈ arg(a) ⇒ a[x i ] = m i }.

Proposition 1 (characterisation of the affectation) A

′ = A[x 1 := m 1 , . . . , x c := m c ] is the biggest refined subset of A, in the sense of inclusion, verifying ∀i ≤ c, ∀m, m |= A ′ x i ⇒ m = m i . Proof. Let B = {a ∈ A : ∀i, x i ∈ arg(a) ⇒ a[x i ] = m i }
|= A ′ x i ⇒ m |= B x i ⇒ m = m i . 2. A ′ is refined. Let A 1 and A 2 be two refined subsets of B. Let m ∈ M(A 1 ∪ A 2 ) = M(A 1 ) ∪ M(A 2 ). If m ∈ M(A 1 )
, there is a variable x such that m |= A1 x, so m |= A1∪A2 x, and symmetrically for m ∈ M(A 2 ). Therefore A 1 ∪ A 2 is refined. As the number of subsets of B is finite, A ′ is refined.

3.

The fact that it is the biggest such set. Let A ′′ be a refined subset of A verifying ∀i, m |=

x i ⇒ m = m i . Then A ′′ is a refined subset of B. By definition of the affectation, A ′′ ⊂ A ′ .
In regard to the graph of the rule, the axiom set induced by a rule with a rank of 2 can be modeled by a similar graph of molecules. In this new kind of graphs, the vertices are the molecules of the solution, and an edge, labeled with p, links two molecules m 1 and m 2 if there is an axiom a in the axiom set such that P(a) = p and M(a) = {m 1 , m 2 }. Under this formalism, a valid maximal affectation is a subgraph of the graph of molecules that is isomorphic to the graph of the rule, with respect to the labels of the edges.

This problem is known to be NP-complete, as it contains the detection of a clique. This property gives clues on the intrinsic complexity of the reactants searching problem, as we detail now.

NP-hardness of the reactants searching problem

Under the assumption of the termination of the evaluation of reaction conditions -in a time then necessarily independent of both n and k -the reactants searching problem is in class NP, as doing all tests between molecules non-deterministically can solve it in a polynomial time. We will therefore study the NP-hardness of the problem, which, to our knowledge, has never been formally established before.

Proposition 2 (NP-hardness)

The reactants searching problem is NP-hard, depending on its arity.

Proof.

Suppose we know an algorithm that solves the reactants searching problem in a time polynomial both in k, the arity of the rule and n, the number of molecules in the solution.

Let G = (V, E) be a graph. We want to know if G contains a clique of k nodes.

Consider the instance of our problem with R

= replace x 1 , . . . , x k by P (x 1 , . . . , x k ) if k-1 i=1 k j=i+1 (x i , x j ) ∈ E and M = V .
Under our hypothesis, this chemical program finds a clique in G in polynomial time, which reduces our problem to that of finding a clique in a graph, establishing the property.

The rule used to show Property 2 has a rank of 2, and its graph is a clique. In other words, it can be considered as a complicated rule since its reaction condition has as many literals as there are pairs of variables. We should therefore find a way to characterise the complexity of a rule. This is what we do in the next section, which defines the stature of a rule.

Efficient Reactants Searching

In this section, we first precisely characterise the complexity of the reactants searching problem for a rule, given its stature, which is defined first. Then, we present the M2JA algorithm, which leverages this characterisation to allow a better complexity than the basic O(n k ) case, most of the time. Finally, we show that the M2JA algorithm can be parallelised easily.

Stature of a rule

Definition 9 (juncture of a rule) Considering that the variables of a rule are completely sorted by an order ⋖, we define J ⋖ (R), the juncture of R for the order ⋖, as the set of variables for which several predicates contain other variables smaller than them:

J ⋖ (R) = {x ∈ V(R) : |{p ∈ P(R) : ∃y ⋖ x, {x, y} ⊂ arg(p)}| ≥ 2} . ( 10
)
Definition 10 (stature of a rule) The stature of a rule is the size of its smallest juncture among all possible orders.

S(R) = min ⋖ |J ⋖ (R)|. ( 11
)
A juncture J ⋖ (R) such that |J ⋖ (R)| = S(R) is said to be minimal.
Let us illustrate the two previous definitions. As rules with a rank of 2 represents a very important part of chemical programs in practice, we will discuss the stature of some rules having this rank.

Inria 2 ր ց 1 2 5 1 5 ւ ց ր ց ւ տ ւ ց ւ ց 2 ↓ 3 1 7 4 3 8 6 ց ւ ց ր տ ւ ց ր ց ւ 4 3 6 4 9 ց ր 7
Figure 2: The bridge, the height and the lattice with their minimal juncture.

• The stature of a tree is 0. By definition, each node of a tree has a single parent except the root which is an orphan. Therefore, following the topological order, we find no node in the potential juncture.

• The stature of a cycle is one. On one hand, regardless of the order chosen, the greatest element has necessarily two smaller neighbours: its predecessor and its successor in the cycle, making the stature is at least equal to 1. On the other hand, for an order that follows the cycle, all the other elements have 0 or 1 parent, so the stature is at most one.

• Other examples of graphs (see Figure 2) include the bridge and the eight, with a stature of 1, as well as the lattice, which has a stature of 2.

Let us now compare the stature of a rule to its arity. Again, the case of the rules of rank 2 is interesting to start with. Proposition 3 (growth of S) For all rules R 1 and R 2 , if V(R 1 ) ⊂ V(R 2 ) and P(R 1 ) ⊂ P(R 2 ), then S(R 1 ) ≤ S(R 2 ).

Proof.

The principle of the proof is to start by removing the elements of P(R 2 ) \ P(R 1 ) one by one and then those of V(R 2 ) \ V(R 1 ) to get from R 1 to R 2 . We show that each step decreases the stature.

• We first show by induction the next property (P (n)) for all n ≥ 0: every pair of rules

(R 1 , R 2 ) satisfying V(R 1 ) = V(R 2 ), P(R 1 ) ⊂ P(R 2 ) and |P(R 2 )| -|P(R 1 )| = n verify S(R ′ 1 ) ≤ S(R ′ 2 ). -For n = 0, R 1 = R 2 so S(R 1 ) = S(R 2 ).
-Assume now that P (n) is true for a n ≥ 0.

Let

R 1 and R 2 such that |P(R 2 )| -|P(R 1 )| = n + 1. Let p ∈ P(R 2 ) \ P(R 1 ). Consider the rule R such that V(R) = V(R 1 ) = V(R 2 ) and P(R) = P(R 1 ) ∪ {p}. We have |P(R 2 )| -|P(R)| = n then by induction, S(R) ≤ S(R 2 ).
Let ⋖ be an order on V(R) such that J ⋖ (R) is minimal. For all variable x,

{p ∈ P(R 1 ) : ∃y ⋖ x, {x, y} ⊂ arg(p)} ⊂ {p ∈ P(R) : ∃y ⋖ x, {x, y} ⊂ arg(p)} so J ⋖ (R 1 ) ⊂ J ⋖ (R), and S(R 1 ) ≤ |J ⋖ (R 1 )| ≤ |J ⋖ (R)| = S(R) ≤ S(R 2 )
, which shows P (n + 1).

By induction, for all rules R 1 , R 2 with the same variables and P(R 1 ) ⊂ P(R 2 ), S(R 1 ) ≤ S(R 2 ).

• Let us now remove variables. Let R 1 and R 2 such that P(R 1 ) = P(R 2 ) and V(R 1 ) ⊂ V(R 2 ). Regardless of the order ⋖, for any x ∈ V(R 1 ), (x

∈ J ⋖ (R 1 )) ⇒ (x ∈ J ⋖ (R 2 )), so J ⋖ (R 1 ) ≤ J ⋖ (R 2 ), so S(R 1 ) ≤ S(R 2 ).
• Let us finally conclude with any rules R 1 and R 2 such that V(R 1 ) ⊂ V(R 2 ), P(R 1 ) ⊂ P(R 2 ). Let R be the rule with the variables V(R 2 ) and the predicates P(R 1 ). We have:

S(R 1 ) ≤ S(R) ≤ S(R 2 ).
Theorem 1 (upper bound on the stature of a rule with a rank of 2) Let R be a rule of rank 2. Then the following inequality is verified, with equality if and only if the graph of R is a clique.

S(R) + rk(R) ≤ |V(R)| (12)
Proof.

Let R be a rule whose rank is 2, and ⋖ an order on V(R). The two smallest variables x and y have, each, at most one smaller neighbour. Therefore, they are not part of

J ⋖ (R), so S(R) ≤ |V(R)| -2 = |V(R)| -rk(R).
If the graph of R is a clique, any variable z different from x and y, has at least x and y as smaller neighbours, so z

∈ J ⋖ (R), then |J ⋖ (R)| = |V(R)| -rk(R), so S(R) = |V(R)| -rk(R).
Conversely, if the graph of R is not a clique, then there exists x and y that are not connected. Let R ′ be a rule with the same variables as R and a predicate connecting all pairs of variables, except (x, y). By Proposition 3, S(R) ≤ S(R ′ ). Let z be a variable different from x and y, and ⋖ an order in which the three largest items are z ⋖ x ⋖ y in that order. Then J

⋖ (R ′ ) = k -3, so S(R) ≤ S(R ′ ) ≤ k -3.
Let us remark that the previous property is false for a rule with a rank greater than 2, and that even the property that S(R) ≤ |V(R)| -2 is false. This is because two variables can be connected by multiple predicates. This is the case for the rule (3) whose graph is illustrated on Figure 1. Note that this rule with 4 arguments has a stature equal to 3, which contradicts Theorem 1. This shows that in the case of higher rank rules, it is necessary to relax this theorem. This can be done using the notion of equivalence of rules presented before.

Theorem 2 (upper bound on the stature for a rule with a rank greater than 2) For any rule R, we can find

R ′ ≡ R such that S(R ′ ) + rk(R ′ ) ≤ |V(R ′ )|.
Proof. Let R be a rule and R ′ be the rule with the same variables as R but only one predicate

(C, {x 1 , • • • , x k })
, where C is the reaction condition of the rule. As only one predicate is present, the stature is 0. We have rk

(R ′ ) ≤ |V(R ′ )|, so S(R ′ ) + rk(R ′ ) ≤ |V(R ′ )|.
Remark 1 The rule R ′ defined in the proof is usually not the one that minimises S(R ′ )+rk(R ′ ). It is simply given as an illustration.

Inria

Proof. Suppose there are m c+1 , . . . , m k such that A k is valid, then

A k = A[x 1 := m 1 , . . . , x k := m k ] = A[x 1 := m 1 , . . . , x c := m c ][x c+1 := m c+1 , . . . , x k := m k ] ⊂ A[x 1 := m 1 , . . . , x c := m c ] = A c . As A k is valid, ∀x ∈ J ⋖ (R), |{m ∈ A c : m |= x}| ≥ |{m ∈ A k : m |= x}| ≥ 1, so A c is valid.
Conversely, suppose A c is valid. We choose one molecule per variable, following the order ⋖. When choosing a molecule for x i , three cases can occur, depending on in(x i ), the number of predicates that contain other variables smaller than x i , which is defined as:

in(x i ) = |{p ∈ P(R) : ∃x j ⋖ x i , {x i , x j } ⊂ arg(p)}| 1. in(x i ) = 0: If x i has no smaller neighbour, there is a molecule m i |= Ac x i , as A c is valid.
2. in(x i ) = 1.: Let p be this single predicate. There is an axiom a such that

P(a) = p in A c because A c is valid. We choose m i = a[x i ].
3. in(x i ) ≥ 2: In this case, x i ∈ J ⋖ (R) (see equation 10), so x i is already affected in A c . As this assignment is valid, m i ∈ M(A c ) and is suitable with all the already chosen variables.

A k made of these choices is maximal because a choice is made for each variable, and valid as we showed the existence in all three cases.

The M2JA algorithm

In this section, we present the algorithm M2JA (Maximisation of the Minimal Juncture Affectation) that solves the reactants searching problem. Algorithm 1 shows the algorithm M2JA for a rule R, with an arity of k and a juncture J ⋖ (R) = {x 1 , . . . , x c }.

It takes as argument a set of axioms A of type AxiomSet organised like a graph in rank 2. A gives access to all the molecules sorted by variables. Each molecule gives access to a set of references to the axioms in which it is an argument, sorted by predicate. According to this representation, the methods called on A in Algorithm 1 respond in constant time, apart from clone() whose complexity is linear in |A|. By convention, the indices of an array tab[] vary between 1 and tab.size.

The algorithm is based on Theorem 3, that suggests to test the validity of all the affectations of a juncture to detect inertia. It is composed of a main loop (starting on line 2), executed for every tuple of molecules (m 1 , . . . , m c ) that defines the affectation of J ⋖ (R). The loop is composed of two parts:

1. In lines 3-26, the affectation

A ′ = A[x 1 ← m 1 , . . . , x c ← m c ] is computed. If it is not valid,
the computation is dropped. The affectation is computed by a successive elimination of molecules. In the first loop, the molecules that were not chosen are removed. In the second loop, the set is refined by removing all the molecules that cannot be in any refined subset of A ′ . For each molecules m and predicate p for which m can be an argument, we keep an index on the first unremoved axiom corresponding for p and containing m, m.first(p) initialised to 1, in order to check efficiently if a molecule must be removed. (See line 18.)

Algorithm 1: Reactants searching for R with V(R) = {x 1 , . . . , x k } and J ⋖ (R) = {x 1 , . . . , x c }. 1 Molecule[] : findReactants(AxiomSet A) : /* Chose values for J ⋖ (R) = {x 1 , . . . , x c } */ 2 forall m 1 |= x 1 , . . . , m c |= x c do /* Compute A ′ = A[j 1 ← m 1 , . . . , j c ← m c ] */ 3 AxiomSet A ′ ← A.clone(); /* A ′ ← {a ∈ A : ∀i ≤ c, x i ∈ arg(a) ⇒ a[x i ] = m i } */ 4 for i ← 1 to c do 5 forall Molecule m ∈ A ′ .
28 for i ← 1 to k do // x σ(1) ⋖ • • • ⋖ x σ(k) 29 switch p ∈ P(R) : ∃x σ(j) ⋖ x σ(i) , {x σ(i) , x σ(j) } ⊂ arg(p) do //
34 case {p} : // x σ(i) ⋗ x σ(j) ∈ arg(p) 35 A ′′ [i] ← A ′′ [j].axioms(p)[A ′′ [j].first(p)].get(x i ); break; 36 otherwise // x σ(i) ∈ J ⋖ (R) 37 A ′′ [i] ← m i ; 38 return A ′′ ;
Proposition 4 (Computation of the affectation) In Algorithm 1 after line 26, A ′ is valid and its molecules that are not marked as removed are M(A[x 1 := m 1 , . . . , x c := m c ]).

Proof.

We will first use Proposition 1 to show that A ′ has the same molecules as the affectation. We need to show the three following points.

1. The property ∀i, ∀m, m |= A ′ x i ⇒ m = m i . We never put the remove field of a molecule to false between lines 3 and 26. So if a molecule m = m i satisfies x i in A ′ line 26, it also does line 9. But that is not possible because it would have been removed on line 6.

2. A ′ is refined. Line 26, valid is true, so changed is false, so the condition of line 19 was always false in the last iteration, which means exactly that A ′ is refined.

3. The fact that it is the biggest such set. Let B be a refined subset of A such that ∀i, ∀m, m |= A ′ x i ⇒ m = m i . Suppose there is a molecule in M(B) \ M(A ′ ). Let m be the first molecule from B removed from A ′ in Algorithm 1. It was after line 9 because B verifies the first condition, so it was removed on line 21. But m is at least in the same axioms there as in B, as it is the first molecule of M(B) \ M(A ′ ) removed. Therefore, m can not be in B as B is refined. This is a contradiction, so M(B) ⊂ M(A ′ ).

We will now show that A ′ is valid. If note, at least one variable would be associated to zero molecule in A ′ . A ′ is not modified during the last loop, so the same property is true at the beginning of the last loop. Consequently, the variable valid is set to f alse on line 13, and line 26 can not be reached.

Proposition 5 (Correctness) Algorithm 1 returns a valid maximal affectation if there is one and ⊥ otherwise.

Proof.

If the algorithm does not return ⊥, it ends on line 38. In this case, on line 26, A ′ = A[x 1 := m 1 , . . . , x c := m c ] is valid by Proposition 4, then by Theorem 3, there is a valid maximal affectation, which is the one returned by the algorithm, the three cases in the switch corresponding to the three cases in the proof of the theorem.

Conversely, if the algorithm returns ⊥, let us show by contradiction that there is no maximum valid affectation. Suppose that there is a valid affectation A[x 1 := m 1 , . . . , x k := m k ]. The affectation A[x 1 := m 1 , . . . , x c := m c ] was tested by the algorithm, but no value has been returned, since line 39 was reached. We deduce that A[x 1 := m 1 , . . . , x c := m c ] was invalid, which is contradictory since A[x 1 := m 1 , . . . , x k := m k ] is valid.

Proposition 6 (Complexity) We give two estimations of C(R, A), the complexity in the worst case of Algorithm 1, one more precise and one simpler:

C(R, A) = O   rk(R) × |A| × x∈J⋖(R) |{m |= x}|   (13) C(R, A) = O |M(A)| S(R)+rk(R) (14)
Proof. If the solution is inert, which is the worst case, there are exactly x∈J⋖(R) |{m |= x}| executions of the main loop. We will show that the complexity of one iteration of the main loop is proportional to the size of the axiom set.

• In the first part, where the molecules that do not fit with the choice are removed, each axiom is considered at most once by argument, leading to a complexity rk(R) × |A|.

• The complexity of the refinement part is defined by the complexity of its two inner loops. On line 18, the field m.f irst(p) can only grow, so each axiom is checked only once by molecule in its arguments. On Line 22, an axiom can only be removed once by argument.

• ) , and the complexity of the algorithm can be expressed as O |M(A)| S(R)+rk(R) .

• |{m |= x}| ≤ |M(A)|, so x∈J⋖(R) |{m |= x}| ≤ |M(A)| |J⋖(R)| , • |A| ≤ |M(A)| rk(R) ≤ |M(A)| rk(R) rk(R)! . Finally rk(R) × |A| × x∈J⋖(R) |{m |= x}| ≤ 1 (rk(R)-1)! |M(A)| S(R)+rk(R
If the order ⋖ defines a minimal juncture, and if we choose R so that Theorem 2 is verified, the worst-case complexity of the algorithm is O n S(R)+rk(R) for n molecules in the multiset. By Theorem 2, S(R) + rk(R) ≤ |V(R)|. This establishes that the algorithm proposed has a complexity which is either similar or better than O(n k ). Note also that this gain can be evaluated at compile time (it depends on the rank and the stature of a rule, both obtained by static analysis of the rule, matter of Section 2).

Parallelisation

Algorithm 1 is composed of a main loop over all the possible affectations of a juncture, in which the affectation is computed and, if it is valid, reactants are exhibited.

The iterations of the main loop are totally independent, i.e. the value of one affectation does not influence the value of another one. Hence, they can be done in parallel, which leads to the following corollary about the time complexity of a parallel execution of the algorithm: Corollary 1 The reactants searching problem for a rule R and n molecules can be solved by less than n S(R) tasks in O n rk(R) operations each.

Let us discuss the implementation of a such a parallelisation. Going further, let us notice that wo sub-problems must be addressed in order to make a fully decentralised chemical machine:

Axioms induction: the construction of the set of induced axioms is easy in the case of a sequential algorithm, as it is just necessary to make all the necessary tests between molecules, according to the predicates of the rule. This problem becomes more difficult in a distributed system. On the one hand, it is required to make each test at least once, to make sure that no reaction is forgotten. On the other hand, making a single test more than once might have a high impact on the results, as the time needed to evaluate a predicate only depends on the chemical program itself and is unpredictable. Thus, a good system would evaluate each predicate exactly once.

Inria

Affectations parallelisation: in practice, the set of axioms might not be directly copied on the correct number of nodes, each of them knowing of which affectation it has to take care. The set of axioms might even not be storable on a single machine. It is necessary to distribute the choice of the affectations. For the same reasons as above, each affectation must be treated at least once, and should be treated at most once.

Related work

We have seen that, for a rule of rank two, the reagents searching problem could be reduced to the subgraph isomorphism problem. This problem has been well documented [START_REF] Conte | Thirty years of graph matching[END_REF] since the first algorithm proposed by Ullman in 1976 [START_REF] Ullmann | An algorithm for subgraph isomorphism[END_REF], which was based on backtracking. In [START_REF] Messmer | Subgraph isomorphism in polynomial time[END_REF], a decision tree is built, allowing the search for isomorphic subgraphs in polynomial time, after a pre-treatment in exponential time. This work cannot be applied in our case, as the pre-treatment would be made on the graph of the molecules, which constantly change depending of the reactions arising in the computation. As shown in [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF], it is possible to reformulate the subgraph isomorphism problem as a Constraint Search Problem (CSP). The distributed version of CSP, disCSP [START_REF] Ismel | Synchronous, asynchronous and hybrid algorithms for DisCSP[END_REF], is very close to the reactants searching problem. The disCSP problem is generally treated by an exhaustive exploration of the nodes of the network [START_REF] Doniec | A discsp solving algorithm based on sessions[END_REF][START_REF] Zivan | Message delay and discsp search algorithms[END_REF], possibly with an optimisation thanks to backtracking [START_REF] Zivan | Parallel backtrack search on discsp[END_REF].

Naturally, the search for molecules has only been treated through exhaustive search, in particular in distributed settings, as studied in works like [START_REF] Gladitz | Shared Memory Implementation of the Gamma-Operation[END_REF] in a shared memory architecture, or more recently, in the architectural framework in [START_REF] Obrovac | On the Feasibility of a Distributed Runtime for the Chemical Programming Model[END_REF]. The problem of allocating molecules in concurrent settings have been addressed in [START_REF] Bertier | Adaptive atomic capture of multiple molecules[END_REF].

Conclusion

Concurrent multiset rewriting requires the search for elements satisfying a reaction condition. Finding such reactants efficiently is a difficult algorithmic problem, so much that previous studies assumed that it was impossible to solve it in less than O(n k ) operations, for a solution with n molecules and a rule having k arguments.

In this paper, we showed that, through a static analysis of the reaction conditions, this problem can be solved in a time bounded by n rk(R)+S(R) , with rk(R) + S(R) ≤ k. For rules whose rank is 2, we were able to characterise the case of equality. We then showed an algorithm for this, and that it is well suited for a decentralised chemical machine, as it can be parallelised.

Apart from the already exposed problems inherent to the decentralisation, many issues remain to be solved in order to have a complete decentralised chemical machine.

The first issue concerns compilation. So far we have only presented the effective search of molecules knowing the graph of a rule. We used logics arguments to express that compilation was indeed possible. It would be interesting to find efficient algorithms, for example, to choose an optimal order of the variables in the rule. Choosing an order for which juncture is minimal allows a choice at compile time, while providing guarantees on the complexity of the algorithm. But the algorithm is more flexible and works regardless of the juncture. The complexity can be further reduced by a finer choice of the order, with respect to the quantity of molecules for each variable.

An other issue is the integration of other aspects of HOCL, including higher order, that changes the behaviour of the program over time, sub-solutions, and of course the ability to have multiple rules in the same solution. In this last aspect, it is possible to parallelise the search upon the rules.

Finally, let us mention that we could go further in the analysis of the reaction conditions. For example, many literals are of the form f (x) • g(y) where x and y are variables, f and g are functions to the same of set E and • is an order or equivalence relation on E. In this case, the search can be improved, for example through a sorting on the values of f (m) and g(n) for the suitable molecules m and n.

Figure 1 :

 1 Figure 1: Graphical representations of the rules (2), of arity 4 and rank 2, and (3), of arity 4 and rank 3.

  be the set introduced to define the affectation (Definition 8), and A ′ = A[x 1 := m 1 , . . . , x c := m c ]. Three points must be demonstrated to obtain the property:

Theorem 3 (

 3 maximisability of the juncture's affectation) Let R be a rule of arity k whose variables are ordered by ⋖ with {x 1 , . . . , x c } = J ⋖ (R), and A be a set of axioms. For all m 1 , . . . , m c ∈ M(A), A c = A[x 1 := m 1 , . . . , x c := m c ] is valid if and only if there are m c+1 , . . . , m k ∈ M such that A k = A[x 1 := m 1 , . . . , x n := m k ] is maximal and valid.

  The complexity of the maximisation is only |V(R)| ≪ |A|. The complexity of the main loop does not exceed rk(R)×|A| and is executed x∈J⋖(R) |{m |= x}| times, so the complexity of the whole algorithm is O rk(R) × |A| × x∈J⋖(R) |{m |= x}| . We can simplify the writing:

  molecules(i) s.t. m = m i do

	6	m.removed ← true;
	7	forall Axiom a ∈ m.axioms( * ) do
	8	a.removed ← true;
		/* Refine A ′ : compute a biggest fix point	*/
			/
		/* From this point, we know there is a reaction, let us find it	*/

9 Boolean valid, changed; 10 repeat 11 changed ← false; 12 for i ← 1 to k do 13 valid ← false; 14 forall Molecule m ∈ A ′ .molecules(i) s.t. ¬m.removed do 15 valid ← true; 16 forall Predicate p ∈ x i .predicates() do 17 Axiom at[ ] ← m.axioms(p); 18 while m.first(p) ≤ at.size ∧ at[m.first(p)].removed do m.incrementFirst(p); 19 if m.first(p) > at.size then // m.axioms(p) is empty 20 changed ← true; 21 m.removed ← true; 22 forall Axiom a ∈ m.axioms( * ) do 23 a.removed ← true; 24 if ¬valid then break; 25 until ¬valid ∨ ¬changed; 26 if ¬valid then continue ;// A ′ is not valid: try next affectation (line 2) /* A ′ is valid: there is a maximal valid affectation A ′′ ⊂ A ′ . *27 Molecule A ′′ [k];

RR n°8408

According to Theorem 3, a valid affectation of the jointure can be maximised into a valid maximal affectation. Lines 27-38 follow exactly the three cases of the proof of this theorem.If no affectation of the juncture is valid, the solution is inert and the algorithm returns ⊥ on line 39.RR n°8408

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

return ⊥;// No affectation is valid: there is no solution Inria