Oleksandr Zinenko

Lorenzo Chelini

Tobias Grosser

Declarative Transformations in the Polyhedral Model

Keywords: polyhedral model, loop nest optimization, tree rewriting

Despite the availability of sophisticated automatic optimizers, performance-critical code sections are in practice still tuned by human experts. Pragma-based languages such as OpenMP or OpenACC are the standard interface to apply such transformations to large code bases and loop-transformation pragmas would be a straightforward extension to provide fine-grained control over a compilers loop optimizer. However, the manual optimization of programs via explicit sequences of directives is unlikely to fully solve this problem as expressing complex optimization sequences explicitly results in difficult to read and non-performance-portable code. We address this problem by presenting a novel framework of composable program transformations based on the internal tree-like program representation of a polyhedral compiler. Based on a set of tree matchers and transformers, we describe an embedded transformation language which provides the foundation for the development of program optimization tactics. Using this language, we express core building blocks such as loop tiling, fusion, or data-layout-transformations, and compose them to higher-level transformations expressing algorithm-specific optimization strategies for stencils, dense linear-algebra, etc. We expect our approach to simplify the development of polyhedral optimizers and integration of polyhedral and syntactic approaches.

Transformations Declaratives dans le Modèle Polyédrique

Résumé : Malgré l'existence d'outils sophistiqués d'optimisation automatique, les parties des programmes dont la performance est cruciale sont toujours optimisées manuellement par des humains experts. Les langages basés sur des directives "pragma", tels que OpenMP ou OpenACC, sont une interface typique pour exprimer les transformations sur des grandes bases de code source. Telles directives pour transformer des nids de boucles seraient une extension naturelle permettant de contrôler l'optimiseur de boucles d'une manière précise. Pourtant l'optimisation manuelle des programmes à travers les séquences des directives de transformation n'est pas toujours souhaitable car ces séquences longues et complexes produisent des programmes peu lisibles et ne beneficiant pas de la portabilité de performance entre les différentes architectures matérielles.

Nous proposons une nouvelle approche pour définir les transformations composables des programmes basée sur la représentation interne d'un compilateur polyédrique sous forme de l'arbre. Grâce à un ensemble des "motifs" et "transformateurs" des arbres, nous décrivons un langage de transformation sur lequel nous basons le développement des tactiques d'optimisation. Avec ce langage, il est possible d'exprimer les transformations basiques, telles que le tuilage, la fusion ou la transposition de données, ainsi que la composition de ces transformations afin de définir une stratégie d'optimisation pour les grandes classes des programmes, telles que les pochoirs, les contractions de tenseurs, etc. Notre approche pourrait simplifier le développement des optimiseurs polyédriques et l'integration des transformations polyédriques et syntaxiques.

Mots-clés : modèle polyédrique, transformations des nids de boucles, réécriture des arbres

Introduction

Software-based solutions address an increasingly large number of significantly complex real-world problems. With combined effect of increased expressivity of high-level programming languages, increased domain-specificity, and diminishing returns of hardware scaling, the quality of optimizing compilers becomes important more than ever to keep the cost of software abstractions manageable. Compilers now have to address multiple, often hardly compatible, goals from code size to energy efficiency, from reducing program runtime to reducing the compilation turnaround time so as to improve the programmer's user experience.

Classical optimizing compiler typically makes use of several internal representations of the program it considers at different stages of the compilation process. Optimization passes are usually based on some form of a graph: abstract syntax trees (often with types and other attributes), control-flow or dependency graphs, etc. These passes often consist in performing tree transformations or graph reductions using a set of pre-defined or user-supplied rules.

Polyhedral compilation, on the other hand, transforms (parts of) the program into a different representation based on integer sets, better amenable to mathematical optimization [START_REF] Feautrier | Polyhedron Model[END_REF]. While this representation enables one to perform a complex combination of loop optimizations at once directly on the mathematical representation, it is rarely used in complete isolation from the graphbased syntactic representations since only specific constrained parts of the input program can be modeled. In practice, polyhedral optimization techniques are not only used together with conventional program transformation approaches, but also rely on classical graph and tree-based algorithms internally [START_REF] Zinenko | Modeling the conflicting demands of parallelism and temporal/spatial locality in affine scheduling[END_REF]. For example, Pluto loop scheduling algorithm operates on a dependence graph whose edges are annotated with integer sets denoting iteration-wise dependences obtained by the polyhedral analysis [START_REF] Bondhugula | A Practical Automatic Polyhedral Parallelizer and Locality Optimizer[END_REF].

The schedule tree structure [START_REF] Verdoolaege | Schedule Trees[END_REF] used by the popular isl library [START_REF] Verdoolaege | Isl: An integer set library for the polyhedral model[END_REF] to represent loop schedules is the glaring example of the practical intertwining of polyhedral and syntax-based techniques. It combines integer set-based or "polyhedral" iteration sets and schedules with the explicit loop nesting structure and coarser-grain relations between dynamic executions of program statements. It also includes ways to precisely control code generation, which are usually not considered by the polyhedral transformations themselves [START_REF] Grosser | Polyhedral AST generation is more than scanning polyhedra[END_REF]. This codeorientation of schedule tree representation has been leveraged to implement advanced transformations and extensions to the polyhedral model that would have require a separate, syntactic description otherwise.

Practical applications of polyhedral compilation techniques often perform only a part of the transformation in the polyhedral representation itself. Loop tiling is often performed as a separate step, controlled by an imperative algorithm that uses the results of polyhedral analyses. Loop fusion is often done outside the main schedule model. Device mapping, parallelization, vectorization or loop unrolling require transformations that cannot be expressed purely in terms of polyhedral schedules and are delegated to subsequent algorithmic steps.

We propose to combine polyhedral analyses and program transformations with more conventional compiler construction techniques based on tree rewriting. We introduce the pattern-matching framework for the schedule tree representation and the integer relations that are at the core of the polyhedral representation. We leverage this framework to build a declarative polyhedral program transformation engine and demonstrate how it can be used to easily express existing and new polyhedral-based optimizations through case studies. Unlike existing script-based polyhedral transformers, our search-and-match approach is not specific to a particular program but allows one to build reusable program transformation tactics in a declarative way, similarly to production optimizing compilers.

Background: Polyhedral Compilation

The polyhedral model is a unified framework to model and transform loops in imperative programs [START_REF] Feautrier | Polyhedron Model[END_REF]. After over three decades of active research and attempts to integrate into production compilers [START_REF] Bondhugula | A model for fusion and code motion in an automatic parallelizing compiler[END_REF][START_REF] Pop | Graphite: Polyhedral analyses and optimizations for gcc[END_REF][START_REF] Grosser | Polly-polyhedral optimization in llvm[END_REF], it has recently gathered new attention as a method of choice for generating efficient code targeting both CPUs and accelerators from domain-specific languages [START_REF] Vasilache | Tensor comprehensions: Framework-agnostic high-performance machine learning abstractions[END_REF][START_REF] Mullapudi | Polymage: Automatic optimization for image processing pipelines[END_REF].

The key power of the polyhedral framework is also a major obstacle to its wider adoption: instead of operating on conventional syntax-or instructionbased intermediate representations, it transforms parts of the program into a linear-algebraic form amenable to mathematical optimization. Specific algorithms and tools have been developed to convert high-level languages [START_REF] Bastoul | Putting polyhedral loop transformations to work[END_REF][START_REF] Verdoolaege | Polyhedral extraction tool[END_REF] and intermediate complier representations [START_REF] Grosser | Polly-polyhedral optimization in llvm[END_REF] as well as algorithms to generate the transformed code back [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF][START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF][START_REF] Chen | Polyhedra scanning revisited[END_REF] have been developed to make the model applicable in practice.

The internal operation of a polyhedral optimizer is widely regarded as a black-box although tools exist to describe the performed transformations in terms of loop transformation primitives [START_REF] Bagnères | Opening Polyhedral Compiler's Black Box[END_REF].

Iteration Domains and Access Relations

The algebraic representation is based on integer sets and relations. Generally, the model is applicable to static control parts (SCoPs) that consist of loops whose boundaries are affine functions of the surrounding loop iterators and parameters (values of which are unknown at compile time but guaranteed to be constant). Every iteration of these loops is identified by an integer vector in a k-dimensional space, where k is the number of nested loops, with coordinates corresponding to the values of loop induction variables. Each statement other than control flow construct has an associated symbolic name and a set integer vectors describing at which iterations it is executed, commonly referred to as iteration domain. When the static control conditions are respected, the iteration domain can be concisely denoted by a system of affine inequalities. For exam-Inria ple, the iteration domain of the gemm kernel in Listing 1 can be expressed as (N ) → {S 1 (i, j) | 0 ≤ i, j < N ; S 2 (i, j, k) | 0 ≤ i, j, k < N } in the tagged-tuple notation introduced by iscc [START_REF] Verdoolaege | Counting affine calculator and applications[END_REF]. Note that N is a parameter in this expression. Individual executions of statements inside loops are called statement instances.

for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) { S1:

D[i][j] = beta * C[i][j]; for (int k = 0; k < N; ++k) S2: D[i][j] += alpha * A[i][k] * B[k][j]; } Listing 1: Generalized Matrix Multiplication (gemm) kernel.
Internal computation is usually expressed as arithmetic expressions and function calls with array elements as arguments. The conditions on the array subscripts are essentially the same as those on control flow conditions. Array accesses can be thus precisely encoded as relations between vectors in the iteration spaces and vectors in the array space whose coordinates are values of the accessed subscripts. These relations are most often defined by piece-wise quasiaffine functions. In our running example from Listing 1, the access relation for S1

is (N ) → {S 1 (i, j) → beta() | 0 ≤ i, j < N ; S 1 (i, j) → D(d 1 , d 2 ) | d 1 = i ∧ d 2 = j; S 1 (i, j) → C(c 1 , c 2 ) | c 1 = i ∧ c 2 = j}.
Polyhedral access relations allow for fine-grain exact data-flow analysis [START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF] and enable loop transformations [START_REF] Bastoul | Mapping Deviation: A Technique to Adapt or to Guard Loop Transformation Intuitions for Legality[END_REF][START_REF] Girbal | Semi-automatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF].

Schedules and Schedule Trees

The iteration domain defines which statement instances should be executed but not in which order. The latter is defined by a schedule that maps points in the iteration space to points in the time space. Statement instances are executed in the lexicographical order of their coordinates in the time space.

While it is possible to express schedules as relations or piece-wise quasi-affine functions, it is often undesirable. This is because statements are likely to share part of their schedule due to the inherent nested strcuture of the generated code. It is also necessary to reflect the relative syntactic order of loops and statements within them, which is often achieved through auxiliary dimensions whose values are always constant for the given statement. 1To address these challenges, Verdoolaege et.al. proposed the schedule tree structure as a way to represent schedules in the polyhedral model [START_REF] Verdoolaege | Schedule Trees[END_REF]. In this structure, statements that share a common partial schedule share an ancestor that describes this partial schedule, removing duplication. The loop nesting maps naturally to parent/child relationships whereas syntactic ordering maps to the left-to-right order of sibling nodes. Beyond simple schedules, schedule trees have been successfully used to support multi-target code generation with offloading mechanisms [START_REF] Grosser | Polyhedral AST generation is more than scanning polyhedra[END_REF][START_REF] Verdoolaege | Polyhedral parallel code generation for cuda[END_REF] and advanced hybrid polyhedral/syntactic transformation techniques [START_REF] Grosser | Hybrid hexagonal/classical tiling for gpus[END_REF][START_REF] Zhao | A polyhedral compilation framework for loops with dynamic data-dependent bounds[END_REF].

Nodes in schedule trees have one of the numerous types. Let us briefly present node types relevant to our examples.

domain-the iteration domain, always located at the root of the tree; band -partial schedule of one or multiple loops (the name refers to the notion of tilable band of loops [START_REF] Bondhugula | A practical and fully automatic polyhedral program optimization system[END_REF]);

filter -restricts the subtree to a subset of the iteration domain;

sequence-imposes the order between its children (statements or loops); extension-introduces new statement instances local to the subtree with respect to its prefix schedule.

All nodes except sequence have no more than one child. For the vectors in the iteration domain, one can reconstruct the relational form of the schedule by taking a range product of all partial schedules in the tree.

Listing 2 presents a simple 2D stencil program. The corresponding schedule tree is depicted in Figure 1. It uses two band nodes: the outer one provides a partial schedule representing the individual time steps. The inner one provides a partial schedule enumerating all data-points within a given time step. Finally, the individual statement types are enumerated within a sequence node. When combined across all dimensions, the schedule tree defines an execution order identical to the original execution order.

The tree structure of a schedule tree allows for the easy application of local transformations. Such transformations include the combination of nested band nodes into a single band, the splitting of a multi-dimensional band node into individual bands, but also more complex transformations such as tiling of nodes or even the application of affine transformations.

for (int t = 0; t <= T; t++) for (int i = 1; i < 1023; i++) for (int j = 1; j < 1023; j++) if (t % 2 == 0) S:

A[i][j] += B[i-1][ j] + B[ i][j+1] + B[ 0][ 0] + B[ i][j-1] + B[i+1][ j]; else T: B[i][j] += A[i-1][ j] + A[ i][j+1] + A[ 0][ 0] + A[ i][j-1] + A[i+1][ j];
Listing 2: A simple 2D stencil kernel This schedule can be rewritten in the relational form as (T ) → {S(t, i, j) → (t, i, j, 0); T (t, i, j) → (t, i, j, 1)}. Note that duplication of (t, i, j) and the presence of auxiliary trailing dimensions.

domain: S(t, i, j) | 1 ≤ i, j ≤ 1022 ∧ 0 ≤ t ≤ T ∧ t mod 2 = 0; T (t, i, j) | 1 ≤ i, j ≤ 1022 ∧ 0 ≤ t ≤ T ∧ t mod 2 = 1
band: S(t, i, j) → (t); T (t, i, j) → (t) band: S(t, i, j) → (i, j); T (t, i, j) → (i, j) sequence filter: S(t, i, j) filter: T (t, i, j)

Declarative Schedule Tree Transformations

Unlike conventional compilation approaches centered around syntax trees, polyhedral compilation techniques typically operate on some representation of a schedule disconnected from any syntactic representation. While schedule representations are often based on algebraic objects such as sets or relations in vector spaces offering large sets of operations and consistency guarantees, they do not easily lend themselves to classical tree rewriting compiler passes. Thanks to the nature of the schedule tree representation that combines schedule and syntactic aspects in a tree shape, it becomes possible to apply tree matching and rewriting techniques without leaving the polyhedral representation.

Structural Schedule Tree Matchers

As first step, we define the schedule tree matchers to enable declarative description of pattern in the schedule tree. A structural matcher essentially replicates the node type-based structure of the schedule tree with additional filtering and wildcarding capabilities. A node matcher consists of an expected (possibly unspecified) node type, a list of child node matchers and an optional filter. In addition to all schedule tree node types, a matcher may expect any node type or a non-empty list thereof. Since schedule trees are often deeper than wider due to only two nodes types being allowed to have multiple children, matching a single node of any type is often sufficient to recognize large portions of a schedule tree. Tree leaves can be matched using a special node type.

Matching Procedure

The tree matcher is specified by the top node it must match, referred to as relative root. The matching procedure starts at the specified node in the schedule tree and performs simultaneous depth-first preorder traversal of the schedule tree and the tree formed by descendants of the relative root matcher. If a mismatch is detected until the entire matcher is traversed, it is immediately reported and the traversal stops. Our choice of explicitly marking tree leaves in the matcher allows us to easily transform the matcher into a sequential form if necessary. Once this transformation is performed, it becomes possible to apply more efficient string-based algorithms for matching against a set of strings or finding all substrings by treating each combination of node type and filter as a unique symbol in some abstract alphabet.

Programming Interface

The API to construct schedule tree matchers is designed to visually resemble the structure of the tree itself in a declarative way. Named variadic functions correspond to node types, and the argument lists enumerate child nodes. This approach enables static checking of most tree invariants, e.g., some node types only allowed to have one child or children of a specific type. Leading arguments include optionally a filtering function and a reference to the "placeholder" node.

The filtering function allows the caller to control the matching more precisely by considering non-structural aspects of the schedule tree. For example, matching only one-dimensional band nodes requires both structural and non-structural properties. The references are used to capture certain nodes in the matched subtree, similarly to captured groups in regular expressions, for future use by the caller. If the captured nodes are non-empty, the underlying subtree is guaranteed to have the structure described by the matcher. The matcher shown in Listing 3 matches all the sub-trees starting at a sequence node with exactly two filters as children, the first of which has a permutable band as a child (checked via a user-provided function isPermutable) while the second has no children.

Extensibility of the Approach: Traversal Matchers

Practically, filtering functions are not limited to matching decisions and can be used to extend the tree matching engine. For example, they can be used to include arbitrary tree traversals in the matching procedure. It suffices to provide a higher-order function that produces the necessary filtering functions, Inria typically closures over the traversal procedure and the nested matcher. Listing 4 illustrates the approach by defining a "descendant" matcher for coincidencefeaturing leaf bands and by using it while looking for sequence nodes that can have such descendants.

auto isCoincidentLeaf = band(hasCoincidence, leaf()); auto matcher = sequence(hasDescendant(isCoincidentLeaf), anyForest());

Listing 4: The matcher API leverages functional programming techniques to support extensions. Here, hasDescendant is a higher-order function.

Schedule Tree Builders

The imperative-style schedule tree construction interface does not allow for declarative tree construction. Instead, we provide schedule tree builders whose programming interface is close to that of the tree matchers. Named variadic functions specify the type of the node, the nesting of the function calls reflects the structure of the tree and the optional leading arguments accept functions that are used to build the non-structural properties of each specific node type (partial schedules for band nodes, conditions for filter nodes and so on). Once again, declarative-style interface requires other components of functional programming to be expressive enough, in particular, to enable lazy or delayed initialization of non-structural properties. The builder itself is merely a description of the subtree to build (one may think of it as an alternative, simplified schedule tree structure). It may be transformed into a standalone tree if possible, i.e., if it is rooted at a domain or an extension node, or grafted at a leaf of the existing tree.

Putting it All Together

The aforementioned lazy initialization supports the declarative find-and-replace procedure on schedule trees. First, a pattern is described using schedule tree matchers, optionally with callback functions to check additional node properties. Then, the transformed tree structure is described using the builders with delayed initialization of the node properties. Finally, one can traverse the tree in any order, trying to match a subtree rooted at the node being traversed and, in case of a match, replace the subtree with the one produced by the builder. Listing 5 demonstrates how simple rectangular loop tiling transformation can be declared. The matcher looks for permutable band nodes anywhere in the tree and captures both the node and its child subtree. The builder splits the band into two nested bands by taking the integer division and the modulo parts of the schedule while keeping the child subtree intact. It can be trivially extended to support more advanced techniques like full/partial tile separation or to limit tiling to the deepest bands without loosing a significant portion of the transformation code clarity. In isl notation, individual relations in a union are uniquely identified by the space in which they live.

A polyhedral relation matcher is formed by a non-empty set of groups of relation property descriptors. For a relation to be considered for a match, it must have all properties described in a group. The connection between groups in the same matcher is controlled by the user. In particular, they may be allowed or not allowed to match the same relation in the union.

Beyond checking the properties, we are interested in capturing some information about the matching relation, e.g. the coefficients of the linear access function or whether a relation that was matched thanks to its surjectivity is also bijective. The capturing mechanism operates through placeholders, each of which has two data components: a constant pattern and a variable candidate. Each relation is checked against the pattern and, in case of a match, it may yield one or more candidates. The description of what constitutes a match and how the candidates are generated is external to the matching engine and can, e.g., be provided by the user. An example of a pattern that yields multiple candidates is "access relation with some output dimensions fixed to a non-parametric constant". In this case, a candidate may be generated for each output dimension fixed to a constant or, for other use cases, for each value of the constant.

Inria

A match is an assignment of candidates to placeholders. One union of relations may yield zero or more matches against the given matcher. In many cases, candidates assigned to different placeholders are required to be distinct. We require this by default, considering other cases invalid assignments. At the same time, if a placeholder is reused within the same matcher expression, we also require that all appearances get assigned the same candidate. In any case, the candidate comparison does not take into account the differences between spaces of the relations, only the candidate descriptions. This behavior allows the matcher to connect different relations in the union with more precision. For example, it captures the fact that exists two relations that are either both bijective or both non-bijective. To support edge cases, the caller can override the definition of valid assignment, e.g., to allow the same candidate to be assigned to different placeholders.

Matcher Specification, Extension and Composition

The implementation of the matching procedure, or the matching engine, is a template definition in C++ sense. An instantiation of the matching engine is specified by data structures for the pattern and the candidate. In addition, it implements functions to define a set of candidates for a given access relation and whether a candidate assignment to the matchers constitutes a valid match. Attempts to use different types of patterns in the same matcher expression are disallowed by construction. If several properties need to be combined to find candidates, the caller must defined how exactly two properties should be combined together, e.g. both present at the same time, at least one, and so on. We argue that, at a price of small code overhead, it allows to support an arbitrarily large combinatorial space of patterns combined freely from the simple ones.

Matching Procedure

The matching procedure is decoupled into two stages, providing sufficient flexibility without sacrificing performance.

First, the engine traverses all relations one by one and defines, for each placeholder, the set of suitable candidates using the user-provided function. If at least one of the placeholders has no suitable candidates, the absence of match is reported immediately and the procedure stops. Then the engine traverses the space of all possible assignments of candidates to the placeholders and checks whether the assignment is valid using another user-provided function. Since the space of possible assignments is combinatorially large, we opt for a branchand-cut traversal approach. Partially-formed assignments are passed to the validation function before adding one more placeholder-candidate pair to the assignment. If the partial assignment is reported invalid, further exploration is not performed. This approach can be easily transformed into a branch-andbound if the assignment needs to be optimal in some sense, or altered to change the exploration/validation ratio if the validation itself is expensive. In cases where invalid partial assignment does not mean that full assignment cannot be valid, the validity filter can always return "valid" for partial assignments.

Programming Interface by Example: Access Matchers

Let us illustrate relation matchers as applied to polyhedral access relations. In most cases, these relations are in fact affine functions in vector spaces mapping the schedule vector to the subscript vector. We provide pattern and candidate descriptions for simple affine expressions (ω = k * ι + c where k and c coefficients form the pattern while ω and ι define a candidate by matching one of the output and input dimensions respectively) as well as for fixed non-parametric strides in the relation range. When targeting access relation unions, we can assume that the source space of all relations is the same, e.g., the schedule space, so it is convenient to only operate on the relation range.

The programming interface is similar in spirit to that of schedule tree matchers, except for the tree parent-child relations, and is based on nested function calls progressively constructing the matcher. Listing 6 illustrates how one can identify if the same 2D array is accessed directly and with transposition. The example above uses additional syntactic sugar to operate separately on output dimensions of the access relation. Contrary to the generic property matchers that are applicable to entire relations, for accesses, it is often required to separately consider individual output dimensions. In the transposed access example, the caller would like to know that, in (. . .

) → {(i, j) → A(a 1 = i, a 2 = j); (i, j) → B(b 1 = j, b 2 = i)}
, the first dimension of the output space in one relation (a 1 ) is equal to the second dimension of the output space in another relation (b 2 ). We handle such situations by augmenting the pattern description with the expected position in the output space, projecting the access relation onto that dimension and continuing to check the properties on a relation with one-dimensional output space. Such per-dimensional behaviour is particularly Inria useful for, e.g., detecting the presence of temporal or spatial locality in array accesses. However, it is not applicable beyond relations that are also affine functions in some vector space.

Relation Transformers

Since relation matchers are significantly less constrained than tree matchers (due to relations being also less structured), we opted for an approach that expresses directly the transformation of a relation rather than building it from scratch. Each pattern-candidate pair may optionally provides a function that transforms a relation given the candidate assignment. If provided, this function will be repeatedly applied to the matching relation for each candidate assignment. A sanity check verifies that the same relation does not get transformed through different placeholder groups (it may be matched several times as long as only one of the matches has an associated transformation function) and that each placeholder has a unique assignment within each relation.

For access relations, it is often convenient to rewrite some dimensions, i.e. array subscripts, keeping others without modification. Relation transformers can achieve that by using the same per-dimension separation mechanism as the relation matchers.

Listing 7 exemplifies the use of relation transformers to transpose accesses to 2D and 3D arrays. Note that each "replace" call describes an individual pattern that gets checked separately during the same pass of the replacement procedure.

isl::ctx ctx = /*...*/ ; auto _i = placeholder(ctx); auto _j = placeholder(ctx); auto _k = placeholder(ctx); auto result = findAndReplace(reads, replace(access(_j, _i), access(_i, _j)), replace(access(_k, _j, _i), access(_k, _i, _j)));

Listing 7: Access relation transformation transposing the two innermost subscripts of 2D and 3D arrays.

Case studies

In this section we demonstrate the applicability and flexibility of our framework considering three use cases. We start by showing how our framework improves upon state-of-the-art user-driven optimizers such as Clay, decoupling the transformation strategy from the code structure. Subsequently, with the second use case we demonstrate how we can achieve the same level of optimization of a newly introduces compiler optimization in the LLVM polyhedral optimizer in a more flexible and concise way, effectively reducing code requirements and unlocking the opportunity for a rule-based optimization engine. We conclude with the third use case showing how we can implement a data-layout transformation to reduce stream alignment issue in stencil patterns, namely dimension-lifted transposition.

Case: Generalizing Directive-Based Loop Transformations

Several tools were proposed to perform classical loop transformations, such as loop fusion or tiling [START_REF] Wolfe | High Performance Compilers for Parallel Computing[END_REF], in the polyhedral framework, including AlphaZ [START_REF] Yuki | Alphaz: A system for design space exploration in the polyhedral model[END_REF], CHiLL [START_REF] Chen | Chill: A framework for composing highlevel loop transformations[END_REF][START_REF] Rudy | A Programming Language Interface to Describe Transformations and Code Generation[END_REF], Clay [START_REF] Bagnères | Opening polyhedral compiler's black box[END_REF] among others. All these tools expose some language to specify program transformations applicable to loops. Although the loop transformations are abstracted to their common names, the languages can be seen as imperative in a sense that they require to specify which loops are targeted using external tags or language-level annotations. Consider how tiling a simple loop nest is expressed in Clay, as shown in Listing 8. Each transformation directive starts with a "beta-prefix" that identifies the loop it applies to, followed by two target loop depths (where to place the tiles loop and the points loop in the nest) and, finally, by the requested tile size. Adding, for example, a time loop around this transformation, or an initialization statement for C[i][j] would break the transformation script immediately since the "beta-prefixes" and/or the loop depths would change. Worse, the user must keep track of the transformation effects on the loop structure to properly spell subsequent transformations. Also note that Clay and the related Chlore algorithm internally transform beta-prefixes into a tree structure, similar to the combination of sequence and filter nodes in the schedule tree [START_REF] Zinenko | Interactive Program Restructuring[END_REF].

/* Clay

tile([0,0,0],1,1,32); tile([0,0,0,0],3,2,32); tile([0,0,0,0,0],5,3,32); */ for (int i = 0; i < 1024; i++) for (int j = 0; j < 1024; j++) for (int k = 0; k < 1024; ++k)

C[i][j] += alpha * A[i][k] * B[k][j];
Listing 8: Core statement of the matrix-matrix multiplication.

With our approach, it is possible to make the transformations more declarative. Instead of binding it to specific loops, we can declare how a compati-Inria ble schedule tree should look like and express the transformations in terms of matched nodes and relations between them. Both the matcher and the builder are expressed as compilable C++ code in Listing 9, offering virtually unlimited flexibility as to controlling program transformations. Listing 9: Declarative transformation for tiling 3D loops with a tile size of 32.

The search-and-replace algorithm can now be applied to any arbitrarily complex schedule tree, repeatedly if necessary, tiling three-dimensional loops in the entire program. Note that we are using helper functions to define properties of the nodes in the transformed trees, such as tileSchedule and pointSchedule computing partial schedules of the bands as shown in Listing 9. A number of this function is shipped with the framework and others can be easily derived using isl .

Case: Hand-Tuned Optimization for GEMM-like Kernels

Generalized matrix multiplication (gemm BLAS kernel) is one of the important computation patterns and is the most optimized kernel in history [START_REF] Low | Analytical modeling is enough for high-performance blis[END_REF]. Preoptimized versions of gemm are commonplace, yet they often under-perform on matrix sizes different from those they were designed [START_REF] Vasilache | Tensor comprehensions: Framework-agnostic high-performance machine learning abstractions[END_REF]. In such cases, even state-of-the-art compilers only achieve a fraction of the theoretical machine performance for a simple textbook style implementation [START_REF] Gareev | High-performance generalized tensor operations: A compiler-oriented approach[END_REF].

A recent improvement in Polly loop optimizer introduced a custom transformation for gemm-like kernels that is controlled outside of the main affine scheduling mechanism [START_REF] Gareev | High-performance generalized tensor operations: A compiler-oriented approach[END_REF]. This transformation applies to a generalized case of tensor contraction of the form

C[i][j] = E(A[i][k], B[k][j], C[i][j])
where the dimension k is contracted and E is some operation between tensors. The matrix-matrix multiplication from Listing 8 is an example of such contraction with E = (×, +). To qualify for the transformation, the kernel must respect the following: perfectly nested loop satisfying the requirements of the polyhedral model; contains three non-empty one dimensional for loops with induction variables incremented by one; contains an innermost statement C πC(IJ) = E(A πA(IK) , B πB(KJ) , C πC(IJ) ) where A πA(IK) , B πB(KJ) , C πC(IJ) and πA(IK), πB(KJ), πC(IJ) are access to tensors A, B, C and permutations of the enclosed indexes respectively. E is a generic expression that contains at least three reads from tensors A, B and C. the interchange of I and J is valid, while for P is interchangeable iff P contains only one element or an associative operation is used to update C.

The candidate loops are found implementing the aforementioned conditions as schedule tree and access relation matchers. In particular, we are looking for a three-dimensional permutable band with a single statement (leaf) featuring specific access patterns: at least three two-dimensional read accesses to different arrays; an permutation of indices that satisfies the placeholder pattern The kernel transformation is derived from [START_REF] Gareev | High-performance generalized tensor operations: A compiler-oriented approach[END_REF] as follows: First, we re-arrange the band dimensions such that j will be the outermost dimension, followed by k and i. Then we apply multiple level tiling and loop interchange to create macro and micro-kernel. Specifically, the builder for the macro-kernel reported Inria in Listing 11 performs L2-cache tiling and interchanges the newly created point loops. After applying the first builder, we define and apply the second builder to create the micro-kernel tiling the points loops so that they fit into vector registers and fully unroll the new innermost loops to simplify subsequent vectorization. Listing 12: Further tiling and unrolling to enable vectorization.

[i, j] → [i, k][k, j].
Overall, the combination of matcher and builder holds in dozens of lines whereas Polly's implementation of the same transformation needs several hundred lines.

We believe that our declarative approach can greatly simply the optimizer implementation through collection of pattern-based optimizations. It also unleashes the opportunity for a rule-based optimization engine where expert user can try out new optimization strategies, and plug them in the optimizer by declaratively restricting the parts of the program to which they apply.

Case: Stencil Vectorization through Data-Layout Transformation

Stencil patterns occupy a special place in the loop optimization literature. As they typically involve accesses to multiple adjacent array elements along multiple dimensions inside nested loops with (mostly) static control flow, they fit the polyhedral model. Furthermore, many important practical use cases can be implemented as stencils, for example difference schemes for systems of differential equations or convolutional neural networks. However, polyhedral compilation based on affine scheduling are often counterproductive for stencil patterns. Direct attempts to minimize dependence distances lead to effective serialization of computation and/or substantially complex control flow overhead, making the transformed code slower than the orig-inal [START_REF] Verdoolaege | Consecutivity in the isl polyhedral scheduler[END_REF][START_REF] Zinenko | Unified Polyhedral Modeling of Temporal and Spatial Locality[END_REF]. Several techniques, more or less closely related to the polyhedral compilation, address specifically stencil parallelization and vectorization [START_REF] Grosser | Hybrid hexagonal/classical tiling for gpus[END_REF]. Since these techniques are often not adapted to non-stencil cases, one has to first find the stencil-like part of the program and then apply the transformation.

Declarative schedule tree and relation patterns simplify the matching of stencils in the polyhedral representation. Let us illustrate how a stencil vectorization technique based on data layout transformation that has not been defined in terms of polyhedral model [START_REF] Henretty | Data layout transformation for stencil computations on short-vector simd architectures[END_REF] can be expressed in our framework.

Candidate Loops

The data-layout transformation (DLT) technique [START_REF] Henretty | Data layout transformation for stencil computations on short-vector simd architectures[END_REF] applies only to innermost loops with no loop-carried dependences. In schedule tree nomenclature, we should start by finding all band nodes that don't have other band nodes as children and that have the last coincidence flag set (assuming coincidence was computed based on all dependences). band(capturedBand, and_(

Note that if the band is permutable, it does not matter whether the coincident flag is set for the last dimension or any other of the band member since the loops can be trivially permuted. This extension to the original technique and the corresponding tree transformation are easy to propose and implement in our approach but were not considered by the original paper.

Access Strides

The transformation validity is further restricted to statements that access either the same or the consecutive elements of an array in consecutive iterations of the innermost loop. This can be easily matched against by computing the range stride of the access relation, i.e. the constant distance between elements in its range. From the relation matcher perspective, the pattern includes the expected stride along with the algorithm to compute it while the candidate is just the space of the matching relation. This can be represented by saying all the last (Python-style index -1) access must have stride either 0 or 1. match(rw, access(dim(-1, stride(ctx, 0)))).size() + match(rw, access(dim(-1, stride(ctx, 1)))).size() == match(rw, access(any())).size() Vector Lane Conflicts Finally, DLT is only necessary when there exist vector lane conflicts that cannot be removed through loop shifting, in particular, if the same value is accessed through different references in different iterations of the innermost loops. This condition can be transformed into a sequence of set operations forming a system of constraints followed by an emptiness check.

Inria

While it can be used inside the relation matcher to create a list of candidates before checking that all accesses do match, it is arguably more pragmatic to perform the operations directly.

Injecting Data Layout Transformation The data layout transformation consists in placing hitherto adjacent array elements at distance L from each other, where L is the number of vector lanes. This is expressed by a union of affine functions defined over disjoint subdomains s∈[0,L),s∈Z {(

ι, i) → ( α, a) | α = ι ∧ a = s -B L + L(i -s•(B U -B L +1) L ) ∧ s B U -B L +1 L ≤ i -B L < (s + 1) B U -B L +1
L } where B L and B U are the lower and the upper inclusive bound on the accessed elements. If this transformation had been performed on the iteration space rather than on the subscript space, it would have corresponded to loop strip-mining followed by loop interchange and coalescing. Once such union of affine functions is constructed, it can be used to declare a schedule tree builder that injects the transformation itself as presented in Listing 13. First, it introduces the new statements for the copies to and from the transformed array through an extension node. Below it, copies "to" are scheduled before the main computation, which itself is scheduled before the copies "from" using a combination of sequence and filter nodes. Finally, partial schedules are specified for the copies and the original subtree is replicated for the main computation. In practice, the transformation is only applicable to loop iterations that fully fit into vector lanes2 and requires additional edge case otherwise. While these cases can also be expressed declaratively as tree builders, they are not essential to our presentation so we omit them for the sake of clarity. The code below features the corresponding builder construction, which clearly indicates the intended subtree structure.

This example lays the foundation for most data-layout or memory-related transformation that rely on copying the data from one array to another and back. Most such transformations follow the extension/sequence/3-filter pattern with only the properties of nodes differing. Copies "to" or "from" additional memory may be omitted in cases where the initial or the final values are irrelevant to the task, creating two other possible tree transformation patterns.

Access rewriting

The final step of DLT consists in changing the original stencil computation to access the transformed array instead of the original one and, optionally, emitting vectorization pragmas. To enable vectorization, it is necessary to access the transformed array in the sequential order and the data required by each iteration now requires different subscripts. In particular, hitherto adjacent elements are now located at distance L. This can be easily expressed using our relation rewriting mechanism. The pattern constitutes the single access subscript (last range dimension of the access relation) and the candidates are affine functions that define that subscript Listing 13: Base case of data layout transformation for stencil vectorization.

Other data motion transformations follow the same tree structure.

can be decomposed into the linear and the constant part, and the latter needs to be scaled by L in order to reflect the data layout transformation.

auto aff = affine(ctx); findAndReplace(stride1Accesses, replace(access(dim(-1, aff)), access(dim(-1, aff.payload_.linear + L * aff.payload_.constant)));

Only the stride-1 accesses should be transformed. They can be easily extracted by reusing a part of the polyhedral relation matcher used during the detection phase. Depending on the underlying implementation, this rewriting can happen directly on the access relations or during code generation.

Even including the pattern and candidate descriptors, the entire code for implementing this advanced vectorization technique fits into a couple of hundred lines and can be easily combined with other techniques expressed in the same way.

Case: Mapping to GPUs and Accelerators

Recent evolutions in hardware for high-performance computing led to the increasing popularity of the accelerator-based architecture. GPUs and other accelerators now power world's fastest supercomputers. Most accelerators share a common programming model that consists of: one or several accelerator programs using a specialized language or a subset of existing language (typically C++), referred to as kernels; Inria a host program that manages the acclerator device, explicitly transfers the data between host and accelerator, and schedules the kernels for execution.

To make computation efficient, accelerators impose stringent conditions on the operations that can be executed by the kernels, especially when it comes to parallel execution. Therefore, the first task most automatic device mappers perform is to find parts of the program that can be executed on the accelerator.

Let us consider GPU mapping as an illustrative example using the algorithm implemented in ppcg that offers a good trade-off between performance and simplicity [START_REF] Verdoolaege | Polyhedral parallel code generation for cuda[END_REF]. ppcg operates on schedule trees and starts by traversing the tree from the top to find band nodes that have both the permutability and at least one coincident flag. Subtrees rooted at such nodes are mapped to the GPU entirely and therefore not traversed in search of other bands. This can be expressed declaratively as "permutable coincident bands that don't have permutable coincident bands as ancestors", which can be trivially transformed into a schedule tree matcher as shown in Listing 14.

isl::schedule_node node, child; auto mappable = and_(isPermutable, hasCoincidence); auto matcher = band(node, and_(mappable, not_(hasAscendant(band(mappable, anyTree())))), anyTree(child)); Listing 14: Schedule Tree Matcher mimicing ppcg's behavior.

When all such nodes are found, the basic mapping transformation is performed as follows. First, a special context node introduces thread and block identifiers to the program. They can only appear below this node, which separates the top part of the tree that is transformed to host code from the bottom subtree that is transformed to the GPU device code. Below the context a extension node introduces calls to the data motion functions that ensure the data is transferred to and from device before and after the kernel, respectively. Depending on the host API being synchronous or not, explicit device synchronization calls can be emitted as well. Data transfers and synchronizations are ordered using a sequence node with corresponding filters. Finally, the original band node is tiled with the resulting tiles loops being mapped to blocks and points loops being mapped to threads using filter nodes. 4 All these steps can be concisely expressed as a schedule tree builder in Listing 15.

Alternatively, the GPU mapping can be composed from different builders applied in sequence, hence reusing their definitions. First, a data-copy builder similar to that of Listing 13 performs data copies. Then, a tiling builder similar to Listing 9 kicks in. Finally, a thread-and block-mapping builder, or even the Once the schedule tree is transformed, a full implementation of an automatic GPU mapper must also provide a customized code generator that emits either CPU or GPU code depending on its position in the schedule tree. This can be done, for example, by keeping track of the points that belong to a subtree that got mapped to the GPU, which can be obtained when the matcher is called.

While this example demonstrates only the minimal GPU mapping strategy, it shows how easily one can devise new strategies using our framework. For example, one can chose to exploit finer-grain parallelism by mapping the innermost mappable loops instead of the outermost. This only requires to replace hasAscendant with hasDescendant in the matcher. More complex strategies like the one proposed in Tensor Comprehensions [START_REF] Vasilache | Tensor comprehensions: Framework-agnostic high-performance machine learning abstractions[END_REF] can be implemented as a series of matchers and builders. Finally, intra-kernel optimizations such as explicit copies to fast shared or private memory can be easily obtained by modifying the data-copy builder of Listing 13. In practice, data-copy transformations only differ by the mapping function they use, while the schedule tree structure remains essentially the same. Thanks to tree node properties being decoupled from its structure in schedule tree builders, it is possible to make the copy builder parameterizable with the mapping and derive the remaining properties from that. Similarly, the entire device mapping builder can be adapted to different architectures, including the upcoming ones for, e.g, near-memory computing, by simple changes to node properties and code generation. [START_REF] Kelly | A framework for unifying reordering transformations[END_REF]. Subsequently, Girbal et.al. proposed the URUK framework to apply loop transformations for cache hierarchy and parallelism based on unimodular schedules [START_REF] Girbal | Semi-automatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF]. Yuki et.al. developed Al-phaZ, a framework that allows the users to express the program as a set of equations based on the Alpha language and manipulate it using script-driven transformations [START_REF] Yuki | Alphaz: A system for design space exploration in the polyhedral model[END_REF]. In their implementation, the authors also support memory (re)-allocation and explicitly represent reductions. Similarly, Yi et.al. presented POET, an interpreted program transformation language designed for applying and explore complex code transformations in different programming languages [START_REF] Yi | Poet: A scripting language for applying parameterized source-tosource program transformations[END_REF]. Donadio et.al. introduced Xlanguage, an embedded DSL based on C/C++ pragma-based that allows the users to generate multi-version programs specifying the type of transformations to apply as well as the transformation parameters [START_REF] Donadio | A language for the compact representation of multiple program versions[END_REF]. Bagnères et.al provided feedback from a polyhedral compiler by expressing it as a sequence of loop transformation directives [START_REF] Bagnères | Opening polyhedral compiler's black box[END_REF]. Their input language, Clay, and the related Chlore algorithm allow the users to examine, refine or freeze sequences of loop transformation directives. Chen et.al. introduced CHiLL, a high-level transformation and parallelization framework that makes use of a model-driven empirical optimization engine to generate and evaluate different code variants [START_REF] Chen | Chill: A framework for composing highlevel loop transformations[END_REF]. It also allows the user to define a transformation script specifying how the code should be optimized. The framework has also been extended with high-level constructs to generate GPU code [START_REF] Rudy | CUDA-CHiLL: A programming language interface for GPGPU optimizations and code generation[END_REF].

The vast majority of polyhedral directive-based loop transformation engines require the user to target specific loop nests or statements within them, as well as to track the changes to the loop structure throughout the transformation process. Our approach on the other hand is to declaratively describe parts of the program that are suitable for a particular optimization making it applicable to arbitrarily complex programs without modifying the transformation script. Since the schedule tree also captures the syntactic information via sequence and filter nodes, one can also use our framework to target specific syntax-driven parts of the program rather than their properties.

Code Offloading for Accelerators

Hseih et.al. proposed to offload code to accelerators based on a simple cost function [START_REF] Hsieh | Transparent offloading and mapping (tom): Enabling programmer-transparent near-data processing in gpu systems[END_REF]. The key idea is to statically identify the code blocks with greater potential in bandwidth saving. If the bandwidth saving in offloading the code block is higher than the cost of initiate and complete it the offload block is marked as "beneficial" to be offloaded. Pattnaik et al. developed an affinity prediction model to understand where to execute a given kernel, i.e., main core or near-memory accelerator [START_REF] Pattnaik | Scheduling techniques for gpu architectures with processing-in-memory capabilities[END_REF]. The prediction model is based on five prediction metrics (i.e. memory intensity) and it is computed offline. Nair et.al. pro-posed a code offloading mechanism based on OpenMP 4.0 user annotations [START_REF] Nair | Active memory cube: A processing-in-memory architecture for exascale systems[END_REF]. Their compiler uses OpenMP pragmas to identify code section that will be executed on the near-memory accelerator, and identify the data region accessed. CAIRO relies on a LLC cache profiler and analytical model to decide potential offloading candidates [START_REF] Hadidi | Cairo: A compiler-assisted technique for enabling instruction-level offloading of processing-in-memory[END_REF]. Nai et.al. proposed to offload all the graph operation on the graph properties on the near-memory accelerator mapping host atomic instructions directly into near-memory atomic instructions using uncacheable memory support available in modern architectures [START_REF] Nai | Graphpim: Enabling instruction-level pim offloading in graph computing frameworks[END_REF]. Ahn et.al. proposed an ISA extension on the host processor allowing the programmer to decide what to offload [START_REF] Ahn | Pim-enabled instructions: A low-overhead, locality-aware processing-in-memory architecture[END_REF].

All these algorithms are complementary to our approach as they mostly control when to perform the offloading while featuring a hard-coded implementation of how to perform it. They can be used in the matching procedure to find subtrees suitable for offloading while the transformation is performed by a combination of schedule tree builder, relation rewriter and code generator.

Conclusion

We presented a flexible framework to perform pattern matching and declaratively express transformations on schedule trees, an internal representation of the polyhedral compiler. Unlike the majority of polyhedral approaches, pattern matching tree-rewriting transformations in our framework can be implemented and augmented using classical compiler construction technology applicable to other tree-like structures in the internal representation of a compiler. We believe our approach can simplify the integration of polyhedral optimizers into production compilers.

Our approach allowed us to easily express and compose a set of program transformations that use the polyhedral analyses but do not require the conventional combination of affine loop transformations, usually obtained by solving linear optimization problems. Complementary to existing loop transformation directives, it now becomes possible to declare, refine and reuse transformations for, e.g., data layout or device mapping. Future work will focus on serializing and sharing the match and build patterns for both the schedule trees and the expressions, as well as on more efficient matching algorithms and whole-program uses. A domain-specific declarative polyhedral language can also be considered to avoid the syntactic hurdles of using declarative lazy-evaluated style in C++. Finally, reusable code polyhedral transformation patterns expressed as matcher/builder pairs can be shipped together with high-performance libraries that rely on generated code or with compilers for novel architectures.

Inria

Inria

Figure 1 :

 1 Figure 1: Schedule tree representation of the code in Listing 2.

Listing 3 :

 3 Schedule tree matchers declaratively describe the structure of a tree.

  isl::schedule_node node, continuation; auto matcher = band(node, isPermutable, anyTree(continuation)); auto builder = band([&](){ return tileSchedule(node); }, band([&]() { return pointSchedule(node); }, subtree(continuation))); replaceDFSPreorderOnce(schedule, matcher, builder);Listing 5: Declarative specification of rectangular loop tiling. Functions tileSchedule and pointSchedule divide and take modulo by tile sizes, respectively.4 Polyhedral Relation Matchers 4.1 General ConceptsLet us first introduce the concepts of access relation matchers before continuing with the matching algorithm and examples. Polyhedral relation matchers allow the caller to identify relations that have certain properties in a union of relations.

  std::vector<Candidate> candidates(isl::map relation, Pattern pattern); bool isValidMatch(std::vector<Placeholder>, std::vector<Candidate>);

  isl::ctx ctx = /*...*/ ; auto _i = placeholder(ctx); auto _j = placeholder(ctx); auto _A = arrayPlaceholder(); auto matcher = allOf(access(_A, _i, _j), access(_A, _j, _i)); auto result = match(readsAndWrites, matcher); some_call(result[_i]); // extracts the candidate for _i Listing 6: Access relation matcher for transposed accesses.

  isl::map transformMap(isl::map map, const Candidate& candidate, const Pattern& pattern);

  This can be expressed by using a relation matcher and a set of callback functions as shown in Listing 10. auto isGemmLike = [&](isl::schedule_node) { auto _i = placeholder(ctx); auto _j = placeholder(ctx); auto _k = placeholder(ctx); auto _A = arrayPlaceholder(); auto _B = arrayPlaceholder(); auto _C = arrayPlaceholder(); /* Restrict the access relations to this subtree */ reads = reads.intersect_domain(node.domain()); writes = writes.intersect_domain(node.domain()); auto mRead = allOf(access(_A, _i, _j), access(_B, _i, _k), access(_C, _k, _j)); auto mWrite = allOf(access(_A, _i, _j)); return match(reads, mRead).size() == 1 && match(writes, mWrite).size() == 1; }; auto matcher = band(_and(is3D, isPermutable, isGemmLike)), leaf()); Listing 10: Combined Schedule Tree and Access Matcher for a simple Tensor Contraction

  isl::schedule_node node = /* obtain node, e.g., from a matcher */ ; auto macroKernel = band([&]() { return tileSchedule(node, /*L2-specific sizes*/ ); }, band([&]() { return swapDims(pointSchedule(node), -2, -1); } )); Listing 11: Tiling for L2-cache locality. Negative indexes in swapDims have Python semantics.

/

  *Apply the previous builder and get the child band.*/ node = macroKernel.insertAt(node.cut()).child(0); auto microKernel = band([&]() { return tileSchedule(node, /*Vectorization sizes*/ ); }, band([&]() { return unrollAll(pointSchedule(node)); }));

Inria 6

 6 Related Work 6.1 Directive-based Transformation Engines Kelly et.al. proposed arguably the first framework to expose high-level transformation on top of the polyhedral model

  3 . These affine functions

	auto s = [&]() { return capturedBand
		.get_prefix_schedule_union_map(); };
	auto dlt = /* construct the DLT function */ ;
	auto from = dlt.set_tuple_id(isl::dim::range, "from");
	auto to	= dlt.set_tuple_id(isl::dim::range, "to");
	extension([&](){ return from.apply_domain(s()).unite(
		to.apply_domain(s()));},
	sequence(
	filter([&](){ return to.range();},
	band([&](){ return to.range().affine_hull()
		.wrapped_domain_map();}))
	filter([&](){ return s().domain();},
	subtree(capturedBand))
	filter([&](){ return from.range();},
	band([&](){ return from.range().affine_hull()
		.wrapped_domain_map();})));

  Listing 15: Schedule Tree Builder mimicing ppcg's behavior. same builder with different identifiers applied twice can be used to introduce thread and block identifiers. First-class support for builder combination is left for future work.

	auto builder =
	context(/*introduce block/thread ids*/ ,
	extension("{[...] -> copy_to[...];"
	" [...] -> copy_from[...];"
	" [...] -> dev_sync[]",
	sequence(
	filter("{copy_to[...]}"),
	filter([&]() { return node.domain(); },
	filter(/*block mapping*/ ,
	band([&]() { return tileSchedule(node); },
	filter(/*thread mapping*/ ,
	band([&]() { return pointSchedule(node); },
	subtree(child)))))),
	fitler("{copy_from[...]}"),
	filter("{dev_sync[]}"))));

RR n°9243

For example, all instances of S 1 (i) → (0, i) are scheduled before any instance of S

(i) → (1, i).RR n°9243

This can be ensured with, e.g., full/partial tile separation given tile size is equal to the number of vector lanes.

In the polyhedral model, subscripts should be affine.RR n°9243

GPU accelerator code is expressed in terms of one thread rather than parallel loops, so filters restrict the generated code to that one thread.RR n°9243

Inria

rue Simone Iff -CS 42112 75589 Paris Cedex 12 Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

This work was partially supported by the European Commission Horizon 2020 programme through the MNEMOSENE grant agreement, id. 780215 and the NeMeCo grant agreement, id. 676240 as well as through Polly Labs (Xilinx Inc, Facebook Inc, and ARM Holdings) and the Swiss National Science Foundation through the Ambizione program.


Pedro Magalhães

Xavier Cristol
email: xavier.cristol@fr.thalesgroup.com

Cornel Ioana
email: cornel.ioana@gipsa-lab.fr

Jerome Mars
email: jerome.mars@gipsa-lab.fr

Hausdorff distance: A new target localization using single hydrophone

Cet article traite du problème de la localisation de source avec un seul capteur en acoustique sous marine. Une nouvelle mesure de distance appelée Hausdorff Distance (HD) est introduite pour comparer différents ensemble de TDOA (Time difference of Arrival) simulés (TDOAs) aux TDOA reçus (TDOAr). Cette méthode utilisée en traitement d'image pour mesurer le degré de similitude entre deux objets est adaptée et appliquée ici dans un contexte complètement nouveau. Deux extensions (la distance maximum de Hausdorff et la distance médiane de Hausdorff) sont présentées sur un cas réel pour trouver la position correcte, (en portée et en profondeur) d'une cible. Les résultats en termes de localisation, de précision sont montrés dans une expérience réelle en cuve.

Introduction

In this paper, we investigate the problem of underwater source detection (target localization) and 2D localization, in range and deep, using a matching technique based on a new metric called Hausdorff Distance. Localization of acoustic sources is classically done by measurement of the time difference of arrival (TDOA) [START_REF] Sayed | Network-based wireless location: challenges faced in developing techniques for accurate wireless location information[END_REF]. Nowadays several solutions exist, based on cross correlation processing [START_REF] Porter | Shallow-water tracking in the Sea of Nazare[END_REF], Bartlett estimator [START_REF] Le Touze | Source localization on a single hydrophone[END_REF] as well as modal propagation.

The new metric Hausdorff Distance was very briefly used recently (2016) in underwater domain by Mours et al. [START_REF] Mours | Target-depth estimation in active sonar: Cramer-Rao bounds for a bilinear soundspeed profile[END_REF]. In our work, we show that, with only one sensor, it is possible to detect and give position based on TDOA analysis. In order to obtain these TDOA, we first estimate the time delay from cross correlation between received signal and transmitted signal, and; then, by finding the best match between real data and simulated data, it is possible to give target position in range and deep. The second section gives a basic understanding of underwater propagation in order to express the context (Ray path propagation). In Sec. 3, we present the experimental facility (tank) developed at GIPSA-LAB, and we explain how the experimental data have been recorded to search the target position. In Sec. 4, we will provide the description of Hausdorff Distance Technique and its extensions and the way it can be used for localization purpose. Finally in Sec. 5 we show the experiments results and the conclusion.

Underwater Propagation

Firstly, it is necessary to express the underwater propagation. Classically, propagation models can be classified into two groups [START_REF] Etter | Underwater Acoustic Modeling and Simulation[END_REF]:

First one is Range dependent model, where variation of environmental parameters (speed of sound and bathymetry) is considered not only as a function of depth, but also in terms of the range and azimuth, Second one is Range independent model, where horizontal stratification of the oceans is assumed, and then variation of environmental parameters is only a function of depth. In these geometries, because only one propagation path is considered, called one-way, there is no oceanographic characteristics that cause the way back (incoming wave) and therefore the solution is based only on the divergent wave (outgoing wave).

In this work, we ignore the way back (incoming wave) and we consider a range independent propagation because the variation of environmental parameters are disregarded.

Ray propagation

Now, we must define briefly the Raypath propagation model. Since early 1960s, it is used commonly for high frequencies and deep water, because it is only valid if the magnitude distances order involved are much greater than the wavelength [START_REF] Brekhovskikh | Fundamentals of Ocean Acoustics[END_REF]. The ray propagation theory is mainly described in [START_REF] Kuperman | Computational Acoustics[END_REF]. Figure 1 shows the variation of the Munk's sound speed profile for small variations of the ray's curvature, according to Snell's law. The received signal is composed of several arrivals that are time shifted and attenuated versions of the emitted signal. In this case, we will consider three different groups only, which differ by the reflection number at the bottom (figure 2) . The limit to 3 is due to high attenuation caused by those reflections at the bottom principally. On the first group -zero bottom reflection, second group -one bottom reflection and third group -two bottom reflections. For the first simulation (a 13 km range), it can be shown in the figure 3, a total of 10 multipaths : 2 for the first group (red color at 1 km range), 4 for the second (green color at 1 km range) and 4 for the third group (blue color at 1 km range). After 1km, only 8 multipaths are still presented due to the fact that the direct path does not appear anymore due to the ray curvature. It can be noted that, each variation in distance results in different TDOA, and the latter reduces with increasing distance. Another conclusion is that, for distance localization, the TDOA intergroup information are more significant, (being the reason whereby they have been used for localization in range). After this presentation, we give the objective of this paper: the TDOA of the received signal is compared with the set of simulated one via the Hausdorff Distance and through this process, we can identify the distance of the target (range and depth) with only one sensor.

For the second simulation, the figure 4 shows the total of 10 multipaths for all depths, 2 for the first group (red color at 1 km range), 4 for the second (light blue color at 1 km range) and third (dark blue color at 1 km range) group each. It can be noted that with increasing depth the TDOA intragroup become more significant due to the dispersion, being TDOA intragroup the reason whereby they have been used for depth localization.

The table 1 shows how the simulations were done and which parameters have considered for each case. To prove and give information on this new localization technique, a set of experiments have been conducted in the experimental facility of GIPSA-LAB at University of Grenoble Alpes. It is composed by a water tank whose size is 1.5 meters length by 1 meter width by 1 meter height, (shown in Figure 56). In our experiment, the signal was sampled at 25MHz with a resolution of 16 bits, and we considered the medium as homogeneous, used for approximation of the linear wave equation. The Doppler effect has not been considered. The sensors, either the transmitter, Tx, and the receiver, Rx, shown at "Fig 1", have 7 degrees' beam spread at 1MHz for 0.5-inch diameter. This allows us to disregard the tank's side considering a narrow beam spread.

The assembled system was static with 3 different distances, and same depths. We chose fixed distance and depth in order to show the localization results.

For the transmission, the sensor was located at 0.1 m deep. For the reception, the sensor was located at 0.3 m deep and 0.7 m range, centered horizontally in the tank (figure 5). The table 2 expresses the signal transmitted in this experiment. We were interested to obtain the impulse response from the

Hausdorff Distance Technique

Maximum Hausdorff Distance

As expressed, the comparison of TDOA was done by using the Hausdorff distance (HD), proposed by Huttenlocher [START_REF] Huttenlocher | Comparing images using the Hausdorff distance[END_REF]. This technique is quiet popular in image processing to measure the degree of similarity among different objects, giving an interesting measure of their mutual proximity. Our merit is to prove that matching technique performed through the minimization of the algorithm between two sets of point is efficient in underwater source localization.

Let consider two sets of points, first being the received TDOA R T = {r 1 , r 2 , r 3 , … , r x } , second the simulated TDOA S T = {s 1 , s 2 , s 3 , … , s y }, the Maximum Hausdorff distance between them is defined as:

        T T T T T T R , S h , S , R h max = S , R H (1)
where:

        x , … 1,2,3, = i min max = S , R h R r T T s r i   (2)         y , … 1,2,3, = j min max = R , S h S s T T r s i   (3)
Maximum Hausdorff Distance is small when every point in one of the sets is near to some point in the other. The Euclidean distance was used to evaluate the norm ‖ . ‖, and selects the farthest which is the largest discrepancy between the two sets of points or maximal distance.

Mean Hausdorff Distance

An extended version of the technique was done by Dubuisson, in 1994, called Modified Hausdorff Distance (MHD) with the Mean Hausdorff Distance [START_REF] Dubuisson | A modified Hausdorff distance for object matching[END_REF]. The main difference is that, all points contribute to measure the average of the distance, which ensures that more measurement points closely resemble the model. In other words this can soften the problems with false alarm or fake point at the received. MHD is given by:

        x , … 1,2,3, = i min X 1 = S , R h X 1 i T T     s r i R r (4)         y , … 1,2,3, = j min Y 1 = R , S h Y 1 j T T     r s i S s (5)

Results and Conclusion

In our experimental case, the final received signal is the combination of different signals with different takeoff angle for each ray path. This means that for each ray, we need to mechanically rotate the sensors and align with the correct angle in order to record the high energy of the transmitted signal, due to the high directivity of the sensor. The received package contains 9 different take off angles, shown at figure 7 Once the signal for each ray is recorded, we estimate the TDOA for each path. We first apply the sum in time domain of the signal (figure 7) and after we crosscorrelate it with the transmitted signal. The maximum local value for each ray is then estimated (figure 8). The arrival time for each ray (TDOAr) is compared with the TDOAs simulated using the ray path propagation simulator. We use the Hausdorff techniques to find the best matching in order to find the correct position. The signal noise ratio for this measurement is 23.69dB. In order to validate the results, we present, in figure 9 and 10, the HD and MHD results (normalized and with the same color mapping in order to direct compare both images in term of accuracy). The dark blue value is the minimum value. It can be considerered as the real position of the target. In both measurements, we can well locate the target (dark blue). In our experiment, we can also notice that the Mean Hausdorff Distance is more accurate, due to its small dark blue area, which is related to its variance of error, compared with the Maximum Hausdorff Distance. The target is located in the red line intersection, at 0.3 m deep and 0.7 m range. To conclude, we prove on this work, that it is possible localize target, with a good precision, using a Hausdorff approach (Classical or Mean) with only one sensor at the receiver. Both techniques have shown to be very robust for underwater localization. The different experimental setup made in a tank proved that is possible to localize real signals with good accuracy.

Figure 1 -Figure 2

 12 Figure 1 -Munk's sound speed profile

Figure 3 -

 3 Figure 3 -Arrival time with distance variation.

Figure 4 -

 4 Figure 4 -Arrival time with depth variation.

Figure 5 -

 5 Figure 5 -Measurement setup.

Figure 6 -

 6 Figure 6 -Experimental facility.

  . The amplitude of the received signal for each ray was combined with the 2 amplifiers (TX and RX). The synchronization of the different Transmitted signal type Time (s) Frequency (Hz) Chirp 100µs 500KHz-1.5MHz signals recorded, with different take off angles, was done by with transmitted signal.

Figure 7 -

 7 Figure 7 -TDOA recorded for 9 take off angles.

Figure 8 -

 8 Figure 8 -Received signal from sum of different take off angles.

Figure 9 -

 9 Figure 9 -Error of Maximum Hausdorff Distance

Table 1 -

 1 Set of simulations.

Table 2 -Transmitted signal.

 2


David Orive-Miguel
email: david.orive-miguel@cea.fr

Lionel Hervé

Jerome I Mars

Laurent Condat

Jean-Marc Dinten

Sophie Morales

Jérôme Mars

Computational Performance of Time-Resolved Diffuse Optical Tomography in a Two-Layer Brain Model

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Traumatic brain injury (TBI) is one of the leading causes of disability and injury-related deaths. Although magnetic resonance imaging or computerized x-ray tomography are standard techniques in TBI treatment protocols there is still a need of a technology that could allow to perform bedside, continuous and non-invasive brain monitoring. Timeresolved diffuse optical tomography (TR-DOT) is an emerging photonics technology that could potentially permit to continuously measure important brain physiological values in TBI cases such as blood oxygenation levels.

Even though TR-DOT has already been applied, for example, to flap viability assessment [START_REF] Di Sieno | Time-resolved diffuse optical tomography for non-invasive flap viability assessment: pre-clinical tests on rats[END_REF] with very promising results there is still a lack of accurate optical properties estimation for brain imaging cases. Nonlinear Monte-Carlo fitting methods have been used in computational analysis [START_REF] Selb | Comparison of a layered slab and an atlas head model for Monte Carlo fitting of time-domain near-infrared spectroscopy data of the adult head[END_REF] but their performance depends significantly on the knowledge of the geometry being scanned and is time consuming. In this work we computationally evaluate an algorithm based on Mellin-Laplace moments and we analyze its accuracy and adaptation to different layer thicknesses.

Method

The computational phantom consists of a two-layer model (8×8×5 cm 3 ) containing an extracerebral and cerebral layer with one light source and four detectors placed at 1 cm distance at the top boundary, see Figure 1a. Monte-Carlo simulations [START_REF] Fang | Mesh-based Monte Carlo method using fast ray-tracing in Plcker coordinates[END_REF] were performed on the phantom, signals were convolved with the instrumental response function (IRF) of single-photon avalanche diode (SPAD) and silicon photomultipliers (SiPM) detectors, and were corrupted with Poisson noise.

For the reconstruction of the optical properties, Mellin-Laplace moments algorithm [START_REF] Zouaoui | Quantification in time-domain diffuse optical tomography using Mellin-Laplace transforms[END_REF] with twelve moments and p = 9 was used. This algorithm is based on finite element method and computes the moments of time-point spread functions in order to achieve a good balance between accuracy and computational time. The algorithm can be used on experimental data by performing measurements in a reference phantom; moreover, knowledge of the shape of the IRF is not necessary.

Results and Discussion

The first experiment consists in evaluating the Mellin-Laplace algorithm performance for estimating the absorption at cerebral layer. Monte-Carlo simulations were done with a constant absorption of µ a,1 = 0.1 cm -1 at the extracerebral layer and an increasing absorption from µ a,2 = 0.1 to 0.5 cm -1 on steps of 0.05 cm -1 for the cerebral layer. The scattering coefficient was maintained homogeneous along both layers because of its little sensitivity at brain layer [START_REF] Selb | Comparison of a layered slab and an atlas head model for Monte Carlo fitting of time-domain near-infrared spectroscopy data of the adult head[END_REF]. Extracerebral and cerebral layer thicknesses were 1 cm and 4 cm respectively; homogeneous absorption through each layer was assumed. An homogeneous phantom of µ a = 0.1 cm -1 was taken as reference measurement. Results are plotted in Figure 1b for different types of detectors; they show that Mellin-Laplace algorithm reconstructs linear absorption estimations with low relative error. It is also seen that IRFs influence the results by yielding an overestimation at absorption values much higher than the reference measurement making the relative error to increase twice with respect to free IRF data. No significant differences are seen in the performance with SPAD or SiPM detectors.

In the second experiment we aim to evaluate whether the algorithm is able to detect the thickness of the two layers; that is, whether the algorithm has the flexibility to adapt to the two-layer geometry when thickness of each layer is not given in advance. The domain was segmented in 0.3 cm-thick layers each of them with an homogeneous absorption. The two-layer geometry interface was estimated as the layer interface with largest absorption difference between above and below layers. Monte-Carlo simulations were performed for cases were extracerebral layer thickness ranged from 0.5 to 1.1 cm on steps of 0.3 cm; cerebral layer absorption ranged from µ a,2 = 0.2 to 0.5 cm -1 on steps of 0.1 cm -1 . Results in Figure 1c show that accurate extracerebral layer thickness is estimated in 0.5 and 0.8 cm cases. Nevertheless, thickness is underestimated in 1.1 cm case because the algorithm has less sensibility to deeper layers. Higher absorption values tend to underestimate the thickness because the algorithm tries to compensate the high absorption by decreasing the extracerebral layer thickness.

Conclusions and Future Work

We studied the performances of a Mellin-Laplace algorithm for two-layer brain models with data that mimics realistic conditions. Accurate results are obtained if each layer thickness is given in advance although the algorithm is also capable of doing good estimations without that a priori information.

In the near future, we will extend the analysis by performing Monte-Carlo simulations on ATLAS head models and doing reconstructions of the absorption at the brain layer.

Acknowledgements

This project has received funding from the European Union's Horizon 2020 Marie Skodowska-Curie Innovative Training Networks (ITN-ETN) programme, under grant agreement no 675332 BitMap. The authors thank to the photonics lab at Politecnico de Milano for providing the instrumental response functions of the detectors.

  Absorption reconstruction for simulated data with and without IRF. Confidence interval of 75%. Layer thickness estimation without using a priori data. Confidence interval of 75%.

Fig. 1 :

 1 Fig. 1: Results using Mellin-Laplace moments algorithm.