Sylvain Lefebvre

Samuel Hornus

Anass Lasram

HA-Buffer: Coherent Hashing for single-pass A-buffer

Keywords: A-buffer, real-time, transparency, spatial hashing Hash

Identifying all the surfaces projecting into a pixel has several important applications in Computer Graphics, such as transparency and CSG. These applications further require ordering, in each pixel, the surfaces by their distance to the viewer. In real-time rendering engines, this is often achieved by recording sorted lists of the fragments produced by the rasterization pipeline. The major challenge is that the number of fragments is not known in advance. This results in computational and memory overheads due to the necessary dynamic nature of the data-structure. Similarly, many fragments which are not useful for the final image-due to opacity accumulation for instance-have to be stored and sorted nonetheless, negatively impacting performance. This paper proposes a novel approach which records and simultaneously sorts all fragments in a single geometry pass. The storage overhead per fragment is typically lower than 8 bits per record, and no pointers are involved. Since fragments are progressively sorted in memory, it is possible to assess during rendering whether a new fragment is useful. Our approach combines advantages of previous approaches at similar levels of performance, and is implemented in a single fragment shader of 24 lines of GLSL.

HA-Buffer: application du hachage cohérent à la création de A-buffer en une passe

Résumé : Plusieurs applications en synthèse d'image nécessitent le calcul de l'ensemble des surfaces visibles au travers d'un pixel. Citons le dessin correct de surfaces transparentes ainsi que le dessin de modèles CSG. Ces applications nécessite également de trier les surfaces, pour chaque pixel, selon leur distance au point de vue.

Pour les applications en temps-réel, ce sont les fragments produits par l'étape de rasterisation qui sont triés et stockés en mémoire vidéo. Le nombre de ces fragments n'étant pas connu à l'avance, il est nécessaire d'utiliser de coûteuses techniques de gestion de la mémoire. De plus, tous les fragments sont traités même si une fraction non négligeable d'entre eux peut être inutile au dessin de l'image finale (grâce, par exemple, à l'accumulation de l'opacité de plusieurs surfaces combinées).

Nous proposons une technique simple pour trier les fragments d'un même pixel au moment de leur rasterisation, sans utiliser de liste chainée (et donc de pointeur). Puisque la liste des fragments pour un pixel est toujours triée, il est possible de déterminer, au moment de sa rasterisation, si un fragment contribuera ou pas à l'image finale, et de le rejetter le cas échéant. La technique combine les avantages de plusieurs approches existantes pour un niveau de performance similaire. Elle a l'unique avantage d'être très simple à coder : 24 lignes de GLSL.

Mots-clés :

A-buffer, temps-réel, transparence, hachage spatial 1 Introduction

Transparency and CSG operations have always been challenging for real-time rasterization. A major difficulty comes from the fact that the on-screen color of a pixel can only be determined after all surfaces influencing the pixel are known and sorted in depth-order. During rasterization, each triangle generates a number of fragments. Each fragment corresponds to a screen pixel -it is a small surface element potentially visible through this pixel. In a classical rasterizer only the fragment closest to the viewer is kept: The rasterizer blindly rejects all fragments which are further away than the current closest, using the depth-buffer algorithm.

Instead, algorithms dealing with transparency or CSG have to produce ordered lists of all the fragments falling into each pixel. This is typically performed in two stages: First, lists of fragments sorted by depth are gathered in each pixel. Second, rendering is performed by traversing the lists, either accumulating opacity and colors (transparency), or applying boolean operations to determine which fragment is visible (CSG). The data-structure is recreated at every frame, and therefore has to be extremely efficient and integrate well with the rasterizer. A large body of work has been dedicated to this problem (see Section 2).

Contributions This paper introduces a novel approach which is both efficient and simple to implement in available graphics APIs. In particular, our approach has the following advantages:

• It requires only a single rasterization pass to both store and sort all fragments, without depth-complexity limitations. • The memory overhead is low, 8 bits per stored fragment in our implementation.

• The implementation is simple, the entire fragment insertion and sort code fitting in less than 24 lines of GLSL (see also Figure 1). Access requires a simple loop without dependent memory reads. • Our approach enables a conservative early rejection of fragments that will have no or little influence on the final image. This is made possible by recent advances in spatial hashing. Our work introduces a hashing mechanism for depth-sorted fragment lists based on coherent hashing [START_REF] García | Coherent parallel hashing[END_REF]. We exploit properties of this hashing scheme that make it uniquely well suited to the problem of storing and sorting per-pixel fragments.

While the main focus of the paper is on transparency, our technique may be applied to other applications such as CSG and voxelization (in the manner of layered depth images [START_REF] Shade | Layered depth images[END_REF]). Our analysis shows that performance is similar to state of the art approaches with unique advantages: sort and store in a single geometry pass, conservative early culling, simplicity of implementation, arbitrary depth complexity and memory-performance tradeoff.

Related work

Most techniques for fragment accumulation implement a form of A-buffer [START_REF] Carpenter | The A-buffer, an antialiased hidden surface method[END_REF]. The A-buffer stores in each pixel the list of fragments that cover that pixel. The fragments are sorted by depth and the size of the list is called the depth-complexity. While the A-buffer is primarily designed for anti-aliasing it inspired many similar approaches for transparency, CSG and voxelization.

We review below the approaches for A-buffer construction most related to ours. Since our work is based on spatial hashing, we recall the necessary background in Section 2.5.

Depth peeling techniques

Depth-peeling [START_REF] Mammen | Transparency and antialiasing algorithms implemented with the virtual pixel maps technique[END_REF][START_REF] Everitt | Interactive order-independent transparency[END_REF][START_REF] Hable | Blister: GPU-based rendering of boolean combinations of free-form triangulated shapes[END_REF] is a multi-pass algorithm leveraging the Z-buffer to process all fragments in depth-order. It mimics selection-sort: after a first rasterization of the scene, the Z-buffer contains only the visible fragments, closest to the viewer. A second rasterization is performed, keeping the closest fragments which are farther than the Z-buffer of the previous iteration. This extracts a second layer of fragments. This process is iterated until no fragment remains.

The obvious drawback of the algorithm is that it requires as many rasterization passes as the depth complexity of the scene. An advantage, however, is that fragments are processed in front-to-back order without having to store all the fragments at once. This affords for early culling of fragments that are hidden, for instance, by opacity accumulation.

Per-pixel fixed-size arrays

A number of approaches pre-allocate fixed-sized buffers per pixel. A per-pixel counter is maintained to accumulate fragments in the buffers. This offers good performance, but imposes a hard limit on the per-pixel depth-complexity for a single rasterization while requiring large amounts of memory. Crassin [START_REF] Crassin | Fast and accurate single-pass A-buffer using OpenGL 4.0+[END_REF] describes an efficient implementation of this approach on mainstream GPUs.

Stencil routed A-buffer [START_REF] Myers | Stencil routed A-buffer[END_REF] captures up to K fragments per pixel by rendering to a multisampled buffer with multisampling antialiasing disabled. A stencil mask is applied to write only one sample at a time. This approach is limited by the maximum number of MSAA samples.

Several approaches propose a mechanism to deal with buffer overflows: For instance merging fragments with closest z values [START_REF] Jouppi | Z 3 : an economical hardware technique for high-quality antialiasing and transparency[END_REF] or relying on stochastic blending [START_REF] Enderton | Stochastic transparency[END_REF]. This is possible for effects which support approximations such as transparency but not for CSG or voxelization. Our algorithm is compatible with such approximations, but can also improve them since the sorted list of fragments gathered so far is available at all times. In particular, our algorithm can keep fragments up to a given opacity threshold, rather than summarily rejecting supernumerary ones.

FreePipe [START_REF] Liu | Freepipe: a programmable parallel rendering architecture for efficient multi-fragment effects[END_REF] allocates fixed size arrays and uses an atomicMin to sort the fragments as they are inserted. If K is the size of the per-pixel buffers, then only the K closest fragments are available after a single pass. Our approach also provides sorting in a single pass but without limiting the number of fragments per-pixel.

Per-pixel linked lists

The original A-buffer proposed an implementation based on per-pixel linked lists. A number of recent work describe how to implement this mechanism on modern rasterization hardware.

Inria

These techniques proceed in two stages: First, unordered lists of fragments are gathered in a geometry pass. Then the lists are sorted.

Yang et al. [START_REF] Yang | Real-time concurrent linked list construction on the gpu[END_REF] describe a linked-list A-buffer implementation in GLSL. A first buffer stores a list-head pointer per pixel. A second buffer stores fragment records (next-pointer, depth and data). When a fragment is inserted, a global counter is incremented to allocate the next fragment record. The new record is then linked to the pixel list by updating the head pointer atomically. When the geometry pass is complete, the linked lists are sorted. This technique is improved by allocating small blocks of fragments instead of single elements, thereby reducing contention on the allocation counter at the expense of additional memory usage [START_REF] Crassin | OpenGL 4.0+ A-buffer v2.0: Linked lists of fragment pages[END_REF].

A major advantage of these approaches is their simplicity of implementation, and their versatility: There is no other constraint on depth complexity than the maximum global number of fragments. There is a unique geometry pass, although a sort pass is required. In particular, fragments having no impact on rendering can only be detected after the sort. The per-fragment pointers incur a memory overhead.

Our approach brings similar advantages but in a single geometry pass (no additional sort), and without pointers. It enables early rejection of fragments. We further compare with these approaches in Section 6.2.

Count, gather, sort approaches

Another category of work relies on parallel stream processing techniques to allocate just enough memory for storing all the fragments seen from a given viewpoint. A first geometry pass counts the fragments covering each pixel. The obtained set of counters is then converted to prefix sums. These sums indicate where to start storing the fragments of each pixel. A fragment buffer of the required size is allocated and a second geometry pass is used to fill it. An optimized implementation of this technique is described in [START_REF] Maule | Memory-efficient order-independent transparency with dynamic fragment buffer[END_REF].

This approach gracefully adapts to large variations in the total number of fragments. It requires however two geometry passes as well as an efficient parallel scan implementation (prefix sum). All fragments have to be gathered, disabling conservative early culling.

Background in parallel hashing

Since our technique is based on parallel spatial hashing we recall here the necessary background to present our work. For the sake of brevity we do not describe spatial hashing techniques outside the scope of this paper [START_REF] Lefebvre | Perfect spatial hashing[END_REF][START_REF] Alcantara | Real-time parallel hashing on the GPU[END_REF] but focus on the technique of García et al. [START_REF] García | Coherent parallel hashing[END_REF] which inspired our technique.

A brief introduction to coherent hashing

Notations Coherent spatial hashing takes as input a set of integer keys K ∈ U , where U is called the universe and is typically the set of 32 bits integer coordinates. Each key is associated with a data record, for instance a 32 bits RGBA color. K is the set of defined keys while U \ K is the set of empty keys. Generally |K| |U |. The hashing scheme encodes K in a hash table H and afford for fast random queries to discover whether a key k ∈ U belongs to K and retrieve its associated data.

The table H stores both keys and data packed together in integer records. Let key(H[x]) and data(H[x]) be respectively the key and the data record stored in H[x]. Storing the keys alongside data is mandatory, and incurs a significant overhead: This doubles the required size of the hash table. In addition, not all entries of H are used: The load factor of a hash table is measured by the ratio |K| |H| . Fortunately, coherent hashing reaches load factors of up to 0.99.

Coherent hashing is based on an open addressing mechanism:

Each key k is associated with a sequence of possible insertion locations (h i (k)) i≥0 . If the key k is inserted at index h i (k) in the table, for some i ≥ 0, we say that the age of that key is i and define age(H[h i (k)]) = i. Coherent hashing keeps track of the age of the keys inserted in the hash during construction.

Parallel construction of the hash table

The keys and their associated data are inserted in the table in parallel. The thread responsible for key k explores its sequence of insertion locations. At each location two cases occur: 1) the location is empty in which case the key is immediately inserted, or 2) the location is occupied by another key q, in which case the algorithm checks whether to evict the key. If q is evicted then it is exchanged with k in the table and q becomes the key to be inserted: the thread resumes with insertion until finding an empty location in the table. Coherent hashing uses the Robin Hood strategy, which consists in evicting keys with a smaller age than the current key. This reduces the maximum age of the keys across the entire table to O(log N ) for a table storing N keys [START_REF] Celis | Robin hood hashing (preliminary report)[END_REF], thereby ensuring that all inserted keys can be accessed quickly.

Since many threads compete to insert keys, care must be taken to preserve the consistency of the hash table. Coherent hashing relies on an atomic instruction to both test for eviction and exchange the keys if needed. To achieve this, age, key and data are packed together into a 64 bits integer record. The age is stored in the most significant bits (MSB), so that an atomicMax operation can be used to both check for eviction and store the record. In this setting, the age is usually small and requires only 3-4 bits from the key and data. The full insertion code of García et l.3 MAX AGE is an upper bound on the maximal age of a key. l.5 Age, key and data are packed into an integer record. Age is in the MSB. l. [START_REF] Enderton | Stochastic transparency[END_REF] The hash relies on an atomicMax operation to automatically evict keys with a smaller age than the current key. l.7 This tests whether the atomicMax inserted the key. l.8-10 If a key was evicted it becomes the current key. KEY and AGE extract information from the packed integer record.

Coherent hash function The code described above is generic and could use any hash function.

Coherent hashing is obtained by generating the sequence of insertion locations using the following Inria hash function:

h i (k) = k + o i mod |H|
where o i is a predefined random offset that does not depend on k. If two neighboring threads process neighboring keys they access coalesced memory addresses at each step i. This greatly improves memory access efficiency during construction and access.

As we explain later this hash function has important properties for us, besides its improved memory access pattern.

Hash access Accessing a key amounts to following its sequence of locations until either it is found or MAX AGE is reached: access ( uint key )

for

( i = 0; i < MAX_AGE ; i ++ ) { h_k_i = hash ( key , i ) ; if ( KEY ( H [ h_k_i ] ) == key ) { return DATA ( H [ h_k_i ] ) ; } } return null ;

Overview

Our approach has two main stages: A gather+sort stage, where the geometry is rasterized, fragments are stored and sorted. A render stage, where the image is produced on screen from stored, depth-sorted fragments.

The critical stage is the gather+sort. We rasterize the scene using a specific fragment shader that shades each input fragment and inserts it in the hash table: In this shader, v Pos is the transformed and projected surface point, shadeFragment() computes the fragment color and opacity as a 32 bits integer, pos are the surface point coordinates remapped in [0, 1] 3 . V and MAX DEPTH define the screen resolution and depth resolution. Typical values are V = 1024 (for simplicity we assume a square viewport), MAX DEPTH = 2 24 . Note that we reverse depth so that large values correspond to fragments closer to the view. The reason for this will become apparent later.

During this pass nothing appears on screen: All fragments are discarded from the frame buffer after being inserted in the hash table.

The second pass is the render stage. A screen-covering quad is rendered with a shader accessing the sorted list of fragments for each pixel and computing the final color: void main () vec2 pos = ( u_Pos . xy * 0.5 + 0.5 ) ; uvec2 screen_pos = uvec2 ( pos . xy * V ) ; renderFragList ( screen_pos ) ;

These two shaders produce the final image. No other steps are required. Section 4 discusses how to implement insert using a variation of coherent hashing-this is the entire gather+sort. Section 4.1.4 explains how to retrieve data during the render stage (renderFragList). Section 5 presents a modification of the algorithm for early-culling of fragments.

The HA-Buffer

This section describes how we modify the coherent hash to obtain an efficient mechanism for simultaneously storing and sorting fragments. There are different aspects to our approach:

• We modify the definition of the keys to exploit collisions and form a list of fragments in each pixel (Section 4.1).

• We show that we can avoid storing the keys together with the data, resulting in significant memory savings (Section 4.1.1). • From the above we deduce a major property of the hash which lets us sort the keys as we insert them.

Hashing fragments

In order to insert a fragment in the hash table, we need to compute a key from information that describes the fragment (screen coordinates, depth, color). García et al. store and access voxels in the hash table by computing a key from their x, y, z coordinates, for instance with linear addressing. The case of fragments is different: We seek to store and retrieve efficiently a list of fragments for each screen pixel.

We propose to exploit collisions to this end: the idea is to form keys (to be hashed) from the screen coordinates only, ignoring the fragment depth. As a consequence, all fragments covering screen coordinates p collide with each other; they all appear along the access sequence (h i (p)) i≥0 . This is illustrated Figure 1.

This automatically produces lists of fragments. However, the downside is that the maximum age is at least the depth complexity of the pixel. To avoid incurring a large penalty in all pixels we implement the maximum age table mechanism described by García et al.: A 2D table covers the screen and records locally the maximum age reached by inserted fragments. This table has a lower resolution than the screen, so as to avoid divergence in neighboring threads during access: one cell of the maximum age table typically covers a block of 8 × 8 pixels.

This first step already provides hashing for fragment lists. However, the main benefits of using the hash are yet to be described.

Reducing memory overheads

When hashing with screen coordinates, fragments of a same pixel will collide with each other as well as with fragments from other pixels. We therefore still need to store the keys alongside the data, so as to remember by which pixel the fragment was produced.

Inria

This is a major overhead, typically 24 bits per fragment for encoding screen coordinates. This overhead is similar to the pointers used by linked-list approaches, making our approach less beneficial from a memory standpoint. Fortunately, we can improve this.

García et al. essentially considers the case where |H| |U |. However, our scenario is very different. The set of keys is limited to the screen coordinates. Assuming a full-screen rendering, we can expect at least one fragment per screen pixel. As a consequence the size of the hash table is expected to be larger than the screen:

|H| > |U |.
This has an important property: Consider the set of keys that can collide at age i in location for

x ∈ [0, |H| -1] of the table: {k ∈ U | k + o i mod |H| = x}.
( i = 0; i < MAX_AGE ; i ++ ) { x = h ( p , i ) ; if ( AGE ( H [ x ] ) == i ) { return DATA ( H [ x ] ) ; } } return null ;
Notation For a fragment covering pixel x at depth z inserted in the hash table at index h a (x) we note (a, z) x the integer record stored in H[h a (x)]. Only the age, depth and data are stored. For clarity we indicate the source pixel in subscript and ignore the data record. By the age equivalence property, we have:

H[h a (x)] = (a, z) q ⇒ q = x
We describe next a consequence of the age equivalence property which greatly facilitates sorting of the fragments in depth order.

Sorting fragments by depth

Recall that collisions are resolved with an atomicMax. This implements the Robin Hood strategy provided that the age is stored in the most significant bits of the integer records. From now on, we also assume that the depth immediately follows the age in the integer record: (a, z) x = MSB a z ...

LSB

. A natural question is whether this mechanism is also sorting fragments in decreasing depthorder, as they are inserted in the hash. It turns out this is not a trivial question. First, key ages and depths are uncorrelated: a key of any depth can be at any age. Second, many threads are working simultaneously, and some fragments may not be compared with each others due to evictions. We now prove that this in fact works as expected, i.e., that we obtain depth-sorted lists as a by-product of using coherent hashing.

Proof We start with a monotonicity property: " The hash table cells are exclusively updated with the atomicMax function. This implies that the value of a cell of H can only increase."

We then prove Proposition A = " At any time, if a cell of the hash table contains (a, z)

x then ∀b ∈ [0..a -1], H[h b (x)] ≥ (b, z) x .":
Assume that cell h a (x) of H contains (a, z) x in H. If a = 0 then proposition A is vacuously true. If a > 0 then consider b ∈ [0..a -1]. The fragment currently stored in h a (x) has been tested for insertion in all indexes h b (x): a record (b, z) x was compared to the value at index p = h b (x). Since it was not inserted at, or was evicted from that position, then at that time the integer record stored at index p was larger than (b, z) x . Because of the monotonicity property stated above, H[p] still contains a value larger than (b, z) x , which proves proposition A.

Finally, to prove our claim, assume that the table stores two fragments covering the same pixel x with depth z at index h a (x) and w at index h b (x), such that b < a. From proposition A it comes: (b, w

) x = H[h b (x)] ≥ (b, z) x , which implies w ≥ z.
Simply by encoding the depth after the age in the integer records, without any change to the insertion algorithm, we obtain lists of fragments ordered by decreasing depth. To obtain a front-to-back ordering we reverse the depth values so that largest values are closest to the viewer.

Discussion From the point of view of sorting, the algorithm exhibits a complexity of Ω(k 2 ) for k fragments in a pixel list. As k increases we can expect this to become inefficient. It is interesting, however, to note that this bears similarity with the sorting scheme of FreePipe [START_REF] Liu | Freepipe: a programmable parallel rendering architecture for efficient multi-fragment effects[END_REF] and the depth-peeling technique. Depth peeling, although very different in nature, also exhibits such a quadratic behavior for front-to-back enumeration of the fragments covering a pixel.

However, in our approach the unique geometry pass does not only sort but simultaneously hashes the fragment for storing them in memory. Therefore, only considering the sort complexity would be misleading: The cost is amortized by the simultaneous storage of the fragments.

In addition, the resulting algorithm maps very well to embarrassingly parallel architectures: Each thread independently inserts keys through a simple atomic operation, and thanks to the coherent scheme memory is accessed in a coalesced fashion. Elaborate gather and sort techniques tend to make such parallelism more difficult to achieve, to implement and to optimize.

Gather+sort algorithm

The complete pseudo-code for the gather+sort mechanism is: We store the fragment records in 64 bits, from MSB to LSB: 8 bits for the age, 24 bits for the depth, 32 bits for the data (RGBA color). PACK forms the 64 bit integer record. H, V and M are respectively the hash table size, the screen size and the max age table size. Note that the stored age starts at 1: We use 0 as a special value to indicate an empty cell in the hash (the hash table is cleared at the start of every frame).

This algorithm is implemented directly in GLSL, using the NV shader buffer store extension. The only complication on current hardware is the 64 bits atomicMax which is not yet available (it will be available on the GK110 NVidia GPU). We emulate it with a 64 bits atomicCAS (Compare And Swap). The resulting GLSL code is only 24 lines of code.

32 bits version Some applications do not require a data record. This is for instance the case of voxelization, where only the fragment depth is necessary. In this case a native 32 bits atomicMax can be used, resulting in an increased performance (Section 6.2).

Insertion failures

The insertion algorithm may fail in two situations: when the age of an inserted key does not fit on 8 bits, or when the hash table is full.

The maximum reachable value for the age of a key depends on the worst depth-complexity and the load factor (the eviction mechanism of Robin Hood hashing is activated more often at higher load factors). A complex scene could require a depth complexity larger than 256, in which case the age can be encoded on more bits. However, as we describe in Section 5 some applications such as transparency allow for an early-cull mechanism which reduces the need for a large maximum age.

A full hash table would also trigger insertion failures. In addition, performance varies with the hash table load factor (Section 6.2 analyzes this trade-off). We therefore dynamically adapt the size of the hash table. It is initially set to the screen resolution, which provides a reasonable estimate for viewing an object centered on screen. We monitor the number of fragments stored in the hash and dynamically reallocate when reaching a high density (80%). Similarly, we decrease the hash table size if a low density is reached (30%). We always increase/decrease the table size by a the same ratio (25%). To count fragments efficiently, we exploit the clear mechanism. The clear is performed by rendering a quad covering the hash table. An occlusion query monitors the number of drawn fragments during the clear. If the fragment being cleared was previously empty, it is discarded. Otherwise, it is drawn and will be counted by the occlusion query. This counts fragments with a one-frame latency. This mechanism has a very small impact on performance.

Offsets for the hash function Recall that the sequence of insertion positions for key k is (h i (k) = (k +o i ) mod |H|) i≥0 where (o i ) i≥0 is a sequence of offsets in the hash table: ∀i, o i ∈ H. The good properties of Robin Hood hashing have been analyzed in a probabilistic setting, as most other hashing schemes. Therefore, in order to ensure a well behaved insertion in table H at all times, it is recommended to regenerate the offset sequence o at random when the hash table is cleared, i.e., before each frame.

Since inserted fragments are sorted in correct front-to-back order in the hash table, our algorithm offers a unique opportunity for conservative early culling. That is, we can check during insertion whether a fragment is still useful for the final rendering based on the sorted list of previously inserted fragments. If we determine that the fragment is no longer useful we stop trying to insert it. The decision is conservative: We may insert too many fragments, but we will never improperly reject them. This is in effect very similar to the early depth cull of the Z-buffer: The efficiency of the mechanism depends on the order with which surfaces will be sent for rasterization. However, in complex scenes it results in a significant savings in the amount of fragments that need to be stored, as illustrated Figure 2.

The GLSL code for early culling follows. It is a simple modification of the earlier insertion algorithm that we specialize here for opacity accumulation. Other applications for which a conservative culling criterion can be devised should benefit from early culling as well. (translucent shadow map + view) at varying opacity. Early-culling directly benefits the frame rate even in this challenging case (opacity threshold: 99%). The insets show the color-coded number of culled fragments. The total number of fragments without culling is 15 420 969 -this is roughly doubled when shadowing is enabled.

This issue only affects the 64 bits version of our code. We therefore use the 32 bits version to estimate the overall performance penalty. Figure 5 measures performance of our reference 32 bits implementation as well as a 32 bits implementation modified to include similar penalties as the 64 bits version: The atomicMax is emulated with atomicCAS and the same additional atomic memory fetch is added. We find that the ratio between both is 1.37 on average. In the remainder of the paper we report the performance measured on our 64 bit implementation and in brackets the projected performance.

Transparency

Transparency effects directly benefit from the HA-buffer. Once fragments have been inserted in the hash table, it suffices to traverse the lists and blend the fragments in front-to-back order. Let F i rgb be the color of fragment i and F i α its opacity:

C i+1 rgb = C i rgb + F i rgb (1 -C i α ) C i+1 α = F i α (1 -C i α )
with C 0 α = 0, C 0 rgb = 0. A last opaque fragment with any chosen background color F bkg rgb is implicitly added. The final pixel color is C K+1 rgb where K the number of fragments for the pixel. We optionally add a translucent media effect inside surfaces, tracking in/out pairs of fragments. This effect is visible Figure 3.

Performance vs load factor Our approach offers a tradeoff between performance and memory: Performance varies depending on the hash table load factor. Figure 5 illustrates this tradeoff. The object of Figure 7 is displayed and rotates on screen, with a flat shading (simple alpha blending of constant color fragments). The viewpoint is chosen to fill the hash up to about 95%, and then the hash table size is gradually increased. As can be seen performance slowly decreases as the load factor increases, with an abrupt drop after 80%: the high load factor decreases the probability to find an empty cell in the hash table during fragment insertion.

We next analyze the performance and memory requirements with respect to other techniques. Analysis: Memory We consider a screen of size V 2 , fragments of 56 bits (depth:24 RGBA:32), and a total number of fragments F . Two existing techniques offer compact storage of the fragments (see Section 2): count-gather-sort and linked-lists. The count-gather-sort approaches only store a per-pixel pointer, but require two geometry passes (count+gather) and a parallel scan on all pixels to complete. The required storage is 32V 2 + 56F bits.

The per-pixel linked-lists store one head pointer per-pixel, as well as a pointer per fragment. The required storage is 32V 2 + (56 + 32)F bits. Paged linked-lists store one head pointer perpixel, a per pixel counter (8 bits) as well as several pages for the fragments. Each page contains K fragments and a pointer. The technique requires at least (32 + 8)V 2 + F/K(56K + 32) bits and at most 56(K -1)V 2 more.

Using our HA-buffer the storage is the hash table and the maximum age table. Each entry of the max age table covers a block of 8 × 8 pixels. The age is encoded on 8 bits (could be less in many situations), resulting in 64 bits per fragment in the hash. The storage is a function of the load factor L: V 2 /8 + 64F/L bits.

Table 1 summarizes memory occupancy for a variety of settings. As can be seen our technique offers compact storage at high load factors. Of course, this impacts performance and we recommend 80% for memory intensive applications, and to target 50% for best performance. Even at this lower load factor the memory consumption remains comparable to that of other schemes thanks to the low overhead per fragment.

Inria

We next discuss performance of these methods.

Analysis: Speed We compare the speed of our technique to several available A-Buffer implementations. In each case we use the original code from the authors, only modified to use a simple flat shading with alpha blending. We measure performance of our 64 bits implementation (without early culling). Results are summarized in table 2 for rendering the view shown in the inset.

Our objective with this analysis is to verify that our technique, despite the simplicity of its algorithm, offers similar levels of performance as the best available approaches. It is worth mentioning that performance measures for A-buffers depend on many factors: Balance between geometric complexity and number of fragments, opacity, early culling, etc. Hardware design also has a huge impact: For instance in table 2 the performance of paged perpixel linked lists increases by 57% when using textures instead of global memory, resulting in the fastest approach. On a Fermi architecture both were equivalent, and slower than our approach.

The table does not include the count-gather-sort of [START_REF] Maule | Memory-efficient order-independent transparency with dynamic fragment buffer[END_REF] but this paper reports performance in the range of paged per-pixel linked lists (see Figure 7 of [START_REF] Maule | Memory-efficient order-independent transparency with dynamic fragment buffer[END_REF]) on a benchmark similar to our test scene. We however expect this approach to outperform all others when the fragment complexity is much higher than the geometric complexity. The strengths of this approach are its tight memory allocation and its efficient construction of contiguous lists of pixels via a parallel scan. Its weaknesses are that it always requires two full-resolution geometry passes -making it less efficient when geometric complexity increases for a same number of fragments -and that it cannot reject fragments based on accumulation effects. In contrast we perform all computations in a single geometry pass, and support rejection of fragments based on the current sorted lists of fragments.

The second fastest approach is paged per-pixel linked lists. As shown in table 2, using 64 bits the HA-buffer is slower than paged linked-lists using texture memory -however the expected performance boost of native atomicMax64 support gives an equivalent performance. We also modified the code of paged linked-list to store only depth in 32 bits textures and compared it to our 32 bits approach which does not suffer from penalties. In this scenario we also achieve a similar frame rate. In conclusion our approach exhibits similar performance than state of the art approaches, with similar memory requirements. Note however that there is a speed-memory tradeoff: At compact memory usage our performance is lower, but applications that can afford more memory directly see a performance increase.

Inria

Shadows The HA-buffer can be used for translucent shadow effects: A HA-buffer is rendered from the light view and used as a shadow map. Since it contains all fragment colors and opacities ordered by depth, translucent effects may be applied. When the HA-buffer for the view is rendered, the HA-buffer shadow map is used for shading fragments. The HA-buffer shadow map is only updated when the light moves. Figures 2, 7 and 6 illustrates this.

Early culling Figure 2 illustrates early culling on a complex object. It is worth noting that the fill-rate remains constant: This outlines the efficiency of early culling since the FPS proportionally increases with the reduction in number of fragments.

Analysis: Stress test Figure 8 is a stress-test of our approach: We render random quads of a same size at various depth in orthographic projection. The quads are then slowly brought together, until they all exactly stack up. The number of fragments remains constant through this process, but the depth complexity in pixels not only increases but becomes more and more uniform. As expected, once all pixels have a same large amount of fragments performance degrades. One reason is that parallelism cannot hide the cost of long lists by processing many small lists simultaneously: The work load is the same everywhere. The other reason is that the sort complexity is in Ω(k 2 ) (Section 4.1.2) and therefore the cost increases significantly with larger number of fragments. Of course, such a catastrophic case is highly unlikely in practice, and as shown Figure 8 using early culling alleviates the issue even in the worst cases.

Constructive Solid Geometry

We implement a simple CSG renderer by walking along the ordered lists of fragments, keeping track of the entry/exit points of the primitives. A bit field with one bit per primitive indicates whether we are inside/outside. As we walk along the list of fragments, the bit field and the boolean expression are updated in order to check whether the current fragment is on the boundary of the desired CSG model. A sample CSG rendering is shown Figure 9. The main challenge with CSG is the high depth complexity and the large number of fragments per pixel involved in even simple shapes. Our implementation is quite naive and different approaches and optimizations should be experimented with, inspired by state of the art work [START_REF] Kirsch | Rendering techniques for hardware-accelerated imagebased CSG[END_REF][START_REF] Hable | CST: Constructing solid trimming for rendering BReps and CSG[END_REF].

Limitations, future work

As we have seen the main limitation of our approach is related to the increase in complexity when long fragment lists appear in all pixels. This is fortunately a rare situation, and it can be alleviated through the use of early culling in transparency applications.

For future work we would like to consider conservative tests regarding CSG. It seems interesting to revisit the works that have been done in the context of depth peeling and see how these would apply to our approach.

Conclusion

When looking at the implementation, the hash mechanism we describe seems a simple modification of coherent hashing to deal with fragment lists. However, the important contribution of this work is to show that this mechanism offers unique properties regarding storage and sorting of fragments. It would be easy to miss these if hashing was naively applied to an A-buffer construction.

The good performance of our approach is due to the simplicity of its implementation and the coherence properties of the hash function. We however do not think this is the fastest approach Inria for A-buffer. It will be especially efficient on scenes of moderate fragment complexity with high geometric complexity. Besides performances, we would like to re-emphasize the properties of our approach: A unique geometry pass is immediately followed by the rendering pass, both utilizing a coherent memory access pattern. No separate sorting is necessary and conservative early culling is simple to add. There are no pointers and the access loop does not even involve dependent reads in memory. The construction and access code fit within 50 lines of GLSL code. Together, these properties cover most of the advantages of prior methods. These benefits make the approach versatile: it takes little time to implement a variety of applications or mix several viewpoints with sorted fragments for shadowing effects.

We hope our approach will facilitate the wide adoption of transparency and CSG as standard effects in real-time rendering engines.

Figure 1 :

 1 Figure 1: Left: The A-buffer stores lists of fragments in each pixel. Instead of lists, we use a single hash table wherein each pixel p is assigned a fixed permutation i = p + o i of locations where to store the fragments covering p. A simple mechanism is used to detect collisions and skip over them. Middle: Transparent rendering with translucent shadows (118 FPS, 1024 2 viewport). Right: The complete algorithm for inserting and simultaneously sorting fragments, called from the fragment shader during rasterization.

  void main () uint data = shadeFragment () ; vec3 prj = v_Pos . xyz / v_Pos . w ; vec3 pos = vec3 ( prj . xy * 0.5 + 0.5 , 1 -( v_Pos . z + ZNear ) /( ZFar + ZNear ) ) ; uvec2 screen = uvec2 ( pos . xy * V ) ; uint depth = pos . z * MAX_DEPTH ; insert ( screen , depth , data ) ; discard ;

Figure 2 :

 2 Figure 2: The hairball model rendered at 1024 2 resolution and two HA-buffers of size 4096 2(translucent shadow map + view) at varying opacity. Early-culling directly benefits the frame rate even in this challenging case (opacity threshold: 99%). The insets show the color-coded number of culled fragments. The total number of fragments without culling is 15 420 969 -this is roughly doubled when shadowing is enabled.

Figure 3 :

 3 Figure 3: Horses and chariot rendered with our HA-buffer using a GeForce GTX 680 and FXAA. 891929 frag., load: 54.4% (12.6 MiB) max depth: 22, FPS: 216.7.

Figure 4 :

 4 Figure 4: Technical models and intricate shapes rendered with our HA-buffer using a GeForce GTX 680 and FXAA. Left: Robotic hand model, 1487942 frag., load: 58.1% (19.6 MiB) max depth: 20, FPS: 91.1 Middle: Joystick model, 539236 frag., load: 51.4% (8.06MiB) max depth: 16, FPS: 348 Right: Dendrite model, 855296 frag., load: 52.2% (12.6 MiB) max depth: 16, FPS: 202.7

Figure 5 :

 5 Figure 5: Millions of fragments per second as a function of hash table load factor. Performance is averaged over multiple runs.

Figure 6 :

 6 Figure 6: Two viewpoints of a glass sculpture rendered with translucent HA-buffer shadows. (Fixed light, V = 1024 2 H = 1450 2 , FXAA enabled). Left: 1021046 frag., load: 48.6% (16.1 MiB) max depth: 16, FPS: 118.0 (shadows on). Right: 1001412 frag., load: 47.6% (16.1 MiB) max depth: 18, FPS: 122.5 (shadows on).

Figure 7 :

 7 Figure 7: Left: The shadow projected by the baby dinosaur on the back-lit egg shell is encoded in a second HA-buffer, used as a shadow map. This back-lit shadow is only possible thanks to the translucent shadow map. View: 744208 frag., load: 71% (8 MiB) max depth: 24, FPS: 119.5 (shadows on). Right: Front view with different lighting and opacity. View: 701822 frag., load: 67% (8 MiB) max depth: 20, FPS: 158 (shadows on). (Fixed light, V = 1024 2 , H = 1024 2 , FXAA enabled).

  al. is given below, followed by a short description:

	insert ( uvec3 key , uint data )
	age = 0;
	while ( ++ age < MAX_AGE ) {
	h_k_i	= hash ( key , age ) ;
	record = PACK ( age , key , data ) ;
	prev	= atomicMax ( & H [ h_k_i ] , record ) ;
	if ( record > prev ) {
	if ( EMPTY ( prev ) ) ) {
	return SUCCESS ;
	} else {
	key = KEY ( prev ) ;
	age = AGE ( prev ) ;
	}
	}
	}
	return FAILURE ;

  In our case this set contains a unique coordinate, which follows from |H| > |U |: the screen does not fold onto itself when mapping into H. This coordinate is the unique pixel coordinate that will send a fragment at age i into location x. It is easily computed from the age as p = xo i mod |H|. No other pixel can produce a fragment that maps to x at age i. We refer to this as the age equivalence property. As a consequence, a fragment inserted in the table belongs to pixel p if and only if age(H[h i (p)]) = i. Thus, we can safely remove from the hash table the pixel screen coordinates. The age is enough to identify the source pixels of the fragments. The code for accessing a fragment covering pixel p becomes: access ( uvec2 p )

Table 1 :

 1 Estimated memory usage (in MiB) of the different A-buffer techniques for a screen size V = 1024 2 and a varying number of fragments F .

	Technique	F = 0.85 M 1.5 M 10 M
	count-gather-sort	9.7	14.0	70.8
	linked-list	12.9	19.7 108.9
	paged, K = 4, lower bound	11.5	16.4	81.3
	paged, K = 4, upper bound	32.5	37.4 102.3
	HA-buffer, L = 90 %	7.2	12.7	84.8
	HA-buffer, L = 80 %	8.1	14.3	95.4
	HA-buffer, L = 50 %	13.0	22.9 152.6

Table 2 :

 2 Performance of different A-buffer techniques measured in our benchmark. The view contains 938174 fragments, and a maximum depth complexity of 10. The model has 871414 triangles (Stanford Dragon model). (*) For the 64 bits version we give in parenthesis projected performance without atomicMax64 emulation and compiler issue.

	Technique	Frame rate
	Dual depth peeling (8 passes)	140 FPS
	Depth peeling (16 passes)	82 FPS
	Stencil routed K-Buffer (16 samples)	226 FPS
	Paged linked lists (global)	157 FPS
	Paged linked lists (texture)	246 FPS
	HA-buffer 50% (64 bits)	187 FPS (257*)
	HA-buffer 80% (64 bits)	83 FPS (114*)
	Paged linked lists (32 bits, global)	180 FPS
	Paged linked lists (32 bits, texture)	290 FPS
	HA-buffer 50% (32 bits)	298 FPS
	HA-buffer 80% (32 bits)	181 FPS

RR n°8282

rue du Jardin Botanique CS20101 54603 Villers-lès-Nancy Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments This work was supported by ERC grant ShapeForge (StG-2012-307877). The authors thank NVIDIA for hardware donation.

In practice the set of defined keys is often cluttered due to spatial coherence. We found that in this situation the set of offsets o i can be made slightly more efficient. We use iterated farthest point sampling [START_REF] Teofilo | Clustering to minimize the maximum intercluster distance[END_REF] to generate offsets which are as far as possible from each others. In addition, the first offsets have a greater distance between them, making it more likely to find a free slot in the hash table at small age values. We seed this process randomly and generate different offsets at every start of the application -this is unfortunately too slow for regeneration at every frame. A drawback is that biases due to the structure in the offsets are more likely to occur. All renderings in the paper use such offsets, while synthetic test cases use random offsets to avoid bias in measurements.

Traversing the list of fragments covering a given pixel

The GLSL code for traversing the front-to-back list of fragments covering a pixel at position screen follows. Haddr and Maddr compute a linear address from the 2D coordinates in row-major order. If the condition in line 16 is verified, then proposition A (Section 4.1.2) guarantees that we have already seen all the fragments covering the query pixel.

Early cull mechanism

The algorithm as we described it so far gathers all fragments per-pixel. However, there are some cases were many fragments end up not being used in the final rendering. This a common situation when the accumulation of semi-transparent surfaces results in an almost opaque media: deeper fragments have no impact on the final image.

Few methods have the ability to efficiently detect these cases: In most other approaches the fragments are gathered first and sorted in a separate pass (see also Section 2). Depth peeling offers such a functionality but due to multiple geometry passes is much slower than our approach (see Section 6.2). Fixed-size arrays methods propose mechanisms to deal with overflows, relying on approximate techniques to merge fragments in excess. We could easily implement a similar mechanism, keeping only the K closest fragments. However, what we propose here is different: The most important part of the algorithm are lines 44-53, which are executed when a fragment f is not inserted. If the age is the same as the stored fragment s (line 44), then both fragments cover the same pixel. Fragment f was not inserted because it is located after s along the view ray. We therefore accumulate the opacity of s (lines 46-48) and test whether f will still be visible by comparing the accumulated opacity to a user-defined threshold.

The other change to the base algorithm are lines 32-38, executed after a fragment e has been evicted by the inserted fragment f . It is important to verify whether e covers the same pixel as f . If it does not, the opacity accumulator is reset (line 33) since it is no longer valid. Note that e may not be the first fragment for its pixel, but since we cannot know the accumulated opacity at e we simply set it to 0. This is conservative since opacity monotonically increases. If e and f cover the same pixel, we do not increment the age of e so that the opacity of f can be accumulated in the next iteration.

Results

In this section we further analyze the performance of our technique. We focus on transparency and only briefly present results on CSG to illustrate the versatility of our technique.

Test setup

For performance results we use a NVidia GeForce GTX 680 (Kepler GK104), NVidia driver 310.90 and OpenGL 4.3. The GK104 lacks native support for atomicMax64, contrary to the soon to be released GK110. We therefore emulate atomicMax64 with atomicCAS64. Unfortunately, on Kepler this triggered a known GLSL compiler issue regarding race conditions and atomics, preventing proper execution of our shader (the same shader runs correctly on a Fermi GPU). To circumvent this issue we have to insert a useless atomic memory fetch in the main loop of our shader. This has a negative impact on our performance but prevents the compiler from generating the problematic code. The squares have 0.5 opacity and are pulled together until exactly stacking, an unfavorable situation for our algorithm. The middle row shows the frame rate, and the bottom row the maximum age in the table. The red curve is our standard scheme. Note the large increase towards the end, as all pixels get the exact same workload. Fortunately, when enabling early culling (green curve) the negative effects are compensated by the large number of culled fragments (threshold for culling: 99% opacity).


Vladimir Reinharz

Yann Ponty
email: yann.ponty@lix.polytechnique.fr

J Ér Ôme Waldisp Ühl

J Ér Ôme

Waldisp Ühl

A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotides distribution

Keywords: RNA, secondary structure, design, weighted sampling, GC-content

Motivations:

The design of RNA sequences folding into predefined secondary structures is a milestone for many synthetic biology and gene therapy studies. Most of the current software uses similar local search strategies (i.e. a random seed is progressively adapted to acquire the desired folding properties) and more importantly do not allow the user to control explicitly the nucleotide distribution such as the GC-content in their sequences. However, the latter is an important criterion for large-scale applications as it could presumably be used to design sequences with better transcription rates and/or structural plasticity. Results: In this paper, we introduce IncaRNAtion, a novel algorithm to design RNA sequences folding into target secondary structures with a predefined nucleotide distribution. IncaRNAtion uses a global sampling approach and weighted sampling techniques. We show that our approach is fast (i.e. running time comparable or better than local search methods), seed-less (we remove the bias of the seed in local search heuristics), and successfully generates highquality sequences (i.e. thermodynamically stable) for any GC-content.

To complete this study, we develop an hybrid method combining our global sampling approach with local search strategies. Remarkably, our glocal methodology overcomes both local and global approaches for sampling sequences with a specific GC content and target structure.

INTRODUCTION

At the core of the emerging field of synthetic biology resides our capacity to design and re-engineer molecules with target functions. RNA molecules are well tailored for such applications. The ease to synthesize them (they are directly transcribed from DNA) and the broad diversity of catalytic and regulation functions they can perform enable to integrate de-novo logic circuits within living cells [START_REF] Rodrigo | De novo automated design of small RNA circuits for engineering synthetic riboregulation in living cells[END_REF] or re-program existing regulation mechanisms [START_REF] Chang | Synthetic RNA switches as a tool for temporal and spatial control over gene expression[END_REF]. Future advances and applications of these * to whom correspondence should be addressed techniques in gene-therapy studies will strongly rely on efficient computational methods to design and re-engineer RNA molecules.

Most of RNA functions are, at least partially, encoded by the three-dimensional molecular structures, which are themselves primarily determined by the secondary structures. The development of efficient algorithms for designing RNA sequences with pre-defined secondary structures is thus a milestone to enter the synthetic biology era. RNAinverse pioneered RNA secondary structure design algorithms. It has been developed and distributed with the Vienna RNA package [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF]. However, only posterior experimental studies revealed the potential and practical impact of these techniques. Thereby, during the last 6 years many improvements and variants of RNAinverse have been proposed. Conceptually, almost all of existing algorithms follow the same approach. First a seed sequence is selected, then a local search strategy is used to mutate the seed and find, in its vicinity, a sequence with desired folding properties. Using this strategy, INFO-RNA [START_REF] Busch | INFO-RNA-a fast approach to inverse RNA folding[END_REF], RNA-SSD [START_REF] Aguirre-Hernández | Computational RNA secondary structure design: empirical complexity and improved methods[END_REF] and NUPACK:Design [START_REF] Zadeh | Nucleic acid sequence design via efficient ensemble defect optimization[END_REF] significantly improved the performance of RNA secondary structure design algorithms. More recent research studies aimed to include more constraints in the selection criteria. RNAexinv focused on the design of sequences with enhanced thermodynamical and mutational robustness [START_REF] Avihoo | RNAexinv: An Extended Inverse RNA Folding from Shape and Physical Attributes to Sequences[END_REF], while Frnakenstein enables to design RNA with multiple target structures [START_REF] Lyngsø | Frnakenstein: multiple target inverse RNA folding[END_REF].

We recently introduced with RNA-ensign a novel paradigm for the search strategy of RNA secondary structure design algorithm [START_REF] Levin | A global sampling approach to designing and reengineering rna secondary structures[END_REF]. Instead of a local search approach, we proposed a global sampling strategy of the mutational landscape based on the RNAmutants algorithm [START_REF] Waldispühl | Efficient Algorithms for Probing the RNA Mutation Landscape[END_REF]. This methodology offered promising performances, but suffered from prohibitive runtime and memory consumption. Following our work, Garcia-Martin et al proposed RNAiFOLD [START_REF] Garcia-Martin | RNAiFold: A constraint programming algorithm for RNA inverse folding and molecular design[END_REF], an alternate methodology that uses constraint programming techniques to prune the mutational landscape. While also suffering from prohibitive running times, it is worth noting that this latter algorithm also proposes a seed-less approach to the RNA secondary structure design problem.

In this paper, we introduce IncaRNAtion, a RNA secondary structure design algorithm that benefits of our recent algorithmic advances [START_REF] Reinharz | A linear inside-outside algorithm for correcting sequencing errors in structured rna sequences[END_REF] to expand our original RNA-ensign algorithm [START_REF] Levin | A global sampling approach to designing and reengineering rna secondary structures[END_REF]. IncaRNAtion addresses previous limitations of RNA-ensign and offers new functionalities. First, while our previous program had a running time complexity of O(n 5 ), IncaRNAtion now runs in linear-time and space complexity, allowing it to demonstrate similar speeds as any local search algorithm. Next, IncaRNAtion is seed-less. Unlike RNA-ensign, it does not require a seed sequence to initiate its search. Finally, IncaRNAtion implements a novel algorithm based on weighted sampling techniques [START_REF] Bodini | Multi-dimensional Boltzmann Sampling of Languages[END_REF] that enables us to control, for the first time, explicitly the GC-content of the solution. This functionality is essential because wild-type sequences within living organisms often present medium or low GC-content, presumably to offer better transcription rates and/or structural plasticity. Previous programs do not allow to control this parameter and tend to output sequences having high GC-contents [START_REF] Lyngsø | Frnakenstein: multiple target inverse RNA folding[END_REF].

We demonstrate the performance of our algorithms on a set of real RNA structures extracted from the RNA STRAND database [START_REF] Andronescu | RNA STRAND: the RNA secondary structure and statistical analysis database[END_REF]. To complete this study, we develop an hybrid method combining our global sampling approach with local search strategies such as the one implemented in RNAinverse. Remarkably, our glocal methodology overcomes both local and global approaches for sampling sequences with a specific GC content and target structure.

METHODS

We introduce a probabilistic model for the design of RNA sequences with a specific GC-content and folding into a predefined secondary structure. For the sake of simplicity, we choose to base this proof-ofconcept implementation on a simplified free-energy function E(•), which only considers the contributions of stacked canonical basepairs. We show how a modification of the dynamic programming scheme used in RNAmutants allows for the sampling of good and diverse design candidates, in linear time and space complexities.

Definitions

A targeted secondary structure S * of length n is given as a noncrossing arc-annotated sequence, where S * i stands for the basepairing position of position i in S * if any (and, reciprocally, S * S * i = i), or -1 otherwise. In addition, let us denote by #gc(s) the number of occurrences of G and C in an RNA sequence s.

Simplified energy model

We use a simplified free-energy model which only includes additive contributions from stacking base-pairs. Using individual values from the Turner 2004 model (retrieved from the NNDB [START_REF] Turner | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF]). Given a candidate sequence s for a secondary structure S, the free-energy of any sequence s of length |S| is given by

E(s, S) = (i,j)→(i ,j )∈S stacking pairs E β s i s j →s i s j
where E β ab→a b is set to 0 if ab = ∅ (no base-pair to stack onto), the tabulated free-energy of stacking pairs (ab)/(a b ) in the Turner model if available, or β ∈ [0, ∞] for non-Watson-Crick/Wobble pairs (i.e. not in {GU, UG, CG, GC, AU or UA}). This latter parameter allows one to choose whether to simply penalize invalid base pairs (β > 0), or forbid them altogether (β = +∞). Positionspecific sequence constraints can also be enforced at this level (details omitted for the sake of clarity) by assigning to E a +∞ penalty (leading to a null probability) in the presence of a base incompatible with a user-specified constraint mask.

2.1.2 GC-weighted Boltzmann ensemble and distribution In order to counterbalance the documented tendency of sampling methods to generate GC-rich sequences [START_REF] Levin | A global sampling approach to designing and reengineering rna secondary structures[END_REF], we introduce a parameter x ∈ R + , whose value will influence the GCcontent of generated sequences. For any secondary structure S, the GC-weighted-Boltzmann factor of a sequence s is B

[x] S (s) such that B [x] S (s) = e -E(s,S) RT • x #gc(s) (1)
where R is the Boltzmann constant and T the temperature in Kelvin.

Summing the GC-weighted-Boltzmann factor over all possible sequences of a given length |S|, one obtains the GC-weighted partition function Z

[x]

S , from which one defines the GC-weighted Boltzmann probability P S (s) of each sequence s, respectively such that

Z [x] S = |s|=n B [x]
S (s)

and

P [x] S (s) = B [x] S (s) Z [x] S . (2)
2.2 Linear-time stochastic sampling algorithm for the GC-weighted-Boltzmann ensemble

Let us now describe a linear-time algorithm to sample sequences at random in the GC-weighted Boltzmann distribution. This algorithm follows the general principles of the recursive approach to random generation [START_REF] Wilf | A unified setting for sequencing, ranking, and selection algorithms for combinatorial objects[END_REF], pioneered in the context of RNA by the SFold algorithm [START_REF] Ding | A statistical sampling algorithm for RNA secondary structure prediction[END_REF]. The algorithm starts by precomputing the partition function restricted to each substructure occurring in the target structure, and then performs a series of recursive stochastic backtracks, using precomputed values to decide on the probability of each alternative.

2.2.1

Precomputing the GC-weighted partition function Firstly, a dynamic programming algorithm computes Z [a,b] N,S the GCweighted partition function (the dependency in x is omitted here for the sake of clarity) for a structure S, assuming its (previously chosen) flanking nucleotides are a and b respectively, either forming a closing base-pair (N = T) or not (N = F). Remark that the empty structure only supports the empty sequence, having energy 0, so one has

Z [a,b] T,ε = Z [a,b] F,ε = e -0/RT = 1. (3)
The general recursion scheme consists in three different terms, depending on the first position in S: Case 1. First position is unpaired (S = •S ):

Z [a,b] T,•S = Z [a,b] F,•S := a ∈B x #gc(a ) • Z [a ,b] F,S ; (4)
Case 2. First position is paired with last position (S = (S )), stacking onto a pre-existing exterior pair (N = T):

Z [a,b] T,(S ) := a ,b ∈B 2 x #gc(a .b ) • e -E β ab→a b RT • Z [a ,b ] T,S ; (5) Algorithm 1: SBx (a, b, N, S) r ←Random Z [a,b] N,S // Random real in [0, Z [a,b] N,S [ switch do case S = ε return ε; // Empty structure case S = • S // First position is unpaired for a ∈ B do r ← r -x #gc(a ) • Z [a ,b] F,S if r < 0 then return a .SBx(a , b, F, S ) case S = ( S ) and N = T // Extremities are involved in stacking base pair for (a , b ) ∈ B × B do r ← r -x #gc(a .b ) • e -E β ab→a b /RT • Z [a ,b ] T,S if r < 0 then return a .SBx(a , b , T, S ).b otherwise // First position is paired without a stacking pair // S = ( S ) S for (a , b ) ∈ B × B do r ← r -x #gc(a .b ) • e -E β ∅→a b RT • Z [a ,b ] F,S • Z [b ,b] T,S if r< 0 then return a .SBx (a , b , T, S ) .b .SBx (b , b, F, S )
Case 3. First position is involved in a base-pair (S = (S )S ), which is not stacking onto an exterior base-pair (N = F or S = ε):

Z [a,b] N,(S )S := a ,b ∈B 2 x #gc(a .b ) • e -E β ∅→a b RT • Z [a ,b ] T,S • Z [b ,b] F,S .
(6) Remark that the number of combinations of a, b and N remains bounded by a constant, thus the complexity of computing Z [a,b] N,S mainly depends on the values taken by S upon subsequent recursive calls. Such values are entirely determined by S at any given step of the recursion, and their dependency can be summarized in a tree having Θ(|S|). Therefore, the computation of Z [a,b] N,S * requires Θ(n) time and space using dynamic-programming.

Stochastic backtrack

Once the GC-weighted partition functions have been computed and memorized, a stochastic backtrack starts from the target structure S * with any exterior bases [a, b] and no nesting base-pair, corresponding to a call SBx (∅, ∅, F, S * ) to Algorithm 1. At each step, a suitable assignment for one or several positions is chosen, using probabilities derived from the precomputation, as illustrated by Figure 1. One or several recursive calls over the appropriate substructures are then performed. On each recursive call, the algorithm assigns at least one nucleotide to a -previously unassigned -position. Moreover, the number of executions of each loops is bounded by a constant. Consequently, the complexity of Algorithm 1 is in Θ(n) time and space.

Self-adaptive sampling strategy

Let us remind that our goal is to produce a set of sequences whose GC-content matches a prescribed value gc. An absolute tolerance κ may be allowed, so that the GC-content of any valid sequence must fall in [gc -κ, gc + κ]. Since sequences of arbitrary GC-content may be generated by Algorithm 1, we use a rejection-based approach [START_REF] Bodini | Multi-dimensional Boltzmann Sampling of Languages[END_REF], previously adapted by the authors in a similar context (Waldispühl and Ponty, 2011). This gives an algorithm which generates k valid sequences in expected time

Θ(k • n √ n) when κ = 0 (or Θ(k • n) when κ is a positive constant) and memory in Θ(k • n).
A complete analysis of the rejection process can be found in an earlier contribution [START_REF] Waldispühl | An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure[END_REF], but let us briefly outline the approach, and the main arguments used to establish its complexity.

As summarized by Figure 2, our adaptive sampling approach simply generates sets of sequences by repeatedly running the stochastic backtrack algorithm. The average GC-content induced by the current value of the x parameter, can then be adequately estimated from the sample, or computed exactly using recent algorithmic advances [START_REF] Ponty | A Combinatorial Framework for Designing (Pseudoknotted) RNA Algorithms[END_REF]. The set of sequences is filtered to only retain valid sequences. The value of the parameter x is then adapted to match the average GC-content (induced by the value of x) with the targeted one. It can be shown that the expected GCcontent is a continuous and strictly increasing monotonic function of x, whose limits are 0 when x = 0 and n when x → +∞. Consequently, for any targeted GC-content gc ∈ [0%, 100%], there exists a unique value xgc such that generated sequences feature, on the average, the right GC-content. In practice, a simple binary search [START_REF] Waldispühl | An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure[END_REF] is used in our implementation, and typically converges after very few iterations. An optimal value for x can also be derived analytically using interpolation after Θ(n) evaluations of Z [a,b] i,j for different candidate values of x, as previously noted [START_REF] Waldispühl | An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure[END_REF] and could be implemented using the Fast-Fourier Transform [START_REF] Senter | Using the fast fourier transform to accelerate the computational search for RNA conformational switches[END_REF].

Overall complexity

S

Case 3: First position in paired to some position, but not involved in a stacking pair.

a b b S S p a ,b = x #gc(a .b ) • e -E β ∅→a b RT • Z [a ,b ] T,S • Z [b ,b] F,S /Z [a,b] N,S
Fig. 1: Stochastic backtrack procedure for a given substructure S: Either the first position is left unpaired (top), a base-pair is formed between the two extremities, stacking onto an exterior base-pair (middle), or paired without creating a stacking, defining two regions on which subsequent recursive calls are needed (bottom). For the empty structure (omitted here), the empty sequence is returned. Positions indicated in red are assigned at the current stage of the backtrack.

constants µx and σx such that the distribution of GC-content asymptotically converges towards a normal law having expectation in µx

• n • (1 + o(1)) and standard deviation in σx • √ n • (1 + o(1)
). Furthermore, the distribution of GC-content is highly concentrated, as asserted by its limited standard deviation, therefore the expected number of attempts required to generate a valid sequence when κ = 0 (resp. κ ∈ Ω(1/ √ n)) grows like Θ( √ n) (resp. Θ(1), i.e. a constant), leading to the announced complexities. Formally, since a suitable weight x must be recomputed for each targeted structure and GC-content, then the number M of iterations required for the converge can be accounted for explicitly, leading to time complexities in Θ((M +

√ n) • k • n) (if κ = 0, i.e. without any tolerance) and Θ(M • k • n) (if κ > 0).

Postprocessing unpaired regions: A local/global (glocal) hybrid approach

Due to our simplified energy model, unpaired regions are not subject to design constraints other than the GC-content, leading to modest probabilities for refolded design candidates to match the targeted structure. To improve these performances and test the complementarity of our global sampling approach with previous contributions based on local search, we used the RNAinverse software to redesign unpaired regions. We specified a constraint mask to prevent stacking base-pairs from being modified and, whenever necessary, reestablished their content a posteriori, as RNAinverse has been witnessed to take some liberties with constraints masks. As shown in Table 1 (Supplementary material), this postprocessing does not drastically alter the GC-content, so the glocal approach reasonably addresses the constrained GC-content design problem.

RESULTS

Implementation

Our software, IncaRNAtion, was implemented in Python 2.7.

We used RNAinverse from the Vienna Package 2.0 [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF]. All time benchmarks were run on a single AMD Opteron(tm) 6278 Processor at 2.4 GHz with cache of 512 KB. The penalty β, associated with invalid base-pairs, was set to 15.

Dataset

To evaluate the quality of our method, we used secondary structures from the RNA STRAND database [START_REF] Andronescu | RNA STRAND: the RNA secondary structure and statistical analysis database[END_REF]. Those are known secondary structure from a variety of organisms. We considered a subset of 50 structures selected by [START_REF] Levin | A global sampling approach to designing and reengineering rna secondary structures[END_REF], whose length ranges between 20 and 100 nucleotides. To ease the visualization of results, we clustered together structures having similar length, stacks density and proportion of free nucleotides in loops, leading to distributions of structures shown in Figure 4.

Design

We ran our method as follows. First, we sampled approximately 100 sequences per structure. Then, we use these sequences as seed in RNAinverse. Finally, we computed the MFE with the RNAfold program from the Vienna Package 2.0 [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF].

Before starting our benchmark, we asses the need for our methods and performed an analysis of the GC-content drift achieved with state-of-the-art software. Using our dataset of 50 structures, we generated 100 samples per structure with classical softwares who do not control the GC-content. Namely, RNAinverse, INFO-RNA, NUPACK:Design and Frnakenstein. We show the distribution of the GC-content of the sequences produced with these softwares in Fig. 5 those distributions.

As anticipated, we observe a clear bias toward high GC-contents and a complete absence of sequence with less than 30% of GC. This striking results motivates a need for methods that enable to explicitly control the GC-content and more precisely that enable to design sequences with low GC-content (i.e. 30% or less). In order to provide a complete overview of the performance of IncaRNAtion, we provide additional statistics for these software in the supplementary material.

Success rate

We started by estimating the success rate of our methodology and computed the percentage of sequences with a MFE structure identical to the target secondary structure. Figure 6 shows our results. We clearly see that before the post-processing step (i.e. RNAinverse) the sequences sampled by IncaRNAtion have a low success rate (first row). As mentioned earlier, this could be explained by the fact that no selection criterion has been at this stage applied to unpaired nucleotides. Remarkably, after the local search optimization (with RNAinverse) of nucleotides in unpaired regions (second row), we observe a dramatic improvement of our success rate. As expected, we observed that length is, in general, not a good predictor for the hardness of designing a structure. Instead, a high number of free nucleotides in the structure seems to be a good measure of the hardness of its design. Similarly, these data also show that designing sequences with low GC-content is challenging for all types of targets.

We investigated further the quality of the sequences generated by IncaRNAtion. In particular, we estimated the capacity of our methods to generate "good" sequences with desired folding capabilities regardless of the property to fold exactly into the target structure. In Figure 7, we show the ratio of well predicted base pairs in the MFE structure of our sampled sequences. As above, we can observe that, in all cases, the sequences that are the hardest to design are those with an extremely low GC-content. Indeed, the energetic contribution of the base pairs to the stability of the structure is weaker. Interestingly, we also notice that the most accurate sequences yield a GC-content of 70 ± 10%. Overall, we observe that all our samples have good folding properties, and that there is a correlation between the "precision" of the samples and the hardness of the design.

We noticed a highly decreased structural sensitivity for the sequences with 15% free nucleotides in the loops. However, one must remain careful interpreting this observation, as the structures within this class all originate from the PDB, and are relatively small (for the complete STRAND DB, the average length is ∼ 526nts, compared to ∼ 38nts around 15% unpaired bases).

Properties of designed sequences

In this section, we further analyze the generated sequences with a MFE structure that folds into the target structure.

A desirable feature in sequence design, is to produce samples with a high sequence diversity and stable secondary structure. Therefore, in the following we will use two useful measures which are the sequence identity of the samples, and the Boltzmann probability of the target structure in the low energy ensemble.

The sequence identity is defined over a set S of aligned sequences (in our case, all sequences have the same length and can be trivially aligned) as :

s 1 ,s 2 ∈S×S     1 |s 1 | i s 1 i ≡s 2 i 1     Seq. identity (7)
where si is the nucleotide at position i in sequence s. Intuitively, this measure captures the diversity of sequences generated by a given method. Next, the Boltzmann frequency is defined, for a structure S where Z s is the partition function of sequence s. This measure tells us how dominant is a structure S in the Boltzmann ensemble of structures over a sequence s. A high value implies a stable structure. We compute this frequency with RNAfold from the Vienna Package 2.0 [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF].

Figure 8 shows the number of solutions generated (i.e. sequences with a MFE structure identical to the target structure). Here, we note that low GC-contents have a strong (negative) influence on the number of sequences generated, and in parallel also affect negatively the sequence diversity. This observation emphasizes the difficulty to design sequences with low GC-content. Once again, large percentages of free nucleotides increase the difficulty of the task.

The thermodynamical stability of the target structure on the designed sequence is another important property when estimating the performance of RNA design algorithms. We estimate the quality of our solutions in Figure 9. First, we observe a slow decline of the structure stability (i.e. the frequency) when the target structure increases in size. Yet, for an average GC-content, the frequency stays over 10% even at size of 100 nucleotides. Next, we note that for the most difficult target structures (i.e. the longer ones or those with high percentages of unpaired nucleotides in loops) the GCcontent have a limited (almost null) influence on the stability of the target structure on the designed sequence. By contrast, this is less true for easiest and small structures with only few free nucleotides in internal loops.

Global sampling vs Local search vs Glocal approach

To conclude this study, we estimate the impact of the design methodology on the performances. More precisely, we aim to determine the merits of a global sampling approach (IncaRNAtion), compared to a glocal procedure (IncaRNAtion + RNAinverse) and a local search methodology (RNA-SSD). To our knowledge, RNA-SSD, beside IncaRNAtion, is the only software that implements an explicit control of the GC-content.

Here, we compare the running time and the sequence diversity of the solutions produced by each software. In addition, we focus on the design of sequences with low GC-contents (30% and less) as they are almost impossible to design with classical software (See Figure 5).

Figure 3 shows the running time of each software. These data demonstrate the efficiency and scalability of our techniques. In particular, this figure suggests that our strategy has the potential to be applied efficiently for designing sequences on long (and difficult) target secondary structures at low GC-content-A task that could have not been achieved before due time requirements.

Next, we show in Figure 10 the average sequence identity achieved by the various methods. Our results show that at extremely low GC-contents (i.e. 10%), IncaRNAtion slightly outperforms RNA-SSD while this advantage becomes less evident when the GCcontent increases. Our experiments on higher GC-contents (i.e. 50% and above) showed that our glocal strategy and the local search approach perform similarly. Similarly, we did not find any clear evidence that a global, local or glocal approach outperforms others when we compare at the thermodynamical stability of the target structure (data not shown).

CONCLUSION

In this article, we described a novel algorithm, IncaRNAtion, for the RNA secondary structure design problem, i.e. the design of an RNA sequence adopting a predefined secondary structure as its minimal free-energy fold. Implementing a global sampling approach, it optimizes affinity towards the target secondary structure, while granting the user full control over the GC-content of the resulting sequences. This extended control does not necessarily induce additional computational demands, and we showed the linear complexity of both the preprocessing stage and the generation of candidate sequences for the design, allowing for the design of larger and more complex secondary structures in a matter of minutes on a single processor (e.g. ∼28 mins for 100 candidate sequences for a ∼1500nts 16s rRNA). We evaluated the method on a benchmark composed of target secondary structures extracted from the RNA STRAND database. We observed good overall success rate, with the notable exception of very low targeted GC-content (10%), and a good to excellent entropy within designed candidates. Finally, we implemented an hybrid approach, using the RNAinverse software as a post-processing step for unpaired regions. This approach greatly increased the success rate of the method, allowing for the design of highly diverse candidates for almost all of the structures in our benchmark, while largely preserving the targeted GC-content.

In the future, we would like to complement this study by further investigating the potential of hybrid local/global -or glocalapproaches. A global sampling approach would capture the positive aspects of design, optimizing affinity towards a given structure while allowing the specification of expressive systems of constraints. Designed sequences would serve as a seed for a restricted local approach which, by breaking unwanted symmetries, would perform the negative part of the design, while ideally maintaining obedience to the constraints. Another perspective of this work is the incorporation of the full Turner energy model, which should in principle yield better designs for unpaired regions.

SUPPLEMENTARY DATA

Benchmark other softwares

To evaluate the performances of IncaRNAtion, we benchmark a set of classical softwares lacking GC-content control. Those are RNAinverse, INFO-RNA, NUPACK:Design and Frnakenstein. We present in Fig. 11 the average sequence identity and frequency for sequences generated them.

Benchmark IncaRNAtion +RNAinverse

To emphasize the usefulness of processing IncaRNAtion sequences with RNAinverse, we present the number of structures for which at least one sequence was generated with the desired MFE in Figure . 12 5.3 Limited impact on GC of local-search postprocessing of IncaRNAtion output

Since local search approaches tend to experience a bias towards GC-rich regions, it could be expected that our glocal approach, by postprocessing unpaired regions using a local search algorithm, would suffer from such a drift. However, as summarized in Table 1, we observed that the local search heuristic used to design nucleotides in loop regions has a very limited impact on the GC-content. For each class of GC-content, we reported the observed GC-content in the sequence initially generated by IncaRNAtion, and the observed GC-content after the RNAinverse postprocessing (as defined in Section 2.3). Our results show that the GC-content is relatively well conserved (less than 6% variation), with a general tendency of the postprocessing step to bring the GC-content back to 50%.

  Fig. 2: General workflow of our adaptive sampling algorithm (Waldispühl and Ponty, 2011).

  It was previously established[START_REF] Waldispühl | An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure[END_REF] that, for each value of x, there exists x #gc(a ) • Z[a ,b] F,S /Z[a,b] Extremities are paired, surrounded by another base-pair, forming a stacking base-pair. b = x #gc(a .b ) • e

Fig. 3 :

 3 Fig. 3: Average time in seconds to generate one sequence for IncaRNAtion and RNAinverse.

Figure 3

 3 Figure 3 presents the average times spent running IncaRNAtion +RNAinverse to generate one sequence with the required GCcontent. As expected, the time grows linearly in function of the length of the structures for IncaRNAtion.

Fig. 5 :

 5 Fig. 5: Overall GC-content distribution for sequences designed using RNAinverse, INFO-RNA, NUPACK:Design and Frnakenstein folding in the desired structure.

Fig. 4 :

 4 Fig. 4: Number of secondary structures per bin, according to our three clustering criteria.

Fig. 7 :

 7 Fig. 7: Structural sensitivity (i.e. # well predicted base pairs / # base pairs in target) of the sampled sequences MFE.

Fig. 8 :

 8 Fig. 8: Number of solutions generated with IncaRNAtion +RNAinverse on the first row and their average sequence identity on the second.

Fig. 10 :

 10 Fig. 10: Sequence identity of IncaRNAtion and RNAinverse for 10 and 30% of GC.

Fig. 11 :

 11 Fig. 11: The average sequence identity and frequency for softwares without GC-content control.

Fig. 12 :

 12 Fig.12:The first row shows the number of structures for which one generated sequence has the structure as MFE when only using IncaRNAtion. The second row shows when we process IncaRNAtion results with RNAinverse.

Table 1 .

 1 Observed GC-content of solutions returned by IncaRNAtion (2nd column) and after the application of the local search postprocessing (3rd column).

		GC-content (%) of designed sequences
	Target GC-content (%)	IncaRNAtion IncaRNAtion + RNAinverse
		(Global)	(Glocal)
	10%	15%	21%	6%
	30%	30%	33%	3%
	50%	48%	49%	1%
	70%	71%	69%	2%
	90%	83%	78%	5%