Di Ge
email: di.ge@irccyn.ec-nantes.fr

Eric Le Carpentier
email: eric.le-carpentier@irccyn.ec-nantes.fr

Dario Farina

Unsupervised Bayesian decomposition of multiunit EMG recordings using Tabu search

Unsupervised Bayesian Decomposition of Multiunit EMG Recordings Using Tabu Search

Di Ge, Eric Le Carpentier, and Dario Farina * , Senior Member, IEEE Abstract-Intramuscular electromyography (EMG) signals are usually decomposed with semiautomatic procedures that involve the interaction with an expert operator. In this paper, a Bayesian statistical model and a maximum a posteriori (MAP) estimator are used to solve the problem of multiunit EMG decomposition in a fully automatic way. The MAP estimation exploits both the likelihood of the reconstructed EMG signal and some physiological constraints, such as the discharge pattern regularity and the refractory period of muscle fibers, as prior information integrated in a Bayesian framework. A Tabu search is proposed to efficiently tackle the nondeterministic polynomial-time-hard problem of optimization w.r.t the motor unit discharge patterns. The method is fully automatic and was tested on simulated and experimental EMG signals. Compared with the semiautomatic decomposition performed by an expert operator, the proposed method resulted in an accuracy of 90.0% ± 3.8% when decomposing singlechannel intramuscular EMG signals recorded from the abductor digiti minimi muscle at contraction forces of 5% and 10% of the maximal force. The method can also be applied to the automatic identification and classification of spikes from other neural recordings.

Index Terms-Bayesian analysis, electromyography (EMG) signal decomposition, Tabu search.

I. INTRODUCTION

I NTRAMUSCULAR electromyography (EMG) signals are the sum of series of action potentials discharged by the motor units (MUs) within the detection volume of the recording electrodes. The interference intramuscular EMG signals can be decomposed into the constituent MU action potentials (MUAPs) to extract the discharge patterns of the motor neurons innervating the muscle. This procedure allows the analysis of MU control properties [START_REF] De | Control scheme governing concurrently active human motor units during voluntary contractions[END_REF] and thus provides a unique insight into the neural drive to muscles. EMG signal decomposition is usually performed by semiautomated methods, often based on serial detection and classification of action potentials.

Since the pioneer work completed by De Luca and coworkers both in the development of detection needles, processing techniques and in the use of EMG signal decomposition for understanding fundamental physiological issues [START_REF] De | Control scheme governing concurrently active human motor units during voluntary contractions[END_REF]- [START_REF] Lefever | A procedure for decomposing the myoelectric signal into its constituent action potentials: Part II. Execution and test for accuracy[END_REF], many techniques have been further developed for the EMG decomposition (for review, see [START_REF] Stashuk | Decomposition and quantitative analysis of clinical electromyographic signals[END_REF]). In classical approaches for intramuscular EMG decomposition, the MUAPs are identified from the interference signal with a threshold and classified as belonging to the active MUs. The challenge in the decomposition is mainly the separation of MUAPs that overlap in time. The global optimization of overlapped MUAPs is indeed an nondeterministic polynomial-time (NP) hard problem, i.e., a problem that cannot be solved by polynomial complexity algorithms. The existing methods either perform on the restrained search spaces [START_REF] Lewicki | Bayesian modeling and classification of neural signals[END_REF], [START_REF] Atiya | Recognition of multiunit neural signals[END_REF] to reduce the complexity or are based on recursive algorithms [START_REF] Mcgill | Optimal resolution of superimposed action potentials[END_REF], [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF] with specific trial strategies and residual threshold estimations. Atiya [START_REF] Atiya | Recognition of multiunit neural signals[END_REF] proposed a robust approach for decomposing overlaps of action potentials in neural recordings by comparing all possible combinations of up to two action potentials (restrained search space). Within the family of the spike sorting algorithms, the dynamic programming method [START_REF] Lewicki | Bayesian modeling and classification of neural signals[END_REF], which uses the fast exploration technique of a k-d tree, is also limited by the memory space necessary to generate such data structure, resulting in an equally restrained search space of up to two overlapping action potentials. This restraint is not justified in several applications, for example, when decomposing intramuscular EMG signals comprising many active MUs. Other approaches are not limited in the number of overlapping sources, however, they usually require the interaction with an operator. For example, the EMGLAB algorithm [START_REF] Mcgill | Optimal resolution of superimposed action potentials[END_REF], [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF] recursively matches the templates (reference MUAPs) and reevaluates the residual errors, similar to the greedy algorithm, but requires an expert intervention when the algorithm fails to reach the residual lower bound [START_REF] Mcgill | Optimal resolution of superimposed action potentials[END_REF] after a number of trials.

In this study, we propose a fully automatic EMG decomposition algorithm that exploits both the signal model likelihood and the regularity of the MU discharge patterns in a similar way as it is done during the manual decomposition. Unlike previous methods, which indirectly use the information on MU discharge regularity in an interactive procedure involving human intervention [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF], the proposed method integrates this information as prior law in a Bayesian framework. The NP-hard problem associated to the decomposition of overlapping MUAPs is solved with the Tabu metaheuristic without limiting the search space. Moreover, contrary to methods for the multichannel surface EMG decomposition [START_REF] Holobar | Correlation-based decomposition of surface EMG signals at low contraction forces[END_REF]- [START_REF] Gazzoni | A new method for the extraction and classification of single motor unit action potentials from surface EMG signals[END_REF], the proposed method requires one single recording channel. The method was tested on simulated and experimental signals.

II. MODEL AND METHOD

A. Forward Model

An interference EMG signal comprising the contributions of I MU discharge sources can be described as

z = I i=1 h i * s i + ǫ ( 1
)
where z is the recorded EMG signal of length N , modeled as a mixture of convolutions of the impulse trains s i ,i =1,...,I (the discharge patterns) and the linear filters h i ,i =1,...,I (MUAPs). The shorthand (•) i replaces the whole set {•,i = 1 ...,I} without ambiguity, thus (s i ) i = {s i ,i =1...,I}.

The following statistical assumptions are made on this mixture of sources. 1) Each impulse source s i =1 x i (the discharge pattern) is modeled as an arriving-time process with uniform amplitudes. The subscript of the indicator function x i is a vector containing the impulse coordinates for the ith source (arriving time of impulses). We also denote n i =dim(x i ), the number of discharges.

2) The MU discharge patterns (1 x i ) i are supposed mutually independent.

3) The MUAP shapes (h i ) i vary slowly during the EMG recording (i.e. variations in shape of the MUAPs occur during time intervals of seconds as order of magnitude). 4) The observation z is corrupted by an additive white Gaussian noise with unknown variance σ 2 ǫ . In relation to assumption 2), it should be noted that two types of correlations exist between MU discharge patterns. The first consists of the common modulation of the time-varying discharge rates between pairs of MUs (common drive) [START_REF] De | Control scheme governing concurrently active human motor units during voluntary contractions[END_REF]. The second is characterized by simultaneous or near-simultaneous discharges of two or more MUs (including but not limited to phase locked and entrained (1 x i ) i ) that occur more often than would be expected with independent processes (short-time synchronization [START_REF] Dietz | Correlation between the dischanges of two simultaneously recorded motor units and physiological tremor[END_REF], [START_REF] De Luca | Physiology and mathematics of myoelectric signals[END_REF]). Models exist (e.g., see [START_REF] Jiang | A simulation method for the firing sequences of motor units[END_REF]) to simulate the interdependent discharge patterns, by introducing both MU common drive and short-time synchronization among MUs. However, these two types of correlation between MU discharge patterns cannot be directly coded into the prior laws of (1 x i ) i . This ap r i o r iinformation is thus ignored in the present study to yield a simple model structure. It is also as a general rule in Bayesian inference problems that independency assumptions often imply lack of tractability to exploit the data dependency properties rather than proved scientific observations.

The data-generating process, given the source parameters

(x i , h i ) i ,σ 2
ǫ , obeys the Gaussian law

P z|(x i , h i ) i ,σ 2 ǫ = 1 2πσ 2 ǫ N exp - ||z -i 1 x i * h i || 2 2σ 2 ǫ . (2)
The following additional assumptions are made on the discharge pattern of each MU.

1) The interspike interval (ISI) T ij = x i,j+1x i,j between two consecutive impulses for a given source s i is larger than a threshold value T R , which indicates a physiological constraint on the ISI (e.g., T R m a yb es e ts i m i l a rt ot h e refractory period of muscle fibers). 2)

(T ij -T R ) follows a Gaussian-shaped distribution (T ij - T R ) ∼N (m i ,σ 2 i ).
3) The ISI variability is smaller than a threshold: σ i /m i < Th σ i ; this constraint expresses the regularity of the discharge patterns and the threshold Th σ i controls the variability of the ISI. According to the previous assumptions, the ISI follows a truncated Gaussian distribution, in which

P (T ij )=0,T ij <T R P (T ij -T R ) ∝ g(m i ,σ 2 i ),T ij ≥ T R . ( 3
)
This assumption is in accordance with the histograms of experimentally observed discharge patterns [START_REF] Moritz | Discharge rate variability influences the variation in force fluctuations across the working range of a hand muscle[END_REF], i.e. approximately Gaussian-shaped with a normalized variance σ i /m i < Th σ i .On a discrete grid, the ISI probability in (3) is well-defined up to a normalization factor. The source-generating process is thus the following:

P x i | m i ,σ 2 i = 1 4 erfc x i,1 -T R -m i √ 2σ i erfc N -x i,n i -T R -m i √ 2σ i 1 2πσ 2 i (n i -1)/2 exp   - n i -1 j =1 (x i,j+1 -x i,j -m i -T R ) 2 2σ 2 i   (4)
for a given configuration of a source s i =1 x i . The two terms 1/2erfc(•) in (4) evaluate the first impulse ISI with probability P (T i,1 >x i,1 ) and the (n i +1) (th) impulse ISI with probability P (T i,n i +1 >N-x i,n i ) respectively. For simplicity, the endeffects by integration are neglected in the following, so that we obtain

P x i | m i ,σ 2 i ≈ 1 2πσ 2 i (n i -1)/2 × exp   - n i -1 j =1 (x i,j+1 -x i,j -m i -T R ) 2 2σ 2 i   . (5)
With the independence assumption of MU discharge patterns (x i ) i , we also obtain

P (x i ) i | (m i ,σ 2 i ) i = I i=1 P x i |m i ,σ 2 i . (6)
In the following, P ((x i ) i ) will denote the prior law for the discharge patterns to simplify the notation.

B. Posterior Probability Law

In a Bayesian framework, we establish the posterior distribution for the unknown variables (x i ,m i ,σ 2 i , h i ) i ,σ 2 ǫ |z :

P (x i ,m i ,σ 2 i , h i ) i ,σ 2 ǫ |z ∝ P z|(x i , h i ) i ,σ 2 ǫ I i=1 P x i |m i ,σ 2 i P (m i ) P σ 2 i P (h i ) P σ 2 ǫ . (7)
The decomposition task is accomplished by jointly maximizing (7) w.r.t. the discharge patterns (x i ) i , their statistics (m i ,σ 2 i ) i , the MUAP shapes (h i ) i and the background noise variance σ 2 ǫ . Equations ( 2) and (4) provide the first two terms of the righthand side of [START_REF] Atiya | Recognition of multiunit neural signals[END_REF]. The remaining terms can be expressed as conjugate priors with noninformative hyperparameters

(α i ,β i ) i , (α s ,β s ,µ 0 ,σ 2 0 ,σ 2 h ) P (m i ) ∼N(µ 0 ,σ 2 0 ),P (σ 2 i ) ∼IG(α i ,β i ) (8)
P (h i ) ∼N(h (0)
i ,σ 2 h ),P (σ 2 ǫ ) ∼IG(α s ,β s ) (9) where IG stands for the inverse Gamma distribution.

For example, the inverse Gamma IG(1, 1) yields an almost uniform distribution that covers a wide range of parameter values. The choice of µ 0 ,σ 2 0 expresses the prior information on the discharge rate of each MU and the certainty of this information. Typically, the Gaussian prior of m i can be interpreted as conveying no specific prior information under the condition that [25 ms, 200 ms] ⊂ [µ 0 -3σ 0 ,µ 0 +3σ 0 ] (i.e. most of the physiological discharge rates ranging from 5 to 40 Hz are considered, corresponding to a mean value of 200 and 25 ms, respectively). The same remark is applied to σ 2 h that controls the uncertainty of the initialized shapes (h (0) i ) i , which are the approximated MUAP obtained from the preprocessing phase (see Section III-A).

III. MAXIMIZATION ALGORITHM

This section describes the main features of the two-phase EMG signal decomposition algorithm based on the maximum a posteriori estimator (MAP) applied to the forward model described earlier; the joint posterior distribution of ( 7) is maximized w.r.t. the unknown parameters

Θ = arg max P (Θ|z) with Θ= (x i , h i ,m i ,σ 2 i ) i ,σ 2
ǫ . The Section III-A presents the overall structure of the method, divided into a preprocessing phase and a complete EMG separation phase that iteratively maximizes [START_REF] Atiya | Recognition of multiunit neural signals[END_REF]. For fixed

(x i ) i , maximization w.r.t. (h i ,m i ,σ 2 i ) i ,σ 2 ǫ
leads to a closed-form solution. On the other hand, the problem of estimating the discharge pattern (x i ) i cannot be practically solved by an exhaustive exploration due to computational time [START_REF] Lewicki | Bayesian modeling and classification of neural signals[END_REF], [START_REF] Lewicki | A review of methods for spike sorting: The detection and classification of neural action potentials[END_REF]. The decomposition problem consists of determining the MU label set and the discharge instants that maximize the posterior distribution expressed by [START_REF] Atiya | Recognition of multiunit neural signals[END_REF]. The complete search space of the discharge pattern (x i ) i , containing all possible alignments of all possible combinations of MU label sets, is exponential [START_REF] Gut | High-precision EMG signal decomposition using communication techniques[END_REF]. In this paper, we propose to use the Tabu search algorithm [START_REF] Glover | Tabu search-Part I[END_REF] to deal with the maximization w.r.t. the discharge pattern (x i ) i .

We clarify the main application issues such as the cost function to minimize (equivalent to the maximization of the marginal posterior law for (x i ) i ), the neighborhood definition (the set of moves allowed in one iteration), and the Tabu list structure.

A. Overall Structure

The decomposition method consists of two phases. First, the EMG signal is segmented and a few representatives of the MUAPs (h (0) i ) i are extracted for initialization of MUAP shapes [see [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF]]. It is implemented as a classic preprocessing for MUAP detection in intramuscular recordings and comprises bandpass filtering, thresholding for detection [START_REF] Lewicki | A review of methods for spike sorting: The detection and classification of neural action potentials[END_REF], [START_REF] Gut | High-precision EMG signal decomposition using communication techniques[END_REF], and isolated MUAP shape classification [START_REF] Mcgill | Automatic decomposition of the clinical electromyogram[END_REF]. In this study, we applied a level threshold proportional to the background noise standard deviation estimate σ ǫ .Theactive segments ({Seg k }) are then detected where either isolated or overlapped MUAPs occur. False positives are preferred over false negatives since (x i ) i = φ in one active segment is an allowed solution of the MAP estimation (the second phase).

The algorithms used in the first phase, inspired from the methods in [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF], [START_REF] Mcgill | Automatic decomposition of the clinical electromyogram[END_REF], are as follows.

1) Bandpass filtering of the raw EMG signal to enhance the MUAP's fast rising edges while suppressing the low-frequency background activities. The first-order differentiator in [START_REF] Mcgill | Automatic decomposition of the clinical electromyogram[END_REF] is adopted for its simplicity, as in [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF]:

z f n = z n +1 -z n -1
, where z and z f represent the raw EMG data and its filtered output, respectively. In the following, z ← z f while we keep the notation z for filtered EMG.

2) Segmenting the prefiltered EMG signal into temporal intervals that contain either overlapped or separated MUAP activities. A level threshold Th Seg is calculated based on a quick estimation of σ ǫ by solving the implicit equation

Th Seg =4 T 1 n =1 z 2 n I(Th Seg ,n) T 1 n =1 I(Th Seg ,n) (10)
in which the indicator function I(p, n)=1if |z n | <pand 0 otherwise. We fix T 1 as the first 1-s interval of z, i.e. 10 000 samples if the sampling frequency is 10 kHz and the coefficient 4 in (10) corresponds to detect only MUAP activities whose absolute amplitudes exceed 4σ ǫ . The solution of ( 10) is of complexity O(T 1 ) (code given in EMGLAB [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF]). We then detect all the local extrema of z whose absolute values exceed Th Seg (empirically believed to be MUAP peaks). A new segment is allocated when the distance between two adjacent extrema exceeds 2P +1, P being the length of (h i ) i . An example of two segments is provided in Fig. 2 for which P is set to 5 ms.

3) Extracting recurrent isolated MUAP forms from the segments to obtain (h (0) i ) i . Canonically registered discrete Fourier transform (CRDFT) coefficients (see [START_REF] Mcgill | Automatic decomposition of the clinical electromyogram[END_REF] and [START_REF] Mcgill | High-resolution alignment of sampled waveforms[END_REF], and references therein for implementation details) are calculated on segments with lengths smaller than 3P to exclude those potentially containing overlapped waveforms and a quadratic metric is adopted on these CRDFT coefficients to make a one-pass classification operation. The progressively averaged CRDFT coefficients are considered as the DFT of the representative templates of (h (0) i ) i , used as initialization in the second phase. The segment length depends on the signal properties and the parameters used for segmentation, such as the sparsity of z (basically depending on the number of active MUs) and the chosen MUAP length P . Since the optimization of (7) w.r.t. (1 x i ) i is only performed for intervals of the active segments, the segmentation step reduces the computational burden of the second phase. The first phase of the decomposition, as described earlier, is completely automated and requires a negligible computational load compared to the second phase, that solves the problem of MUAP overlaps.

The novel part of the decomposition method consists of the second phase. After segmentation, the maximizationdecomposition phase is applied serially to each active segment {Seg k } to obtain: 1) the complete discharge patterns (x i ) i for each MUAP; 2) the MUAP shape estimates (h i ) i ; and 3) the background noise variance σ 2 ǫ , with an MAP criterion. For each segment, the joint posterior distribution [see [START_REF] Atiya | Recognition of multiunit neural signals[END_REF]] is maximized over the complete search space of the discharge patterns (x i ) i in an automatic algorithm. The kernel of this problem (MUAP overlapping) is dealt with by a metaheuristic algorithm, called Tabu search, first introduced by Glover [START_REF] Glover | Tabu search-Part I[END_REF] in the area of operational research. The generic Tabu search algorithm and its application are described in Section III-C. The optimization strategy is designed to overcome limitations such as early convergence into local optima in NP-hard problems and is thus adapted to solve the combinatorial problem of MUAP overlaps.

The MAP optimization are summarized as follows. 1) Initialize (h (0) i ) i (MUAP shapes) using the information from phase 1 (segmentation). 2) Initialize (m i ,σ 2 i ) i by their conjugate prior laws, and σ ǫ by σ ǫ from (10).

3) Initialize x i = φ (1 x i = 0) for all i. 4) Iterate: a) For each Seg k , maximize:

P ((x (k ) i ) i |(x (-k ) i , h i ,m i ,σ 2 i ) i ,σ 2 ǫ , z) w.r.t. (x (k ) i ) i using Tabu search b) Maximize (7) w.r.t. continuous parameters (h i ,m i ,σ 2 i ) i ,σ 2 ǫ (see Section III-B); until convergence. In the previous steps, x (k ) i = x i ∩ Seg k ,i =1,.
..,I,denote the discharge coordinates vector of each MU within the given Seg k and the superscript (-k) in the following indicates variables belonging to segments ∪ j =k Seg j .

The calculations in the step 4(b) are straightforward maximization problems, for which we provide detailed expressions in Section III-B. Step 4(a), however, is a typical NP-hard combinatorial problem whose complete space amounts to 2 I dim(Seg k ) configurations for the kth active segment, where I denotes the number of classes (MUs) and dim(Seg k ) the length of the segment. Its maximization algorithm is explained in Section III-C using the Tabu search metaheuristic. Step 4(a) is solved by an iterative Tabu search that runs until convergence is detected either by exhausting all valid neighbors or by reaching the max-imum allowed iteration counts M as stated in the algorithm in Section III-C. The external loop is run for both the steps 4(a) and (b) until the global convergence of Θ. Jointly, the steps 4(a) and (b) bear a certain resemblance to a Gibbs simulation scheme [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF]. However, while a Gibbs simulation scheme iteratively generates samples asymptotically distributed to [START_REF] Atiya | Recognition of multiunit neural signals[END_REF], the proposed algorithm iteratively maximizes the posterior probability in [START_REF] Atiya | Recognition of multiunit neural signals[END_REF], and finally, builds an estimate of Θ in an MAP sense.

B. Maximization on Continuous Parameters

We describe first the maximization procedures in step 4(b). 1) Maximization w.r.t. m i : The maximization problem can be written as: for all i, m i = max

m i P (Θ|z) = max m i P (x i |m i ,σ 2 i )P (m i ) = max m i   exp    - 1 2σ 2 i n i -1 j =1 (S ij -m i ) 2    P (m i )   where {S ij = x i,j+1 -x i,j -T R } is the set of Gaussian sam- ples. P (x i |m i ,σ 2 i
) is given in [START_REF] Stashuk | Decomposition and quantitative analysis of clinical electromyographic signals[END_REF]. Since its prior law P (m i ) is also a Gaussian and the product of two Gaussian is still Gaussian, the maximum is found at

m i = (µ 0 /σ 2 0 +( n i -1 j =1 S ij )/σ 2 i ) (1/σ 2 0 +(n i -1)/σ 2 i ) . ( 11
)
2) Maximization w.r.t. σ 2 i : The maximization problem can be written as: for all i,

σ 2 i = max σ 2 i P (Θ|z) = max σ 2 i P (x i |m i ,σ 2 i )P (σ 2 i ) = max σ 2 i σ -(n i -1) i exp    - 1 2σ 2 i n i -1 j =1 (S ij -m i ) 2    P (σ 2 i ).
Since the inverse Gamma distribution is the conjugate prior for Gaussian likelihood with known mean, their product is still inverse Gamma with the maximum found at

σ 2 i = β i +( n i -1 j =1 (S ij -m i ) 2 )/2 α i +1+(n i -1)/2 . ( 12
)
To take into account the regularity constraint that σ i /m i < Th σ i , we impose that σ 2 i ← min{σ 2 i , (Th σ i m i ) 2 }. This constraint can be tuned depending on the expected regularity of the discharge pattern.

3) Maximization w.r.t. (h i ) i : We propose to jointly maximize (h i ) i by creating a concatenated column vector of size IP × 1: h =[h 1 ,...,h I ] and the maximization problem can be written as

h = max h P (Θ|z) = max h P (z|(x i ) i , h,σ 2 ǫ )P (h) = max h exp - 1 2σ 2 ǫ z -Mh 2 P (h) with M =[M 1 ,...,M I ],
where each submatrix M i is the Toeplitz matrix of convolution such that

M i h i =1 x i * h i .
It follows that z -Mh = zi 1 x i * h i . Since P (h) ∼ N (h (0) ,σ 2 h I), the maximum of the marginal posterior distribution is found at

h = σ -2 h I + σ -2 ǫ M t M -1 σ -2 h h (0) + σ -2 ǫ M t z . ( 13
)
We note that the matrix σ -2 h I + σ -2 ǫ M t M is PI × PI positive definite; thus, its inversion operation can be replaced by a Cholesky factorization followed by matrix left division to reduce the calculation load. A minimum degree preordering can also be used since M t M is sparse for any MU source. The motivation behind the joint maximization of (h i ) i is to ensure one-step convergence to the joint Gaussian mode, instead of a coordinate-by-coordinate maximization.

4) Maximization w.r.t. σ 2 ǫ : The maximization problem can be written as follows:

σ 2 ǫ = max σ 2 ǫ P (Θ|z) = max σ 2 ǫ P (z|(x i , h i ) i ,σ 2 ǫ )P (σ 2 ǫ ) = max σ 2 ǫ σ -N ǫ exp - 1 2σ 2 ǫ C P (σ 2 ǫ )
where C = z -Mh 2 denotes the energy of the reconstruction error. And the maximum is found at

σ 2 ǫ = β s + C/2 α s +1+(N/2) . (14)

C. Maximization on Combinatorial Sets

A Tabu search approach [START_REF] Glover | Tabu search-Part I[END_REF] is applied to solve step 4(a) as the following optimization problem:

Minimize: c(y) ∈R: y ∈ Y. ( 15
)
The objective function c(y) should be a real-valued cost function, while Y is a discrete set containing all combinations of y. Associated with each y ∈ Y is the set Neighbor(y) composed of all possible moves (trial solutions) from a current configuration y. A subset T ⊂ Y ,theTabu list, comprising prohibited moves is managed using a first-in-first-out (FIFO) structure. All allowed trials from a given state y belongs to the set Neighbor(y)\T .The Tabu search tackles the optimization problem in (15) as follows:

First, the nonimproving moves are made possible through step 4 since we do not necessarily have c(y m ) ≤ c(y m -1 ):a best move rather than an improving move is chosen. Second, visited local optima are forbidden for future moves by the FIFO mechanism of T that remembers to a certain degree (depending on the FIFO length) local optima to avoid. We begin by reformulating the step 4(a) from an optimization point of view. Let

n (k ) i =dim(x (k ) i ) denote the discharge counts of MU i in Seg k . It follows that arg max (x i ) ( k ) i P ((x i ) (k ) i |(x (-k ) i , h i ,m i ,σ 2 i ) i ,σ 2 ǫ , z) = arg max (x i ) ( k ) i exp      -z (k ) -I i=1 h i * 1 (k ) x i 2 (2σ 2 ǫ )      I i=1 σ -n ( k ) i i exp   - n ( k ) i j =0 x (k ) i,j+1 -x (k ) i,j -m i -T R 2 2σ 2 i    = arg min (x i ) ( k ) i 1 2σ 2 ǫ z (k ) - I i=1 h i * 1 (k ) x i 2 + log σ i I i n (k ) i + I i    n ( k ) i j =0 x (k ) i,j+1 -x (k ) i,j -m i -T R 2 2σ 2 i    (16)
where for MU i, x

i,0 and x

(k ) i,n ( k ) i

+1

denote the last discharge in the (k -1)th segment and the first discharge in the (k + 1)th segment, respectively. The first term in the criterion of ( 16) can be identified as the likelihood of the EMG signal z given (h i , 1 x i ) i and the last two terms as the prior information (regularity) of the discharge patterns. From the Bayesian point of view, they are naturally combined through the marginal posterior distribution on (x i ) i in the equation.

In the following, the superscript (k) is omitted without ambiguity since only (x i ) (k ) i is considered for each segment. The following elements are defined to apply the Tabu search:

1) y =(x i ) i and Y =({0, 1} dim(Seg k ) ) I define a state and the search space respectively; 2) the cost function c(y) is expressed in (16); 3) a cyclic fixed-length T (Tabu list) is updated to keep track of previous solutions; each state is characterized by MUAP classes and discharge instants; 4) for a given state, Neighbor((x i ) i ) is defined as:

(x * i ) i i 1 x * i -1 x i 1 ≤ 2, i |n * i -n i | < 2
in which |•| 1 defines the L1-norm of the vector (or the Hamming distance in the case of two binary sequences). It can be verified that such a neighborhood allows either adding a spike at any unoccupied position of the same MUAP, removing or changing the position of an existent spike. Furthermore, the neighborhood is both selfincluded ((x i ) i ∈ Neighbor((x i ) i )) and reflective ((

x i ) * i ∈ Neighbor((x i ) i ) ⇀ ↽ (x i ) i ∈ Neighbor((x i ) * i )).
It is also easily verified that the cardinality of Neighbor((x i ) i ) amounts to (I + i n i ) • dim(Seg)i n 2 i , or approximately (I + i n i ) • dim(Seg) for dim(Seg) ≫ n i . Because of the constraint on ISI (for example, indicating the absolute refractory period of muscle fibers), the number of discharges n i for each MU in a given active segment is limited. The complexity of each Tabu iteration is therefore linear w.r.t. the number of sources I and the data length dim(Seg), in contrast to the NP-hardness of the space.

D. Maximization on Two Segments

The segment-by-segment maximization method, as described earlier, often converges to local rather than global maximal configurations of (x i ) i . It was experimentally observed that this phenomenon occurs when the following two conditions concur (see Section V-B for illustration): 1) two or more different MUs have similar action potential shapes; and 2) at least two consecutive segments are separated by a time interval smaller than T R . Under these conditions, an error in the estimation of (x i ) i in one segment is susceptible to yield a systematic estimation error on the adjacent segment. A solution to the early convergence problem is to maximize jointly on two consecutive segments using a similar algorithm structure of Tabu search as in Section III-C, if the two segments are separated by less than T R .

For this purpose, we define the discharge patterns in two connected segments as a concatenated vector:

x (k,k+1) i = [x (k ) i ; x (k +1) i ],
where k and k +1 denote the segments index. z (k,k+1) covers the interval of both segments and the in-between area of noises and n

(k,k+1) i = n (k ) i + n (k +1) i
. Thus, the optimization function reads arg max

(x i ) ( k,k+1) i P ((x i ) (k,k+1) i |(x -(k,k+1) i , h i ,m i ,σ 2 i ) i ,σ 2 ǫ , z) = arg min (x i ) ( k,k+1) i 1 2σ 2 ǫ z (k,k+1) - I i=1 h i * 1 (k,k+1) x i 2 + I i n ( k,k+1) i j =0 x (k,k+1) i,j+1 -x (k,k+1) i,j -m i -T R 2 2σ 2 i + log σ i I i n (k,k+1) i . ( 17
)
Similarly to [START_REF] Moritz | Discharge rate variability influences the variation in force fluctuations across the working range of a hand muscle[END_REF] , allowing either to add, remove or relocate an existent MUAP across two segments. The neighborhood space is enlarged to (I + i n

(k,k+1) i )(dim(Seg k )+dim(Seg k +1 )) from (I + i n (k ) i )dim(Seg k ).
Note that for the two-segment approach the iterative maximization process is also carried out for each value of k.

With these modifications, a false label-switched configuration in one segment is no longer stable, since relocating an MUAP to the adjacent segment would further maximize [START_REF] Atiya | Recognition of multiunit neural signals[END_REF]. Comparisons are made in Section V-B on experimental EMG recordings to confirm that the two-segment maximization technique outperforms the one-segment version.

IV. VALIDATION METHODS

Validation tests were conducted using a PC equipped with a 3.2-GHz Intel Pentium 4 CPU processor and a 3 G RAM.

A. Simulation Model

An EMG signal generation tool designed by Farina et al. [START_REF] Farina | A model for the generation of synthetic intramuscular EMG signals to test decomposition algorithms[END_REF] was used to simulate intramuscular EMG signals. In order to test the performance of the proposed decomposition algorithm, the following model parameter values were fixed (see [START_REF] Farina | A model for the generation of synthetic intramuscular EMG signals to test decomposition algorithms[END_REF] for details on the model).

1) Four to five MUs were activated in a single channel with ∆E max (the maximum normalized energy difference between MUAPs of different classes) set to 0.7. 2) Discharge rates of the active MUs were set in the range of 10 -15 Hz. 3) Maximum degree of overlapping between action potentials of different MUs was set as α % ,max =1(complete overlap). 4) Percentage of MUAPs in the interference signal which overlapped in time was set to N % s = 70%. 5) SNR ratio was SNR dB =10. Five signals with the earlier properties were generated to test the proposed method. The sampling frequency of the simulated signals was 10 kHz. For both simulation and experimental tests in the following, the value of T R was fixed to 10 ms, i.e. 100 samples for a 10 kHz sampling frequency. Simulated signals had duration of 10 s.

B. Experimental Recordings

Five healthy men (age, mean ± SD =25.3 ± 4.5 years) participated in the study. The study was conducted in accordance with the Declaration of Helsinki, approved by the local Ethics Committee, and written informed consent was obtained from all subjects prior to participation.

The EMG measurements were obtained from the abductor digiti minimi muscle. The fifth finger was fixed in a custommade brace to record the force exerted during an isometric contraction of the muscle (Politecnico di Torino, Torino, Italy). The subjects performed three maximal voluntary contractions (MVCs) with the abductor digiti minimi, with each trial separated by a resting period of 2 min. The greatest force was used as the reference MVC for the other contractions. A pair of wire electrodes made of Teflon coated stainless steel (A-M Systems, Carlsborg, WA; diameter 50 µm) was inserted into the muscle with a 25-G needle. The insulated wires were cut to expose only the cross section at the tip. The needle was inserted to a depth of a few millimeters below the muscle fascia and removed to leave the wire electrodes inside the muscle. The intramuscular EMG signals were amplified and provided one bipolar recording (Counterpoint EMG, DANTEC Medical, Skovlunde, Denmark) that was bandpass filtered (500 Hz-5 kHz) and sampled at 10 kHz. A reference electrode was placed around the wrist. Each subject performed two 60-s contractions at forces of 5% and 10% MVC in random order and with 15 min of rest in between.

As reference, the MUAPs were identified from the intramuscular recordings with a validated decomposition algorithm (EMGLAB) [START_REF] Mcgill | Optimal resolution of superimposed action potentials[END_REF], [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF]. This interactive algorithm includes a user interface for manually editing and verifying the results [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF]. The EMGLAB displays a segment of the EMG signal, the templates of the action potentials of the identified MUs, the discharge patterns, and a close-up of the signal for resolving missed discharges and overlaps. Full, regular patterns provided confidence that the decomposition was correct, whereas gaps, extra discharges or uneven intervals provided an indication of possible decomposition errors. To assist in identifying missed discharges the program displays bars in the signal panel that indicate the expected discharge times of each MU. The signal portion can then be viewed in the close-up panel which displays the signal at an expanded scale and allows matching MU templates to be selected. The close-up panel also displays overlaps at an expanded scale which allows visual verification of the result and different sets of templates to be selected and adjusted to find the correct fit. Commands are also available for undoing identifications and deleting or merging templates. The acquired EMG signals were manually inspected with this user interface by an expert operator. The decomposition with EMGLAB has been previously cross-validated by decomposing on pairs of simultaneous EMG signals containing the same MU activities [START_REF] Mcgill | Validation of a computer-aided EMG decomposition method[END_REF]. The discharge patterns estimated in this way were used as reference. To evaluate the performance of the proposed methods, we followed the accuracy criterion as introduced in [START_REF] De Luca | Decomposition of surface EMG signals[END_REF] for the ith MU train

A(i)= N Dis (i) -N FP (i) -N FN (i) N Dis (i) × 100%
where N Dis (i) is the number of "true discharges" of the MU and N FN (i) and N FP (i) are, respectively, the number of false negatives and the number of false positives produced by the unsupervised decomposition algorithms for the ith MU when compared to the reference from manual decomposition (the true discharges were assumed to be those extracted with manual decomposition with EMGLAB). An MU discharge identified by the proposed method was considered a correct estimate if it was detected within a window of 1 ms centered at the time of occurrence of the true discharge. It was also shown in [START_REF] De Luca | Decomposition of surface EMG signals[END_REF] that by adopting

A = 1 I I i=1 A(i)
as the overall decomposition accuracy, any MU train is of the same significance regardless of the MUAP shape, duration and number of discharges. σ i /m i < 0.3 corresponded to the limit of regularity apriori imposed (this limit can be changed by the user if necessary). Thus, among the MUs identified by the proposed method, those with σ i /m i > 0.3 were automatically excluded from the analysis. The overall decomposition accuracy A was thus averaged only on those A(i) values corresponding to MUs under the condition that σ i /m i < 0.3. Results for experimental signals were obtained from the decomposition of the first 10 s of the experimental recordings.

V. R ESULTS

A. Simulations

Fig. 1(a) shows a portion of a simulated signal. The MUAP shapes are illustrated in Fig. 1(c). The labels assigned to each MUAP are also displayed. In this example, the decomposition process took less than 5 min per iteration as described in steps 4(a) and (b) of the proposed method on an EMG signal of duration 10-s sampled at 10 kHz. Two iterations were run in this test and in all following results. Fig. 1(d) shows a segment with four overlapping MUAPs, in which the positive lobe of MUAP #3 is compensated by the negative lobes of MUAP #1 and MUAP #2 (destructive overlap).

Table I shows the decomposition results on five simulated EMG signals. The MU discharge patterns from the decomposition were validated if σ i /m i < 0.3, and the overall efficiency was calculated on validated MU sources. The one-segment maximization algorithm was sufficient for accurate decomposition of all simulated EMG signals.

B. Experimental Signals

In the experimental EMG signals recorded for 10% MVC force, the total length (duration) of the active segments corresponded to approximately 50% of the total recording duration when 8 MUs were active. Among the active segments, those containing single MUAP (length inferior to 3P ) accounted for approximately 50% of the total segments.

Contrary to the simulated signals, the one-segment maximization method was not optimal in the case of decomposition of the experimental EMG signals. Fig. 2 compares the two versions of the automatic decomposition method (one-and two-segment maximization) in the case of a 10% MVC EMG signal comprising seven active MUs. The MU labels in the upper panels of both Fig. 2(a) and (b) was obtained by trained experts using EMGLAB, while the middle panels show the reconstructed EMG signal based on the automatic decomposition method. In this example, the label switched pair {#5, #7} in Fig. 2(a) appears across the two segments that are sufficiently close to each other (within 5 ms); their MUAP shape estimates given in Fig. 2 would yield an unusually small ISI and penalize the posterior distribution through the prior laws on the discharge pattern of x 7 . The joint optimal solution on two segments is more "global" in a sense that the corresponding residual level is lowered, while the ISI regularity is almost unaltered, yielding a higher posterior probability [see Fig. 2(b)].

TABLE II COMPARISON OF TWO VERSIONS OF OPTIMIZATION METHOD

To better illustrate the efficiency of the two-segment joint maximization method in comparison to the one-segment maximization method, decomposition statistics for each MU source are listed in Table II for one representative signal. Three MU discharge patterns are validated out of five for a two-segment joint maximization algorithm while only one of them is validated for the one-segment version. Some A(i) are not cal-

TABLE III

ACCURACY OF THE TWO-SEGMENT OPTIMIZATION ON EXPERIMENTAL SIGNALS culated as N FP (i)+N FN (i) exceeds the grand total of true discharges for some sources in the first case. From these examples, it is evident that the decomposition quality of the twosegment maximization is superior to that of the one-segment maximization, while the computational time is less than doubled. The results also show that the two-segment optimization method corresponds to a better tradeoff between computation time and accuracy as a function of segment length for adjacent segments separated by less than 10 ms apart from each other as in Fig. 2. The accuracy would still improve by jointly optimizing the marginal posterior law w.r.t. all the segments (the entire EMG signal); however, the Tabu search computational load would be unacceptable in this case because the neighborhood complexity is dependent of the segment length (see discussions in Section III-C and III-D); at the other extreme, a segmentby-segment strategy requires less computational time but is less reliable (see Fig. 2). The following results on experimental signals were all obtained with the application of the two-segment optimization.

Fig. 3 reports an example of decomposition of a segment which contains seven MUs. This example shows that the method proposed is not limited by the number of overlapping MUAPs. as identified by an expert operator, the number of validated MU discharge patterns and the overall accuracy rate A.N o less than half of the MU sources were validated using the criterion σ/m < 0.3 and A was greater than 85% in all experimental results. The average accuracy was 91.5% ± 4.4% and 88.4% ± 2.8% for the signals recorded during 5% and 10% MVC contractions, respectively. The computational time required for the experimental signal decomposition (10-s duration) was in the range of 15-30 min, longer than for simulated signals mainly because of the higher discharge rates in experimental recordings.

The instantaneous discharge rates for the two validated MUs from the signal EMG #5 at 10% MVC (see Table III) are reported in Fig. 4(a) and (b), as estimated by the proposed automatic decomposition method. The estimated MUAP shapes are reported in Fig. 4(c) and (d). False-positive and false-negative statistics are also given in Fig. 4.

VI. DISCUSSION

We have proposed a method for full decomposition of singlechannel intramuscular EMG signals that does not require any manual interaction. The attained accuracy approximates 90%, compared to the results obtained by an expert operator with semiautomatic decomposition methods, for the case of contraction levels up to 10% MVC and at most eight concurrently active MUs. In practical applications, this accuracy would be sufficient for the analysis of MU properties, such as the discharge rate and ISI variability. However, it may be insufficient to investigate physiological mechanisms for which a very precise decomposition is required, such as MU synchronization [START_REF] Nordstrom | Estimating the strength of common input to human motoneurons from the crosscorrelogram[END_REF]. A manual decomposition may still be required in this case.

The proposed MAP estimator provides an efficient and fully unsupervised approach to cope with the decomposition of multiunit EMG signals. The Bayesian statistical model on the EMG data generation process allowed the inclusion of available prior information, such as the Gaussian-like distribution of ISIs, the refractory period of muscle fibers and the regularity in the spike discharge trains. The discharge pattern regularity has long been viewed as either a hint for the expert interaction or a posteriori evaluation criterion of a given decomposition method [START_REF] Mcgill | Optimal resolution of superimposed action potentials[END_REF], [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF], [START_REF] Lewicki | A review of methods for spike sorting: The detection and classification of neural action potentials[END_REF], and is now introduced into the completely automated decomposition method through the use of prior laws on 1 x i . Although only general aspects of MU discharges are considered in this paper, further information could be integrated through the prior laws. The maximization of the posterior probability was achieved iteratively w.r.t. the continuous parameters and the NP-hard combinatorial discharge pattern. In addition, an extension of the segment-by-segment optimization method was introduced to treat the early convergence problem.

Some aspects of the method can be improved in future research. For example, the configurations of the Tabu algorithm remain an open problem, as to how to fix the FIFO memory length and the maximum number of iterations M , that in this paper has been empirically fixed to 20. Another issue concerns the design of neighborhood of a given state. The computational load for each iteration grows as more configurations are included in the neighborhood, with better chances of avoiding local optima. The search space may also be increased by including optimal alignment of the action potentials [START_REF] Mcgill | High-resolution alignment of sampled waveforms[END_REF] which is now limited by the sampling period. Finally, the proposed Bayesian decomposition approach is based on the modeling of a single-channel EMG recording. It is expected that a multichannel extension will improve the performance by exploiting the interchannel inference.

Concerning the statistical assumption on the independence of discharge patterns, we note that it does not constitute a limitation in the applicability of the approach but implies that the potential correlation between discharge patterns is not in-cluded in the prior laws. Thus, the method can still be applied to cases with high correlation between MU discharge patterns, such as tremor. In these cases, the correlation information is simply not exploited in the statistical model on which the method is based. In fact, the experimental signals were recorded from a muscle that presents relatively high degree of short-term synchronization between MUs [START_REF] Farina | Amplitude cancellation of motor-unit action potentials in the surface electromyogram can be estimated with spike-triggered averaging[END_REF]; thus, the assumption was not strictly valid in the analyzed experimental signals.

In conclusion, a new method for the decomposition of intramuscular EMG signals was proposed. It is fully automatic and provided an accuracy of approximately 90% on simulated signals and experimental single-channel recordings. The method does not make any assumptions on the particular shape of the action potentials or their spectral properties and can thus be applied to the spike classification problems of intraneural or intracortical recordings, for example.

Fig. 1 .

 1 Fig. 1. Decomposition example using one-segment maximization on simulated EMG data. (a) Simulated EMG recordings and reference discharge pattern. (b) Reconstructed EMG and discharge pattern estimated by the decomposition method. No errors are present in this decomposition example. (c) Estimates of the MUAP shapes. (d) Four MUAPs overlapped in the first segment analyzed in this example on an expanded time scale.

Fig. 2 .

 2 Fig. 2. Comparison of Tabu search within the one-and two-segment optimization approach. (a) Upper panel: EMG signal portion in which two segments are shown, delimited by square brackets. The reference MU labels are indicated with numbers as obtained by semiautomatic decomposition with an expert operator. Middle panel: The EMG signal reconstructed on the basis of the proposed automatic decomposition method applied to one-segment at a time. The estimatedMU labels are shown and it can be observed that MU labels #5 and #7 are switched. Lower panel: The reconstruction error. (b) Same as in (a) but with the application of the two-segment optimization strategy. The MU labels are in this case all correctly estimated. Lower reconstruction error is observed while the MU discharge patterns remain regular. (c) and (d) MUAP shapes of the switched pair {#5, #7}.

  (c) and (d) show similar energy level. A local maximum of (x i ) (k ) i includes MU#7 instead of MU#5 in Fig. 2(a) and this local configuration error is propagated to Seg k +1 since placing an MU#7 in x (k +1) 7

Fig. 3 .

 3 Fig. 3. (a) Upper panel: Segment of EMG and the MU labels decomposed by an expert operator. Middle panel: Signal reconstructed by the automatic decomposition method and the estimated MU labels. The proposed method estimated all labels correctly in this example. Lower panel: The reconstruction error. (b) Estimated MUAP shapes.

Fig. 4 .

 4 Fig. 4. (a) and (b) Instantaneous discharge rates obtained by the decomposition results on EMG of 10% MVC. (c) and (d) Corresponding MUAP shapes. (a) FP = 0, FN = 3 / Total = 171. (b) FP = 1, FN = 12 / Total = 137. (c) MUAP #1. (d) MUAP #2.

  , x

	(k,k+1) i,0 discharge in Seg k -1 and the first discharge in Seg k +2 respec-and x (k,k+1) i,n ( k,k+1) denotes the last +1 i tively. The neighborhood definition can be easily generalized
	(see definition in Section III-C) on a concatenated version
	(k,k+1) (x i ) i

TABLE I DECOMPOSITION

 I RESULTS ON SIMULATED EMG SIGNALS(10 S)

Authorized licensed use limited to: Aalborg Universitetsbibliotek. Downloaded on March 01,2010 at 09:31:04 EST from IEEE Xplore. Restrictions apply.


Pavel Labath

Joachim Niehren

A Uniform Programming Language for Implementing XML Standards ⋆

Keywords: Xml transformations, database queries, functional programming languages, compilers

We propose X-Fun, a core language for implementing various Xml standards in a uniform manner. X-Fun is a higher-order functional programming language for transforming data trees based on node selection queries. It can support the Xml data model and XPath queries as a special case. We present a lean operational semantics of X-Fun based on a typed lambda calculus that enables its in-memory implementation on top of any chosen path query evaluator. We also discuss compilers from Xslt, XQuery and XProc into X-Fun which cover the many details of these standardized languages. As a result, we obtain in-memory implementations of all these Xml standards with large coverage and high efficiency in a uniform manner from Saxon's XPath implementation.

Introduction

A major drawback of query-based functional languages with data trees so far is that they either have low coverage in theory and practice or no lean operational semantics. Theory driven languages are often based on some kind of macro tree transducers [START_REF] Neumann | Locating matches of tree patterns in forests[END_REF][START_REF] Hakuta | XQuery Streaming by Forest Transducers[END_REF][START_REF] Frisch | Streaming XML transformation using term rewriting[END_REF], which have low coverage, in that they are not closed under function composition [START_REF] Fülöp | Syntax-Directed Semantics -Formal Models based on Tree Transducers[END_REF] and thus not Turing complete (for instance type checking is decidable [START_REF] Maneth | XML type checking with macro tree transducers[END_REF]). The W3C standardised languages XQuery [START_REF] Robie | XQuery 3.0: An XML Query Language[END_REF] and Xslt [START_REF] Kay | XSL Transformations (XSLT) Version 3.0. W3C Last Call Working Draft[END_REF], in contrast, have large coverage in practice (string operations, data joins, arithmetics, aggregation, etc.) and in theory, since they are closed by function composition and indeed Turing complete [START_REF] Kepser | A simple proof for the Turing-Completeness of XSLT and XQuery[END_REF]. The definitions of these standards, however, consist of hundreds of pages of informal descriptions. They neither explain how to a build a compiler in a principled manner nor can they be used as a basis for formal analysis.

A second drawback is the tower of languages approach, adopted for standardised Xml processing languages. What happened in the case of Xml was the development of a separate language for each class of use cases, which all host the XPath language for querying data trees based on node navigation. Xslt serves for use cases with recursive document transformations such as Html publishing, while XQuery was developed for use cases in which Xml databases are queried. Since the combination of both is needed in most larger applications, the Xml pipeline language XProc [START_REF] Walsh | XProc: An XML Pipeline Language[END_REF][START_REF] Walsh | XML Calabash[END_REF][START_REF] Zergaoui | Innovimax: QuiXProc[END_REF] was developed and standardised again by the W3C. This resulted in yet another functional programming language for processing data trees based on XPath.

For resolving the above two drawbacks, the question is whether there exists a uniform core language for processing data trees that can cover the different Xml standards in a principled manner. It should have a lean and formal operational semantics, support node selection queries as with XPath and it should be sufficiently expressive in order to serve as a core language for implementing XQuery, Xslt, and XProc in a uniform manner.

Related work. An indicator for the existence of a uniform core language for Xml processing is that the omnipresent Saxon system [START_REF]Saxonica: SAXON 9.5: The XSLT and XQuery Processor[END_REF] implements Xslt and XQuery on a common platform. However, there is no formal description of this platform as a programming language, and it does not support the Xml pipeline language XProc so far. Instead, the existing implementations of XProc, Calabash [START_REF] Walsh | XML Calabash[END_REF] and QuiXProc [START_REF] Zergaoui | Innovimax: QuiXProc[END_REF], are based on Saxon's XPath engine directly.

The recent work from Castagna et. al. [START_REF] Castagna | A Core Calculus for XQuery 3.0[END_REF] gives further hope that our question will find a positive answer. They present an XPath-based functional programming language with a lean formal model based on the lambda calculus, which thus satisfies our first two conditions above and can serve as a core language for implementing a subset of XQuery 3.0. We believe that relevant parts of Xslt and XProc can also be compiled into this language, even though this is not shown there. The coverage, however, will remain limited, in particular on the XPath core (priority is given to strengthening type systems). Therefore, our last requirement is not satisfied.

Contributions.

In this paper, we present the first positive answer to the above question based on X-Fun. This is a new purely functional programming language. X-Fun is a higher-order language and it supports the evaluation of path-based queries that select nodes in data trees. The path queries are mapped to X-Fun expressions, whose values can be computed dynamically. In contrast to most previous interfaces between databases and programming languages, we overload variables of path queries with variables of X-Fun. In this manner, the variables in path queries are always bound to tree nodes, before the path query is evaluated itself. We note in particular, that path queries are not simply mapped to X-Fun expressions of type string.

The formal model of the operational semantics of X-Fun is a lambda calculus with a parallel call-by-value reduction strategy. Parallel evaluation is possible due to the absence of imperative data structures. The main novelty in X-Fun admission of tree nodes as values of type node. Which precise nodes are admitted depends on a tree store. New nodes can be created dynamically by adding new trees to the tree store. The same tree can be added twice to the store but with different nodes. How nodes are represented internally can be freely chosen by the X-Fun implementation and is hidden from the programmer. X-Fun can serve as a uniform core language for implementing XQuery, Xslt and XProc. In order to do so, we have developed compilers of all three languages into X-Fun. We also discuss how to implement X-Fun in an in-memory fashion on top of any in-memory XPath evaluator. Based on our compilers, we thus obtain new in-memory implementations of XQuery, Xslt and XProc with large coverage. Our implementation has very good efficiency and outperforms the most widely used XProc implementation by a wide margin.

Outline. In Section 2 we introduce our general model of data trees, alongside its application to Xml documents. The syntax and type system of the X-Fun language is introduced in Section 3. The applications of X-Fun to Xml document transformation is studied in Section 4, where we discuss compilers from other Xml processing languages into X-Fun. Section 5 contains our notes on the implementation of X-Fun and the results of our experiments.

Preliminaries

We introduce a general concept of data trees which will be used in the X-Fun language. We also show how to instantiate the trees to the Xml data model.

Data values and data trees

We fix a finite set Char whose elements will be called characters. A data value "c 1 • • • c n " is a word of characters for c 1 , . . . , c n ∈ Char . We define String = Char * to be the set of all data values, and nil="" to be the empty data value.

Next, we will fix a natural number k ≥ 1 and introduce data trees in which each node contains exactly k data values with characters in Char .

A node label is a k-tuple of data values, i.e., an element of (String) k . The set of data trees T of label size k over Char is the least set that contains all pairs of node labels and sequences of data trees in T . That is, it contains all unranked trees t with the abstract syntax t ::= l(t 1 , . . . , t n ), where n ≥ 0 and l ∈ String k . It should be noticed that the set of node labels is infinite, but that each node label can be represented finitely.

XML data model

For Xml data trees, we can fix k = 4 and Char the set of Unicode characters, and restrict ourselves to node labels of the following forms, where all v i are data values:

("element", v 1 , v 2 , nil) ("attribute", v 1 , v 2 , v 3 ) ("comment", nil, nil, v 3 ) ("processing-instruction", v 1 , nil, v 3 ) ("document", nil, nil, nil) ("text", v 1 , nil, nil)
An element ("element", v 1 , v 2 , nil) has three non-nil data values: its type "element", a name v 1 and a namespace v 2 . An attribute has four data values: its type, a name v 1 , a namespace v 2 , and the attribute value v 3 . A text node contains its type and its text value v 3 . Besides these, there are comments, processing instructions and the rooting document node.

Language X-Fun

In this section, we introduce X-Fun, a new functional programming language for transforming data trees. X-Fun can be applied to all kinds of data trees with a suitable choice of its parameters. We will instantiate the case of data trees satisfying the Xml data model concomitant with XPath as a query language.

We start with introducing the types and values of X-Fun (Section 3.1). Then we explain how to map path queries to X-Fun values, by using particular X-Fun expressions with variables (Section 3.2). The general syntax of X-Fun expressions is given in Section 3.3. Some syntactic sugar and an example of an X-Fun program are given in Sections 3.5 and 3.6. Discussion of the typing rules for X-Fun's type system and the formal semantics of X-Fun can be found in the research report [START_REF] Labath | A Uniform Programming Language for Implementing XML Standards[END_REF].

Types and Values

The X-Fun language supports higher-order values and expressions with the following types:

T ::= none | node | tree | number | bool | char | T 1 × . . . × T n | [T ] | T 1 → T 2 | T 1 ∪ T 2
A value of type char is an element of Char , a value of type tree is an element of T . A value of type number is a floating point number, while the values of type bool are the Boolean values true and false. A value of type node will be a node of the graph of one of the trees stored by the environment of the X-Fun evaluator. The precise node identifiers chosen by the evaluator are left internal (to the mapping from trees to graphs).

As usual, we support list types [T ] which denote all lists of values of type T , product types T 1 × . . . × T n whose values are all tuples of the values of types T i , and function types T 1 → T 2 whose values are all partial functions of values of type T 1 to values of type T 2 . Besides these, we also support type unions in the obvious manner.

A data value "c Since XPath can return sequences of items of different types, we define the type pathresult as node ∪ number ∪ string ∪ bool. The result of evaluating a path expression will then be of type [pathresult]. To be able to specify path expressions, we define the type path as [char ∪ pathresult ∪ [pathresult]], i.e., as list of characters, individual items returned by a path expression, and whole sequences of those items.

1 • • • c n " ∈ String is considered as a list of characters of type string = [char]. A node label is considered a k-tuple of strings, i.e.,

XPath queries as X-Fun expressions

We will consider XPath expressions as values of our programming language. This is done in such a manner that the variables in XPath expressions can be bound to values of the programming language. For instance, if we have an XPath expression $x//book[auth=$y] then one might want to evaluate this expression while variable x is bound to a node of some tree and variable y to some data value. In X-Fun, the above query will be represented by the following expression of type path, where x is a variable of type node and y a variable of type string:

x :: '/' :: '/' :: 'b' :: 'o' :: 'o' :: 'k' :: '[' :: 'a' :: 'u' :: 't' :: 'h' :: '=' :: y :: ']' :: nil

The concrete syntax of X-Fun supports syntactic sugar for values of type path, so that the above expression can be defined as:

"$x//book[auth=$y]"
In order to enable the evaluation of path expressions, X-Fun supports a builtin function evalPath of type path → [pathresult]. In an implementation of X-Fun, this function can be mapped straightforwardly to existing XPath evaluators.

Syntax of X-Fun expressions

X-Fun is a purely functional programming language whose values subsume higherorder function, trees, strings, numbers and Boolean values. The evaluation strategy of X-Fun is fully parallel, which is possible since no imperative constructs are permitted.

The syntax of X-Fun programs E is given in Figure 1. All expressions of X-Fun are standard in functional programming languages, so we only briefly describe different kinds of subexpressions of X-Fun programs.

A variable x is evaluated to the value of the corresponding type. The constant expression c returns the respective constant, which can be a Boolean value, a number or a character from Char . The list constructor E 1 :: E 2 prepends an element to a list, while the tuple constructor (E 1 , . . . , E n ) constructs tuples.

The match expression match E { P 1 → E 1 , . . . , P n → E n } selects one of the branches E i based on the patterns P i , which are matched against the value of E. The pattern x : T captures a matched value of type T into a variable. The pattern !(E) matches the value against the value of expression E. Here, the matching of functional values, or lists/tuples that contain functions is not permitted. Pattern P 1 :: P 2 matches a list if P 1 and P 2 match its head and tail, while the pattern (P 1 , . . . , P n ) matches tuples.

A function expression fun x:T 1 → T 2 { E } returns a new function, with the argument x : T 1 and the return value of type T 2 obtained by the evaluation of the function body E. The expression E 1 (E 2 ) applies a function to a value. X-Fun also supports exception handling, where exceptions are values of type string.

Expressions

E ::= x | c | E 1 :: E 2 | (E 1 , . . . , E n ), n ≥ 2 | match E { P 1 → E 1 , . . . , P n → E n } | fun x:T 1 → T 2 { E } | E 1 (E 2 ) | try E 1 catch(x) E 2 | raise(E) Patterns P ::= x : T | !(E) | P 1 :: P 2 | (P 1 , . . . , P n ), n ≥ 2
Fig. 1. Syntax of X-Fun's expressions

Builtin operators

At the beginning of the evaluation, the environment contains bindings of the global variables given in Figure 2. The first block contains three functions, whose semantics are parameters of the language, and depend on the query language and data model. For a label l and a sequence of trees h, the function application makeTree(l, h) returns the data tree l(h), if l(h) is a well-formed data tree (e.g., in the Xml data model attributes cannot have children) and raises an exception otherwise. The function evalPath(p) evaluates a path expression p. Whenever p is not well-formed (e.g., with respect to the XPath 3.0 specification) an error is raised. Note that path expressions are X-Fun values, which means they can be computed dynamically by the X-Fun program using information from the input data tree. We will also define functions evalPath T , on top of evalPath, for T = [node], [string], etc. These functions verify (using a match expression with a typecase) that the result of the path call is of type T and raise an exception otherwise.

Parameters

The next four operators are generic and do not depend on the specific kind of data trees. The variable nil refers to the empty list. A function application subtree(v) returns the subtree rooted at node v, while a function application label(v) returns the label of the node. The function addTree returns the identifier of the root node of the tree, and is used for storing the graph of the tree in the environment. This function can be used to access nodes of newly generated trees by starting path navigation from their root.

Syntactic sugar

In the X-Fun snippets in the rest of the paper we shall employ some syntactic shortcuts, which enable us to express more succinctly some X-Fun constructs: list concatenation We shall use the binary operator * to concatenate two lists. simplified patterns When the type of a capture variable can be deduced from the matched expression we shall omit the ": T " in the capture pattern. This happens when the match expression is used to decompose lists and tuples instead of doing a typecase. For example, we shall simply write match E { h :: t → E 1 , e → E 2 } to get the head and tail of a list. let-declarations We shall use the syntax let . This is unambiguous since tuples always have at least two members.

x 1 = E 1 , . . . , x n = E n in E instead of match (E 1 , . . . , E n ) { (x 1 , . . . , x n ) → E }

Example

In Figure 3 we illustrate a transformation that converts an address book into Html. The address fields are assumed to be unordered in the input data tree, while the fields of the output Html addresses should be published in the order name, street, city and, phone. An X-Fun program defining this transformation is given in Figure 4. Starting at the root it first locates all address records, and applies the function convert address to each of them. For each address record, the program first extracts the values of the fields name, street, and city located at some children of x. These values are then bound to variables named alike and later output as text nodes. The example program uses the standard map function, which can be defined in X-Fun for every T and T ′ as follows

map T →T ′ = fun x : (T → T ′ ) × [T ] → [T ′ ] { match x { ( f , head : : t a i l ) → f ( head ) : : map T →T ′ ( f , t a i l ) o t h e r → nil } }
and the functions element and text, which are wrappers around makeTree which facilitate creation of nodes of the correct kind. ( "$x/ c h i l d : : name/ t e x t ( ) " ) , s t r e e t = evalPath [node] ( "$x/ c h i l d : : s t r e e t / t e x t ( ) " ) , c i t y = evalPath [node] ( "$x/ c h i l d : : c i t y / t e x t ( ) " ) i n element ( " l i " , element ( "p" , map node→tree ( name , subtree ) ) : : element ( "p" , map node→tree ( s t r e e t , subtree ) ) : : element ( "p" , map node→tree ( c i t y , subtree ) ) : : map string→tree ( fun x : string → tree { element ( "p" , text ( "Phone : " * x ) : : n i l ) } , evalPath In this section, we briefly sketch translations from the standard Xml processing languages, Xslt XQuery and XProc. A more thorough treatment of this topic can be found in the full version of the paper. [START_REF] Labath | A Uniform Programming Language for Implementing XML Standards[END_REF] By implementing these three compilers, we obtain a uniform implementation of the whole Xml processing stack based on a single X-Fun evaluator.

Xslt. Each template in the Xslt stylesheet is translated to a function in X-Fun. Furthermore, for each mode, we produce an additional function which implements the selection of the correct template from the set of templates associated with that mode according to their match patterns. The call-template and apply-templates instructions are translated as calls to the template or mode functions respectively. In the copy-of instruction, the nodes returned by the XPath expression are copied to the output using the subtree function and strings and numbers are converted to a new text node with a call to makeTree. The instructions constructing elements, attributes and other Xml nodes translate to corresponding calls to makeTree. The for-each instruction translates to a call to map, where the list to map over is produced by a call to evalPath and the mapping function is the body of the for-each instruction. Other Xslt instructions like if and choose can be translated similarly.

XQuery. The feature that most distinguishes XQuery is the Sql-like Flwor expression. It enables the programmer to create a stream of tuples using the for and let clauses, filter them with a where clause and then reorder them using the order by clause. There is no single expression in X-Fun which covers this functionality, but it is easy to build it piecewise. Using several evalPath calls we can construct the list of tuples which corresponds to the tuple stream of XQuery. Sorting and filtering of a list are functions easily definable in a functional language, and the functionality of where and order by is translated to calls to these functions. The sort and filter conditions are given again by calls to evalPath with the appropriate XPath expression. Translation of other XQuery constructs like the if expressions and functions proceeds in a straight forward manner.

XProc. By encapsulating each processing step in a function, X-Fun can easily express the multi-stage processing which is inherent in XProc. The pipelines then become simple function compositions. XProc steps which invoke XQuery or Xslt processing are handled by defining a function whose body is the translation of the respective program. Simple XProc steps like split-sequence, which splits a sequence of documents into two based on an XPath criterion are defined as normal X-Fun functions and provided as a library. The pipeline them simply calls these functions to do the required processing. The rest of the constructs like choosing among alternative subpipelines (choose) or looping over documents in a sequence are compiled to match and map expressions in X-Fun.

Implementation and Experiments

We have implemented a proof-of-concept X-Fun language evaluator in the Java programming language. We have instantiated X-Fun with the Xml data model, using standard Java libraries for manipulating Xml trees. We have used XPath as the path language, as implemented by Saxon. We have used standard techniques for implementing functional languages, using the heap to store the values and the environment of the program and a stack for representing recursive function calls. We reduce an expression in all possible positions in an arbitrary order.

We have attempted to interface our implementation with Tatoo, a highly efficient evaluator of an XPath fragment based on [START_REF] Arroyuelo | Fast in-memory xpath search using compressed indexes[END_REF]. Unfortunately, the penalty of crossing the language barrier (Tatoo is implemented in OCaml) shadowed all performance gains from a faster implementation, so we could not perform any significant experiments. To see the difference in performance in using a faster XPath implementation, we would need to implement X-Fun in OCaml as well.

We have also implemented the compilers of Xslt and XQuery into X-Fun. In order to support real-world Xslt and XQuery, they need support for additional features, like modules and various optional attributes of expressions in these languages (e.g., grouping with the group-starting-with attribute, etc.). However, none of these limitations are fundamental and they are not implemented because of their volume. The supported fragment is wide enough to run all queries from the XMark [START_REF] Schmidt | XMark: A benchmark for XML data management[END_REF] benchmark.

We don't have an XProc compiler implementation, but for the purposes of testing we have run X-Fun on manually translated programs.

Experiments

To evaluate the performance of our implementation, we have compared it with the leading industry tool, the Saxon Xslt and XQuery processor. To compare our performance on XProc pipelines, we have used Calabash, the most frequently used XProc processor, as baseline. The tests were run on a computer with an Intel Core i7 processor running at 2.8 GHz, with 4GB of RAM and a SATA hard drive, running 64-bit Linux operating system.

First, we have compared the running time of our implementation on XQuery programs. We used the queries from the XMark benchmark, and the results are in Figure 5. The tests show that the running time of both tools is comparable. X-Fun is faster in case of simple queries (Q6, Q7, Q15, which contain just a simple loop), while Saxon is faster on queries involving joins (e.g., Q8, Q9, Q11). On the rest of queries our implementation of X-Fun is at most 20% slower that the competition, which we consider a good result as Saxon is a highly optimised industry tool, while we have not spent much time optimising the performance of our X-Fun implementation.

For the Xslt test, we used a transformation publishing an address book to Html. The transformation in question is a more elaborate version of the program in Figure 4, and it includes about 40 XPath expressions. The tests show that Saxon is about 4 times faster than our tool (for example, 15.7 vs. 63 seconds on a 200 MB document) and that the time of both tools scales linearly with the document size.

In the XProc comparison, we have a simple pipeline consisting of 4 steps. First, it selects subtrees from the input document, splits the resulting sequence into two based on the presence of some node. The documents from the two sequences are then joined into pairs and these pairs are concatenated to form a single document again. We have compared the performance of Calabash with Query X-Fun Saxon Q1 13.5 our implementation of the pipeline in X-Fun. Both implementations show linear scalability with respect to size of the input and the pipeline, as can be seen in Figures 6 and7 (for scaling the pipeline size, we simply composed the described pipeline with itself). However, our own implementation is consistently at least two times faster, and for the larger pipelines the difference is even more apparent. While the relatively low processing speed per megabyte can be explained by the need to create many small documents (the element per megabyte density is much higher compared to the previous tests), it is surprising to see an implementation specifically designed for processing XProc be outperformed by our unoptimised implementation of the pipeline steps. We have presented X-Fun, a language for processing data trees and shown that can serve as a uniform programming language for Xml processing and as a uniform core language for implementing XQuery, Xslt, and XProc on top of any existing XPath evaluator. Our implementation based on Saxon's inmemory XPath evaluator yields surprisingly efficient implementations of the three W3C standards, even there is a lot of space left for optimisation. We have obtained results which are a match for the Saxon's XQuery and Xslt evaluators and in the case of XProc, first results show that we are already faster than Calabash.

Our prime objective in future is to build streaming implementations of X-Fun, and thus of XQuery, Xslt, and XProc. The main ideas behind it are described in a technical report [START_REF] Labath | A Functional Language for Hyperstreaming XSLT[END_REF]. These streaming implementation will serve in the tools called QuiXQuery, QuiXslt, and QuiXProc. A first version of QuiXslt is freely available for testing on our online demo machine [START_REF]Quix tools suite[END_REF] while streaming is not yet available for our current QuiXProc implementation.

  as a value of type label = string k . Hedges are considered as lists of trees of type hedge = [tree].

Fig. 2 .

 2 Fig. 2. Builtin operators of X-Fun

  as a more familiar way to declare variables. tuple arguments We shall allow tuple arguments to functions to be written without an extra pair of parentheses. I.e., f (a, b) instead of f ((a, b))

Fig. 3 .

 3 Fig. 3. Publication of an address book in Html except for secret entries

  fun book : t r e e → t r e e { l e t b o o k r o o t = addTree ( book ) i n l e t c o n v e r t a d d r e s s = fun x : node→ t r e e { l e t name = evalPath [node]

Fig. 4 .

 4 Fig. 4. X-Fun program converting address books to Html

Fig. 5 .

 5 Fig.5. Running time in seconds of X-Fun and Saxon on queries from the XMark benchmark on a 500 MB document. The three queries marked with '*', due to their complexity, were run on a 300 MB document.

Fig. 7 .

 7 Fig. 7. Performance of X-Fun and Calabash on a 2 MB document with varying pipeline size

Acknowledgement. We would like to thank Guiseppe Castagna and Kim Nguyen for their helpful discussions about the type system of X-Fun.

⋆ This research was supported in part by the grant VEGA 1/0979/12