Jose E Morales

Guillaume James

Arnaud Tonnelier

Solitary waves in the excitable Burridge-Knopoff model

I. Introduction

• The Burridge-Knopoff model (BK) has been originally introduced to capture certain nonlinear mechanisms of earthquake fault dynamics through a simple nonlinear lattice model. The BK model describes more generally the interaction of two elastic media under frictional contact, one being pulled at constant speed V . This medium is discretized as a chain of blocks (see Fig. Each block is subject to a friction force F (v), where v denotes the sliding velocity.

• Friction experiments with metals, rocks, hydrogels, polymers suggest the existence of non-monotonic friction forces increasing at high velocities (see Fig. 2).

• For non-monotonic friction forces, the BK model behaves like an elastic excitable medium and supports travelling waves.

• Overall purpose: Study the existence of solitary waves (i.e. localized travelling waves) depending on model parameters.

II. The model

• The BK model is expressed by the following equations:

γ du n dt = k∆y n -F (V + u n ) -y n , dy n dt = u n
where ∆y n = y n+1 -2y n + y n-1 denotes the discrete Laplacian, γ a mass parameter, k the coupling parameter, F (v) the friction function, u n the velocity variable and y n the position variable.

• A class of non-monotonic friction forces was explored to probe solitary wave existence. Fig. 3 and Fig. 4 present two of the main friction functions considered, a cubic friction function and an idealized piecewise linear friction function, respectively. • Our interest is to find solutions u n (t) = ϕ(n -ct) and y n (t) = ψ(n -ct)

F (v) v 1 1 -α b a F 0 (v) v 1 1 -α a
that describe solitary waves. The functions ϕ, ψ satisfy the advance-delay differential equation:

-cγϕ (ξ) = k [ψ(ξ + 1) -2ψ(ξ) + ψ(ξ -1)] -F (V + ϕ(ξ)) -ψ(ξ), -c ψ (ξ) = ϕ(ξ).

III. Travelling pulse solutions

The piecewise linear friction function in Fig. 4 is defined by

F (v) = v a -αH(v -a)
where H(x) is the Heaviside function. For this friction force, we look for a solution ϕ(ξ) that satisfies:

(H) :

       ϕ(0) = ϕ(ξ 1 ) = -β ϕ(ξ) > -β if ξ ∈] -∞, ξ 1 [∪]0, +∞[ ϕ(ξ) < -β if ξ ∈]ξ 1 , 0[, where β = V -a and |ξ 1 | = solitary wave width.
The resulting linear advance-delay differential equation

c 2 γϕ - c a ϕ + ϕ = k∆ϕ + ac(δ ξ 1 -δ 0 )
can be solved by Fourier transform:

ϕ(ξ; c, ξ 1 , k) = R e i2πλξ g c,ξ 1 ,k (λ)dλ,
where:

g c,ξ 1 ,k (λ) = αc(1 -e -2πiλξ 1 ) 4π 2 c 2 x 2 γ -4k sin 2 (πx) + 2πicx a -1 -1
We solve system (H) in the overdamped regime, for a large enough wave velocity c, small coupling constant k and V ≈ a. This leads to the following:

Existence theorem for solitary waves

Assume γ < 1 4a 2 . Fix c > b aγ ln 1+b 1-b with b = 1 -4a 2 γ.
Then for all k small enough, for V = V (k) = a + O(k) > a, ∃ a solitary wave solution of the BK model taking the following form when k → 0: • A solitary wave solution is shown in Fig. 6 supported by the existence theorem.

ϕ(ξ) = αc[K(ξ -ξ 1 ) -K(ξ)] + O(k)
• Solitary waves satisfying (H) exist only above a minimal coupling strength k min = O(V -a). Fig. 7 presents the coupling strength k as a function of the wave speed c of the solitary wave.

IV. Dynamical simulations

• The propagation of a stable solitary wave is shown in Fig. 8. There exist solitary waves with oscillatory tails for V → a as shown in Fig. 9. Certain nonmonotonic friction forces also lead to bistability between the ground state and periodic oscillations, see Fig.

V. Conclusions

• In the excitable Burridge Knopoff model, we report for the first time the existence of solitary waves for a wide range of non-monotonic cubic-like friction forces and parameters. For a piecewise-linear and discontinuous friction force, we proved existence of solitary waves at small coupling k for V ≈ a.

• Besides the classical earthquake fault interpretation, the excitable BK equations may find interesting applications for the description of friction at the interface of patterned surfaces, where the blocks correspond to discrete surface elements. Other applications involve mechanical oscillator networks and nonlinear electric transmission lines.

Fig. 1

 1 Fig.1 Schematic view of the Burridge-Knopoff model

Fig. 2 :

 2 Fig. 2: Friction of a gel

Fig. 3 :

 3 Fig. 3: Smooth friction

  Fig. 5: Shape of function K(ξ)

Fig. 8 :

 8 Fig. 8: Velocities u n (t).Fig. 9: Velocity ϕ(n -ct) (solid), position ψ(n -ct) (dashed). Fig. 10: Velocities u n (t). Fig. 11: Velocity u n (t) at n = 17500


A Limem

G Delmaire
email: gilles.delmaire@lisic.univ-littoral.fr

M Puigt

G Roussel
email: gilles.roussel@lisic.univ-littoral.fr

D Courcot

Non-negative matrix factorization under equality constraints -a study of industrial source identification

Keywords: Non-negative Matrix Factorization, Quadratic Optimization, Air quality

This work is devoted to the factorization of an observation matrix into additive factors, respectively a contribution matrix G and a profile matrix F which enable to identify many pollution sources. The search for G and F is achieved through Non-negative Matrix Factorization techniques which alternatively look for the best updates on G and F .

These methods are sensitive to noise and initialization, and-as for any blind source separation method-give results up to a scaling factor and a permutation. A Weighted Non-negative Matrix Factorization extension has also been proposed in the literature, so that different standard deviations of the data matrix components are taken into account. However, some estimated profile components may be inconsistent with practical experience. To prevent this issue, we propose an informed Non-negative Matrix Factorization, where some components of the profile matrix are set to zero or to a constant positive value. A special parametrization of the profile matrix is developed

Introduction

Approximate matrix factorization techniques are of great importance in several areas, such as pattern recognition, antenna array processing, or environmental data processing. Non-negative matrix Factorization (NMF) aims to factorize a non-negative matrix as the product between two matrices, such that each entry of the latter matrices is null or positive. It emerged under the name Positive Matrix Factorization (PMF) with Paatero and Tapper's contributions [START_REF] Paatero | Least squares formulation of robust non negative factor analysis[END_REF][START_REF] Paatero | Positive matrix factorization: a non negative factor model with optimal utilization of error estimates of data values[END_REF]. It then appeared in the field of signal and image processing with the work by Lee and Seung [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. Some of the most well known approaches have been proposed by Cichocki and Zdunek [START_REF] Cichocki | NMFLAB for signal processing[END_REF] and Lin [START_REF] Lin | Projected gradients methods for non-negative matrix factorization[END_REF], and are based on a projected gradient method. Recently, some weighted versions of NMF were investigated and took into account specific standard deviations associated to each data point. Inconsistencies in practical results lead to investigate additional assumptions such as sparseness constraints [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraint[END_REF][START_REF] Berry | Algorithms and applications for approximate non-negative matrix factorization[END_REF], fixed row and column sums constraints [START_REF] Ho | Non-negative matrix factorization with fixed row and column sums[END_REF], or orthogonality constraints [START_REF] Yoo | Non-negative matrix factorization with orthogonality constraints[END_REF][START_REF] Choi | Algorithms for orthogonal non-negative matrix factorization[END_REF].

Over the last decade, matrix factorization has been extensively investigated in the field of atmospheric sciences, and particularly for receptor models with the aim to re-construct the impacts of emissions from different pollutant sources, e.g., airborne particles [START_REF] Hopke | PM source apportionment and health effects: Intercomparison of source apportionment results[END_REF]. Such source apportionment studies enable to identify the main emission sources and their relative contribution at different monitoring sites. However, Viana et al. [START_REF] Viana | The effects of particle-size distribution and chloride depletion of sea-salt aerosols on estimating atmospheric deposition at a coastal site[END_REF] reported that, depending on the applied NMF method, these techniques may provide highly different solutions. Moreover-and to the best of our knowledge-NMF tools fail in discriminating some natural sources (i.e., crustal matter from local or regional re-suspension and road dust), as well as sources containing sulfate and nitrate components. Moreover, the separation of mineral particles emitted from both industrial dust emissions and natural sources appears as a difficult task [START_REF] Escrig | Application of optimally scaled target factor analysis for assessing source contribution of ambient PM10[END_REF]. In this paper, we propose to take into account some partial a priori knowledge on some source profiles in order to improve the estimation to the unknown components.

The remainder of the paper is organized as follows. In Section 2, we present the concepts of some state-of-the-art NMF methods. Section 3 introduces our proposed approach while its performance is investigated in Section 4. We conclude and discuss future directions of this work in Section 5.

A short review of Non-negative Matrix Factorization

NMF aims to estimate two non-negative matrices whose product approximates the observed one. The non-negativity constraints make the decomposition purely additive whereas other ones like Principal Components Analysis (PCA) and Independent Components Analysis (ICA) can provide negative components. Informed NMF is located between NMF methods-which only assume non-negativity of the matrices-on one hand, and on the other hand, regression models which assume that one matrix is exactly known. In this paper, we introduce an informed NMF method, since matrices may be partially known.

Receptor modelling

This model is very generic since it is often used in signal and image processing. In environmental field, the factorization model (1) below is best known as receptor modelling and links the data matrix and the active sources. The factorization enables to approximate the n × m data matrix X by the product of two matrices,

X ≈ G • F, (1)
where

• X is the n × m data matrix, where n is the number of samples and m is the number of species. In environmetrics, x ij -the (i, j) th element of X-accounts for the concentration, expressed in ng/m 3 , of the j th chemical species coming from the i th sample. A row gathers concentrations of all registered species for the current sample.

• G stands for the n × p contribution matrix, where p is the pre specified number of sources and n is the number of samples. In environmetrics and chemometrics, the (i, k) th element of G-denoted g ik -is referred to as the massic contribution from Source k to Sample i, expressed in µg/m 3 .

• F is a p × m matrix of profiles. Its (k, j) th term f kj stands for a mass ratio (in ng/µg) corresponding to a percentage of the j th species with respect to the whole mass of the Source k.

The number p of sources is usually chosen such that np + pm nm. The a priori choice for the appropriate value of p mostly depends on the nature of the investigated data [START_REF] Guillamet | Introducing a weighted non-negative matrix factorization for image classification[END_REF] and on the expert's knowledge of the number of sources. In the absence of such information, p can be estimated in a preprocessing stage by applying some information theory techniques [START_REF] Fishler | Detection of signals by information theoretic criteria: general asymptotic performance analysis[END_REF] such as AIC or BIC or extended versions.

Non-negative matrix factorization

Non-negative matrix factorization has received a lot of attention since Lee and Seung [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] published their popular multiplicative algorithm. In fact, NMF looks for an approximate factorization of the data matrix according to receptor modelling (1) under the additional non-negativity constraint. The basic NMF solves min

G,F ≥0 n i=1 m j=1 ((x ij -(GF ) ij )) 2 = min G,F ≥0 X -G • F 2 f , (2)
where G, F ≥ 0 means that all components of G, F are non-negative and . f stands for the Frobenius norm1 . The problem in ( 2) is convex with respect to G and F separately, but not for both at the same time. As a consequence, most algorithms are only able to yield local minima. Classical NMF algorithms thus alternatively search for F (respectively G) when G is fixed (respectively F ). As an alternative, Kim et al. proposed a general framework [START_REF] Kim | Algorithms for nonnegative matrix and tensor factorizations: a unified view based on block coordinate descent framework[END_REF], called Block Coordinate Descent, which includes several methods described hereafter. Different iterative approaches to solve (2) were proposed in the literature:

• Most commonly used ones are the multiplicative updates rules. Lee and Seung [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] mentioned that this update is a good trade-off between speed of convergence and ease of implementation. The update of F (respectively G) is made through a multiplicative update coefficient which depends on the data matrix X and the contribution matrix G (respectively the profile matrix F ):

F ← F • (G T X) (G T GF ) , G ← G • (XF T ) (GF F T ) , (3)
where X • Y and X Y respectively account for component-wise product and element-wise division between two matrices. According to the update rule (3), the Frobenius norm is not increasing and becomes constant if a limit point is reached.

A limit point (G, F ) is said to be stationary if it checks the Karush-Kuhn Tucker (KKT) conditions, i.e.,

F ≥ 0, G ≥ 0, G T • (G • F -X) ≥ 0, (G • F -X) • F T ≥ 0, F • G T • (G • F -X) = 0, G • (G • F -X)) • F T = 0.
Stationarity is only a necessary condition of a local minimum. Moreover, Lin [START_REF] Lin | On the convergence of multiplicative update algorithms for non-negative matrix factorization[END_REF] reported that some limit points of (3) which are not stationary may exist, especially if some components of F and G are initialized to zero.

• As an alternative, Alternating Non-negative Least Squares (ANLS) consists of alternatively estimating G (respectively F ) by setting F (respectively G). The mathematic formulation of the problem at iteration r + 1 may be divided into two sub-problems depending on the results obtained at iteration r [START_REF] Kim | Nonnegative matrix factorization based on alternating nonnegativity constrained least squares and active set method[END_REF]:

F r+1 = arg min F ≥0 X -G r F 2 f , G r+1 = arg min G≥0 X -GF r+1 2 f . (4)
Each step of (4) may be seen as a set of independent non-negative least squares problems involving one column or one row of the data matrix. Let x i be the i th column of the matrix X and f r+1 j be the j th column of the matrix F r+1 . Searching for F amounts to find separately each column of F ,

f r+1 j = arg min f j ≥0 x i -G r f j 2 f ∀j ∈ {1...m} ∀i ∈ {1...m}.
The resolution of this problem is generally more time consuming at each iteration than the multiplicative updates but the number of iterations is expected to be lower. This method may also be considered as a special case of the Block Coordinate Descent framework in the case of two blocks [START_REF] Kim | Algorithms for nonnegative matrix and tensor factorizations: a unified view based on block coordinate descent framework[END_REF]. Grippo and Sciandrone [START_REF] Grippo | On the convergence of the block nonlinear gauss-seidel method under convex constraints[END_REF] state that every limit point obtained from ( 4) is a stationary point.

• Lin [START_REF] Lin | Projected gradients methods for non-negative matrix factorization[END_REF] proposed to develop projected gradient-based techniques for bound constrained optimization. The iterative update of the matrix is done through a projected gradient approach on the feasible region:

F r+1 ← P F r -α r ∇ F r X -G r F 2 f ,
where the projection P enables to map the candidate back into the bounded region if it lies outside:

∀(i, j) ∈ {1...n}×{1...m}, (P (X)) i,j =    X i,j if m i,j ≤ X i,j ≤ M i,j , M i,j if X i,j ≥ M i,j , m i,j if X i,j ≤ m i,j .
Many articles-such as [START_REF] Bertsekas | Nonlinear Programming[END_REF] for example-are devoted to the research of a step α r which checks a sufficiently fast decreasing for each step. This stage is also the most time consuming one.

Even if most work-including the above methods and our contributionfocus on approximated solutions of the factorization of X, please note the existence of recent investigations on exact factorization. For example, Vavasis [START_REF] Vavasis | On the complexity of nonnegative matrix factorization[END_REF] showed the NP-hardness of computing such a factorization. Recently, Arora et al. proposed a near-separable NMF algorithm which provably runs in polynomial time under some additional assumptions on the input matrix X [START_REF] Arora | Computing a nonnegative matrix factorization -provably[END_REF]. Near-separability has been furtherly explored in, e.g., [START_REF] Bittorf | Factoring nonnegative matrices with linear programs[END_REF][START_REF] Gillis | Fast and robust recursive algorithms for separable nonnegative matrix factorization[END_REF].

Weighting the Non-negative Matrix Factorization

Recently, some authors have proposed an uncertainty weighted version of NMF, expressed either as a sequence of columns updates [START_REF] Delmaire | Une version pondérée de la factorisation matricielle non négative pour l'identification de sources de poussières. Application au littoral de la mer du nord[END_REF] or as a whole update [START_REF] Ho | Non-negative matrix factorizations algorithms and applications[END_REF]. The function to be minimized derives from the classical NMF cost function and involves a weight matrix, i.e., min

G,F ≥0 n i=1 m j=1 (x ij -(GF ) ij ) σ ij 2 ,
where σ ij stands for the standard deviation obtained from the measurement of the i th sample and the j th species, gathered in the uncertainty matrix Σ = {σ ij }. A whole weight matrix W may be expressed:

W = 1 n×m Σ•Σ
, where each element of the matrix 1 n×m is equal to 1.

Ho [START_REF] Ho | Non-negative matrix factorizations algorithms and applications[END_REF] derived the Karush-Kuhn Tucker conditions in the weighted case, which are outlined as

F ≥ 0, G ≥ 0, G T • (W • (G • F -X)) ≥ 0, (W • (G • F -X)) • F T ≥ 0, F • G T • (W • (G • F -X)) = 0, G • (W • (G • F -X)) • F T = 0.
Ho also introduced a unified expression of an iterative solution of this criterion using the maximization minimization theorem [START_REF] Ho | Non-negative matrix factorizations algorithms and applications[END_REF]. The main advantage lies in the whole iterative update of each iteration, i.e.,

F ← F • G T (W • X) G T (W • (GF )) , G ← G • (W • X)F T (W • (GF ))F T , (5)
In order to prevent division by 0, a small quantity ε is added to the denominator of the relation (5), i.e.,

F ← F • G T (W • X) G T (W • (GF )) + ε1 p×m , G ← G • (W • X)F T (W • (GF ))F T + ε1 n×p .
Using weights is particularly important in the case of physical measurements, where data are obtained with radically different uncertainties. Note that the update rules in (5) may be viewed as an extension of the classical multiplicative update described in (3) to weighted data. Indeed, by setting W to the matrix 1 n×m , the updates rules in (3) are exactly those in [START_REF] Choi | Algorithms for orthogonal non-negative matrix factorization[END_REF]. From this remark, weighted NMF (WNMF) behaves as the classical NMF for convergence aspects. In fact, some non-stationary limit points may exist with weighted multiplicative updates, especially in the case when some entries of F and G are initialized to zero [START_REF] Lin | On the convergence of multiplicative update algorithms for non-negative matrix factorization[END_REF].

Other criteria may also be investigated, notably the minimization of the weighted Kullback-Leibler divergence. For this case, the expression of weighted multiplicative updates is available in [START_REF] Ho | Non-negative matrix factorizations algorithms and applications[END_REF]. To conclude, the main interest for this kind of algorithm lies in the direct update of both matrices and its ease of implementation.

Solving a constrained problem

Practically, the profiles recovery is never completely blind. In some cases, a priori information on components are available and some values can be set to zero if some species are absent from a source profile. This kind of knowledge should be included as constraints in our algorithms. Up to now-and to the best of our knowledge-no state-of-the-art method supports weighted cost functions with constraints.

Introduction of equality constraints

In many cases, some known components of the source profile may be provided by experts. The aim of this section is to find a parametrization that naturally integrates this knowledge. Constraints may be formulated as follows. Let Ω be a p × m binary matrix which specifies the presence or the absence of constraints on each element of the profile matrix F ,

∀(i, j) ∈ {1...p} × {1...m}, Ω ij = 1 if F ij has to be set, 0 otherwise.
Φ is the p × m matrix of set values. By construction, Φ ij is equal to zero when Ω ij = 0. It may thus be easily checked that

Φ • Ω = Φ . ( 6
)
Using the definition of matrices (Φ, Ω), known components of F may easily be set to zero or to other positive values according to the following relationship:

F • Ω -Φ = 0 . (7)
Let:

• f i be the i th column of the F matrix,

• ϕ i be the i th column of the Φ matrix,

• M i be the (l i × p) constraint matrix made up with the i th column of the matrix Ω (i.e., containing l i constraints),

• δ i be the vector of the constraint values extracted from Φ.

Constraints may be expressed for the i th column vector of the profile matrix as

M i f i = M i ϕ i = δ i , ∀i ∈ {1...m}. (8)
As an example, let us consider a case with, say p = 5 sources and l i = 2 constraints for all i. The latter are here assumed to be set on the second and the fourth position, and to be respectively equal to, say 80 and 30. The previous matrices then read:

M i = 0 1 0 0 0 0 0 0 1 0 , ϕ T i = 0 80 0 30 0 , δ T i = 80 30 , ∀i ∈ {1...m}.
Let Γ i be the (p × (p -l i )) matrix whose columns span the supplementary space to the rows of M i such that rank M i Γ T i = p. Γ i checks the following normalization relation

M i Γ i = 0 l i ×(p-l i ) , Γ T i Γ i = I (p-l i )×(p-l i ) , ∀i ∈ {1...m},
where 0 l i ×(p-l i ) and I (p-l i )×(p-l i ) are respectively the null matrix and the identity matrix. By noticing in ( 8) that (f i -ϕ i ) has to lie in the space spanned by the columns of Γ i , it turns out that

∆f i f i -ϕ i = Γ i θ i ∀i ∈ {1...m}, (9)
where θ i is the ((p -l i ) × 1) vector of free parameters. The expression of one column of the profile matrix may then be deduced

f i = ϕ i + Γ i θ i ∀i ∈ {1...m}. (10)
The expression in (9) reads in matrix form,

∆F = F -Φ, ( 11
)
where ∆F is the matrix gathering each column ∆f i , i = 1...m. It may be easily stated by using ( 6) and ( 7) that

∆F • Ω = 0. ( 12
)
From [START_REF] Fishler | Detection of signals by information theoretic criteria: general asymptotic performance analysis[END_REF], we can express F with respect to ∆F and Φ as,

F = ∆F + Φ. ( 13
)
By combining the relationship [START_REF] Gillis | Fast and robust recursive algorithms for separable nonnegative matrix factorization[END_REF] with the fact that 1 = Ω + (1 -Ω), we obtain

F = Ω • Φ + Ω • ∆F + (1 -Ω) • ∆F + (1 -Ω) • Φ. ( 14
)
Introducing the above properties ( 6) and ( 12) into ( 14) leads to the simplified expression of the profile matrix

F = Ω • Φ + (1 -Ω) • ∆F. (15)
Since (1 -Ω), Ω, Φ, ∆F are non-negative matrices, the profile matrix verifies F ≥ Φ, according to Relation [START_REF] Grippo | On the convergence of the block nonlinear gauss-seidel method under convex constraints[END_REF]. This matrix form is also equivalent to its vectorial expression [START_REF] Evans | Effect of sea salt calcium carbonate interactions with nitric acid on the direct dry deposition of nitrogen to Tampa Bay, Florida[END_REF]; both may be used in the following sections.

Constrained problem formulation

The aim of this section is to find new update rules which verify the defined above constraints. Constraints defined through [START_REF] Cichocki | NMFLAB for signal processing[END_REF] have to be incorporated and-similarly to (5)-multiplicative rules are investigated. The whole minimization problem may be formulated as

min G,F Tr (X -G • F ) T ((X -G • F ) • W ) s. t. G ≥ 0, F ≥ Φ, F = Ω • Φ + (1 -Ω) • ∆F. (16)
The proof of the KKT conditions are detailed in Appendix A for the case Φ = 0. The extension of these conditions should be easily derived when Φ is different from the null matrix. For the sake of clarity, let us define

Ω (1 -Ω) and F Ω • F + Φ.
The KKT conditions (see Appendix A for details) for the constrained weighted problem remain unchanged, i.e.,

F ≥ 0, G ≥ 0, Ω • G T ((G F -X) • W ) ≥ 0, (W • (G F -X)) F T ≥ 0, F • Ω • G T ((G F -X) • W ) = 0, G • (W • (G F -X)) F T = 0. ( 17
)
The profile matrix (respectively the contribution matrix) should be split into columns (respectively rows). The problem ( 16) is only solved for the profile matrix since the same technique may be applied to the contribution matrix. To solve [START_REF] Guillamet | Introducing a weighted non-negative matrix factorization for image classification[END_REF], the criterion may be split into independent parts which are only related to one column of the profile matrix. The problem may be formulated as a collection of quadratic optimization sub-problems,

J (f i ) = (x i -Gf i ) T D w i (x i -Gf i ) s. t. M i f i = δ i and f i ≥ ϕ i , ∀i ∈ {1...m}, (18)
where x i stands for the i th column of the X matrix. Remember that the matrix D w i is defined by D w i diag(w i ), where w i accounts for the i th column of the weight matrix W .

The sub-problem (18) shall be investigated for two different cases.

The case

δ i = 0
Here, we are interested in the solution of the special case when δ i = 0 which leads to ϕ i = 0. Substituting ϕ i = 0 into (10) yields a simplified expression, which tells that f i is a linear combination of column vectors from Γ i :

f i = Γ i θ i ∀i ∈ {1...m}, (19)
where θ i ((p -l i ) × 1 ) stands for the free parameters vector extracted from f i . By substituting [START_REF] Hopke | PM source apportionment and health effects: Intercomparison of source apportionment results[END_REF] into [START_REF] Ho | Non-negative matrix factorization with fixed row and column sums[END_REF], it turns out to be a weighted non-negative least square function, depending on the free parameters vector θ i :

K (θ i ) = (x i -(GΓ i )θ i ) T D w i (x i -(GΓ i )θ i ) s. t. θ i ≥ 0 ∀i ∈ {1 . . . m}. (20)
It appears to be a WNMF (as detailed in [START_REF] Ho | Non-negative matrix factorizations algorithms and applications[END_REF]) applied to the search of the vector θ i and to a new contribution matrix GΓ i . Ho [START_REF] Ho | Non-negative matrix factorizations algorithms and applications[END_REF] expressed the vector at iteration (k + 1) as a function of the vector at iteration k by using a maximization minimization strategy. As a consequence, we can respectively replace f i and G in his contribution by θ i and GΓ i . The solution may be obtained also by splitting the first part of (5) into vectorial form, and we obtain the following update rule for θ i :

θ k+1 i ← θ k i • Γ T i G T (D w i x i ) Γ T i G T (D w i GΓ i θ k i ) ∀i ∈ {1...m}. (21)
Referring to [START_REF] Hopke | PM source apportionment and health effects: Intercomparison of source apportionment results[END_REF], the expression of the i th column of the profile matrix may be deduced and reads

f k+1 i ← f k i • Γ i Γ T i G T (D w i x i ) Γ i Γ T i G T (D w i Gf k i ) ∀i ∈ {1...m}.
Given that Γ i Γ T i = diag(1 p×1 -ω i ), with ω i the i th column of the matrix Ω, it may be noticed that this operator only selects active components among the profile vector in [START_REF] Choi | Algorithms for orthogonal non-negative matrix factorization[END_REF]. The previous expression may be summarized into matrix formulation, i.e.,

F ← F • (1 p×m -Ω) • G T (W • X) G T (W • (GF )) . (22)
In ( 22), (1 p×m -Ω) has to be seen as a constraint mask applied to the components of the profile matrix. As a result, ( 22) is equivalent to WNMF since an appropriate choice of the initial profile (F 0 • Ω = F 0 ) leads exactly to the same iterations.

The case δ i = 0

Back to [START_REF] Evans | Effect of sea salt calcium carbonate interactions with nitric acid on the direct dry deposition of nitrogen to Tampa Bay, Florida[END_REF], let us use it in the case when, δ i is different from zero. By replacing it into (18), it appears as a non-negative least square cost function

K (θ i ) = (x i -Gϕ i -(GΓ i )θ i ) T D w i (x i -Gϕ i -(GΓ i )θ i ) s. t. θ i ≥ 0 ∀i ∈ {1...m}. (23)
It should be noticed that ( 20) and ( 23) are similar, except that x i in (20) has been replaced by (x -Gϕ) i in [START_REF] Laversin | Tracing the atmospheric particulate emissions from an iron and steel works[END_REF]. From [START_REF] Kim | Nonnegative matrix factorization based on alternating nonnegativity constrained least squares and active set method[END_REF], it results in a new expression of the free parameters vector at iteration (k + 1):

θ k+1 i ← θ k i • Γ T i G T D w i (x i -Gϕ i ) Γ T i G T (D w i GΓ i θ k i ) ∀i ∈ {1...m}.
Expression [START_REF] Fishler | Detection of signals by information theoretic criteria: general asymptotic performance analysis[END_REF] enables to update the i th column of the profile matrix:

f k+1 i -ϕ i ← f k i -ϕ i • Γ i Γ T i G T D w i (x i -Gϕ i ) Γ i Γ T i G T (D w i G(f k i -ϕ i ))
∀i ∈ {1...m}.

The previous equation may be written in matrix form, in a similar way as in [START_REF] Kim | Algorithms for nonnegative matrix and tensor factorizations: a unified view based on block coordinate descent framework[END_REF], by noticing that the matrix Γ i Γ T i acts as a component selector operator equivalent to the constraint mask (1 p×m -Ω):

F -Φ ← (F -Φ) • (1 p×m -Ω) • G T (W • (X -GΦ)) G T (W • (G(F -Φ))) . ( 24
)
This equation is the global expression which enables to update iterates of the profile matrix. It is the generalization of ( 22) and has to be compared to Expression [START_REF] Choi | Algorithms for orthogonal non-negative matrix factorization[END_REF], which stands for the unconstrained version. It can be seen as the constrained WNMF applied to (X -GΦ) instead of X (the non-negativity constraint is required on the equivalent data). Moreover, it can be easily checked that any limit point called F verifies the complementary slackness condition (last left condition from ( 17)). However, this method suffers from the same convergence problems as the NMF multiplicative updates proposed by Lee and Seung2 [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. Indeed, Gonzales and Zhang pointed the existence of numerical examples showing the fail of this algorithm to reach a stationary point [START_REF] Gonzales | Accelerating the Lee-Seung algorithm for nonnegative matrix factorization[END_REF]. As our method can be reduced to the same approach by setting W = 1 n×m and by imposing Ω = Φ = 0 p×m , the same comment here applies. Non stationarity may typically occur when some profile components are initialized to zero.

In the same way, the contribution matrix is updated according to unconstrained WNMF:

G ← G • (W • X)F T (W • (GF ))F T . ( 25
)
However, an unconstrained WNMF on the contribution matrix does not ensure the non-negativity of (X -GΦ). A projection step of equivalent data is thus performed before updating the contribution matrix.

Outline of the algorithm

From the previous results, we derive the algorithm whose steps are detailed in Algorithm 1. The normalization is made after each computation step by normalizing each row of the profile matrix and by updating the corresponding column of the contribution matrix. As the normalization has no effect on the product G • F , the value of the trace in ( 16) remains unchanged before or after the normalization. However, it makes the constraints unverified within iterations. From a practical point of view, if a limit point is Algorithm 1 Update rules of the CWNMF method.

Reading Φ and Ω.

Initialization of G and F Normalize F and G while the stopping rule is not checked do Check for the non-negativity of (X -GΦ) and project it if necessary Search for F at constant G according to [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] Search for G at constant F according to [START_REF] Lin | On the convergence of multiplicative update algorithms for non-negative matrix factorization[END_REF] Normalize F and G end while reached (i.e., F k+1 = F k ), constraints are observed to be verified asymptotically. Despite the fact that-as for other NMF methods [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF][START_REF] Ho | Non-negative matrix factorizations algorithms and applications[END_REF]-we cannot guarantee the stationarity when a limit point is reached, we never faced convergence issues-i.e., non-stationarity-in the numerous experiments we conducted on our dataset. It should be noticed that limit points tending to zero should be checked with great care. Indeed, although the third KKT condition is always satisfied, the sign of the gradient in the second KKT condition may be negative.

Lastly, we did not meet any kind of issue with the initialization of the matrices, as reported in Appendix B. It must be noticed that the equality constraints limit the number of free parameters to initialize, which probably tends to make the informed NMF method more robust to initialization than blind NMF, as shown in our experiments described in Appendix B.

Application to industrial source identification

Introduction

Data relative to airborne particles with size lower than 10 µm collected at Dunkerque were considered. This site is a medium city located in Northern France on the coast of the North Sea (Figure 1). There is an important industrial activity, especially with the presence of an integrated steelwork complex. Sampling was performed on the roof of the Research Institute and the sampling point is directly under the influence of emissions of the industrial sector when the wind blows from direction 260 degres to 320 degres. In this work, we have selected 92 samples collected under a wider wind sector (i.e. from 250 degres to 60 degres) to consider other influences such as marine, crustal and long-range transport origins. Concentrations of 19 elements and ions listed in Table 1 have been determined for each sample. In our previous work [START_REF] Delmaire | Une version pondérée de la factorisation matricielle non négative pour l'identification de sources de poussières. Application au littoral de la mer du nord[END_REF], we focused our attention on the identification of main sources contributing to atmospheric particulate background. Four particles types were evidenced: sea salts, aged sea-salts, secondary inorganic aerosols and crustal particles. Their respective profiles were in good agreement with results from literature [START_REF] Evans | Effect of sea salt calcium carbonate interactions with nitric acid on the direct dry deposition of nitrogen to Tampa Bay, Florida[END_REF][START_REF] Ten-Harkel | The effects of particle-size distribution and chloride depletion of sea-salt aerosols on estimating atmospheric deposition at a coastal site[END_REF][START_REF] Zhao | Acidic species and chloride depletion in coarse aerosol particles in the US east coast[END_REF]. We can assume that the composition of such emissions does not vary with time. In the present work, we aim to identify the industrial emission sources which additionally influence the level of atmospheric particulate at our sampling site. The integrated steelwork complex gathers processes devoted to coke making, iron ores sintering, steel production and slag treatments. A ferromanganese alloys production plant concerned by atmospheric emissions is also located next to the above industrial site. The particulate emissions have different chemical compositions from one process to another and do not always occur simultaneously. For these reasons, it is not adequate to use a single profile to track the impact of this industrial site on atmospheric particulate level. In this study, it is proposed that 5 industrial sources are expected to contribute significantly to airborne particulate matter at a local scale [START_REF] Laversin | Tracing the atmospheric particulate emissions from an iron and steel works[END_REF].

Major source identification

We assume that 4 background sources are completely known whereas 5 industrial sources are partially unknown. The matrices related to the constraints-i.e., Ω and Φ-are respectively given in Tables 1 and2. The first four rows represent background sources while the others account for industrial sources. When a value is set to 1, it means that the corresponding species is constrained. For example, copper for the sixth source is set to 0.

Al Ca Cr Cu Fe K Mg Mn Na Ni Pb Sn Ti V Zn Cl -NO - 3 SO 2- 4 NH + 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 0 1 1 1 1 0 1 1 0 1 0 0 0 1 0 1 0 0 1 0 1 0 0 0 1 1 1 0 1 0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 1 1 1 1 0 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 1 1 0 1 1 0 1 0 1 0 1 1 1 1
Table 2: Matrix Φ of set profile values.

Al Ca Cr Cu Fe K Mg Mn Na Ni Pb Sn Ti V Zn Cl -NO - 3 SO 2- 4 NH + 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 421 316 263 367 140 .2 .1 147 136 64 2.5 122 .1 .1 0 15 .3 .3 3 2.4 .5 .5 0 10 0 0 0 9 30 0 253 0 0 0 0 0 0 200 280 217 0 0 12 0 0 0 12 37 0 312 0 0 0 0 0 0 550 0 78 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

To sum up, 55 parameters are free in the profile matrix, they are located in the 5 last rows corresponding to the 5 industrial sources. No constraint is applied on the contribution matrix. As a result, it contains 92 × 9 free parameters. The estimation of these parameters is ordered according to sources in Figures 23456. In each figure, the top plot shows the contribution of one source (corresponding to one column of the contribution matrix) versus samples, while the middle one gives a rate of the whole mass emitted. The bottom graph shows the source profile (corresponding to one row of the profile matrix). The fifth source denoted blast furnace / steel plant sources is represented on Figure 2 while the steel slag source is given in Figure 3. Two sintering plant sources are detailed on Figures 4 and5 while a ferromanganese plant source is reported in Figure 6.

By comparing these profiles with experimental data defined for each source [START_REF] Laversin | Tracing the atmospheric particulate emissions from an iron and steel works[END_REF], the profile issued from Figure 2 is referred to as iron rich particles emitted by blast furnaces or steel plant. These two sources have got very similar chemical identities. These features have also been reported for similar processes located on other sites all over the world [START_REF] Machemer | Characterization of airborne and bulk particulate from iron and steel manufacturing facilities[END_REF][START_REF] Oravisjarvi | Source contributions to PM2.5 particles in the urban air of a town situated close to a steel works[END_REF]. The second profile (see Figure 3) is characterized by a high proportion of calcium and iron and then various metallic impurities. Unambiguously, this profile can be attributed to fine particles which gather The third profile (see Figure 4) highlights species in the following order: Fe> Ca> Al> Mg> Mn> K> Na. It is ascribed to ores sintering plant source and such particles correspond to fugitive emissions depending on all handling performed on sintered ores. It comes from the heating of raw materials. However, a lack of potassium and sodium may be noted. The fourth source profile (Figure 5) is expected to be mainly composed with Cl -, K, and Fe. It stands for the sintering chimney point source. But it can be noticed that there is a lack of Fe particles and also NH + 4 is over estimated. Finally, a Ferromanganese plant source may be recognized on Figure 6, it can be noted that the expected profile is very well recovered.

A brief summary of significant sources show that 3 sources (blast furnaces source, steel slags source, sintering chimney point source) present a few peaks over 3 µg/m 3 . The steel slag source and the ferromanganese plant source are minor contributors to particles pollution. Hence, these results enable to evidence the impact of point source emissions on atmospheric particulate levels and particularly provide for the first time a relative contribution of fugitive emissions from this site.

Conclusion

We have introduced in this article the use of basic equality constraints and have derived theoretical expressions of constrained WNMF. We provided global multiplicative rules which take into account at the same time linear equality constraints and weights. The second interest stands in the choice of the application which aims to find natural and industrial profiles. Partial knowledge of these profiles enables to focus the search on species which are completely unknown. As a consequence, according to chemists, our results allow them to better understand the industrial profiles. This shows the relevance of our proposed approach. Upcoming prospects concern the integration of inequality constraints in the algorithms.

The differentiation of H 2 and H 3 then results in

∂H 2 ∂F = Ω • (G T (X • W )), ∂H 3 ∂F = 2[G T (W • (G(Ω • F )))] • Ω. (A.1)
According to (A.1), the differentiation of the function H(.) gives

∂H ∂ F = 2Ω • G T ((G[Ω • F ] -X) • W ) . (A.2)
The Lagrangian function may be expressed as

L(G, F ) = 1 2 H( F ) -Tr(Γ G • G) -Tr(Γ F • F • Ω).
Cancelling the derivatives of this Lagrangian function with respect to F yields the following relation,

∂L ∂ F = Ω • G T ((G[Ω • F ] -X) • W ) -Γ F = 0, (A.3)
leading to an expression of Γ F . Karush-Kuhn Tucker conditions for the minimization of ( 16) may be summed up as

F ≥ 0, G ≥ 0, ∂H ∂ F ≥ 0, ∂H ∂G ≥ 0, F • Γ F = 0, G • Γ G = 0. (A.4)
For the sake of clarity, let us define F Ω • F . By replacing the relations (A.2) and (A.3) into (A.4) and by using the same analogy with the matrix G, the final KKT conditions for the constrained weighted problem are measures-similar to SNRs-on the contribution and the profile matrices.

F ≥ 0, G ≥ 0, Ω • G T ((G F -X) • W ) ≥ 0, (W • (G F -X)) F T ≥ 0, F • Ω • G T ((G F -X) • W ) = 0, G • (W • (G F -X)) F T = 0.
• When applied to the columns of G, the chosen measure is named Mixing-Error Ratio (MER) [START_REF] Vincent | The 2008 signal separation evaluation campaign: A community-based approach to large-scale evaluation[END_REF]. We express each estimated column ĝj of the matrix G as ĝj = g coll j + g orth j ,

where g coll j and g orth j are respectively collinear and orthogonal to the true vector g j . The MER (in dB) associated with this estimated column is then expressed as MER j = 10 log 10 g coll j 2 / g orth j 2 .

This definition provides one MER for each source-respectively denoted as S 1 , S 2 , and S 3 -so that p MERs are obtained.

• When applied to the lines of F -or equivalently to the columns of F Twe measure a Signal-to-Interference Ratio (SIR) which is classically met in blind source separation [START_REF] Vincent | The 2008 signal separation evaluation campaign: A community-based approach to large-scale evaluation[END_REF]. The SIR is computed for each estimated source. From a practical point of view, we apply the above MER algorithm on the columns of F T as an SIR measure.

The sensitivity of the algorithm is studied by choosing different initializations of the profile matrix F init . We exclusively use initial profile matrices which satisfy the constraints defined in Tables B. [START_REF] Bertsekas | Nonlinear Programming[END_REF] and B.4. For each initialization, we run the algorithm during 20 thousand iterations.

We performed numerous tests in this configuration and all the obtained enhancement provided by our NMF approach was consistent. For the sake of readibility, we only gather 12 significant initial profile matrices in Table B. [START_REF] Cichocki | Generalized alpha-beta divergences and their application to robust non-negative matrix factorization[END_REF]. This table provides three performance indices:

1. the first three columns show the SIR obtained for each source when each line of F is set to its initialization. We see that the performance is pretty low, as it ranges from 1.15 to 24.07 dB. While the SIR can be really low for Source 2 and Source 3, we could not reach an SIR lower than 10 dB with Source 1 because the mass to be distributed over the free parameters is quite low, as discussed above. 2. The last three columns show the SIRs obtained after running the algorithm. We notice that, whatever the initialization, the SIR for each source keeps almost constant and is high, i.e., between 54.61 and 81.66 dB. This shows that the algorithm converges to the same stationary point. 3. Such a result is seen again in the three middle columns. They show the MERs obtained after applying the NMF algorithm. Here again-and as for the SIR-the measure obtained for each source remains almost constant over the trials. It can be noticed that SIRs-estimated over the profile matrix-are greater than MERs which estimated over the contribution matrix. It may be due to the fact that the constraints are applied on the profile matrix only.

Other experiments have been conducted with less constraints. We noticed that relaxing some of the above constraints made the algorithm more sensitive to the initialization. This was expected as less constraints means more parameters to estimate.

To conclude, these experiments show the relevance of our informed NMF method. It not only provides a good performance using the additional knowledge on the profile matrix F , but it is also less sensitive to the initialization than its blind counterpart.

Figure 1 :

 1 Figure 1: Location of the area under study.

Figure 2 :

 2 Figure 2: Blast furnaces source.

Figure 3 :

 3 Figure 3: Steel slag source.

Figure 4 :

 4 Figure 4: Ores sintering plant non point source.

Figure 5 :

 5 Figure 5: Sintering chimney point source.

Figure 6 :

 6 Figure 6: Ferromanganese plant source.

Table 1 :

 1 Contraint matrix Ω.

Table B .

 B 5: Real profile matrix F .

	Fe	Ca	SO 4	Zn	Mg	Al	Cr
	700 100	80	60	40	20	0
	300 400	5	0	200	75	20
	400 200	80	0	120 200	0

  Table B.6: Performance of the NMF approach for different initializations. Performance criteria: SIR (in dB) and MER (in dB).

	trial	SIR (on F init )		MER (on G)			SIR (on F )
	number	S 1	S 2	S 3	S 1	S 2	S 3	S 1	S 2	S 3
	1	24.07 18.56 16.50 29.89 30.06 32.44 54.61 78.18 54.81
	2	17.31 18.56 16.50 29.89 30.05 32.44 54.78 81.66 54.61
	3	13.96 18.56 16.50 29.89 30.06 32.44 54.64 80.19 54.78
	4	24.07	3.87	16.50 29.89 30.06 32.44 54.61 78.18 54.81
	5	24.07	6.11	16.50 29.89 30.06 32.44 54.61 78.18 54.81
	6	24.07 18.56	3.19	29.89 30.06 32.44 54.61 78.18 54.81
	7	24.07 18.56	1.15	29.89 30.06 32.44 54.61 78.18 54.81
	8	17.31 18.56	1.15	29.89 30.06 32.44 54.61 78.18 54.81
	9	13.96 18.56	1.15	29.89 30.06 32.44 54.61 78.18 54.81
	10 13.96	6.11	1.15	29.89 30.06 32.44 54.61 78.18 54.81
	11 13.96	3.87	1.15	29.89 30.06 32.44 54.61 78.18 54.81
	12 10.47	3.87	1.15	29.89 30.06 32.44 54.61 78.18 54.81

Please note that other cost functions are also widely used, e.g., Kullback divergence[START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] or flexible parametric divergences such as generalized alpha-beta divergences[START_REF] Cichocki | Generalized alpha-beta divergences and their application to robust non-negative matrix factorization[END_REF].

Please note the existence of alternative algorithms which fix this issue[START_REF] Gillis | Accelerated multiplicative updates and hierarchical als algorithms for nonnegative matrix factorization[END_REF].

Acknowledgements

The authors would like to thank the anonymous reviewers for their helpful comments and valuable suggestions to improve the quality of the paper. This research was financed by the integrated steel and mining company Arcelor-Mittal.

Appendix A. Expression of the KKT conditions in the constrained problem

Karush-Kuhn Tucker conditions are carried out without loss of generality in the case when Φ = 0. It may be easily extended to the general case. For the sake of concision, let assume that Ω (1 p×m -Ω). The minimization of the criterion [START_REF] Guillamet | Introducing a weighted non-negative matrix factorization for image classification[END_REF] leads to consider a function H(.), defined by

Appendix B. Sensitivity of the algorithm to the initialization

The proposed NMF algorithm has been tested in the framework of pollutant source separation. As explained in Subsection 3.3-and as it is the case for classical NMF approaches [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF][START_REF] Ho | Non-negative matrix factorizations algorithms and applications[END_REF]-the robustness of our method to initialization is not theoretically guaranteed (but was shown not to be an issue with the considered application). In this appendix, we aim to investigate this behaviour: we generate synthetic data according to industrial source profiles available in the literature [START_REF] Oravisjarvi | Source contributions to PM2.5 particles in the urban air of a town situated close to a steel works[END_REF]. Seven chemical species-listed in Table B.3-are under consideration. Three industrial sources coming from steel industry (blast furnaces source, steel slag source, ores sintering plant) shown in this order in Table B.5 are assumed to be active.

The data matrix consists of 50 samples and 7 species, with a known uncertainty measure σ ij -provided by a chemical expert-associated to each data point x ij . The uncertainty measure allows us to compute the weight matrix W , according to w ij = σ -2 ij . A uniform noise ranging in [-min{λ σ ij , x ij }; λ σ ij ] is added while keeping non-negativity of the data. Please note that λ is related to an input Signal-to-Noise Ratio (SNR). In the experiments below, the input SNR is equal to 30 dB. The actual profile matrix is given in Table B. [START_REF] Choi | Algorithms for orthogonal non-negative matrix factorization[END_REF].

In our method, we can select which components of the profile matrix are to be set, making it highly flexible. In the tests reported in this appendix, 9 constraints are taken into account and are shown in Table B.3, where the 1s indicate the location of the constraints in the profile matrix. The matrix Φ is provided in Table B.4. Note that 70 percent of the first source profile is set, so that it only subsists 30 percent of it to be distributed over the remaining species. To assess the quality of the fit, we apply different performance


Alice Allen

Cecilia Aragon

Christophe Becker

Jeffrey C Carver

Andrei Chis

Benoit Combemale

Mike Croucher

Kevin Crowston

Daniel Garijo

Ashish Gehani

Carole Goble

Robert Haines

Robert Hirschfeld

James Howison

Kathryn Huff

Caroline Jay

Daniel S Katz

Claude Kirchner

Kateryna Kuksenok

Ralf Lämmel

Oscar Nierstrasz

Matthew Turk

Rob V Van Nieuwpoort

Matthew Vaughn

Jurgen Vinju

Lightning Talk: "I solemnly pledge" A Manifesto for Personal Responsibility in the Engineering of Academic Software

Keywords:

Software is fundamental to academic research work, both as part of the method and as the result of research. In June 2016 25 people gathered at Schloss Dagstuhl for a week-long Perspectives Workshop and began to develop a manifesto which places emphasis on the scholarly value of academic software and on personal responsibility. Twenty pledges cover the recognition of academic software, the academic software process and the intellectual content of academic software. This is still work in progress. Through this lightning talk, we aim to get feedback and hone these further, as well as to inspire the WSSSPE audience to think about actions they can take themselves rather than actions they want others to take. We aim to publish a more fully developed Dagstuhl Manifesto by December 2016.

I. INTRODUCTION

Software is fundamental to academic research work, both as part of the method and as the result of research. With the advent of artifact evaluation committees of conferences, journals that include source code and running systems as part of the published artifacts, as well as the increasing push to reproducibility, we foresee that software will only increase in importance as part of the academic process.

In June 2016, 25 people gathered at Schloss Dagstuhl for a weeklong Perspectives Workshop [START_REF]Dagstuhl Perspectives Workshop 16252[END_REF] to produce a roadmap towards future professional software engineering for softwarebased research instruments and other software produced and This work is licensed under a CC-BY-4.0 license. used in an academic context (i.e., for research, not administration). The group was carefully picked to be broad in its range of disciplines (including Astronomy, Social Sciences, Biology, Chemistry, Computer Science, Physics, and Humanities), roles (including computer science researchers, general and specialist research software engineers and systems administrators) and career stages (from PIs and institute heads to PhD students and postdocs). Despite its ambiguous title, "Engineering Academic Software," the workshop was on the Engineering of Academic Software (not software for engineers).

A Dagstuhl Perspectives Workshop results in a Manifesto. The open science and research software communities have been very active in creating manifestos: the Science Code Manifesto [2], Karlskrona Manifesto for Sustainability Design [START_REF]Karlskrona Manifesto for Sustainability Design[END_REF], the Reproducibility Manifesto [START_REF]Reproducibility manifesto[END_REF], and Principles for Software Citation [START_REF] Smith | FORCE11 Software Citation Working Group[END_REF], FAIR data [START_REF] Wilkinson | The FAIR guiding principles for scientific data management and stewardship[END_REF], and so on. Why do we need another Manifesto?

First, our manifesto is to be less about what others should do, and more about what we, as individuals, should do. That is, it is more in the style of a personal responsibility pledge like those on open access [7] and open peer review [8]. It is easy for us to declare "the community should do X", "funding panels should do Y" and "promotion committees should do Z", while conveniently forgetting that we are the community, we are the panelists, we are the committee members.

Second, our manifesto places emphasis on the scholarly value of academic software. For some of our group, engineering academic software is chiefly a means to an endto pro-duce robust and reliable software as an instrument as part of a wider research investigation. For others the software is the end in itselfthe software is the research. In both cases the software has scholarly merit.

II. THE DRAFT PERSONAL MANIFESTO

Currently at 20 there are too many pledges so we are actively working on reducing the number and simplifying the message. At the WSSSPE4 meeting, through this lightning talk, we aim to get feedback and contribute to this process.

Table I presents our pledge list organized into three broad areas: (1) recognition of academic software, (2) academic software development processes and (3) the intellectual content of academic software. Each pledge should be actionable by an individual. Each pledge has a story that will be developed in the full manifesto.

TABLE I. MANIFESTO DRAFT PLEDGES

A. Recognition of academic software 1 I will properly cite software used to produce my research results 2 I will point out improper or missing citations to software when I am reviewing publications. 3 I will make explicit how to cite the software I make available. 4 I will recommend software experts for funding agencies to include in their review processes. 5 I will invite developers of software that enables my research to be coauthors on my papers. 6 I will recognize software contributions in hiring and promotion within my institution. 7 I will recognize software contributions at conferences, e.g. dedicated sessions, and prizes. 8 I will support and publish in journals that recognise software contributions. 9

I will contribute to sustaining the software I rely on for my research.

B. Academic software development processes

10 I will develop software as open source right from the start whenever possible. 11 I will document my academic software for users with instructions and examples. 12 I will package, release and archive versions of my software 13 I will consider and document the sustainability of my research software. 14 I will publish how I organize and run my software projects 15 I will match software engineering practices I recommend to the needs and resources of projects. 16 I will help scientists improve the quality of their software without passing judgment.

C. The intellectual content of academic software

17 I will acknowledge that source code is a legitimate part of the academic discourse. 18 I will publish the intellectual contributions of my research software. 19 I will distinguish the intellectual contribution of my software from its service contribution. 20 I will examine the source code of academic software contributions and encourage others to do so as well.

The pledges A and B (1-16) are targeted at both developers and users of academic software; pledges C (17-20) are more focused on developers. Pledges to recognize software are selfevidently responsibilities to be borne by users. Perhaps more implicit are the benefits to users embedded in the software development processes. Open source (10) and documentation (11) enables feedback and potential community engagement during critical design stages, and examples and instructions by their definition must embed an understanding of use cases and usability pitfall. Improving transparency of process, and analytic consideration of future use cases (12, 13, 14) are valid usercentered approaches and (15, 16) directly speak to assessing user needs and crafting appropriate interventions.

III. NEXT STEPS

This is still work in progress. We recognize that 20 pledges is roughly twice as many as desirable, and that the pledges need to be succinct and easy to understand and adopt, more in the style of the Reproducibility Manifesto [START_REF]Reproducibility manifesto[END_REF]. We are currently (1) looking to merge pledges where feasible and (2) exploring the use of roles to structure them; for example "When I develop software ...When I write research papers ...When I evaluate colleagues ..." and so on. At WSSSPE we intend to run an online vote on the pledges structured on roles.

We aim to develop the full manifesto by December 2016 and publish it as a Dagstuhl Manifesto [9]. We intend that community groups promote its contents among researchers and research software engineers, and we use it to influence decision makers to enable our respective communities to execute these pledges with moral and financial support.

ACKNOWLEDGMENT

We thank Schloss Dagstuhl -Leibniz-Zentrum für Informatik GmbH for supporting and hosting Workshop 16252.


X Garbet

J H Ahn

S Breton

P Donnel

D Esteve

R Guirlet

H Lütjens

T Nicolas

Y Sarazin

C Bourdelle

O Février

G Dif-Pradalier

P Ghendrih

V Grandgirard

G Latu

J F Luciani

P Maget

A Marx

A Smolyakov

This paper investigates the interplay of neoclassical, turbulent and MHD processes, which are simultaneously at play when contributing to impurity transport. It is shown that these contributions are not additive, as assumed sometimes. The interaction between turbulence and neoclassical effects leads to less effective thermal screening, i.e. lowers the outward flux due to temperature gradient. This behavior is attributed to poloidal asymmetries of the flow driven by turbulence. Moreover sawtooth crashes play an important role to determine fluxes across the q = 1 surface. It is found that the density profile of a heavy impurity differs significantly in sawtoothing plasmas from the one predicted by neoclassical theory when neglecting MHD events. Sawtooth crashes impede impurity accumulation, but also weaken the impurity outflux due to the temperature gradient when the latter is dominant.

Synergetic effects of collisions, turbulence and sawtooth crashes on impurity transport.

Introduction

Impurity transport is a subject of renewed interest due to several facts. First plasma facing components made of tungsten are currently tested in several tokamaks, in view of a future implementation in Iter. This has triggered an active research on means to prevent influxes of heavy impurities. Furthermore, burning plasmas will produce helium, which should be expelled from the plasma core and pumped. Finally impurities with medium charge number are commonly injected in the scrape-off layer to cool down the edge plasma. These impurities should be prevented from reaching the plasma core to avoid dilution and excessive radiative losses. These constraints raise questions on the mechanisms that underlie impurity transport for various mass and charge numbers. It turns out that collisional, turbulent and MHD processes are simultaneously at play and contribute to impurity transport. It is shown here that these contributions are not additive, as assumed sometimes. Moreover the interplay between turbulence, neoclassical and sawtooth crash effects leads to less effective thermal screening, i.e. a weakening of the outward flux due to temperature gradient. It is reminded that the neoclassical impurity flux is given by the following expression

Γ Zψ = -D nc Z ∂N Z ∂r + V nc Z N Z
where N Z is the impurity density. The ratio of the pinch velocity

V nc Z to the diffusion coef- ficient D nc Z is of the form V nc Z D nc Z = Z ∂ ln N i ∂r + HZ ∂ ln T i ∂r ,
where N i and T i are the ion density and temperature. The number H is the thermal screening factor and is equal to -1 2 when the impurity is in the Pfirsch-Schlüter regime, while the main ion species is in the bananaplateau regime. A multi-species collision operator has been implemented in the GYSELA gyrokinetic code, in order to compute neoclassical and turbulent particle fluxes on an equal footing [1]. It has been benchmarked against the theory of neoclassical transport [2] and the NEO code [3,4]. The agreement is found to be satisfactory for both the diffusion coefficient and the pinch velocity for various collisionality regimes (see Fig. 1 and Fig. 2).

The upgraded version of the GYSELA code has been run above the instability threshold of Ion Temperature Gradient driven modes (ITG) [5]. It is found that neoclassical and turbulent processes contribute to impurity transport in a complex manner. This question has been assessed as follows. Three simulations are run for each impurity :

1. Purely neoclassic (called "neoclassic" hereafter): in this case, all non-axisymmetric toroidal modes are filtered out at each time step (i.e. all Fourier modes with n = 0 are set to zero, with n the toroidal mode number). The neoclassical flux is then the sum of both v E and magnetic drift v D contributions:

Γ neo z ψ = d 3 v F z (v D,z + v n=0 E,z ) • ∇r) ψ (1)
2. Mainly turbulent (called "turbulent" hereafter): in this case, single-species collisions only are retained (namely ν ii and ν zz ) , so that momentum or energy exchange be-tween species are not taken into account. Retaining intra-species collisions is important and sufficient to account, among others, for the collisional damping of zonal flows, which are known to contribute efficiently to turbulence self-regulation. In the case of trace impurities, the collision operator for the main ions only is important. Indeed, impurity-impurity collisions are negligible in this case. The turbulent flux is then governed by the electric drift:

Γ turb z ψ = d 3 v F z (v E,z • ∇r) ψ (2)
3. Full (called "total" hereafter): in this case, no simplification is made to the collision operator, nor any filtering applied to the electric potential. More precisely, all terms of the collision operator are retained, involving both intra-and inter-species collisions, in the turbulent regime. As an example, the case of tungsten is considered hereafter. The same trend is observed for other impurities. In Fig. 3 the total tungsten flux (red) from the self-consistent simulation is compared to the sum (black) of the neoclassical (blue) and turbulent (green) fluxes, that come from the reduced simulations. It appears that the total flux differs from the sum of the neoclassical and turbulent fluxes, by more than a factor two at some radial locations. This result shows that neoclassical and turbulent impurity transports are not additive. One explanation for this synergy comes from the impact of turbulence on the magnitude and radial shape of in/out poloidal asymmetries. The pattern of the impurity density in the poloidal plane is a signature of this multiscale dynamics. The left panel of Fig. (4) shows the level of n = 0 fluctuations of impurity density, while the right panel shows non axisymmetric n = 0 small scale fluctuations. A similar structure is observed for the electric potential. Small scale fluctuations are clearly modulated by n = 0 perturbations. Conversely, it is expected that turbulence generates n = 0, m = ±1 perturbations (m is the poloidal wave number), that produces a significant contribution to neoclassical transport. The interplay via large scale flows is expected to be a key mechanism to explain the synergistic effects that are observed. Moreover the coefficient of thermal screening, which quantifies the outward pinch due to the temperature gradient, is above the neoclassical value H = -1 2 in presence of turbulence, and can become positive. This is an unfavorable configuration, as it may lead to impurity accumulation. This behavior is attributed to the same mechanism, i.e. poloidal asymmetries of the flows driven by turbulence, via a mechanism similar to the effect of poloidal asymmetries (e.g. centrifugal force or potential asymmetries produced by RH heating). Such asymmetries are able to greatly modify neoclassical coefficients, up to one order of magnitude, as predicted theoretically [6,7]. The interplay of neoclassical transport and sawtooth events has been assessed by implementing a set of fluid equations in the XTOR-2F code [8], which model the impurity flows and are consistent with the gyrokinetic model in the collisional regime. The impurity density and parallel momentum equation read

∂ t N z + ∇ • (N z V z ) + ∇ × B B 2 . ∇P z e z = S z (3)
N z m z ∂ t V ,z + (V z • ∇)V ,z = -∇ P z + N z e z E + R ,zi (4)
where N z and V z are the impurity density and fluid velocity, e z and m z the charge and mass, S z is the impurity density source, and R ,zi is the parallel friction force given by

R ,zi = -m z N z ν zi (V ,z -V ,i ) - 2 5 C 0 m z N z ν zi Q ,i P i ( 5
)
where V ,i and Q ,i are the main ion parallel velocity and heat flux. Turbulent diffusion and pinch velocities can be added via an ad-hoc diffusion/convection transport model in Eq.( 3). Banana-plateau neoclassical transport coefficients can be added in the same manner. However the Pfirsch-Schlüter neoclassical flux should not be added in that way since it is already accounted for in the fluid equations Eqs. (3,4) . The expression of the collisional drag force Eq.( 5) is essential in that regard. Using the Braginskii value of 0 = 1 would be a mistake, since no thermal screening is found when all species (main ion and impurity) are in the Pfirsch-Schlüter regime [10]. The appropriate choice of C 0 to recover the correct value

H = -1 2 is C 0 = 3 2
, which is the one that has been implemented in the XTOR code. It has been verified that the profiles of impurity density agree with the neoclassical predictions in MHD quiet plasmas.

Once the Eqs. (3,4) have been implemented and verified, the XTOR-2F code has been run in a regime of parameters where a steady sawtoothing regime sets on [11]. It is found that the impurity density behaves differently with or without sawteeth [9]. Typically the impurity density profile is less peaked than the one computed without sawtooth collapses. This is due to the flattening effect of sawtooth crashes. Neoclassical relaxation processes are negligible during a crash, since they are too slow. Conversely neoclassical transport is the prominent transport channel during the recovery phase. Hence the interplay comes from modifications of the profiles during the crash, which persist during the recovery phase, thus modifying neoclassical fluxes. When the main ion density gradient is large chosen such that the impurity neoclassical flux is outward due to thermal screening. It appears clearly that the impurity profile after 9 crashes is less hollow than the one obtained without sawteeth, i.e. with neoclassical transport alone.

|∇N i /N i | ≥ |∇T i /2T i |,
For practical purposes, it is interesting to compare the post-crash profiles with those predicted by the Kadomtsev reconnection model [14,15,16]. The final impurity profile differs from the one predicted by the Kadomtsev reconnection model when the pre-crash impurity density profile is slightly depleted. This is due to sawteeth induced corrugations of the density, which are clearly seen in XTOR-2F simulations after each crash (see Fig. 6). Density perturbations also appear on the main ion post-crash density, and have been observed by fast-sweeping reflectometry tomography on Tore Supra [12]. On the contrary, the post-crash impurity profile is consistent with the Kadomtsev prediction when the pre-crash profile is hollow (Fig. 7), as already found in [13]. The agreement is not perfect though.

During the crash, the impurity transport dynamics is quite complex. In particular, it appears that impurities are concentrated after the crash within an annulus in between the inversion and mixing radii, while fine structures develop inside the inversion radius [13,9]. The flux can be decomposed in two components: one due to the perpendicular fluid E × B drift velocity (the diamagnetic component is negligible for high Z impurities), and the other one due to magnetic flutter. These fluxes are defined as follows

Γ z,δB = N z V ,z δB B (6)
and

Γ z,E = N z V E (7)
It is found that the radial flux due to E × B drift is much larger than the contribution of magnetic flutter in most places (see Fig. 8 and Fig. 9). Hence fine structures are due to convection, not to magnetic flutter.

Conclusion

In summary, it appears that collisions, turbulence and MHD processes contribute synergistically to impurity transport, thus questioning the usual assumption of additivity. The interplay between turbulent and neoclassical transport has been studied with the GYSELA gyrokinetic code, which has been equipped with a comprehensive collisional operator. It is found that the impurity flux found in a self-consistent simulation differs from the sum of neoclassical and turbulent fluxes calculated separately. This behavior is attributed to poloidal convective cells generated by turbulence, which modify neoclassical fluxes. Hence synergy is the result of poloidal asymmetries and interaction between disparate scales. Moreover the thermal screening effect gets weaker. The two fluid non linear MHD code XTOR-2F has been modified to model the dynamics of heavy impurities. The collisional drag force has been adjusted such that the Pfirsch-Schlüter impurity flux agrees with the one predicted by neoclassical theory, including the thermal screening term. Simulations of sawtoothing plasmas indicate that the time average impurity profile is flatter than the one predicted by neoclassical theory alone. This is equivalent to a weakening of the thermal screening term in situations where the outward pinch velocity due to temperature gradient dominates over the contribution of the density gradient responsible for accumulation. The post-crash profiles are reasonably in agreement with those predicted by the Kadomtsev reconnection model, when pre-crash profiles are sufficiently contrasted. Fine structures of the impurity density develop during the crash and persists some time. These structures results from small scale convective scales associated with the crash dynamics. The overall impurity profile relaxation is controlled by the E × B convection.

Figure 1 :

 1 Figure 1: Neoclassical diffusion coefficients vs ν * Z calculated with the GYSELA code, NEO code and analytical theory

Figure 2 :

 2 Figure 2: Neoclassical pinch velocities vs ν * Z calculated with the GYSELA code, NEO code and analytical theory

Figure 3 :

 3 Figure 3: Neoclassical (blue) and turbulent (green) fluxes calculated separately, and their sum (black) compared to the self-consistent total flux (red)

Figure 4 :

 4 Figure 4: Poloidal map of tungsten density fluctuations in a GYSELA gyrokinetic simulation. left panel: n = 0 axisymmetric fluctuations. right panel: n = 0 non-axisymmetric fluctuations

Figure 5 :

 5 Figure 5: Impurity profile with (blue line) or without sawtoothing activity (red line). The impurity neoclassical flux is outward due to thermal screening.

1 Figure 6 :Figure 7 :

 167 Figure 6: Comparison of the final impurity profile to Kadomtsev model for an initial slightly hollow profileFigure 7: Comparison of the final impurity profile to Kadomtsev model for an initial hollow profile

Figure 8 :

 8 Figure 8: Radial flux due to E × B drift velocity during a sawtooth crash

Figure 9 :

 9 Figure 9: Radial flux due to magnetic flutter during a sawtooth crash

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium and the French Research Federation for Fusion Studies and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. We benefited from HPC resources from GENCI, and IFERC.