Fabian Gruber

Manuel Selva

Diogo Sampaio

Christophe Guillon

Antoine Moynault

Fabrice Rastello

Data-Flow/Dependence Profiling for Structured Transformations

Keywords: Performance Feedback, Polyhedral Model, Loop Transformations, Compiler Optimization, Binary, Instrumentation, Dynamic Dependence Graph

Profiling feedback is an important technique used by developers for performance debugging, where it is usually used to pinpoint performance bottlenecks and also to find optimization opportunities. Assessing the validity and potential benefit of a program transformation requires accurate knowledge of the data flow and dependencies, which can be uncovered by profiling a particular execution of the program.

In this work we develop poly-prof, an end-to-end infrastructure for dynamic binary analysis, which produces feedback about the potential to apply complex program rescheduling. Our tool can handle both inter-and intraprocedural aspects of the program in a unified way, thus providing interprocedural transformation feedback.

Introduction

The most effective program transformations for improving performance or energy consumption are typically based on rescheduling instructions so as to expose data locality and/or parallelism. Optimizing compilers typically attempt, via static analysis, to build a representation of the program precise enough to enable useful program transformations. A key issue faced when analyzing general-purpose languages is the ambiguities introduced by the language itself: for example the use of pointers typically restricts the ability to precisely characterize the data being accessed, in turn triggering conservative and approximate dependence analysis [START_REF] Griffith | GCC: the complete reference[END_REF][START_REF] Lattner | LLVM: A Compilation Framework for Lifelong Program Analysis & Transformation[END_REF][START_REF] Trifunovic | GRAPHITE Two Years After: First Lessons Learned From Real-World Polyhedral Compilation[END_REF]. Such frameworks therefore rely on conservative assumptions, limiting the ability to reason on the legality of complex program restructuring. In this work we specifically target structured transformations, that is, a (possibly interprocedural) reordering of operations involving complex sequences of multidimensional loop transformations, including (partial) loop fusion/fission, loop tiling [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF], etc.

When a region of the source program fits specific syntactic and semantics restrictions, such as avoiding function calls, using only arrays as data structures, and very simple conditional statements with no indirections [START_REF] Collard | Fuzzy array dataflow analysis[END_REF][START_REF] Doerfert | Spolly: speculative optimizations in the polyhedral model[END_REF], transformation frameworks such as the polyhedral model [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF] showed that multidimensional loop nest rescheduling can be successfully implemented, leading to significantly improved performance [START_REF] Grosser | Polly -Performing polyhedral optimizations on a low-level intermediate representation[END_REF][START_REF] Pop | GRAPHITE: Polyhedral analyses and optimizations for GCC[END_REF][START_REF] Trifunovic | GRAPHITE Two Years After: First Lessons Learned From Real-World Polyhedral Compilation[END_REF]. The input program can be massaged manually via function inlining, loop rewriting/normalization, etc. to enable such static analyses to succeed, but this is a rare scenario: In full programs and in particular those relying on external binaries visible only to the linker in compiled form, often data allocation and even the full call stack is inaccessible to static analysis.

Dynamic analysis frameworks [START_REF] Beyls | Discovery of Localityimproving Refactorings by Reuse Path Analysis[END_REF][START_REF] Faxén | Embla-data dependence profiling for parallel programming[END_REF][START_REF] Manuel | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF][START_REF] Simbürger | PolyJIT: Polyhedral Optimization Just in Time[END_REF] address this limitation by reasoning on a particular execution of the program, through the analysis of its execution trace. That is, disambiguation of dependence information and accessed data is addressed by monitoring the stream of operations and accessed addresses. In general, the result is only valid for that particular execution. Such systems provide some feedback to the user, pinpointing the probable absence of dependencies along some loop, thus highlighting potential parallelism [START_REF] Ketterlin | Profiling Data-Dependence to Assist Parallelization: Framework, Scope, and Optimization[END_REF][START_REF] Tournavitis | Semi-automatic extraction and exploitation of hierarchical pipeline parallelism using profiling information[END_REF][START_REF] Vandierendonck | The Paralax infrastructure: automatic parallelization with a helping hand[END_REF][START_REF] Wang | Integrating profile-driven parallelism detection and machine-learning-based mapping[END_REF] or SIMD vectorization [START_REF] Holewinski | Dynamic trace-based analysis of vectorization potential of applications[END_REF].

The contributions of this work are: 1. The development of an inter-procedural intermediate representation that compacts program traces effectively while being amenable to polyhedral optimization, and all associated algorithms and implementation to make the process scalable to full programs. 2. poly-prof, a tool that provides optimization feedback using binary translation and instrumentation. poly-prof is the first framework that provides feedback on structured transformation potential for arbitrary binaries.

Overview of poly-prof

In contrast to, for instance, basic loop parallelization, structured transformations can be found even in the presence of data dependencies between loop iterations, dramatically broadening the class of loop transformations that can be applied on the program. This, however, requires detailed information about such dependencies. That is, instead of just proving the absence of dependencies one must capture and represent them in order to be able to find patterns and reason about their structure.

An overview of poly-prof is shown in Fig. 1. Taking a carefully generated execution trace of the program, the first objective of poly-prof is to construct a useful and analyzable representation of the trace. Building such a representation constitutes our first key set of contributions, and covers the first three stages below. We then extended a polyhedral compiler. The goal here is not to provide automated transformations, but instead to assist the user in figuring out where optimizations should be implemented, and importantly which ones. This corresponds to the fourth stage below that forms the second contribution of this work. Interprocedural control structure To eventually reconstruct a compact and useful representation of the program trace, the first stage starts by dynamically computing the intraprocedural control flow graphs (CFGs) and the interprocedural call graph (CG) from an optimized binary. This information is then used to construct the interprocedural loop context tree which combines the notions of loop-nesting-trees and their equivalent for call graphs, the recursive-componentset, both of which are input to the next stage. Their structure and construction are detailed in section 3. Dynamic Dependence Graph This stage generates the actual dynamic dependence graph (DDG) [START_REF] Nethercote | Redux: A dynamic dataflow tracer[END_REF], a trace representing both the dynamic execution of a program's instructions as well as their data dependencies. It uses the interprocedural loop context tree to construct dynamic interprocedural iteration vectors, a new interprocedural representation of the program execution which unifies two abstractions: the intraprocedural schedule tree [START_REF] Kelly | A unifying framework for iteration reordering transformations[END_REF][START_REF] Verdoolaege | Schedule trees[END_REF] and the interprocedural calling-context-tree [START_REF] Ammons | Exploiting Hardware Performance Counters with Flow and Context Sensitive Profiling[END_REF]. Both the DDG and the dynamic interprocedural iteration vector are described in section 4. Compact polyhedral DDG The third stage compacts, or folds, the DDG into a polyhedral program representation where, essentially, sequences of individual operations in the trace are merged together into sets of points. Intuitively, when a loop is executed it generates one execution of a statement for each loop iteration and therefore one point in the trace for each execution. This stage amounts to folding those points back into a loop-equivalent representation, which itself can be the input of a polyhedral compiler. However we proceed in a general case, where folding occurs across multi-level loops but also, possibly recursive, procedure calls, leading to interprocedural folding. Details are presented in [START_REF] Gruber | Building of a Polyhedral Representation from an Instrumented Execution: Making Dynamic Analyses of non-Affine Programs Scalable[END_REF], and an outline of the process is presented in Sec. 5. DDG polyhedral feedback The last stage of poly-prof consists of a customized state-of-the-art polyhedral toolchain, for the purpose of analyzing very large polyhedral programs. It analyzes the polyhedral DDG and suggests sequences of program transformations needed to expose multilevel parallelism and/or improved data locality, along with numerous statistics on the original and potentially transformed program including SIMDization potential. Section 6 describes this last step as well as how poly-prof provides human interpretable feedback to the user.

Section 7 illustrates this feedback through a few case studies. Section 8 thoroughly evaluates poly-prof on the Rodinia benchmarking suite. Related work is discussed in Section 9 before concluding.

Interprocedural Control Structure

As most of the execution time of a program is usually spent in loops, poly-prof has been designed to find loop-based transformations. In practice, however, as illustrated in Ex. 1 of Fig. 3, interesting loop nests are often spread across multiple functions over long call chains, obfuscating them to traditional static analysis. poly-prof is designed to be able to represent both classical loops as well as function calls and returns in a unified way with polyhedra [START_REF] Feautrier | Polyhedron model[END_REF][START_REF] Michael | A data locality optimizing algorithm[END_REF].

DDG Polyhedral Feedback

poly-prof has been designed to handle several scenarios. 1. Nested loops containing function calls with side effects: here calls can be viewed as atomic instructions, but profiling the storage locations they access or modify and compute the dependencies between each call and the rest of the loop body is required. 2. Loops containing function calls that themselves contain loops, as in the Ex. 1 of Fig. 3: here the objective is to view the whole interprocedural region as a multidimensional loop nest. 3. Recursive functions, as in Ex. 2 of Fig. 3: here the primary objective is to avoid the depth of the generated representation to grow proportionally with the depth of the call stack; A secondary objective is to detect any underlying regular loop structure amenable to polyhedral compilation (after recursion-to-iteration conversion).

This section explains how to go from the execution of a program to a representation of its interprocedural control structure.

The first step of poly-prof is to extract the inter-and intraprocedural static control structure from a program in compiled binary form. This static control structure of the program will then be used in a later stage of poly-prof to transform the stream of raw events (jump/call/return) gathered during execution into loop events (entry/iterate/exit).

poly-prof performs this through instrumentation of jump, call, return and a few other instructions to dynamically trace any transfer of control during execution. It thus "dynamically" extracts the control-flow-graph (CFG) of each function and then proceeds to statically build the loop-nestingtree. It also dynamically extracts the call-graph (CG) of the whole program and builds the counterpart of loop-nestingtrees for call-graphs, the recursive-component-set. An advantage of our approach is that only the part of a program that is actually executed will be analyzed. For large programs with only a small performance-critical part this can substantially reduce the amount of work for the analyzer.

Control-flow-graph and loop-nesting-tree

For the profiled CFG, the loop detection algorithm used by poly-prof is the one described by Havlak [START_REF] Havlak | Nesting of Reducible and Irreducible Loops[END_REF]. As formalized by Ramalingam in [START_REF] Ramalingam | On Loops, Dominators, and Dominance Frontiers[END_REF], a loop-nesting-forest can be recursively defined (the actual algorithm is almost linear) as follows: 1. Each strongly connected component (SCC) of the CFG containing at least one cycle constitutes the region of an outermost loop; 2. for each loop, one of its entry nodes is designated the header node of the loop; 3. all edges within the loop that target its header node are called back-edges; 4. "removing" all back-edges from the CFG allows one to recursively define the next level of sub-loops.

An example CFG and its corresponding loop-nesting-tree is given in Figures 2a &2b. Here, the CFG can be partitioned into one SCC (giving rise to loop L 1 with B as its header) plus two separated basic-blocks A and E. The region of L 1 , once its back-edge (D, B) has been removed, contains one sub-SCC formed by the loop L 2 and the basic-block B. Among the two entry nodes C and D of L 2 only one, C, is selected to be its header. As depicted in Alg. 1, those notions (region, header, back-edge) are important as they allow associating loop events (E: entry, I: iterate, and X: exit) with the stream of raw control events (jump, call, return) that are produced by our instrumented execution. Here, generation of loop events is driven by the control event "jump" that is emitted by our instrumentation. Whenever we detect the start of a new iteration of a given loop (Line 8), all live sub-loops are considered exited, that is, we emit an "exit" event (Line 4). This is especially important for recursive loops as further explained in Alg. 2.

Call-graph and recursive-component-set

To uniquely identify each dynamic execution of an instruction, also called a dynamic instruction, poly-prof uses interprocedural iteration vectors. This is described in more detail in Section 4. Note, that the modeling of programs containing recursive function calls with calling-context-paths is Algorithm 1 CFG-loop events generated from a jump event.

A B Header L1 C Header L2 D E Back-edge L2 Back-edge L1 L 1 L 2 Return (a) r A L 1 :B E B L 2 :C C D (b) L 1 Main B C D E (c) components = {L 1 } L 1 .entries = {B} L 1 .headers = {B, C} (d)
Input:

• event: Branch event.

• inLoops: memory inefficient since the length of the paths is proportional to the depth of the recursion. While recursive function calls are found to be very uncommon in performance critical code 1 , we do need to handle them to ensure robustness. poly-prof handles recursiveness in an elegant way that we believe to be useful beyond the restricted scope of this paper, for example to detect properties of tree-recursive calls. The main data structure for treating recursive control flow is the recursive-component-set which is for the call-graph what the loop-nesting-tree is for the control-flow-graph. Before providing more details, we go through the illustrating Ex. 2 of Fig. 3. Here, the edge from B 3 to B 0 (in orange) is not a CFG-edge but a recursive call to B from call site B 3 . The recursive-component-set computed from the CG contains one SCC with a cycle made up of a single function B. This example raises several questions: 1. Should C 0 be part of the associated recursive loop? It actually depends on the context: It should, when called from B 1 , while it should not when called from D 0 . 2. What about B 5 ? One should observe that B 5 will be executed as many times as the number of recursive calls to B: In other words, B 5 should be part of a loop.

To conclude, while CG-cycles clearly represent potential dynamic loop structures, a CG-edge does not have the same semantic as a CFG-edge. In particular, a call will never exit a recursive loop. For poly-prof, the recursive loop "L 1 " of Ex. 2 is a dynamic notion defined as follows: 1. Entering L 1 is characterized by a call to B (step 1). 2. L 1 is exited when this (first) call gets unstacked, that is, when the number of returns reaches the number of calls (step 22). 3. Everything executed in between is part of the loop and iteration, and the corresponding increment of induction variables takes place whenever a call/return to/from B occurs (steps 10,15,20,and 21). As one can see, once the recursive-component-set has been computed from the CG, the only relevant information dynamically used to devise recursive-loop events corresponds to the header functions, here B.

As already stated, the recursive-component-set is for the CG what the loop-nesting-tree is for the CFG. In a similar way it can be recursively defined as follows:

1. Find all the top-level SCCs with at least one cycle in the CG. Each gives rise to a distinct recursive-component. 2. For each component, determine all its entry nodes. 3. Repeat the following phases until no more cycles exist: a. For a given SCC, choose an entry node and add it to the headers-set of the recursive-component (top-level SCC) it belongs to. b. Remove all edges inside this SCC that point to this node. Havlak's loop-nesting-forest construction algorithm can easily be adapted to build the recursive-component-set in almost linear time. The end result of the algorithm is a possibly empty set of recursive-components where each recursivecomponent has a non-empty set of headers plus a non-empty set of entries associated with it. Alg. 2 uses this data structure to associate loop events (entry, iterate, and exit) to control events (call and return). Here: 1. entering a recursive loop is characterized by a call to a component's entry function (Line 4); 2. a new iteration is started whenever a call/return to/from one of the components' header occurs (Line 6); 3. an exit occurs when all the iterating calls to the headers have been unstacked, that is, when the number of returns equals the number of calls, and we are returning from the original function that entered the loop (Line 18). Tracking the state of the call stack is done with the following two data structures: L.stackcount represents for a recursive-component L a counter of the number of calls-to a header minus the number of returns from it; L.entry represents the function through which L was entered.

DDG: Dynamic dependence graph

The objective of the second stage of poly-prof ("Instrumentation II") is to profile the dynamic dependence graph of a given execution, that is, to build a graph that has one vertex for every dynamic instruction and one edge for every data dependence. Because we want to enable feedback with structured loop transformations, we need to map this graph to a "geometric" space that reflects the structural properties of the program. To this end, we tag each dynamic instruction with its iteration vector (IV). The IVs uniquely identify each Algorithm 2 Different recursive-loop events generated from a call or a return event.

Input:

• event, inLoops: same as for Alg. 1 Emitted events: • E C (L, B): call-to a function header of recursive-component L and entry into the corresponding loop. B is the current basic-block after the call. dynamic instruction and naturally span a geometric space. A data dependency is then simply represented as the pair of the IVs of the producer and the consumer.

• I C (L, B) / I R (L, B): call-to / return-
To handle interprocedural programs we also need a notion of calling context that is scalable in the presence of recursive calls. Our dynamic interprocedural iteration vector (dynamic IIV) described in this section addresses those objectives by unifying two notions: 1. Kelly's mapping which describes intraprocedural IVs and is used by the polyhedral framework [START_REF] Kelly | A unifying framework for iteration reordering transformations[END_REF]; 2. calling-context-paths used by profiling feedback tools. We first briefly recall those two notions. Kelly's mapping For a given function, Kelly's mapping can be explained using a form of schedule tree [START_REF] Verdoolaege | Schedule trees[END_REF] as shown in Fig. 4. Here a schedule tree is nothing else than a decorated loop-nesting-forest. The first decoration consists of associating a "static" index to each loop and basic-block: Recall the recursive characterization of loops given by Ramalingam in the previous section. For a given loop region (e.g. L j in the schedule-tree of the fused version in Fig. 4), its sub-regions (once back-edges have been removed -here statements S and T ) form a directed-acyclic-graph (reduced DAG represented in dashed in Fig. 4) that can be numbered using a topological order. This numbering is used to index the corresponding nodes at every level of the loop-nesting-tree (e.g. S(0) and T (1) for the fused schedule or L i (0) and L i ′ (1) for the fissioned one). The second decoration consists of associating a canonical induction variable, that is, an induction variable that starts at value 0 and increments by 1, to each loop-vertex. As an example the loop-vertex associated with L j is labeled

A 0 A 1 A 2 A 3 B B B B M 0 Main A B L 1 L 2 (a) Example Main A B (b) r(Main) M 0 r(A) r(B) B 0 L 2 :B 1 B 3 A 1 A 2 B 1 B 2 A 0 L 1 :A 1 A 3 (c) step : trace event dynamic IIV (⊥) 1 : M 0 N(M 0 ) (M 0 ) 2 : A 0 C(A 0 ) (M 0 /A 0 ) 3 : A 1 E(L 1 , A 1 ) (M 0 /L 1 , 0, A 1 ) 4 : B 0 C(B 0 ) (M 0 /L 1 , 0, A 1 /B 0 ) 5 : B 1 E(L 2 , B 1 ) (M 0 /L 1 , 0, A 1 /L 2 , 0, B 1 ) 6 : B 2 N(B 2 ) (M 0 /L 1 , 0, A 1 /L 2 , 0, B 2 ) 7 : B 1 I(L 2 , B 1 ) (M 0 /L 1 , 0, A 1 /L 2 , 1, B 1 ) 8 : B 2 N(B 2 ) (M 0 /L 1 , 0, A 1 /L 2 , 1, B 2 ) 9 : B 3 X(L 2 , B 3 ) (M 0 /L 1 , 0, A 1 /B 3 ) 10 : A 2 R(A 2 ) (M 0 /L 1 , 0, A 2 ) 11 : A 1 I(L 1 , A 1 ) (M 0 /L 1 , 1, A 1 ) . . . . . . . . . (d) L 1 (1), (i 1 ) M 0 (0) r A 3 (2) A 0 (0) A 1 (0) A 1 (0) A 2 (1) L 2 (1), (i 2 ) B 3 (2) B 0 (0) B 1 (0) B 2 (1)
(e) with L j as well as the static index 0 followed by its canonical variables j, resulting in L j (0), (j). For any given statement, an IV with Kelly's mapping is nothing else than the vector given by the path from the root to its corresponding leaf, called its iteration vector. Fig. 4 shows this mapping both in its textual form using region names and numerical form using indices. As one can see, an interesting property of the numerical form of this mapping is that the scheduling of the original code is given by the lexicographical order of the so obtained iteration vectors (unique per dynamic instance).

M 0 M 1 Main B 0 B 1 B 2 B 3 B D 0 D C 0 C B 4 B 5 (f) Example 2 Main B C D L 1 (g) Main B(M 1 ) D(M 0 ) C(D 0 ) C(B 1 ) B(B 3 ) C(B 1 ) B(B 3 ) C(B 1 ) B(B 3 ) (h) step : trace call-stack event dynamic IIV M (⊥) 1 : M 0 M N(M 0 ) (M 0 ) 2 : D 0 M /D C(D 0 ) (M 0 /D 0 ) 3 : C 0 M /D/C C(C 0 ) (M 0 /D 0 /C 0 ) 4 : M 1 M R 2 (M 1 ) (M 1 ) 5 : B 0 M /B E C (L 1 , B 0 ) (M 1 /L 1 , 0, B 0 ) 6 : B 1 M /B N(B 1 ) (M 1 /L 1 , 0, B 1 ) 7 : C 0 M /B/C C(C 0 ) (M 1 /L 1 , 0, B 1 /C 0 ) 8 : B 2 M /B R(B 2 ) (M 1 /L 1 , 0, B 2 ) 9 : B 3 M /B N(B 3 ) (M 1 /L 1 , 0, B 3 ) 10 : B 0 M /B/B I C (L 1 , B 0 ) (M 1 /L 1 , 1, B 0 ) 11 : B 1 M /B/B N(B 1 ) (M 1 /L 1 , 1, B 1 ) 12 : C 0 M /B/B/C C(C 0 ) (M 1 /L 1 , 1, B 1 /C 0 ) 13 : B 2 M /B/B R(B 2 ) (M 1 /L 1 , 1, B 2 ) 14 : B 3 M /B/B N(B 0 ) (M 1 /L 1 , 1, B 3 ) 15 : B 0 M /B/B/B I C (L 1 , B 0 ) (M 1 /L 1 , 2, B 0 ) 16 : B 1 M /B/B/B N(B 1 ) (M 1 /L 1 , 2, B 1 ) 17 : C 0 M /B/B/B/C C(C 0 ) (M 1 /L 1 , 2, B 1 /C 0 ) 18 : B 2 M /B/B/B R(B 2 ) (M 1 /L 1 , 2, B 2 ) 19 : B 4 M /B/B/B N(B 4 ) (M 1 /L 1 , 2, B 4 ) 20 : B 5 M /B/B I R (L 1 , B 5 ) (M 1 /L 1 , 3, B 5 ) 21 : B 5 M /B I R (L 1 , B 5 ) (M 1 /L 1 , 4, B 5 ) 22 : ⊥ X R 2 (L 1 , ⊥) (⊥) (i) M 0 (0) r D 0 (0) C 0 (0) B 1 (1) B 3 (3) B 0 (0) B 1 (1) B 2 (2) M 0 (0) M 1 (0) M 1 (0) L 1 (0), ( i
Calling-context-tree Differently from the schedule tree, the calling-context-tree [START_REF] Ammons | Exploiting Hardware Performance Counters with Flow and Context Sensitive Profiling[END_REF] (CCT) is only enumerative (it does not contain any loop indices) and reflects a dynamic behavior. In other words, it encodes the dynamic-call-tree in a "compact" way. The calling-context-tree of the example in Fig. 3f is reported in Fig. 3h. This figure is slightly different from the "original" CCT, but corresponds to the current practice: to differentiate two calls to a common function from different basic-blocks, callees are labeled with call sites (in purple and under parenthesis in Fig. 3h). In this example, a calling-context-path (for example M 1 /B 3 . . . B 3 /B 1 /C) can be as long as the number of recursive calls to B, but the calling context is fully encoded making it possible to differentiate the different contexts within which basic-block C 0 is executed. Obviously, looking at this example, one wants to fold all the repeated calls from B 1 to C.

Dynamic IIV The dynamic IIV is basically a combination of Kelly's mapping and the CCT. Similarly to Kelly's mapping, the dynamic IIV alternates between context-ids and canonical induction variables. Differently from Kelly's mapping, but more like the CCT, each context-id in a dynamic IIV is a, possibly empty, stack of calling contexts and the identifier for a basic block.

Examples of IIVs are shown in Fig. 3a: Function A contains a loop L 1 that contains a call to function B, itself containing a loop L 2 . In this interprocedural example one clearly wants to see that basic-block B 1 in loop L 2 belongs to a twodimensional nested loop. This is reflected by our dynamic IIV (see steps 5,7 of Fig. 3d) which will be

(M 0 /L 1 , i 1 , A 1 /L 2 , i 2 , B 1 )
where i 1 (resp. i 2 ) are the canonical induction variables of loop L 1 (resp. L 2 ). Here, our context-ids are loop-ids (e.g. L 2 ) or statement-ids (e.g. B 0 ) and each context-id is decorated with a, possibly empty, call stack (e.g. A 1 ).

1 for (i =0; i <n;i ++) 2 for (j =0; j <= i ; j ++) 3

{ S; T } 1 for ( i =0; i < n ; i ++) 2 for ( j =0; j <= i ; j ++) 3 { S } 4 for ( i ' =0; i ' <n ; i ' ++) 5 for ( j ' =0; j ' <= i ' ; j ' ++) 6

{ T } (a) Nested loop before and after fission L i (0), (i)

L j (0), (j)

S(0) T(1) r L i (0), (i)
L j (0), (j) Similar to Kelly's mapping one can also construct a dynamic schedule tree from the dynamic IIVs of a program execution. Also, note that the schedule tree is for the dynamic IIVs what the calling-context-tree is for the set of calling context paths. The schedule tree for our example is shown in Fig. 3e. The relationship between these three structures is summarized in Fig. 5. As described in more details further, poly-prof exposes the schedule tree to the programmer in the form of a flame-graph [START_REF] Gregg | The flame graph[END_REF] where the root of the tree is on the bottom and leaves on the top.

S(0) T(0) r L j' (0), (j') L i' (1), (i') (b) Corresponding schedule trees (reduced DAG in dashed) S → [L i , i, L j , j, S ] → [0, i, 0, j, 0] T → [L i , i, L j , j, T ] → [0, i, 0, j, 1] S → [L i , i, L j , j, S ] → [0, i, 0, j, 0] T → [L i ′ , i ′ , L j ′ , j ′ , T ] → [1, i ′ , 0, j ′ , 0] (
Ex. 3f illustrates how the recursive-component-set is used to "fold" calling-context-paths in the presence of recursive calls. Here, M 1 is the caller site (step 4) before entering the recursive loop L 1 (step 5). Looking at C 0 , its multiple instances are indexed by the corresponding recursive-loop induction variable i 1 (respectively 0,1, and 2 at steps 7, 12, and 17). Within the loop it gets executed when C is called from B 1 . The associated dynamic IIV is (M 1 /L 1 , i 1 , B 1 /C 0 ). As already mentioned, B 5 is also part of the loop: Indeed there are as many instances of B 5 as there are calls to B from B 3 (2 calls at steps 10,15). Observe that the value of i 1 does not reflect the size of the call stack: It does not go up and down. It keeps increasing. The main reason for doing so is related to our polyhedral back-end: 1. we want our indexing to be lexicographically increasing; 2. we want to match to the polyhedral formalism where iterators can be represented using simple canonical induction variables. To do so, any return (except the last that exits the recursive loop) associated with a call to B (steps 20,21) leads to incrementing i 1 .

As previously described, loop events allow updating the dynamic IIV. This is done as detailed in Alg. 3. Here, diiv corresponds to the dynamic IIV, where the rightmost is the innermost dimension. Each dimension contains two fields, the induction variable IV, followed by the context variable CTX: As an example, I(L 1 , A 1 ) applied to (M 0 /L 1 , 0, A 2 ) increments the IV of L 1 and updates the CTX to A 1 , leading to (M 0 /L 1 , 1, A 1 ). CTX.last represents the last element of the context variable: As an example, if CTX = M 0 /D 0 /C 0 , CTX.last = C 0 . CTX.push(B) appends B to the context variable: As an example C(C 0 ) applied to (M 1 /L 1 , 0, B 1 ) leads to (M 1 /L 1 , 0, B 1 /C 0 ). CTX.pop() does the reverse: As an example R(M 1 ) applied to (M 0 /D 0 ) pops D 0 and updates the last element to M 1 leading to a CTX of (M 1 ). addDimension(0, B) adds those two fields to the vector: As an example E C (L 1 , B 0 ) applied to (M 1 ), appends L 1 to CTX and adds a dimension with innermost context set to B 0 , leading to (M 1 /L 1 , 0, B 0 ). removeDimension() does the reverse: For example, X(L 2 , B 3 )

applied to (M 0 /L 1 , 0, A 1 /L 2 , 1, B 2 ) leads to (M 0 /L 1 , 0, A 1 /B 3 ).

Algorithm 3 Updating of dynamic IIV (diiv)

Input:

• event: Branch event.

• diiv: Current dynamic IIV. Output:

• Updated dynamic IIV 1 if event == C ( B ) or E C (L , B ): 2 diiv . innermost . CTX . push ( L ) 3 if event == E (L , B ) or E C (L ,
I1 → I2 I2 → I4 I4 → I4 IV Label IV Label IV Label (cj,ck) (cj',ck') (cj,ck) (cj',ck') (cj,ck) (cj',ck') (0,0) (0,0) (0,0) (0,0) (0,1)
(0,1) (0,1) (0,1) (0,1) (0,0) . . . . . . . . . . . . . . . . . . Table 1. Dependency input stream from example in Fig. 6 main properties are outlined here. Our folding algorithm borrows ideas from trace compression [START_REF] Ketterlin | Prediction and Trace Compression of Data Access Addresses Through Nested Loop Recognition[END_REF] with the notable difference that is based on a geometric approach and is able to perform over-approximations as briefly described below.

Folding interface The output of the second stage of polyprof (Instrumentation II) that feeds the folding stage is a stream of IIVs and labels. The IIV representation used at this stage slightly differs from the previous section. IIVs are split into two parts: 1. The context corresponds to the non-numerical part of the vector; 2. The coordinates correspond to the numerical part. Folding is then performed for each context separately. The folding algorithm takes as input for each dynamic instruction and each dynamic dependence: 1. ì I : its iteration vector, which uniquely identifies it (i.e., its unique coordinates). 2. Label: a vector a( ì I ) of associated integer values. And as an output it produces: 1. a union of polyhedra that represent the set of all ì I . 2. for each polyhedron P, an affine function A such that for all ì I ∈ P, A( ì I ) = a( ì I ). For dynamic instructions, a( ì I ) are the integer and pointer values produced by the instruction, if any. For dependencies, a( ì I ) are the IV of the producer instruction. Dependency recognition Some of the streams compressed by the folding stage are dependencies. To illustrate its functionality on dependencies, consider bpnn_layerforward, a kernel of the backprop benchmark from the Rodinia benchmark suite. Part of the input stream and corresponding folded output dependencies can be seen in Tables 1 and2. Note that iterators used by poly-prof are canonical ones (here cj and ck) computed on the fly through program instrumentation, that do not necessarily match the ones present in the original code (e.g. j of the outer loop starts at one). SCEV recognition At the machine code level, even very regular programs contain a large amount of "unimportant" code, such as instructions that increment loop counters or calculate addresses relative to a base pointer. Those are identified by compilers as scalar evolution expressions (SCEVs) [

I1 → I2 0 ≤ c j ≤ 15, 0 ≤ ck ≤ 42 c j ′ = c j + 0ck, ck ′ = 0c j + ck I2 → I4 0 ≤ c j ≤ 15, 0 ≤ ck ≤ 42 c j ′ = c j + 0ck, ck ′ = 0c j + ck I4 → I4 0 ≤ c j ≤ 15, 1 ≤ ck ≤ 42 c j ′ = c j + 0ck, ck ′ = 0c j + ck -1
Table 2. Output of the folding algorithm for the dependencies stream shown in Table 1 69], as they can be expressed as functions of the canonical induction variables. Detecting SCEVs is important for two reasons: First of all, the chains of dependencies associated to their computation shall be ignored, as otherwise it greatly and unnecessarily constrains possible code transformations; Second, it allows to detect an important class of memory access patterns, the strided accesses, that are used to evaluate the profitability of some loop transformations.

If we detect during folding of an instruction that all the Label values of a polyhedron can be expressed using a SCEV, that is, if the folding succeeds in constructing an affine function to express them, the corresponding points in the folded polyhedra of IIVs as well as the associated dependencies are removed from the DDG. This happens for example for instructions I5 and I8. For instruction I5, the folding algorithm finds a SCEV with value a(cj, ck) = 0cj + 1ck + 1.

Over-approximations Even in programs where the hot region is affine such as in PolyBench [START_REF] Pouchet | Polybench: The polyhedral benchmark suite (version 4.2)[END_REF], profiling the entire benchmark reveals a large amount of non-regular parts. Keeping an accurate representation of non-regular parts causes scalability issues both in the profiling part and in the polyhedral back-end of poly-prof. Our idea to handle such non-regular parts in a scalable fashion is through overapproximation. This over-approximation, detailed in [START_REF] Gruber | Building of a Polyhedral Representation from an Instrumented Execution: Making Dynamic Analyses of non-Affine Programs Scalable[END_REF], concerns labels for non-affine dependencies or non-affine SCEVs and polyhedra for dependencies and instructions with iteration domains with holes.

DDG Polyhedral Feedback

An essential motivation for folding DDGs into polyhedral structures is to enable the use of advanced polyhedral compilation systems, which are capable of finding a schedule that maximizes parallelism, finds tiling opportunities, etc. [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF].

Polyhedral compilation of folded-DDGs Typically, a polyhedral compiler is applied to small numerical kernels made of a handful of statements [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF][START_REF] Girbal | Semi-Automatic Composition of Loop Transformations for Deep Parallelism and Memory Hierarchies[END_REF][START_REF] Mehta | Improving Compiler Scalability: Optimizing Large Programs at Small Price[END_REF][START_REF] Pradelle | Scalable hierarchical polyhedral compilation[END_REF]. Polyhedral schedulers suffer scalability challenge for larger programs [START_REF] Mehta | Improving Compiler Scalability: Optimizing Large Programs at Small Price[END_REF]: their complexity typically grows exponentially with the number of statements in the input program. Our DDG folding and over-approximation techniques allow going from programs with thousands of statements (vastly exceeding the typical program scale these schedulers can handle) to only a few hundreds.

Numerous customizations for scalability of the polyhedral compiler have been implemented, ranging from constraining the space of legal schedules to accelerate the scheduling process to approximation of the code generation process to quickly output a decorated simplified AST describing the program structure after transformation. For example, the presence of large integer constants causes combinatorial blow up in the ILP solver used to solve the scheduling problem [START_REF] Papadimitriou | Combinatorial optimization: algorithms and complexity[END_REF]. We implemented a parameterization of iteration domains, to replace those constants by a parameter (an unknown, but constant integer value). That is, a domain The reader may refer to the available implementation in PoCC [START_REF] Pouchet | the PoCC polyhedral compiler collection[END_REF] for further details about the simplifications implemented, computation of profitability metrics is implemented in the PolyFeat module. Final output The main visual support used for reporting aggregated feedback is the dynamic schedule-tree described in Sec. 4, along with a simplified AST that shows the code structure after the application of the suggested structured transformation. This AST embeds various metrics on instruction count, loop properties (parallelism and tilability) and the list of statements surrounded by every loop. This lets the user visualize the potential effort in manually writing the code structure corresponding to the application of the suggested transformation.

{[i] : 0 <= i < 1024} is replaced by [n] → {[i] : 0 <= i < n ∧ n ≥ 1024}
To better visualize the various program regions, metrics can weight each tree-node and be rendered with flame-graphs [START_REF] Gregg | The flame graph[END_REF]. An example discussed in Sec. 7 is reported in Fig. 7. Here, the "surface" of a region in the flame graph is proportional to its estimated computation weight. Such flame graphs are SVG files, which can be clicked to obtain detailed information on each region/box. In the example of Fig. 7, loop/call nodes are marked with the label loop/call, respectively. Colors and gradients are used, that is, non-interesting regions can be grayed out. In Fig. 7, grayed regions are non-affine and blacklisted (initialization and extensive calls to libc) ones. Node width is used to highlight hotness of regions. In our example, functions adjustweight and layerforward take substantially more space than other regions.

Numerous additional extensions have been developed, to provide useful feedback to the user: various metrics about the size and structure of polyhedral regions in the DDG, before (original program) and after (transformed program) applying an optimizing schedule, the possible code structure including the precise fusion/distribution scheme, a simplified sequence of core loop transformations that should be applied to achieve the target transformation, etc. Given the extensive textual length of the feedback we provide, an example is shown only in the supplementary document.

Case Studies

This section aims to present case studies that illustrate both the type of feedback provided by poly-prof as well as the different transformations it can suggest. Current QEMU [START_REF] Bellard | QEMU, a Fast and Portable Dynamic Translator[END_REF], on which poly-prof is based, cannot handle newer AVX instructions. Feedback was obtained from binaries compiled Case study I This case study illustrates a simple (structured) feedback: pinpointing that dependencies live within the first quartant [START_REF] Utpal | Dependence Analysis for Supercomputing[END_REF] allows exposing data locality and fine-grain parallelism through a loop interchange. For this study, we selected backprop, a supervised learning method used to train artificial neural networks extracted from the Rodinia benchmark suite [START_REF] Che | A characterization of the Rodinia benchmark suite with comparison to contemporary CMP workloads[END_REF]. In addition to what was outlined in the previous case study, poly-prof exploits its ability to associate an iteration vector with each memory access, and to match pointer values with scalar evolution. Since the scalar evolution expressions for a memory access describe the addresses it reads/writes, they can be used to detect strided accesses along any dimension, which in turn can be used to detect vectorization potential. The flame-graph provided by poly-prof is reported in Fig. 7 where the regions of interest are: 1. the first call (of two) to bpnn_layerfoward with a constant value of n2 = 16 ; 2. the last call (of two) of bpnn_adjust_weights with a constant value of ndelta = 16. Both functions contain a 2Dnested-loop, called L layer and L adjust respectively, for which poly-prof suggests an interchange and SIMDization. For example, see L layer on Tab. 3 for which the two loop dimensions are reported to be on lines 253 (outer) and 254 (inner) of file backprop.c: The loop nest is fully permutable, that is, interchange is possible; only the outermost loop is parallel; and there are more stride 0/1 accesses along the outermost dimension (100%) than along the innermost (67%).

To enable the suggested transformations, one needs to specialize the two interesting function calls and array expand the scalar sum. The other calls to bpnn_layerfoward and bpnn_adjust_weights are left unchanged, since they (a) take up much less of the overall program runtime (b) have different arguments and do not profit from the transformation. Applying the transformation improved, in our case, bpnn_layerfoward from 0.5 GFlop/s to 2.8 GFlop/s and bpnn_adjust_weights from 0.3 GFlop/s to 5.1 GFlop/s.

Case study II This case study illustrates an advanced feedback: providing exact dependence vector "directions" allows exposing data locality and coarse-grain parallelism through loop skewing and tiling. For this study, we selected GemsFDTD, a finite difference time domain method from the SPEC CPU2006 benchmark suite [START_REF] Henning | SPEC CPU2006 Benchmark Descriptions[END_REF] written in Fortran90. Analyzing Fortran code is not a problem for poly-prof as it works at the binary level. However, the compiler we used (gfortran-8.1.1) messes up the debug information, making it necessary for the user to shift the line numbers for the provided code references by hand. This case study fully exploits poly-prof's ability to model (and compress in a polyhedral form) the data dependencies, instead of simply checking their existence/absence. This knowledge about the structure of dependencies allows poly-prof to check for tiling opportunities. First, polyprof detects that four functions from the benchmark are fat (i.e. they execute a large amount of the program's total number of dynamic instructions): updateH_homo, updateE_homo, UPML_updateH, and UPML_updateE; Inside the first two of those functions are the five hottest loop nests, so we focus on them. As reported in Fig. 4, poly-prof annotates all five loops as fully parallel and tilable. So to obtain a speedup we tile each loop along all dimensions with a tile size of 32 and mark the outermost loop parallel with an OMP PARALLEL DO directive. We recall that tiled code can always be also coarsegrain parallelized using wavefront parallelism, as exploited by the Pluto polyhedral scheduler [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF]. Tiling and parallelizing the loops increased performance in updateE_homo from 1.3 GFlop/s to 2.7 GFlop/s and updateH_homo from 1.3 GFlop/s to 3.7 GFlop/s

Experiments

The goal of this section is to demonstrate that poly-prof can be systematically applied on a full benchmark suite, and find potential for optimization. Note that the output of poly-prof for each benchmark is extensive: flame graph, statistics on each sub-region, potential structured transformation(s), simplified annotated AST after the application of the transformation, complete AST, etc. Consequently, we only illustrate here the application of poly-prof on Rodinia by using aggregate metrics. These metrics are not meant to be used as is by the end user: instead, the user is expected to work on one benchmark at a time, and navigate the feedback we provide.

Experiment I In addition to the case studies, we evaluate below our tool-chain using the latest revision (3.1) of the Rodinia benchmark suite [START_REF] Shuai | Rodinia: A benchmark suite for heterogeneous computing[END_REF][START_REF] Che | A characterization of the Rodinia benchmark suite with comparison to contemporary CMP workloads[END_REF]. As poly-prof does not support multithreaded applications yet, each benchmark is manually modified to run in a single thread, and compiled using the same compiler and flags as the case studies (GCC 8.1.1, -g -O2 -msse3).

The dynamic analysis is based on the QEMU-plugin instrumentation interface [START_REF] Guillon | Program Instrumentation with QEMU[END_REF]. We have extended the interface itself to support multiple interacting C/C++ plugins. Plugins primarily work at the level of the generic QEMU compiler instructions, making them CPU architecture agnostic. The choice of using QEMU for instrumenting the code is orthogonal to our contribution. There are some other candidates (e.g. Valgrind [START_REF] Nethercote | Valgrind: A Framework for Heavyweight Dynamic Binary Instrumentation[END_REF] or pin [START_REF] Luk | Pin: building customized program analysis tools with dynamic instrumentation[END_REF]) for implementing a dynamic binary analysis tool, each having their own advantages and disadvantages. Since we use a shadow memory to track data dependencies, dynamic analysis obviously does not come for free. As an example, the total CPU time (summing for all cores) on our server required by the first three stages of poly-prof to analyze the full Rodinia benchmark suite is 3h 6' (the full execution including libc was instrumented).

We have built our polyhedral feedback pass described in Sec. 6 in the PoCC compiler [START_REF] Pouchet | the PoCC polyhedral compiler collection[END_REF] that implements the PluTo scheduler [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF], along with a new polyhedral DDG analysis and optimization in PoCC/PolyFeat [START_REF] Pouchet | the PoCC polyhedral compiler collection[END_REF]. Extensions for scalability have been implemented in the schedulers and code generator, and in PolyFeat.

Experiment II To show problems static approaches encounter with the Rodinia suite we also ran the static LLVMbased [START_REF] Lattner | LLVM: A Compilation Framework for Lifelong Program Analysis & Transformation[END_REF] polyhedral compiler Polly [START_REF] Grosser | Polly -Performing polyhedral optimizations on a low-level intermediate representation[END_REF] over the entire suite. We used Polly version 7.0.1 and the flags -O3 -ffastmath -polly-process-unprofitable. Kernels that span multiple functions were inlined to allow Polly to see the same code region as poly-prof. Calls to functions from libc or the OpenMP runtime were not inlined. Where such calls are present this usually results in Polly being unable to analyze the kernel, though it can handle calls to simple functions such as exp or sqrt. Polly was unable to build a polyhedral model of the whole region of interest for any of the 19 benchmarks. It was able to model some smaller subregions, 1D or 2D loop nests, in most benchmarks, but in nearly all cases its own profitability metric decided not to optimize them. Two notable exceptions are the heartwall and lud benchmarks. In heartwall Polly was able to model a sequence of nine 2D loop nests which accounts for roughly two thirds of the code in the body of the kernel. For lud it managed to model the whole inner 3D loop nest of the kernel, but not the outermost loop. In all benchmarks the inability to model the outermost loops blocks Polly from exploiting the thread level parallelism inherent in the Rodinia benchmark suite. for heartwall, hotspot, and lud is the consequence of not supporting lattices at folding time: These programs contain hand linearized nested loops whose bounds use modulo expressions and so are not recognized as fully affine. Note that, even when parts of a benchmark are not affine, we can still find affine over-approximations for those regions, and potentially find transformations for the program as a whole.

Summary statistics.

Based on the statistics provided by poly-prof, the biggest region for which the optimizer suggests a transformation has been selected by hand. The code reference is reported in the column Region. We considered a region to be interprocedural (column interproc.) if inlining was required to perform the transformation or if it contained a call to libc or the OpenMP runtime. Column %ops reports the percentage of dynamic operations of the program executed inside the region, while %Mops and %FPops reports the percentage of memory (resp. floating point) operations of the region itself. Note that the sum of %Mops and %FPops can be greater than 100% since on x86 a single instruction can both load and store to/from memory and perform an operation.

The column Reasons why Polly failed lists the reasons why LLVM Polly was unable to model the whole region as an affine program. They are coded as: R. unhandled function call; C. complex CFG (break/return); B. non-affine loop bound or non-affine conditional statements; F. non-affine access function (includes pointer indirection); A. unhandled possible pointer aliasing; P. base pointer not loop invariant.

The next group of metrics shows what can be achieved via semantics-preserving structured transformations. skew displays whether skewing is used in the proposed transformation, we tend to avoid skewing unless it really provides improvements in parallelism and tilability. %||ops gives the percentage of dynamic operations that can be parallelized using OpenMP parallel pragmas. If a non-inner loop dimension is detected as parallel, then all its operations are considered to be parallelizable. As a loop has at least two iterations, at least two parallel blocks can be exposed when a loop is reported parallel. Similarly, %simdops reports the percentage of operations that occur in parallel innermost loops.

The %reuse/%Preuse metrics report space locality that is available in the program: %reuse is the percentage of load/stores that are stride-0 or stride-1 in the existing innermost loops in the program, while %Preuse reports the maximal percentage of load/store operations that can be made stride-0 or stride-1 via a sequence of loop permutations.

We report the maximal loop depth of the region in the source code (ld-src) and in the binary code (ld-bin). This shows whether the compiler performed any transformation that modifies the loop depth (e.g., full loop unrolling for cfd). Next the maximal tiling depth (TileD) is reported, along with %Tilops, the percentage of operations that can be tiled.

As soon as a region can be tiled, coarse-grain (wavefront) parallelism is possible, and data reuse could be improved. poly-prof does not currently provide feedback on temporal locality potential, but as illustrated in the backprop case study %reuse allows to evaluate spatial locality improvements through tiling/interchange.

Finally, metrics C/Comp./fusion outline the complexity of the loop fusion/distribution structure that originates from the structured transformation proposed, and is an indication of the difficulty to manually implement a transformed code. Any outermost loop with more than 5% of the total region operation counts as one "component". For example, if the region is made of two consecutive loop nests executing each half of the operations, then 2 components will be reported. C reports the number of components in the binary code; Comp. the number of components after applying the proposed structured transformation, using the fusion heuristic reported in fusion (M for maximal loop fusion, and S for smartfuse, a somewhat balanced fusion/distribution strategy).

Note that streamcluster exhausted memory at the scheduling stage, and therefore no result is displayed.

Related Work

This section describes previous work related to profiling, performance debugging, and trace compression.

Performance feedback tools Using profiling feedback is a widely used technique both for performance debugging and for guiding optimization [START_REF] Michael | Overcoming the Challenges to Feedbackdirected Optimization (Keynote Talk)[END_REF] (FDO). Tools such as Perf [START_REF] Carvalho De Melo | The new linux'perf'tools[END_REF] or Intel VTune [START_REF] Reinders | VTune performance analyzer essentials[END_REF] gather hardware counter values and, using sampling, provide statistics (such as instruction throughput or cache miss rate) at the granularity of machine instruction. Various interfaces such as HPCToolKit [START_REF] Adhianto | HPCTOOLKIT: tools for performance analysis of optimized parallel programs[END_REF] allow navigating through the gathered statistics thanks to a besteffort mapping from binary to source code. Compiler based profiling tools, such as gprof [START_REF] Graham | Gprof: A call graph execution profiler[END_REF], make it possible to provide complementary statistics (such as branch probabilities). While these allow pin-pointing important parts of the code where the programmer or compiler should focus his attention on, the metrics they report are very low level and consequently they leave the role of finding optimizing code transformations to the programmer or compiler heuristics [START_REF] Sopeju | AutoSCOPE: Automatic Suggestions for Code Optimizations Using PerfExpert[END_REF]. Recent work such as in MAQAO [START_REF] Andres S Charif-Rubial | CQA: A code quality analyzer tool at binary level[END_REF][START_REF] Djoudi | Maqao: Modular assembler quality analyzer and optimizer for itanium 2[END_REF], AutoSCOPE [START_REF] Sopeju | AutoSCOPE: Automatic Suggestions for Code Optimizations Using PerfExpert[END_REF] or MIAMI [START_REF] Marin | MIAMI: A framework for application performance diagnosis[END_REF] combines these metrics with static analysis of binary/source code. Intel IACA [START_REF]Intel Architecture Code Analyzer -User's Guide[END_REF], a purely static analysis tool, uses a precise machine model, which the vendor does not publish, to pin-point resource bottlenecks and predict performance of snippets of binary code.

Dynamic data-flow analyses have been presented to provide useful performance feedback. The detection of parallelism (such as vectorization) along canonical directions has particularly been investigated [START_REF] Ao | ParaInsight: An Assistant for Quantitatively Analyzing Multi-granularity Parallel Region[END_REF][START_REF] Butt | Runtime analysis of application binaries for function level parallelism potential using QEMU[END_REF][START_REF] Faxén | Embla-data dependence profiling for parallel programming[END_REF][START_REF] Ketterlin | Profiling Data-Dependence to Assist Parallelization: Framework, Scope, and Optimization[END_REF][START_REF] Kim | Prospector: A dynamic data-dependence profiler to help parallel programming[END_REF][START_REF] Li | DiscoPoP: A profiling tool to identify parallelization opportunities[END_REF][START_REF] Tournavitis | Semi-automatic extraction and exploitation of hierarchical pipeline parallelism using profiling information[END_REF][START_REF] Vandierendonck | The Paralax infrastructure: automatic parallelization with a helping hand[END_REF][START_REF] Wang | Integrating profile-driven parallelism detection and machine-learning-based mapping[END_REF], as it requires only relatively localized information. Another use case is the evaluation of effective reuse [START_REF] Berg | Fast data-locality profiling of native execution[END_REF][START_REF] Beyls | Discovery of Locality-Improving Refactorings by Reuse Path Analysis[END_REF][START_REF] Liu | Pinpointing data locality problems using data-centric analysis[END_REF][START_REF] Marin | MIAMI: A framework for application performance diagnosis[END_REF] with the objective of pinpointing data-locality problems.

The APOLLO [START_REF] Jimborean | VMAD: An Advanced Dynamic Program Analysis and Instrumentation Framework[END_REF][START_REF] Manuel | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF] framework uses an interesting technique that corresponds to detecting dynamic regularities (affine expressions) of memory accesses in nested loops, and use this information to perform speculative loop transformation, which differs in scope and challenges from the profiling feedback tool developed in this paper. Polyhedral compilation Integer linear algebra is a natural formalism for representing the computation space of a loop nest. The polyhedral framework [START_REF] Feautrier | Polyhedron model[END_REF] leverages, among others, operators on polyhedrons (e.g. union, intersection, projection), enumeration (for code generation [START_REF] Bastoul | Generating loops for scanning polyhedra: Cloog users guide[END_REF]), and parametric integer linear programming [START_REF] Feautrier | Parametric integer programming[END_REF] (for dependence analysis [START_REF] Collard | Fuzzy array dataflow analysis[END_REF]). Historically, it has been designed to work on restricted programming languages, and was used as a framework to perform source-to-source transformations. More recently, efforts have been made to integrate the technology in mainstream compilers (e.g. Graphite [START_REF] Pop | GRAPHITE: Polyhedral analyses and optimizations for GCC[END_REF][START_REF] Trifunovic | GRAPHITE Two Years After: First Lessons Learned From Real-World Polyhedral Compilation[END_REF] for GCC [START_REF] Griffith | GCC: the complete reference[END_REF] and Polly [START_REF] Grosser | Polly -Performing polyhedral optimizations on a low-level intermediate representation[END_REF] for LLVM [START_REF] Lattner | LLVM: A Compilation Framework for Lifelong Program Analysis & Transformation[END_REF]). The set of loop transformations (known as affine transformations) that the polyhedral model can perform is wide and covers most of the important ones for exposing locality and parallelism to improve performance [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF]. Tools such as PoCC [START_REF] Pouchet | the PoCC polyhedral compiler collection[END_REF] provide a convenient interface to the numerous existing state-of-the-art scheduling heuristics and polyhedral libraries such as ISL [START_REF] Verdoolaege | isl: An Integer Set Library for the Polyhedral Model[END_REF]. Dynamic dependence graph Shadow memory [START_REF] Zhao | Umbra: Efficient and scalable memory shadowing[END_REF] records a piece of information for each storage location used in a program. For dependency tracking this is usually the last statement or dynamic instruction that modified that location. But shadow memories are also a core component of memory error debugging tools [START_REF] Nethercote | Valgrind: A Framework for Heavyweight Dynamic Binary Instrumentation[END_REF][START_REF] Serebryany | AddressSanitizer: A Fast Address Sanity Checker[END_REF].

Like poly-prof, Redux [START_REF] Nethercote | Redux: A dynamic dataflow tracer[END_REF] builds a data-flow graph from binary level programs. It has, however, no notion of loops or calling contexts and does not try to compress the produced graph and is consequently only able to handle very small programs. A few techniques exist to address the high overhead of monitoring data dependencies. The most common is the use of static analysis so as to remove redundant instrumentation [START_REF] Li | Fast data-dependence profiling by skipping repeatedly executed memory operations[END_REF], or to avoid monitoring must-dependencies [START_REF] Ketterlin | Profiling Data-Dependence to Assist Parallelization: Framework, Scope, and Optimization[END_REF]. Another approach that leads to over-approximation, consists of reducing the size of the shadow memory through the use of signature based addressing [START_REF] Li | An efficient datadependence profiler for sequential and parallel programs[END_REF]. The most sophisticated technique is the one developed in SD 3 [START_REF] Kim | SD3: An Efficient Dynamic Data-Dependence Profiling Mechanism[END_REF] that amounts to detecting strided memory accesses to compress the shadowing. A last technique that could be used for our purposes is the parallelization of shadowing as done in [START_REF] Li | Fast data-dependence profiling by skipping repeatedly executed memory operations[END_REF][START_REF] Moseley | Shadow profiling: Hiding instrumentation costs with parallelism[END_REF].

Calling context tree Using calling context trees to disambiguate instructions in different calling contexts is an idea from Ammons, Ball and Larus [START_REF] Ammons | Exploiting Hardware Performance Counters with Flow and Context Sensitive Profiling[END_REF]. But in the presence of recursive functions the size of their CCTs grows proportionally with the depth of the recursion, leading to an unreasonably high memory overhead. Loop-call context trees [START_REF] Sato | On-thefly Detection of Precise Loop Nests Across Procedures on a Dynamic Binary Translation System[END_REF] simply encode the calling context of intraprocedural loops and suffer from the same size problems in recursive programs.

Conclusion and Future Work

In this paper we introduced poly-prof, a profiling-based polyhedral optimization feedback tool. poly-prof tackles the problem of dynamic data dependence profiling and polyhedral optimization from binary programs, while addressing the associated scalability challenges of handling traces with billions of operations, as we demonstrated. Our tool handles non-regular programs with recursive function calls. Numerous technical contributions as presented in this paper were required to enable structured transformation feedback on binary programs, a significant step in complexity compared to prior approaches typically limited to unstructured transformations. We implemented numerous feedback reporting schemes for the user in poly-prof: flame graphs, statistics on each sub-region, proposed potential structured transformation, simplified annotated AST after the application of the transformation, complete AST, etc.

While this represents an important step towards being able to provide polyhedral feedback on full-scale applications there are issues that still need to be addressed. This includes the development of under-approximation schemes in the DDG, as well as overall scalability enhancements, like using approximate (non-optimal) polyhedral scheduling strategies.

Lastly, the process of mapping our results back from the machine code back to the source code is in itself a research topic, and currently we do not provide much more than what objdump does. Our ongoing efforts in this direction leverage polyhedral program equivalence techniques [START_REF] Bao | PolyCheck: Dynamic Verification of Iteration Space Transformations on Affine Programs[END_REF][START_REF] Schordan | Verification of polyhedral optimizations with constant loop bounds in finite state space computations[END_REF][START_REF] Verdoolaege | Equivalence checking of static affine programs using widening to handle recurrences[END_REF].

Figure 1 .

 1 Figure 1. Overview over the poly-prof framework.

Figure 2 .

 2 Figure 2. (a/b) Example CFG and associated loop-nestingtree; (c/d) Example CG and associated recursive-componentsset.

Figure 3 .

 3 Figure 3. (a/f) Code skeletons; Purple is used to represent CFG, while orange is used to represent CG; (b/g) Call-graph and recursive-component-set; (c) Loop-nesting-forests; (h) Calling-context-tree with call-sites information; (d/i) Example trace, associated loop events, and dynamic interprocedural interation vectors. R 2 stands for two consecutive returns; (e/j) Dynamic schedule tree; (k) Folded domains.

Figure 5 .

 5 Figure 5. Dynamic schedule tree

  prior to scheduling. We control the number of parameters introduced by reusing the same parameter for a range of values, that is, if the value x ∈ [1024s, 1024 + s] with s ∈ Z (we typically set s = 20), then we replace x by n + (x -1024).

Figure 7 .

 7 Figure 7. Annotated flame-graph for backprop (see Sec 7). with GCC 8.1.1 and the flags -g -O2 -msse3. Speedup measurements are using Intel icc and ifort compilers (version 18.0.3, flags -Ofast -march=native -mtune=native) on a machine with a Xeon Ivy Bridge CPU with two 6 core CPUs (24 hyperthreads), each running at 2.1Ghz (GFlop/s are averaged over 50 runs).

  Stack of currently live loops. From outermost to innermost.

	Emitted events:
	• E(L, H): entry into loop L; H is the header.
	• I(L, H): iteration of loop L; H is the header.
	• X(L, B): exit of loop L, jumping to basic-block B
	• N(B): local jump to basic-block B.
	1 if event . type is local -jump :
	2	B := event . dstBB
	3	while L := inLoops . peek () and L . isCFG and B not in L:
	4	L . visiting := false ; inLoops . pop (); emit X(L ,B)
	5	if B is loop -header :
	6	H := B ; L := loop ( H )
	7	if L . visiting = false :
	8	L . visiting := true ; inLoops . push ( L ); emit E(L ,H)
	9	else : emit I (L , H )
	10	emit N ( B )
	11 else :
	12	. . . # Alg. 2

  from a function header of recursivecomponent L and iteration of the corresponding loop. B is the current basicblock after the call/return. • X R (L, B): return from a function header of recursive-component L and exit from that loop.

1 if event . type is call : 2 F := event . callee ; B := event . dstBB 3 L := F . recursive_component 4 if F is recursive -component entry and L . entry == undef : 5 L . entry := F ; inLoops . push ( L ); emit E C (L ,B) 6 elif F is recursive -component header : 7 while L ' := inLoops . peek () and L ' in L : 8 L ' . visiting := false ; inLoops . pop (); 9 emit X ( L ' ,B ) 10 L . stackcount ++; emit I C (L , B ) 11 else : emit C (F , B ) 12 if event . type is return : 13 F := event . callee ; B := event . dstBB 14 while L ' := inLoops . peek () and L ' . isCFG and L' in F 15 L ' . visiting := false ; inLoops . pop (); emit X(L' ,B) 16 L := F . recursive_component 17 if F is recursive -component entry and 18 L . stackcount == 0 and L . entry == F : 19 L . entry := undef ; emit X R (L , B ) 20 elif F is recursive -component header : 21 L := F . recursive_component ; L . stackcount --; 22 emit I R (L , B ) 23 else : 24 . . . # Alg. 1

  54,

	Id	Polyhedron (cj,ck)	Label expression f(cj,ck)

Table 3 .

 3 backprop case study. Reported lines (e.g. 52, 253, . . . ) are from debug information. Suggested interchange is represented using a permutation of code lines. Statistics/properties per loop dimension as follows: (outer, inner).

	Fat regions			% ops	interchange+SIMD	parallel permutable % stride 0/1 speedup
	backprop_kernel.c:52 (L layer )	14% backprop.c:(254,253) (yes,no)	(yes,yes)	(100%, 50%)	5.3 x
	backprop_kernel.c:57 (L adjust ) 46% backprop.c:(322,320) (yes,yes)	(yes,yes)	(100%, 50%)	7.8 x
	Fat regions	op	tiling		speedup
	update.F90:106 20%	update.F90:{106,107,121}	2.6 x
	update.F90:240 18%	update.F90:{240,241,244}	1.9 x

Table 4 .

 4 GemsFDTD case study. Reported lines are shifted debug info.

  Table 5 presents summary statistics about the Rodinia 3.1 (CPU only) benchmarking suite, that we computed/aggregated by processing the feedback from poly-prof and Polly on each benchmark.Column %Aff reports the percentage of dynamic operations that are part of a fully affine region without overapproximation. The low proportion of affine code reported

	Bench.	#V	#E	%Aff	Region	%ops	%Mops	%FPops	interproc.	Reasons why	Polly	failed	skew	% | |ops	%simdops	%reuse	%Preuse	ld-src	ld-bin	TileD	%Tilops	C	Comp.	fusion
	backprop	10 M	10 M	85%	facetrain.c:25	67%	30%	76%	Y	A		N	100%	100%	50%	100%	2D	2D	2D	100%	6	4	S
	bfs	2 M	1 M	21%	bfs.cpp:137	55%	66%	18%	N	BF		N	100%	100%	1%	1%	3D	3D	2D	100%	1	1	M
	b+tree	23 M	24 M	49%	main.c:2345	26%	99%	0%	N	BF		N	100%	100%	44%	44%	3D	3D	3D	100%	15	4	S
	cfd	251 M	372 M	98%	*3d_cpu.cpp:480	98%	42%	99%	Y	F		N	100%	61%	18%	42%	5D	4D	3D	100%	1	3	S
	heartwall	4 G	8 G	1%	main.c:536	99%	38%	56%	Y	RCBF	N	100%	100%	0%	0%	7D	6D	5D	100%	1	3	S
	hotspot	11 M	16 M	0%	*_openmp.cpp:318	81%	35%	89%	Y	B		Y	100%	100%	3%	3%	4D	4D	2D	100%	1	1	S
	hotspot3D	210 M	256 M	99%	3D.c:261	49%	28%	81%	N	BF		N	100%	99%	11%	11%	4D	4D	3D	100%	1	1	M
	kmeans	513 M	647 M	97%	*_clustering.c:160	97%	56%	98%	Y	RFA		N	100%	100%	46%	53%	4D	4D	4D	100%	1	3	S
	lavaMD	879 M	1 G	0%	kernel_cpu.c:123	99%	69%	92%	N	BF		N	100%	100%	0%	0%	4D	4D	3D	100%	1	2	S
	leukocyte	4 G	9 G	39%	detect_main.c:51	37%	64%	64%	Y	RCBFAP	N	100%	100%	63%	63%	4D	4D	3D	100%	11	5	S
	lud	42 M	71 M	4%	lud.c:121	99%	44%	70%	Y	BF		N	99%	98%	0%	1%	5D	5D	3D	99%	3	3	S
	myocyte	1 M	866 K	89%	main.c:283	99%	80%	80%	Y	CBA	N	100%	99%	47%	47%	4D	3D	1D	99%	1	3	S
	nn	1 M	2 M	1%	nn_openmp.c:119	31%	42%	71%	Y	RF		N	100%	0%	0%	0%	1D	1D	1D	100%	1	1	M
	nw	80 M	93 M	99%	needle.cpp:308	79%	73%	27%	Y	RF		Y	100%	100%	77%	77%	4D	4D	2D	100%	2	2	S
	particlefilter	628 M	678 M	27%	*_seq.c:593	99%	5%	14%	N	CF		N	99%	100%	55%	55%	3D	3D	2D	100%	22	2	S
	pathfinder	62 M	48 M	67%	pathfinder.cpp:99	31%	83%	16%	N	BP		Y	100%	0%	0%	40%	2D	2D	2D	100%	1	1	M
	srad_v1	1 G	2 G	99%	main.c:241	99%	31%	93%	Y	RF		N	99%	100%	18%	18%	3D	3D	2D	100%	1	1	S
	srad_v2	600 M	864 M	98%	srad.cpp:114	96%	31%	92%	Y	RF		N	100%	100%	14%	14%	3D	3D	2D	100%	1	1	S
	streamcluster	779 M	1 G	97%	*_omp.cpp:1269	99%	6%	13%	Y	RCBFAP	-	-	-	-	-	6D	6D	-	-	52	-	-

Table 5 .

 5 Summary statistics computed from poly-prof's feedback on the Rodinia benchmark suite.

the Rodinia benchmark suite, for example, does not use any recursion

Acknowledgements

This work was supported in part by the U.S. National Science Foundation awards CCF-1645514, CCF-1731612 and CCF-1750399, and by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement d'avenir.


Raphael Abele
email: raphael.abele@st.com

Daniele Fronte
email: daniele.fronte@st.com

Pierre-Yvan Liardet
email: pierre-yvan.liardet@st.com

Jean-Marc Boi
email: jean-marc.boi@lis-lab.fr

Jean-Luc Damoiseaux
email: jean-luc.damoiseaux@lis-lab.fr

Djamal Merad
email: djamal.merad@lis-lab.fr

Autofocus in infrared microscopy

Keywords: Infrared imaging, Microscopy, Focusing, Image analysis, Image quality, Image decomposition, Wavelet transforms, Statistics

Autofocus (AF) is a widely investigated subject in the fields of natural scene images, industrial assembly and biologic microscopy. This paper proposes a new effective AF method for infrared (IR) microscopy in the context of the Integrated Circuit industry (IC). The proposed method operates in the wavelet domain using a custom orthogonal wavelet for the 2D Discrete Wavelet Transform (DWT). The quality criterion of our AF algorithm relies on the standard deviance of the DWT coefficients, computed per subband and per level. Tested on several optical magnifying lenses, our method is robust time-efficient, and usable on-the-fly in the IC location system.

I. INTRODUCTION

Integrated Circuits (ICs 1 ) are electronic components whose applications seem unlimited: they are used in many fields such as wearable technologies and IoT (Internet of objects). In order to protect the privacy of users, some ICs need to be secured, and their security needs to be validated. One way to characterize a secure IC is to study its behavior following a physical disruption. Such disruption may be obtained by a laser shot targeted on the internal structures of the IC. The accuracy of the shot relies on the precision of the laser power calibration and 3D positioning inside the IC. In our study, the Object Of Interest (OOI) is considered as made of three layers of different materials: (1) the silicon, (2) a conductor material forming the conductive tracks, and (3) an electrical insulator, as illustrated in Fig. 1.

The internal OOI structures are visible on the surface of the conductive tracks, which is the Surface Of Interest (SOI). The Project funded by the Association Nationale de la Recherche et de la Technologie (ANRT). 1 A glossary is present at the end of this paper.

Fig. 1: Concept mapping of three materials composing our IC.

3D positioning of the laser on the SOI is done with an infrared (IR) camera that can see through the silicon layer. In our work, the optical system and the laser source are interdependent, so the optical adjustment implies laser position adjustment.

Usually, these adjustments are made manually by skilled men, which raises two problems: time loss and non-reproducibility of the characterization (because of potential imprecision). An autofocus (AF) system could partially solve this problem by automating one dimension of the 3D positioning process. Section (II) shortly introduces some works related to the AF. Then, section III presents our particular work environment. Section IV describes the AF method we propose, and section V details the overall experiment results.

II. AROUND THE AUTOFOCUS

The AF mechanism is a deterministic algorithm that determines the lens position for which the system is "well focused" on the scene/object of interest (see Fig. 2), that is, the position for which the image of the scene or object is sharpest.

Two AF approaches exist: active and passive. The active approach depends on an additional system measuring the distance from the system to the scene/object. For example, ultrasonic sound waves and IR light reflection are two ways of measuring distances. Once the optics-scene distance is known, the correct lens position to obtain a focused image is estimated. Such an active system is not available, so we consider a passive approach relying on image analysis. This analysis determines the most relevant object/scene image, which we then use to find the correct focus. At this point, we are confronted with the issue of Image Quality Assessment (IQA). The Human Visual System (HVS) is the most reliable tool for IQA: the main difficulty is the interpretation of human subjectivity with algorithms following objective rules and criteria [START_REF] Wang | Why is image quality assessment so difficult?[END_REF]. Some subjective methods attempt to emulate HVS using the Mean Opinion Score, but they are time-consuming and generally non-deterministic; these methods are partly based on objective metrics pooled together to predict the quality score [START_REF] Streijl | Mean opinion score (MOS) revisited: methods and applications, limitations and alternatives[END_REF].

Objective IQA methods attempt to be correlated with the HVS. Depending on the amount of available information from the image source, methods are categorized as Full-Reference (FR) based if an original, non-distorted image is available [START_REF] Sheikh | Image information and visual quality[END_REF] or No-Reference (NR) based otherwise [START_REF] Kamble | No-reference image quality assessment algorithms: A survey[END_REF]. If just some of the original image features are known, the method is qualified as Reduced-Reference (RR) based. See reference [START_REF] Kim | Automatic focus control for assembly alignment in a lens module process[END_REF] as an example. RR-IQA methods are generally guided/optimized NR-IQA methods. In the AF context, both of the NR and RR-IQA methods are used, since the quality criterion is hypothetical. Images are ranked according to this criterion, and the best-ranked image is associated with the correct lens position.

Statistics of order 1, 2 and higher are widely used in signal processing since the 1960's. In image analysis, each statistic point out an image feature [START_REF] Kumar | Importance of statistical measures in digital image processing[END_REF]:

• Average Gradient (AG): reflects the contrast and the clarity of the image. It can be used to measure the spatial resolution of a test image, where a larger average gradient indicates better spatial resolution. • Skewness metric: a statistical measure of the direction and extent to which a data set deviates from a distribution. For a standard normal distribution, high skewness indicates asymmetry of the data. In this case, the data contains a greater amount of information. These statistics can be computed on the image data or on its transform, that is, in spatial, frequency or time-frequency domains.

Since AF algorithms rely on specific image analysis, they cannot be generalized. Several methods exist, each adapted to its context such as images of natural scenes [START_REF] Kim | Automatic focus control for assembly alignment in a lens module process[END_REF], [START_REF] Schlag | Implementation of automatic focusing algorithms for a computer vision system with camera control[END_REF]- [START_REF] Zhang | A new auto-focus measure based on medium frequency discrete cosine transform filtering and discrete cosine transform[END_REF], low-contrast images [START_REF] Xu | Robust automatic focus algorithm for low contrast images using a new contrast measure[END_REF] or microscopic images [START_REF] Mateos-Pérez | Comparative evaluation of autofocus algorithms for a real-time system for automatic detection of mycobacterium tuberculosis[END_REF]- [START_REF] Zhang | Evaluation of autofocus algorithm for automatic dectection of caenorhabditis elegans lipid droplets[END_REF]. In paper [START_REF] Zhang | Evaluation of autofocus algorithm for automatic dectection of caenorhabditis elegans lipid droplets[END_REF], the author studied several autofocus methods in a particular case (lipid droplets inside microscopic worms). The focus criterion is clarity-evaluation based. He tested the performance and accuracy of sixteen autofocus methods: histogram, intensity, statistic, derivative and transform-based. He found that the absolute Tenengrad algorithm had the best performance against accuracy for its specifics. Finally, each approach uses well justified metrics. In our proposal we define metrics adapted to IR microscopy for the IC industry, which implies some constraints (see Section III).

III. WORK ENVIRONMENT

A. Available material

Observation of the SOI under the silicon layer is possible thanks to its optical properties. Indeed, around the wavelength 1µm and higher, the silicon absorption coefficient is not significant [START_REF] Green | Optical properties of intrinsic silicon at 300 k[END_REF], [START_REF] Green | Self-consistent optical parameters of intrinsic silicon at 300 k including temperature coefficients[END_REF]. In other words, silicon is "transparent" to lights whose wavelengths are around 1µm and above (such as IR light).

Our motorized optical system provided by AlphaNov 2 is composed of:

• an uncooled IR camera • an optical microscope with 4 magnifying lenses (2.5x, 5x, 20x and 50x) allowing us to observe the OOI (around 5mm 2 ) and its internal structures (micro-metric scale) Three types of IR camera are available on the market: Long-Wave, Mid-Wave and Short-Wave (LWIR, MWIR and SWIR, respectively). On the one hand, MWIR and LWIR sensors detect thermal emissions from objects, and are efficient when these are warmer than their surroundings. On the other hand, SWIR cameras use the reflected light, much like the slightly shorter wavelengths of the visible spectrum (see Fig. 3). In our experiment we worked with the latter type of camera.

With the unique capabilities of our SWIR camera, we are able to see the SOI through the silicon layer. To propose an automated method for focusing on the SOI, our algorithm must take in account some specifics arising from the camera and from the SOI (see Section III-B).

B. Constraints and special features

Concerning the camera: The IR camera sensor is an InGaAs, sensitive to a large IR wavelength range (900nm to 1700nm); the embedded Wide Dynamic Range technology (WDR) dynamically detects the relevant wavelengths to take into account for best image quality. Despite its accuracy, this camera stays sensitive to external thermal interference. Since the sensor is uncooled (thermally not insulated), the thermal instability of the environment produces thermal noise, which affects the image acquisition (see Fig. 4).

Cooled IR camera

Uncooled IR camera Fig. 4: An example of images acquired with cooled and uncooled infrared cameras (source: flir.com).

Concerning the object of interest:

Two distinct surfaces are considered for the autofocus algorithm: the silicon and the SOI. As said in III-A, silicon is transparent to our IR light and should not be visible. However, in many cases, contaminants (e.g. dusts) are ingrained on the silicon, which reveals the surface to any light. This implies a possible conflict of interest, considering the sharpness criterion: the silicon surface with its contaminants may be sharper than the SOI. This point underlines that sharpness is not, in practice, the appropriate criterion to exclusively characterize an SOI. Fig. 5 shows to the notable states during the autofocus process.

We may compare the conductive layer of an IC to a city, and its internal structures to buildings. As shown in Fig. 5: Images acquired at several focal point positions (green points) according to the Z-axis. Component out of focus (left), focused on the silicon surface (middle) and on the conductive track surface (right). Fig. 6, everything is orthogonally disposed, in regard to the rectangular shape of the component. This criterion is a good candidate for our AF algorithm.

However, a difficulty arises from the SOI image texture. An IC is built by microscopic material deposition, which implies a textural grain visible with a high magnifying lens (20x and 50x). Considering the SOI images, this grain can be compared to a noise in the images: structures are clusters of points, more or less dense (illustrated in Fig. 7). In the worst case, detection of the structure boundaries turns out to be difficult and the images could need preprocessing to be well treated.

To sum up, we may state that: a. Sharpness is not a good criterion to make the difference between the silicon surface and the SOI. b. On the SOI, horizontal and vertical structures boundaries are visible. c. Because of textural grain and noise, topological information is not pertinent (e.g. line detection). d. The AF method has to be effective for each magnifying lens (2.5x, 5x, 20x and 50x). Since we know which features are expected in the targeted image (b.), our IQA can be qualified as RR-based. To avoid preprocessing such as image denoising (c.), we propose a transform-based analysis in the next section.

IV. AUTOFOCUS METHOD PROPOSAL

The criterion analyzed for an AF system has to be strong enough to prevent false results. In that sense, for a robust AF algorithm, we need a robust quality metric.

A. Analysis domain proposal

As previously described, image preprocessing should be performed before any image-based analysis (e.g. denoising). But it could break information continuity in our images (e.g. spatial information), and could be time-consuming. To preserve this continuity, our approach is transform-based, which may also save computation time. This transform is described below.

Wavelet transform:

Several image transforms exists. One of the most commonly used is wavelet decomposition, but has limitation due to its orthogonality (i.e. there is little orientation information). Others exist such as the steerable pyramid and curvelet transforms: the steerable pyramid transform overcome this orientation limitation and the curvelet transform is an extension of the wavelet concept, with less redundancy and more orientation information [START_REF] Simoncelli | The steerable pyramid: a flexible architecture for multi-scale derivative computation[END_REF]- [START_REF] Hill | Contrast sensitivity of the wavelet, dual tree complex wavelet, curvelet, and steerable pyramid transforms[END_REF]. However, according to the features described in section III-B, the main information needed for our analysis is the horizontal and vertical components of images. Therefore, we only consider wavelet decomposition for our IQA. Moreover, it may permit skipping some useless information related to material texture or even noise. For this purpose, the wavelet choice is important, considering that each wavelet matches specific signal information.

Wavelet choice:

Wavelets are used to decompose signal information, depending on their characteristics and topology. Here we need to find the best wavelet matching to the interesting information in our IR images, that is, the structures contours. Vu and Chandler [START_REF] Vu | A fast wavelet-based algorithm for global and local image sharpness estimation[END_REF] described a method to construct an image sharpness map based on 2D wavelet decomposition. The implementation of such a sharpness map could allow us to visually evaluate how image information is decomposed by a wavelet.

To that end, we took three photos of the SOI with 5x, 20x and 50x magnifying lenses. For each image, we built a sharpness map considering the following: two orthogonal and two bi-orthogonal wavelets, respectively:

1. Haar wavelet [START_REF] Haar | Zur theorie der orthogonalen funktionensysteme[END_REF], since the image can be interpreted as binary: reflexive versus non-reflexive materials, structures versus background (I). 2. A Custom Orthogonal Wavelet (COW) built from a father wavelet whose scaling coefficients are defined in (1) and and from which filter bank coefficients are computed in Tab. III. See its illustration in Tab. II. 3. Cohen-Daubechies-Fauveau (CDF 9/7) 9/7 wavelet [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF], having a great capacity to extract textures, is illustrated in Tab. IV. 4. An Optimized CDF 9/7 9/7 wavelet proposed in [START_REF] Zhang | A novel 9/7 wavelet filter banks for texture image coding[END_REF] and illustrated in Tab. V. Our goal is to find which wavelet match the most the relevant information in the image, that is, structures salience. Tab. VI shows obtained sharpness maps. We visually estimate that COW seems to match the interesting information, that is, relative to structures' salience. Contrary to CDF 9/7 9/7 wavelet, COW seems to avoids a great part of the image information relative to thermal noise and textural grain.

φ = 1 k 1 k 1 k 1 k (1)

B. Our IQA method

In [START_REF] Vu | A fast wavelet-based algorithm for global and local image sharpness estimation[END_REF], Vu and Chandler also proposed a wavelet-based algorithm for estimating both the global and local image sharpness (Fast Index SHarpness -FISH). Their IQA comes from the log-energy calculation of a three-level separable Discrete Wavelet Transform (DWT) using CDF 9/7 9/7 wavelet. The choice of this wavelet can be justified by its effectiveness to match textures, and thus to extract a large range of data. The energy evaluation brings information about the overall data values. The sharper an image is, the more its DWT contains high frequencies; and the more there are sparse frequencies, the higher the energy is. Here the sharpest image could be either a well focused image, or an unfocused one, because each contains sparse high frequencies.

We generalize a DWT analysis as follows: 1) 2D DWT of grayscale image using a given wavelet, on N levels. Let XY n denote the DWT subbands at each level n, where XY is either LH, HH or HL DWT subbands. 2) Given a statistic function F , F XYn is F computed on each subband XY at each decomposition level n. 3) For each decomposition level n, F n is the pounding sum of each F XYn : where α, β and γ are coefficients with α + β + γ = 1. 4) The total F tot of the DWT is the sum of each F n weighted as follows:

F n = α F LHn + β F HLn + γ F HHn (2)
F tot = N n=1 2 N -n F n (3)
Our algorithm follows these four steps, with the specifics listed below:

• According to the study in section IV-A, COW is used to match the relevant image information. • We consider the Standard Deviation (SD) (4), reflecting the image data variability or diversity, since we cannot rely on statistics based on spatial arrangement information, or on its quantity/quality estimation (see section III-B).

SD XYn = 1 C C c=1 (S XYn (c) -µ) 2 (4)
where µ is the mean of the considered subband, and c is the coefficient number of this subband from 1 to C. • Only the HL and LH subbands (5) are used: horizontal and vertical components are kept, leaving the diagonal information held by HH (see section III-B).

SD n = SD LH + SD HL 2 (5)
• The weighting to obtain SD tot is reversed compared to F tot described below:

SD tot = N n=1 2 n SD n (6)
For each decomposition level, some information are dropped due to the sub-sampling. Then, in the highest level, the root image information remains. With this IQA we are then able to rank images by their estimated quality. Our AF method is based on this criterion: the best ranked image in a video corresponds to the focused image, from which the focus position can be estimated (see next section).

V. EXPERIMENT AND DISCUSSION

Considering a linear movement given to the optical system toward the OOI, we take a video. The initial position is the farthest position from the OOI, the last is the nearest possible (the physical limits). We make such an acquisition for each magnifying lenses (2.5x, 5x, 20x and 50x), at 80 FPS. From a video, images are extracted so that each image correspond to a lens position, and are evaluated with our proposed quality metric. In this way, the highest-ranked image gives the best matching focus according to our quality criterion.

A. Results

We firstly compare our algorithm to FISH [START_REF] Vu | A fast wavelet-based algorithm for global and local image sharpness estimation[END_REF]. Fig. 8 presents the results.

The two approaches have different behaviors. With every magnifying lens, our algorithm effectively catch the SOI focus position (the star in Fig. 8); however, we observe unexpected behaviors around the focus point, which could correspond to optical distortions. The FISH algorithm catch the SOI focus position, but is largely more sensitive to distortions from the optics, the camera and other non-relevant information. This reflects the importance of the textural overcoming the structural information. This sensitivity bias the result, which is inoperable for 2.5x and 5x magnifying lenses. These results show that our method is adapted and robust to our specifics.

In order to better describe the image characteristics, we compute other statistics on the same scheme as our algorithm: entropy, Auto-Correlation (AC), skewness and kurtosis. See Fig. 9.

Considering the interpretation of each statistic applied to images, these measures reveal a complex image information. Dis-focused or not, images are saturated: in natural scenes, dis-focused images are considered as blurred. In this case the image data are smoothed and contains little information. In our study the textural part of the data is strong, and light distortions implied by the optical system add unexpected features to the images. Then, even using the DWT with our COW, spatial information is too noisy to be interpreted without preprocessing. The SD measure has a greater abstraction level, enough to remain a reliable quality criterion.

B. Runtime

As studied in [START_REF] Rioul | Wavelets and signal processing[END_REF] and [START_REF] Shukla | Efficient algorithms for discrete wavelet transform: with applications to denoising and fuzzy inference systems[END_REF], the computational complexity of DWT using filter banks implementation is O(n). Our algorithm threat the wavelets decomposition linearly per subband. Since only two subbands per level are considered, and for a three-level DWT, the global complexity is around:

O(n) + 2O( n 2 ) + 2O( n 4 ) + 2O( n 8 ) ≈ O(n) (7)
We implemented our architecture in Python language, on an 64-bit PC with an Intel Core i5-6200U at 2.30GHz and 8Gb of RAM. The camera resolution is 320x256 pixels. Additional packages are used for matrix operations (NumPy) and for DWT (PyWavelet) using filter banks.

For a video containing 1500 images, our implementation run in around 3.37 seconds. Then, an image quality is computed in 2.24 ms. In our application, this time is not disturbing: the focus research can be reduced if a previous AF is already done. This solution may considerably reduce video length and so the image stack.

C. Evaluation of silicon thickness

As said in Section I, two adjustments are necessary before a laser shot: the 3D positioning, and the laser power. The laser power setting depends on the silicon thickness on the SOI. A simple way to measure the silicon thickness consists in subtracting the positions of the SOI (1) and the silicon surface [START_REF] Streijl | Mean opinion score (MOS) revisited: methods and applications, limitations and alternatives[END_REF]. The SOI position is available with the proposed algorithm. Changing our AF criterion may permit to focus on (2) instead of (1). For example, in the case of contaminated silicon surface, the sharpness could be a good criterion.

VI. CONCLUSION

We proposed a robust passive autofocus algorithm. Our quality criterion allows us to efficiently focus on the conductive track surface inside an integrated circuit, independently of the magnifying lens. This proposal responds to a need for accuracy and reproducibility improvement in the context of the security characterization of integrated circuits using laser shots. The proposed algorithm is based on wavelet decomposition analysis, using a wavelet designed for our specifics.

Other approaches may be adopted to study our IR images, as steerable pyramid transforms, or any image-based approach, but they require some preprocessing. Note that data such as the diagonal details of the DWT decomposition stay untapped, and may also be used for IR images denoising. Despite the optical system limitations arising from the uncooled camera, the proposed approach does not need any image enhancement; its time efficiency makes our autofocus method usable for onthe-fly characterizations.

Fig. 2 :

 2 Fig. 2: Minimalist optical system; the image on plane 2 (green) is well focused on the optical sensor whereas the image on plane 1 is not.

  Fig. 3: Infrared wavelength in the electromagnetic spectrum.

Fig. 6 :

 6 Fig. 6: A simple representation of an IC and its internal structures.

Fig. 7 :

 7 Fig. 7: Two lines delimiting structures boundary (50x).

Fig. 8 :

 8 Fig. 8: Graphs of focus measures of images extracted from videos, based on FISH (red) and on ours (black), for different magnifying factors: 2.5x (a), 5x (b), 20x (c) and 50x (d); the stars mark the real SOI focus positions.

Fig. 9 :

 9 Fig. 9: Graphs of statistics measures of images extracted from videos, following our method, for different magnifying factors: 2.5x (a), 5x (b), 20x (c) and 50x (d); the AC, kurtosis, entropy, AG and SD are in green, magenta, blue, red and black, respectively. The entire graphs are on the left, whereas zoomed in views around the correct focus index are on the right; the stars mark the real SOI focus positions.

TABLE I :

 I Haar wavelet illustration

Scaling function (φ)

Wavelet function (ψ)

TABLE II :

 II Custom Orthogonal Wavelet (COW) illustration

	Scaling function φ	Wavelet function ψ

TABLE III :

 III COW filter bank

	k	Analysis LPF	Analysis HPF	Synthesis LPF	Synthesis HPF
	0	0.35355339	-0.35355339	0.35355339	0.35355339
	1	0.35355339	0.35355339	0.35355339	-0.35355339
	2	0.35355339	-0.35355339	0.35355339	0.35355339
	3	0.35355339	0.35355339	0.35355339	-0.35355339

TABLE IV :

 IV CDF 9/7 9/7 wavelet illustration

	Analysis scaling function φ	Analysis wavelet function ψ

TABLE V :

 V Optimized CDF 9/7 9/7 wavelet illustration

Analysis scaling function φ

Analysis wavelet function ψ

TABLE VI :

 VI Sharpness map in function of magnifying factor and wavelet in use