Keywords: Aspect Oriented Programming (AOP), Stochastic Simulation, Multiple Replications In Parallel (MRIP), GPU, CUDA, Warp-Level Parallelism (WLP)

Stochastic simulations involve multiple replications in order to build confidence intervals for their results, and Designs Of Experiments (DOEs) to explore their parameters set. In this paper, we propose Warp-Level Parallelism (WLP), a GPU-enabled solution to compute Multiple Replications In Parallel (MRIP) on GPUs (Graphics Processing Units). GPUs are intrinsically tuned to process efficiently the same operation on several data, which is not suited to parallelize MRIP or DOEs. Our approach proposes to rely on small thread groups, called warps, to perform independent computations such as replications. This approach has proved to be efficient on three classical simulation models, but originally lacked the transparency users might expect. In this work, we enhance WLP using Aspect Oriented Programming (AOP). Our work describes the way to combine CUDA and AOP, and brings forward the techniques available to exploit AOP in a CUDA-enabled development.

INTRODUCTION

Replications are a widespread method to obtain confidence intervals for stochastic simulation results. It consists in running the same stochastic simulation with different random sources and averaging the results. According to the Central Limit Theorem, the average result is approximated in an accurate enough way by a Gaussian Law, for a number of replications greater than 30. Thus, for a number of replications greater than 30, we can obtain a confidence interval with a satisfactory precision. This average result depends on the stochastic variability of the application. Some applications can settle for fewer replications, but are consequently less keen to take advantage of MRIP.

There are many cases where a single simulation can last for a while, so 30 of them run sequentially may represent a very long computation time. Because of this overhead, 30 replications are hardly run in most simulations. Instead, a good practice is often to run 3 replications when debugging, and 10 replications are commonly used to compute a confidence interval. To maintain an acceptable computation time while running 30 or more replications, many scientists proposed to run in parallel these independent simulations. This approach has been named Multiple Replication in Parallel (MRIP) in the nineties [START_REF] Pawlikowski | Distributed stochastic discrete-event simulation in parallel time streams[END_REF]. As its name suggests, its main idea is to run each replication in parallell [START_REF] Hill | Object-oriented pattern for distributed simulation of large scale ecosystems[END_REF]; [START_REF] Pawlikowski | Towards credible and fast quantitative stochastic simulation[END_REF]. In addition, when we explore an experimental plan we have to run different sets of replications, with different factor levels according to the experimental framework [START_REF] Hill | Object-oriented analysis and simulation[END_REF]; [START_REF] Amblard | MDA compliant design of SimExplorer a software tool to handle simulation experimental frameworks[END_REF]. In this paper, we will not consider any constraints that need to be satisfied when implementing MRIP. One of the main barriers that often prevents simulationists to achieve a decent amount of replications is the lack of knowledge in the parallelization techniques. Another common hindrance is the amount of parallel computing facilities available. Our work tackles this problem by introducing a way to harness the computational power of GPUs (Graphics Processing Units) to process MRIP or DOEs faster than on a scalar CPU (Central Processing Unit).

GPUs deliver such an overwhelming power at a low cost that they now play an important role in the High Performance Computing world. However, this kind of devices display major constraints, tied to its intrinsic architecture. Basically, GPUs have been designed to deal with computation intensive applications such as image processing. One of their well-known limits are the slow memory accesses. Indeed, since GPUs are designed to be efficient at computation, they badly cope with applications frequently accessing memory. Except by choosing the right applications, the only thing we can do to overcome this drawback is to wait for the hardware to evolve in such a way. Since the NVIDIA generation codenamed Fermi, GPUs have shown a move in this way by improving cache memories available on the GPU. This leads to better performances for most applications at no development cost, only by replacing the old hardware by the state-of-the-art one (Kepler at the time of writing).

Now, what we can actually think about is the way we program GPUs. Whatever the programming language or architecture chosen to develop an application with, CUDA (Compute Unified Device Architecture) or OpenCL, the underlying paradigm is the same: SIMT (Single Instruction, Multiple Threads). Thus, applications are tuned to exploit the hardware configuration, which is a particular kind of SIMD architecture (Single Instruction, Multiple Data). To obtain speed-ups, parallel applications must be implemented in an SIMD-compliant way. This point reduces the scope of GPU-enabled applications.

The SIMT paradigm automatically groups into 32-wide bundles called warps. Warps are the base unit used to schedule both computation on Arithmetic and Logic Units (ALUs) and mem-ory accesses. Threads within the same warp follow the SIMD pattern, i.e. they are supposed to execute the same operation at a given clock cycle. If they do not, a different execution branch is created and executed sequentially every time a thread needs to compute differently from its neighbours. The latter phenomenon is called branch divergence, and leads to significant performance drops. However, threads contained in different warps do not suffer the same constraint. They are executed independently, since they belong to different warps.

In this paper, we describe Warp-Level Parallelism (WLP), a paradigm to evaluate the approach of using GPUs to compute MRIP, using an independent warp for each replication. In order make it easier to use, we introduce an Aspect Oriented Programming (AOP) declination of WLP using an handcrafted preprocessor. As a matter of fact, this paper also details the different solutions to couple AOP with CUDA. It shows what AOP can bring to CUDA-enabled applications, in terms of software engineering and extensibility. In the further sections, we will:

• Describe a mechanism to run MRIP on GPU;

• Propose an implementation of our approach: WLP;

• Introduce the different approaches to implement AOP and those compatible with CUDA;

• Detail the AOP version of WLP;

• Benchmark WLP with three different simulation models.

GENERAL CONCEPTS OF GPU PROGRAMMING AND ARCHITECTURE

This section does a brief recall of the major concepts introduced by GPU programming and especially by CUDA. It also basically describes how a GPU architecture is organized, since these aspects are directly tied to our approach.

The Single Instruction Multiple Threads (SIMT) paradigm

SIMT is the underlying paradigm of any CUDA application. It is based on the well-known SIMD paradigm. While using SIMD, the same instruction is executed in parallel on multiple computational units, but takes different data flows in input. Instead of viewing SIMT as a simple SIMD variant, one needs to understand that it has been created to simplify applications development on GPU. The main idea is first to allow developers to deal with a unique function, named kernel, which will be run in parallel on the GPU. Second, developers manipulate threads in SIMT, which are a much more common tool nowadays than traditional vectors enabling SIMD parallelization.

In order to handle SIMT more easily, CUDA introduces different bundles of threads. As a matter of fact, threads are grouped into blocks, which size and 3D-geometry are defined by the user. The whole blocks of a kernel form a 2D grid. Each thread will be uniquely identified in the kernel thanks to an identifier computed from a combination of its own coordinates and of its belonging block's. More precisely, in addition to grid and blocks, CUDA devices automatically split threads into fixed-size bundles called warps. Currently, warps contain 32 threads. They are extremely important in the low-level mechanisms ruling execution on a GPU.

As long as NVIDIA has defined both CUDA-enabled GPUs' architecture and the SIMT paradigm, the latter is not only convenient, it also perfectly fits the device. Its sole purpose is to be used on a GPU architecture, which is quite different from other multi-core architectures, especially from CPU ones, as we will see in the next part.

Basic architecture of a GPU

While a CPU possesses few cores, each of them allowing the execution of one thread at a time, a GPU possesses a small number of Streaming Multiprocessors (SM) (for instance an NVIDIA Fermi C2050 has 14 SMs). Each SM embeds an important number of computational units (there are 32 floating point computational units -called Streaming Processor (SP) -on each SM of a Fermi C2050). In theory, the floating-point computation power of a GPU device is equal to the number of SMs multiplied by the number of SPs. Another figure that needs to be considered in the architecture is the number of warp schedulers. The latter are key elements of CUDA performance. In fact, memory accesses are done per warp. However, because of memory latency, the warps-schedulers select the warps that have their data ready to process. Consequently, the more warps can be scheduled, the more the memory latency can be hidden.

While the first generation of NVIDIA GPUs was issued with a single warp-scheduler per SM [START_REF] Lindholm | NVIDIA tesla : A unified graphics and computing architecture[END_REF], later generations, such as Fermi and Kepler, own respectively two and four warp schedulers per SM [START_REF] Wittenbrink | Fermi GF100 GPU architecture[END_REF]. Figure 1 shows a simplified representation of a SM of the Fermi architecture. Warps are first employed when threads need to be scheduled on the SM they have been assigned to. In fact, threads within a warp also achieve memory accesses in parallel, before processing the same instruction on this data. To sum up, threads are bound to each other, and must execute the same instructions according to SIMD machinery. Meanwhile, warps are the smallest unit that run in parallel on the different SMs of a GPU, and are the smallest GPU element that is able to process independent code sections. Indeed, given that different warps either run on different SMs, or on the same but at different clock ticks, they are fully independent to each other.

Blocks dispatching and warps scheduling on NVIDIA GPUs

Now that we have introduced the basic functioning of CUDA-enabled GPUs, let us detail the particular features that will help us to achieve MRIP on such architectures. We will see in this part how our GPU-enabled MRIP implementation relies on the scheduling features provided by NVIDIA CUDA devices. The first generation of NVIDIA GPUs were only able to run a single kernel at a time. Thus, blocks of threads were dispatched through all the available SMs in a more or less logical way: SMs were activated in turn, striding indices four by four. When every SM had been activated, the process started again.

One of the key features of architectures from Fermi and more recent is the ability to run several kernels in parallel on the same device. To do so, the way blocks of threads are dispatched through the device has been redesigned in a new fashion. Now, every block of threads, no matter which kernel it belongs to, is first handled by a top-level scheduler referred to as the GigaThread Engine. It is supposed to dispatch blocks of threads to the Streaming Multiprocessors (SMs). The point is CUDA has always proposed asynchronous kernel calls to developers. Now that devices can run several kernels in parallel, GigaThread needs to take into account any potential upcoming kernel. Consequently, the dispatcher cannot reserve all the SMs to run a first kernel, given that a second one could be launched at any time. When the second kernel appears, some resources will still be available so that they can be assigned to the new kernel blocks. Moreover, GigaThread enables immediate replacement of blocks on an SM when one completes executing. Since context switching has been fastened with recent architectures, blocks of threads can fully take advantage of the hardware device thanks to GigaThread dispatching capabilities. From an external point of view, and since we do not have the real specifications, we have noted that the dispatching of blocks does not seem to be deterministic. NVIDIA uses a specific way to place blocks on SMs: indeed, SMs will not be enabled in order. SMs bearing non-consecutive identifiers will in fact run consecutively ordered blocks.

A WARP MECHANISM TO SPEED UP REPLICA-TIONS

Two problems arise when trying to port replications to GPU threads, considering a replication per thread. First, we generally compute few replications, whereas we have seen that GPUs needed to achieve large amounts of computations to hide their memory latency. Second, replications of stochastic simulations are not renowned for their SIMD-friendly behaviour. Usually, replications fed with different random sources will draw different random numbers at the same point of the execution. If a condition result is based on this draw, divergent execution paths are likely to appear, forcing threads within a same warp to be executed sequentially because of the intrinsic properties of the device.

The idea that we propose in this paper is to take advantage of the previously introduced warp mechanism to enable fast replications of a simulation. Instead of having to deal with Thread-Level Parallelism (TLP) and its constraints mentioned above, we place ourselves at a slightly higher scope to manipulate warps only. Let this paradigm be called Warp-Level Parallelism (WLP), as opposed to TLP. Now running only one replication per warp, it is possible to have each replication to execute different instructions without being faced to the branch divergence problem.

But to successfully enable easy development of simulation replications on GPU using one thread per warp, two mechanisms are needed.

First, it is necessary to restrict each warp to use only one valid thread. By doing so, we ensure not to have divergent paths within a warp. Moreover, we artificially increase the device's occupancy, and consequently, we take advantage of the quick context switching between warps to hide slow memory accesses. Theoretically, we should use the lowest block size maximizing occupancy. For instance, a C2050 board owns 14 SMs, and can schedule at most 8 blocks per SM. In this case, the optimal block size when running 50 replications would be 32 threads per block. This situation is represented in Figure 2, where we can see two warps running their respective first threads only. The 31 remaining threads are disabled, and will stall until the end of the kernel. Unfortunately, the GigaThread scheduler, introduced in the previous section, does not always enable a kernel to run on every available SM. In addition, SMs' memory constraints might compromise this ideal case by reducing the number of available blocks per SM. Second, there has to be an easy solution to get a unique index for each warp. TLP relies essentially on threads identifiers to retrieve or write data back. Thus, WLP needs to propose an equivalent mechanism so that warps can be distinguished to access and compute their own data.

Thanks to the two tools introduced in this section, it is possible to create a kernel where only one thread per warp will be valid, and where it will be easy to make each valid thread compute different instructions, or work on different data depending on the new index.

Although we could not figure out the real behavior of the GigaThread Engine dispatcher, the characteristics noticed in this part are sufficient to evaluate the performance of the dispatching policy. Furthermore, the new scheduling features introduced of NVIDIA GPUs benefit the overall performance of our warp-based approach, given that it highly relies on warp scheduling and block dispatching.

IMPLEMENTATION

Now that we have defined our solution, we will propose an implementation in this section. To do so, we need to focus on two major constraints: first, we should keep a syntax close to C++ and CUDA, so that users are not confused when they use our approach. Second, we need to propose compile-time mechanisms as much as possible. Indeed, since WLP only exploits a restricted amount of the device's processing units, we have tried to avoid any overhead implied by our paradigm.

This paper intends to prove that our approach is up and running. Thus, this section will only introduce a restricted number of keywords used by WLP. As we have seen previously, we first have to be able to identify the different warps, in the same way SIMT does with threads. One way to obtain the warp identifier is to compute it at runtime. Indeed, we know that warps are formed by 32 threads in current architectures NVIDIA (2011b). Thus, knowing the running kernel configuration thanks to CUDA defined data-structures, the warp identifier can be determined using simple operations only, similarly to what [START_REF] Hong | Accelerating CUDA graph algorithms at maximum warp[END_REF] have done. The definition of a warpIdx variable containing the warp's identifier can be written as in Figure 3 Conceptually, this definition is ideal because warpIdx is declared as a 'constant variable', and the warp identifier does not change during a kernel execution. This formula fits with the CUDA way to number threads, which first considers threads' x indices, then y and finally z, within a block. The same organization is applied to blocks numbering [START_REF] Kirk | Programming Massively Parallel Processors[END_REF]. Please note that the warpSize variable is provided by CUDA. This makes our implementation lasting since warp sizes may evolve in future CUDA architectures.

Although this method introduces superfluous computations to figure out the kernel's configuration, we find it easier to understand for developers. Another way to compute the warp's identifier would have been to write CUDA PTX assemblyNVIDIA (2011a). The latter is the Instruction Set Architecture (ISA) currently used by CUDA-enabled GPUs. CUDA enables developers to insert inlined PTX assembly into CUDA high-level code, as explained in NVIDIA (2011b). However, this method is far less readable than ours, and would not be more efficient since we only compute warpIdx once: at initialization. This warp identifier will serve as a base in WLP. When classical CUDA parallelism makes a heavy use of the runtime-computed global thread identifier, WLP proposes warpIdx as an equivalent.

Now that we are able to figure out threads' parent warp, let us restrain the execution of the kernel to a warp scope. Given that we need to determine whether or not the current thread is the first within its belonging warp, we will be faced to problems similar to those encountered when trying to determine the warp identifier. In fact, a straightforward solution relying on our knowledge of the architecture quickly appears. It consists in determining the global thread identifier within the block to ensure it is a multiple of the current warp size. Once again, the kernel configuration is issued by CUDA intrinsic data structures, but we still need a reliable way to get the warp size to take into account any potential evolution. Luckily, we can figure out this size at runtime thanks to the aforementioned warpSize variable. Consequently, here is how we begin a warp-scope kernel in WLP: We now own the bricks to perform WLP, but still lack a user-friendly API. Indeed, it would not be adapted to ask our users to directly use complex formulas without having wrapped them up before in higher-level calls. A first attempt to do so is implemented though macros. As compiletime mechanisms, macros do not cause any runtime overhead. They are also perfectly handled by nvcc, the CUDA compiler. Our previous investigations result in two distinct macros: WARP_BEGIN and WARP_INIT, which respectively mark the beginning of the warp-scope code portion, and correctly fill the warp identifier variable. When WARP_INIT presents no particularities, except the requirement to be called before any operations bringing into play warpIdx, WARP_BEGIN voluntary forgets the block-starting brace following the if statement. By doing so, we expect users to place both opening and closing braces of their WLP code if needed, just as they would do with any other block-initiating keyword.

To sum up, please note once again that this implementation mainly targets to validate our approach. Still, it lays the foundation of a more complete API dedicated to WLP. Unfortunately, macros are not convenient to use. They do not provide control check until compilation. Macros are also quite hard to debug, because they are inlined in the code. Thus, they are not suited to write production codes. Macros were useful in our case to validate the concept though.

AOP DECLINATION OF WLP

In this section, we study the possibility to take advantage of the inputs brought by Aspect Oriented Programming (AOP) [START_REF] Kiczales | Aspect-oriented programming[END_REF] in GPU programming, and especially when using CUDA.

AOP consists in defining entry points where code is inserted in order to modify the program's behavior. These entry points are called pointcuts in AOP. A traditional example of aspects usage is the insertion of a pointcut to wrap a function call. For instance, the function call can be wrapped around the verification of a pointer, or an exception catching block. The point is this change has a small impact on the code thanks to the way aspects are handled. Indeed, the specific operations are externalized, and usually inserted at compile time. This way, operations can be done without being inserted directly in the code.

AOP generally involves a third-party software that will preprocess the source code in order to actually add the equivalent parts matching the aspect directives. For C++, a preprocessor was released in 2002 and is called AspectC++ [START_REF] Spinczyk | AspectC++: an aspect-oriented extension to the c++ programming language[END_REF]. AspectC++ comes with lots of features but is still experimental, and thus does not fully support the C++ syntax and standard constructs. The aspect has to be defined in an "aspect header" (.ah extension) that will be used by the AspectC++ preprocessor (ac++) to weave the C++ code matching the pointcuts.

Our first attempt has been to try to implement WLP through pointcut matching any function called wlp_*(). In our case, we expected to obtain a similar behavior to what plain macros delivered. The conditional statement, previously achieved by the WARP_BEGIN macro, would have been added with the aspect specific code. The warpIdx variable, previously declared and initialized by WARP_INIT, would have been defined through an inlined function. Unfortunately, such an approach could not be completed for several reasons. First WLP would have been available at a kernel scope only, whereas the original implementation allows to turn each programming block into a WLP-run operation. Second, AspectC++ displays troubles parsing CUDA code and the Thrust library, making it unavailable within a CUDA project.. This led us to abandon AspectC++ and to head towards another aspect implementation.

A classic way to implement handcrafted aspects in C++ is by using templates. This design is close to policy classes [START_REF] Alexandrescu | Modern C++ Design[END_REF] where each new behavior that is intended to be inserted is implemented at the heart of a new template class. This class is seen as a component, which can be plugged into the base class or into other components. In order to support this feature, each involved class must accept a template parameter. Moreover, when using C++03, the template function's prototype has to match the wrapped function's to enhance it. This forces the user to write its own aspect any time he wants to use WLP. This is why we abandoned this method. A new feature from C++11 called variadic templates would enable us to provide aspect facilities through templates without the need for the user to write its own template aspect class. However, nvcc, the CUDA compiler, does not support C++11 at the time of writing. In the same way, a third approach could have been considered by harnessing variadic macros. Yet, the latter are not standard in C++03 either.

Other languages, such as Java for instance, implement aspects through annotations inserted in the code. The latter are understood by a preprocessor that makes them part of the language keywords. As long as the C-style fashion to deal with aspects has shown unavailable or awkward in a CUDA-enabled project, we turned to annotations. This implied to write our own preprocessor, which is designed to adopt the same behavior as Java annotations. The preprocessor itself is a simple Perl script that will rapidly parse the code to find the annotations. Its behavior is simple and focuses on detecting the annotations in the original code, and neither evaluates the whole code, nor builds a full syntax tree. This which prevents issues with unsupported C++/CUDA features that are encountered by more complete pieces of software like AspectC++.

Macros make the code harder to debug because they are inlined. An aspect preprocessor can workaround this issue by taking advantage of the #line pragma. This pragma can be used to point to the compiler what will be the next line number to consider, no matter the previous line number. By doing so, when our preprocessor weaves the code, it also inserts the said pragma. The aspect is then totally transparent for the user in case of a build error or when debugging. He will be warned of problems in his code at the exact line they were before any aspect code was inserted.

Concretely, in order to define a WLP kernel, the "// @warp" annotation just needs to be added before the kernel implementation. Two others annotations have been defined to enable WLP on part of the code only. The code parts that need WLP are delimited by "// @warp: begin" and "// @warp: end" at, respectively, the beginning and end of the block. Equivalent keywords from the macro WLP implementation and the AOP version are summed up in Table 1. Obviously, this preprocessing phase occurs prior to the classical building stages from CUDA.

In a more concrete way, Listings 1 and 2 expose dummy code snippets using aspects to enable WLP. Listing 1 shows an example of code for a whole WLP kernel, while Listing 2 presents a code where WLP is bounded to part of the code only.

Aspect code

Macro equivalence // @warp WARP_BEGIN { WARP_INIT // actual code } // @warp: begin WARP_BEGIN { WARP_INIT // @warp: end } Table 1: Equivalence between original WLP through macros and aspect implementation

RESULTS

In this part, we introduce three well-known stochastic simulation models in order to benchmark our solution. We have compared WLP's performances on a Tesla C2050 board to those of a state-of-the-art scalar CPU: an Intel Westemere running at 2.527 GHz. For all of the three following models, each replication runs in a different warp when considering the GPU, whereas the CPU runs the replications sequentially. The following implementations use L'Ecuyer's Tausworthe three-component PRNG, which is available on both CPU and GPU respectively through Boost.Random and Thrust.Random [START_REF] Hoberock | Thrust: A Parallel Template Library[END_REF] libraries. Random streams issued from this PRNG are then split into several sub-sequences according to the Random Spacing distribution technique Hill (2010).

Description of the models

First, we have a classical Monte Carlo simulation used to approximate the value of Pi. The application draws a succession of random points' coordinates. The number of random points present in the quarter of a unit circle are counted and stored. At the end of the simulation, the Pi approximation corresponds to the ratio of points in a quarter of the unit circle to the total number of drawn points. The output of the simulation is therefore an approximate Pi value. The second simulation is a M/M/1 queue. For each client, the time duration before its arrival and the service time is randomly drawn. All other statistics are computed from these values. The program outputs are the average idle time, the average time in queue of the clients and the average time spent by the clients in the system. Because it did not impact the performances, the parameters of the random distribution are static in our implementation. Only the number of clients in the system and the number of replications, which modify the execution time, can be specified when running the application.

The last simulation is an adaptation of the random walk tests for PNRGs exposed in [START_REF] Vattulainen | Mission impossible: Find a random pseudorandom number generator[END_REF]. The idea is to simulate a walker moving randomly on a chessboard-like map. The original application tests the independence of multiple flows of the same PRNG. To achieve this, multiple random walkers are run with different initializations of a generator on identically configured maps. Basically, each walker computes a replication. In the end, we count the number of walkers in every area of the map. Depending on the PRNG quality, we should find an equivalent number of walkers in each area. When the original version splits the map in four quarters, our implementation uses 30 chunks to put the light on the opportunity of our approach when there are many divergent branches in an application.

Comparison CPU versus GPU warp

As we can see in Figure 5, the CPU computation time of the Monte Carlo application approximating the value of Pi grows linearly with the number of replications. The GPU computation time increases only by steps. This behaviour is due to the huge parallel capability of the device. Until the GPU is fully loaded, adding another replication does not impact the computation time, because they are all done in parallel. So, when the device is full, any new iteration will increase the computation time. This only happens on the 65 th replication because the GPU saved some resources in case a new kernel would have to be computed simultaneously. The same mechanism explains that after this first overhead, a new threshold appears and so on.

Due to this behaviour, GPUs are less efficient than CPUs when the board is nearly empty. When less than 30 replications are used, more than two-thirds of the board computational power is idle. Because sequential computation on CPU is widely faster than sequential computation on GPU, if only a little of the parallel capability of the device is used, the GPU runs slower. But when the application uses more of the device parallel computation power, the GPU becomes more efficient than the CPU.

The pattern is very similar for the second model: the M/M/1 queue (see Figure 6). When the board does not run enough warps in parallel, the CPU computation is faster than the GPU one. But with this model, the number of replications needed for the GPU approach to outperform the CPU is smaller than what we obtained with the previous simple model. The GPU computation is here faster as soon as 20 replications are performed, when it required 30 replications to show its efficiency with the first model. This can be explained by GPUs' architecture, where memory accesses are far more costly than floating point operations in terms of processing time. If the application has a better computational operations per memory accesses ratio, it will run more efficiently on GPU. Thus, the GPU approach will catch with the CPU one faster. This point is very important because it means that depending on the application characteristics, it can be adequate to use this approach from a certain number of replications, or not. A solution is to consider the warp approach only when the number of replications is enough to guaranty that most of the applications will run faster.

Comparison GPU warp versus GPU thread

If the advantages of WLP-enabled replications compared to CPU ones in terms of computation time have been demonstrated with the previous examples, it is necessary to determine if outperforms the classic TLP.

This case study has been achieved using the last model introduced: our adaptation of the random walk. Figure 7 shows the computation time noticed for each approach: CPU, GPU with WLP and GPU with TLP (named thread in the caption). Obviously, CPU and WLP results confirm the previous pattern: the CPU computation time increases linearly when the WLP one increases by steps. TLP follows logically the same evolution shape as WLP. Although it is impossible to see it here because the number of replications is too small, it also evolves step by step, similarly to the warp approach. WLP consumes a whole warp for each replication. In the same time, TLP activates 32 threads per warp. Thus, the latter's steps will be 32 times as long as WLP's. Having said that, we easily conclude that the first step in TLP will occur after the 2048 th replication.

As we can see in Figure 7, the computation time needed by the thread approach is significantly more important than the computation time of the warp approach (about 6 times bigger for the first 64 replications). But WLP catches up with TLP when the number of replications increases. When more than 700 replications are performed, the benefit of using the warp approach is greatly reduced. The best use of the warp approach for this model is obtained when running between 20 and 700 replications. Please note that this perfectly matches our replications amount requirement. It even allows the user to run another set of replications according to an experimental plan, or to run another set of replications with a different high quality PRNG. The latter practice is a good way to ensure that the input pseudo-random streams do not bias the results. These results are backed up by the output of the NVIDIA Compute Profiler for CUDA applications. The latter tool allows developers to visualize many data about their applications. In our case, we have studied the ratio between the time spent accessing global memory versus computing data. Such figures are displayed in Figure 8 for both TLP and WLP versions of the random walk simulation. Our approach obviously outperforms TLP, given that the ratio of overall Global Memory access time versus computation time is about 2.5 times bigger for TLP.

To explain this ratio, let us recall that computation time was lower for WLP. Since the same algorithm is computed by the two different approaches, we should have noticed the same amount of Global Memory accesses in the two cases. In the same way, the profiler indicates significant differences between Global Memory reads and writes for TLP and WLP. These figures are summed up in Table 2:

CONCLUSION

This paper has shown that using GPUs to compute MRIP was both possible and relevant. Having depicted nowadays GPUs' architecture, we have detailed how warp scheduling was achieved on such devices, and especially how we could take advantage of this feature to process codes with a high rate of branch divergent parts. Our approach, WLP (Warp-Level Parallelism), can also help users to easily distribute their DOE experimental plans with replications on GPU.

WLP has been implemented thanks to simple arithmetic operations. Consequently, WLP displays a minimalist impact on the overall runtime performance. In order to validate the approach, the internal mechanisms enabling WLP have first been wrapped in macros. Then, for the sake of user-friendliness, another high level version has been proposed following an aspectoriented approach. Aspects are implemented through annotations, preprocessed by a Perl script to generate the corresponding WLP blocks. At the time of writing, our version is functional and allows users to create blocks of code that will be executed independently on the GPU. Each warp will run an independent replication of the same simulation, determined by the warp identifier figured out at runtime. By doing so, we prevent performances to drop as they would do in an SIMT environment confronted to branch-divergent execution paths. WLP also tackles the GPU underutilization problem by artificially increasing the occupancy.

To demonstrate our approach performances, we have compared the execution times of a sequence of independent replications for three different stochastic simulations. Results show that WLP is at least twice as fast as a thread running on a cutting-edge CPU when asked to compute a reasonable amount of replications, that is to say more than 30 replications. This will always be the case when a stochastic simulation is studied with a design of experiments, where for each combination of deterministic factors, at least 30 replications shall be run, according to the previously mentioned Central Limit Theorem. WLP also overcomes the traditional CUDA SIMT performances by up to 6 to compute the same set of replications. Here, SIMT suffers of an underutilized GPU, whereas WLP takes advantage of the fast scheduling of warps.

Insofar performances of WLP increase with the Fermi architecture compared to Tesla, but it is difficult to forecast the same behaviour on new architectures without adapting the approach.

The lack of information regarding the hardware and especially the schedulers ruling threads execution on a CUDA GPU forces us to perform new experiments and possibly adapt WLP.

Still the aspect-oriented declination of WLP opens new perspective in terms of software engineering for GPU devices. We have shown that AOP could be harnessed for such devices, provided an handcrafted preprocessor is available. Future releases of the CUDA toolkit should even withdraw this constraint by allowing established solutions such as AspectC++ and C++11 to fully match CUDA's specificities. Aspect could then deliver their full potential, and in our case, allow us to implement an OpenCL declination of WLP without changing the way it is used in client source code.

Figure 1 :

 1 Figure 1: Simplified architecture of an NVIDIA Fermi C2050 streaming multiprocessor

  Figure 2: Representation of thread disabling to place the application at a warp-level

  Figure 3: Const-definition of warpIdx

Figure 4 :

 4 Figure 4: Directive enabling warp-scope execution

Listing 1 :

 1 A whole WLP kernel // @warp __global__ v o i d T e s t K e r n e l ( f l o a t * d e v i c e A r r a y ) { d e v i c e A r r a y [ 0 ] = 0 . ; d e v i c e A r r a y [ warpIdx ] = 1 . ; }

Listing 2 :

 2 WLP on part of the code only __global__ v o i d T e s t K e r n e l 2 ( f l o a t * d e v i c e A r r a y ) { d e v i c e A r r a y [ 0 ] = 0 . ; // @warp : b e g i n d e v i c e A r r a y [ warpIdx ] = 1 . ; // @warp : end } This model takes two input parameters: the number of random points to draw and the number of replications to compute.

Figure 5 :

 5 Figure 5: Computation time versus number of replications for the Monte Carlo Pi approximation with 10000000 draws

Figure 6 :

 6 Figure 6: Computation time versus number of replications for a M/M/1 queue model with 10000 clients

Figure 7 :

 7 Figure 7: Computation time versus number of replications for a random walk model with 1000 steps (above: 100 replications, below: 1000 replications)

Figure 8 :

 8 Figure 8: Comparison of TLP and WLP ratio of the overall Global Memory access time versus computation time TLP WLP Reads 225 18 Writes 302 104

Table 2 :

 2 Number of read and write accesses to Global Memory for TLP and WLP versions of the Random Walk

partially funded by the Auvergne Regional Council and the FEDER


H Frankowska

E M Marchini
email: elsa.marchini@polimi.it

M Mazzola
email: marco.mazzola@imj-prg.fr

A RELAXATION RESULT FOR STATE CONSTRAINED INCLUSIONS IN INFINITE DIMENSION

In this paper we consider a state constrained differential inclusion ẋ ∈ Ax + F (t, x), with A generator of a strongly continuous semigroup in an infinite dimensional separable Banach space. Under an "inward pointing condition" we prove a relaxation result stating that the set of trajectories lying in the interior of the constraint is dense in the set of constrained trajectories of the convexified inclusion ẋ ∈ Ax + coF (t, x). Some applications to control problems involving PDEs are given.

Proposition 3.3. Assume that (2.8)-(2.9) hold true with time independent k R , φ ∈ R + ,

Introduction

We study a class of infinite dimensional differential inclusions subject to state constraints. Interest in these kind of equations arises in several contexts. Differential inclusions find a natural application in a research area of great development, the control theory, and the infinite dimensional setting allows to apply our results to control problems involving PDEs. Hence, models describing many physical phenomena such as diffusion, vibration of strings, fluid dynamics, may be included in our analysis.

In this paper we are concerned with the differential inclusion ( with coF (t, x(t)) the closed convex hull of F (t, x(t)). The operator A is the infinitesimal generator of a strongly continuous semigroup S(t) : X → X, X is an infinite dimensional separable Banach space, F : I × X X is a set-valued map with closed non-empty images, I = [0, 1] and t 0 ∈ I. The trajectories of the differential inclusion (1.1) are understood in the mild sense (see [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]) and are subject to the state constraint. Namely given a set K ⊂ X, we restrict our attention to the trajectories satisfying (1.3) x(t) ∈ K, for t ∈ [t 0 , 1] .

In this paper we shall always assume that K is the closure of an open subset of X. When satisfying the constraint, a trajectory x is called feasible. Differential inclusions, and control systems, in presence of state constraints, are largely employed in applied sciences. One of the tools playing a key role in this context consists in approximating feasible trajectories by trajectories lying in the interior of the constraints. It is used for instance to establish regularity properties of value functions, to justify the Date: January 1, 2015. This work was co-funded by the European Union under the 7th Framework Programme "FP7-PEOPLE-2010-ITN", grant agreement number 264735-SADCO.

use of the Maximum Principle in normal form, to prove existence and regularity results of optimal solutions. The classical technique employed to construct the approximating trajectories relies on the possibility of directing the velocity into the interior of the constraint K whenever approaching the boundary of K. To this aim, in the finite dimensional setting, an "inward pointing condition" was proposed by Soner, see [START_REF] Soner | Optimal control with state-space constraints[END_REF], to get continuity of the value function associated to an optimal control problem with dynamics ẋ ∈ F (x) independent of t. Since then, this subject has received considerable attention, a partial list of references includes [START_REF] Bettiol | On trajectories satisfying a state constraint: W 1,1 estimates and counter-examples[END_REF][START_REF] Bettiol | L ∞ estimates on trajectories confined to a closed subset[END_REF][START_REF] Cannarsa | On Bolza optimal control problems with constraints[END_REF][START_REF] Frankowska | Discontinuous solutions of Hamilton-Jacobi-Bellman equation under state constraints[END_REF][START_REF] Frankowska | Filippov's and Filippov-Wazewski's theorems on closed domains[END_REF][START_REF] Frankowska | Existence of neighbouring feasible trajectories: applications to dynamic programming for state constrained optimal control problems[END_REF].

Defining the oriented distance from x ∈ X to K by

d K (x) = inf k∈K x -k X if x / ∈ K -inf k∈(X\K) x -k X otherwise,
the inward pointing condition, in the case of time independent F and state constraints with a locally C 1,1 boundary, takes the following form:

(1.4) min v∈F (x) ∇d K (x), v ≤ -ρ, ∀ x ∈ ∂K ,
cf. [START_REF] Bettiol | On trajectories satisfying a state constraint: W 1,1 estimates and counter-examples[END_REF][START_REF] Frankowska | Filippov's and Filippov-Wazewski's theorems on closed domains[END_REF]. As in many applied models state constraints having nonsmooth boundary are present, a number of papers made extensions of (1.4) to the nonsmooth setting. However, contrary to the smooth case, here some regularity of the dynamics F (t, x) is usually required both in t and in x. On the other hand, it may happen in some applications that the dynamics depends on t in a merely measurable way. In order to extend the theory to this situation, in the recent works [START_REF] Frankowska | Discontinuous solutions of Hamilton-Jacobi-Bellman equation under state constraints[END_REF][START_REF] Frankowska | On relations of the adjoint state to the value function for optimal control problems with state constraints[END_REF] a new inward pointing condition (equivalent to the classical one if K has smooth boundary) is proposed: for any "bad" velocity v pointing outside the constraint, there exists a "good" one v such that the difference v -v points inside in a uniform way, namely:

(1.5)

   ∀ x ∈ ∂K, ∃ ρ > 0, ∀ t ∈ I, ∀ v ∈ F (t, x) with sup ξ∈∂d K (x) ξ, v ≥ 0, ∃ v ∈ F (t, x) satisfying sup ξ∈∂d K (x) ξ, v -v < -ρ .
Here ∂d K (x) denotes the Clarke generalized gradient of d K at the point x ∈ X. Under this assumption, in [START_REF] Frankowska | Discontinuous solutions of Hamilton-Jacobi-Bellman equation under state constraints[END_REF][START_REF] Frankowska | On relations of the adjoint state to the value function for optimal control problems with state constraints[END_REF] some approximation results were proved in order to get uniqueness of solutions for a constrained Hamilton-Jacobi-Bellman equation.

The purpose of the present paper is to perform the analysis in the infinite dimensional setting, the natural framework for many phenomena described by PDEs. Also in this case we need results which permit to approximate feasible trajectories by trajectories staying in the interior of the state constraints. Assuming an inward pointing condition, Theorem 3.2 below guarantees the existence of the required approximation. Notice that, although the literature dealing with infinite dimensional control theory (and infinite dimensional differential inclusions) is quite rich, see e.g. the books [START_REF] Barbu | Analysis and control of nonlinear infinite dimensional systems[END_REF][START_REF] Bensoussan | Representation and control of infinite dimensional systems[END_REF][START_REF] Bensoussan | Representation and control of infinite dimensional systems[END_REF][START_REF] Fattorini | Infinite-dimensional optimization and control theory[END_REF][START_REF] Lasiecka | Control theory for partial differential equations: continuous and approximation theories[END_REF][START_REF] Li | Optimal control theory for infinite-dimensional systems[END_REF], the recent paper [START_REF] Cannarsa | Optimal control for evolution equations with memory[END_REF] and the bibliography therein, to our knowledge, no similar results are known in this setting. As an application, we obtain our main result, a relaxation theorem in infinite dimension (see Theorem 3.1).

We deal with great generality, allowing the state space X to be a separable Banach space. Hence, our analysis applies to some interesting and delicate frameworks as the space of essentially bounded functions and the space of continuous functions. For this reason, in this context, the relaxation theorem is obtained under a version of condition (1.5), requiring some uniformity on a neighborhood of ∂K and with respect to the semigroup. Nevertheless, as illustrated in Section 3, if some compactness conditions are satisfied, a much more simple condition, analogue to the finite dimensional (1.5) is sufficient.

We consider the following inward pointing condition:

∀ x ∈ ∂K, ∃ η, ρ, M > 0 such that ∀ t ∈ I, ∀ x ∈ K ∩ B(x, η), (1.6) ∀ v ∈ coF (t, x) satisfying sup τ ≤η, ξ∈∂d K (z 0 )
ξ, S(τ ) v ≥ 0, for some z 0 ∈ B(x, η),

∃ v ∈ coF (t, x) ∩ B(v, M ) satisfying sup z ∈ B(S(τ )x, η) τ ≤ η, ξ ∈ ∂dK (z) ξ, S(τ ) (v -v) < -ρ.
Notice that condition (1.6) deals with the set-valued map coF , since, in order to prove the relaxation theorem, we must be able to approximate relaxed trajectories by relaxed trajectories lying in the interior of K. However, under additional compactness conditions, the first convex hull can be removed from (1.6). Quite interesting for the applications is the case when X is a Hilbert space, see Section 5. In this framework we will provide an alternative version of condition (1.6), which drastically simplifies the analysis when the set of constraints K is convex. The inward pointing condition needed here involves projections on convex sets rather than generalized gradients of the oriented distance function which belong to the dual space X * . 1.1. Outline of the paper. Section 2 contains a list of notations, definitions, and assumptions in use. The main theorems are stated in Section 3. Some results which allow to simplify the inward pointing condition are also proposed. The Hilbert setting is analyzed in Section 4, while Section 5 is devoted to some applications involving PDEs and integrodifferential equations.

Preliminaries

In this Section we list the notation and the main assumptions in use throughout the paper.

2.1. Notation.

-B(x, r) denotes the closed ball of center x ∈ X and radius r > 0; B is the closed unit ball in X centered at 0; -given a Banach space Y , L(X, Y ) denotes the Banach space of bounded linear operators from X into Y , C(I, X) the space of continuous functions from I to X, L 1 (I, X) the space of Bochner integrable functions from I to X, and L ∞ (I, X) the space of measurable essentially bounded functions from I to X; -•, • stands for the duality pairing on X * × X; -µ is the Lebesgue measure on the real line. We will use the following notion of solution.

Definition 2.1. Let t 0 ∈ I and x 0 ∈ X. A function x ∈ C([t 0 , 1], X) is a (mild) solution of (1.1) with initial datum x(t 0 ) = x 0 if there exists a function f ∈ L 1 ([t 0 , 1], X) such that (2.1) f (t) ∈ F (t, x(t)), for a.e. t ∈ (t 0 , 1) and (2.2) x(t) = S(t -t 0 ) x 0 + t t 0 S(t -s) f (s) ds, for any t ∈ [t 0 , 1],
i.e. f is an integrable selection of the set valued map t F (t, x(t)) and x is a mild solution (see [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]) of the initial value problem

(2.3) ẋ(t) = Ax(t) + f (t), for a.e. t ∈ [t 0 , 1] x(t 0 ) = x 0 .
In order to simplify the notation, for t 0 ∈ I and x 0 ∈ X, we denote by S K (t 0 , x 0 ) the set of mild solutions x of (1.1) satisfying x(t 0 ) = x 0 and (1.3), and by f x the corresponding measurable selection in (2.3).

Notice that, since S(t) is a strongly continuous semigroup, there exists M S > 0 such that

(2.4) S(t) L(X,X) ≤ M S , for any t ∈ I.
The differential inclusion (1.1) is a convenient tool to investigate for example the semilinear control system

(2.5) ẋ(t) = Ax(t) + f (t, x(t), u(t)) a.e. t ∈ [t 0 , 1] u(t) ∈ U
where U is an appropriate separable metric space of controls. Setting F (t, x) = f (t, x, U ), we can reduce (2.5) to (1.1) by applying a measurable selection theorem.

2.2.

Assumptions. In our main theorems, we will assume the following conditions:

-positive invariance of K by the semigroup:

(2.6) S(t) K ⊂ K ∀ t ∈ I;
-∀ t ∈ I and any x ∈ X, F (t, x) is closed, and, for any x ∈ X,

(2.7) the set-valued map F (•, x) is Lebesgue measurable; -F (t,
•) is locally Lipschitz in the following sense: for any R > 0, there exists k R ∈ L 1 (I; R + ) such that, for a.e. t ∈ I and any x, y ∈ RB,

(2.8) F (t, x) ⊂ F (t, y) + k R (t) x -y X B;
-there exists φ ∈ L 1 (I; R + ) such that, for a.e. t ∈ I and any x ∈ X, (2.9)

F (t, x) ⊂ φ(t)(1 + x X )B.

The main results

In this section we state the main results of the paper. The first is a relaxation theorem.

Theorem 3.1. Assume (1.6) and (2.6)-(2.9). Then, for any ε > 0 and any feasible trajectory x of (1.2), (1.3), there exists a trajectory x of (1.1) satisfying

(3.1) x(t 0 ) = x(t 0 ), x(t) ∈ Int K, for any t ∈ (t 0 , 1]
and

(3.2) x -x L ∞ (t 0 ,1;X) ≤ ε.
The key point in the proof of Theorem 3.1, is a result on approximation of feasible trajectories, by trajectories lying in the interior of the constraint K. Theorem 3.2. Assume (2.6)-(2.9) and that

∀ x ∈ ∂K, ∃ η, ρ, M > 0 such that ∀ t ∈ I, ∀ x ∈ K ∩ B(x, η), (3.3) ∀ v ∈ F (t, x) satisfying sup τ ≤η, ξ∈∂d K (z 0 ) ξ, S(τ ) v ≥ 0, for some z 0 ∈ B(x, η), ∃ v ∈ F (t, x) ∩ B(v, M ) satisfying sup z ∈ B(S(τ )x, η) τ ≤ η, ξ ∈ ∂dK (z) ξ, S(τ ) (v -v) < -ρ.
Then, for any ε > 0 and any feasible trajectory x of (1.1), (1.3), there exists a trajectory x of (1.1) satisfying (3.1) and (3.2).

In the following propositions pointwise versions of the inward pointing condition (3.3) are proposed, see the applications in Section 5. In the next proposition the convexity of values of F is needed on the boundary of K.

F (•, x) is continuous for any x ∈ X, (3.
Proposition 3.4. Assume that (3.4) and (2.8)-(2.9), with time independent k R , φ ∈ R + , hold true, that X is reflexive, for any x ∈ ∂K and t ∈ I, F (t, x) is convex, and the set-valued map ∂d K (•) is upper semicontinuous in x, and ∂d K (x) is compact. Notice that when d K is C 1 on a neighborhood of ∂K, then condition (3.6) is satisfied. In the proof of Theorem 3.1, we need to approximate relaxed trajectories by relaxed trajectories lying in the interior of K. This is the reason why the inward pointing condition (1.6) required in this case involves the set-valued map coF . By the way, the first convex hull in (1.6) can be removed in some special cases, as indicated in the next proposition. Proposition 3.5. Suppose that for every t ∈ I and x ∈ ∂K, coF (t, x) is closed and the set-valued map

(3.7) [0, 1] t co v ∈ F (t, x) : sup ξ∈∂d K (x)
ξ, v ≤ 0 is upper semicontinuous with closed values. Assume that (3.4) and (2.8)-(2.9), with time independent k R , φ ∈ R + , hold true. Further, suppose that either (3.5) is satisfied, or that X is a reflexive space and (3.6) is satisfied. Then, the following assumption: for any x ∈ ∂K, there exists ρ > 0 such that ∀ t ∈ I and any v ∈ F (t, x) satisfying sup

ξ∈∂d K (x) ξ, v ≥ 0, (3.8) ∃ v ∈ co F (t, x) satisfying sup ξ∈∂d K (x) ξ, v -v < -ρ ; implies (1.6).
In the following remark, the special case of affine forcing terms is analyzed, providing further simplification.

Remark 3.6. If d K is C 1 on a neighborhood of ∂K and, for a subset U ⊂ Y , (3.9) F (t, x) = f 0 (t, x) + g(t, x)U,
where Y is a Banach space,

f 0 : I × X → X and g : I × X → L(Y, X),
then the classical inward pointing condition implies (1.5) with F and also with co F whenever either U is compact or Y is reflexive. Namely, assume (3.4), and (2.8), (2.9) for time independent k R , φ ∈ R + , with F replaced by f 0 and g.

If either U is compact or Y is reflexive, then co F (t, x) = f 0 (t, x) + g(t,
x)co U . Then the inward pointing condition:

∀ x ∈ ∂K, ∀ t ∈ I, ∃ ū ∈ U such that ∇d K (x), f 0 (t, x) + g(t, x)ū < 0 (3.10)
implies (1.5) both with F and with co F . Indeed, by compactness of [0, 1] and continuity of f 0 (•, x) and g(•, x), assumption (3.10) yields:

∀ x ∈ ∂K, ∃ ρ > 0, ∀ t ∈ I, ∃ ū ∈ U satisfying ∇d K (x), f 0 (t, x) + g(t, x)ū < -ρ. (3.11)
Let t ∈ I and u ∈ co U be so that ∇d K (x), f 0 (t, x) + g(t, x)u ≥ 0. Thus

∇d K (x), f 0 (t, x) ≥ -∇d K (x), g(t, x)u .
Then, taking ū as in (3.11), we obtain

∇d K (x),f 0 (t, x) + g(t, x)ū -f 0 (t, x) + g(t, x)u = ∇d K (x), g(t, x)ū -g(t, x)u ≤ ∇d K (x), g(t, x)ū + f 0 (t, x) < -ρ,
yielding (1.5) with F and also with co F .

Under the same assumptions and F given by (3.9), let us consider two examples, where condition (3.10) can be further simplified. Case 1: 0 ∈ U . If ∇d K (x), f 0 (t, x) < 0, for any x ∈ ∂K and t ∈ I, then (3.10) holds for ū = 0. Case 2: U is the unit sphere in Y . Here, if (3.12) ∇d K (x), f 0 (t, x) < g(t, x) * ∇d K (x) Y = 0, for any x ∈ ∂K and t ∈ I, then (3.10) holds for

ū = - g(t, x) * ∇d K (x) g(t, x) * ∇d K (x) Y , (3.13)
where g(t, x) * is the adjoint of g(t, x). Indeed, in this case, for any x ∈ ∂K, ∇d K (x), f 0 (t, x) + g(t, x)ū = ∇d K (x), f 0 (t, x) -g(t, x) * ∇d K (x) Y < 0, yielding (3.10).

The case of Hilbert spaces

Here we analyze the case when the state space X is Hilbert. In this setting, we show that if the state constraint is convex then the inward pointing condition can be drastically simplified by involving projections on convex sets instead of generalized gradients of the oriented distance function which do belong to the dual space X * . This turns out to be very useful in the applications, as we will show in Section 5.

Let •, • X be the scalar product in X and let K be a proper closed subset of X such that K = Int K. Denote by Z the set of points z ∈ X ∂K admitting a unique projection P ∂K (z) on ∂K. This set is dense in X (see [START_REF] Preiss | Differentiability of Lipschitz functions on Banach spaces[END_REF]). For every z ∈ Z, set

n z = z -P ∂K (z) z -P ∂K (z) X sgn(d K (z)) .
A new inward pointing condition involving n z is proposed in this Hilbert framework in order to obtain results analogous to those from Section 3.

Theorem 4.1. Assume (2.6)-(2.9). Then, (i) the assertions of Theorem 3.1 are valid under the following inward pointing condition:

∀ x ∈ ∂K, ∃ η, ρ, M > 0 such that ∀ t ∈ I, ∀ x ∈ K ∩ B(x, η), (4.1) ∀ v ∈ co F (t, x) satisfying sup τ ≤η, z∈Z∩B(x,η) n z , S(τ ) v X ≥ 0, ∃ v ∈ co F (t, x) ∩ B(v, M ) satisfying sup τ ≤η, z∈Z∩B(S(τ )x,η) ξ, S(τ ) (v -v) X < -ρ.
(ii) the assertions of Theorem 3.2 are valid under the following inward pointing condition:

∀ x ∈ ∂K, ∃ η, ρ, M > 0 such that ∀ t ∈ I, ∀ x ∈ K ∩ B(x, η), (4.2) ∀ v ∈ F (t, x) satisfying sup τ ≤η, z∈Z∩B(x,η) n z , S(τ ) v X ≥ 0, for some z 0 ∈ B(x, η), ∃ v ∈ F (t, x) ∩ B(v, M ) satisfying sup τ ≤η, z∈Z∩B(S(τ )x,η) n z , S(τ ) (v -v) X < -ρ.
Again, these conditions can be simplified when the data satisfy some compactness assumptions.

Proposition 4.2. Assume that (2.8)-(2.9) with time independent k R , φ ∈ R + and (3.4) hold true. Further suppose that either (3.5) is valid, or F (t, x) is convex for any t ∈ I and x ∈ ∂K, and

∀ x ∈ ∂K, ∃ r > 0 such that the set n z : z ∈ Z ∩ B(x, r) is pre-compact. (4.3)
Then, the following assumption: for any x ∈ ∂K, there exists ρ > 0 such that for any t ∈ I and v ∈ F (t, x) satisfying inf ε>0 sup z∈Z∩B(x,ε) n z , v X ≥ 0, (4.4)

there exists v ∈ F (t, x) satisfying inf ε>0 sup z∈Z∩B(x,ε) n z , v -v X < -ρ . implies (4.2).
Remark 4.3. The proof of Proposition 4.2 provided in Section 6 implies that it is still valid if (4.3) is replaced by the following less restrictive assumption: for x ∈ ∂K define N (x) :=Limsup z→x, z∈Z {n z } (the Kuratowski upper limit) and assume that for all x ∈ ∂K the set N (x) is compact and for every ε > 0 there exists δ > 0 such that

n z ∈ N (x) + εB ∀ z ∈ Z ∩ B(x, δ).
In particular, if ∂K is of class C 1 , then the above holds true.

4.1. Convex state constraints. If K convex, then the inward pointing conditions (4.1), (4.2), and (4.4) can be weakened by replacing Z with K c := X \ K. Indeed, any z ∈ K c admits a unique projection on ∂K and, as proved in the following proposition, for any z ∈ Int K ∩ Z we can find an element w ∈ K c such that n z = n w .

Proposition 4.4. Let K be a closed convex set such that K = Int K. Then, for any z ∈ Int K ∩ Z, there exists w ∈ K c such that z -P ∂K z = P ∂K w -w. In particular, n z = n w .

Proof. Let z ∈ Int K ∩Z and P ∂K (z) be its unique projection on ∂K. By the Hahn-Banach theorem, there exists p ∈ X such that p X = 1 and p, P ∂K (z) X ≤ p, k X , for any k ∈ K.

Let M + = x ∈ X : p, x -P ∂K (z) X ≥ 0 ⊇ K and M = ∂M + = x ∈ X : p, x -P ∂K (z) X = 0 .
Then M is a closed hyperplane in X and there exists a unique projection P M (z) of z on M. Actually, since P ∂K (z) ∈ M,

z -P M (z) X ≤ z -P ∂K (z) X ,
and, since K ⊂ M + and z lies in the interior of K,

z -P M (z) X ≥ z -P ∂K (z) X ,
we deduce that P M (z) = P ∂K (z). Take

w = z + 2(P ∂K (z) -z). As z ∈ Int K ⊂ Int M + , we have p, w -P ∂K (z) X = p, P ∂K (z) -z X < 0, yielding w ∈ X \ M + ⊂ K c . Further, for any x ∈ M, 0 = z -P M (z), x -P M (z) X = z -P ∂K (z), x -P ∂K (z) X = P ∂K (z) -w, x -P ∂K (z) X .
This implies that P ∂K (z) = P M (w). Finally, since M + is a closed convex set, w admits a unique projection P M + (w) = P M (w) = P ∂K (z). So, for any k ∈ K ⊂ M + , w -P ∂K (z), k -P ∂K (z) X ≤ 0, implying P ∂K (z) = P ∂K (w). This ends the proof.

Examples

The examples analyzed in this section describe some classical models involving partial differential equations and integrodifferential equations, to which we may apply our abstract results. In all the examples the state constraints satisfy the positive invariance property (2.6).

5.1.

A one-dimensional heat equation. The first example is a one-dimensional parabolic equation describing the heat flux in a cylindrical bar, whose lateral surface is perfectly insulated and whose length is much larger than its cross-section. The Neumann boundary conditions are assumed, corresponding to the requirement that the heat flux at the two ends of the bar is zero. For x = x(t, s) : [0, 1] × [0, 1] → R we consider the following inclusion (we omit the variable s in the sequel)

(5.1) ẋ(t) ∈ Ax(t) + F (t, x(t)), t ∈ [0, 1] x(0) = x 0 .
The state space is X = H 1 (0, 1) and the linear operator acting as Ax = x -x with domain

D(A) = x ∈ H 2 (0, 1; R) : x (0) = x (1)
= 0 is the infinitesimal generator of a strongly continuous semigroup S(t) : X → X, see e.g. [START_REF] Temam | Infinite-dimensional dynamical systems in mechanics and physics[END_REF]chapter II]. (The notation prime stands for the distributional derivative.) Classical results in PDEs ensure that, if the initial datum x 0 takes nonnegative values, then the solution x to ẋ(t) = Ax(t), x(0) = x 0 takes nonnegative values. The reader is referred to [START_REF] Alikakos | An application of the invariance principle to reaction-diffusion equations[END_REF] or [START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF], containing a number of examples of sets invariant under the action of the semigroup associated with A. In particular, if the state constraint is the cone of nonnegative functions:

x(t) ∈ K = x ∈ X : x ≥ 0 ,
then the invariance property (2.6) is satisfied. Moreover, K is convex and Int K = ∅. The state space is endowed with the scalar product

x, y X = x(0)y(0) + x , y L 2 (0,1) , for any x, y ∈ X, whose associated norm

x 2 X = |x(0)| 2 + x 2 L 2 (0,1)
for any x ∈ X, turns out to be equivalent to the usual one • H 1 (0,1) . We show next that the set (5.2)

n z : z ∈ K c is pre-compact.
Since K is a closed convex cone, then any z ∈ K c can be uniquely represented as

z = P ∂K (z) + b(z),
with b(z) ∈ K -, here K -is the negative polar cone to K. By [START_REF] Wierzbicki | Projection on a cone, penalty functionals and duality theory for problems with inequaltity constraints in Hilbert space[END_REF],

K -= p ∈ X : p is nondecreasing and p(0) ≤ p (s) ≤ 0, for a.e. s ∈ [0, 1] ,

see also [START_REF] Outrata | An augmented Lagrangian method for a class of convex continuous optimal control problems[END_REF] where an explicit formula for b is provided. To prove (5.2), notice that

n z = b(z) b(z) X : z ∈ K C ⊂ Q := p p X : p ∈ K -, p = 0 ⊂ ∂B.
Any y ∈ Q satisfies y nondecreasing and -1 ≤ y(0) ≤ y (s) ≤ 0, for a.e. s ∈ [0, 1].

So, taking a sequence {y n } in Q,

y n W 1,∞ (0,1) := y n L ∞ (0,1) + y n L ∞ (0,1) ≤ 3,
implying that y n (0) → y(0) (up to a subsequence). Since y n is nondecreasing, Helly's selection theorem, see [START_REF] Kolmogorov | Introductory real analysis[END_REF], allows to deduce that, (again up to a subsequence), y n (s) → g(s), for a.e. s, with g ∈ L ∞ (0, 1), and, applying Lebesgue dominated theorem, we deduce that y n → g in L 2 (0, 1). Further,

y n (s) = y n (0) +
n z , v -v X < -ρ .

5.2.

Fourier's problem of the ring. In the second example we consider the temperature distribution in a homogeneous isotropic circular ring with diameter small in comparison with its length and perfectly insulated lateral surfaces. This problem can be modeled by a one-dimensional equation with periodic boundary conditions

(5.3) ẋ(t) ∈ Ax(t) + F (t, x(t)), t ∈ [0, 1] x(0) = x 0 , where x = x(t, s) : [0, 1] × [0, 1] → R (s is omitted as in the previous example), the state space is X = H 1 per (0, 1) := x ∈ H 1 (0, 1; R) : x(0) = x(1)
, the linear operator acting as Ax = x with domain D(A) = H 2 (0, 1; R) ∩ H 1 per (0, 1) is the infinitesimal generator of a strongly continuous semigroup S(t) : X → X, see e.g. [START_REF] Butzer | Semi-groups of operators and approximation[END_REF]. As before we supplement inclusion (5.3) with the state constraint

x(t) ∈ K = x ∈ X : x ≥ 0 .
Then, K satisfies condition (2.6), see for instance [START_REF] Martin | Invariant sets for perturbed semigroups of linear operators[END_REF] dealing with invariant sets for semigroups. Again, K is a closed and convex set with non empty interior. Hence, by the results contained in subsection 4.1, the inward pointing conditions (4.1), (4.2), and (4.4) hold true with Z replaced by K c .

5.3.

A model for Boltzmann viscoelasticity. The last example deals with the phenomena of isothermal viscoelasticity. An integrodifferential inclusion is involved, since, as outlined in the seminal works of Boltzmann and Volterra [START_REF] Boltzmann | Zur Theorie der elastischen Nachwirkung[END_REF][START_REF] Boltzmann | Zur Theorie der elastischen Nachwirkung[END_REF][START_REF] Volterra | Sur les équations intégro-différentielles et leurs applications[END_REF][START_REF] Volterra | Leçons sur les fonctions de lignes[END_REF], a correct description of the mechanical behavior of elastic bodies requires the notion of memory. The key assumption in this theory is that both the instantaneous stress and the past stresses influence the evolution of the displacement function y = y(x, t) : Ω × R → R. Here Ω ⊂ R 3 , a bounded domain with smooth boundary ∂Ω, represents the region occupied by the elastic body. Omitting the variable x in the sequel, we study the following inclusion

(5.4) ÿ(t) + A y(t) - ∞ 0 µ(s)y(t -s) ds ∈ F(t, y(t)), t > 0,
where, A = -∆ with domain D(A) = H 2 (Ω) ∩ H 1 0 (Ω), according to the assumption that the body is kept fixed at the boundary of Ω, the memory kernel µ, taking into account the viscoelastic behavior, is supposed to be a (nonnegative) nonincreasing and summable function on R + , with total mass κ = ∞ 0 µ(s)ds ∈ (0, 1), piecewise absolutely continuous, and thus differentiable almost everywhere with µ ≤ 0. Equation (5.4) is supplemented with the following initial condition

y(0) = y 0 , ẏ(0) = z 0 , y(-s) |s>0 = φ 0 (s),
for some prescribed data y 0 , z 0 , φ 0 , the latter taking into account the past history of y. Applying Dafermos' history approach, see [START_REF] Dafermos | Asymptotic stability in viscoelasticity[END_REF], we can write (5.4) as a differential inclusion of type (1.1). To this aim, we first introduce an auxiliary variable which contains all the information about the unknown function up to the actual time η t (s) = y(t) -y(t -s), t ≥ 0, s > 0 and we recast problem (5.4) as the system of two variables y = y(t) and η = η t (s) with initial conditions y(0) = y 0 , ẏ(0) = z 0 , η 0 = η 0 = y 0 -φ 0 .

Here the operator T is the infinitesimal generator of the right-translation semigroup on the memory space M = L 2 µ (R + ; H 1 0 (Ω)), namely, T η = -η with domain dom(T ) = η ∈ M : η ∈ M, η(0) = 0 .

The notation prime standing for the distributional derivative, and η(0) = lim s→0 η(s) in H 1 0 (Ω). In [START_REF] Cannarsa | Optimal control for evolution equations with memory[END_REF], details and related bibliography can be found, jointly with some applications of this model to optimal control problems. Now, defining the linear operator A on the state space X = L 2 (Ω) × H x(t) = (y(t), z(t), η(t)),

x 0 = (y 0 , z 0 , η 0 ), F (t, x(t)) = 0, F(t, y(t)), 0 , we view (5.5) as the following problem in X:

(5.6) ẋ(t) ∈ Ax(t) + F (t, x(t)), for t ∈ I x(0) = x 0 .

The operator A generates a strongly continuous semigroup of contractions S(t) : X → X whose first component is a solution of (5.4) for F = 0. Here, taking as a state constraint in (5.6) K = B, we deduce that (2.6) is satisfied. Further, let U, f 0 and g be as in Remark 3.6. Since K = B in a Hilbert space, an appropriate adaptation on the basis of Section 4 of the inward pointing condition (3.10) reads: for any x ∈ X with x = 1 and any t ∈ I, there exists ū ∈ U such that x, f 0 (t, x) + g(t, x)ū < 0. In particular in the case when U is the unit sphere in R N , the condition x, f 0 (t, x) < g(t, x) * x R N = 0, implies the above inward pointing condition, for ū defined by ū = -g(t,x) * x g(t,x) * x R N . As discussed in Remark 3.6, Proposition 3.5 holds in this case implying the validity of Theorems 3.1 and 3.2.

  4) and F (t, x) is compact, for any t ∈ I and any x ∈ ∂K. (3.5) Then, assumption (1.5) implies (3.3). Consequently, if (1.5) holds true with F replaced by co F , then (1.6) is satisfied.

(3. 6 )

 6 Then, assumption (1.5) implies (3.3) and (1.6).

s 0 y

 0 n (τ )dτ → y(s) := y(0) + s 0 g(τ )dτ, for any s ∈ [0, 1].Hence g = y , y ∈ W 1,∞ (0, 1) ⊂ H 1 (0, 1), yielding the required pre-compactness. Let F satisfy assumptions (2.8)-(2.9) with time independent k R , φ ∈ R + , (3.4), and let F (t, x) be convex, for any x ∈ ∂K and t ∈ I. Taking into account Proposition 4.2 and the results in subsection 4.1, the inward pointing condition (4.2) is implied by the following assumption: for any x ∈ ∂K there exists ρ > 0 such that for any t ∈ I and v ∈ F (t, x) satisfying inf ε>0 sup z∈K c ∩B(x,ε) n z , v X ≥ 0, there exists v ∈ F (t, x) satisfying inf ε>0 sup z∈K c ∩B(x,ε)

  ) + A (1 -κ)y(t) + ∞ 0 µ(s)η t (s)ds ∈ F(t, y(t)), ηt = T η t + ẏ(t).

  1 0 (Ω) × M, acting asA(y, z, η) = z, -A (1 -κ)y + ∞ 0 µ(s)η(s) ds , T η + z with domain dom(A) = (y, z, η) ∈ X z ∈ H 1 0 (Ω) , η ∈ dom(T ), (1 -κ)y + ∞ 0 µ(s)η(s)ds ∈ H 2 (Ω) ∩ H 1 0(Ω) and setting

Appendix

This section contains some technical results needed in the course of the investigation. The first is an infinite dimensional version of the Filippov Theorem, see [START_REF] Frankowska | A priori estimates for operational differential inclusions[END_REF]Theorem 1.2], modified here in a suitable form for our scopes.

Lemma A.1. Let δ 0 ≥ 0 and t 0 ∈ I. Assume (2.7)-(2.8), let y be a solution to (2.3) in [t 0 , 1], for some f ∈ L 1 (t 0 , 1; X). Set R = 1 2 max t∈[t 0 ,1] y(t) X , γ(t) = dist f (t), F (t, y(t))

and

then, for any y 0 ∈ y(t 0 ) + δ 0 B and any β > 1, there exists a solution x to (1.1)

and

where f x ∈ L 1 (t 0 , 1; X) is so that (2.1) and (2.2) hold true for f = f x .

Proof. The proof proceeds exactly as in [START_REF] Frankowska | A priori estimates for operational differential inclusions[END_REF]Theorem 1.2]. The only difference is in the first line of page 109, while applying Lemma 1.3 of [START_REF] Frankowska | A priori estimates for operational differential inclusions[END_REF]. The assumption β > 1 is needed to ensure that, for a.e. t, the set

Indeed, if γ(t) = 0, the definition of γ ensures that f (t) ∈ F (t, y(t)), while, if γ(t) > 0, since β > 1, from the very definition of distance and the measurable selection theorem we get that there exists a measurable selection

Notice that in finite dimension we can take β = 1, recovering the original Filippov Theorem.

The second result is a version of the mean value theorem for the oriented distance d K in Hilbert spaces. Here we make use of the notations introduced in Section 4.

Lemma A.2. Let (X, , X ) be a Hilbert space. For every x, y ∈ X we have

Proof. Let x, y ∈ X and fix ε > 0. It is enough to consider the case x = y. For every s ∈ [0, 1] set γ(s) = x + s(y -x). Then d K (γ(•)) is absolutely continuous. It is sufficient to prove that for almost every s ∈ [0, 1]

Indeed, then it follows that

The maps d K (γ(•)) and d 2 K (γ(•)) are almost everywhere differentiable. Denote by D ⊂ (0, 1) the set on which both functions are differentiable and fix s ∈ D. We distinguish three cases.

Case 1:

Consequently, for all h > 0 small enough,

Then we obtain d ds

Case 2: γ(s) ∈ Int K. Similarly to Case 1, for every 0 < h < min{ε, s} and every z ∈ Z ∩ (B(γ(s), h 2 ) K c ) we have (A.1) and

Then we obtain d ds It remains to verify that in the case where d ds d K (γ(s)) = 0, we have sup z∈Z∩B(γ(s),ε) n z , y -x X ≥ 0 .

Assume by contradiction that there exists α > 0 such that sup z∈Z∩B(γ(s),ε)

2 ) ∩ Int K converge to γ(s) as i → +∞ and set, for every 0 ≤ r ≤ h, γ i (r) = ξ i + r(y -x). Fix i ∈ N and define

n z , y -x X ≤ -α for a.e. r ∈ [0, h i ] .

Then, γ i (h i ) ∈ Int K and, consequently, h i = h. Summarizing, we have obtained that for every i ∈ N d dr d K (γ i (r)) ≤ sup z∈Z∩B(γ i (r), ε 4 )

n z , y -x X ≤ -α for a.e. r ∈ [0, h] .

Therefore, for every i ∈ N and h ∈ (0, ε 4 y-x X ) we have d K (γ i (h)) -d K (γ i (0)) ≤ -α h . Taking the limit as i → +∞, we obtain d K (γ(s + h)) -d K (γ(s)) ≤ -α h . This is impossible, since d ds d K (γ(s)) = 0.


Pablo Ruiz
email: pruiz@vicomtech.org

Montse Cuadros
email: mcuadros@vicomtech.org

Thierry Etchegoyhen
email: tetchegoyhen@vicomtech.org

Normalización Léxica

Lexical Normalization of Spanish Tweets with Preprocessing Rules, Domain-Specific Edit Distances, and Language Models

Keywords: Spanish microtext, lexical normalization, Twitter, edit distance, language model microtexto, español, castellano, normalización léxica, Twitter, distancia de edición, modelo de lengua

published or not. The documents may come

Introduction and objectives

Studies on the lexical normalization of Spanish microtext are scarce, e.g. Armenta et al. 2003, which predates Twitter and focuses on SMS. The lack of resources and tools for the normalization of Spanish tweets lead us to develop a baseline system, identifying sources of error and means of improvement. The system comprises data resources to model the domain, as well as analysis modules.

Evaluating the system, we identified sources of error in dictionary and entity data, or in candidate-ranking elements like edit cost estimation and language model querying.

The Spanish tweet normalization system presented in this paper achieved above-average performance among a set of 13 competing systems at SEPLN's 2013 Tweet-Norm task 1 .

The system's architecture and components are presented in Section 2, resources employed in Section 3, and settings and results-evaluation 1 www.sepln.org/?news=workshop-tweet-norm &lang=en in Section 4. Conclusions and future work are discussed in Section 5.

Architecture and components

The system's architecture and components are shown in Figure 1 and explained in following.

Preprocessing

The preprocessing module was based on regexes and custom lists.

A set of case-insensitive regexes detected emoticons and delengthened OOVs with repeated characters, as well as mapping OOVs to DRAE2 onomatopoeias. Repeated letters were reduced to a single letter, unless a word with a repeated letter was found in Aspell's Spanish inflected form dictionary (v1.11.3)3 . E.g. vinoo was preprocessed to vino, but creeeen gives creen.

Custom lists were used to identify abbreviations, and expand them if needed. Lists were also used to resegment tokens commonly written together in microtext. Microtext expressions such as RT or HT were considered in-vocabulary (IV).

Candidate generation

A minimum edit distance technique was used to obtain candidate corrections [START_REF] Damerau | A technique for computer correction of spelling errors[END_REF]. Up to two character edits (insertions, deletions, substitutions) were performed on preprocessed OOVs. Variants found in Aspell's dictionary were accepted as correction-candidates. The OOV itself was part of the candidate-set, since it's necessary to determine whether to keep the OOV (e.g. for proper nouns) or to edit it.

Candidate ranking

A candidate's rank prior to named-entity detection reflected the weighted combination of language model (LM) scores and edit distance.

For LM scores, the content of the n-gram looked up in the LM was configurable: either the trigram ending in the candidate, or a string with the candidate as the middle token, and up to nine tokens long. The LM score was the n-gram logprob returned by model lookup.

For distance scores, the Levenshtein distance between the OOV and each candidate was obtained. Each edit was assigned a cost. Costs were domain-specific, determined by surveying common errors in Spanish microtext. E.g. editing k as q (like in one of the edits from kiero to quiero) cost less than uncommon edits.

Costs were also inspired by spelling error frequencies for Spanish, reported by Ramirez and López, 2006. E.g. they found that 51.5% of errors were accent omissions. Accordingly, a cost model was created where replacing a nonaccented character with its accented variant cost less than other substitutions.

Candidates where edit cost was higher than a threshold could be demoted or filtered out (configurable). A weighted sum of LM and distance scores was used to rank candidates.

Only the highest-ranked candidate moved forward to the next step in the workflow.

Entity detection

This step determines whether the highestranked candidate so far (as per LM and editdistance scores) should be considered final, or if an entity-candidate should be proposed instead.

Exploiting NERC models included in FreeLing 3.0 4 was assessed. The models 4 nlp.lsi.upc.edu/freeling/ showed a strong tendency to label uppercaseinitial tokens as entities. This was not optimal, since microtext does not reliable follow standard casing rules for proper-noun identification. FreeLing's confidence scores were also not a reliable indication of a token's proper vs. common noun status.

Given the difficulties in employing NERC, the following heuristics were implemented. For the best-ranked candidate according to the LM and edit distance scores, an uppercase-initial and an all-caps variant were created (e.g. Messi and MESSI for OOV messi). These variants were looked up in entity dictionaries. In case of a match, two factors determined whether to keep the non-entity candidate, or to accept an entity candidate. First, edit distance between the non-entity candidate and the original OOV. If larger than a threshold, the non-entity candidate was demoted. Second, segmental traits of the entity vs. non-entity candidate. E.g. if the entity contained sequences uncommon in Spanish (e.g. uu), it was demoted. Accented characters, on the other hand, promoted either candidate.

Postprocessing

The original OOV's case could undergo decasing via regex application or candidate-set generation. The selected candidates were recased to match the original OOV's case.

Tokens were also uppercased when they followed a sentence delimiter. 5

Resources

As a known-word (IV) dictionary, DRAE was approximated by using Aspell's dictionary. Entity lists were obtained from the JRC Names6 database. A list of named entities manually annotated in the SAVAS corpus7 , which consists of 200 hours of Spanish news broadcasts from 2012 was also used. The SAVAS corpus was useful since it covers entities from current events, often discussed on Twitter. An entity-list for internal use at Vicomtech-IK4 was also employed.

Normalization does not require entity classification or linking, but merely identifying whether a token belongs to an entity or not. Accordingly, in our entity lists multiword entities were split into their tokens. Tokens for which a lowercase variant exists in Aspell's dictionary were filtered out.

For measuring candidate distance, a cost matrix was created. Costs were estimated by surveying the frequency of character substitutions in Spanish tweets, and inspired by [START_REF] Ramírez | Spelling Error Patterns in Spanish for Word Processing Applications[END_REF]. Table 1 provides example costs. Using the table, editing alli to allí costs 0.5; kiero to quiero costs 1.5. Other cost models were also assessed (see Section 4).

Error

Correction Cost (each) a, e, i, o, u, n á, é, í, ó, ú, ñ 0.5 k, null q, u 0.75 p, a, z m, u, k 1

Table 1: Edit Costs

We created three 5-gram case-sensitive language models with Kenlm8 , using an unk token. The corpora for the LMs contained one million sentences each, and vocabulary size was around 140,000. The Subtitles corpus contained film and documentary subtitles. The Tweets corpus contained Spanish tweets with no hashtags, usernames, or URLs. We also used a corpus extracted from Europarl [START_REF] Koehn | Europarl: A Parallel Corpus for Statistical Machine Translation[END_REF].

For the Tweets corpus, tweets with language value es and European time zones were collected in the spring of 2013. Only tweets for which Hunspell (v1.3.2) detected no errors were accepted. In order to decrease false positives, Hunspell dictionaries were enriched with entity lists. Tweet tokenization largely treated emoticons, URLs and repeated punctuation as single tokens. For tweets where there was at least 70% of token-overlap with other tweets, only one exemplar was accepted.

Settings and evaluation

Accuracy with the baseline settings (delivered for Tweet-Norm) was 0.610 for the dev corpus (475 sentences, 654 annotated OOVs) and 0.606 for the test corpus (564 sentences, 632 annotated OOVs). Edit distance and LM scores were weighted equally, and the LM had been trained on the Subtitles corpus (see Section 3).

Results with an LM trained on the Tweets corpus were comparable (0.599 for the test corpus and 0.594 for dev). Results with the Europarl LM dropped 5%. This suggests that LMs for tweet normalization can be trained correctly with off-domain texts containing short sentences and showing colloquial vocabulary and style, such as film subtitles.

The main sources of error for the baseline settings, based on a sample of 200 errors, are shown in Table 2. Data problems account for 45% of errors. Adding abbreviations, and regexes for delengthening or mapping tokens to their DRAE form, would improve results. Other problems are gaps in dictionaries (e.g. missing Twitter slang, or DRAE word pachanga missing from Aspell), incorrect data (e.g. bieber in lowercase in JRC Names), or ambiguous data, such as Dutch name Ruben (from JRC Names) competing with the Spanish form Rubén. Adding a domain dictionary and accessing DRAE instead of Aspell would help solve these issues, as well as exploiting JRC Names' metadata to disallow foreign names in conflict with a Spanish name.

Issues with candidate scoring and ranking via the LM and distance metrics account for most of the remaining problems. Correction Model issues in Table 2 refer to cases where the distance score between the correct candidate and the OOV is too low to compete with candidates that have a better LM score in the context. The distance model in the baseline system is not context sensitive, e.g. replacing k with q and deleting u cost 0.75 each, regardless of context. We tested adding context sensitivity at the character level to the model. We defined a cost of 0.5 for corrections that involve two edits, but repair common errors in the domain, e.g. rendering ki as qui, x as ch or wa as gua. Distance-1 variants common in the domain that were given a cost of 0.5 are p for pe, t for te or k for ka, and others following the pattern of using a character as a shorthand for the sequence that sounds like the character's name. Corrections for colloquial generic variants like ao for ado were also given a cost of 0.5. Accuracy improved 1.06 % in the test corpus and 0.61 % in the dev corpus. The numbers are small, but no overcorrection occurred as a result of the modifications, which suggests a positive trend.

As for LM and LM Querying errors, the high proportion of OOVs in microtext makes it difficult to find relevant contexts to rank candidates. The following example illustrates a common pattern in the genre. The underlined tokens are OOV in @Idoia_LIA buenoa dias mi vida. Our system wrongly rendered buenoa as buena instead of buenos. The distance is 1 for both candidates. The token following OOV buenoa is OOV dias. Given that buenos días is a strong collocation, the corrected token, días, would be a useful feature for the LM to promote the correct candidate, buenos. However the tokens adjacent to the OOV under consideration (buenoa) are also OOV. Thus the LM scores are restricted to the unigram probability of candidates for the OOV to be normalized, without the benefit of context.

These difficulties to exploit LM information led us to test prioritizing distance scores over LM scores. By weighting distance scores at 70% and LM scores at 30% accuracy improved by 1.36% in the test set and by 0.76% in the development set.

Issues with Ranking and Selection Criteria in Table 2 are also related to the interaction between LM and distance scores. We implemented rules to modify the ranking for candidates that are missing an accent, but are otherwise exactly like the OOV. It is very likely for such candidates to be correct; accentomissions cover 51.5% of the errors in [START_REF] Ramírez | Spelling Error Patterns in Spanish for Word Processing Applications[END_REF]. When candidates whose only error is an accent omission were ranked second, we required the difference between their LM score and the LM score of the first-ranked candidate to be higher than a threshold. By setting this threshold at 1.5, accuracy improved by 2.07% in the test set and by 0.76% in the development set.

Another way to improve candidate selection was the following: In the baseline, candidates at distance 2 were demoted, but allowed to compete with each other and with the OOV. This caused some overcorrection of IV tokens and of correct proper nouns, e.g. hombrecillos for pobrecillos. By filtering out candidates above a distance of 1.5, accuracy increased by 0.6 % in the test set and 2.29 % in the dev set.

Conclusions and future work

We presented a baseline system for the normalization of Spanish tweets. The system uses rules to preprocess OOVs into forms closer to IV tokens. Edit-distance scores and language model probabilities rank correction candidates. Edit-costs were estimated taking into account common errors in tweets. Candidate assessment via the language model posed difficulties, given the high frequency of OOVs in the candidates' context.

We also showed configurable settings that improve results by 5.09 % (test corpus) and 4.42 % (dev corpus).

In terms of future work, data problems caused 45% of the system's errors: additional preprocessing rules, IV entries, domain-specific entries like Twitter slang, or more accurate named entities would improve results

Context modeling is also an area to improve. [START_REF] Han | Lexical normalisation of short text messages: makn sens a #twitter[END_REF] devised a context model with features based on non-contiguous IV tokens, and used fuzzy matching between candidates and the context model.

Figure 1 :

 1 Figure 1: System Architecture 5 The delimiters considered were . ! ? " …

Table 2 :

 2 System's Errors

Spanish Academy dictionary, www.rae.es

aspell -l es dump master | aspell -l es expand

optima.jrc.it/data/entities.gzip

www.fp7-savas.eu/savas_project

kheafield.com/code/kenlm/