Francesco Komauli
email: francesco.komauli@studenti.unimi.itmomigliano@di.unimi.it

Alberto Momigliano

Property-Based Testing of Abstract Machines: an Experience Report

Contrary to Dijkstra's diktat, testing, and more in general, validation, has found an increasing niche in formal verification, prior or even in alternative to theorem proving. Validation, in particular, propertybased testing (PBT) is quite effective in mechanized meta-theory of programming languages, where theorems have shallow but tedious proofs that may go wrong for fairly banal mistakes. In this report, we abandon the comfort of high-level object languages and address the validation of abstract machines and typed assembly languages. We concentrate on Appel et al.'s list-machine benchmark [2], which we tackle both with αCheck, the simple model-checker on top of the nominal logic programming αProlog and the PBT library FSCheck for F#. This allows us to compare the relative merits of exhaustive-based PBT in a logic programming style versus the more usual randomized functional setting. We uncover one major bug in the published version of the paper plus several typos and ambiguities thereof. This is particularly striking, as the paper is accompanied by two full formalizations, in Coq and Twelf. Finally, we do a bit of mutation testing on the given model, to asses further the trade-off between exhaustive and randomized data generation. Spoiler alert: the former performs better.

Introduction

Does this sound familiar? You are in the middle of a long but supposedly straightforward formal proof trying to drag your favourite proof assistant to confirm the blindingly obvious truth of the result, when you get stuck in an unprovable part of the derivation; you realize that the statement needs to be adjusted or more commonly that there is something afoul, typically a banality, in your specification. If only you had a way to realize that the theorem was unprovable before wasting precious time in a doomed proof attempt, perhaps some kind of testing . . . This is where, with due respect to Dijkstra and his "thou shall not test" commandment, validation prior to theorem proving comes in: as a matter of fact, such tools have been available for some time in the leading proof assistants: to cite only two of them, the NitPick/QuickCheck combination [START_REF] Christian Blanchette | Automatic Proof and Disproof in Isabelle/HOL[END_REF] (and descendants) in Isabelle/HOL and QuickChick [START_REF] Paraskevopoulou | Foundational Property-Based Testing[END_REF] in Coq.

Our angle is not verification vs. validation in general, but in a particular domain: the mechanization of the meta-theory (MMT) of programming languages (PL) and related calculi. As any practitioner may testify, the main properties of interest are well known and have mathematically shallow proofs. The difficulties lies in the potential magnitude of the cases one must consider and in the trickiness that some encodings require, typically those having to deal with binding signatures. Here, very minor mistakes in the specification, even at the level of what we would consider a typo, may severely frustrate the verification effort, to the point to make it not cost-effective. Further, we aim to support the "working semantics" in her work in developing, and eventually formally proving correct, such calculi, rather than concentrating on models that we already know to be correct.

The model checking technique we have adopted is property-based testing (PBT): while it originated as a testing techniques for concrete (functional) programs, it is also applicable at the meta-theoretical level, especially in the PL domain, as shown by the success [START_REF] Casey | Run your research: on the effectiveness of lightweight mechanization[END_REF] of tools such as PLT-Redex (https:// redex.racket-lang.org) and, "si parva licet componere magnis", αCheck [START_REF] Cheney | αCheck: A mechanized metatheory model checker[END_REF]. PBT's data generation strategy comes in several flavours: random [START_REF] Claessen | QuickCheck: a lightweight tool for random testing of Haskell programs[END_REF], exhaustive [START_REF] Runciman | Smallcheck and lazy smallcheck: automatic exhaustive testing for small values[END_REF], a combination of the two [START_REF] Duregård | Feat: functional enumeration of algebraic types[END_REF], with an ever increasing emphasis on coroutining the generation and testing phases [START_REF] Fetscher | Making Random Judgments: Automatically Generating Well-Typed Terms from the Definition of a Type-System[END_REF][START_REF] Lampropoulos | Generating good generators for inductive relations[END_REF]. An analysis of some of those strategies as applied to MMT was carried out in [START_REF] Fachini | Validating the Meta-Theory of Programming Languages (Short Paper)[END_REF].

In this paper, we abandon the comfort of high-level object languages, which have been investigated extensively, and address the validation of abstract machines and typed assembly languages. That seems more challenging, since instructions operating at a low level provide less "structure" for counterexamples, which then tend to be substantially more complex. Further, it is hard to generate sequences of machine states that yield meaningful machine runs. Similar remarks have appeared in the context of random generation of C programs [START_REF] Yang | Finding and understanding bugs in C compilers[END_REF]. Not coincidentally, the authors of [START_REF] Hritcu | Testing noninterference, quickly[END_REF] suggest that "[counterexamples to properties such as non-interference may be] well beyond the scope of naive exhaustive testing" (pag. 12).

We take on Appel et al.'s list-machine benchmark [START_REF] Appel | A List-Machine Benchmark for Mechanized Metatheory[END_REF], dubbed by the authors CIVmark, standing for Compiler Implementation Verification. We concentrate mostly on its version 1.0, as originally presented at LFMTP'06, but we also touch on 2.0, see Section 6.1. CIVmark was conceived as a benchmark for "machine-checked proofs about real compilers", spurred by a civilized criticism of POPLMark as being too biased towards issues of binders. It consists of a basic pointer machines with instructions such as car, cdr, cons and (un)conditional jump, endowed with a standard SOS, but with a reasonably sophisticated type system: the latter, being able to make static predictions about lists being empty or not, guarantees that a well typed machine's run does not get stuck. Two implementations, including sketch of proofs of type soundness in Twelf and Coq come with the paper.

A version of this model is also part of the benchmarks distribution of PLT-Redex [START_REF] Findler | The PLT-Redex list-machine model[END_REF] and as such being subjected to some PBT. All the more reasons to tackle with two very distinct PBT's approaches, which fits well with the existing formalizations:

1. An encoding in the nominal logic programming αProlog [START_REF] Cheney | Nominal logic programming[END_REF], to be tested with its model checker αCheck [START_REF] Cheney | αCheck: A mechanized metatheory model checker[END_REF]; 2. A functional encoding with the PBT library FSCheck for F#. This allows us to compare the relative merits of exhaustivity-based PBT in a logic programming style versus the more usual randomized functional setting polarized by QuickCheck. It is also a stress test for αCheck, since CIVmark does not exploit any of the features offered by nominal logic that makes αProlog effective, see the case studies at https://github.com/aprolog-lang/checker-examples, which are all devoted to high-level languages, including some studied in [START_REF] Casey | Run your research: on the effectiveness of lightweight mechanization[END_REF]; instead, the CIVmark benchmark brings to the fore the naivete of αProlog's search strategy.

A testing approach is "good" only if it uncovers bugs, and so we did: we falsified the type preservation property as presented in the paper, caused by an incorrect specification of typing for values. This is particularly striking, considering the paper comes with two formalized proof of type soundness. The mystery disappears once we realized that the Coq implementation of that judgment was different from the paper and coincided with the definition that we had reverse-engineered from the counter-example to preservation 1 . We also found several typos and ambiguities in the typing rules in the published paper, but this is not unusual where there is no connection between the text in the paper and the model verified by the proof assistant.

That was encouraging, but to asses further the trade-off between exhaustive and randomized data generation, we resorted to some mutation testing [START_REF] Pezzè | Software testing and analysis -process, principles and techniques[END_REF] on the given debugged model. As we will see in Section 5, exhaustive generation, in all its naivety, tends to be more cost effective than pure random PBT. Finally, we report on some ongoing work on the list-machine 2.0 and on replaying the first section of [START_REF] Hritcu | Testing noninterference, quickly[END_REF], which deals with an abstract machine for dynamic information-flow control. Here again, exhaustive PBT shows its colours.

The list machine

We present here the syntax, static and dynamic semantics and finally meta-theory of CIVmark [START_REF] Appel | A List-Machine Benchmark for Mechanized Metatheory[END_REF].

The plumbing of the machine

The list-machine, as the name suggests, operates over an abstraction of lists, where every value is either nil or the cons of two values, in other terms S-expressions with only nil as atom.

value a ::= nil | cons(a 1 , a 2 )
Given a set of variables and labels, the machine features the following set of instructions:

ι 1 , ι 2 , . . . : I instructions jump l : I jump to label l branch-if-nil v l : I if v = nil then jump to l fetch-field v 0 v : I fetch the head of v into v fetch-field v 1 v : I fetch the tail of v into v cons v 0 v 1 v
: I make a cons cell in v halt : I stop executing ι 1 ; ι 2 : I sequential composition

We rely on several association lists: programs map labels to instructions and stores variables to values.

program p ::= end | p, l n : ι store r ::

= { } | r[v → a]
We use a functional notation such as p(l) for look-up into such structures. The notation r[v → a] assumes that the variable v is not in the domain of r and we use r[v := a] = r for functional update.

Dynamic and static semantics

The operational semantics is summarized in the inference rules in Fig. 1: given a fixed program, the small-step relation (r, ι) p → (r , ι ) works on store/instruction configurations in a continuation-passing style. Then, a big-step semantics is defined as the Kleene closure of the small-step relation, with the instruction halt that signals the end of a program execution. A program p is said to run if it runs in the big-step relation, starting from the instruction at p(l 0 ) with an initial store containing only the binding v 0 = nil, until a halt instruction is reached.

The type system assigns to each variable a list type that is then refined to empty and nonempty lists, to guarantee safety of certain operations, e.g. fetch-field.

type τ ::= nil | list τ | listcons τ (r, (ι 1 ; ι 2 ); ι 3 ) p → (r, ι 1 ; (ι 2 ; ι 3 )) step-seq r(v) = cons(a 0 , a 1 ) r[v := a 0 ] = r (r, (fetch-field v 0 v ; ι)) p → (r , ι) step-fetch-field-0 r(v) = cons(a 0 , a 1 ) r[v := a 1 ] = r (r, (fetch-field v 1 v ; ι)) p → (r , ι) step-fetch-field-1 r(v 0 ) = a 0 r(v 1 ) = a 1 r[v := cons(a 0 , a 1 )] = r (r, (cons v 0 v 1 v ; ι)) p → (r , ι) step-cons r(v) = cons(a 0 , a 1 ) (r, (branch-if-nil v l; ι)) p → (r, ι) step-branch-not-taken r(v) = nil p(l) = ι (r, (branch-if-nil v l; ι)) p → (r, ι ) step-branch-taken p(l) = ι (r, jump l) p → (r,
p → (r , ι ) (p, r , ι ) ⇓ (p, r, ι) ⇓ run-step (p, r, halt) ⇓ run-halt [v 0 = nil] = r p(l 0 ) = ι (p, r, ι) ⇓ p ⇓ run-prog
Figure 1: Small and big step operational semantics

The type system includes therefore the expected subtyping relation and a notion of least common supertype τ τ . A type environment Γ is a mapping between variables and types and subtyping is extended to environments width and depth-wise. A program typing Π is an association list of labeled environments, where Π(l) = Γ represents the types of the variables when entering a block labeled with l.

typing env Γ ::= { } | Γ, v : τ program typing Π ::= { } | Π, l : Γ
Type-checking is stratified in several judgments following the structure of a program as a labeled sequence of blocks. At the bottom, instruction typing Π instr Γ{ι}Γ transforms a precondition Γ into post-condition Γ under the program typing Π. The check-block relation handles terminal instructions halt and jump. Block typing is extended to all blocks of a program p with the check-blocks relation Π blocks p, whereby for each label l in program p, the associated block ι = p(l) is type-checked with its environment Γ = Π(l). Finally, the top-level type-checking relation |= prog p : Π states that a program p type-checks with a program typing Π when blocks and environments have matching labels, every block type-checks with its corresponding environment, and the initial environment is Γ 0 = (v 0 : nil).

Π instr Γ{ι 1 }Γ Π instr Γ {ι 2 }Γ Π instr Γ{ι 1 ; ι 2 }Γ check-instr-seq Γ(v) = list τ Π(l) = Γ 1 Γ[v := nil] = Γ Γ ⊂ Γ 1 Π instr Γ{branch-if-nil v l}(v : listcons τ, Γ ) check-instr-branch-list Γ(v) = listcons τ Π instr Γ{branch-if-nil v l}Γ check-instr-branch-listcons Γ(v) = nil Π(l) = Γ 1 Γ ⊂ Γ 1 Π instr Γ{branch-if-nil v l}Γ check-instr-branch-nil Γ(v) = listcons τ Γ[v := τ] = Γ Π instr Γ{fetch-field v 0 v }Γ check-instr-fetch-0 Γ(v) = listcons τ Γ[v := list τ] = Γ Π instr Γ{fetch-field v 1 v }Γ check-instr-fetch-1 Γ(v 0 ) = τ 0 Γ(v 1 ) = τ 1 (list τ 0 ) τ 1 = list τ Γ[v := listcons τ] = Γ Π instr Γ{cons v 0 v 1 v}Γ check-
Π instr Γ{ι 1 }Γ Π; Γ block ι 2 Π; Γ block ι 1 ; ι 2 check-block-seq Π(l) = Γ 1 Γ ⊂ Γ 1 Π; Γ block jump l check-block-jump

Properties

One of the few drawbacks of PBT as a testing technique for programming units is coming up with meaningful properties. In MMT, this is a non-issue, as PL calculi come equipped with the meta-theorems they must satisfy, when they are not variants of standard results. As far as the list-machine is concerned, we are spoiled with choices: the paper (pages 467-8) lists more than a dozen theorems ranging from basic properties of subtyping to type soundness. We report the top-level ones:

Progress: given a well-typed instruction and a well-typed store, the machine steps or halt.

|= prog p : Π Π instr Γ{ι}Γ r : Γ step-or-halt(p, r, ι)
Preservation: if a well-typed block steps, there is an environment type-checking the next instruction.

|= prog p : Π env Γ r : Γ Π; Γ block ι (r, ι) p → (r , ι ) ∃Γ . env Γ ∧ r : Γ ∧ Π; Γ block ι
Soundness: a well-typed program, after a finite amount of steps, either halts or can make another step.

|= prog p : Π Γ = Π(l 0 ) ι = p(l 0 ) (r, ι) p → * (r , ι ) step-or-halt(p, r , ι )
Note that all the above properties are existential and while this is a non-issue for (constructive) theorem proving, it may be a challenge for testing, not only in the functional setting, but also when logic programming is concerned, as we shall see in Section 3.1.

The Twelf implementation adds more properties, mostly linked to its own peculiar meta-theory. So, what shall we test? Past experience with αCheck [START_REF] Cheney | αCheck: A mechanized metatheory model checker[END_REF] suggests that testing intermediate lemmas is beneficial, while systems such as PLT-Redex tend to go for the bull's eye. We report in Section 5.2 some partial answers to this dilemma.

Mutation testing

How do we know that PBT is effective in catching bugs in MMT? There are very few F#/Prolog implementation of PL calculi in the wild that we can sink our metaphorical teeth in. In truth, PLT-Redex has collected a small set of bugged models, but for one, you really have to fancy Lisp2 , and secondly, it is not clear how meaningful those (very few) bugs are. Case in point, their list-machine model comes with only 3 (three) bugs that are caught by counterexamples with at most two instructions.

One way to asses the effectiveness of PBT is via mutation testing [START_REF] Pezzè | Software testing and analysis -process, principles and techniques[END_REF]. As well known, the latter is a form of white box testing, whereby a program is changed in a localized way by introducing a fault. The resulting program is called a "mutant" and the aim of a testing suite is to fail the faulted code, which is known as "killing" the mutant.

While mutation testing is widely applied to programming languages [START_REF] Jia | An Analysis and Survey of the Development of Mutation Testing[END_REF], the fact is that we do not have a theory for mutations of PL models. Standard mutation operators from (imperative) programming languages largely do not apply; hence we took some inspiration from operators for the closest we can get to models specified by derivation rules, Prolog [START_REF] Juliano | Applying Mutation Testing in Prolog Programs[END_REF]; this resulted in the following mutation operators, adapted to a strongly-typed setting:

Clause mutations: predicate deletion and swap, replacement of conjunction by disjunction.

Operator mutations: arithmetic and relational operator mutation.

Variable mutations: replacing variable by (anonymous) variable and vice versa

Constant mutations: replacing constant by constant, by (anonymous) variable and vice versa.

We thus manually produced two dozens mutations of the list-machine, half of which are specific to a relational implementation of part of the model. In Fig. 3 we show two sample mutations: the first one, ported from PLT-Redex, is a variable mutation, by which the step-cons rule fails to update the store forgetting r . The second one is a constant mutation where, removing the list type constructor, v gets the wrong type.

While we readily acknowledge that a manual approach to mutation testing is seriously limited, we point out that this is common in the field and that a general theory and implementation of a tool for mutation testing of PL calculi is beyond the scope of this paper. 3 αProlog implementation αProlog [START_REF] Cheney | Nominal logic programming[END_REF] is a logic programming language particularly suited to encoding PL calculi and related systems due to its support for nominal logic. However, since this case study is purposely first order only, we did not get to use any of those goodies and the encoding is many sorted pure Prolog code.

r(v 0 ) = a 0 r(v 1 ) = a 1 r[v := cons(a 0 , a 1 )] = r (r, (cons v 0 v 1 v ; ι)) p → ( r , ι) step-cons* Γ(v) = listcons τ Γ[v := r r list τ] = Γ Π instr Γ{fetch-field v 1 v }Γ check-instr-fetch*
In fact, it follows quite closely the reference Twelf implementation by Appel [START_REF] Appel | A List-Machine Benchmark for Mechanized Metatheory[END_REF]. The only place one would naturally try and use name types, and hence inherit α-equivalence, is encoding of variables and labels. Yet, the machine model calls for distinguished (initial) variable and label v 0 and, l 0 and this suggests an encoding based on an enumeration of constants -remember, we are in the business of bounded model checking, so we tend to avoid using infinite types as integers. The only downside is the need for an explicit inequality predicate. We use association lists for all the environments the listmachine uses (note the Haskell like syntax for type abbreviations) and standard Prolog predicates for related operations such as looking up or updating (non-destructively) such a list. Do not be fooled by the functional notation, it is flattened at compile time to the equivalent relation. The complete code can be found at https://bitbucket.org/fkomauli/list-machine.

= value. var_lookup([(V,A)|R],V) = A. var_lookup([(V,A)|R],V1) = A1 :- not_same_var(V,V1), A1 = var_lookup(R,V1).
The small-step semantics is an unsurprisingly translation of the rules in Fig. 1 Note that check_instr(Pi,G,I,G') will instantiate G' during execution. All other checks expect fully grounded inputs, i.e. no type inference is possible.

PBT with αCheck

αCheck [START_REF] Cheney | αCheck: A mechanized metatheory model checker[END_REF] is a tool for checking desired properties of formal systems implemented in αProlog. The idea is to test properties/specifications of the form H 1 ∧ • • • ∧ H n → A by searching exhaustively (up to a bound) for a substitution θ such that θ (H 1 ), . . . , θ (H n ) all hold but the conclusion θ (A) does not. In this paper, we identify negation with negation-as-failure, but the tool includes also another strategy [START_REF] Cheney | Advances in Property-Based Testing for αProlog[END_REF]. The concrete syntax for a check is

#check "name" depth: G => A.
where G is a goal and A an atom or constraint 3 . As usual in Prolog the free variables are implicitly universally quantified and depth is the user-given bound. The above pragma is translated to the formula

∃ X : τ. G ∧ gen τ ( X) ∧ ¬A (1)
where gen τ are type directed exhaustive generators automatically compiled by the tool to ground X, the intersection of the free variables of G and A, needed to make the use of NF sound. The user can, alternatively, specify her own generators as αProlog code, if she feels she needs to implement a smarter generation strategy, as we will see shortly and in Section 6.2.

A query such as (1) amounts to a derivation-first approach, which generates all "finished" derivations of the hypothesis G up to a given depth, considers all sufficiently ground instantiations of variables, and finally tests whether the conclusion finitely fails for the resulting substitution. αCheck implements a simple-minded iterative deepening search strategy over a hard-wired notion of bound, which roughly coincides with the number of clauses that can used in the derivation of each of the premises Let's look at a check to be more concrete, progress. While αProlog does not have first class existentials and disjunction, it is a basic Prolog exercise to code it: pred step_or_halt (program,store,instr). step_or_halt(P,R,instr_halt). step_or_halt(P,R,I) :-step(P,R,I,R',I').

#check "progress" 10: check_program(P,Pi), check_block(Pi,G,I), store_has_type(R,G) => step_or_halt(P,R,I).

Keeping in mind that the tool add generators for P,R,I before trying to refute step_or_halt(P, R, I), informal mode analysis for step tells us that R',I' will be ground when the partial proof tree for the check is built.

The same approach will not work for preservation:

pred exists_env (program_typing,store,instr). exists_env(Pi,R,I) :-store_has_type(R,G), env_ok(G), check_block(Pi,G,I).

#check "pres" 20 : check_program(P,Pi), step(P,R,I,R',I'), env_ok(G), store_has_type(R,G),check_block(Pi,G,I) => exists_env(Pi,R',I').

This because store_has_type(R,G) expects G to be ground and we are in effect trying to use αProlog for an impossible type inference task. The solution is to write a specialized generator build_env for type environments (and recursively, for types and variables). The peculiarity is that we need to add a local depth bound, so that smallish environments can be built independently from the hard-wired bound, which is additively distributed along all the atoms in the check.

pred exists_env_b (int,program_typing,store,instr). exists_env_b(N,Pi,R,I) :-N > 0, exists_env_b(N -1,Pi,R,I). exists_env_b(N,Pi,R,I) :-N >= 0, build_env(N,G), store_has_type(R,G), env_ok(G), check_block(Pi,G,I).

#check "pres_b" 20 : ... => exists_env_b(4,Pi,R',I').

Please see Section 5 for empirical results. Even discounting our obvious bias, specifying and validating properties ("spec'n'check") in αCheck is dead simple, requires very little effort and more than often turns out to be pretty useful.

F# implementation

The granddaddy of PBT being QuickCheck, it is natural to look into a functional implementation of the benchmark to establish a baseline for αCheck's efficacy. We chose FsCheck (https://fscheck. github.io/FsCheck), partly because we were familiar with it, but mostly because, differently from other libraries, we hoped we could leverage FsCheck's ability to provide automatic data generation for any datatype via .Net reflection to start spec'n'check without further ado. How wrong we were.

The functional implementation of the machine is unremarkable: we use F#'s maps for all sorts of environments and integers for variables and labels. The step function is total, since type-checking see to that. The type-checking function follows the imperative specification in the paper (Section 8.2).

This is what a check such as progress looks like in FsCheck's DSL, keeping in mind that ==> operationally means: (lazily) pass to the post-condition all the tests that satisfies the pre-conditions and discard the other up to a (configurable) limit:

(typecheck-program pi p) && (typecheck-block pi g i) && (store-has-type r g) ==> lazy (step-or-halt p r i)

Having carefully read [START_REF] Hritcu | Testing noninterference, quickly[END_REF], we were expecting low coverage for checks as the above, since uniform distributions are not likely to find data that has to satisfy severe constraints; still, we were not quite ready for the number we got: zero. FsCheck, as typical in adaptations of QuichCheck, provides a monadic language to write your own generators and so we did. In [START_REF] Hritcu | Testing noninterference, quickly[END_REF], the emphasis was writing generators so that the machine would run longer without getting stuck; here we need to produce well-typed programs, and this is far from immediate, since we have no type-inference whatsoever. For the progress property, we need to generate simultaneously a program p, a program typing pi that type-checks with p, a store r compatible with a type environment g, a label l that belongs to program p and the instruction i associated to label l.

Rather than showing the code, which is available in the repo https://bitbucket.org/fkomauli/ list-machine/branch/fsharp, we sketch some of the strategy. The starting point is giving more bias to cons instructions, so as to populate the store giving a chance for branching and fetching to do something interesting. We choose an instruction only if we know it is safe to execute in the given type environment. With the instruction itself, the generators produces the environment updated after the execution of the instruction. We then build sequence of non-terminal instructions. Jumps are delicate: a jump is always directed forward to a random label with id greater than the current block. If the chosen landing block has an incompatible typing environment or there is no label to jump forward to, then an halt instruction is generated instead. And this is just part of the reasoning behind it.

The smart generators developed for this benchmark consists of over a hundred lines of very dense F# code. Even if it sounds like a small number, the time spent behind it was dozens of hour dedicated to studying data correlation and distribution to generate meaningful objects. The process continued by improving distributions and fixing bugs that were found out by analyzing coverage. Though not a formal proof, we also validated with FsCheck itself the soundness of generators, as far as providing data that satisfies the constraints we imposed. Completeness is out of the question, but even in systems such as QuickChick [START_REF] Paraskevopoulou | Foundational Property-Based Testing[END_REF] where such proofs are expressible, this is far from automatic.

Shrinkers are a necessary evil of random generation, transforming large counterexamples into smaller ones that can be understood and acted upon. The QuickCheck philosophy across the board is to put in the hand of the user all the hard choices and implementing shrinkers is one. While shrinking blocks can be achieved by removing non-terminal instructions, we must safeguard several sources of correlation between data: for example, programs and program typings cannot be shrunk independently, as we must ensure that they still type-check. For the latter we wrote a type inference function to generate them together from scratch. This strategy worked well enough to reduce counterexamples to sizes sufficiently close to the minimal ones easily found with the iterative deepening approach of αCheck.

Experimental results

Reasons of space and the desire not to bore the reader silly suggest to present only a selection of all the experiments that we have carried out. Full details can be found in [START_REF] Komauli | Property-Based Testing Abstract Machines[END_REF].

A cautionary tale

A testing approach is any good only if it uncovers bugs, but we were not expecting to find any in the model presented in the paper, which came equipped with a type soundness proof formalized in two proof assistants. So, imagine our surprise when FsCheck, and later on αProlog came up with this counterexample to type preservation: p = (l 0 : cons(v 0 , v 0 , v 0 ); jump l 1 ); (l 1 : fetch-field(v 0 , 0, v 0 ); jump l 2 ); (l 2 ; halt)

Π = (l 0 : [v 0 → nil]); (l 1 : [v 0 → listcons nil]); (l 2 : [v 0 → nil]) r = [v 0 → cons(cons(nil, nil), nil)] i = fetch-field(v 0 , 0, v 0 ); jump l 2 (instruction at l 1 )
which after a single step yields:

r = [v 0 → cons(nil, nil)] i = jump l 2
However, there can be no Γ satisfying the postconditions: any such type environment must contain either the binding [v 0 → listcons τ] or [v 0 → list τ] to accommodate r , but both would not be compatible with what is required by the jump. After a fair amount of soul searching, we zeroed on the encoding of the the judgment value-has-ty, (page 481 of [START_REF] Appel | A List-Machine Benchmark for Mechanized Metatheory[END_REF]). nil : nil nil : list τ cons(a 0 , a 1 ) : listcons τ a : listcons τ a : list τ This is an essential component of the definition of r : Γ, which can be stated as for all v ∈ dom(r), it holds r(v) : Γ(v). The listcons case sounds fishy, since it makes no assumption about the types of a 0 and a 1 . Once we changed that case to: a 0 : τ a 1 : list τ cons(a 0 , a 1 ) : listcons τ the counterexample failed to show up. This change was also confirmed by inspecting the Coq implementation (it is left as an exercise in the Twelf one), which defines value-has-ty exactly like that.

Less dramatically, the check-instr-branch-listcons typing rule in the original paper contains additional preconditions similar to those in the check-instr-branch-list rule, regarding the jump target environment; however in the Coq implementation they are not present, and they could not be, considering the proof of equivalence with the algorithmic version of type checking, where they are notably absent. Several typos also present, viz. rule check-instr-cons0 (Section 8.1) and in the type checking algorithm (instruction typecheck instr Π Γ ι). Typos were spotted through formalization, not PBT.

The moral is, of course, that if there is no formal connection between a formalization and the paper reporting it, Isabelle/HOL and literate Agda being the precious exceptions, you may want to be skeptical of the latter. Similar findings appear in [START_REF] Casey | Run your research: on the effectiveness of lightweight mechanization[END_REF]. We used our "home-baked" mutation testing to assess the efficacy of PBT suites implemented with αCheck and FsCheck w.r.t. the 11 mutations applicable to both implementation. Checks are divided in three groups: top-level theorems as in our Section 2.3, intermediate lemmas collected in Section 8 of [START_REF] Appel | A List-Machine Benchmark for Mechanized Metatheory[END_REF] and auxiliary checks regarding low-level details of the machine specification. We have also ported from PLT-Redex to both αProlog and F# few dozens unit tests, as an additional baseline for mutant killing. In doing, so we have made some of them parametric, so that exhaustive and random generation of those parameters could make those tests more far-reaching.

Mutation analysis

Checks that took more than a minute to find a counterexample are considered timed-out 4 . We do not compare times to find counterexamples in αCheck vs. FsCheck, because it is a serious case of oranges and apples. We did take timing information's separately, and again we refer to [START_REF] Komauli | Property-Based Testing Abstract Machines[END_REF] for the full breakdown. In summary, the FsCheck implementation of PBT is quick (as promised, although not as much as you would expect), not only because of the efficiency of the host language, but also thanks to the flexibility of the configuration of how checks are executed. The whole FsCheck BPT suite completes in around 5 minutes: roughly 10 seconds are taken by top-level checks and about 40 seconds by the intermediate lemmas, while the remaining time is taken by the auxiliary checks. Instead, execution time of the αCheck suite is not as predictable: we manually set the depth bound so as to make the execution stay withing 10 seconds. Exception to this discipline are the progress and preservation checks, whose search space size tend to explode several minutes to complete. Unit tests, even in a parametric fashion, execute in just few seconds in both implementations. What we can reasonably compare is the rate of mutants killed by the two testing approaches (Fig. 4) and the αCheck implementation comes marginally ahead. Lemmas and low-level checks did not perform well, while top-level theorems were capable of killing most mutants. In fact, the soundness property was the one which killed most, as we can see in Fig. 5 and it is a good countermeasure against possible errors contained in the model. This is consistent to the PLT-Redex model, where soundness is the only property checked.

Progress

While we do not want to read too much in such a limited experimentation, it is safe to say that αCheck keeps its ground, considering how inexpensive it is to set up.

Extensions

Here we very briefly report on extending the list machine model to its 2.0 version, see http://www.cs. princeton.edu/ ~appel/listmachine/2.0/ and on replicating some of the results in [START_REF] Hritcu | Testing noninterference, quickly[END_REF]. In both cases, we fixed αCheck as the testing approach: in the first case, to showcase how nice it is not having to adapt generators and shrinkers and in the second one because the random approach has already been thoroughly investigated by much smarter people than us.

List-machine 2.0

In the final part of the paper [START_REF] Appel | A List-Machine Benchmark for Mechanized Metatheory[END_REF], the model is extended to cater for indirect jumps. This entails some small but far-reaching changes: labels are now coerced into values and the nil value is replaced by the initial label l 0 . We generalize the jump instruction and add a way to get the current label. This yields the following changes to the operational semantics:

r(v) = l p(l) = ι (r, jump v) p → (r, ι ) step-jump r[v := l] = r (r, get-label l v; ι) p → (r, ι) step-get-lab
We modify the type system for the new instruction(s):

Π(l) = Γ 1 Γ ⊂ Γ 1 v ∈ dom(Γ 1 ) Π; Γ block get-label l v; jump v check-block-jump2
Updating the αProlog implementation to reflect those changes is a matter of half an hour and we refer to the online documentation for details (https://bitbucket.org/fkomauli/list-machine/2.0). As the paper does not detail the soundness proof, the first thing we did was running our PBT suite to validate our implementation of the extension.

Preservation as formulated beforehand fails! This time, it is not a specification error, but the fact that the statement of the theorem has to be generalized. In fact, in the new typing rule both instructions have to occur in sequence and if we take a single step after get-label, no typing rule will match the remaining jump, thus making the typing requirement in the existential conclusion fail. One solution is to modify the latter so as to try another optional step to "get over" the jump. After this modification, αCheck did not find any other counterexamples w.r.t. the 2.0 model.

Testing non-interference, not so quickly

It is natural to consider, in addition to type soundness, more intensional properties, such as dynamic secure information-flow control (IFC). This choice is not casual, since this is the topic of [START_REF] Hritcu | Testing noninterference, quickly[END_REF], which uses classic QuickCheck. The paper claims that, for the hardest-to-find bugs, minimal counterexamples are well beyond the scope of naive exhaustive testing. Thus, we want to test this. We mainly concentrate on Section 1-4 of [START_REF] Hritcu | Testing noninterference, quickly[END_REF] and refer the reader to the paper for full details of the model and their testing outcome.

The starting point is a simple, but not simplistic abstract stack-and-pointer machine with instructions such as push, pop, load and store. A machine state consists of a program counter, a stack, a randomaccess memory, and a fixed instructions sequence. In dynamic IFC, security levels (labels, here only public ⊥ and private ) are attached to run-time values and propagated during the execution, making sure that private data does not leak. In this model, the values manipulated by the machines are labeled integers x@L. The property that the abstract machine must satisfy is end-to-end noninterference: if we call indistinguishable two machine states if they differ only in their private values, non interference guarantees that in any execution starting with indistinguishable states and ending in a halted state, the final states are indistinguishable as well.

In [1, Section 2.3], the authors present an operational semantics for the machine, which, while intuitive, turns out to be bugged; then they detail (Section 4) the sophisticated testing strategies they had to program to catch the four bugs inserted. To give an idea, on the left we report the buggy rule for Load it is natural to try and make the search strategy more modular, possibly using the notion of hookable disjunction from TOR [START_REF] Schrijvers | Tor: Modular search with hookable disjunction[END_REF]. Orthogonally, model-checking is all about pruning the search space and in the context of MMT, ideas from [START_REF] Roberson | Efficient software model checking of soundness of type systems[END_REF] could be valuable.

Finally, it is folklore that linear logical framework are well suited to represent in an elegant and concise way PL calculi with imperative features, as initially suggested by [START_REF] Cervesato | A Linear Logical Framework[END_REF]. However, when it comes to reason about them, the offering is quite meagre, with a few exceptions [START_REF] Momigliano | A formalization of an Ordered Logical Framework in Hybrid with applications to continuation machines[END_REF]. This is due to the (current) lack of meta-reasoning facilities in the leading (well, the only) sub-structural logical framework, Celf. An alternative worth pursing hence is validation via PBT in such a setting. In [START_REF] Miller | Property-Based Testing via Proof Reconstruction: Work-in-progress[END_REF], the authors have presented a proof theoretical reconstruction of PBT for relational specifications: this leverages the Foundational Proof Certificate framework [START_REF] Chihani | Foundational Proof Certificates in First-Order Logic[END_REF] to describe both the generation and the testing phase in terms of focused search. The idea translates directly to linear logic as well and can be immediately prototyped by encoding what it boils down to very simple Lolli-like meta-interpreters in a host language such as λ Prolog. The conjecture is that a linear encoding would help particularly exhaustive data generation by making structures such as stores and type environments implicit and turning them into (atomic) logical assertions.

Figure 2 :

 2 Figure 2: Instruction and block typing

Figure 3 :

 3 Figure 3: Sample mutations

Figure 4 :

 4 Figure 4: Mutation analysis for all properties and unit tests

Figure 5 :

 5 Figure 5: Mutation analysis for top-level theorems

  into a predicate step (program, store, instr, store, instr), where the first three arguments are inputs. Type checking is slightly more interesting: we summarize the main type definitions:

	pred check_instr	(program_typing,env,instr,env).
	pred check_block	(program_typing,env,instr).
	pred check_blocks (program_typing,program).
	pred check_program (program,program_typing).

The Twelf "proof" is left as an exercise and the above typing judgment undefined.

It is not just a question of syntax: the code for the only check in[START_REF] Findler | The PLT-Redex list-machine model[END_REF], namely progress, is so procedural to make it hard to relate to the formal statement in the paper.

Since in this paper we make no use of nominal features, this means syntactic equality.

Checks executed on machine with Intel Core i5-4-200U CPU,a clock speed of 1.60GHz and 8GB of RAM, with 64 bit Ubuntu 17.10 Artful and Linux kernel 4.13.0

and on the right the fixed one, where i(pc) = Load and the box signals where the issue is. pc, (x@L : s), m =⇒ pc, ( m(x) : s), m Load * pc, (x@L : s), m =⇒ pc, ( m(x)@L : s), m Load OK

QuickCheck managed to locate the first two bugs with a relatively naive generation strategy by which indistinguishable states where generated together by randomly creating the first and modifying the second in their secret part. The other two bugs necessitated a far more complex strategies for generating meaningful sequences of instructions and addresses, akin to the ones we used for the FsCheck implementation of the list machine.

αCheck found the first two bugs in less than a minute without any setup. The third one required writing a generator that yields more structured programs, such as sequences of push and store so that values in memory can change during execution, a necessary condition to find distinguishable states. This generator is a simplification of the weighted and sequence strategy in [START_REF] Hritcu | Testing noninterference, quickly[END_REF]. The fourth bug seemed, however, out of reach of αCheck, until we got it instead using a plain αProlog query for non-interference, rather than a check, the difference being the search strategy: depth-first vs. un-optimized iterative deepening. After the last bug has been fixed, the query did not find any counterexample and completes in about 5 minutes.

Section 2 in the paper ends by introducing three additional and more outlandish bugs. We got them all with the previous techniques (vanilla check, check with a generator and query with a generator) with queries being the more efficient one. Additionally, we have also carried out some experiments with a second version of the machine that includes a jump instruction [1, Section 5], but the results are too preliminary to draw any conclusion.

Conclusions and future work

Our experience suggests that off-the-shelf PBT tools are already quite useful in validating the metatheory of PL models, may they be already formalized as the list machine and their evolution, or the more if under under development. PBT helps in finding errors in the specifications and also in adjusting the statements of theorems when the model changes, as in the case of list-machine 2.0. From the costs/benefits perspective, exhaustive generation, even in the naive way we have considered, seems to be a winner over the random approach. Validating low-level languages brings in more challenges, but those can be handled with the tools we have and some additional work. Clearly, there are many other TAL models, see Chapter 4 of [START_REF]Advanced Topics in Types and Programming Languages[END_REF], that could confirm this conclusion.

We are keenly aware that FsCheck and αProlog are the extremes of a small number of PBT-tools that we could have used for this case study, starting from running PLT-Redex on our mutations, to other logic programming model-checkers such as Bedwyr [START_REF] Baelde | The Bedwyr System for Model Checking over Syntactic Expressions[END_REF]. The present paper is not meant to be an exhaustive (sorry for the pun) survey of any applicable tool. Still, a gap that should be filled is replaying the benchmark with approaches that goes beyond pure generate-and-test and try to automatically derive (random) generators that intrinsically satisfy certain pre-conditions. The obvious candidate is the new QuickChick [START_REF] Lampropoulos | Generating good generators for inductive relations[END_REF], but also Bulwahn's smart generators in Isabelle/HOL [START_REF] Bulwahn | The New Quickcheck for Isabelle -Random, Exhaustive and Symbolic Testing under One Roof[END_REF] are a possibility.

αCheck performs better than we hoped for, considering that its implementation is nothing more than a OCaml interpreter for nominal logic programming, exploring the full search space in an iterative deepening way up to a bound. In other terms, it is probably orders of magnitude slower than standard Prolog and it makes no effort to prune the search space or to explore it in more flexible ways. As the preservation example shows, a hard-wired additive bound is, to say the least, inconvenient and therefore


Zhao Li

C.-J Richard Shi
email: cjshi@ee.washington.edu

An Efficiently Preconditioned GMRES Method for Fast Parasitic-Sensitive Deep-Submicron VLSI Circuit Simulation*

We propose an efficiently preconditioned generalized minimal residual (GMRES) method for fast SPICE-accurate transient simulation of parasitic-sensitive deep-submicron VLSI circuits. First, when time step-sizes vary within a predefined range, the preconditioned GMRES method is applied to solve circuit matrix equations rather than LU factorization. The preconditioner we use comes directly from the previously factorized L and U matrices. Second, to keep using the same preconditioner during nonlinear iteration, the successive variable chord method is applied as an alternative to the Newton-Raphson method. An improved piecewise weakly nonlinear definition of MOSFETs is adopted and the low-rank update technique is implemented to refresh the preconditioner efficiently. With these techniques, the number of required LU factorizations during transient simulation is reduced dramatically. Experimental results on power/ground networks have demonstrated that the proposed method yields SPICE-like accuracy with an about 18X overall CPU time speedup over SPICE3 for circuits with tens of thousands elements.

Introduction

In modern deep-submicron VLSI circuit design, parasitic effects are no longer ignorable with the higher operation frequency, lower supply voltage and smaller device feature size [START_REF] Phillips | Simulation Approaches for Strongly Coupled Interconnect Systems[END_REF]. Accurate post-layout verification requires full-chip simulation of large-scale circuits together with massive extracted parasitic elements, which come from substrate, power/ground networks, interconnects, etc. For such kind of circuits, the per-iteration cost with SPICE [START_REF] Nagel | SPICE: A Computer Program to Simulate Semiconductor Circuits[END_REF] during transient simulation is dominated by costly LU factorizations. New approaches are required for fast-yetaccurate simulation of parasitic-sensitive circuits. For example, several efficient methods [START_REF] Chen | Efficient Large-Scale Power Grid Analysis based on Preconditioned Krylov-subspace Iterative Methods[END_REF][4][5] [START_REF] Zhao | Hierarchical Analysis of Power Distribution Networks[END_REF] have been proposed to speedup the simulation of power/ground networks. These methods, although orders of magnitude faster than SPICE, are mainly tailored for purely linear circuits and have difficulties with nonlinear circuits incorporated for full-chip simulation. Recently, [START_REF] Li | A Coupled Iterative/Direct Method for Efficient Time-Domain Simulation of Nonlinear Circuits with Power/Ground Networks[END_REF] proposed to couple nonlinear circuits together with power/ground networks based on the Gauss-Seidel style relaxation [START_REF] Mccalla | Fundamentals of Computer-Aided Circuit Simulation[END_REF]. However, when nonlinear and linear circuits are strongly coupled, the number of nonlinear iterations could be very large.

A key approach to improving the efficiency of SPICEaccurate simulation of large-scale parasitic-sensitive circuits is to reduce the number of LU factorizations during time-domain simulation. In [START_REF] Acar | TETA: Transistor-Level Waveform Evaluation for Timing Analysis[END_REF], a fixed time step-size is assumed and the successive chord method [START_REF] Mccalla | Fundamentals of Computer-Aided Circuit Simulation[END_REF] is applied to linearize nonlinear devices. Thus, only one LU factorization is required during the whole transient simulation. However, the application of this method is limited since excessive nonlinear iterations might be required with only a single chord defined for the entire operating region of nonlinear devices. Further, the assumption of a fixed time step-size is not adequate for circuits with widely distributed time constants.

Recently, we proposed SILCA [START_REF] Li | SILCA: Fast-Yet-Accurate Time-Domain Simulation of VLSI Circuits with Strong Parasitic Coupling Effects[END_REF] for VLSI circuits containing strong parasitic coupling effects. Two linearcentric ideas were used to help keep a circuit matrix constant during transient simulation with variable time step-sizes: 1) Semi-implicit iterative integration scheme to keep equivalent conductance of capacitor/ inductor companion models constant for a large time interval; 2) Successive variable chord method to keep linearized conductance of nonlinear devices constant for a large voltage/current range. With these two ideas, the number of LU factorizations can be reduced dramatically. However, as pointed out in [START_REF] Li | SILCA: Fast-Yet-Accurate Time-Domain Simulation of VLSI Circuits with Strong Parasitic Coupling Effects[END_REF], the absolute stability region of the iterative integration corrector is related to the number of iterations. Consequently, the stability of the iterative trapezoid formula could become worse than that of the standard trapezoid formula in practice. Therefore, heuristic measures have been incorporated in SILCA to ensure the A-stability at the cost of more iterations and/or more LU factorizations.

The proposed preconditioned GMRES method [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] in this paper borrows the ideas in SILCA to reduce the number of LU factorizations during variable time step-size transient simulation. Instead of using semi-implicit integration predictor and iterative integration corrector in SILCA, standard integration formulas are applied here. Thus, the stability problem encountered in SILCA is naturally avoided. Furthermore, the number of nonlinear iterations with the preconditioned GMRES method is smaller than that with SILCA and could be comparable to that with SPICE.

It is well known that the key to fast convergence of the GMRES method is to design an efficient preconditioner [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. A good preconditioner should be as close to the inverse of a circuit matrix as possible and easy to derive. In this paper, we reuse the previously factorized L and U matrices as the preconditioner in transient simulation. The details are described as below, 1) When time step-sizes h n vary within a predefined range of the basis time step-size h, the preconditioned GMRES method rather than LU factorizations is applied to solve circuit matrix equations. The preconditioner comes directly from the previously factorized L and U matrices base on the basis time step-size h. 2) To keep using the already factorized L and U matrices as the preconditioner during nonlinear iteration, we apply the successive variable chord method [START_REF] Li | SILCA: Fast-Yet-Accurate Time-Domain Simulation of VLSI Circuits with Strong Parasitic Coupling Effects[END_REF] as an alternative to the Newton-Raphson method and improve the piecewise nonlinear definition of MOSFETs. The low-rank update technique is implemented to refresh the preconditioner efficiently. Further, incomplete L and U matrices can be derived to act as the preconditioner for the better performance. With these techniques, the GMRES method is able to converge in a small number of iterations and the number of required LU factorizations is reduced dramatically.

This paper is organized as follows. Section 2 presents the improved piecewise nonlinear definition of MOSFETs. The new preconditioned GMRES method is proposed in Section 3. Section 4 summarizes experimental results on general nonlinear circuits and power/ground network examples. Finally, conclusions are drawn in Section 5.

PWNL definition of MOSFETs

In SILCA, a heuristic piecewise weakly nonlinear (PWNL) definition of MOSFETs was proposed. The purpose was to keep the circuit matrices constant when nonlinear iterations are performed within a PWNL region. Therefore, the number of LU factorizations can be reduced. In this section, we start from discussing the convergence property of the successive variable chord (SVC) method. Then systematic rules to generate the PWNL regions of MOSFETs are described. Finally, we briefly review the low-rank update technique.

Suppose that nonlinear iterations are performed within a PWNL region of a nonlinear function f(x) to solve f(x)=0, as shown in Fig. 1, nonlinear iteration can be expressed by,

g x f x x i i i ) ( 1 (1)
where g is the chord for this PWNL region. Let the exact solution be x * = x i + i = x i+1 + i+1 . Subtracting x * from both sides of Eq. (1) gives,

g x f i i i ) ( 1 (2)
By the Taylor expansion of f(x) at x i , we obtain the following error estimation,

g x f g x f i i i i i 2 ) ( ' ' ) ) ( ' 1 ( 2 1
(3) Eq. ( 3) shows clearly that a quadratic convergence rate is achieved if g is equal to f'(x i ), which is the Newton-Raphson method. Otherwise, the convergence rate is reduced to be linear, which is the case for the successive variable chord method. We observed that, on one hand, the smaller the |1f'(x i )/g| is, the closer to the quadratic convergence rate Eq. ( 3) is. On the other hand, the larger the |1f'(x i )/g| is, the larger the range of a PWNL region could be. Apparently, there exists a tradeoff between the convergence rate and the range of a PWNL region. We define the following condition with a parameter <1,

g x f i ) ( ' 1 (4)
For a PWNL region as shown in Fig. 1, it can be derived from Eq. ( 4) that, It should be noted that the above analysis is done in the context that nonlinear iterations are performed within a PWNL region. In case that nonlinear iterations run across two or more PWNL regions, such as the example shown in Fig. 2 where the exact solution resides at the boundary of two PWNL regions, the following condition should be satisfied to achieve convergence, 0 2

x x ( 6 ) Thus, g 1 and g 2 should satisfy the following inequality,

a x x x f x f a g a g 0 1 0 1 2 1 ) ( ) ( ) ( ) ( (7)
In our experiments, to satisfy both Eq. ( 5) and Eq. ( 7), the chord is chosen to be the maximum derivative in each PWNL region. With the knowledge of device model behaviors, such as monotonicity, the maximum derivative for each PWNL region can be computed. It should be noted that PWNL regions of a nonlinear function is equivalent to piecewise constant (PWC) regions of the first-order derivatives of the nonlinear function. The following three rules can be used to generate PWNL regions for the MOSFET model.

1) The maximum voltages of V ds and V gs are predefined. In our experiments, we use V dd as the maximum voltage for both of them. Given model parameters, the maximum g ds and g m can be calculated. 2) With a predefined <1, the PWC regions for g ds and g m are calculated as below, 3) A lower bound of g ds and g m is predefined, so that the rule (2) will stop whenever g ds and g m are less than the predefined lower bound. This is necessary to avoid a PWC region for g ds and g m being too narrow. With the above rules, the PWC regions for g ds and g m of the MOSFET level 1 model in the two-dimensional V ds and (V gs -V th ) plane are shown in Figs. 3 and4, respectively. There are five PWC regions for both g ds and g m .

It is clear from Figs. 3 and4 that g ds and g m reach their maximum values in different PWNL regions. It should be noted that effects due to V bs have been incorporated into V th . For the MOSFET level 1 model, g mbs has a simple relationship with g m [START_REF]Star-Hspice Manual[END_REF], For the simplification purpose, we always use the maximum in our experiments. Therefore, the PWC regions for g sb th dV dV / mbs are the same as those for g m . The rules of generating PWNL regions can be applied to complicated MOSFET models such as BSIM3 [START_REF]Star-Hspice Manual[END_REF], as well.

As mentioned in [START_REF] Fujisawa | A Sparse Matrix Method for Analysis of Piecewise-Linear Resistive Networks[END_REF], the low-rank update technique is an efficient method to update the factorized L and U matrices when only a few nonlinear devices change their PWNL regions. When a MOSFET switches its operating PWNL region, the contribution of this MOSFET to the circuit matrix changes as follows,

Preconditioned GMRES method

The transient simulation flow of the proposed preconditioned GMRES method is shown in Algorithm I described below. It is clear from Algorithm I that LU factorizations are only performed when time step-sizes vary out of the predefined h n /h range. In other cases, the L and U matrices are either kept constant or updated by the low-rank update technique when nonlinear devices switch their piecewise nonlinear regions. During the whole process, the L and U matrices are used for forward/ To examine the efficiency of the proposed GMRES method, a power/ground network example as shown in Fig. 5 is simulated, which is similar to that used in [START_REF] Li | SILCA: Fast-Yet-Accurate Time-Domain Simulation of VLSI Circuits with Strong Parasitic Coupling Effects[END_REF]. The power and ground supply networks are modeled as two RCL mesh layers (parasitic coupling capacitors are not shown in Fig. 5). In our example, between these two layers is a 20-stage inverter chain, representing nonlinear circuits. Furthermore, RCL loads are added for each inverter to model interconnect lines between the adjacent stages. The size of two RCL meshes can be changed to vary the number of elements.

Tables II, III and IV summarize the simulation results for power/ground network examples using the GMRES method with the full LU preconditioner, the type I ILU preconditioner and the type II ILU preconditioner, respectively. The error tolerance for the GMRES method is set to 1e-8. SPICE3 simulation results are also included in Table II. For clarity, the run time comparison is shown in Fig. 6. It is expected that more speedup could be achieved for larger power/ground networks. It is seen that the GMRES method with the type II ILU preconditioner achieves the best speedup over SPICE3 for the largest power/ground network -18.02X. The reason is that the number of matrix elements in the type II ILU preconditioner is much less than those in the LU preconditioner and the type I ILU preconditioner, especially for large matrices. For the power/ground network example with 4002 elements, the histograms of the number of L and U matrix elements during transient simulation are shown in Fig. 7. The number of matrix elements in full L and U matrices is 3116915 for this example. It is observed that the number of L and U matrix elements in the type II ILU preconditioner is reduced to about 1/10~1/5 of that in the full LU preconditioner.

With the error tolerance set to 1e-8, the average number of GMRES iterations in each GMRES solving process (#GMRES Iter / #GMRES) is about 3.20 to 3.35 for the GMRES method with the full LU preconditioner as shown in Table II. It increases to about 3.75 to 4.50 and 4.05 to 4.75 for the GMRES method with the type I ILU preconditioner (Table III) and the type II ILU preconditioner (Table IV), respectively. It is clear that the proposed preconditioner is efficient for the GMRES method during time-domain VLSI circuit simulation. In our experiments, when the error tolerance is further decreased to 1e-10 for higher accuracy, the average number of GMRES iterations in each GMRES solving process increases to about 6.60 to 8.35 for the GMRES method with the type II ILU preconditioner. As a result, the GMRES method requires more run time when the error tolerance is made smaller. Therefore, there exists a tradeoff between the accuracy and efficiency. Figure 8 shows the output waveform of the inverter chain between the power and ground networks. It is seen that the low-level voltage of the output is larger than the ideal ground voltage due to the IR drop and L*dI/dt effects. It is clear in Fig. 8 that the accuracy with the proposed GMRES method is comparable to that with SPICE3.

It is worthy noting that the number of required LU factorizations (#Tran LU) is reduced dramatically with the preconditioned GMRES method compared to SPICE3 (#Tran LU = #Tran Iter). Furthermore, compared to the simulation results with SILCA, a similar speedup over SPICE3 has been achieved by the preconditioned GMRES method. However, the proposed preconditioned GMRES method is free of numerical stability problems as encountered in SILCA. Although the number of iterations (#Tran Iter) with the preconditioned GMRES method is increased to about 1.5X due to the PWNL definition of MOSFETs, it is still much less than that with SILCA.

Conclusion

In this paper, an efficiently preconditioned GMRES method is presented to speedup the transient simulation of parasitic-sensitive deep-submicron VLSI circuits. The cost of LU factorizations is minimized since they are only performed if time step-sizes vary violently. When time step-sizes change within a predefined range, the GMRES method is invoked with the preconditioner coming from the previously factorized L and U matrices. An improved PWNL definition of MOSFETs is also proposed to reduce the number of nonlinear iterations. With these techniques, the cost of required LU factorizations has been reduced dramatically. Orders of magnitude speedup has been achieved on power/ground network examples with the SPICE-like accuracy. The speedup could be further boosted by the parallel implementation of the GMRES method [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF].

Figure 1 .

 1 Figure 1. SVC method for f(x) within a PWNL region.

r e g io n 2 P W N L r e g io n 1 Figure 2 .

 212 Figure 2. SVC method for f(x) near the boundary of two PWNL regions.

  the above representation, the low-rank update technique can be applied.

Figure 3 .

 3 Figure 3. PWC regions of g ds in two-dimensional V ds -(V gs -V th ) plane.

Figure 4 .

 4 Figure 4. PWC regions of g m in two-dimensional V ds -(V gs -V th ) plane.

Figure 5 .

 5 Figure 5. The power/ground network example.

Figure 6 .Figure 7 .

 67 Figure 6. Run time variation with the number of elements in P/G network examples

Figure 8 .

 8 Figure 8. The output waveform of the power/ground network example

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

* This research was supported by DARPA NeoCAD Program under Grant No. N66001-01-8920, an NSF CAREER Award under Grant No. 9985507, and SRC/NSF's Joint Mixed-Signal Initiative under Contract CCR0120371.

backward substitution (FBS) and act as the preconditioner for the GMRES method.

Three types of preconditioners are tested in our experiments: 1) The full L and U matrices. 2) Type I incomplete L and U (ILU) matrices approximated from the full L and U matrices -a matrix element a(i,j) in the L or U matrix is removed if |a(i,j)|<c•|a(i,i)| and |a(i,j)|<c•|a(j,j)|. 3) Type II incomplete L and U matrices approximated from the full L and U matrices -a matrix element In Algorithm I, the PWNL definition of MOSFETs described in Section 2 is used for the preconditioner. However, it is not necessary for the GMRES method to use the same PWNL definition to construct the circuit matrix equations. Instead, the GMRES method could still use standard MOSFET models so that accurate model information is included. This is especially useful for incorporating nonlinear capacitors, since nonlinear capacitors are generally simplified to linear ones when used for the preconditioner. Furthermore, by using standard MOSFET models for the GMRES method, the number of nonlinear iterations with the preconditioned GMRES method could be close to that with SPICE.

Experimental results

General nonlinear circuit examples

To verify the proposed GMRES method on general nonlinear circuits, several digital, analog and RF circuits have been tested and results are shown in Table I. The preconditioner for the GMRES method is the full L and U matrices. From Table 1, we see that the number of LU factorizations is reduced dramatically compared to that with SPICE. For the simplification purpose, we used the PWNL definition of MOSFETs for both the preconditioner and circuit matrix equations solved by the GMRES method. Therefore, the number of nonlinear iterations is generally increased to less than 2X of that with SPICE. Compared to SILCA, the proposed GMRES method generally requires less number of nonlinear iterations.


Romain Couillet

Gilles Wainrib

Harry Sevi

Tiomoko Hafiz

Ali

Hafiz Tiomoko Ali

hal-01633450v2 TRAINING PERFORMANCE OF ECHO STATE NEURAL NETWORKS

Keywords: Neural networks, ESN, random matrix theory

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

INTRODUCTION

In the realm of both artificial and biological neural networks, one usually differentiates long-and short-term memory networks. While the former are inherently based on a dedicated rewiring of the network during training phase (i.e., the network edges adapt to the input-to-output training by means of backpropagation) and work in a feed-forward manner, the latter do not change their connectivity matrix but maintain past-input information within the network via self-connection (thereby being recurrent networks rather than feed-forward networks) so that, for energy conservation reasons (inputs are continuously fed into the networks), the past-input memory decays exponentially fast in this case. These networks, less considered than their feed-forward counterparts, have recently been reinstated by Jaeger [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks-with an erratum note[END_REF] under the name of echostate networks (ESN), who defends their performance superiority in short-term memory settings [START_REF] Jaeger | Short term memory in echo state networks[END_REF].

A striking feature of ESN's is that the (fixed) network connectivity matrix can be chosen as a random matrix and that choice is even a very satisfactory one. This observation and intuition from random matrix theory suggests the possibility for a much expected theoretical analysis of such networks. The objective of the present article is to lay down such theoretical bases to propose a first study of the training performance of these networks.

Precisely, we shall consider an echo-state network of n nodes and will prove that, as both the training sequence length T and n grow large simultaneously, the mean-square error Couillet and Tiomoko Ali's work is supported by the ANR RMT4GRAPH ANR-14-CE28-0006. performance of the task training procedure tends to be deterministic in the limit, for all well-behaved connectivity matrices W ∈ R n×n . Then, specifying W to belong to specific classes of random matrices, we will further show that the aforementioned deterministic behavior takes on a very simple expression that exhibits the salient features of the ESN. Many consequences can be drawn of this analytical formula, some of which will be discussed here.

MAIN RESULTS

Generic W

Consider a neural network of n nodes having at each time t a joint state x t ∈ R n . We assume the following state evolution

x t+1 = W x t + mu t+1 + ηε t+1 (1)
for some initial (say empty) state at large negative t, with W ∈ R n×n the connectivity matrix, . . . , u -1 , u 0 , u 1 , . . . ∈ R the network (scalar) inputs, m the input-to-network connectivity vector, and ηε t ∈ R n an additional in-network noise of amplitude η > 0. A visual representation of such a network is depicted in Figure 1.

We seek here to evaluate the performance of the so-called training task consisting in relating the input sequence {u t } to a pre-determined output sequence r 0 , . . . , r T -1 from a given linear combination of the network states x 0 , . . . , x T -1 for a duration T . This is traditionally performed using a leastsquare regression. Precisely, we shall define

ω ≡ X X T X -1 r , T ≤ n XX T -1 Xr , T > n
where r = [r 0 , . . . , r T -1 ] T and X = [x 0 , . . . , x T -1 ] ∈ R n×T , and consider X T ω as the least-square estimate of r with mean-square estimation error

E η (u, r) ≡ 1 T r -X T ω 2 . ( 2
)
This quantity is clearly identically zero if T ≤ n and we shall thus only consider the non-trivial case where T > n.

Our objective is to understand the behavior of E η (u, r) defined in (2) in the limit where n, T → ∞, for an arbitrary (and then for a random) matrix W . The underlying idea behind our approach is that the noise ε t , which regularizes the otherwise unstable network, will tend to concentrate as n, T grow large, thereby leading to a (much desired) more deterministic behavior of the network as a whole. Considering X as a random matrix through these ε t 's, we shall provide an asymptotically consistent estimate of E η (u, r).

To ensure proper conditioning in the large n, T limit, the following assumption is needed.

Assumption 1 (Large n, T limit) As n, T → ∞,

1. lim n n T = c ∈ [0, 1) 2. lim sup n W < 1 3. lim sup n AA T < ∞ where A = M U , M = [m, W m, . . . , W T -1 m] ∈ R n×T and U = {u i-j } 1≤i,j≤T .
Item (2) is a (often too stringent) sufficient condition to ensure network stability, while Item (3) is merely a nottoo demanding technical conditioning. The aforementioned matrix U will play a key role in what follows. Note that its columns constitute delayed versions of the sequence u -(T -1) , . . . , u T -1 and thus can be seen as successive snapshots of the (T -step behind) network memory as time elapses.

With Assumption 1 in place, we have our first main result.

Theorem 1 (Asymptotic MSE) Let Assumptions 1 hold and let r ∈ R T be of O( √ T ) Euclidean norm. Then, with E η (u, r) defined in (2), as n → ∞, E η (u, r) - 1 T r T Qr → 0 almost surely, where Q ≡ (I T + R + η -2 A T R-1 A) -1 , with (R, R) a solution to the implicit system R = c 1 n tr S i-j R-1 T i,j=1 R = ∞ q=-∞ 1 T tr J q (I T + R) -1 S q
and [J q ] ij ≡ δ i+q,j and S q ≡ k≥0 W k+(-q) + (W k+q + ) T (a + = max(x, 0)).

Sketch of Proof 1

The idea behind the proof consists in applying the so-called deterministic equivalent method (see e.g., [START_REF] Couillet | Random Matrix Methods for Wireless Communications[END_REF]Chapter 6]) of random matrix theory on the random matrix 1 T X T X. That is, defining Qγ = ( 1 T X T X -γI T ) -1 for all γ > 0, the resolvent of 1 T X T X, the method consists in determining a matrix Qγ such that, as n, p → ∞, a T ( Qγ -Qγ )b almost surely vanish for all deterministic and bounded norm a, b vectors. This is performed via the Gaussian integrationby-parts and Nash-Poincaré inequality framework devised by Pastur in [4, Chapter 2]. Once this is achieved, we show that the aforementioned convergence still holds when γ = 0, which is valid as long as c < 1. Calling Q the limit of Qγ as γ → 0 provides our final result. The complete details of the proof are provided in the article [START_REF] Couillet | The asymptotic performance of linear echo state neural networks[END_REF].

Although seemingly intractable, the form taken by Q is quite interesting in itself. First, note that R and R only depend on W so that the terms R and M T R-1 M contain all the information about the W found in Q. The noise level η 2 then trades the need for regularization of Q and the need for emphasizing the part A T R-1 A versus I T + R; the former matrix indeed has its columns living in the span of the columns of U with weights imposed by the network matrix M T R-1 M . Thus, we expect the memory performance of the ESN to relate strongly to both η 2 and the effect of R and R.

An interesting corollary is found when c = 0 (i.e., n/T → 0) in which case one easily shows that R = 0 and R = S 0 . And this thus brings the almost sure convergence

E η (u, r) - 1 T r T I T + 1 η 2 U T DU -1 r → 0 with D the matrix with entry D ij = m T (W i-1 ) T S -1 0 W j-1 m (recall that S 0 = k≥0 W k (W k ) T ).
We recognize here the diagonal entries of D to be the quantity

D ii = m T (W i-1 ) T   k≥0 W k (W k ) T   -1 W i-1 m
known in the ESN literature as (the value i of) the Fisher memory curve [START_REF] Jaeger | Short term memory in echo state networks[END_REF][START_REF] Ganguli | Memory traces in dynamical systems[END_REF], an abstract mesure of the ability of the ESN to maintain an i-step old input in memory.

Having an insight on the couple (R, R) for non trivial values of c is however quite involved, and it is interesting to consider specific cases where this expression simplifies. In particular, one may choose W to be a given snapshot of a random matrix model. Since n, p are assumed large, by means of random matrix identities, this one realization will have an almost sure deterministic behavior in the limit, thereby leading to explicit approximations for E η (u, r).

We shall consider next two classical random matrix models for W modelling directed and undirected random graphs.

Non-Hermitian random W

The first model assumes a non-Hermitian structure for W but with statistical invariance when (left-or right-) multiplied by orthogonal matrices. In this case, R merely becomes c 1-c I T while R is essentially (1 -c)S 0 and we have the following corollary of Theorem 1.

Corollary 1 (Non-Hermitian Random W ) Let the Assumptions of Theorem 1 hold and take W to be random with left and right orthogonal invariance, and m of unit norm independent of W . Then

E η (u, r) -(1 -c) 1 T r T I T + 1 η 2 U T DU -1 r → 0
almost surely, where D is diagonal with

D ii ≡ 1 n tr W i-1 (W i-1 ) T S -1 0 .
In particular, if W = σZ with Z a random orthogonal and orthogonally invariant matrix, 1) .

1 then D ii = (1 -σ 2 )σ 2(i-
Corollary 1 provides much expected insights on the behavior of ESN's in the case of non-Hermitian random W . In particular, from the diagonal structure of D and the fact that the columns of U are delayed versions of the input sequence {u t }, the performance of the ESN will depend on the accuracy of the representation of r as a linear combination of weighted delays of {u t }. Note interestingly that these delays satisfy i≥1 D ii = 1, so that the ESN effectively distributes its memory capacity along the successive delays (in an exponentially decaying manner). The specific choice of the singular values of W impact these delays.

Figure 2 provides a performance comparison between Monte Carlo simulations and our theoretical results for the training of an interpolation task of chaotic data (here the Mackey-Glass model [START_REF] Glass | A simple model for phase locking of biological oscillators[END_REF]). An accurate match is observed between theory and practice with increased precision as n, T grow large, consistently with our results. Despite the apparent accuracy on this example, it is nonetheless important to stress that Theorem 1 is only valid for fixed η > 0 and growing n, T . In particular, the approximation for small η's may dramatically fail in some specific cases that we simulated. 1 This is often referred to as a Haar matrix.

Hermitian random W

A second natural example is to assume that the ESN is now an undirected graph, so that W is symmetric. In this case, the matrix R is no longer diagonal but is found to be solution of a simple fixed-point equation. Precisely, assuming W Hermitian orthogonally invariant with normalized empirical eigenvalue distribution (i.e., the measure n -1 i δ λi(W ) with λ i (W ) the eigenvalues of W ) converging to µ, we have

R ab = c t |a-b| µ(dt) q∈Z t |q| 1 T tr(J q (I T + R) -1 )
for all a, b ∈ {1, . . . , T }. Since R is a Toeplitz matrix with exponentially decaying first column entries, this equation only involves a few parameters and is easily solved. An interesting case is when µ is symmetrical (i.e., µ(-t) = µ(t)) in which case all values of R ab for a -b odd vanish. Figure 3 provides a visual comparison of R when W is a Gaussian matrix with i.i.d. entries with or without symmetry.

As opposed to the case of non-Hermitian W , the ESN output is not merely mapped to a delayed version of past inputs but a more intricate combination of them. This intuitively induces performance losses when it comes to fulfilling pure delay task. This is confirmed in Figure 4, where ESN's with or without Hermitian symmetry are trained on a τ -delay memory task on the Mackey-Glass dataset. While the performance decay incurred by increasing τ for the non-symmetric case is moderate, it is instead dramatic for symmetric matrices. As such, undirected random graphs are seen to perform poorly on related simple memory tasks.

CONCLUDING REMARKS

Several interesting remarks can be drawn from the theoretical results in Section 2. Consider for instance the scenario where W = σZ for Z a Haar matrix and in which the sought-for output r is defined as r = √ T U T b for some vector b ∈ R k and k small. Then, from Corollary 1,

E η (u, r) 1 -c b T U   I T + i≥0 (1 -σ 2 )σ 2(i-1) η 2 U T •,i U •,i   -1 U T b
which can be shown to converge to zero as η → 0 and may, for every η, be numerically minimized over σ. In particular, letting b i = α i-1 for some α ∈ (-1, 1), one can show that σ 2 = |α| minimizes E η (u, r). Thus, the network should (as one would expect) align as much as possible to the dependence of r t in u t-i for each i. For more general than expo-nential decaying relations between r t and u t-i , and in particular to account for heterogeneous memory dynamics, one may appropriately consider W = diag(σ 1 Z 1 , . . . , σ L Z L ) for different values of σ and (different sized) independent Haar Z , parametrizable upon the application context. For k not small however, or for more general models of r as a function of u (especially non linear models), the convergence of E η (u, r) to zero as η → 0 is not guaranteed. In this regime of small η's, our analysis however theoretically breaks so that the instability shortcomings already observed by Jaeger in [START_REF] Jaeger | Short term memory in echo state networks[END_REF] cannot be resolved at this point and remain an open question. Nonetheless, it clearly appears through our results that large enough η's (it can be proved enough to have η 2

n -1 2 ) induce system stability.

Generalizations of the present work encompass the extension to the testing (as opposed to training) performance of ESN's (see [START_REF] Couillet | The asymptotic performance of linear echo state neural networks[END_REF] for preliminary results) along with the challenging consideration of non-linear activation functions, i.e., generalizing (1) to x t+1 = S(W x t + mu t+1 + ηε t+1 ) for some (pointwise) sigmoid function S. The theoretical difficulties incurred by the latter may be circumvented by ideas from mean field dynamics.

Fig. 1 .

 1 Fig. 1. Echo-state neural network.

10 - 5 Fig. 2 .

 52 Fig. 2. (Normalized) MSE for the Mackey Glass one-step ahead task, W = σZ for Z Haar, n = 200, T = 400. Error bars indicate one standard deviation of noise realizations.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Upper 9 × 9 part of R for c = 1/2 and σ = 0.9 for W with i.i.d. zero mean Gaussian entries [left] and W Gaussian symmetric (Wigner) [right].


Central Message

A novel TAVR device needs to take the next steps towards success: enhanced image guidance and demonstration of reproducibility.

List of abbreviations

Text

Implanted for the first time in a human patient in 2014, the J-Valve transcatheter aortic valve replacement (TAVR) concept has scored significant advancements. In perspective, it can be an attractive option for TAVR in aortic regurgitation, with a good deal of potential for success. The current paper [START_REF] Zhu | Transapical transcatheter aortic valve replacement with a novel TAVR system -J-Valve TM in high-risk patients with severe aortic valve diseases[END_REF] represents the first report of early and follow-up outcomes in a significant number of patients treated with this device. Unquestionably, this experience remains limited in terms of number of implants and amount of patient-years of follow-up. Yet, some prominent aspects deserve to be underlined. First, this technology has been developed through the contribution of cardiac surgeons, who also currently are its only users, testifying once again the vitality and the propensity towards evolution of our discipline. Second, the separation of locators and valve -and their consequent deployment in differed times -represents a typical feature and, potentially, an effective mean to improve immediate results (perivalvular leakage, anatomically coherent positioning) especially in aortic regurgitation cases.

Nonetheless, the rate of perivalvular leakage reported herein remains important and, not surprisingly, greater in aortic stenosis cases. It appears that the effectiveness of this TAVR concept relies on optimal landing of the locators at the nadir of each native aortic cusp: irregular calcifications in aortic stenosis reasonably represent an anatomical obstacle in several cases. Missed positioning of the locators -as acknowledged by the Authors -is at the root of many instances of paravalvular leakage in both regurgitation and stenosis cases. To such respect, bicuspid valve anatomy can be considered as a foe for the success of this TAVR system, including in the aortic regurgitation cases for which it was originally intended. As it is a crucial issue, the landing of the locators at the cusps' nadirs cannot be based on tactile

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 3 Anselmi, Verhoye
feedback alone, but requires effective image-based guidance. Real-time threedimensional transesophageal echocardiography (3D TEE) may prove to be the best suited modality -and even essential -to such scope. 3D TEE has demonstrated far greater accuracy and fineness than fluoroscopy and angiography in the intraprocedural evaluation of the mechanism of aortic regurgitation after deployment of a TAVR device. Real-time 3D TEE provides information about the degree of elliptical expansion of the valve device, about the anatomical relationship between the prosthetic and native commissures (antianatomical expansion) [START_REF] Shibayama | 3D Assessment of Features Associated With Transvalvular Aortic Regurgitation After TAVR: A Real-Time 3D TEE Study[END_REF]. All of these information can be discriminating in intraprocedural decision-making. Also, TEE has been already proposed as the primary imaging modality to guide the deployment of a TAVR device [START_REF] Verhoye | Association of transaortic approach and transoesophageal echocardiography as the primary imaging technique for improved results in transcatheter valve implantation[END_REF].

In transcatheter valve therapy, success is achieved through effectiveness and reproducibility. The pathway towards success of this promising TAVR concept will probably demand not only technological advancements (miniaturization of the 27-Fr introducer, availability for retrograde delivery routes) but also full exploitation of intraprocedural 3D TEE capabilities, and -most importantly -proof of reproducibility in larger international trials. Cooperation with allied heart team professionals (anesthesiologist, sonographer cardiologist) will be therefore essential.

In the end, this pathway further stimulates us to take advantage of the surgical knowledge of cardiac anatomy, physiology and pathology, in order to continue being the cutting edge in the innovation for valve therapy.

  Amedeo ANSELMI, MD PhD, Jean-Philippe VERHOYE, MD PhD Division of Thoracic and Cardiovascular Surgery, Pontchaillou University Hospital, Rennes, France Corresponding author: Amedeo ANSELMI, MD PhD Division of Thoracic and Cardiovascular Surgery, Pontchaillou University Hospital 2, rue Henri Le Guilloux 35 033 Rennes, France Tel: (0033) 06 37 06 03 79; Fax: (0033) 02 99 28 24 96 Email: amedeo.anselmi@alice.it Conflicts of interest: None.