Benoît Boyer

Axel Legay

Louis-Marie Traonouez

A Formalism for Stochastic Adaptive Systems

published or not. The documents may come

L'archive ouverte

A Formalism for Stochastic Adaptive Systems

Benoît Boyer, Axel Legay, and Louis-Marie Traonouez Inria / IRISA, Campus de Beaulieu, 35042 Rennes CEDEX, France Abstract. Complex systems such as systems of systems result from the combination of several components that are organized in a hierarchical manner. One of the main characteristics of those systems is their ability to adapt to new situations by modifying their architecture. Those systems have recently been the subject of a series of works in the software engineering community. Most of those works do not consider quantitative features. The objective of this paper is to propose a modeling language for adaptive systems whose behaviors depend on stochastic features. Our language relies on an extension of stochastic transition systems equipped with [START_REF] Corp | Advanced laboratory on embedded systems[END_REF] an adaptive operator that allows to reason about the probability that a system has to adapt its architecture over time, and (2) dynamic interactions between processes. As a second contribution, we propose a contract-based extension of probabilistic linear temporal logic suited to reason about assumptions and guarantees of such systems. Our work has been implemented in the Plasma-Lab tool developed at Inria. This tool allows us to define stochastic adaptive systems with an extension of the Prism language, and properties with patterns. In addition, Plasma-Lab offers a simulation-based model checking procedure to reason about finite executions of the system. First experiments on a large case study coming from an industrial driven European project give encouraging results.

Context

Critical systems increasingly rely on dynamically adaptive programs to respond to changes in their physical environments. Reasoning about such systems require to design new verification techniques and formalisms that take this model of reactivity into account [START_REF] Cheng | Software engineering for self-adaptive systems: A research roadmap[END_REF].

This paper proposes a complete formalism for the rigorous design of stochastic adaptive systems (SAS), whose components' behaviors and environment changes are represented via stochastic information. Adding some stochastic feature to components' models is more realistic, especially regarding the environment aspect, e.g. the probability of hardware failure, the fire frequency in a forest or a growing city population. . . We view the evolution of our system as a sequence of views, each of them representing a topology of the system at a given moment of time. In our setting, views are represented by a combination of Markov chains, and stochastic adaptive transitions that describe the environment evolution as transitions between different views of the SAS (e.g. adding or removing components). Each view thus associates the new environment behaviour and a new system configuration ( Global Verif.

Fig. 1. Illustration of SAS Methodology topology, addition or suppression of system components. . . ). The incremental design is naturally offered to the system architect who can extend easily an existing model by creating new views.

Properties of views can be specified with Bounded Linear Temporal Logic (BLTL) that allows to reason about finite execution. To reason about sequences of view, we propose Adaptive BLTL (A-BLTL) that is an extension of BLTL with an adaptive operator to reason about the dynamic change of views. We also show that the formalism extend to contracts that permit to reason about both assumptions and guarantees of the system.

Consider the system described in Figure 1. This system is composed by three different views linked by adaptive transitions (represented by dashed arrows). Each view contains some system components in a particular topology denoting a configuration of the SAS. Each local property p 1 , p 2 , p 3 is attached to one or more views. There is also global property Φ. The SAS is initially designed by View 1 and View 2 and the black dashed arrows. Properties p 1 , p 2 are validated for the corresponding views and Φ is validated against this complete initial version of the system. To fit with the upcoming settings of the system, the system architect updates the model by adding View 3 with new adaptive transitions (in grey). This requires to only validate p 1 and p 3 against View 3 and to validate again the global property Φ against the system including all the three views.

We propose a new Statistical Model Checking (SMC) [START_REF] Younes | Verification and planning for stochastic processes with asynchronous events[END_REF][START_REF] Sen | On statistical model checking of stochastic systems[END_REF] algorithm to compute the probability for a SAS to satisfy an A-BLTL property. SMC can be seen as a trade-off between testing and formal verification. The core idea of SMC is to generate a number of simulations of the system and verify whether they satisfy a given property expressed in temporal logics, which can be done by using runtime verification approaches [START_REF] Havelund | [END_REF]. The results are then used together with algorithms from the statistical area in order to decide whether the system satisfies the property with some probability. One of the key points to implement an SMC algorithm is the ability to bound a priori the size of the simulation, which is an important issue. Indeed, although the SAS can only spend a finite amount of time in a given view, the time bound is usually unknown and simulation cannot be bounded. To overcome the problem, we expand on the work of Clarke [START_REF] Younes | Statistical verification of probabilistic properties with unbounded until[END_REF] and consider a combination of SMC and model checking algorithm for untimed systems.

As a second contribution, we propose high-level formalisms to represent both SAS and A-BLTL/contracts. The formalism used to specify SAS relies on an extension of the Reactive Module Language (RML) used by the popular Prism toolset [START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF]. Properties are represented with an extension of the Goal and Contract Specification Language (GCSL) [START_REF] Arnold | Contracts and behavioral patterns for systems of systems: The EU IP DANSE approach[END_REF] defined in the DANSE IP project [START_REF] Danse | Designing for adaptability and evolution in sos engineering[END_REF]. This language offers English-based pattern to reason about timed properties without having to learn complex mathematics inherent to any logic.

Finally, as a last contribution, we have implemented our work in Plasma-Lab [START_REF] Boyer | Plasma-lab: A flexible, distributable statistical model checking library[END_REF] -a new powerful SMC model checker. The implementation has been tested on a realistic case study defined with industry partners of DANSE.

Modeling Stochastic Adaptive Systems

In this section, we present the formal model used to encode behaviors of adaptive systems. In Section 2.1, we introduce Markov chains (MC) to represent individual components of a view. Then, in Section 2.2, we show how to describe views as well as relations between views, i.e., adaptive systems.

Discrete and Continuous Time Markov Chains

Definition 1. A (labelled) transition system is a tuple T = (Q, q, Σ, →, AP, L) where Q is a set of states, q ∈ Q is the initial state, Σ is a finite set of actions, →: Q × Σ × Q is the transition relation, AP is a set of atomic propositions, and L : Q → 2 AP is a state labelling function that assigns to each state the set of propositions that are valid in the state.

We denote by q a -→ q ′ the transition (q, a, q ′ ) ∈→. A trace is a finite or infinite alternating sequence of states and time stamps ρ = t 0 q 0 t 1 q 1 t 2 q 2 . . . , such that ∀i.∃a i ∈ Σ.q i ai -→ q i+1 . Time stamps measure the cumulative time elapsed from a time origin. In discrete time models delays are integer values (i.e., t 0 = 0, t 1 = 1, t 2 = 2) and therefore they can be omitted. In continuous time models they are real values. We denote by |ρ| the length of a trace ρ. If ρ is infinite then |ρ| = ∞. A trace is initial if q 0 = q and t 0 = 0. We denote by trace n (T ) (resp. trace(T )) the set of all initial traces of length n (resp. infinite traces) in T . Let 0 ≤ i ≤ |ρ|, we denote ρ |i = t 0 q 0 t 1 q 1 . . . t i-1 q i-1 the finite prefix of ρ of length i, ρ |i = t i q i t i+1 q i+1 . . . the suffix of ρ that starts at position i, and ρ[i] = q i the state at position i.

We now extend transition systems with probabilities to represent uncertainty of behaviors or of material on which the system is running. We present two semantics, either with discrete or continuous time, that are both compatible with our setting. A discrete time Markov chain (DTMC) is a state-transition system in which each transition is labelled by a probability P(s, s ′ ) to take the transition from state s to state s ′ . Definition 2. A (labelled) DTMC is a tuple D = (Q, q, Σ, →, AP, L, P) where:

-(Q, q, Σ, →, AP, L) is a labelled transition system, -P : Q×Q → [0, 1] is a transition probability matrix, such that q ′ ∈Q P(q, q ′ ) = 1 for all q ∈ Q, -→ is such that q a -→ q ′ iff P(q, q ′ ) > 0, and for each state q there is at most one action a ∈ Σ such that q a -→ q ′ for some q ′ . In continuous time Markov chains (CTMCs) transitions are given a rate. The sum of rates of all enabled transitions specifies an exponential distribution that determines a real value for the time spent in the current state. The ratio of the rates then specifies which discrete transition is chosen. Definition 3. A (labelled) CTMC is a tuple C = (Q, q, Σ, →, AP, L, R) where:

-(Q, q, Σ, →, AP, L) is a labelled transition system, -R : Q × Q → R ≥0 is a transition rate matrix, -→ is such that q a -→ q ′ iff
R(q, q ′ ) > 0, and there is a unique a ∈ Σ such that q a -→ q ′ . In our setting, a view of a system is represented by the combination of several components. We can compute the parallel composition C 1 C 2 of two DTMCs (resp. CTMCs) defined over the same alphabet Σ. Let (Q 1 , q 1 , Σ, → 1 , AP 1 , L 1 ) and (Q 2 , q 2 , Σ, → 2 , AP 2 , L 2 ) be the two underlying transition systems. We first compute their parallel composition, which is a labelled transition system (Q, q, Σ, →, AP, L), where Q = Q 1 × Q 2 , q = (q 1 , q 2 ), AP = AP 1 ∪ AP 2 , L(q) = L 1 (q 1 ) ∪ L 2 (q 2 ) and the transition relation → is defined according to the following rule:

q 1 a -→ 1 q ′ 1 q 2 a -→ 2 q ′ 2 (q 1 , q 2 ) a -→ (q ′ 1 , q ′ 2 ) (1)
Then, in case of DTMCs, the new transition probability matrix is such that P((q 1 , q 2 ), (q ′ 1 , q ′ 2 )) = P(q 1 , q ′ 1 ) * P(q 2 , q ′ 2 ), and in case of CTMCs the new transition rate matrix is such that R((q 1 , q 2 ), (q ′ 1 , q ′ 2 )) = R(q 1 , q ′ 1 ) * R(q 2 , q ′ 2 ). DTMCs with different alphabets can also be composed and they synchronize on common actions. However, if both DTMCs can perform a non synchronized action, a uniform distribution is applied to resolve the non determinism. In case of CTMCs, the two actions are in concurrence, such that if

q 1 a -→ 1 q ′ 1 with a ∈ Σ 2 , then (q 1 , q 2 ) a -→ (q ′ 1 , q 2 ) and R((q 1 , q 2 ), (q ′ 1 , q 2 )) = R(q 1 , q ′ 1 )
. In what follows, we denote by

Sys = C 1 C 2 • • • C n the DTMC (resp. CTMC) that results from the composition of the components C 1 , C 2 , . . . , C n .

Stochastic Adaptive Systems (SAS)

An adaptive system consists in several successive views. It starts in an initial view that evolves until it reaches a state in which an adaptation is possible. This adaptation consists in a view change that depends on a probability distribution that represents uncertainty of the environment. Definition 4. A SAS is a tuple (∆, Γ, S, sys, ) where:

-∆ = {C 1 , C 2 , . . . , C n }
is a set of DTMCs (resp. CTMCs) that are the components of the SAS. -Γ is the set of views of the SAS, such that each view is a stochastic system obtained from the parallel composition some components from ∆. -sys ∈ Γ is the initial view.

-S is the set of states of the SAS. S is the union of the states of each view in Γ , i.e., for each state s ∈ S there exists

{C 1 , C 2 , . . . , C k } ⊆ ∆ such that s ∈ Q 1 × Q 2 × • • • × Q k (where ∀i, 1 ≤ i ≤ k, Q i is the set of states of C i ). - ⊆ S × [0, 1] S is a set of adaptive transitions.
Observe that the number of components per state may vary. This is due to the fact that different views may have different components. Observe also that it is easy to add new views to an existing adaptive system without having to re-specify the entire set of views. An element (s, p) ∈ consists in a state s from a view sys ∈ Γ and a probability distribution p over the states in S. When s = s ′ , we denote s s ′ if there exists p such that (s, p) ∈ and p(s ′ ) > 0, which means that state s can be adapted into state s ′ with probability p(s ′ ).

A trace ρ in a SAS is either a finite combination of n traces ρ = ρ 0 ρ 1 . . . ρ n , such that for all 0 ≤ i ≤ n -1, ρ i = t 0i s 0i t 1i s 1i . . . t li s li is a finite trace of sys i ∈ Γ , and s li s 0i+1 , and t 0i+1 = 0, and ρ n is a finite or infinite trace sys n ∈ Γ . Otherwise ρ may be an infinite combination of finite traces ρ = ρ 0 ρ 1 . . . that satisfy for all i the same constraints.

A logic for SAS properties

Probabilistic Adaptive Bounded Linear Temporal Logic

We consider quantitative verification of dynamic properties that are expressed via a quantitative extension of the Adaptive Linear Temporal Logic (A-LTL) proposed in [START_REF] Zhang | Model-based development of dynamically adaptive software[END_REF]. Our logic, which we call Adaptive Bounded Linear Temporal Logic (A-BLTL), relies on an extension of Bounded Linear Temporal Logic (BLTL) combined with an adaptive operator. Although the logic is not strictly more expressive than BLTL, it is more suitable to describe properties of individual views, as well as global properties of the adaptive system, and it allows to develop specific algorithms for these properties. In the last part of the section, we also show how one can define contracts on such logic, where a contract [START_REF] Meyer | Applying "design by contract[END_REF] is a pair of assumptions/guarantees that must be satisfied by the system.

We first introduce BLTL, a logic used to express properties on individual views. In BLTL, formulas are built by using the standard Boolean connectors ∧, ∨, =⇒ , ¬, and the temporal operators G, F , X, U borrowed from Linear Temporal Logic (LTL). The main difference between BLTL and classical LTL is that each temporal modality is indexed by a bound k that defines the length of the run on which the formula must hold. The semantics of a BLTL formula is defined in Table 1 for finite executions of CTMC/DTMC ρ = t 0 s 0 t 1 s 1 t 2 s 2 . . . , In case formulas are nested, the value of k adapts incrementally. We now generalize BLTL to adaptive systems. For doing so, we introduce an adaptive operator in the spirit of [START_REF] Zhang | Using temporal logic to specify adaptive program semantics[END_REF]. The new logic A-BLTL is an extension of BLTL with an adaptive operator Φ Ω = ⇒ ≤k Ψ , where Φ is a BLTL formula, Ψ is an A-BLTL formula, Ω is a predicate over the states of different views of the SAS, and k is a time bound that limits the execution time of the adaptive transition. We will also consider unbounded versions of the adaptive operator.

ρ |= X ≤k Φ ≡ ∃i, i = max{j | t0 ≤ tj ≤ t0 + k} and ρ |i |= Φ ρ |= Φ1 U ≤k Φ2 ≡ ∃i, t0 ≤ ti ≤ t0 + k and ρ |i |= Φ2 and ∀j, 0 ≤ j ≤ i, ρ |j |= Φ1 ρ |= Φ1 ∧ Φ2 ≡ ρ |= Φ1 and ρ |= Φ2 ρ |= ¬Φ ≡ ρ |= Φ ρ |= P ≡ P ∈ L(s0) ρ |= true ρ |= f alse
Definition 5 (A-BLTL semantics). Let Φ Ω = ⇒ ≤k Ψ be an A-BLTL formula and ρ = t 0 s 0 t 1 s 1 t 2 s 2 . . . be an execution of the SAS:

ρ |= Φ Ω = ⇒ ≤k Ψ ≡ ∃i, i = min{j | t 0 ≤ t j-1 ≤ t 0 + k ∧ ρ |j |= Φ ∧ s j-1 s j ∧ Ω(s j-1 , s j )} ∧ ρ |i |= Ψ (2)
The property is unbounded if k = ∞, and in that case we write

Φ Ω = ⇒ Ψ .
According to Definition 5, an execution ρ satisfies an adaptive formula Φ Ω = ⇒ ≤k Ψ if and only if there exists a minimal prefix of ρ that satisfies Φ and reaches a state s j-1 , such that Ω(s j-1 , s j ) is satisfied, and such that the suffix of ρ from state s j satisfies Ψ . Therefore to satisfy this formula it is necessary to observe an adaptation compatible with Ω. We relax this constraint by introducing a new operator Φ Ω -→ ≤k Ψ , for which an adaptation is not necessary but triggers a check of Ψ when it happens. It is equivalent to the following formula:

Φ Ω -→ ≤k Ψ ≡ (Φ Ω = ⇒ ≤k true) =⇒ (Φ Ω = ⇒ ≤k Ψ ).
We finally introduce stochastic contracts, that are used to reason about both the adaptive system and its environment via assumptions and guarantees.

Definition 6 (Contracts for SAS).

A contract is defined as a pair (A, G), where A and G are respectively called the Assumption and the Guarantee. A SAS M satisfies the contract (A, G) iff ∀ρ, ρ |= A ⇒ ρ |= G, where ρ is a trace of M and ρ |= A (resp. G) means the trace ρ satisfies the assumption A (resp. the guarantee G). In that case we write M |= (A, G).

Verifying SAS Properties using SMC

SMC [START_REF] Younes | Verification and planning for stochastic processes with asynchronous events[END_REF][START_REF] Sen | On statistical model checking of stochastic systems[END_REF] is an alternative to model checking [START_REF] Baier | Principles of Model Checking[END_REF][START_REF] Clarke | Model checking[END_REF] that employs Monte Carlo methods to avoid the state explosion problem. SMC estimates the probability that a system satisfies a property using a number of statistically independent simulation traces of an executable model. By using results from the statistic area, SMC decides whether the system satisfies the property with some degree of confidence, and therefore it avoids an exhaustive exploration of the state-space of the model that generally does not scale up. It has already been successfully experimented in biology area [START_REF] Clarke | Statistical model checking in biolab: Applications to the automated analysis of t-cell receptor signaling pathway[END_REF][START_REF] Jha | A bayesian approach to model checking biological systems[END_REF][START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF], software engineering [START_REF] Clarke | On simulation-based probabilistic model checking of mixed-analog circuits[END_REF] as well as industrial area like aeronautics [START_REF] Basu | Statistical abstraction and model-checking of large heterogeneous systems[END_REF].

The basic algorithm used in SMC is the Monte Carlo algorithm. This algorithm estimates the probability that a system Sys satisfies a BLTL property P by checking P against a set of N random executions of SyS. The estimation p is

given by p = N 1 f (ρ i ) N
where f (ρ i ) = 1 if ρ i |= P , 0 otherwise. Using the formal semantics of BLTL, each execution trace ρ i is monitored in order to check if P is satisfied or not. The accuracy of the estimation increases with the number of monitored simulations. This accuracy can be controlled thanks to the Chernoff-Hoeffding bound [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF]. It relates N to δ and ε, that are respectively the confidence and the error bound of p:

P r(|p -p| < ε) ≥ 1 -δ if N ≥ ln( 2 δ ) 2ε 2 (3)
According to this relation, the user is able to trade off analysis time in return for accuracy of the result using the parameters δ and ε.

Knowing that there exists techniques to monitor BLTL properties [START_REF] Havelund | Synthesizing monitors for safety properties[END_REF], the model checking of A-BLTL is rather evident. Given an A-BLTL property Φ 1 Ω = ⇒ ≤k Φ 2 , the monitor will first check whether the run satisfies Φ 1 using classical runtime verification techniques. If no, then the property is not satisfied. If yes, then one checks whether there is a pair of two successive states between t 0 and t 0 + k that satisfies Ω. The latter is done by parsing the run. If this pair does not exist, then the property is not satisfied. Else we start a new monitor from the suffix of the run starting in the second state of the pair in order to verify Φ 2 .

Verifying unbounded SAS Properties using SMC

We propose a method inspired by [START_REF] Younes | Statistical verification of probabilistic properties with unbounded until[END_REF] to check unbounded A-BLTL properties. The principle is to combine a reachability analysis by model checking the underlying finite-state machine, with a statistical analysis of the stochastic model using the algorithms introduced previously.

We consider an A-BLTL property Φ Ω = ⇒ Ψ , where Φ and Ψ are BLTL formulas. We first consider the reachability problem with objective G = {s | ∃s ′ .s s ′ ∧ Ω(s, s ′ )}. The preliminary to the statistical analysis is to compute the set Sat(Reach(G)), that is all the states of the SAS that may eventually reach a state in G. This can be computed using classical model checking algorithms for finite-state machines. Only the underlying automata of the DTMCs or CTMCs of the SAS is used for this analysis. As stochastic quantities are ignored, efficient symbolic techniques exist to speed up this process [START_REF] Burch | Symbolic model checking: 1020 states and beyond[END_REF].

1: procedure Check(Φ Ω = ⇒ Ψ, ρ) 2:
if ¬Check(Φ, ρ) then return false 3:

end if 4:

i ← 0, prec ← null, curr ← null 5:

while i < |ρ| do 6:

prec ← curr, curr ← ρ[i] 7:

if curr ∈ Sat(Reach(G)) then return false 8:

end if 9:

if Ω(prec, curr) then return Check(Ψ, ρ |i ) 10:

end if 11:

i ← i + 1 12:

end while 13: end procedure Fig. 2. Algorithm to monitor unbounded A-BLTL properties Once this preliminary computation is performed, the Check algorithm from Figure 2 is used to monitor the runs of the SAS. The algorithm takes as input the A-BLTL property and a run ρ. The run should be in general infinite as there is no bound on the length of the runs that satisfied an unbounded A-BLTL property. In that case the states would be generated on-the-fly. The algorithm returns true or false, whether the run satisfied the property. We also denote Check(Φ, ρ) as the monitoring of the BLTL property Φ. Then the first step on line 2 is to monitor Φ on the run ρ. If the result is false then the property Φ Ω = ⇒ Ψ is not satisfied. Otherwise the algorithm searches through ρ for two states prec and curr such that Ω(prec, curr) is true. This is possible if curr belongs to the precomputed set Sat(Reach(G)). If Ω is satisfied the last step on line 10 is to monitor Ψ from the current position in ρ.

For homogeneous Markov chains (with constant probability matrices, as it is the assumption in this paper), the algorithm almost surely (with probability 1) terminates, since it either reaches a state where Ω is unreachable, or the probability to reach two states that satisfy Ω is not null. It can be iterated to check sequences of adaptive operators, that is to say properties where Ψ is also an unbounded A-BLTL.

A software engineering point of view

In this section, we propose high level formalisms to specify both adaptive systems and their properties. Then, we define semantics of those formalisms by exploiting the definitions introduced in the previous sections. This gives us a free verification technology for them. The situation is illustrated in Figure 3. PRISM toolset [START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF]. Due to space limit, the syntax common to RML and A-RML is only briefly described here 1 .

Adaptive RML systems as a high level formalism for SAS

We represent adaptive systems with Adaptive Reactive Module Language (A-RML), an extension of the Reactive Module Languages (RML) used by the

The RML language is based on the synchronisation of a set of modules defined by the user. A module is declared as a DTMC or CTMC, i.e., some local variables with a set of guarded commands. Each command has a set of actions of the form λi:ai where λ i is the probability (or the rate) to execute ai. A-RML extends RML such that each module can have some parameters in order to define its initial state. .

. . endmodule

The optional channel identifiers prefixing commands are used to strongly synchronise the different modules of a RML system. A module is synchronised over the channel chan if it has some commands prefixed by chan. We say a command is independent if it has no channel identifier.

In A-RML a system is a set of modules and global variables. The modules synchronise on common channels such that the system commands are the independent commands of each module and the synchronised commands. A command synchronised on chan forces all the modules that synchronised over chan to simultaneously execute one of their enabled commands prefixed by chan. If one module is not ready, i.e. it has no enabled commands for chan, the system has no enabled command over chan. Similarly to a module, if the system reaches a state with a non-deterministic choice, it is solved by a stochastic behaviour based on a uniform distribution. This solution allows to execute the A-RML system in accordance with the DTMC/CTMC models. system SYS_NAME ( < Parameters >) < global_variables > . . . < module_declarations > . . . endsystem An adaptive system consists in a set of different views, each represented by an A-RML system, and a list of adaptations represented by adaptive commands. The adaptive environment is used to specify which one is the initial view and what adaptations are possible. . . . endadaptive An adaptive command is similar to module command. It has a guard g k that applies to the current view SYS NAME, and a set of actions λi:ai where λ i is the probability (or the rate) to execute action a i = SYS NAME'(e 0 , . . . , e m ). This action defines the next view SYS NAME' after performing the adaptation. This view is determined according to the states of the previous view by setting the parameters of SYS NAME' with the expressions e 0 , . . . , e m evaluated over SYS NAME. The execution of the adaptive command is done in accordance with the SAS semantics.

Theorem 1. The semantics of A-RML can be defined in terms of SAS.

The proof of the above theorem is a direct consequence of the fact that semantics of RML is definable in terms of composition of MC, and that the definition an adaptive command can also be represented as a MC.

A contract language for SAS specification

The Goal and Contract Specification Language (GCSL) was first proposed in [START_REF] Arnold | Contracts and behavioral patterns for systems of systems: The EU IP DANSE approach[END_REF] to formalise properties of adaptive systems in the scope of the DANSE project. It has a strong semantics based on BLTL but it has a syntax close to the hand written English requirements. Dealing with formal temporal logic is often an issue to formalise correctly the initial English requirements. Most of the time the formalisation frequently contains some mistakes, which is due to the nesting of the temporal operators. The difficulty for correctly specifying properties is enough to make the overall methodology useless.

The GCSL syntax combines a subset of the Object Constraint Language (OCL) [START_REF]Ocl v2.2[END_REF] (used to define state properties, i.e., Boolean relations between the system components) and English behavioural patterns used to express the evolution of these state properties during the execution of the system. The usage of OCL is illustrated in Example 1.

Example 1. We consider a SAS describing the implementation of an emergency system in a city. The city area is divided as a set of districts where each district may have some equipment to fight against the fire, e.g. some fire stations with fire brigades and fire fighting cars. Each district is also characterised by a risk of fire and the considered damages are mainly related to the population size of each district. The requirement "Any district cannot have more than 1 fire station, except if all districts have at least 1" ensures the minimal condition for the equipment distribution in the city. We use syntactic coloration to make the difference between the parts of the language used in the property: the words in red are identifiers from the model, the blue part is from OCL, like collection handling, and the black words are variables:

City.itsDistricts→exists(district | district.ownedFireStations > 1) implies City.itsDistricts→forAll(district | district.ownedFireStations ≥ 1)
GCSL patterns are used to specify temporal properties. In this section we only present a subset of such patterns that is considered to be general enough to specify properties of a large set of industry-examples from the DANSE project. After having read this section, the user shall understand that the set can be easily increased. Each pattern can nest one or more state properties, denoted in the grammar by the non-terminals <OCL-prop> and <arith-rel>, that respectively denote a state property written in OCL or an arithmetic relation between the identifiers used in the model. The non-terminal <int> denotes a finite time interval over which the temporal pattern is applied, and <N> is a natural number. The patterns can be used directly or combined with OCL: applying a pattern to a collection of system components defines a behavioural property that is applied to each element of the collection. We present below an excerpt of the complete GCSL grammar available in [START_REF] Arnold | Contracts and behavioral patterns for systems of systems: The EU IP DANSE approach[END_REF]: We now propose A-GCSL, a syntax extension for GCSL that can be used to describe adaptive requirements of SAS. A-GCSL extends the GCSL grammar by adding a new pattern that allows to express adaptive relations as done with the two adaptive operators defined in Section 3. The first pattern of <dyna-spec> is equivalent to the operator Ω = ⇒ and the second one denotes Ω -→. Any adaptive requirement has three elements (A, Ω, G) that are called assumption, trigger and guarantee, respectively. The assumption and guarantee are specified in GCSL, whereas the trigger is in OCL. The syntax allows to compose the patterns by specifying the guarantee with an adaptive pattern. For instance, a composed requirement of the form if Φ 1 holds and for all rule that satisfies Ω then (if Φ 2 holds and for all rule that satisfies Ω ′ then Φ 3 holds) holds is equivalent to the property Φ 1 The attribute underFire denotes the number of districts in which a fire has been declared. If there are more than three fires in the city, then the fire stations change their usual emergency management into a crisis one. When such management is activated, the firemen have 50 hours to fix the problem. The requirement can be translated in A-GCSL using the following formula:

Ω = ⇒ Φ 2 Ω ′ = ⇒ Φ 3 . The A-GCSL
Φ6 = `underF ire = 0 ´underF ire≥3 ---------→ ≤10000 di :district " ¬di.decl =⇒ G ≤10000 (di.decl =⇒ F ≤50 di.f ire = 0) "
In [START_REF] Arnold | Contracts and behavioral patterns for systems of systems: The EU IP DANSE approach[END_REF], we have showed that any GCSL pattern can be translated into a BLTL formula. The result extends as follows.

Theorem 2. Any A-GCSL pattern can be translated into an A-BLTL property.

This result is an immediate consequence of the definition of the adaptive pattern.

Experiments with SAS

Our work has been implemented in a new statistical model checker named Plasma-Lab [START_REF] Boyer | Plasma-lab: A flexible, distributable statistical model checking library[END_REF], a platform that includes efficient SMC algorithms, a graphical user interface (GUI) and multiple plugins to support various modelling languages. The tool is written in Java that offers maximum cross-platform compatibility. The GUI allows to create projects and experiments, and it implements a simple and practical mean of distributing simulations using remote clients. Plasma-Lab also provides a library to write new plugins and create custom statistical model checkers based on arbitrary modelling languages. Developers will only need to implement a few simple methods to build a simulator and a logic checker, and then beneficiate from Plasma-Lab SMC algorithms.

Plasma-Lab can be used as a standalone application or be instantiated within other softwares. It has already been incorporated in various performancecritical softwares and embedded hardware platforms. The current plugins allow to simulate biologic models, models written in RML and A-RML, but it has also the capabilities to drive an external engine to perform the simulations, like MATLAB, Scilab, or DESYRE a simulator for adaptive systems developed by Ales [START_REF] Corp | Advanced laboratory on embedded systems[END_REF].

CAE model

Together with our industrial partners in the DANSE project, we have developed the Concept Alignment Example (CAE). The CAE is a fictive adaptive system example inspired by real-world Emergency Response data to a city fire. It has been built as a playground to demonstrate new methods and models for the analysis and visualization of adaptive systems designs.

The CAE describes the organization of the firefighting forces. We consider in our study that the city is initially divided into 4 districts, and that the population might increase by adding 2 more districts. Different and even more complex examples can be built using the components of this design.

A fire station is assigned to the districts, but as the fire might spread within the districts, the system can adapt itself by hiring more firemen. We can therefore design a SAS with three views as described in Figure 4. Adaptive transitions exist between these views to reflect changes in the environment and adaptations of the system. Initially in View 1, the system can switch to View 2 when the population of the city increase. This change models an uncertainty of the environment, and for the purpose of this study we fix its probability to 0.01 Then, if the number of fires becomes greater than 2, the system adapts itself by switching to View 3. If the number of fires eventually reduces and becomes lower than 2, the system might return to View 2. Again, as this change is uncertain, we fix its probability to 0.8. We design several A-RML models of the system that consist in two types of modules: District and FireStation, both based on a CTMC semantics. First, we study a model AbstractCAE that is an abstract view of the SAS. In this model, the District module, presented below, is characterized by a constant parameter p, that determines the probability of fire, and by two Boolean variables decl and men, that respectively defines if a fire has been declared and if the firemen are allocated to the district. The module fireStation has one constant parameter distance i for each module of the system. This parameter determines the probability to react at a fire, such that the greater the distance, the lower the probability. However a fire station can only treat one fire at a time, which is encoded with a Boolean variable allocated. The fire stations and the districts synchronize on channels allocate and recover, that respectively allocate firefighters to the district and bring them back when the fire is treated. The different views are constructed by instantiating and renaming the modules presented above. -ConcreteCAE Full is the full model of ConcreteCAE, with the 3 views and all the adaptive transitions between them.

Checking requirements

The requirements are expressed in A-GCSL and translated to A-BLTL. We first check the model AbstractCAE against A-BLTL properties with adaptive operators. Our goal is to verify that the transitions between the different views of the system occurs and satisfy some properties. The first property, if [true] holds then there exists a rule satisfying [underfire ≤ 1] and Always [!maxfire], checks that when the system is in View 1, it eventually switches to View 2 when the number of districts that have declared a fire (underfire) is still lower than 1, and that as a result the system remains safe for a limited time period, i.e., the number of districts that have declared a fire is not maximum (maxfire is false). To check this property we limit the analysis to the model AbstractCAE 1 with only View 1 and View 2. The A-GCSL property is translated in an A-BLTL formula: Φ 1 = true underfire≤1 =======⇒ G ≤1000 !maxfire, and the results in Table 2 show that the probability to satisfy the property is only 50%. This justify the need to add a second fire station, as in View 3.

The second property, if [true] holds then there exists a rule satisfying [true] and Always [!maxfire], checks that from View 2 a second fire station is quickly added, which switches the system to View 3, and that then the system is safe. The property is checked on the model AbstractCAE 2 using the A-BLTL formula : Φ 2 = true true = = ⇒ ≤100 G ≤10000 !maxfire. Finally, with the property if [true] holds then there exists a rule satisfying [true] and [true], we check that from View 3 the system eventually returns to View 2. Therefore we use the model AbstractCAE 3 that starts in View 3 and we check the A-BLTL formula Φ 3 = true true = = ⇒ ≤100 true. The AbstractCAE models are simple enough to be able to perform reachability analyses and check the unbounded A-BLTL properties presented above using Algorithm 2. In a second step we consider the models ConcreteCAE to better evaluate the safety of the system. The state spaces of these models contain several millions of states, and therefore, they can only be analyzed by purely SMC algorithms. We verify the two following properties:

-Always !maxfire, to check that the maximum of fire intensity of 10 is never reached in any district. This corresponds to Φ 4 = G ≤10000 !maxfire. -Whenever [fire > 0] occurs [fire = 0] within [50 hours], to check that a fire in a district is totally extinct within 50 hours. This corresponds to Φ 5 = G ≤10000 d6.fire > 0 =⇒ F ≤50 d6.fire = 0 .

These two properties are first checked for each view of the system. The results in Table 2 show that while View 1 and View 3 are surely safe, View 2 is frequently unsafe. But when we check these properties on the complete adaptive model ConcreteCAE Full, with the three views, we can show that the system remains sufficiently safe. It proves that after a change of the environment (the increase of population) the system is able to adapt itself to guaranty its safety.

In the last experiment of Table 2 we check the A-GCSL property presented in Example 3. This bounded adaptive A-GCSL property is checked using the full ConcreteCAE model.

We have performed each experiment in Plasma-Lab with a confidence δ = 0.01 and an error bound ε = 0.02. The results in Table 2 give the probabilities estimation and the time needed to perform the computation.

Conclusion

This paper presents a new methodology for the rigorous design of stochastic adaptive systems. Our model is general, but the verification procedure can only reason on a finite and known set of views. Our formalism is inspired from [START_REF] Zhang | Using temporal logic to specify adaptive program semantics[END_REF], where both the stochastic extension and high level formalisms are not considered.

In future work, we will extend this approach to purely dynamic systems. Another objective is to extend the work to reason about more complex properties such as energy consumption. Finally, we shall exploit extensions of SMC algorithms such as CUSUM [START_REF] Page | Continuous inspection schemes[END_REF] which permits to reason on switches of probability satisfaction. This would allow us to detect emergent behaviors.

Fig. 3 .

 3 Fig. 3. SAS verification flow

  module MOD_NAME ( < Parameters >) < local_vars > . . . [ chan ] g k -> (λ0 : a0 ) + . . . + (λn : an );

  adaptive init at SYS_NAME ( < Initial values >) . . . { SYS_NAME | g k } -> λ0 :{ a0 } + . . . + λn :{ an };

  <GCSL> ::= <OCL-coll>->forAll(<variable>| <pattern>) | <OCL-coll>->exists(<variable>| <pattern>) | <OCL-prop> | <pattern> <pattern> ::= whenever [<prop>] occurs [<prop>] holds during following [<int>] | whenever [<prop>] occurs [<prop>] implies [<prop>] during following [<int>] | whenever [<prop>] occurs [<prop>] does not occur during following [<int>] | whenever [<prop>] occurs [<prop>] occurs within [<int>] | [<prop>] during [<int>] raises [<prop>] | [<prop>] occurs [<N>] times during [<int>] raises [<prop>] | [<prop>] occurs at most [<N>] times during [<int>] | [<prop>] during [<int>] implies [<prop>] during [<int>] then [<prop>] during [<int>] <prop> ::= <OCL-prop> | <arith-rel> Example 2. Consider the following requirement about the model described in Example 1: "The fire fighting cars hosted by a fire station shall be used all simultaneously at least once in 6 months". This requirement uses both GCSL and OCL patterns: City.itsFireStations→forAll(fStation | Whenever [fStation.hostedFireFighting-Cars → exists(ffCar | ffCar.isAtFireStation)] occurs, [fStation.hostedFireFight-ingCars→forall(ffCar | ffCar.isAtFireStation = false)] occurs within [6 months])

Example 3 .

 3 grammar is the following: <dyna-spec> ::= if [<GCSL>] holds and for all rule that satisfies [<prop>] then ( <GCSL> | <dyna-spec> ) holds | if [<GCSL>] holds then there exists a rule satisfying [<prop>] and ( <GCSL> | <dyna-spec> ) holds Consider again the system in Example 1 and the following A-GCSL requirement: if [ City.underFire = 0 ] holds and for all rule such that rule satisfies [ City.underFire ≥ 3 ] then [ City.itsDistricts→forall(district | district.decl = false =⇒ whenever [ district.decl = true ] occurs, [ district.fire = 0 ] occurs within [50 hours] ) ] holds

Fig. 4 .

 4 Fig. 4. Components and Views in the CAE model

  module District ( const int p ) decl : bool init false ; men : bool init false ; [] ! decl -> p /1000: ( decl '= true ); [ allocate ] decl & ! men -> ( men '= true ); [ recover ] decl & men -> 1/ p : ( decl '= false ) & ( men '= false ); endmodule We refine this model to better encode the behaviour of the SAS. In this new model ConcreteCAE a new variable fire of module District ranges from 0 to 10 and grades the intensity of the fire. The fire stations can now assign several cars (from 0 to 5) to each districts. Therefore the variables men and allocated becomes integers. module District ( const int p ) fire : [0..10] init 0; decl : bool init false ; men : [0..5] init 0; [] fire =0 -> p /1000: ( fire '=1); [] fire >0 & fire <10 -> p /((1+ men )*100): ( fire '= fire +1); [] fire >0 & ! decl -> ( fire * fire )/10: ( decl '= true ); [ allocateSt1 ] decl & fire >0 -> ( fire * fire )/10: ( men '= men +1); [ allocateSt2 ] decl & fire >0 -> ( fire * fire )/10: ( men '= men +1); [] men >0 & fire >0 -> men /10: ( fire '= fire -1); [ recover ] decl >0 & fire =0 -> 1000: ( men '=0)&( decl '= false ); endmodule From the two models we can consider several subparts composed by one or several views of the SAS. Adaptive commands are used to model the transitions between the different views. -AbstractCAE 1 consists in View 1 and 2 from model AbstractCAE. -AbstractCAE 2 consists in View 2 and 3. -AbstractCAE 3 has the same views as AbstractCAE 2 but is initiated in View 3 instead of View 2. -ConcreteCAE 1 only consists in View 1 from model ConcreteCAE. -ConcreteCAE 2 only consists in View 2. -ConcreteCAE 3 only consists in View 3.

Table 1 .

 1 Semantics of BLTL with |ρ| ≥ k. If |ρ| < k, it is extended by duplicating the last state enough times.

Table 2 .

 2 Experiments on CAE models

	Property	CAE Model	Estimation interval Consumed Time
	Φ1	AbstractCAE 1 View 1, View 2	[0.53, 0.56]	1351s
	Φ2	AbstractCAE 2 View 2, View 3	[0.84, 0.86]	11s
		AbstractCAE 3
	Φ3	AbstractCAE 2 starting	[0.98, 1]	1363s
		from View 3
	Φ4	ConcreteCAE 6 4 dist. 1 sta.	[0.95, 0.99]	11s 9s
	Φ4	ConcreteCAE 2	[0.46, 0.5]	15s
	Φ5	6 dist. 1 sta.	[0.21, 0.25]	13s
	Φ4	ConcreteCAE 3	[0.98, 1]	30s
	Φ5	6 dist. 2 sta.	[0.98, 1]	31s
	Φ4	ConcreteCAE Full	[0.89, 0.93]	25s
	Φ5	4-6 dist. 1-2 sta.	[0.82, 0.86]	42s
	Φ6	ConcreteCAE Full 4-6 dist. 1-2 sta.	[0.47, 0.51]	109s

The full syntax can be found at http://project.inria.fr/plasma-lab/


Hrishikesh Deshpande

Pierre Maurel

Christian Barillot

Detection of Multiple Sclerosis Lesions using Sparse Representations and Dictionary Learning

The manual delineation of Multiple Sclerosis (MS) lesions is a challenging task pertaining to the requirement of neurological experts and high intra-and inter-observer variability. It is also time consuming because large number of Magnetic Resonance (MR) image slices are needed to obtain 3-D information. Over the last years, various models combined with supervised and unsupervised classification methods have been proposed for segmentation of MS lesions using MR images. Recently, signal modeling using sparse representations (SR) has gained tremendous attention and is an area of active research. SR allows coding data as sparse linear combinations of the elements of over-complete dictionary and has led to interesting image recognition results. The dictionary used for sparse coding plays a key role in the classification process. In this work, we have proposed to learn class specific dictionaries and developed a new classification scheme, to automatically detect MS lesions in 3-D multi-channel MR images.

Introduction

Multiple sclerosis is a chronic, autoimmune disease of the central nervous system, characterized by structural damages of axons and their myelin sheathes. During progression of the disease, certain areas of brain develop MS lesions. The evolution of MS lesions is highly variable and is not fully known. MS is more common in North America and Europe and is more prevalent in young adult population, causing non-traumatic disabilities.

Magnetic Resonance Imaging (MRI) holds the capability of detecting abnormalities in 95% of the patients with MS and is the best paraclinical method for imaging MS [START_REF] Grossman | Perspectives on Multiple Sclerosis[END_REF]. These images are analyzed to find the number and spatial patterns of the lesions, appearance of new lesions and the total lesion load, which are key parameters in the current MS diagnostic setup. However, manual segmentation of MS lesions is a laborious and time consuming process and is prone to high intra-and inter-expert variability. Therefore, there is a need for fully automated MS lesion detection methods that can handle large variety of MR data and which can provide results that correlate well with expert analysis [START_REF] Polman | Diagnostic criteria for multiple sclerosis: 2010 Revisions to the McDonald criteria[END_REF].

Over the past years, various approaches for semi-automatic and automatic segmentation of MS lesions have been proposed. In these methods, different image features, classification methods and models have been tried, but they usually suffer from high sensitivity to the imaging protocols and so usually require tedious parameter tuning or specific normalized protocols [START_REF] García-Lorenzo | Review of automatic segmentation methods of multiple sclerosis white matter lesions on conventional magnetic resonance imaging[END_REF]. More recently, sparse representation has evolved as a model to represent an important variety of natural signals using few elements of an overcomplete dictionary. Many publications have demonstrated that sparse modeling can achieve state-of-the-art results in image processing applications such as denoising, texture segmentation and face recognition [START_REF] Li | An Efficient Dictionary Learning Algorithm and Its Application to 3-D Medical Image Denoising[END_REF][START_REF] Wright | Robust Face Recognition via Sparse Representation[END_REF]. In [START_REF] Wright | Robust Face Recognition via Sparse Representation[END_REF], given multiple images of individual subjects under varying expressions and illuminations, the images themselves were used as dictionary elements, for classification. Such a method uses dictionary learning to analyze image as a whole. Mairal et al [START_REF] Mairal | Discriminative learned dictionaries for local image analysis[END_REF] proposed to learn discriminative dictionaries better suited for local image descrimination tasks. In medical imaging, local image analysis is of prime importance and it could be interesting to see the performance of sparse representation and dictionary learning based classification methods in the context of disease detection. Some researchers have reported works on segmentation of endocardium and MS lesions using dictionary learning [START_REF] Weiss | Multiple Sclerosis Lesion Segmentation Using Dictionary Learning and Sparse Coding[END_REF][START_REF] Rosas-Romero | Segmentation of endocardium in ultrasound images based on sparse representation over learned redundant dictionaries[END_REF]. Weiss et al. proposed an unsupervised approach for MS lesion segmentation, in which a dictionary learned using healthy brain tissue and lesion patches is used as basis for classification [START_REF] Weiss | Multiple Sclerosis Lesion Segmentation Using Dictionary Learning and Sparse Coding[END_REF].

Our approach differs from this method in several ways. In [START_REF] Weiss | Multiple Sclerosis Lesion Segmentation Using Dictionary Learning and Sparse Coding[END_REF], authors use only FLAIR MR images for analysis of clinical data. However, MS lesions appear in different intensity patterns in various MR sequences, which include T1 (T1-w MPRAGE) and T2-weighted, (T2-w) and Proton Density (PD). The complementary information in these MR images can further assist in classifying MS lesions. We build our analysis using above mentioned MR sequences. Our major contribution is however that we learn class specific dictionaries for healthy brain tissues and lesions that promote the sparse representation of healthy and lesion patches. The lesion patches are well adapted to its own class dictionary, as opposed to the other. Thus, we can use the reconstruction error derived from sparse decomposition of test patch on to these dictionaries for classification. In the dataset, the healthy class patches outnumber lesion patches and exhibit more variability. Thus we use different dictionary lengths for modelling individual class patches. In this manner, we take into consideration the data variability and class imbalance in healthy and lesion classes. Finally, supervised approach for detection results in omission of tuning of one parameter as mentioned in [START_REF] Weiss | Multiple Sclerosis Lesion Segmentation Using Dictionary Learning and Sparse Coding[END_REF], making this method fully automatic. In the following sections, we describe our new approach and provide its evaluation using clinical images.

Methodology

As shown in Figure 1, we first preprocess MR images for noise removal and then extract the image patches of predefined size using brain mask. These patches are normalized and are divided into the training and test sets for healthy brain tissue and lesion classes, with the help of manual segmentation images. Using training signals, we derive different classification approaches by either learning single dictionary or two separate dictionaries for both the classes. Finally, for a given test patch, the reconstruction error based classification method is developed, followed by voxel-wise classification and lesion detection. The following subsections briefly describe these steps.

Patch Extraction and Training Set

We divide the intracranial MR volume into several 3-D patches and flatten them into one dimensional concatenated vectors representing intensities of T1-w MPRAGE, T2-w, PD and FLAIR images. Keeping the computational complexity of further analysis in mind, we extract a patch every M voxels in each direction. As described earlier, we develop supervised approach by labelling these patches as belonging to either healthy or lesion class. If, in a patch, the number of voxels manually labelled as lesions exceeds a threshold T L = 6mm 3 , it is included in a lesion set, or in healthy set otherwise. For every subject, we obtain around 1.5 × 10 6 patches for healthy and 10 3 to 10 5 patches for lesion class, depending on the lesion load for each patient. These patches are finally normalized to limit their individual norms below or equal to unity, as per constraint imposed by dictionary learning.

Sparse Representations and Dictionary Learning

Sparse representation of the data allows the decomposition of signal into linear combination of few basis elements in an overcomplete dictionary. Consider a signal x ∈ R N and an overcomplete dictionary D ∈ R N ×K . The sparse coding problem can be stated as min a a 0 s.t. x = Da or x -Da Basis pursuit algorithm solves the convex approximation of the problem above by replacing l 0 norm with l 1 norm that also results in sparse solution [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. Thus, the sparse coding problem can be given by min

a x -Da 2 2 + λ a 1 , (1)
where λ controls the trade-off between representation error and sparsity.

The fixed dictionaries like wavelets can be efficient if a background analytical model can be inferred. On the other hand, the dictionary learning from underlying data has produced exciting results with greater data adaptibility and has replaced the use of generic models. For a set of signals {x i } i=1,.,m , the dictionary learning problem is to find D such that each signal can be represented by sparse linear combination of its atoms. This can be stated as the following optimization problem min

D,{ai} i=1,..,m m i=1 x i -Da i 2 2 + λ a i 1 .
(

) 2
The optimization is carried out as two-step process involving the sparse coding step with fixed D and the dictionary update step with fixed a.

Classification

(a) Using Single Dictionary : In the context of MS lesion classification, the simplest idea, similar to [START_REF] Weiss | Multiple Sclerosis Lesion Segmentation Using Dictionary Learning and Sparse Coding[END_REF], could be to use a single dictionary learned from healthy and lesion class patches. As the lesions are outliers with respect to the healthy brain intensities, the decomposition of lesion patch using this dictionary would result in higher representation error than that for the healthy tissue patch. For a given test patch, we calculate the sparse coefficients and reconstruction error, and assign it to the lesion class if this error is greater than chosen threshold. The threshold is selected by observing the histogram of the error map.

(b) Using Class Specific Dictionaries (Same Length) : Here, we learn class specific dictionaries D 1 and D 2 for healthy and lesion classes, respectively. Given a test patch x ∈ R N , classification is performed in two steps: In the first step, sparse coefficients a i are obtained using Eq (1) for each class i=1 (Healthy) and 2 (Lesion). The test patch is then assigned to class c such that

c = argmin i x -D i a i 2 2 . (3)
(c) Using Class Specific Dictionaries (Different Lengths) : The dictionaries learned using above mentioned approach does not take into account the data variability between two-classes. The size of the dictionary plays a major role in the data representation. For healthy class data with more variability and number of training samples than that for the lesion class, we allow larger dictionary length for healthy class data and study its effect on MS lesion classification.

Voxel-wise Classification and Lesion Detection

As already stated, there is some overlap between patches. However, to obtain voxel-wise classification, each voxel needs to be assigned to either of the classes. This is achieved using majority voting, in which, the voxel under consideration is classified as healthy or lesion, using majority votes of all patches which contain that voxel.

The voxelwise classification image is further processed to obtain the lesion based detection image. A lesion is said to be detected if

R D T R GT R GT ≥ T O ,
where R D and R GT are respectively the candidate regions in the classification image and the ground truth, whereas T O is the threshold indicating overlap between them as a fraction of ground truth lesion.

Dataset and Preprocessing

The proposed approach was validated on MRI volumes of 14 MS patients acquired by Verio 3T Siemens scanner. T1-w MPRAGE, T2-w, PD and FLAIR MR modalities were chosen for the experiment. The volume size for T1-w MPRAGE and FLAIR is 160 × 256 × 256 and voxel size is 1 × 1 × 1 mm 3 , whereas for T2-w and PD scans, the volume size is 192 × 256 × 44 and voxel size is 1 × 1 × 3 mm 3 . The manual segmentation images obtained from neurological experts are referred to as ground truth lesion masks.

For MR images of each patient, the imaging artifacts are corrected by denoising using non-local means and Intensity Inhomogenity Correction (IIH). The images so obtained are then registered with respect to T1-w MPRAGE volume and are processed further to extract the intra-cranial region.

Results and Discussions

We implemented our method using MATLAB and Python. The packages AN-IMA and N3 ITK were used for denoising, registration and IIH correction, respectively [START_REF] Coupé | An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[END_REF][START_REF] William | Multi-modal volume registration by maximization of mutual information[END_REF][START_REF] Tustison | N4ITK: Improved N3 Bias Correction[END_REF]. We used the neuroimaging software Brain Extraction Tool (BET) for brain extraction [START_REF] Smith | Fast robust automated brain extraction[END_REF]. For dictionary learning and sparse coding, we used SPArse Modeling Software (SPAMS) package [START_REF] Mairal | Online Dictionary Learning for Sparse Coding[END_REF].

We performed the experiments on 14 subjects using Leave-One-Subject-Out-Cross-Validation. Different parameters have been tested for the methods. It was found that image patch of size 5 × 5 × 5, with a patch every 2 voxels in each direction, was optimal with respect to the classification efficiency. The dictionary length of 5000 and sparsity parameter λ = 0.95 were optimal selections for dictionary learning method. For voxel-wise classification method, we then recorded the number of voxels that belong to True Positives (TP), False Negatives (FN), False Positives (FP) or True Negatives (TN) and the classification methods were finally validated by calculating sensitivity= T P T P +F N and Positive Predictive Value (PPV) = T P T P +F P .

In the first method, we studied the classification by learning single dictionary with the help of both healthy brain tissue and lesion patches. We chose sparse penalty factor λ = 0.85 in the sparse coding step and performed the classification for various threshold values on the histogram of error map, as explained previously. We then selected the threshold for which the best voxelwise classification results were obtained in terms of both sensitivity and PPV. It was observed that the method suffered with a very large number of false positive detections.

Next, we learned class specific dictionaries for healthy and lesion classes, each. We used dictionary lengths of 5000 for signal representation of each class. The mean sensitivity and PPV obtained using this approach were 91.5% and 7.5%. This method performs better than the previous method but still contains many false positives. The primary reason behind this can be the difference in the data variability of each class signals. The healthy class patches have more variability in terms of representation of white matter (WM), gray matter (GM) and cerebrospinal fluid (CSF), as compared to the variations in the representation of lesions. Hence, we adopted different dictionary lengths for representation of these classes. We used dictionary lengths of 5000 and 1000 respectively, for healthy and lesion classes. Table 1 summarizes the results of the voxelwise classification for the three methods described above. It can be seen that using class specific dictionaries with the same dictionary length improves both sensitivity and PPV, as compared to the first method. But PPV in the second method is still low, indicating that there are still large number of false positives, which can explain higher sensitivity. Using different dictionary lenghts, as implemented in third method, drastically reduces the number of false positives, which can be seen by the significant increment in PPV, while keeping the sensitivity in the acceptable limit.

The mean PPV and sensitivity for lesion detection with class specific dictionaries of different lengths are shown in Table 2 for various overlap thresholds T O . To be consistent with the threshold T L incorporated in learning stage (Refer Section 2.1), we ignore very small lesions with volumes less than T L . It can be seen that we detect 61% of the lesions with the overlap threshold of 1% . Moreover, in 49% of the lesions detected, at least 40% of the voxels are correctly classified by the method. In Figure 2, we show the results for patient 8, for all the methods discussed above. The detection image is superimposed on FLAIR MR image. It can be observed that methods (a) and (b) have large number of false positives. We get the best classification results using class specific dictionaries with different dictionary lengths. But, in terms of voxelwise classification, there are still few false positives and true negatives around the actual lesion. This does not pose a major problem for lesion detection as long as significant portion of the actual lesion is being classified correctly. There are, however, some false positive lesion detections.

We are aware that we do not have a very large population for training. Hence we investigated the incorporation of longitudinal database into our analysis by considering MR sequences at 3 time points (M0, M3 and M6) for all the patients. As the lesions evolve over the course of time, it is fair to consider that each new dataset will enrich our learning model. Thus, we modified the training data, for each patient, in two ways: (1) Data at time-points M0 and M3, with 26 datasets and (2) Data at time-points M0, M3 and M6, with 39 datasets. However, the lesion detection experiments for the same test subjects, as in previous experiments, using class specific dictionaries with the lengths of 5000 and 1000 for healthy and lesion class respectively, did not show any significant improvement in the sensitivity and PPV. This suggests that the population for training the dictionaries earlier was sufficient and the dictionaries should be adapted to learn 1.

more specific structures viz. WM, GM and CSF versus lesions to help improve the detection.

Conclusion

In this paper, we have proposed a new supervised approach to automatically detect multiple sclerosis lesions using dictionary learning. We investigated the performance of three methods which either use one dictionary, treating lesions as outliers, or use class specific dictionaries for healthy and lesion classes, wherein the underlying data for each class is represented by the dictionary and sparse coefficients. We further studied the effect of using different dictionary lengths, allowing larger dictionaries to represent the complex data and concluded that such method minimizes the false positive detections in the classification. Although the method using class specific dictionaries follows supervised approach, contrary to the single dictionary based classification method, which does not necessarily require training data, it is worth mentioning that the former method eliminates one parameter: threshold on error map. This crucial paramater is not easy to tune and could lead to worse classification results for small errors in the brain extraction procedure.

To further improve the results, it would be interesting not to learn only one dictionary for healthy brain tissues, but derive dictionaries more specific to WM, GM and CSF, in addition to the lesions. One could also study the role of sparse coefficients in addition to the reconstruction error, in the classification step.

Fig. 1 .

 1 Fig. 1. Flowchart of MS Lesion Detection using Dictionary Learning

Fig. 2 .

 2 Fig. 2. Classification results for Patient 8. For illustration purpose, one slice has been arbitrarily selected. True Positives are in red, False Positives are in cyan, False Negatives are in green. Methods (a), (b) and (c) are the same as in Table1.

Table 1 .

 1 Voxel-wise classification results using: (a) Single Dictionary (SD), with 5000 atoms learned using healthy and lesion class data, (b) Class Specific Dictionaries with Same Lengths (CSD SL): 5000 atoms each and (c) Class Specific Dictionaries with Different Lengths (CSD DL): 5000 atoms for healthy class and 1000 atoms for lesion class. Sensitivity and Positive Predictive Value (PPV) (%) are given for each method.

		(a)		(b)		(c)
	Patient	SD		CSD SL	CSD DL
		Sens.	PPV	Sens.	PPV	Sens.	PPV
	1	42	1	97	3	53	31
	2	74	1	98	2	66	41
	3	73	1	91	2	63	27
	4	91	2	98	17	57	68
	5	61	1	95	10	54	65
	6	91	7	89	29	38	55
	7	78	1	85	3	20	32
	8	72	1	98	3	69	21
	9	80	1	36	2	4	9
	10	66	1	97	9	61	52
	11	89	2	98	12	66	41
	12	75	1	99	8	52	36
	13	78	1	100	3	77	31
	14	59	1	100	2	78	17
	Average	73.5	1.57	91.5	7.5	54.14	37.57

Table 2 .

 2 Performance analysis for lesion detection using Class Specific Dictionaries with Different Lengths (CSD DL) for each class, with 5000 atoms for healthy class dictionary and 1000 atoms for lesion class dictionary.

		TO = 0.01 TO = 0.1 TO = 0.2 TO = 0.3 TO = 0.4
	PPV (%)	61.67	58.41	56.53	54.31	49.40
	Sensitivity (%)	60.97	57.58	56.67	54.56	49.94

2 ≤ ε, where a 0 is l 0 norm of the sparse coefficient vector a ∈ R K and ε is error in representation.