Adnan Bouakaz

Pascal Fradet

Alain Girault

Symbolic Buffer Sizing for Throughput-Optimal Scheduling of Dataflow Graphs

The synchronous dataflow model is widely used to design real-time streaming applications which must assure a minimum quality-of-service. A benefit of that model is to allow static analyses to predict and guarantee timing (e.g., throughput) and buffering requirements of an application. Performance analyses can either be performed at compile time (for design space exploration) or at run-time (for resource management and reconfigurable systems). However, these algorithms, which often have an exponential time complexity, may cause a huge runtime overhead or make design space exploration unacceptably slow. In this paper, we argue that symbolic analyses are more appropriate since they express the system performance as a function of parameters (i.e., input and output rates, execution times). Such functions can be quickly evaluated for each different configuration or checked w.r.t. many different non-functional requirements. We first provide a symbolic expression of the maximal throughput of acyclic synchronous dataflow graphs. We then perform an analytic and exact study of the minimum buffer sizes needed to achieve this maximal throughput for a single parametric edge graph. Based on these investigations, we define symbolic analyses that approximate the minimum buffer sizes needed to achieve maximal throughput for acyclic graphs. We assess the proposed analyses experimentally on both synthetic and real benchmarks.

I. INTRODUCTION

Synchronous dataflow (SDF) graphs [START_REF] Lee | Synchronous data flow[END_REF] are widely used to design digital signal processing and concurrent real-time streaming applications. This model comes with static analyses that guarantee the boundedness and liveness of an application as well as predictable performances (e.g. throughput, latency, memory requirements).

Performance analyses of SDF graphs check whether nonfunctional requirements are met. They can be performed both at design time and at run-time. At design time, it is a crucial step in the development of embedded applications. Many decisions and settings of the system need to be explored (e.g. hardware/software partitioning, memory allocation, granularity and different implementations of tasks, processor speeds, etc.) and the best options that satisfy the non-functional requirements can be chosen. At run-time, performance analysis is performed either for resource management or to cope with the dynamic behavior of parametric versions of SDF. Indeed, in response to the increasing complexity of systems, many parametric dataflow models have been proposed (e.g., PSDF [START_REF] Bhattacharya | Parameterized dataflow modeling of DSP systems[END_REF], SPDF [START_REF] Fradet | SPDF: a schedulable parametric data-flow MoC[END_REF], BPDF [START_REF] Bebelis | BPDF: a statically analyzable dataflow model with integer and boolean parameters[END_REF], πSDF [START_REF] Desnos | PiMM: parametrized and interfaced dataflow meta-model for MPSoCs runtime reconfiguration[END_REF], SADF [START_REF] Damavandpeyma | Parametric throughput analysis of scenario-aware dataflow graphs[END_REF], etc.) in which the graph (e.g., its rates or channels) may change at run-time.

Throughput is one important timing constraint of realtime systems. For example, a video decoder should decode a minimum number of frames per second. Using a throughputoptimal scheduling policy, such as self-timed scheduling, allows the designer to guarantee such timing requirements, but also to use dynamic voltage and frequency scaling (DVFS) techniques to reduce energy consumption in case the maximum throughput is larger than the desired quality-of-service [START_REF] Stuijk | Throughput-constrained DVFS for scenario-aware dataflow graphs[END_REF]. Furthermore, for embedded systems, it is primordial to determine the minimum buffer sizes that allow such scheduling.

In this paper, we propose symbolic analyses of dataflow graphs. Our analyses consider the rates and execution times of actors as parameters. Most non-functional properties of the application can be described as a function of these parameters. By evaluating these functions for specific values of parameters, the properties/performances of specific configurations can be obtained. We propose two symbolic analyses of acyclic graphs under self-timed scheduling to answer the following questions:

Q1. What is the maximum throughput of the application? Q2. What are the minimum buffer sizes that allow the application to achieve its maximum throughput? Although our symbolic analyses may give only approximate (but safe) results, they are quite useful in many cases (see Fig. 1) :

1. At early design stages, the SDF model of the application is usually partially specified and many system settings remain to be explored. For instance, the execution time of an actor depends on whether it is implemented in software or in hardware, or which algorithm will be used to implement its functionality, etc. Therefore, design space exploration may require analyzing of a potentially huge number of configurations (path 3 ). Symbolic analyses are a big advantage in this case. Formulas are generated only once and simply evaluated for each possible configuration (i.e., set of parameters) (path 4 ).

2.

Similarly, non-functional requirements of parametric dataflow models at compile time can be expressed symbolically 1 as parametric formulas. Then, the requirements can be either checked by evaluating formulas for all potential configurations (path 6b ) or, better, by an analytic proof (path 6a ) to ensure, for instance, that the buffering requirements are less than some threshold whatever the configuration of the application.

3.

For dynamic models and runtime resource management, appropriate settings have sometimes to be chosen dynamically. Consider, for instance, a parametric application where frequency scaling is used to guarantee a specific throughput and minimize power consumption. In such case, frequency must be adjusted at each parameter change. Instantiating the graph (path 5 ) and performing a numerical analysis is far too costly at run time. Consequently, fast analyses, like the evaluation of symbolic formulas, are required (path 6b ).

4.

Finally, even for completely static SDF models, many analyses have an exponential complexity. In practice, the exact computation of throughput and latency is acceptable at compile time (path 1 ). However, exact algorithms for minimal buffer sizes are too expensive even for small graphs. Moreira et al. [START_REF] Moreira | Buffer sizing for rate-optimal single-rate data-flow scheduling revisited[END_REF] shows that this problem is NP-complete for homogeneous SDF (HSDF) graphs. Moreover, SDF-to-HSDF conversion may lead to an exponential growth of the size of the graph. Our symbolic analysis (path 2 ) is much more efficient and its approximate solution can also be considered as a starting point to prune the parameter space and hence improve the performance of the exact algorithm.

Contributions. The contributions of our paper are threefold: a. We present and prove a duality theorem and show how it can be used to prune many cases.

b. We give a new analytic characterization of the parametric data-dependency A p q --→ B, called enabling patterns, which can be used to build symbolic analyses, e.g., for buffer sizing. c. We describe exact and approximate polynomial-time symbolic analyses to answer the above questions (Q1 and Q2).

The article is organized as follows. Section II introduces synchronous dataflow graphs and the needed definitions. Section III presents throughput analysis and a new generic result required to solve the second question. We present in Section IV an exact analysis of the minimum buffer size for a single edge SDF graph A p q --→ B. These results are extended to acyclic graphs using approximate analyses in Section V. We evaluate our technique on both synthetic and real benchmarks in Section VI. Finally, we review related work in Section VII and conclude in Section VIII. The interested reader will find proofs and additional details in a companion paper [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF].

II. BACKGROUND

In this section, we first introduce the application model and the scheduling policy. Then, we review some useful properties.

Application model: An SDF graph G = (V, E) consists of a finite set of actors V and a finite set of edges E which can be seen as FIFO channels. The atomic execution of a given actor (called firing) consumes data tokens from all its incoming edges (its inputs) and produces data tokens to all its outgoing edges (its outputs).

The number of tokens consumed (resp. produced) at a given port at each firing is called the consumption (resp. production) rate. An actor can fire only when all its input edges have enough tokens (i.e., at least the number specified by the corresponding rate). An edge may contain some initial tokens (also called delays). The execution time of an actor X is denoted by t X . For instance, Fig. 2(a) shows an SDF graph with two actors A and B, with execution times t A = 20 and t B = 7 respectively. The production and consumption rates are 8 and 5 respectively. The top edge in Fig. 2 The state of a dataflow graph is the vector of number of tokens present at each edge (i.e., buffered in each FIFO channel). Each edge carries zero or more tokens at any moment. The initial state of the graph is specified by the number of initial tokens. The initial state of the graph of Fig. 2

(b) is [i AA =1, i AB =0, i BB =1, i BA =15
]. An iteration of an SDF graph is a non empty sequence of firings that returns the graph to its initial state. For the graph in Fig. 2(a), firing actor A five times and actor B eight times forms an iteration. The repetition vector z = [z A =5, z B =8] indicates the number of firings of actors per iteration. If such vector exists the graph is said to be consistent [START_REF] Lee | Synchronous data flow[END_REF]. We denote by z X the number of firings of actor X in the iteration.

Homogeneous SDF (HSDF) is a restriction of SDF where all the production and consumption rates are equal to 1. HSDF graphs are particularly useful because (i) their throughput can be computed as the inverse of the maximal cycle mean (MCM); and, (ii) any SDF graph can be converted into an equivalent HSDF graph. The cycle mean of a cycle is equal to the sum of execution times of the actors in the cycle divided by the number of delays (i.e., initial tokens) in the channels of this cycle. This provides a way to compute the throughput of any SDF graph. Unfortunately, the translation from SDF to HSDF may lead to an exponential increase of the number of nodes.

Scheduling policy: In this paper, we focus on as soon as possible (ASAP) scheduling of consistent graphs without autoconcurrency (i.e., two firings of the same actor cannot overlap). In such self-timed executions [START_REF] Sriram | Embedded multiprocessors: scheduling and synchronization[END_REF], an actor fires as soon as if it is idle (no auto-concurrency) and has enough tokens on its input channels. We assume there are sufficient processing units, e.g., there are as many processors as actors or actors are implemented in hardware. ASAP scheduling allows to reach the maximal throughput. Such schedules are naturally pipelined and composed of a prologue followed by a steady state that repeats infinitely. Fig. 6 shows the ASAP schedule of the graph in Fig. 2(a).

The multi-iteration latency of the first n iterations of the execution of a graph G, written L G (n), is equal to the finish time of the last firing in all firings composing the first n iterations (assuming timing starts at the very first firing). The period of the execution (denoted by P G ) is the average length of an iteration and formally defined as

P G = lim n→∞ L G (n)
n The throughput of the execution, T G , is the number of iterations per unit of time, and hence equals 1 P G .

The FIFO channel A p q --→ B with bounded size d can be modeled by adding a backward channel B q p --→ A with d initial tokens, as shown in Fig. 2(b). This conservative modeling, assumed in most works, enforces that an actor can start firing only if there is enough space on its output channels. Fig. 2(b) also shows how to make the prevention of auto-concurrency explicit by adding self-edges with rates equal to 1 and one initial token. For a channel A p q --→ B with d initial tokens, the i th firing of B (denoted B i ) is enabled if and only if the number of produced tokens is larger than iq. Hence, B has to wait for the j th firing of A (denoted A j ) such that jp + d ≥ iq. We thus have:

B i ≥ A j with j = iq -d p (1)
This equation characterizes the data-dependency between A and B. However, because of the ceiling function, this equation is not suited to symbolic manipulations. We propose in Section IV a new characterization that is more intuitive and suitable to reason about buffer sizes.

III. THROUGHPUT AND DUALITY

In this section, we first determine the exact maximal throughput for acyclic SDF graphs. Then, we present a new property that will be used to answer the minimum buffer sizes question.

Property 1 (Throughput). The maximal throughput of an acyclic SDF graph G = (V, E) is equal to

T G = 1 max A∈V {z A t A } (2)
Hence, the minimal period is

P G = max A∈V {z A t A }.
Proof: This is easily shown by considering the MCM analysis of the corresponding HSDF graph. The only cycles in HSDF(G), the HSDF graph equivalent to an acyclic graph G, are those used to represent the infinite firings of the same actor (see Fig. 3(a)). For each actor A, its corresponding cycle contains one delay and z A instances of A. Thus, the cycle mean is equal to z A t A , the MCM is equal to max A∈V z A t A and denotes the inverse of the maximal throughput of HSDF(G) and G.

We say that actor A imposes a higher load than actor B when z A t A > z B t B . The throughput and period of an acyclic graph is defined by the actor which has the highest load. This implies that this actor never gets idle once the execution enters the steady state. The dual of an SDF graph G, denoted G -1 , is obtained by reversing all edges of G.

Theorem 1 (Duality theorem). Let G be any (cyclic or not) live graph and G -1 be its dual, then

T G = T G -1 and ∀i. L G (i) = L G -1 (i).
Proof: The detailed proof can be found in [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF]. We first prove that the graph HSDF(G) is the dual (after renaming actors) of HSDF(G -1 ) (see Fig. 3). Therefore, both HSDF graphs have the same MCM and, by Eq. ( 2) T G = T G -1 . We prove that ∀i. L HSDF(G) (i) = L HSDF(G -1 ) (i) by unfolding both HSDF graphs for i iterations.

Fig. 4 illustrates the duality theorem with the SDF graph G a of Fig. 2(b). The latency of the first iteration of the ASAP execution of that graph is equal to the latency of the first iteration of its dual i.e., L Ga (1) = L Ga -1 (1).

We use the transformation of a graph to its dual as well as the associated theorem at several occasions during the analysis of minimal buffer sizes.

IV. THE PARAMETRIC GRAPH

A p q --→ B
This section focuses on the simplest SDF graph made of a single edge: G = {A p q --→ B}. It provides exact answers for the minimum buffer sizes needed to achieve maximal throughput. Graph G is parametrized by production and consumption rates 1 p, q ∈ N + as well as execution times t A , t B ∈ R + . This section shows that the symbolic analysis, even for such simple graphs, is quite involved.

The balance equation z A p = z B q entails that the repetition vector of this graph is:

[z A = q/ gcd(p, q), z B = p/ gcd(p, q)]
and, according to Property 1, its period is:

P G = max(z A t A , z B t B ) (3)

A. Enabling patterns

We introduce enabling patterns which characterize the data-dependency between a producer and a consumer. Compared to Eq. ( 1), they are more intuitive and suitable to reason about buffer sizes.

Enabling patterns between the producer A and consumer B are defined by the following grammar:

P ::= A i B j | [P ] x=1..k | P 1 ; P 2
where i, j, k evaluate to a positive integer.

An enabling pattern P is either a basic pattern (A i B j ), a repetition for k times ([P ] x=1..k ), or a sequence of enabling patterns P 1 ; P 2 , The expressions i, j or k are arithmetic

Case A. p ≥ q Let p = kq + r with 0 ≤ r < q Case A.1. r = 0 A B k (4) Case A.2. q ≤ 2r A B k ; A B k+1 αj j=1•• q-r gcd(p,q) (5) Case A.3. q > 2r A B k βj ; A B k+1 j=1•• r gcd(p,q) (6)
where αj = jr q-r -(j-1)r q-r and βj = jq r

-(j-1)q r -1. Case B. p < q Let q = kp + r with 0 ≤ r < p Case B.1. r = 0 A k B (7) Case B.2. p ≥ 2r A k+1 B ; A k B γj j=1•• r gcd(p,q) (8) Case B.3. p < 2r A k+1 B λj ; A k B j=1•• p-r gcd(p,q) (9)
where γj = jp expressions made of integers, parameters or pattern variables defined by enclosing repetition patterns.

The semantics of an enabling pattern is defined w.r.t. two counters c A and c B representing the number of completed firings of A and B (initially 0). The pattern A i B j ; P w.r.t. (c A , c B ) means that:

• c A firings of A have produced enough tokens to fire actor B exactly c B times; • then, if A is not fired at least i times more then B cannot be fired; otherwise B can be fired j times; • the subsequent pattern P is considered with the new values (c A + i, c B + j).

A repetition

[P ] x=1..k is a sequence of k patterns P . The pattern [P ] x=1..k is also written [P ] k if the pattern variable x is not used in P .
A correct enabling pattern must cover an entire iteration: at the end of the pattern, we should have c A = z A and c B = z B . For instance, the enabling pattern of A --→ B is:

A B ; A B 2 ; A B ; [A B 2 ] 2
which is illustrated in Fig. 6. It can also be written as:

A B ; [A B 2 ] i i=1••2
This representation is very useful when the length and shape of enabling patterns depend on the parameters of the graph.

B. Enabling patterns of A p q

--→ B Depending on the production and consumption rates p and q, there are six possible enabling patterns. Proof: We show here the proof of cases (A.1) and (A.2). The remaining cases can be found in [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF].

• Case (A.1) If r = 0 then p = kq and the repetition vector is z = [1, k]. Each firing of A enables k firings of B and the enabling pattern is A B k .

• Case (A.2)(p = kq + r and q ≤ 2r) The first firing of A enables only k firings of B, yielding the pattern A B k . The number of remaining tokens in the channel after sequence AB k is equal to r. Since 2r ≥ q, the second firing of A enables k + 1 firings of B, yielding the aggregated pattern A B k ; A B k+1 . It follows that pattern A B k+1 can be repeated a number of times denoted α 1 , after which there will be r + α 1 (r -q) tokens left. So, α 1 is the largest integer such that r + α 1 (r -q) ≥ 0. Hence, α 1 = r q-r . The next firing of A will only enable k firings of B (A B k ). By the same reasoning as above, this will be followed by [A B k+1 ] α2 where α 2 is the largest integer such that r+α 1 (r-q)+r+α 2 (r-q) ≥ 0. Hence, α 2 = 2r q-r -r q-r . This process is repeated until all firings of A and B of the iteration take place, yielding the complete pattern

A B k ; [A B k+1 ] αj j=1••m
for some m that remains to compute. The number of firings of A in this pattern is equal to c A = m j=1 1 + α j , while the number of firings of B is equal to c B = m j=1 k+(k+1)α j . A correct pattern should cover the entire iteration, i.e., c A = z A and c B = z B , which implies (details omitted, see [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF]) m = q-r gcd(p,q) . The enabling patterns of Fig. 5 can be used to derive a quasi-static schedule for the parametric graph A p q --→ B.

C. Minimum buffer size of

A p q --→ B
We now use enabling patterns to compute the minimum size of the buffer A p q --→ B (denoted θ A,B ) such that the ASAP execution achieves the maximal throughput or, equivalently, the minimal period given by Eq. ( 3). The buffer size is modeled by adding a backward edge with θ A,B initial tokens. We distinguish two cases: Enabling patterns are the key to find the minimum number of tokens that allows this behavior. A trivial case is (A.1) where p = kq and the enabling pattern is A B k . In order to perform the first two firings of A consecutively, the backward edge should have at least 2p tokens. Since qt A ≥ kqt B (hence t A ≥ kt B ), the k firings of B complete before the third firing of A which still needs 2p initial tokens in order to fire again immediately. Hence, the minimum buffer size is 2p.

• Case z A t A ≥ z B t B (i.
• Case z A t A < z B t B (i.e., qt A < pt B ): Actor B has the highest load and should fire consecutively for maximal throughput. However, in general all firings of B cannot be consecutive since initially, there is no token on the edge. The previous approach can still be followed using the duality theorem. Since the graph G and its dual G -1 have the same throughput, we can apply the former reasoning on G -1 where B is the producer and has the highest load.

Property 3. If z A t A ≥ z B t B ,
the minimum buffer sizes of A p q --→ B for maximal throughput are given by the symbolic formulas of Fig. 7.

Case I. Case I.1. A.1 ∨ ((A.2 ∨ A.3) ∧ t A ≥ (k + 1)t B ) θ A,B = 2p + q -gcd(p, q) ( 10
) Case I.2. B.1 ∨ ((B.2 ∨ B.3) ∧ t B ≤ kt A ) θ A,B = p + q -gcd(p, q) + t B t A p (11)
Case II.

Case II.1.

(A.2 ∧ r ≥ r q-r r q-r +1 t B ) ∨ (A.3 ∧ r ≥ 1 q r t B ) where r = t A -kt B θ A,B = 2p + q -gcd(p, q) + t B -r r r (12)
Case II. (jr mod (q-r)) [START_REF] Sriram | Embedded multiprocessors: scheduling and synchronization[END_REF] where n is the smallest positive integer such that Proof: Here, we only outline the proof and focus on one case ((I.1) with (A.2)). A detailed proof can be found in [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF].

Let δ j be the minimum number of tokens in the backward edge (representing the buffer) such that the j th firing of A can occur immediately after the (j -1) th firing of A. By definition of θ A,B , we have θ A,B = max j δ j . Let x j denote the number of firings of B that have finished by the start of the j th firing of A. Hence, δ j = jp -x j q.

The three cases of Fig. 7 should be read as (I) else (II) otherwise (III). CASE (I): At any given enabling point (i.e., any in the enabling pattern), all newly enabled firings of B complete their execution before the next enabling point.

In terms of the enabling pattern cases identified in Fig. 5, case (I) must be split into two exclusive subcases, (I.1) when p ≥ q and (I.2) otherwise. The conjunction of the condition for case (I) with p ≥ q yields the condition A.1 ∨ ((A.2 ∨ A.3) ∧ t A ≥ (k + 1)t B ), which is shown in Fig. 7. The other conditions are obtained similarly.

Let us prove the result for case (I.1) and pattern (A.2), i.e., p = kq + r ∧ q ≤ 2r. According to the enabling pattern (Eq. ( 5)), at most (k+1) firings of B can be enabled at a given point. They all run in parallel with one firing of A. Case (I) requires that t A ≥ (k + 1)t B .

We first expand the enabling pattern of Eq. ( 5) into an infinite pattern in order to compute the minimum buffer size θ A,B over the infinite execution of the graph:

A B k ; A B k+1 αj j∈N + (15)
For the sake of the proof, since α j ≥ 1, we can rewrite Eq. ( 15) into the equivalent infinite pattern:

A B k ;A B k+1 prologue ; [A B k+1 ] αj -1 ;A B k ;A B k+1 block j j∈N + (16)
The minimum of tokens to enable the first firing of A must be δ 1 = p. Then, δ 2 = δ 1 + p = 2p because no B has finished before the second A starts. Then,

δ 3 = δ 2 + (p -kq) = 2p + r. Subsequently, δ 4 = δ 3 + (p -(k + 1)q) = δ 3 + (r -q),
and all α 1 subsequent values of δ i are obtained from δ i-1 by adding (r -q) which is negative. The sequence (α j ) is defined in Fig. 5. Hence, δ 1+(α1+1)+1 = δ 3 +α 1 (r-q) = 2p+r+α 1 (r-q). This ends at the last of the α 1 patterns A B k+1 , so the next value δ 1+(α1+1)+2 is obtained by adding p -kq = r because of pattern A B k , yielding 2p + r + α 1 (r -q) + r. The computation of the infinite sequence (δ j ) is illustrated in Fig. 8. Within each block j, the subsequence (δ h ) is strictly decreasing because (r -q) is negative, so its maximum value is the value of the entry point, which we denote by j = 1 + j i=1 (α i + 1) + 2. We thus have:

δ j = 2p+r+ j i=1
(α i (r-q)+r) = 2p+r+jr+(r-q) jr q -r It follows that the maximum value of the infinite sequence (δ j ) is:

θ A,B = max j∈N δ j = max j∈N (2p + r + jr + (r -q) jr q -r ) = 2p + r + (q -r) max j∈N ( jr q -r - jr q -r )
As a conclusion, θ A,B = 2p + r + (q -r -gcd(p, q)) = 2p + q -gcd(p, q). CASE (II): Case (I) is not satisfied, but, for any block (e.g.,

[[A B k+1 ] αj ; A B k ] in case (A.2)
), all firings of B during this block complete their execution before the first enabling point in the next block.

This case is illustrated in Fig. 9. A block is of the form

[A B 2 ] αj ; A B .
The maximum value of α j is 2. Therefore, five firings of B have to run in parallel with three firings of A. So, we must have 5t B ≤ 3t A . The computed sequence (δ j ) in case (II) is similar to that of case (I) but with small shifts. block block Fig. 9: An ASAP execution (p = 8, q = 5, t A = 13, t B = 7).

CASE (III):

Otherwise. This is the most complicated case to solve since the sequence (x j ) does not follow the enabling patterns. Our solution is based on the following observations. We define a catch-up sequence as a sequence of consecutive firings of B (i.e., without gaps) that may spread over many blocks. Fig. 10 illustrates a catch-up sequence over two blocks.

catch-up sequence

Enable: Finish before: The key observation is the following. For the firings of A inside a catch-up sequence, the number of firings of B that finish before firings of A actually follows the enabling pattern of graph A t A t B ----→ B, i.e., as if time was produced and consumed instead of tokens. Furthermore, the maximum of sequence (δ j ) occurs inside the maximal (in terms of blocks) catch-up sequence. For instance, n in Eq. ( 14) represents the length of the maximal catch-up sequence.

In Case (III), when actors A and B impose the same load (i.e., z A t A = z B t B ), the catch-up sequence spreads all over the iteration. In this worst-case scenario, all four cases (III.A.2, III.A.3, III.B.2 and III.B.3) give the same upper bound:

θ u A,B = 2(p + q -gcd(p, q)) (17)
This bound is also tight, in the sense that for all p, q, there exist t A and t B such that θ A,B given in Eq. ( 14) is equal to θ u A,B . This upper bound does not depend on the execution times of the actors. Therefore, it can be used as a safe buffer size if the execution times of actors are unknown.

Our last result concerns the minimum buffer size of the graph A p q --→ B when there are d initial tokens.

Property 4. If If z

A t A ≥ z B t B and the channel A p q
--→ B contains d initial tokens, then the minimum buffer size θ A,B that allows the maximum throughput is

θ A,B = max{0, θ A,B -d + d mod gcd(p, q)} (18)
with θ A,B as defined in Fig. 7 Proof: See [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF].

If z A t A < z B t B the minimum buffer sizes are computed identically but on the dual graph.

V. ACYCLIC GRAPHS Exact symbolic analyses for a single edge are already so complex that they seem to be out of reach for arbitrary (even acyclic) graphs. This section shows how to use the previous results to obtain approximate analyses for the minimum buffer sizes of general acyclic dataflow graphs. To achieve this, we make use of a technique that linearizes the firings of actors. --→ B}, where, as illustrated in Fig. 6, the firings of A are consecutive while those of B are neither consecutive nor uniformly distributed. Let f B (i) denote the finish time of the i th firing of actor B. In order to derive formulas that can be composed (e.g., when dealing with a chain of actors), we want to transform B into a fictive actor B u that fires consecutively as many times as B and such that

A. Linearization of graph A

∀i. f B (i) ≤ f B u (i)
Actor B u has a starting time t 0 B u and an execution time t B u , and since it fires consecutively f B u (i) = it B u + t 0 B u . We present one linearization method, called Stretch , illustrated in Fig. 11. Method Stretch determines B u by increasing the execution time of B in order to fill the gaps. We distinguish two cases:

• Case qt A ≥ pt B : We take t B u = qt A
p . The starting time can be shown (see the proof in [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF]) to be t 0

B u = t A +t B -gcd(p,q) p t A . So, ∀i. f B u (i) = qt A p i + t A + t B - gcd(p, q) p t A (19)
Method Stretch may advance the starting of some firings, but always postpone their endings. • Case qt A < pt B : Firings of B are consecutive in the steady state. Therefore, we can take t B u = t B and using the duality theorem, we have L G (n) = L G -1 (n) = nz B t B + ∆ B,A and we can take t 0 B u = ∆ B,A , where ∆ B,A is computed on the schedule of the first iteration of the dual graph B q p --→ A as the difference between the end of the last firing of A and the end of the last firing of B (see [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF] for details).

In both cases, it can be shown that this linearization is tight in the sense that ∃i. f B (i) = f B u (i). We can see in Fig. 11 that both 5th firings of B and B u finish at the same time.

Thanks to this linearization, the chain A p q ---→ B p q ----→ C can be treated by first scheduling the subgraph A p q --→ B, then linearizing the firings of B if they are not consecutive, then scheduling the subgraph B u p q ---→ C, and finally combining the two schedules.

B. Minimum buffer sizes for maximal throughput

We give approximate minimum buffer sizes that allow general acyclic graphs to reach their maximal throughput. We first present formulas to compute safe upper bounds for general acyclic graphs, then we present a heuristic that improves this bound for chains, trees, and in-trees (DAGs with only joins). These special kinds of graphs, especially chains, are common in streaming applications.

1) Safe upper bounds:

We first present a negative result. Let G be an acyclic graph and let the size of each channel A p q --→ B be equal to θ A,B as defined in Section IV-C. These buffer sizes do not always permit maximal throughput 2 . A simple counterexample is the graph G b = {A Locally, these buffer sizes allow the producers to run freely without any constraint from the consumers. However, when they are put together, the maximal throughput, where actor C fires consecutively, cannot be achieved (see Fig. 12). The computation of θ B,C assumes that the execution time of actor B is t B = 11. However, as illustrated in Fig. 12, there are gaps between the firings of B due to the data-dependency A → B. The global execution proceeds as if the execution times of B were sometimes longer than 11. --→ B in G is at least 2(p + q -gcd(p, q)), then the ASAP execution of the graph achieves the maximal throughput.

Proof: We first present the proof for chains. Let G be the chain {A

1 p1 q1 ----→ A 2 p2 q2 ----→ A 3 → • • • → A n }, according to Eq. (2), the minimal period of G is P G = max i=1..n {z Ai t Ai }.
2 They do however allow maximal throughput in some specific cases described in the next section

The period and therefore the throughput remain the same if the execution time of each actor A i is considered to be P G z A i . Let G = be the version of G where all actors have the same load as the maximum load in G. Then G and G = have the same period and throughput.

If the size of each buffer

A i pi qi ---→ A i+1 in G = is θ u Ai,Ai+1 = 2(p i + q i -gcd(p i , q i ))
, then G = still achieves the maximal throughput. Indeed, size 2(p 1 + q 1 -gcd(p 1 , q 1 )) for the first channel allows both A 1 and A 2 to run consecutively in the steady state (see Eq. ( 17)). Similarly, size 2(p 2 + q 2 -gcd(p 2 , q 2 )) for the second channel allows both A 2 and A 3 to run consecutively, and so on.

Since graph G = with these buffer sizes achieves the maximal throughput, reducing the execution times of actors in G = to their original values will never decrease the throughput of the graph thanks to the monotonicity of the self-timed execution. Hence, graph G with these buffer sizes achieves the maximal throughput.

The proof for general graphs with no undirected cycle can be found in [START_REF] Bouakaz | Symbolic analysis of dataflow graphs (extended version)[END_REF].

Note 1: Since the minimum buffer sizes below which the graph is definitely not live are equal to p+q -gcd(p, q) [START_REF] Battacharyya | Software synthesis from dataflow graphs[END_REF], Property 5 provides a first solution which is less than twice the exact one. For parametric dataflow models, the upper bound θ u A,B can actually be reached for some configurations. However, if the system supports dynamic reallocation of memory, it is still useful to evaluate the minimal buffer sizes to adjust buffers after each configuration change.

Unfortunately, Property 5 does not hold for general acyclic graphs that contain undirected cycles. A counterexample is the graph

G c = {A 4 3 ---→ B 3 8 ---→ D, A 1 3 ---→ C 3 2
---→ D} with t A = 4, t B = 3, t C = 12 and t D = 8. The repetition vector is z = [6, 8, 2, 3] and all actors impose the same load (i.e., ∀X. z X t X = 24). The ASAP execution when all buffer sizes are equal to their upper bound 2(p + q -gcd(p, q)) is shown in Fig. 13. Actor A does not fire consecutively so the throughput is not maximal. The reason is that the chain A C D imposes an earliest start time for D that is after the earliest start time imposed by the chain A B D. More precisely, the first firing of actor D is delayed by actor C (i.e., by the second chain), which delays the 7 th firing of B which in turn delays the 8 th firing of A. Let f u D,1 (resp. f u D,2 ) denote the linear upper bound on finish times of actor D following the first (resp. second) chain. We have

f u D,1 (i) = 8i + 16 and f u D,2 (i) = 8i + 28, hence f u D,2 (i) > f u D,1 (i).
In order to prevent the second chain from impacting the schedule of the first chain, we must increase the size of buffer B D so that B may fire during 28 -16 time units. Since B produces 3 tokens and t B = 3 the size of the B D buffer must be increased by 28-16 3 3 = 12.

Fig. 13: ASAP execution of G c

In the general case, for a chain A 1

p1 q1 ----→ A 2 → • • • pn-1 qn-1
--------→ A n where all actors have the same load, we have

f u An (i) = t An i + z An t An n-1 i=1 p i + q i -gcd(p i , q i ) p i z Ai (20)
If the actors do not have the same load, we consider, as in the proof of Property 5, the chain where all actors have the same load i.e., the maximum load of the original chain.

Property 6. Let two different chains from A 1 to A n such that f u An,1 (i) = t An i + x 1 , f u An,2 (i) = t An i + x 2 and x 1 < x 2 .
To prevent the second chain from disturbing the schedule of the first one, it suffices to increase the size of the last channel

A n-1 p q --→ A n of the first chain by ζ with ζ = x 2 -x 1 t An-1 p (21)
Proof: The proof is based on the following observation. Let A p q --→ B be a graph with two actors such that z A t A = z B t B . Hence, θ A,B = 2(p + q -gcd(p, q)). Even if the first firing of B starts only at time t A p (p + q -gcd(p, q)), actor A can still fire consecutively provided that the buffer size is θ u A,B (the proof uses Eq. ( 19), details omitted). Hence, if all actors in the first chain are delayed as indicated by Eq. ( 20), then the sizes 2(p i + q i -gcd(p i , q i )) still allow the first actor to fire consecutively. Note that after introducing these delays, no actor gets idle once it starts executing. Let x 1 and x 2 denote the start time of the last actor in the two different chains from A 1 to A n . If x 2 > x 1 , the extra delay (x 2 -x 1 ) imposed by the second chain will hinder the previous property. Let A n-1 p q --→ A n be the last channel in the first chain. Increasing the size of this channel by x2-x1 t A n-1 p will allow actor A n-1 to fire consecutively during the extra delay (x 2 -x 1 ). That is, the impact of the second chain on the first one is avoided.

Note 2:

The value of ζ, like the value of θ u A,B , does not actually depend on execution times. Indeed, Eq. (20) shows that x 1 and x 2 can be expressed as z An t An k 1 and z An t An k 2 . Since z An t An = z An-1 t An-1 , term x2-x1 t A n-1 can be expressed as z An-1 (k 2 -k 1 ) which does not depend on execution times. This approach can be extended to deal with any acyclic graph. The final property should be that for any two different chains from A 1 to A n in the graph such that

f u An,1 (i) = t An i + x 1 , f u An,2 (i) = t An i + x 2 and x 1 < x 2 the size of the last channel A n-1 p q --→ A n of the "fastest" chain is at least equal to 2(p + q -gcd(p, q)) + x2-x1 t A n-1 p.
This property can be enforced by computing, for each node X and each chain from a source actor to X, the function f u X,j (i) = t X i + x j . Then, as in Property 6, we add to each final edge of these chains the value ζ computed w.r.t. the maximum computed x j . Such an approach is safe but not always needed (e.g., when the predecessors of a node do not have common ancestors). We do not describe the refined and optimized algorithm here.

2) Improving the upper bounds:

In this section, we improve the minimum buffer sizes for chains, trees, and in-trees, starting with chains. We say that a chain is monotone if each actor imposes a higher load than its successor or if each actor imposes a lower load than its successor.

Definition 1. The chain

A 1 → • • • → A n is monotone if and only if (∀i. z Ai t Ai ≥ z Ai+1 t Ai+1 ) ∨ (∀i. z Ai t Ai ≤ z Ai+1 t Ai+1 )
Let θ Ai,Ai+1 be the size of the buffer between A i and A i+1 as computed in Section IV-C. This size allows the single edge graph A i pi qi ---→ A i+1 to reach its maximal throughput. Property 7. A monotone chain A 1 → • • • → A n where the size of each buffer A i → A i+1 is at least θ Ai,Ai+1 achieves its maximal throughput.

Proof: Suppose that the chain has the descending order ∀i. z Ai t Ai ≥ z Ai+1 t Ai+1 . The size θ A1,A2 allows actor A 1 to run consecutively, which is the same behavior as when there is no constraint on buffer sizes. Then, the size θ A2,A3 allows actor A 2 to fire as soon as it is enabled by actor A 1 . Actually, if there were no dependence from A 1 to A 2 , the size θ A2,A3 would allow actor A 2 to run consecutively. The same reasoning shows that all actors are fired as it there were no buffer size constraints. Therefore, the chain achieves its maximal throughput. The case of ascending order can be easily dealt with using the duality theorem.

Note that Property 7 is only a sufficient condition simply because θ Ai,Ai+1 allows actor A i to fire consecutively. However, as soon as when i > 1, it is not always necessary in order to achieve the maximal throughput. E.g., let

G d = {A 1 2 2 ---→ A 2 4 3
---→ A 3 } such that t A1 = 28, t A2 = 20 and t A3 = 15. Though θ A2,A3 = 12, the minimum size of channel A 2 → A 3 is actually 9. Indeed, as illustrated in Fig. 14, even if this size delays firings of A 2 (see the 4 th firing), the introduced delay does not prohibit A 1 from firing consecutively. Property 7 also holds for non monotone chains made of an ascending sub-chain followed by a descending one. We say that those chains are of the form . The computed buffer sizes on both wings allow the actors at the top to run consecutively.

Unfortunately, as illustrated in Fig. 12 in the previous section, Property 7 does not hold for any chain. Our solution is hence to put any chain on the form by using the same approach as proof of Property 7; i.e., by increasing the execution times of some actors (without exceeding the maximum load P G ), then computing the buffer sizes as in Property 7, and finally restoring the original execution times. Fig. 15 illustrates this solution.

Any chain on the form obtained by increasing the load of the actors of the original chain is a valid solution. For example, the chain A B . . . K can be transformed in the red chain which is of the form . Actually, any chain of this form inside the gray area is a valid solution. However, we can choose one that minimizes the buffer sizes. According to Section IV-C, the channel I → J is in case (I.1) (Fig. 7) and the size of the channel is 4 whatever the chosen load. The channel J → K is in case (III) and increasing the load of J may put it in case (I.1) where the buffer size is smaller. Increasing loads serve to load Fig. 15: Transformation of a chain to a form.

put the channel in a less expensive case as described in Fig. 7. Furthermore, as long as the load remains in the gray area, several choices can be made e.g., to facilitate the treatment of subsequent channels. There are many heuristics using this idea to minimize buffer sizes. They involve other criteria such as the expected size gain that can be used to prioritize the treatment of channels. We do not describe them here but one of them is evaluated in the experiments section (see Section VI).

The case of trees is solved in the same way. If the tree does not contain any sub-trees (i.e., it consists of a set of chains originating from the same root node), then the load of the root node is first increased to be equal the maximum of all loads in the tree and then the previous method can be applied on every chain composing the tree. This is correct because the computed buffer sizes will allow the root actor to run consecutively. If the tree contains sub-trees, the same process is first applied recursively on sub-trees and then it proceeds with each subtree replaced by its root node.

VI. EXPERIMENTS

In this section, we compare the results of our heuristic presented in Section V-B2, with the safe upper bounds (Section V-B1) and the exact minimum buffer sizes using many randomly generated SDF graphs and some real benchmarks. First, we compare the results of the heuristic with the safe upper bounds 2(p + q -gcd(p, q)) using two million randomly generated chains of 10 actors where production and consumption rates (resp. execution times) are uniformly distributed over the interval [START_REF] Battacharyya | Software synthesis from dataflow graphs[END_REF]20] (resp. [START_REF] Battacharyya | Software synthesis from dataflow graphs[END_REF]200]). The number of firings per iteration, X z X , of every generated chain has been bounded by 6 × 10 3 . For each graph, we compute the ratio of the total buffers sizes obtained by our approach to the sum of the lower bounds, p + q -gcd(p, q) 3 . Each black dot in Fig. 16 represents the obtained ratio for one graph, while the 3 Without lost of generality, tokens are assumed to have the same size.

red line represents the average of ratios. Fig. 16 shows that, in average, our heuristic reduces the total buffer sizes by 20% compared to the upper bounds.

Secondly, we compare the results of the heuristic with the exact minimum buffer sizes. The exact solution is obtained using a dichotomic search that first checks (using symbolic simulation) whether there is a channel capacities distribution, whose total size is at the middle between the lower bound and the already obtained solution by the heuristic, which allows the graph to achieve its maximal throughput. Depending on the answer of this request, the algorithm proceeds to either the top or bottom part of the search space, and so on. Due to the exponential complexity of the minimum buffer sizes problem, we evaluate our approach on only 10 4 randomly generated chains of four actors where production and consumption rates (resp. execution times) are uniformly distributed over the interval [START_REF] Battacharyya | Software synthesis from dataflow graphs[END_REF][START_REF] Ghamarian | Parametric throughput analysis of synchronous data flow graphs[END_REF] (resp. [START_REF] Battacharyya | Software synthesis from dataflow graphs[END_REF]100]). The blue (resp. red) line in Fig. 17 represents the average of the ratios of the exact (resp. approximate) solution to the lower bound. So, in average, our heuristic over-approximates the exact solution by 25%. Furthermore, Fig. 17 also shows that the exact solution is 30% above the lower bound. Hence, by extrapolating this result, our heuristic would over-approximate the exact solutions of Fig. 16 by 35%. Finally, we evaluate the heuristic using five real applications: the H.263 decoder, the data modem and sample rate converter from the SDF 3 benchmarks [START_REF] Stuijk | SDF 3 : SDF for free[END_REF], the fast fourier transformer (FFT), and the time delay equalizer (TDE) from the StreamIt benchmarks [START_REF] Thies | An empirical characterization of stream programs and its implications for languages and compiler design[END_REF]. All these graphs have a chain structure. Table I shows some characteristics of these applications together with the obtained results. Our approach improves better the upper bounds in case of chains with a form (H.263 decoder and FFT). It comes close to the upper bound for the sample rate converter since the two actors with the highest loads are the right and left ends of the chain; increasing the loads of the other actors to get a monotone chain results in a size of almost 2(p + q -gcd(p, q)) for every channel. Load: Few symbolic results about SDF graphs can be found in the literature. In this section, we present the most relevant ones. Consistency can easily be checked analytically. The repetition vector can be computed symbolically as is it done in most dynamic parametric SDF models (e.g., [START_REF] Bebelis | BPDF: a statically analyzable dataflow model with integer and boolean parameters[END_REF], [START_REF] Fradet | SPDF: a schedulable parametric data-flow MoC[END_REF]).

There is no exact analytical solution to check the liveness of a graph with buffers with fixed bounds. In [START_REF] Bebelis | BPDF: a statically analyzable dataflow model with integer and boolean parameters[END_REF], [START_REF] Bempelis | Boolean Parametric Data Flow: Modeling -Analysis -Implementation[END_REF], the authors apply Eq. ( 1) transitively (which leads to nested ceilings) on edges of each cycle in the graph. Then, the obtained equations are linearized by over-approximating the ceiling function (i.e., x < x + 1). Yet, this is a conservative liveness analysis. As proved in [START_REF] Battacharyya | Software synthesis from dataflow graphs[END_REF], the minimum buffer size for which the simple graph A p q --→ B is live is equal to p+q -gcd(p, q) 4 . This however does not imply that any graph whose channels are sized this way is live. Still, this analytical equation is used in many buffer sizing algorithms to compute a lower bound as a starting solution [START_REF] Bempelis | Boolean Parametric Data Flow: Modeling -Analysis -Implementation[END_REF], [START_REF] Stuijk | Exploring trade-offs in buffer requirements and throughput constraints for synchronous dataflow graphs[END_REF].

Let s i denotes the token timestamp vector, where each entry corresponds to the production time of tokens in the i th iteration of the graph. Then, the max-plus algebra can be used to express the evolution of the token timestamp vector: s i = M s i-1 . It has been proved that the eigenvalue of matrix M is equal to the period of the graph. In case of parametric rates, it is sometimes possible to extract a max-plus characterization of the graph with a parametric matrix [START_REF] Skelin | Worst-cas throughput analysis for parametric rate and parametric actor execution time scenario-aware dataflow graphs[END_REF]. However, this works only in cases where Eq. ( 1) can be somehow simplified to get rid of the ceiling function (e.g., when p=1).

[10] presents a parametric throughput analysis for SDF graphs with bounded parametric execution times of actors but constant rates. Since rates and delays are non-parametric, the SDF-to-HSDF transformation is possible and the throughput analysis is based on the MCM. Therefore, all cycle means are linear functions in terms of parametric execution times. By using these linear functions, the parameter space is thus divided into a set of convex polyhedra called "throughput regions", each with a throughput expression. This approach has been extended in [START_REF] Damavandpeyma | Parametric throughput analysis of scenario-aware dataflow graphs[END_REF] to the case of scenario-aware dataflow (SADF) graphs.

A different analytic approach to estimate lower bounds of the maximum throughput is to compute strictly periodic schedules instead of ASAP schedules (e.g., [START_REF] Bodin | Periodic schedules for cyclo-static dataflow[END_REF]). This approach is similar to our Stretch linearization method. The advantage of the strictly periodic scheduling approach is its capability to handle cyclic graphs. However, not all cyclic graphs have strictly periodic schedules (it depends on the number of initial tokens). Furthermore, experiments on real-life benchmarks show that these approaches result in huge over-approximations (sometimes 7 times the exact value) [START_REF] Bodin | Periodic schedules for cyclo-static dataflow[END_REF]. In theory, the overapproximation is not even bounded.

VIII. CONCLUSION

We have studied analytically the different cases of the execution of a completely parametric single edge dataflow graph. Then, we have presented the exact symbolic solutions for the minimum buffer size needed by a single edge to achieve its maximal throughput. Using these results and a linearization technique, we have provided safe upper bounds of buffer sizes of acyclic graphs for maximal throughput. Furthermore, we 4 The equation is slightly different when there are initial tokens.

have proposed a heuristic to improve these bounds for graphs with a chain or a tree structure. Experimental results show that our heuristic improves the upper bound by 20% in average and can give the optimal solution for some real applications. We are following the same approach for exact and approximate symbolic evaluations of the latency of parametric graphs. Future work will concern the extension of these analysis to deal with general (i.e., possibly cyclic) dataflow graphs.

Fig. 1 :

 1 Fig. 1: Symbolic and numerical analyses.

  (b) has 15 initial tokens (represented by the black dot).

Fig. 2 :

 2 Fig. 2: (a) An SDF graph; (b) The same graph with channel size constraint and auto-concurrency disabled.

3 Fig. 3 :

 33 Fig.3: SDF-to-HSDF transformation of a graph and its dual.

Fig. 4 :

 4 Fig. 4: Illustration of the duality theorem.

r-- 1 .Fig. 5 :

 15 Fig. 5: Enabling patterns.

3 6 -

 36 -→ B is A 2 B ; i.e., after every two firings of actor A, one firing of B is enabled. The enabling pattern of A 8 5

Fig. 6 :

 6 Fig. 6: An ASAP execution (p = 8, q = 5, t A = 20, t B = 7).

Property 2 .

 2 Fig. 5 gathers all possible enabling patterns of a channel A p q --→ B.

  e., qt A ≥ pt B ): Actor A has the highest load and should fire consecutively for maximal throughput.

  2. (B.2 ∧ r ≤ 1 p r t A ) ∨ (B.3 ∧ r ≤ r p-r r p-r +1 t A ) where r = t B -kt A θ A,B = p + 2q -gcd(p, q) + r t A -r (p -r) (13) Case III. Case III.1. A.2 θ A,B = 2p+q+r-gcd(p, q)+ n-1 max j=1

  nr t B -r ≥ nr q-r and r = t A -kt B . Cases III.(A.3), III.(B.2), III.(B.3) see [6].

Fig. 7 :

 7 Fig. 7: Minimum buffer size θ A,B when z A t A ≥ z B t B .

prologue block 1 Fig. 8 :

 18 Fig. 8: Sequence (δ j ) in case (I.1) and pattern (A.2).

Fig.

  Fig. An ASAP execution (p = 8, q = 5, t A = 23, t B = 14).

  Consider the graph G = {A p q

Fig. 11 :

 11 Fig. 11: Linearization of B for the graph A p q --→ B (p = 8, q = 5, tA = 20, tB = 7). Case qt A ≥ pt B .

  C} with t A = 16, t B = 11 and t C = 12. The repetition vector is z = [4, 3, 6]. Actor C imposes the higher load, hence the minimal period of this graph is P G b = z C t C = 72. We have θ A,B = 9 and θ B,C = 6.

Fig. 12 :

 12 Fig. 12: ASAP execution of G b Property 5. Let G be a graph without any undirected cycle, if the buffer of every channel A p q

Fig. 14 :

 14 Fig. 14: ASAP execution of G d .

Fig. 16 :

 16 Fig. 16: Heuristic vs. upper bounds.

Fig. 17 :

 17 Fig. 17: Heuristic vs. exact solution.

TABLE I :

 I Experimental results for real benchmarks.

	graph	# actors	A z A	Upper bound	Optimal size	Heuristic
	(a) modem	6	37	32	20	31
	(b) sample con.	6	612	60	34	57
	(c) H.263 dec.	4	1190	2378	1257	1257
	(d) FFT	11	94	992	504	808
	(e) TDE	27	2867	7328	3680	5384

Mathematically speaking, the results of this section can be generalized to rates defined as positive real numbers.


B Colombet

M Woodman

C G Bénar

J M Badier

AnyWave: a cross-platform and modular software for visualizing and processing electrophysiological signals

Keywords: multi-platform software, signal processing, electrophysiology, MATLAB, Python, EEG, MEG

BACKGROUND

The importance of digital signal processing in clinical neurophysiology is growing steadily, involving clinical researchers and methodologists. There is a need for crossing the gap between these communities by providing efficient delivery of newly designed algorithms to end users.

We have developed such a tool which both visualizes and processes data and, additionally, acts as a software development platform.

NEW METHOD

AnyWave was designed to run on all common operating systems. It provides access to a variety of data formats and it employs high fidelity visualization techniques. It also allows using external tools as plug-ins, which can be developed in languages including C++, MATLAB and Python.

RESULTS

In the current version, plug-ins allow computation of connectivity graphs (non-linear correlation h2) and time-frequency representation (Morlet wavelets). The software is freely available under the LGPL3 license.

COMPARISON WITH EXISTING METHODS

AnyWave is designed as an open, highly extensible solution, with an architecture that permits rapid delivery of new techniques to end users.

CONCLUSIONS

We have developed AnyWave software as an efficient neurophysiological data visualizer able to integrate state of the art techniques. AnyWave offers an interface well suited to the needs of clinical research and an architecture designed for integrating new tools.

We expect this software to strengthen the collaboration between clinical neurophysiologists and researchers in biomedical engineering and signal processing.

Introduction: the need for a platform linking clinicians and methodologists

Signal processing techniques play an increasingly important role in clinical neurophysiology. Two communities are involved in the use of these techniques. On the one hand, researchers in clinical neurophysiology have adopted new techniques for working with EEG since the development of personal computer and signal digitization, allowing filtering, re-montaging without change of the original signal, mapping and source localization [START_REF] Rosenblatt | Computerized EEG monitoring[END_REF][START_REF] Wendling | Clinical engineering in neurophysiology: an object-oriented platform for signal visualization and processing[END_REF]. On the other hand, researchers in signal processing have been very active to solve the challenging mathematical problems posed by extracting information from brain signals [START_REF] Baillet | Electromagnetic brain mapping[END_REF]. However, there are currently several practical barriers prevent efficient exchange between the two communities. First, clinical research software is typically closed source and cannot be extended by third parties. Second, researchers deal with different data formats and often redefine basic preprocessing and visualization functions. Finally, new signal processing tools are difficult to distribute to clinical users as they often are not easy to use. Several available software packages have been developed by methodologists as MATLAB (Mathworks, Natick, MA) toolboxes [START_REF] Delorme | EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF][START_REF] Litvak | EEG and MEG data analysis in SPM8[END_REF][START_REF] Oostenveld | FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data[END_REF] or stand-alone software [START_REF] Aguera | ELAN: a software package for analysis and visualization of MEG, EEG, and LFP signals[END_REF][START_REF] Brunet | Spatiotemporal analysis of multichannel EEG: CARTOOL[END_REF]. They are usually open source, and the code can be adapted to the users' needs. However, none of them offers an efficient, easy to use interface that is appropriate for clinical work, where visualization of signals and marking events is important. Therefore, there is a strong need for a common platform which would allow rapid interaction between clinical researchers in the two communities. Our goal was to develop a tool at the interface between methodologists and researchers in clinical neurophysiology. Our first objective was to design a cross-platform software which provides basic visualizations, marking tools and data export functionality. Our second objective was to enable rapid incorporation of signal processing modules in order to facilitate the use of newly developed algorithms by methodologists. We present here the architecture of the AnyWave software, with an emphasis on the signal processing plug-ins, which can be developed in C++, MATLAB or Python. Our software is available for free under the LGPL3 license, along with full documentation and examples at http://meg.univ-amu.fr/wiki/AnyWave . The source code can be obtained at https://github.com/anywave/anywave

Architecture and general principles

In order to combine performance and portability we have chosen the Qt framework which is available under LGPL licensing and uses C++ as a programming language. The VTK library, also in C++ and under the LGPL license, was also used for 3D and 2D data representation, such as mappings of sensors activities for EEG or MEG recordings. A software development kit is also provided for developing plugins from MATLAB and Python scripts and incorporating them into AnyWave. This process will be described in section 4.

The VTK library was chosen because it provides algorithms (Delaunay2D, Delaunay3D, etc.) which are ready to use as well as features implementing 3D renderings. The library is only used to render 2D or 3D mapping of activities. There is no requirement to use the VTK library for developers who would like to implement their own plug-ins.

In order to save time between different reviewing sessions, filtering options, montage and visualization settings are automatically saved along with the data file. Hence, when the user opens the data file again; all previous settings are automatically applied.

Software components

AnyWave is composed of a set of components; some of them are running in different threads of execution (see Figure 1). This allows for a responsive graphic user interface during heavy calculation or when loading a large amount of data. For example, the Data Server component can be launched several times, each time creating a new thread. In AnyWave, data is processed using a client-server mechanism. The Data Server component is the server. Others components or modules can become clients: to do so a connection must be established with the server. Opening a connection creates a new Data Server component in a new thread, providing the client with a dedicated resource from which to request data. The data flow is presented in Figure 2. The data server thread is in charge of loading the requested part of data into memory, applying montage (through the Montage Manager Component), filtering data if necessary and then sending data.

The main part of AnyWave is the Visualization component which is the graphic user interface. This component is a client of Data Server component; therefore the main interface will request data to a server running in the background. When data is loaded from a file, each channel is handled based on its name. Then, the Montage Manager can apply user defined montages, e.g., bipolar montages or re-referenced data, and change the order of channels. User-defined montages can be saved to a file in order to apply that montage to another file subsequently. A data file may also contain events.

Events are called markers in AnyWave and are managed by the Markers Manager component, using a dedicated user interface.

The Process Manager component handles optional signal processing algorithm modules (see next section). It is able to launch, pause or cancel signal processing modules as well as handle their input and output data.

Loading and managing optional plug-ins is the role of the Plug-ins Manager component. Signal processing plug-ins also become clients of Data Server component and then request data without being interrupted by other modules that may request data at the same time.

Plug-ins

One goal in the design of AnyWave was modularity. To achieve that, AnyWave was built with the concept of plug-in that will add features to the software. A plug-in is a C++ module built using the provided software development kit (SDK). There are several types of plug-ins that could be built using the SDK: "reader", "writer", "display" or "signal processing algorithms". Signal processing algorithms can be written in MATLAB or Python, and MATLAB plugins can be compiled to avoid requiring target computers to have MATLAB installed. We are not providing a complete and monolithic solution but rather a modular way of handling input/output and signal processing. With this design, it is possible to support new file formats without requiring a new version of AnyWave, thanks to a new reader or writer plug-in. It is also possible to adapt and include an existing signal processing algorithm by implementing it in a signal processing plug-in. This is particularly easy if the original algorithm is written in MATLAB or Python language.

Reader and writer plug-ins

A reader plug-in is a module that opens a data file and sends data (signals and markers) to AnyWave.

Currently, there are several plug-ins provided that can read several EEG or MEG file formats from different manufacturers: MicroMed (.trc), Neuroscan (.cnt), Brain Products (.vhdr), CTF (.meg4), 4D

NeuroImaging, as well as EDF/BDF format (.bdf, .eeg). Adding new formats only requires writing a new reader plug-in (in C++/Qt). Writer modules are used to export data for further use in other software, and currently exported formats include Brain Vision Analyser (Brain Products GmbH) (.vhdr) and AnyWave binary format (.ades). Note that the ADES format is not the official file format associated with AnyWave. It is simply a convenient binary format with a text header, allowing exchanging data between software and platforms. The documentation is available at http://meg.univ-amu.fr/wiki/AnyWave:ADES As in the case of reader plug-ins, the list of supported file formats available to export can be extended by adding new writer modules that are independent of the main program.

Signal processing and Display plug-ins

Signal processing plug-ins get data from AnyWave and can return computed data directly to AnyWave (see section 4). Each plug-in corresponds to a process.

There are two main kinds of signal processing plug-ins: (i) Those which will run in the background, independently of the visualization components. (ii) Those which will be coordinated with the visualization component.

These processes request data, run an algorithm on the data they receive, and then terminate. They can optionally return some data to the visualization component.

The difference between the two types of signal processing plug-ins (or processes) is that Background launched processes will run once when requested by the user and then terminate. Display processes are linked to what the visualization component is currently displaying. Whenever the displayed content changes, those processes are executed. For example, the time-frequency signal processing plug-in computes the time-frequency representation for one or more channels displayed by the visualization component and sends back their respective time/frequency representations. As it is designed as a signal processing linked to the display, each time the content of the channels is updated (while the user is browsing the data file) the process is run and the time/frequency representations are updated. A Signal processing plug-in (or process) may require data from the Data Server component. This is the case for processes that are designed to run in the background.

Processes designed to be linked to the current display, will not need to request for data as the data is automatically provided by the visualization component. A process can access any piece of data if needed (based on time or markers for example). The data received by the client have been processed through the montage manager and some filtering may occur also, depending on the settings the user set. However, the process can request the Data Server component to provide unfiltered data or ask for a specific filtering to be applied on a particular channel type. Some algorithms can produce data which are not vectors of values, which is the default data type that the visualization component can display. This is the case for example, for time-frequency analyses. This is where Display plug-ins are used to extend the capabilities of the visualization component and allow displaying other data types, like 2D matrices.

Visualization

Signal visualization

AnyWave displays signals in views. By default one main view is used to display the most common channel types (EEG, MEG, SEEG). The user can choose to add other views and can define options to handle the desired type of channels to display and how they are displayed. This permits the visualization of different type of data or different time scales in separate views (see Figure 3). Special care has been dedicated to mapping the signal, which can be acquired at high sampling rate, to pixels on the screen. When a large amount of data is displayed on the screen, the horizontal resolution of a computer screen effectively bins together adjacent time points, and simply downsampling the data for sake of efficiency could lead to aliasing. Instead, for each horizontal bin, we have taken into account the minimum and maximum signal within this bin in order to correctly rasterize the signal. This approach retains high frequency information such as high frequency oscillations. Figure 7 illustrates the rendering of an oscillation at two different time scales.

When the horizontal scaling is changed, signals are displayed again using this technique.

The visualization component is connected to a Data Server thread, as explained in section 2.1.

In this scheme, handling large amount of data is done by only requesting the data part that needs to be displayed. We found that on a computer equipped with an Intel Core I7 at 2.8 GHz, the rendering of 60 channels of data on a 20 s page, sampled at 2034Hz, takes less than 250ms.

Filtering and scaling options

AnyWave can filter data differently depending on the channel type: EEG channels may have different filtering settings than MEG channels. The same principle applies to gain settings, although a global gain factor can be applied to all channels at once. The horizontal scale can be modified by setting the number of centimeters by seconds.

The filtering procedure for each channels, is dispatched on all the available CPU cores. This is accomplished using the QtConcurrent API from the Qt framework. Therefore, CPUs with 4 cores available will filter data from 4 channels at the same time.

Maps of sensors activities

AnyWave can generate 3D or 2D representations of the sensors and compute the mapping of their activities at a given time. Sensors coordinates must be provided in order to compute the mesh representation or the 2D projection (Figure 3). Coordinates can be read from the data file or extracted from an EEG template which is embedded in AnyWave. It is also possible to import coordinates of the sensors from a text file.

Markers

AnyWave uses markers to note particular events in signals (see Figure 4). Markers can also be read from the data file or imported afterwards. The user can also add her/his own markers manually or by using a signal processing algorithm that generates markers. Markers can consist of either selected time points or sections of data, with the option of specifying a subset of channels. AnyWave offers a dedicated user interface for managing markers, which allows the user to edit existing markers, launch signal processing on selected markers and navigate in signals through markers. When dealing with many markers, it is possible to define filtering rules that allow the user to quickly switch between different sets of markers based on their labels, integer values, and other metadata.

The software can handle a large amount of markers. The user can add markers in a convenient way, using keyboard shortcuts to navigate through existing markers for example. It is also possible to prepare a list of possible markers to add depending on the preceding event or marker. Then, when adding a marker, the user can quickly choose the type of marker he/she wants to add.

Displaying markers using filtering rules offers an ergonomic way to present only relevant markers on signal views. Markers can also be saved to a text file for processing in another software tool (MATLAB for example).

Montages

Montages are well known by clinicians who are used to review EEG or SEEG data on a daily basis.

They describe which sensors or signals will be used from the data file and how they will be presented to the user. SEEG data are often visualized using bipolar montages, which can be implemented within AnyWave using the dedicated Montage user interface. A montage is involved in both the visualization component and the data server component: the visualization component applies the defined montage to inform the signal views on what sensors should be displayed on the screen. The data server component makes use of the montage while getting data from the file, before sending them to clients. When opening a data file for the first time, a default montage is automatically defined. This montage contains all the sensors/channels as found originally in the data file. However, if the data has already been open, the last montage used with the file is automatically loaded and set.

Developing signal processing plug-ins

The most basic way of adding a plug-in is to implement the necessary C++ classes (see 2.2.2).

However, MATLAB and Python are common tools for implementing signal processing algorithms;

AnyWave offers a way to run processes written in these languages. To do so, a dedicated software development kit is provided with tutorials (available at http://meg.univ-amu.fr/wiki/AnyWave).

MATLAB and Python scripts

If a full copy of MATLAB is available on the system, AnyWave can run arbitrary MATLAB scripts (using the MATLAB Engine API) as plug-ins with full access to the same data as a C++ plug-in. For Python, AnyWave comes with many scientific libraries that user scripts can use, or if desired, the user can install a distribution that includes the suitable scientific libraries. Both MATLAB and Python scripts can request and process data, and optionally return data or results. On startup, AnyWave will look for plug-ins in the default directory or in the user's home folders and automatically add them to the Processes menu. The only requirement to create a Python or a MATLAB plug-in is to create a folder and place in it at least two required files: the MATLAB main function (main.m) or the Python (__main.py__) main script. The second file is a descriptive text file that will inform AnyWave about the plug-in's properties (name, description, etc.). The format of the descriptive text file is exactly the same for MATLAB and Python plug-ins. Once the two files are completed, the folder containing it must be placed in a location where AnyWave will notice it after a restart.

The scripting method of developing and providing plug-ins is useful because it does not require compiling code and the script that defines the plug-in can be changed on the fly to provide immediate updates to the algorithm, allowing the developer to iterate quickly. MATLAB and Python are multiplatform languages, so developers can easily add and share their own code without targeting a particular platform. Python plug-ins offer the possibility to make a bridge between existing Python solutions as Openelectrophy [START_REF] Garcia | OpenElectrophy: An Electrophysiological Data-and Analysis-Sharing Framework[END_REF] or MNE-Python [START_REF] Gramfort | MEG and EEG data analysis with MNE-Python[END_REF] on the one hand and AnyWave on the other hand. An illustration is shown in the code 2 example below. code 1 example: Code 1: Example of MATLAB Plug-in. This script demonstrates the use of MATLAB code as a plug-in.

The script first obtains 30 seconds of data from AnyWave, with the current filter setting, detects for each channel where the signal exceeds its 99th percentile, and creates markers based on this detection. The results of the execution of this script can be seen on figure 5. code 2 example: import numpy import mne from mne.io.array import create_info, RawArray # get all data & create MNE data structure # using anywave.getData() method: # first parameter is the starting position in seconds. # second parameter is the duration of data in seconds. # Specifying -1 as duration indicates ALL the data available. anywave.getData(0, -1) # getData() will fill an attribute if anywave module called input_channels # input_channels is an array of AwChannel objects.

MATLAB Compiled plug-ins

Python is available freely (on most common platforms), which is not the case for MATLAB.

To run plug-ins without the MATLAB license, AnyWave supports plug-ins which embeds MATLAB compiled code. To build this type of plug-in, both the MATLAB Compiler and a C++ compiler are required. Fortunately, the MATLAB code can be easily ported from a MATLAB script as the MATLAB functions used to communicate with AnyWave have the same names and parameters. The final user will need only the MATLAB Compiler Runtime in order to run the plug-in (freely available at http://www.mathworks.fr/products/compiler/mcr/). Developers have the choice of providing compiled plug-ins or Matlab scripted plug-ins. We encourage this second option. Our github (https://github.com/anywave) can host dedicated branches for new MATLAB scripted plug-ins.

To distribute MATLAB scripted plug-ins to end users who are not using the MATLAB software, we propose to build a compiled version for them, using the AnyWave compatible MCR version.

C++ plug-ins

C++ is the language used to build AnyWave and all its components. Therefore, it is possible to develop all types of plug-ins using C++ where MATLAB and Python can only be used to build Signal Processing Algorithms (or processes). C++ plug-ins required to be compiled for specific hardware architectures, so a Mac OS X plug-in will not run on a Linux platform, for example. However, the C++ API offers, for the moment, more features while developing plug-ins. The detailed documentation on this purpose can be found at http://meg.univ-amu.fr/wiki/AnyWave A possible link between C++ sofwares, such as OpenVibe [START_REF] Renard | OpenViBE: An Open-Source Software Platform to Design, Test and Use Brain-Computer Interfaces in Real and Virtual Environments[END_REF], and AnyWave can be established. For example, developers could use OpenVibe C++ API in their implementation of a C++ signal processing plug-in.

As an example of a C++ signal processing plug-in, Figure 6 shows a time-frequency representation, which uses two plug-ins of different types. The first plug-in is a signal processing algorithm that takes as input one or more signals and computes their time-frequency representation, based on Morlet wavelets [START_REF] Bénar | Consensus Matching Pursuit for multi-trial EEG signals[END_REF]. For each input signal, a new signal is created as output and configured as "custom" type. As this signal does not contain a vector of value but a two dimensions matrix, the visualization component will not be able to natively display it. This is where the second plug-in, a Display plug-in, is used to render the 2D matrix in the visualization component. Note that the signal processing module is linked to display: when the user moves to another part of the file, the corresponding time-frequency representation is automatically re-computed and display is updated. This is the default behavior for signal processing plug-in of type Display.

Concluding remarks and future directions

AnyWave responds to a consistent need to deliver state of the art signal processing and visualization techniques in the context of the clinic and wherever an efficient, user friendly interface is required.

As a software toolbox and platform, it is not alone. In this context, academic free software have been rapidly developing. A few of those are provided as stand-alone programs, for example MNE (http://martinos.org/mne/), Brainstorm (http://neuroimage.usc.edu/brainstorm/), openMEEG [START_REF] Gramfort | OpenMEEG: opensource software for quasistatic bioelectromagnetics[END_REF], Cartool [START_REF] Brunet | Spatiotemporal analysis of multichannel EEG: CARTOOL[END_REF] or Python programs [START_REF] Garcia | OpenElectrophy: An Electrophysiological Data-and Analysis-Sharing Framework[END_REF] but the most popular ones are MATLAB toolboxes [START_REF] Delorme | EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF][START_REF] Oostenveld | FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data[END_REF]. In each case, the scopes of these toolboxes are oriented by original objectives of the developers. Brainstorm is oriented to source localization and presents a userfriendly windows interface although it can be batch scripted [START_REF] Tadel | Brainstorm: a user-friendly application for MEG/EEG analysis[END_REF]. Fieldtrip is a collection of scriptable functions from preprocessing to statistics but requires a certain level of programming to use effectively [START_REF] Oostenveld | FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data[END_REF]. EEGLAB is well designed for processing of continuous signals and evoked activities [START_REF] Delorme | EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF] but its capacities in signal representation and handling markers are limited. Additionally, it can be easily adapted to new requirements. In particular, while algorithms can be written in C++ for full speed, we can also take advantage of the plethora of tools already available in MATLAB and Python toolboxes. Furthermore, the MATLAB and Python frameworks allow a plug-in to save data in other file formats for further processing with other tools if necessary.

We do not want to propose a replacement for existing software but rather a complementary solution with good visualization of signals and a user friendly interface. With AnyWave, it is now possible for neurophysiologists to visualize and process data through a user-friendly interface and for researchers in biomedical engineering and signal processing can easily test and put into practice new algorithms, making the collaboration between both worlds more efficient. To our knowledge, only one development [START_REF] Wendling | Clinical engineering in neurophysiology: an object-oriented platform for signal visualization and processing[END_REF] has had a similar goal. It was, however, developed on NextStep/OpenStep operating system, a platform that no longer exists and didn't allow for the use of multiple programming languages.

In summary, AnyWave is designed as a dynamic solution: plug-ins will add new features, new signal processing algorithms, new supported file format and so on. The software currently supports working with continuous data set. Support for epoched data and a dedicated user interface to manage averaged files and evoked potentials is currently under development. The architecture allows future extensions to support additional types of information, such as MRIs and models for source localization.

FIGURE LEGENDS

#

  create numpy array from input_channels data = numpy.array([ch.data for ch in anywave.input_channels]) # Get the sampling rate for data by reading sampling_rate attribute # of the fist AwChannel object in input_channels. freq = anywave.input_channels[0].sampling_rate # Get the electrode labels labels = [ch.label for ch in anywave.input_channels] # create the MNE data structure types = ['meg'] * len(labels) info = create_info(labels, freq, types) raw = RawArray(data, info) # run ICA & save results ica = mne.preprocessing.ICA(n_components = 50) ica.fit(raw, decim=3) # Saving ica result along with data file using the data file name ica.save(anywave.data_path + '-ica.fif') Code 2: This script illustrates how to use MNE-Python to compute ICA on data coming from AnyWave. AnyWave specific methods or attributes are available through an automatically imported module named anywave. Note that depending on MNE-Python version, the line importing create_info can fail. Try importing from mne.io.meas_info instead.

Figure 1 :

 1 Figure 1: Software components, organized in layers from the operating system to the dynamically loaded plug-ins.

Figure 2 :

 2 Figure 2: Data flow across software components (rectangles) and plug-ins (diamond shapes).

Figure 3 :

 3 Figure 3: Overview of the main window (the visualization component), presenting different type of traces of the same recording and a 3D mapping of sensors activities. Different filters have been applied for the different types of data.

Figure 4 :

 4 Figure 4: Overview of the main window (the visualization component), showing markers. To the left, the markers' user interface is docked and allows quick navigation through markers. Markers are visible on the top view, with their labels and associated values.

Figure 5 :

 5 Figure 5: Signals with markers, added by a scripted plug-in.

Figure 6 :

 6 Figure 6: An example of Time-Frequency representation of a signal.

Figure 7 :

 7 Figure 7: Illustration of the visual rendering of a high frequency oscillation (red arrow) at two different time scales.

Figure 1 .

 1 Figure 1. Software components, organized in layers from the operating system to the dynamically loaded plug-ins.

Figure 2 .

 2 Figure 2. Data flow across software components (rectangles) and plug-ins (diamond shapes).

Figure 3 .

 3 Figure 3. Overview of the main window (the visualization component), presenting different types of traces of the same recording and a 3D mapping of sensors activities. Different filters have been applied for the different types of data.

Figure 4 .

 4 Figure 4. Illustration of the visual rendering of a high frequency oscillation (red arrow) at two different time scales.

Figure 5 .

 5 Figure 5. Overview of the main window (the visualization component), showing markers. To the left, the markers' user interface is docked and allows quick navigation through markers. Markers are visible on the top view, with their labels and associated values.

Figure 6 .

 6 Figure 6. Signals with markers, added by a scripted plug-in.

Figure 7 .

 7 Figure 7. An example of time-frequency representation of a signal.