Olivier Commowick
email: olivier.commowick@inria.fr

Nicolas Wiest-Daesslé

Sylvain Prima

Automated diffeomorphic registration of anatomical structures with rigid parts: Application to dynamic cervical MRI

We propose an iterative two-step method to compute a diffeomorphic non-rigid transformation between images of anatomical structures with rigid parts, without any user intervention or prior knowledge on the image intensities. First we compute spatially sparse, locally optimal rigid transformations between the two images using a new block matching strategy and an efficient numerical optimiser (BOBYQA). Then we derive a dense, regularised velocity field based on these local transformations using matrix logarithms and M-smoothing. These two steps are iterated until convergence and the final diffeomorphic transformation is defined as the exponential of the accumulated velocity field. We show our algorithm to outperform the state-of-the-art log-domain diffeomorphic demons method on dynamic cervical MRI data.

Introduction

In medical image analysis, one is often confronted with the problem of registering anatomical structures containing both hard, rigid (typically, bones) and soft, non-rigid (most other tissues) parts. Such problems are met for instance when following-up spinal cord lesions in MRI for the diagnosis of multiple sclerosis [START_REF] Okuda | Asymptomatic spinal cord lesions predict disease progression in radiologically isolated syndrome[END_REF], or when assessing cervical injuries using dynamic/kinematic MR imaging with positional changes [START_REF] Abitbol | Dynamic MRI of the spine[END_REF]. Many methods have been developed for both fully global rigid registration and fully local non-rigid registration separately [START_REF] Maintz | A survey of medical image registration[END_REF], but the literature on hybrid methods, allowing for adequate registration of the structures depending on the stiffness of their components, is still quite sparse.

The earliest work we know of is that of Little et al. [START_REF] Little | Deformations Incorporating rigid structures[END_REF], who showed how to incorporate rigid structures into a deformation field, using radial basis functions; this was later improved by others to make the field invertible and even diffeomorphic [START_REF] Arsigny | A Fast and Log-Euclidean Polyaffine Framework for Locally Linear Registration[END_REF][START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF][START_REF] Seiler | Geometry-Aware Multiscale Image Registration via OBBTree-Based Polyaffine Log-Demons[END_REF]. However, these methods require the user to specify which structures are rigid, which led to the development of semi-automated methods in which rigidity can be locally favoured/enforced through a regularisation term in the criterion to be minimised [START_REF] Rohlfing | Intensity-based non-rigid registration using adaptive multilevel free-form deformation with an incompressibility constraint[END_REF]. This idea was later improved to allow for this term to be adaptively tuned to the structures to register, through prior segmentation of the rigid parts or design of a stiffness map (typically computed from the image intensities; e.g. bones have high intensities in CT) [START_REF] Loeckx | Nonrigid image registration using free-form deformations with a local rigidity constraint[END_REF][START_REF] Ruan | Nonrigid registration using regularization that accomodates local tissue rigidity[END_REF][START_REF] Staring | A rigidity penalty term for nonrigid registration[END_REF]. Instead of segmenting rigid parts, it was also proposed to define several anchors, to which is attached an unknown polyaffine transformation, which can be subsequently estimated using a modified EM-ICP algorithm [START_REF] Taquet | Spatially Adaptive Log-Euclidean Polyaffine Registration Based on Sparse Matches[END_REF].

In this paper, we propose an iterative two-step method to compute a diffeomorphic non-rigid transformation between images of structures with rigid parts, without any user intervention or prior knowledge on the image intensities (to compute rigid parts or anchors). First we compute spatially sparse, locally optimal rigid transformations between the two images by adopting a new (as opposed to classical, translation-based) block matching strategy, made possible by the use of an efficient numerical optimiser (BOBYQA) (Sec. 2.1). The rationale behind this original strategy is our hope to recover both large rotations and subvoxel displacements. Then we derive a dense, regularised velocity field based on these local transformations using matrix logarithms and M-smoothing (Sec. 2.2). The floating image is then resampled and the two steps are iterated until convergence; the final diffeomorphic transformation is defined as the exponential of the accumulated velocity field. We finally compare our algorithm with the state-ofthe-art log-domain diffeomorphic demons method [START_REF] Vercauteren | Symmetric log-domain diffeomorphic registration: A demons-based approach[END_REF] on dynamic cervical and multiple sclerosis MR images (Sec. 3).

Material and Methods

To compute a diffeomorphism T between a reference image I and a floating image J, we iterate between two steps: computation of a sparse set of locally optimal rigid transformations using block matching between I and J • T l (Sec. 2.1) and computation of a dense velocity field δLT l computed from these locally estimated transformations (Sec. 2.2). Given that the transformation T is initially set to the identity (T 0 = Id), and that the initial velocity field is set to LT 0 = log T 0 = 0, the velocity field is then updated as LT l+1 = LT l + δLT l . This two-step algorithm stops at the iteration l when δLT l is close to 0, and the final diffeomorphism is computed as T = T l = exp(LT l ). The complete algorithm is outlined in Sec. 2.3. For the sake of clarity, we detail the two steps using the simpler notations I, J and δLT (Sec. 2.1 and 2.2).

Computing a sparse set of locally optimal rigid transformations

Classical block matching algorithm. In this approach, that we do not follow, one first defines a set of blocks in each image, before matching each block in the reference image I with the most similar block in the floating image J. Similarity is typically computed using a measure on the voxel intensities, such as the sum of squared differences or the squared correlation coefficient in monomodal problems, or the mutual information or the correlation ratio in multimodal problems. The most common approach to optimise the similarity measure (at least in medical image analysis) is to perform an exhaustive search of the block with the highest similarity in J, within a given neighbourhood of each block in I.

This strategy implicitly assumes that the local motion between the images can be well recovered by a discrete translation (i.e. defined on the discrete grid of the image). Subvoxel displacements and large rotations are thus likely to be missed. It is all the more true when registering piecewise rigid structures, because in this case there exists no single, global rigid movement, that could be corrected before non-rigid registration.

Modified block matching algorithm. Recent advances in nonlinear optimisation allow for testing another strategy. We first define a set of blocks in I (as in the standard strategy), and then we propose to directly compute the rigid transformation best superposing each of these blocks with J, using a similarity measure on the voxel intensities. As opposed to the standard, discrete translation-based strategy, the computation of the similarity measure for a given block in I and a given tested transformation implies resampling to build the block in J. Given that the solution space is no longer finite, this leads to a potentially much more computationally expensive algorithm.

We propose to use the recent BOBYQA algorithm [START_REF] Powell | The BOBYQA algorithm for bound constrained optimization without derivatives[END_REF] to implement this idea. In essence, BOBYQA is a derivative-free, trust-region method which uses successive approximations of the similarity measure by quadratic functions, whose maxima can be computed analytically. It is very similar to the classical NEWUOA algorithm except that bounds must be specified on the variables. We thus end up with a set of blockwise-estimated optimal rigid transformations between I and J. In practice, however, we do not estimate a transformation for the blocks in I having a low variance σ 2 . The set of estimated transformations (R 1 , . . . , R m ) is thus spatially sparse, due to these missing transformations, and also due to the resolution of the grid of blocks in I, which is different from that of I. In addition, we weigh each estimated transformation R i with a weight w i set equal to the similarity measure; here we use the squared correlation coefficient, to be insensitive to local intensity changes, thus 0 ≤ w i ≤ 1.

Estimating a dense velocity field

The set of estimated transformations (R 1 , . . . , R m ) is spatially sparse, but is also noisy and likely to contain outliers (due to the noise in the images to be registered and the potential errors in local registrations). How to estimate a dense (n = card(I)) and smooth velocity field δLT from (R 1 , . . . , R m )? We propose to use the logarithms of these m transformations, defined in the space of 4 × 4 real matrices (M 4 (R)) restricted to those whose last row contains only zeros, and to estimate n intermediate matrices in the same space (that we name log S 1 , . . . , log S n by analogy) as the minimisers of a criterion C:

(log S 1 , . . . , log S n ) = arg min log S1,...,log Sn   n i=1 j∈Vi w j ρ(|| log S i -log R j || 2 )d(|v i -v j | 2 )   ,
where:

ρ : R → R + is a robust error norm, -||.|| is the Frobenius norm in M 4 (R), and |.| is the Euclidean norm in R 3 , -v j is the coordinate of the central voxel of the block where R j was estimated, -v i is the coordinate of the voxel where S i is to be computed, -V i is a neighbourhood around the position v i ; note that the sum over j ∈ V i must read: "the sum over all the points v j where a rigid transformation R j was estimated, and which are inside V i ", -w j is the weight defined in Sec. 2.1, -d(.) : R 3 → R + is a (spatial) error norm.

It must be clear that we do not estimate S i and then its logarithm; we do estimate log S i directly; we use this notation here only as a convention for the sake of simplicity. Solving this minimisation problem is known as local M-smoothing [START_REF] Mrazek | On robust estimation and smoothing with spatial and tonal kernels[END_REF], due to the use of a robust error norm and of spatial neighbourhoods to design C; C can be minimised using gradient descent, where each transformation can be estimated independently of the others. Using a particular adaptive, data-dependent step size leads to an easy-to-interpret update formula for each log S i [START_REF] Mrazek | On robust estimation and smoothing with spatial and tonal kernels[END_REF]:

log S k+1 i = j∈Vi w j ρ ′ (|| log S k i -log R j || 2 )d(|v i -v j | 2 ) log R j j∈Vi w j ρ ′ (|| log S k i -log R j || 2 )d(|v i -v j |
2 ) It can be seen from this formula that ρ ′ acts as a tonal kernel, while d acts as a spatial kernel. After convergence, each finally estimated log S i is a linear combination of the logarithms of the rigid transformations R j ; following Arsigny et al. [START_REF] Arsigny | A Fast and Log-Euclidean Polyaffine Framework for Locally Linear Registration[END_REF], we define the final dense velocity field δLT as δLT (v i ) = log(S i ).v i , ∀i = 1, . . . , n. In practice, we define ρ as the Welsch function, which leads to ρ ′ (a 2 ) = exp(-a 2 /2λ 2 ), and we define d as d(b 2 ) = exp(-b 2 /2θ 2 ); this leads to two similar expressions for the two kernels (with different bandwidths). V i is spherical with radius 2θ (to achieve an approximate 95% confidence interval for a Gaussian law). To initialise the gradient descent algorithm, S 0 i is computed as the solution of the update formula by setting ρ ′ (a 2 ) = 1 (i.e. no tonal kernel).

Complete algorithm

The final estimated transformation is a diffeomorphism [START_REF] Arsigny | A Fast and Log-Euclidean Polyaffine Framework for Locally Linear Registration[END_REF]. We perform all the update calculations on the velocity field, whose exponential is required only once per iteration to resample the floating image. This is the same approximation as that done by Vercauteren et al., who showed experimentally that exp(LT l-1 ) • exp(δLT ) could be approximated by exp(LT l-1 + δLT ) for a small enough velocity field δLT [START_REF] Vercauteren | Symmetric log-domain diffeomorphic registration: A demons-based approach[END_REF].

Implementation details

For the block matching: size of the blocks: 7 voxels; grid step size: 3 voxels; minimal intensity variance in the blocks: 1/4 of the maximum squared intensity; search radiuses within BOBYQA: 2 voxels (translation) and 5 degrees (rotation). Estimate local rigid transformations using block matching (Sec. 2.1): R = (R1, . . . , Rm) ← block-matching(I, J • T l-1 ) 5:

Interpolate a dense velocity field using M-smoothing (Sec. 2.2): δLT ← M-smooth(R) 6:

Increment the velocity field: LT l = LT l-1 + δLT 7:

Regularise (elastic-like) the velocity field: LT l ← Gν * LT l 8:

Compute T l = exp(LT l ) to resample J 9:

until δLT is sufficiently small

For the M-smoothing: kernel bandwidths:

λ 2 = med j =h || log R j -log R h || 2 /2 (
tonal), θ = 4 voxels (spatial); ν = 4 voxels in the elastic-like regularisation. We use a 3-level multiresolution strategy, and the resampling of the floating image is done using trilinear interpolation. Run-time of the algorithm (dual core Xeon 3.0 GHz PC): about 6 min (vs 2 min for the log-domain diffeomorphic demons [START_REF] Vercauteren | Symmetric log-domain diffeomorphic registration: A demons-based approach[END_REF]).

Validation and Results

We propose to assess our algorithm quantitatively on ten patients with traumatic cervical cord injury, who got dynamic cervical MRI (T2-w, size 384 × 384 × 14, voxel size 0.8 × 0.8 × 3 mm 3 ) with two different positions each: either flexion/neutral, or extension/neutral [START_REF] Abitbol | Dynamic MRI of the spine[END_REF]. For each patient, we manually defined landmarks on the cervical/thoracic vertebrae C1-C3-C6-T1 (and T4 when visible), the pontomedullary junction, and the gnathion (lower border of the mandible) on each of the two MRI. We considered the neutral position as the reference image in the extension/neutral setting, and the flexion as the reference image in the flexion/neutral setting. For a given patient, the registration accuracy was evaluated as the root mean square error (RMSE) computed over the homologous landmarks after registration using four different methods: global rigid registration (M 1 ) [START_REF] Ourselin | Block Matching: A General Framework to Improve Robustness of Rigid Registration of Medical Images[END_REF], log-domain diffeomorphic demons (M 2 ) [START_REF] Vercauteren | Symmetric log-domain diffeomorphic registration: A demons-based approach[END_REF], our algorithm (M 3 ), and M 2 initialised using M 3 (M 4 ); both M 2 and M 3 are initialised using M 1 .

We also assessed our algorithm visually on two patients with MS lesions in the spinal cord, and two other patients with tumours in the spinal cord, with two time points each (T1-w, size 256 × 256 × 64, voxel size 1 × 1 × 1 mm 3 ). The box-and-whisker plot in Fig. 1 computed from the ten patients shows our algorithm (M 3 ) to significantly outperform both M 1 (paired t-test: p = 3×10 -4 ) and M 2 (paired t-test: p = 2 × 10 -3 ), with much smaller error and much smaller error dispersal. M 3 is also slightly better than M 4 (paired t-test: p = 5 × 10 -2 ). This suggests that the log-domain diffeomorphic demons performs worse than our algorithm even when properly initialised (using M 3 instead of M 1 ). The results of M 1 , M 2 and M 3 on one of the ten patients are shown in Fig. 2, and that of M 1 and M 3 are shown on one of the MS patients in Fig. 3.

Conclusion & Perspectives

It appears that our strategy for non-rigid registration, based on the computation of locally optimal rigid transformations in the first place, allows us to recover displacements and deformations of piecewise rigid structures (as seen e.g. in dynamic cervical MRI) much better than standard methods which are implicitly based on locally optimal translations, such as the log-domain diffeomorphic demons algorithm. This original strategy was made possible by (i) the use of an up-to-date very efficient optimiser and (ii) the design of a specific regularising procedure on the (sparse) set of locally estimated rigid transformations, based on robust estimation techniques. A future line of research could be to combine our regularisation technique with those previously proposed in this context [START_REF] Loeckx | Nonrigid image registration using free-form deformations with a local rigidity constraint[END_REF][START_REF] Ruan | Nonrigid registration using regularization that accomodates local tissue rigidity[END_REF][START_REF] Staring | A rigidity penalty term for nonrigid registration[END_REF].

Our intuition is also that our algorithm could perform very well in more general problems, without necessarily rigid structures involved, and on other image modalities; we will evaluate this in a near future.

Algorithm

1 :

 1 Initialize T to identity: T 0 ← Id = exp(LT 0 ) and the velocity field to 0: LT 0 ← 0 2: for each pyramid level of the multiresolution scheme, do 3: repeat 4:

Fig. 1 .Fig. 2 .

 12 Fig. 1. Quantitative evaluation of registration accuracy. Box-andwhisker plot of registration errors (unit: millimetres) for the 4 compared methods (M1 to M4). The errors and the error dispersal are much smaller for M3 compared to M1 and M2; M4 is also slightly worse than M3, which suggests that M4 actually degrades the results compared to M3 when initialised with M3.

Fig. 3 .

 3 Fig. 3. Registration results on a patient with MS lesions in the spinal cord. (a) patient image at time point t0; (b) patient image at time point t1 registered with M1 and (e) with M3; (d,g) zoom on (b,e); (h) zoom on (a); (c,f) zoom on the difference (registered minus reference) image. Note that in this case M2 (not displayed) performs as well as M3. These snapshots, and in particular the zoomed difference images, visually show that M3 gives a better result than a simple global rigid registration.

Acknowledgments. We thank Élise Bannier, Jean-Yves Gauvrit and Aurore Esquevin (CHU, University Hospital of Rennes, F-35043 Rennes, France) for discussions on the clinical protocol and for acquiring the MR images. We also thank Charles Kervrann (INRIA, SERPICO Project, F-35042 Rennes, France) for last-minute and life-saving reading suggestions.


Maxime Taquet

Benoit Scherrer

Olivier Commowick

Jurriaan Peters

Mustafa Sahin

Benoît Macq

Simon K Warfield

Benoît Scherrer

Registration and Analysis of White Matter Group Differences with a Multi-Fiber Model

Keywords: Diffusion Imaging, Multi-Fiber Models, Registration, White Matter

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

Introduction

Diffusion magnetic resonance imaging offers the ability to investigate in vivo the white matter microstructure. The representation of the signal by diffusion tensor images (DTI) has proven useful for population analysis in two ways [START_REF] Zhang | High-dimensional spatial normalization of diffusion tensor images improves the detection of white matter differences: an example study using amyotrophic lateral sclerosis[END_REF]. First, scalar features extracted from DTI, such as the fractional anisotropy (FA), may indicate the presence of brain diseases. Second, the use of DTI in registration improves the detection of morphometric differences, compared to scalar images.

The single tensor diffusion model has, however, proven inaccurate for two main reasons. First, it cannot represent the signal arising from multiple fibers with heterogeneous orientations in one voxel. Second, it does not account for the non-monoexponential decay observed when imaging at high b-values. Novel models addressing one or both of these issues have been introduced [START_REF] Minati | Physical foundations, models, and methods of diffusion magnetic resonance imaging of the brain: A review[END_REF] : Qball imaging, spherical deconvolution, 4th order tensors, DOT, and others. Most of them focus on describing the general shape of the diffusion at each voxel.

In contrast, mixture models represent each fiber bundle independently, keeping the interpretability of single fiber models while accounting for crossing fibers. Therefore, scalar quantities such as the fractional anisotropy (FA) can still be computed for each fiber independently. This property makes them very attractive for population analysis.

While the literature on the registration of complex diffusion models is growing (e.g., [START_REF] Barmpoutis | Registration of high angular resolution diffusion mri images using 4th order tensors[END_REF][START_REF] Yap | Sphere: Spherical harmonic elastic registration of hardi data[END_REF]), no method has been developed to register mixture models. This lack of a registration method limits the use of mixture models in population analysis despite their attractiveness. This issue has been previously reported, and has incited researches to register the raw diffusion weighted images instead [START_REF] Bergmann | Two-tensor fiber tractography[END_REF].

The remaining of this paper is organized as follows. Section 2 introduces the diffusion mixture models. Section 3 presents a method to compute weighted average of mixture models. Section 4 develops a similarity metric for diffusion mixture images. Section 5 presents the integration of the developed methods in a registration algorithm and analyzes its complexity. Section 6 presents experimental results on a dataset of 45 subjects. Finally, Section 7 concludes.

Diffusion Mixtures

The basic idea behind multi-fiber models is to fit a single fiber model to each of the fiber bundles present in the voxel. If S i (x) is a suitable model to represent the diffusion process in a single fiber, then,

S(x) = N i=1 f i S i (x) (1)
is a multi-fiber model for N crossing fibers with relative volumetric occupancy given by f i . The assumption behind these models is that the exchange of water molecules between populations of fibers is negligible during the diffusion time [START_REF] Minati | Physical foundations, models, and methods of diffusion magnetic resonance imaging of the brain: A review[END_REF].

The simplest multi-fiber model is the multi-tensor model in which S i = S 0 e -bg T Dig . More complex multi-fiber models have later been introduced [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF]. Potentially, any single fiber model can be extended to a multi-fiber model by means of mixtures. One such model, the biexponential decay model [START_REF] Minati | Physical foundations, models, and methods of diffusion magnetic resonance imaging of the brain: A review[END_REF], represents each fiber bundle by a Gaussian mixture to capture the non-monoexponential decay of the signal. The corresponding multi-fiber model would be a mixture of Gaussian mixtures which is itself a Gaussian mixture. A natural parameterization of diffusion Gaussian mixtures is the set of pairs (fraction, covariance matrix), that we write: {(f 1 , Σ 1 ), ..., (f N , Σ N )}. Alternatively, to connect with the tensor formalism, the inverse of the covariance matrix, D i = Σ -1 i , can be used.

Weighted Combination of Mixtures

Computing weighted combinations of voxel values is at the basics of interpolation (the value in one location is the weighted combination of the values in the neighborhood), smoothing (the value at a grid voxel is replaced by a weighted combination of the values in a neighborhood) and atlas construction (the value at one voxel is the average of the values in the aligned subjects').

Gaussian mixture simplification (GMS) was introduced to efficiently compute weighted combinations of diffusion mixture [START_REF] Taquet | Interpolating multi-fiber models by gaussian mixture simplification[END_REF]. In this section, we underline the important aspects of this method. The idea behind GMS is that computing weighted combinations of mixtures would be trivial if the number of components of the result could be arbitrarily large. Indeed, the linear combination of K mixture models with N components is a mixture models with M = KN components:

M C = K k=1 w k M k = K k=1 w k N j=1 f k j S k j (x) ≡ M i=1 g i S i (x). (2)
We refer to this mixture as the complete mixture. GMS optimizes the parameters of a simplified mixture M S = N j=1 h j R j (x) with N < M components that best approximates M C . The energy function to be minimized is the cumulative differential entropy (the reference to the location x is omitted for clarity):

D(M C , M S ) = N j=1 i:πi=j g i D(S i ||R j ) = N j=1 i:πi=j g i S i (g) log S i (g) R j (g) dg, (3)
where g is the gradient vector and where latent variables π i cluster the components of the complete mixture S i in N clusters each represented by a single component of the simplified mixture, R j ; π i = j means that S i is best represented by R j . Following the recent developments in probabilistic clustering, an EM scheme is used to minimize (3). Banerjee et al showed that both the E-step and the M-step can be solved in closed form for mixtures of exponential distributions [START_REF] Banerjee | Clustering with bregman divergences[END_REF]. For Gaussian mixtures, the E-step consists in optimizing the latent variables π i by computing the Burg matrix divergence between the covariance matrices of each component of M C (Σ S i ) and each component of M S (Σ R j ):

π i = arg min j B Σ S i , Σ R j = arg min j Tr Σ S i Σ R j -1 -log Σ S i Σ R j -1 . (4)
As for the M-step, it sums up to calculating:

Σ R j = i:πi=j f i Σ S i i:πi=j f i and h j = i:πi=j f i . (5)
Alternating ( 4) and ( 5) until convergence provides the parameters (h j and Σ R j ) of the resulting mixture. A log-Euclidean version of this interpolation scheme is obtained by replacing all covariance matrices by their logarithm.

Generalized Correlation Coefficient for Mixtures

The correlation coefficient, invariant under linear transformations of the voxel intensities, is widely used in mono-modal image registration. The inter-subject variability of diffusivity values motivates the introduction of a generalized correlation coefficient, invariant under these differences. In DTI, this variability has been reported and partially accounted for in some registration methods [START_REF] Zhang | High-dimensional spatial normalization of diffusion tensor images improves the detection of white matter differences: an example study using amyotrophic lateral sclerosis[END_REF]. The correlation coefficient between blocks F and G is defined as the scalar product of the normalized blocks:

ρ(F, G) = F -µ F ||F -µ F || , G -µ G ||G -µ G || ,
where µ F is the mean of the image values in the block. It is invariant if F (and/or G) is replaced by aF + b. It has been generalized to vector images by redefining the means µ F and µ G as the projection of the block onto a constant block T [START_REF] Ruiz-Alzola | Nonrigid registration of 3d tensor medical data[END_REF]:

F -µ F = F -F, T T ||T || 2 .
The corresponding generalized correlation coefficient is invariant if F is replaced by aF + bT where T is now any constant vector block. The definition of a scalar product between two blocks of mixture models seems impractical if not impossible. We therefore further generalize the correlation coefficient by substituting the inner product by a more general scalar mapping, m(M f , M g ):

ρ(M f , M g ) = m M f -m(M f , T )T n m (M f -m(M f , T )T ) , M g -m(M g , T )T n m (M g -m(M g , T )T )
,

where n m (M) 2 = m(M, M) is a generalization of the norm. This definition does not guarantee the invariance property of the metric for any scalar mapping. One can show that the invariance is preserved as long as the scalar mapping is linear with respect to the constant block T :

m(aM f + bT , T ) = a m(M f , T ) + b m(T , T ). ( 6
)
To preserve the interpretability of ρ as a similarity metric, it needs to be symmetric, equal to one in case of perfect match and lower than one in any other case. These constraints on ρ translate into the following constraints on m:

m(M f , M g ) = m(M g , M f ) (7)
n m (aM f ) = a n m (M f ) (8) |m(M f , M g )| ≤ n m (M f )n m (M g ). ( 9
)
The latter is a generalized form of the Cauchy-Schwartz inequality for inner products. Conditions [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Taquet | Interpolating multi-fiber models by gaussian mixture simplification[END_REF][START_REF] Banerjee | Clustering with bregman divergences[END_REF][START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF], the choice of T and the definitions of the addition M + T and multiplication by a scalar aM, stand together as a model to define a correlation coefficients in potentially any space. For DTI, if T is an isotropic tensor block (T (x) = DI 3×3 ), m is the log-Euclidean scalar product, and the log-Euclidean algebra is used, then ρ is invariant under linear transformations of the eigenvalues in the log-domain [START_REF] Taquet | A generalized correlation coefficient: Application to dti and multi-fiber dti[END_REF]. For multi-tensor images, we fix T (x) = ( 1 N , DI 3×3 ), ..., ( 1 N , DI 3×3 ) , and we define the addition of T , and the multiplication by a scalar component-wise in the log-domain. The scalar mapping m(M f , M g ) is defined by pairing the tensors in each voxel to maximize the linear combination of pairwise scalar products. Let M f (x) = {(f 1 , F 1 ), ..., (f N , F N )} and M g (x) = {(g 1 , G 1 ), ..., (g N , G N )} defined on a domain Ω, we have:

m(M f , M g ) = x∈Ω max π N i=1 f i g π(i) F i , G π(i) ,
where π is a pairing function associating one tensor of M g to each tensor of M f . This scalar mapping satisfies conditions [START_REF] Assaf | Composite hindered and restricted model of diffusion (charmed) mr imaging of the human brain[END_REF][START_REF] Taquet | Interpolating multi-fiber models by gaussian mixture simplification[END_REF][START_REF] Banerjee | Clustering with bregman divergences[END_REF][START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF]. Interestingly, the resulting generalized correlation coefficient is invariant under any global (within the block) linear transformation of all eigenvalues in the log-domain. This similarity metric is therefore robust to the inter-subject variability of diffusivities.

Implementation and Complexity

The developed methods were integrated in the efficient block matching registration algorithm described in [START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF]. The parameters used are the following: 4 pyramid level, 10 iterations per level, block size: 5 × 5 × 5, outlier removal rate: 20%. The implementation was multi-threaded. On a 8 core workstation, with 220×220×176 two-fiber images, the entire registration takes 1.5 hour. All weighted combinations were computed until complete convergence of the soft clustering. The average number of iterations required for that convergence is 4.

Results

The registration was applied to a clinical dataset of 45 subjects, 13 controls and 32 patients with tuberous sclerosis complex (TSC), a rare genetic disease associated with impaired white matter integrity. A DTI and a multi-tensor model with three components (one isotropic and two anisotropic) were reconstructed for each subject [START_REF] Scherrer | Toward an accurate multi-fiber assessment strategy for clinical practice[END_REF].

Validation

An alternative to the method presented in this paper would be to select one of the two tensors in each voxel (e.g. the one with the highest FA) and to perform single tensor registration on this image. Here, we validate that our method works better than this simple alternative. The quality of the alignment is assessed by the sum of square differences of each eigenvalue after alignment of control subjects. Indeed, while the diffusivities can significantly differ in diseased brain, they are approximately equal for healthy subjects. We performed 26 randomly chosen registrations with four levels of regularization, totalizing 104 registrations.

In each voxel, the eigenvalues were averaged between the two anisotropic components (weighted by their fractions). Results show that multi-fiber registration performs significantly better than single tensor registration (Fig. 1).

Atlas Construction

An atlas was constructed using our registration and alternating three steps: aligning all subjects to the current atlas (initially a randomly chosen subject), averaging the aligned subjects (using the weighted combination of mixtures), applying the mean inverse field to the resulting average [START_REF] Guimond | Average brain models: A convergence study[END_REF]. This atlas remarkably shows areas where multiple fibers are consistently present in all subjects (Fig. 2).

Morphometry

The clinical hypothesis according to which there is substantial white matter shrinkage in TSC subjects was tested by performing a one-tailed two sample ttest on the log-Jacobian of the deformation fields [START_REF] Ashburner | Identifying global anatomical differences: deformation-based morphometry[END_REF]. The subject classes were then randomly permuted 4000 times to assess the null distribution of extreme t-scores. The entire process was repeated with single tensor images. As a result, multi-fiber registration reveals more white matter differences (> 3800 significant voxels) than single tensor (< 1000 voxels) (Fig 2). The entire process was then repeated with DT-REFinD, a state-of-the art DTI registration algorithm [START_REF] Yeo | Dt-refind: Diffusion tensor registration with exact finite-strain differential[END_REF], to test whether the improved detection of differences is truly due to the knowledge brought by multi-fiber models. Again, DT-REFinD did not capture all the differences detected by multi-fiber registration (< 1300 voxels) (Fig 2).

Tract-Based Statistics

Some structural subnetworks are believed to be impaired in TSC patients. To test this hypothesis, we analyzed the FA profile along the median tract of the arcuate fasciculus, generated on the atlas by a probabilistic tractography algorithm [START_REF] Bergmann | Two-tensor fiber tractography[END_REF](Fig. 3(a)). A one-tailed two-sample t-test was performed at every location. A threshold t 0 was then set to the t-statistics and the length of the contiguous supra-threshold segments were recorded. The null distribution of these lengths was assessed by randomly permuting the subjects classes 4000 times. The operation was repeated for a wide range of thresholds (1.5 ≤ t 0 ≤ 4.5) to estimate the robustness of our findings. For t 0 = 2.7 (p 0 = 0.01), with our multifiber registration, four significant clusters, together representing 15% of the tract were detected, indicating a strong impairment of this subnetwork in TSC patients (Fig. 3(b)). These findings were robust to the choice of t 0 for any t 0 ≤ 3.4.

In contrast, single tensor registration only revealed one cluster representing 5% of the fiber, which was not robust outside the range 1.9 ≤ t 0 ≤ 2.8.

Conclusions

This paper introduced a registration and atlas construction method to align multi-fiber models. A proper interpolation method and a robust similarity metric were presented. Results in both morphometry and tract-based statistics demonstrated that multi-fiber registration reveals more group differences than DTI registration. We therefore believe that this registration method opens new doors to understanding brain disorders based on multi-fiber models.

Fig. 1 .

 1 Fig.1. Comparison of the single tensor and multi-fiber registration in terms of the SSD between eigenvalues after alignment, for different regularization parameter values[START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF]. Multi-fiber registration significantly improves the quality of the registration.

Fig. 2 .

 2 Fig. 2. (top)The two-tensor atlas built by means of the developed registration method reveals crossing pathways common to all anatomies. (bottom) White matter volume shrinkage in tuberous sclerosis represented by the p-value maps. Multi-fiber registration reveals more differences than single tensor registration and DT-REFinD[START_REF] Yeo | Dt-refind: Diffusion tensor registration with exact finite-strain differential[END_REF]

Fig. 3 .

 3 Fig. 3. (a) Arcuate fasciculus, a set of fibers involved in language, on which tract based statistics was performed, (b) The FA profile in TSC patients shows significantly disrupted white matter fascicules in different clusters, indicated by the stars.

Acknowledgments MT thanks the FNRS, BAEF and WBI for their financial support. This work was supported in part by NIH grants R01 EB008015,R01 LM010033, R01 EB013248, and P30 HD018655 and by the Boston Children's Hospital Translational Research Program.


Philippe Meton
email: philippe.meton@lss.supelec.fr

Andrea Cozza
email: andrea.cozza@supelec.fr

Florian Monsef

Marc Lambert

Jean-Christophe Joly

Pierre Bruguière

J.-C Joly

Numerical analysis of pulsed local plane-wave generation in a TREC

Keywords: component; wavefront synthesis, test facilities, time reversal, TREC, electromagnetic pulses

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

INTRODUCTION

The ability to generate deterministic wavefronts is fundamental in any test facility aiming at characterizing how a device under test (DUT) responds to electromagnetic energy. This need ranges from EMC applications to radar imaging, passing through antenna tests. The common point in all of these scenarios is the presence of locally-planar wavefronts, with features such as polarization, direction of arrival, bandwidth, etc., that are preliminary defined. The simplest approach implies the use of anechoic environments (anechoic chambers, OATS, TEM cells, etc.) where a source (or a collection) generate a far-field distribution appearing as locally planar in the proximity of the DUT. Though of simple interpretation, the implementation of this technique can be quite cumbersome when a large number of incidence angles are required, particularly in the case of complete 3D characterization, because of the need of complex mechanical displacements. An alternative approach, using fixed sources, but without involving distributed multiple sources is the TREC.

II. THE TIME-REVERSAL ELECTROMAGNETIC CHAMBER

The TREC is a test facility based on the properties of timereversed wavefronts and the propagation of waves in diffusive media [START_REF] Moussa | Experimental demonstration of directive pulsed wavefront generation in reverberation chambers[END_REF][START_REF] Cozza | Emulating an Anechoic Environment in a Wave-Diffusive Medium through an Extended Time-Reversal Approach[END_REF]. In its original form it consisted of a reverberating cavity, driven at frequencies where it can support a diffuse field configuration, where the electromagnetic field can be well represented by a random incoherent field distribution, totally depolarized. Time-reversed signals naturally lead to spatial coherence functions of the field. They thus provide the ability to generate deterministic wavefronts, with a spatial resolution limited by the average spatial-coherence cells of the medium [START_REF] Wolf | Introduction to the Theory of Coherence and Polarization of Light[END_REF], i.e., about half a wavelength. Target wavefronts can thus be defined for a virtual source, and time-reversed to converge onto its phase-center. The main interest of this procedure is that the use of an isotropic diffusive medium implies that any direction of arrival can be generated without needing the line-of-sight conditions required by free-space like environments.

Results from numerical simulations will be provided in this respect, showing how the physical dimensions of the TREC cavity and the bandwidth of the test signals impose a lower bound to the phase precision of the planar approximation, as well as its amplitude uniformity. A scalar electric field configuration will be considered, by means of a quasi-2D cavity. Conclusions will be given with respect to the design of DUT testing in a TREC.

Figure 1 .

 1 Figure 1. An example of pulsed wavefront generated by a TREC, for a central frequency of 1 GHz and a bandwidth of 265 MHz.The dominant parts of the wavefronts are the incoming and outgoing waves crossing a circular array of field probes.


Anne-Marie Kermarrec
email: anne-marie.kermarrec@irisa.fr

François Taïani
email: f.taiani@lancaster.ac.uk

Constellation: Programming decentralised social networks

The past decade has witnessed a dramatic shift in the way the Web is used. The internet has entered our homes, our factories, and our lives like never before. User-generated content services (Flickr, Youtube, Delicious) and social networks (Twitter, FaceBook) have grown exponentially. The digital systems we live with are now composed of hundreds of millions of computing devices, of as many users, and of Terabytes of data, that is dynamically produced, shared, disseminated and searched globally. This data represents a fantastic potential to leverage information about each and every user (their circles of friends, their interests, their activities, the content they generate), and use it to provide better, more personalised, on-line services.

Although traditional social networks have been extremely successful, they today only encompass a tiny part of these many sources of information about users. It is not quite clear whether they should encompass more: Twitter has been known to experience regular outages because of the sheer load of data it must process. Facebook is regularly criticised out of privacy concerns.

These growing difficulties have prompted work on novel and fully decentralised alternatives to traditional social networks [START_REF] Bertier | The gossple anonymous social network[END_REF][START_REF] Pujol | The little engine(s) that could: scaling online social networks[END_REF]. Where traditional social engines fail to provide information that escape them, these decentralised data systems seek the information where it ultimately is: at the user.

To implement these decentralised social networks, gossip protocols appear as a natural solution, as they intrinsically tend to be highly resilient, efficient, and scalable. Existing gossip-based social networks have however so far been limited to simple use-cases, that typically take a uniform view of users, peers, and the data they hold [START_REF] Leroy | Cold start link prediction[END_REF][START_REF] Bai | Gossiping personalized queries[END_REF]. Similarly, the various mechanisms developed for search, recommendation, and personalisation in these systems currently only exist as independent standalone solutions, and lack a clear framework to allow their reuse, and composition, a key step to allow the incremental construction of more complex systems.

To progress to full-fledged decentralised social networks, we posit in this paper that we now need to move to gossipbased social applications that can simultaneously cater for different types of data and services. To design, and evaluate these new approaches in a modular and incremental manner, we further argue that we need specific and dedicated programming technologies. To help in this task, this paper sketches the main ingredients of a new programming language, Constellation, that seeks to simplify the realisation and experimentation with modular social gossip-based applications. Constellation is based on two central observations: (i) future decentralised social applications will need to handle heterogeneous forms of data and self-organisation, and (ii) these applications will need to better support modularity and incremental design, by providing appropriate

Decentralised Social Networks

We first provide some background on decentralised social networks, before presenting the two main mechanisms of Constellation (Section 3.1).

System model

A traditional decentralised social networks uses a peerto-peer architecture in which each peer (a smart phone, a tablet, a laptop) is associated with a user (Fig. 1). (We will discuss cases later in which a peer might also be associated to a web-site, a venue, a storage node.) Peers can connect to each other using point-to-point networking, but they only have a partial and very limited view of the rest of the system: Typically a small-size neighbourhood of other peers.

Each peer maintains some data (or profile) about its user, for instance her friends, or the set of items tagged in a collaborative tagging service (delicious), or the set of the locations visited in a geolocated social networks (foursquare). Peers use this data to organise themselves in a distributed overlay so that similar peers end up connected together. E.g. in Fig. 1, Alice and Bob have similar location profiles, and have therefore been selected to be neighbours of each other (using a mechanism we detail just below).

The resulting overlay can then be exploited to propose a range of personalised services such as search, recommendation, and query extension. The intuition behind this is that if Alice and Bob have similar profiles, then Alice's data can be leveraged to provide a better service to Bob. Here for instance, Alice might know of locations of interest to Bob (here the snowy house), and reciprocally (the circus tent).

Constructing the overlay

To construct its similarity overlay, a gossip-based social network uses a two-layer structure (Fig. 2) [START_REF] Jelasity | T-man: Gossip-based overlay topology management[END_REF][START_REF] Voulgaris | Epidemic-style management of semantic overlays for content-based searching[END_REF]. Both layer maintain an overlay, in which peers have a fixed list of neighbours. In Fig. 2, Alice is connected to Bob, Carl, and Dave in the bottom RPS layer, and to Carl and Bob in the upper layer (clustering). Periodically, each peer selects a user from its view and exchanges information about its neighbours with the final goal of converging to an optimal topology of users with 'similar' profiles in the top layer.

More precisely, the bottom RPS periodically provides each with a random sample of the rest of the network and thus guaranties the convergence of the second layer (clustering), while making the overall system highly resilient against churn and partitions. This is achieved by having peers exchange and shuffle their neighbours list in periodic gossip rounds to maximise the randomness of the RPS graph over time [START_REF] Jelasity | Gossip-based peer sampling[END_REF]. Alice might for instance request Carl's list of RPS neighbours, and randomly decide to replace Dave by Ellie (received from Carl) in her RPS view.

The clustering layer sits on top of the RPS layer, and implements a local greedy optimisation procedure that leverages both neighbours returned by the RPS, and current neighbours from the clustering views [START_REF] Jelasity | T-man: Gossip-based overlay topology management[END_REF][START_REF] Voulgaris | Epidemic-style management of semantic overlays for content-based searching[END_REF]. A peer (say Alice in Fig. 2) will periodically update its list of similar neighbours with new neighbours found to be more 'similar' to her in the RPS layer. This guarantees convergence under stable conditions, but can be particularly slow in large systems. This mechanism therefore is complemented by a swap mechanism in the clustering layer (Fig. 3), whereby two neighbouring peers (here Alice and Bob) exchange their neighbours lists (Step 1), and seek to construct a better neighbourhood based on the other peer's information (Step 2).

In Fig. 3 for instance, Alice has a thing for hearts, and prefers like-minded people. Bob, on the other hand, has a minor interest in hearts, but is much more interested in diamonds. When Alice sends to Bob her current list of neighbours, and Bob sends Alice his (Step 1), both discover new potential neighbours closer to their interests. Alice thus drops Ellie for Carl, and Bob drops Alice for Ellie.

Constellation

Design goals

Constellation's two main design goals are simplicity and composability. This translates into mechanisms and programming constructs that can very concisely represent gossip-based social networks in which different types of peer and clustering approach coexist (Sec. 3.2), at different levels of composition (Sec. 3.3).

Heterogeneous self-organisation

The peer-to-peer system that results from the approach presented in Section 2 is similar to a large set of physical particles (the peers) submitted to a uniform law of attraction (the similarity measure). In this physical metaphor, the RPS layer plays the role as thermal excitation (as in simulated annealing), while the clustering layer provides a local gradient approach to converge to an optimal topology (i.e. one that maximise similarity between nodes).

In Constellation, that kind of uniform clustering sys- tem is described by indicating (i) which data each node holds (e.g. interests below, a simple list of strings denoting the tags used by a user); and (ii) which similarity measure to apply between nodes (with the keyword clusters; here counting the overlap in tags between two users).

RPS

User { List<String> interests clusters potentialFriends { |u| overlap( this.interests, u.interests ) }}

Social applications may however involve different types of node (e.g. users and websites) and different similarities between these nodes, which a uniform approach cannot capture. To help program these more advanced systems, Constellation allows the declaration of multiple similarities between different types of nodes. For instance, in Fig. 4, in-terestingSites links each user with a list of websites this user might be interested in. With an RPS layer adapted to return both a random sample of both User and Website nodes, each user will eventually converge to a list of websites most likely to interest her. A complicating factor here comes from the fact that websites have no neighbourhood, which prevents the direct use of a gradient optimisation in the clustering layer (cf. Fig. 2). To work around this problem, Constellation provides the using keyword (last line in Fig. 4), which indicates here that user nodes should use their list of potential friends to optimise their list of available websites (a form of convergence piggybacking).

Virtualisation

The use of multiple similarity and heterogeneous peers described offer a world in which different 'attraction laws' (different similarities) can act on different particles (users, websites). In this model, all nodes are on the same level and usually correspond to a physical machine in a peer-topeer system. Some social applications, however, can benefit from treating entities that are not associated with individual machines as 'virtually distributed', and thus provide an additional mechanism to compose data and topologies. This applies for instance to the search queries made by a user (in a decentralised search system), to the videos stored by a storage node (in a video storage system), or to the items tagged by a contributor (in a collaborative tagging system).

To support these systems, Constellation allows nodes (called containers) to host "virtual nodes" (vnodes for short), i.e. entities which can maintain a similarity neighbourhood as nodes do (Fig. 5), but are not directly associated with a physical device. The neighbours of a vnode might belong to the same container node, or to distant containers. To allow vnodes to cluster themselves with other similar vnodes, Constellation includes two mechanisms: vertical data-aggregation, and virtualised clustering.

Vertical data aggregation allows a container node to access the data managed by its virtual nodes in the form of a multiset (a 'bag') of data items (red shaded arrow on Fig. 5). For instance, if a node User contains Query vnodes in a variable queries (Fig. 6), a user can access a collection of all her queries with this.queries. If each query maintains in turn a list of tags (as in Fig. 6), all tags contained in the queries of a user can be aggregated using this.queries.tags.

Using vertical aggregation, a Constellation program can define the neighbourhood of a container in terms of its vnode content (green solid arrow on Fig. 5). E.g. in Fig. 6 the similarity relationship similarUsers uses the list of tags of the queries made by a user to find similar users.

Virtualised clustering. Although vertical aggregation can be used on its own, its main use is to couple the similarity neighbourhood of vnodes with that of their containing node. This is useful semantically as vnodes (e.g. queries within users) are often contextually linked to their containing node: Similar queries from similar users are more likely to call for related answers, than similar queries from dissimilar users. Another reason pertains to efficiency: coupling the clustering of vnodes to that of containers allows for efficient realisations that do not overburden individual machines.

The general mechanism of virtualised clustering is shown in Fig. 7 for the code of Fig. 6: a vnode (e.g. a query x) can only maintain links with vnodes that lay in the neighbourhood of its containing node (here Alice). This is as if for x, the rest of the system only consists in Alice and her neighbours Bob and Carl, and the queries (vnodes) they contain (u, v, w, y, and z).

The clustering scope of queries is declared with the set scope keyword (Fig. 6), which indicates which similarity neighbourhood (here similarUsers) to use to scope the clustering of queries. The result in Fig. 7 is a two-level clustering system: users gravitate towards other users hosting similar queries, and queries are linked so similar searches in a user's neighbourhood, in a way similar to that proposed in [START_REF] Bertier | The gossple anonymous social network[END_REF].

Outlook

Constellation offers two other features we have not discussed: the ability for nodes and vnodes to probe the data in their neighbourhood (horizontal data aggregation), and the capability of vnodes to migrate between containing nodes. Although functionally simple, these four mechanisms (virtualisation, virtualised clustering, data aggregation, and migration) can be combined arbitrarily, and we think offer a particularly attractive playing ground to experiment with novel gossip-based social algorithms.

Constellation suggests that social infrastructures might be approached in terms of bio-inspired systems, an aspect we would like to investigate. It might also provide paths to reason about the correctness a systems through static analysis, which we would like to explore in the future.

Figure 1 :

 1 Figure 1: A typical gossip-based social network mechanisms to compose data and topologies.

Figure 2 :

 2 Figure 2: Gossip-based distributed clustering

Figure 3 :Figure 4 :

 34 Figure 3: Peer-to-peer clustering mechanism

Figure 5 :Figure 6 :

 56 Figure 5: Vertical data aggregation

Figure 7 :

 7 Figure 7: Virtualised clustering of Alice's queries (for the code of Fig. 6)

Acknowledgement

This work has been partially supported by the ERC Starting Grant GOSSPLE number 204742.


Anh-Dung Phan

Nikolaj Bjørner

David Monniaux Verimag

Anatomy of Alternating Quantifier Satisfiability (Work in progress)

. We report on experiments on different variants of the auxiliary procedures. So far, there is an edge to applying SMT-TEST proposed in [10], while we found that a simpler approach which just eliminates quantified variables per round is almost as good. Both approaches offer drastic improvements to applying default quantifier elimination.

Introduction

Can formulas with nested quantifiers be checked effectively for satisfiability? Several algorithms exist in the context of Quantified Boolean Formulas that handle alternation of quantifiers [START_REF] Plaisted | A satisfiability procedure for quantified Boolean formulae[END_REF][START_REF] Biere | Resolve and Expand[END_REF]. They are specialized for eliminating variables over Booleans. An algorithm for alternating quantifier satisfiability was given in [START_REF] Monniaux | Quantifier Elimination by Lazy Model Enumeration[END_REF] for the case of linear arithmetic over the Reals. It integrates tightly an All-SMT loop and projection based on Fourier-Motzkin elimination or Chernikov projection. A question arises whether the ideas lift to other projection procedures. Also, are there reasonable alternatives to All-SMT and how do they compare? This ongoing work presents a generalized algorithm of that presented in [START_REF] Monniaux | Quantifier Elimination by Lazy Model Enumeration[END_REF] which abstracts the auxiliary procedures. We instantiate the generalization to projection functions based on virtual substitutions, i.e. substitution methods that replace quantifiers by disjunctions of bounded variables. The specialization is for Linear Integer Arithmetic based on Cooper's procedure and used for formulas that arise from Duration Calculus Model Checker (DCMC).

Linear Integer Arithmetic (LIA) or Presburger Arithmetic, introduced by Mojzaesz Presburger in 1929, is a first-order theory of integer which accepts + as its only operation. A classic example of representing some amount of money by 3-cent coins and 5-cent coins appears in LIA as follows:

∀z (z ≥ 8 → ∃x ∃y (3x + 5y = z))

After Presburger proved decidability of LIA [START_REF] Stansifer | Presburgerś Article on Integer Airthmetic: Remarks and Translation[END_REF], LIA attracted a lot of attention due to applications in different areas. Cooper's algorithm [START_REF] Cooper | Theorem proving in arithmetic without multiplication[END_REF] is a substitution-based decision procedure. The Omega Test is a projection-based decision procedure for LIA and employed in dependence analysis of compilers [START_REF] Pugh | The Omega test: a fast and practical integer programming algorithm for dependence analysis[END_REF]. A variant that integrates elements from both Cooper's method and the Omega Test is implemented in Z3 [START_REF] Bjørner | Linear Quantifier Elimination as an Abstract Decision Procedure[END_REF]. While Z3 can handle non-trivial LIA problems, applications from DCMC also expose limitations of using quantifier elimination alone. The time complexity of all procedures for Presburger Arithmetic is high. Let n denote the length of a LIA formula; running time of any decision procedure is at least 2 2 cn for some constant c > 0 [START_REF] Michael | Super-Exponential Complexity of Presburger Arithmetic[END_REF]. Moreover, Oppen proved a triply exponential upper bound 2 2 2 cn for worst-case running time of Cooper's algorithm [START_REF] Oppen | A 2 2 2 pn Upper Bound on the Complexity of Presburger Arithmetic[END_REF]. This paper is organized as follows. Section 2 presents the generalized algorithm with a few supporting procedures. Section 3 presents different methods for producing candidate formulas. Section 4 instantiates the algorithm with a concrete procedure for virtual substitutions. We discuss implementation details and benchmarks in Section 5 and Section 6 concludes the paper.

Alternating Quantifier Satisfiability

This section develops an algorithm that is an abstraction of the alternating quantifier satisfiability algorithm presented for Linear Real Arithmetic in [START_REF] Monniaux | Quantifier Elimination by Lazy Model Enumeration[END_REF]. The abstraction is formulated such that various quantifier elimination methods can be plugged in, including virtual substitutions.

Definitions

The algorithm being developed relies on two procedures for extrapolation and projection. We first describe the requirements for these procedures and discuss the main algorithm later.

Definition 1 (Extrapolant). Given two formulas A and B, a formula C is an extrapolant of A and B if the following conditions are satisfied:

A ∧ B is unsat then C = f alse A ∧ B is sat then A ∧ C is sat, ¬B ∧ C is unsat We use A, B to denote an extrapolant.
Extrapolation is typically understood as finding new data points outside a set of existing points. Intuitively, C has empty intersection with ¬B and non-empty intersection with A. There are many possible extrapolants for each pair of formulas, and the definition here does not specify how to compute an extrapolant. An example of a trivial extrapolant is described in the following procedure: when A ∧ B is satisfiable we can take B as an extrapolant, otherwise take false.

Definition 2 (Projection πx.(C|M )). Let M and C be quantifier-free formulas where variable x only occurs in C (x / ∈ F V (M ) where F V (M ) denotes the set of free variables in M ). Assume C ∧ M is satisfiable. A projection procedure πx.(C|M ) computes a quantifier-free formula satisfying the conditions:

1. F V (πx.(C|M )) ⊆ F V (C) \ {x} 2. πx.(C|M ) is sat 3. (M ∧ πx.(C|M )) → ∃x C
Similar to extrapolation we only gave the conditions that projection functions have to satisfy for the developing algorithm to be sound. There is a choice of algorithms for implementing πx.(C|M ). A possible way is to use virtual substitutions, where we derive a quantifier-free formulas by substituting variables by one or more disjunctions. Virtual substitutions will be presented with more details in Section 4.

Quantifier Test: algorithm QT

QT is defined as a recursive function in Algorithm 1. The four arguments of QT can be explained in the following way: C is a context formula which reflects collected information at

Algorithm 1: QT (i, C, x, M ) if C ∧ M i is unsat then return (f alse, M ) end if i = n then return ( C, M i , M ) end (C , M ) ← QT (i + 1, C, M i , x, M ); if C = f alse then return ( C, M i , M ) end M k ← M k , ∀ k = i; M i ← M i ∧ ¬(πx i .(C |M i )); return QT (i, C, x, M ); Algorithm 2: QE(x, F n ) M k ← true, ∀ k < n; M n ← F n ; C ← false; C ← true; while C = false do (C , M ) ← QT (1, ¬C, x, M ); C ← C ∨ C ; end return C;
current iteration; x and M are vectors of quantified variables and formulas respectively, and i is the index to access current elements in the above vectors. We use the notation M i to refer to the i-th element of vector M and imply the same notation for other vectors.

Given any nested quantified formula

F 1 in the form of Q 1 x 1 Q 2 x 2 ...Q k x k F where Q i ∈ {∀, ∃} and F is a quantifier-free formula, we can convert F 1 into a form of ∀x 1 ¬∀x 2 ...¬∀x n ¬F n
where F n is also quantifier-free. This leads to a sequence of formulas F i such that: F i ≡ ∀x i ¬F i+1 for i < n. Before QT is called, the initial values of M i are initialized to true for i < n and M n is initialized to F n , The final value of M 1 is false if and only if F 1 is unsatisfiable.

Theorem 1 (Partial Correctness). Assume C is satisfiable. The algorithm QT (i, C, x, M ) returns a pair of the form (C , M ) where C is an extrapolant of C, F i .

We do not provide a detailed proof, but we briefly discuss correctness and the invariant F i M i , that are mutually inductive. The main idea is that F i M i holds for the initialization step for QT, and F i ⇒ ¬(πx i .(C |M i )) also holds for each QT iteration; therefore, we strengthen M i in the end of QT and preserve the invariant at the same time. The if branches of QT return C, M i when M i cannot be strengthened any more. Moreover, QT also ensures that C, M i ∧ F i is unsat. Therefore, the last C, M i (corresponding to the strongest version of M i ) is also an extrapolant of C, F i . Termination: Algorithm QT does not terminate for arbitrary instantiations of projection and extrapolation. The projection and extrapolation operators we examine here are well-behaved (with respect to termination) in the following sense: (1) The extrapolation procedures do not introduce new atoms, so there will be only a finite number of new extrapolants one can make.

(2) The projection procedures are also finitary: they produce only a finite set of projections for the case of linear arithmetic.

Quantifier Elimination: algorithm QE

A quantifier-free version of F 1 is obtained by executing QT until saturation in Algorithm 2. The algorithm initializes a vector of formulas M and strengthens these formulas as much as possible in a loop. The intuition of QE is described as follows:

• Run QT (1, ¬f alse, x, M ), we obtain a formula C 1 where C 1 F 1 .

• Execute QT (1, ¬C 1 , x, M ), we get C 2 , a disjoint formula of C 1 , where C 2 F 1 . • Run QT (1, ¬(C 1 ∨ C 2 ), x, M ), we obtain a next formula C 3 where C 3 F 1 .
• When QT (1, ¬C, x, M ) returns false, we get C as a disjunction of disjoint formulas where

C ≡ C 1 ∨ C 2 ∨ ... ∨ C k and C F 1 .

Algorithm QT by example

We use a small example to illustrate the algorithm QT :

∀y ∃z (z ≥ 0 ∧ ((x ≥ 0 ∧ y ≥ 0) ∨ -y -z + 1 ≥ 0)) (1)
The formulas corresponding to (1) are:

F 1 = ∀y ¬F 2 , F 2 = ∀z ¬F 3 , F 3 = z ≥ 0 ∧ ((x ≥ 0 ∧ y ≥ 0) ∨ -y -z + 1 ≥ 0)
and quantifiers are:

x 1 = y, x 2 = z
Algorithm QT also maintains formulas M 1 , M 2 and M 3 that are initialized as follows:

M 1 = true, M 2 = true, M 3 = F 3 = z ≥ 0 ∧ ((x ≥ 0 ∧ y ≥ 0) ∨ -y -z + 1 ≥ 0)
It maintains the invariants:

F i M i , F V (M i ) ⊆ {x 1 , . . . , x i-1 } for i = 1, 2, 3 (2)
Finally, the algorithm propagates a context formula C i between levels. The context formula is updated during propagation. When C i is propagated from level i to i + 1 or i -1, it results in a formula (C i+1 or C i-1 ) that has non-empty intersection with C i and is implied by M i . This formula is an extrapolant of C i and M i as defined in Definition 1. The extrapolant C i on level i satisfies:

F V (C i ) ⊆ {x 1 , . . . , x i-1 } (3)
It contains only variables that are free above level i.

Let us run QT on the sample formula. In the initial state, vector M = true, true,

F 3 , C 1 = true, i = 1. 1. C 1 ∧ M 1 is true ∧ true. It is satisfiable, so let us choose an extrapolant C 2 for C 1 and M 1 . C 2 := true is an extrapolant because C 1 ∧ C 1 is satisfiable and ¬M 1 ∧ C 2 is unsatisfiable. Set i := 2.
2. C 2 ∧ M 2 is also satisfiable and similarly we set C 3 := true, i := 3.

3. C 3 ∧ M 3 , which is C 3 ∧ F 3 ,
:= ¬(x ≥ 0 ∧ y ≥ 0). 5. C 2 ∧ M 2 , which is ¬(x ≥ 0 ∧ y ≥ 0), remains satisfiable. Let us set C 3 := ¬(y ≥ 0), i := 3.
C 3 satisfies the conditions for being an extrapolant.

6. C 3 ∧M 3 (= C 3 ∧F 3 ) is still satisfiable. We return the extrapolant C 3 := z ≥ 0∧-y-z+1 ≥ 0 and set i := 2.

7. The formula y ≤ 1 implies ∃z.C 3 (they are actually equivalent), so we can update M 2 := M 2 ∧ ¬(y ≤ 1) and maintain the invariant

F 2 M 2 . Let us simplify M 2 , ¬(x ≥ 0 ∧ y ≥ 0) ∧ ¬(y ≤ 1), to x < 0 ∧ y > 1.
8. At this point M 2 implies x < 0. So the next extrapolant C 3 will also imply x < 0.

However, M 3 cannot be satisfiable with C 3 . We are done with level 3 and return false to level 2. In response, level 2 propagates the extrapolant C 2 := (x < 0 ∧ y > 1) up to level 1.

9. Similar to step 4, F 1 ≡ ∀y.¬F 2 implies that ∀y.¬C 2 ≡ ¬∃y.C 2 . So if we take πy.(C 2 |M 1 ) to be x < 0, then we can update M 1 := ¬(x < 0), which is x ≥ 0.

10. At this point let us check the levels below using x ≥ 0. Then M 2 ∧ x ≥ 0 is unsatisfiable, so there is no refinement under this assumption. Return false to level 1. Level 1 is done, and we conclude the formula is satisfiable and an output quantifier-free formula is x ≥ 0.

Extrapolation

A trivial extrapolant has been described in Section 2.1. We will here discuss two other versions of computing extrapolants A, B .

SMT-TEST

The approach used in [START_REF] Monniaux | A Quantifier Elimination Algorithm for Linear Real Arithmetic[END_REF][START_REF] Monniaux | Quantifier Elimination by Lazy Model Enumeration[END_REF] is to enumerate conjunctions of literals that satisfy A ∧ B. Suppose L := 1 , . . . , n are the literals in the satisfying assignment for A ∧ B. An extrapolant is the intersection of L and an unsatisfiable core of L ∧ ¬B. Our implementation of SMT-TEST extrapolation is using a single satisfiability check to extract a (not necessarily minimal) unsatisfiable core.

NNF strengthening

NNF strengthening is a process of deriving a stronger formula by replacing literals by f alse. We start with a formula C which is a transformation of B to NNF. For each literal in C in order, replace that literal by f alse and check the conditions for extrapolation so that ¬B ∧ C is unsat and A ∧ C is sat. The first check is redundant (which holds by construction) and the second is not redundant. An extrapolant C = A, B is computed according to Algorithm 3. NNF strengthening gives us stronger formulas which potentially help reduce the number of iterations in procedure QT. We currently check satisfiability in each round during NNF strengthening. Another approach is to evaluate C[l/ false] using a model for C.

Algorithm 3: NNF extrap- olant C ← N N F (B); foreach literal ∈ C do C ← C[l/ false]; if A ∧ C is sat then C ← C end end return C
4 Projection specialized to Linear Integer Arithmetic

We are here interested in LIA since LIA decision procedure is used as the core of DCMC and plays a central role in feasibility of the model-checking approach [START_REF] Michael | On Tool Support for Duration Calculus on the basis of Presburger Arithmetic[END_REF]. Duration Calculus (DC) is an extension of Interval Temporal Logic with the notion of accumulated durations allowing succinct formulation of real-time problems. Chop ( ) is the only modality appearing in DC; however, the modelchecking problem in DC is transformed to satisfiabilitychecking of a LIA formula in size exponential to the chop-depth [START_REF] Fränzle | Efficient Model Checking for Duration Calculus?[END_REF].

Cooper's algorithm for Presburger Arithmetic corresponds to quantifier elimination using virtual substitutions, and besides SMT-TEST for extrapolation we also consider strengthening formulas in negation normal form (NNF) by replacing literals by false. Virtual substitution methods work directly on formulas in NNF, so SMT-TEST is potentially not required.

Virtual substitutions

Virtual substitutions on LIA are performed by means of Cooper's algorithm. The algorithm removes quantifiers in the inside-out order using the following transformation:

∃x. φ ⇐⇒ φ[ /ax < t, ⊥/ax > t] ∨ δ i=1 ax<t φ[t + i/ax] ∧ δ | t + i
where δ is the least common multiple of all divisors d in divisibility constraints d | t and δ is the least common multiple of all coefficients a in comparison constraints ax t where a > 0 and ∈ {<, ≤, =, ≥, >}.

A quantified formula is transformed to a disjunction by a series of substitution steps. The disjunction can be represented symbolically (it corresponds to an existential quantifier over a finite domain) and may contain redundant disjuncts. The new divisibility constraint δ | t + i could be replaced by many divisibility constraints of small divisors a in each substitution. In this setting, smaller divisors of inner formulas lead to fewer number of case splits for the next quantifier alternation.

Contextual simplification

Our projection procedure πx.(C|M ) admits using a context M when processing C. The simplification depends on the strength of M and it helps to trim down unnecessary cases in a virtual substitution method later. This process is called contextual simplification which is easy to implement in an SMT solver. Algorithm 4 contains a procedure for contextual simplification of formula ϕ. It works with a logical context, asserts that ϕ is unequal to p and recurses over sub-formulas of ϕ to replace them by true or false.

The approach also applies to non-Boolean domains (although this generalization is not required for QT ): the auxiliary algorithm then takes terms of arbitrary types and checks if Algorithm 4: CS(ϕ) if ϕ is unsat then return false end let M be a model for ϕ; let p be a fresh propositional variable; assert p = ϕ; return CS(ϕ, M, p)

Algorithm 5: Auxiliary algorithm CS(ϕ, M, p) if (p = M(ϕ)) is unsat in current context then return M(ϕ) end foreach immediate subformula ψ i in ϕ[ψ 1 , . . . , ψ k ] do push; let p i be a fresh propositional variable; assert ϕ[ψ 1 , . . . , ψ i-1 , p i , ψ i+1 , . . . , ψ k ] = p; ψ i ← CS(ψ i , M, p i ); pop; end return ϕ[ψ 1 , . . . , ψ k ]
the terms are forced equal to the value provided in the model M. Instead of recursing on sub-formulas it can recurse on sub-terms of arbitrary types. When ϕ is represented as a DAG and has exponentially many sub-formulas, the result can in fact get much larger than the input. A practical implementation of CS should therefore impose limits on how many times a sub-formula is traversed.

Contextual simplification is used in STeP [START_REF] Bjørner | STeP: Deductive-Algorithmic Verification of Reactive and Real-Time Systems[END_REF] for simplifying formulas produced from verification condition generators and [START_REF] Dillig | Small formulas for large programs: On-line constraint simplification in scalable static analysis[END_REF] develops a routine that works on formula trees in NNF, where it is used to speed up abstract interpreters that propagate formulas. Algorithm 4 is used in Z3 with the observation that a model M for ϕ can be used to prune checks for value forcing and it applies to subterms of arbitrary type. Z3 also contains cheaper contextual simplification routines that rely on accumulating equalities during a depth-first traversal. While this cheaper algorithm can replace sub-terms by constants, it is not necessarily an advantage to use in context of SMT solving: a non-constant sub-term can be much more useful for learning general purpose lemmas during search.

Evaluation

We implemented different combinations for the instantiated algorithm in the preview version of Z3 4.0. Projection procedures have been used including (A) -full quantifier elimination and (B) -partial quantifier elimination. There are three variants of extrapolation: (0) -trivial extrapolation, (1) -NNF strengthening and (2) SMT-TEST. Furthermore, we also implemented (X) -contextual simplification and (Y) -no contextual simplification. These components in order constitute 12 different combinations which are named in short as aix where a ∈ {A, B}, i ∈ {0, 1, 2} and x ∈ {X, Y }. For example, A0Y denotes a combination of trivial extrapolation, full quantifier elimination and no contextual simplification.

We attempt to compare our algorithm (from now on called AQS algorithm) with Z3's quantifier elimination algorithm. Non-random benchmarks are collected from DCMC. Not only is the huge size (exponential to the chop-depth of DC formulas) of LIA formulas problematic, the nested nature between universal and existential quantifiers makes the problems even harder [START_REF] Michael | On Tool Support for Duration Calculus on the basis of Presburger Arithmetic[END_REF].

We divided benchmarks into two sets. Set 1 consists of 32 easy formulas having from 56 to 94 quantifiers with file sizes ranging from 15KB to 33KB in SMT-LIB format. Z3's quantifier elimination can process these formulas within a few seconds. They are chosen for the purpose of recognizing incompetent candidates in 12 combinations above. Set 2 have 64 hard instances with 69-768 quantifiers and take 50-500KB in SMT-LIB format. They are beyond the scope of Z3's quantifier elimination algorithm. We use them to test scalability of different combinations. Benchmarks and experimental results are available on the Z3 website 1 . The benchmark sets have some specific characteristics: coefficients are quite small and constraints are sparse (consisting of a few variables). These features help limit the number of disjunctions in virtual substitutions. Figure 1 summarizes the running time of different configurations of AQS as well as Z3's quantifier elimination algorithm. Each benchmark is given a 30-second timeout. The graph shows the accumulated running-time for solving all 32 problems. The winner is the configuration A0Y which means running AQS with projection implemented as full quantifier elimination and using trivial extrapolation. This configuration solves all benchmarks within 20 seconds. The configuration A2Y (using SMT-TEST instead of trivial extrapolation) is a close runner-up. It also solves all benchmarks, but requires 35 seconds. This benchmark set is simply too small to draw clear conclusions between these approaches. Partial quantifier elimination (configurations with prefix B) is bad on all configurations. The experiments also suggest that strong context simplification is pure overhead in all configurations. In an earlier prototype outside of Z3, however, strong context simplification was an advantage. We attribute this to how constraints get simplified and subsumed when being passed between AQS and Z3's quantifier elimination procedure. Z3's built-in quantifier elimination procedure has a slower ramp up time and is able to solve all problems within 175 seconds. Experiment 2 was performed on Set 2 for all configurations of AQS with timeout of 300 seconds. Experimental results are shown in Figure 2. We omit running time for configurations with partial quantifier elimination, they solve almost no problems, and we also omit running time for the default quantifier elimination routine that also does not solve any problem. The experimental results indicate that AQS algorithm performs well on formulas with many blocks of nested quantifiers. We have gained an order of magnitude speedup for Duration Calculus application. The A2Y (using full projection and SMT-TEST ) configuration scales well on our benchmarks although A0Y (using full projection and trivial extrapolation) comes intriguingly close.

Conclusions

We presented an anatomy of the algorithm proposed in [START_REF] Monniaux | Quantifier Elimination by Lazy Model Enumeration[END_REF] for checking satisfiability of formulas with alternating quantification. We proposed a set of generalizations, applied them to Presburger Arithmetic, and evaluated the generalizations to benchmarks from a model checker for Duration Calculus. So far the experience has been that the satisfiability algorithms, when instantiated with SMT-TEST (and to some extent trivial extrapolation) perform orders of magnitude better than general purpose quantifier elimination. We are currently investigating additional alternatives to the algorithms presented here. One alternative is to instantiate quantifiers incrementally using virtual substitutions. The idea is similar to how quantifiers are instantiated using E-matching in SMT solvers. SMT solvers create a propositional abstraction of formulas, including quantifiers. If there is a satisfying assignment to the abstracted formula that does not depend on the quantified sub-formulas, then the formula is satisfiable. Otherwise, quantified formulas are model-checked with respect to the current model and only instantiated (by axioms of the form "(∀xϕ[x]) ⇒ ϕ[t]") when the interpretation for free variables cannot be extended to an interpretation that satisfies the quantified formulas.

Figure 1 :

 1 Figure 1: Accumulated running time of AQS vs. Z3 on benchmark set 1

Figure 2 :

 2 Figure 2: Accumulated running time of AQS on benchmark set 2

  is satisfiable as well. In this case we will propagate back to level 2 a formula that intersects with C 3 and implies F 3 . So let us return C . Now F 2 ≡ ∀z.¬F 3 implies that ∀z.¬C 3 ≡ ¬∃z.C 3 . So we can strengthen M 2 with the negation of any formula that implies ∃z.C 3 . This is a projection as denoted by the notation πz.(C 3 |M 2 ) in Definition 2. It can take the current state of M 2 into account. In this case we set πz.(C 3 |M 2 ) to x ≥ 0 ∧ y ≥ 0 and update M 2

3 := z ≥ 0 ∧ x ≥ 0 ∧ y ≥ 0 and update the level i := 2.

4

Pascal Fontaine, Amit Goel (eds.); SMT 2012, pp. 118-128

http://research.microsoft.com/projects/z3/qt2012.zip

Acknowledgments

This work grew out of a course on SMT solving organized by Flemming Nielson and Hanne R. Nielson. We are grateful for valuable comments and permission of using DCMC benchmarks from Michael R. Hansen. Anh-Dung Phan is supported by the IDEA4CPS project granted by the Danish Research Foundation for Basic Research.