Sertaç Karahoda
email: skarahoda@sabanciuniv.edu

Osman Tufan Erenay

Kamer Kaya
email: kaya@sabanciuniv.edu

Uraz Cengiz Türker

Hüsnü Yenigün
email: yenigun@sabanciuniv.edu

Parallelizing Heuristics for Generating Synchronizing Sequences

Synchronizing sequences are used in the context of finite state machine based testing in order to initialize an implementation to a particular state. The cubic complexity of even the fastest heuristic algorithms known in the literature to construct a synchronizing sequence can be a problem in practice. In order to scale the performance of synchronizing heuristics, some algorithmic improvements together with a parallel implementation of these heuristics are proposed in this paper. An experimental study is also presented which shows that the improved/parallel implementation can yield a considerable speedup over the sequential implementation.

Introduction

Model Based Testing (MBT) uses formal models of system requirements to generate effective test cases. Most MBT techniques use state-based models, where the behaviour of the model is described in terms of states and state transitions. There has been much interest in testing from finite state machines (FSMs) (e.g., see [START_REF] Chow | Testing software design modeled by finite-state machines[END_REF][START_REF] Hennie | Fault-detecting experiments for sequential circuits[END_REF][START_REF] Ural | On minimizing the lengths of checking sequences[END_REF][START_REF] Hierons | Reduced length checking sequences[END_REF][START_REF] Petrenko | Testing from partial deterministic FSM specifications[END_REF][START_REF] Simão | On reducing test length for FSMs with extra states[END_REF]). Common to most FSM based testing methods is the need to bring the system under test (SUT) to a particular state. When there is a trusted reset input in the SUT, this is quite easy. However, sometimes such a reset input is not available, or even if it is available, it may be time consuming to apply the reset input. Therefore there are cases where the use of a reset input is not preferred [START_REF] Hierons | Generating a checking sequence with a minimum number of reset transitions[END_REF][START_REF] Schrammel | Chaining test cases for reactive system testing[END_REF][START_REF] Groz | Inferring finite state machines without reset using state identification sequences[END_REF].

A synchronizing sequence 4 for an FSM M is a sequence of inputs such that no matter at which state M currently is, if this sequence of inputs is applied, M is brought to a particular state. Therefore a synchronizing sequence is in fact a compound reset input, and can be used as such to simulate a reset input in the context of FSM based testing [START_REF] Jourdan | Reduced checking sequences using unreliable reset[END_REF].

A synchronizing sequence may not exist for an FSM. However, as the size of the FSM gets larger, there almost always exists a synchronizing sequence [START_REF] Berlinkov | On the probability of being synchronizable[END_REF].

For an FSM M with n states and alphabet size p, checking if M has a synchronizing sequence can be decided in time O(pn 2 ) [START_REF] Eppstein | Reset sequences for monotonic automata[END_REF]. Since a synchronizing sequence will possibly be used many times in a test sequence, computing a shortest one for an FSM is of interest, but this problem is known to be NP-hard [START_REF] Eppstein | Reset sequences for monotonic automata[END_REF]. There exist a number of heuristics, called synchronizing heuristics, to compute short synchronizing sequences, such as Greedy [START_REF] Eppstein | Reset sequences for monotonic automata[END_REF] and Cycle [START_REF] Trahtman | Some results of implemented algorithms of synchronization[END_REF] both with time complexity O(n 3 +pn 2 ), SynchroP and SynchroPL [START_REF] Roman | Synchronizing finite automata with short reset words[END_REF] with time complexity O(n 5 + pn 2 ), and FastSynchro [START_REF] Kudlacik | Effective synchronizing algorithms[END_REF] with time complexity O(pn 4 ). The upper bound for the length of the synchronizing sequence that will be produced by all of these heuristics is O(n 3 ). Although synchronizing sequences are important for testing methods, the scalability of the synchronizing heuristics has not been addressed thoroughly. For practical applications, the use of even the fastest algorithms (Greedy and Cycle) with cubic complexity can be a problem.

In this work we investigate the use of modern multicore CPUs to scale the performance of synchronizing heuristics. We consider the Greedy algorithm to start with, as it is one of the two cheapest synchronizing heuristics (in practice as well [START_REF] Roman | Forward and backward synchronizing algorithms[END_REF]), known to produce shorter sequences than Cycle [START_REF] Roman | Forward and backward synchronizing algorithms[END_REF], and has been widely used as a baseline to evaluate the quality and speed of more advanced heuristics. To the best of our knowledge, this is the first work towards parallelization of synchronizing heuristics. Although, a parallel approach for constructing a synchronizing sequence for a partial machines is proposed in [START_REF] Türker | Parallel algorithm for deriving reset sequences from deterministic incomplete finite automata -submitted[END_REF], the method proposed in [START_REF] Türker | Parallel algorithm for deriving reset sequences from deterministic incomplete finite automata -submitted[END_REF] is not exact (in the sense that it may fail to find a synchronizing sequence even if one exists) and also it is not a polynomial time algorithm.

All synchronizing heuristics consist of a preprocessing phase, followed by synchronizing sequence generation phase. As presented in this paper, our initial experiments revealed that the preprocessing phase dominates the runtime of the overall algorithm for Greedy. Therefore for both parallelization and for algorithmic improvements of Greedy, we mainly focus on the first phase of the algorithm. With no parallelization, our algorithmic improvements alone yield a 20x speedup on Greedy for automata with 4000 states and 128 inputs. Furthermore, around 150x speedup has been obtained for the same class of automata, when the improved algorithm is executed in parallel with 16 threads.

The rest of the paper is organized as follows: In Section 2, the notation is given, and synchronizing sequences are formally defined. We give the details of Eppstein's Greedy construction algorithm in Section 3. The proposed improvements and the parallelization approach together with implementation details are described in Section 4. Section 5 presents the experimental results and Section 6 concludes the paper.

Preliminaries

FSMs are used to describe a reactive behaviour, i.e., when an input is applied to an FSM, it produces an output as a response. However, the output sequence produced by the application of a synchronizing sequence does not play a role. Therefore, in the context of synchronizing sequences, an FSM can simply be considered as an automaton where the state transitions are only performed by the application of an input, and no output is produced.

In this work, we only consider complete deterministic automata. An automaton is defined by a triple A = (S, Σ, δ) where S is a finite set of n states, Σ is a finite set of p input symbols (or simply inputs) called the alphabet. δ : S ×Σ → S is a transition function. If the automaton A is at a state s and if an input x is applied, then A moves to the state δ(s, x). Figure 1 shows an example automaton A with 4 states and 2 inputs. An element of the set Σ is called an input sequence. We use |w| to denote the length of w, and ε is the empty input sequence. We extend the transition function δ to a set of states and to an input sequence in the usual way. We have δ(s, ε) = s, and for an input sequence w ∈ Σ and an input symbol x ∈ Σ, we have δ(s, xw) = δ(δ(s, x), w). For a set of states S ⊆ S, we have δ(S , w) = {δ(s, w)|s ∈ S }.

We use the notation δ -1 (s, x) to denote the set of those states with a transition to state s with input x. Formally, δ -1 (s, x) = {s ∈ S|δ(s , x) = s}.

Let A = (S, Σ, δ) be an automaton, and w ∈ Σ be an input sequence. w is said to be a merging sequence for a set of states S ⊆ S if |δ(S , w)| = 1, and S is called mergable. Any set {s} with a single state is mergable, since ε is a merging sequence for {s}. w is called a synchronizing sequence for A if |δ(S, w)| = 1. A is called synchronizable if there exists a synchronizing sequence for A. For example, the automaton given in Figure 1 is synchronizable, since abbbabbba is a synchronizing sequence for the automaton. Deciding if an automaton is synchronizable or not can be performed in polynomial time based on the following result.

Proposition 1 ( [START_REF] Eppstein | Reset sequences for monotonic automata[END_REF][START_REF] Natarajan | An algorithmic approach to the automated design of parts orienters[END_REF]). An automaton A = (S, Σ, δ) is synchronizable iff for all s i , s j ∈ S, there exists a merging sequence for {s i , s j }.

For a set of states C ⊆ S, let C 2 = { s i , s j |s i , s j ∈ C} be the set of all multisets with cardinality 2 with elements from C, i.e. C 2 is the set of all subsets of C with cardinality 2, where repetition is allowed. An element s i , s j ∈ C 2 is called a singleton if s i = s j , otherwise it is called a pair.

As Proposition 1 makes it explicit, checking the existence of merging sequences for pairs of states is needed to decide if an automaton is synchronizable.

In addition, the heuristic algorithms also make use of the merging sequences for pairs. For both checking the existence of merging sequences and finding a merging sequence (in fact for finding a shortest merging sequence) for pairs of states of an automaton, one can use the notion of the pair automaton, which we define next. Definition 1. For an automaton A = (S, Σ, δ), the pair automaton A of A is defined as A = (S 2 , Σ, ∆), where for a state s i , s j ∈ S 2 and an input symbol x ∈ Σ, ∆( s i , s j , x) = δ(s i , x), δ(s j , x) .

Eppstein's Algorithm

In this section, we explain Eppstein's Greedy algorithm, and we present an observation on the timing profile of the algorithm. This observation guided our work on the improvements and parallelization of the algorithm, which will be explained in Section 4. Greedy (and also all other synchronizing heuristics mentioned in Section 1) has two phases. In the first phase, a shortest merging sequence for each mergable pair of states is found. If all pairs are mergable, these merging sequences are used to construct a synchronizing sequence in the second phase.

For a pair of states s i , s j of an automaton A = (S, Σ, δ), checking the existence of a merging sequence for {s i , s j }, and computing a shortest merging sequence for {s i , s j } can be performed in time O(pn 2 ) by finding a shortest path from the state s i , s j of the pair automaton A to a singleton state in A using Breadth First Search (BFS). Since we will have to check the existence and find merging sequences for all pairs of states, one can instead use a backward BFS, seeded at singleton states of the pair automaton, as explained below.

For an automaton A = (S, Σ, δ), a function τ : S 2 → Σ , is called a pairwise merging function (PMF) for A, if for all s i , s j ∈ S 2 , τ ( s i , s j ) is a shortest merging sequence for {s i , s j } if {s i , s j } is mergable, and τ ( s i , s j ) is undefined if {s i , s j } is not mergable. Note that PMF for an automaton A is not unique, and it is a total function iff A is synchronizable. Algorithm 1 computes such a PMF τ for a given automaton A, where initially τ ( s, s ) = ε for the singleton states in S 2 (line 1), and τ ( s i , s j ) is considered to be "undefined" for pair states in S 2 (line 2). The algorithm iteratively computes the values of τ (.) as it discovers shortest merging sequences for more pairs in S 2 .

Algorithm 1 keeps track of a frontier set F which is initialized to all singleton states at line 3. Throughout the algorithm, R represents the remaining set of pairs with τ ( s i , s j ) still being undefined. In each iteration of the algorithm (lines 5-6), a BFS step is performed by using BFS step F2R given in Algorithm 2. BFS step F2R constructs the next frontier F from the current frontier F , by considering each s i , s j ∈ F (line 2). Lines 4-5 of BFS step F2R identify a pair s i , s j ∈ R such that s i = δ(s i , x) and s j = δ(s j , x) for some x ∈ Σ, and lines 6-7 performs the necessary updates. Since this algorithm considers, in a sense, the reverse transitions of s i , s j in the frontier F to reach to pairs s i , s j in R, we call it as "Frontier to Remaining (F2R)" BFS step.

Algorithm 1 eventually assigns a value to τ ( s i , s j ) if {s i , s j } is mergable. Based on Proposition 1, A is synchronizable iff there does not exist a pair state

Algorithm 1: Computing a PMF τ : S 2 → Σ (F2R based) input : An automaton A = (S, Σ, δ) output: A PMF τ : S 2 → Σ foreach singleton s, s ∈ S 2 do τ ( s, s ) = ε; foreach pair si, sj ∈ S 2 do τ ( si, sj ) = undefined; F ←-{ s, s |s ∈ S}; // all singleton states of A R ←-{ si, sj |si, sj ∈ S ∧ si = sj}; // all pair states of A while F is not empty do F, R, τ ←-BFS step F2R(A, F, R, τ );
Algorithm 2: BFS step F2R input : An automaton A = (S, Σ, δ), the frontier F , the remaining set R, τ output: The new frontier F , the new remaining set R , and updated function τ

F ←-∅; foreach si, sj ∈ F do foreach x ∈ Σ do foreach s i , s j such that s i ∈ δ -1 (si, x) and s j ∈ δ -1 (sj, x) do if τ ( s i , s j ) is undefined then // s i , s j ∈ R τ ( s i , s j ) ←-xτ ( si, sj ); F = F ∪ { s i , s j }; let R be R \ F ;
s i , s j with τ ( s i , s j ) being undefined when Algorithm 1 terminates. We can now present Eppstein's Greedy algorithm based on Algorithm 1.

The Greedy algorithm keeps track of a current set C of states yet to be merged, initialized to S at line 4. A pair s i , s j ∈ C 2 is called an active pair.

In each iteration of the while loop at line 7, an active pair s i , s j ∈ C 2 is found such that it has a shortest merging sequence among all active pairs in C Algorithm 3: Eppstein's Greedy Algorithm input : An automaton A = (S, Σ, δ) output: A synchronizing sequence Γ for A (or fail if A is not synchronizable) compute a PMF τ using Algorithm 1; if there exists a pair si, sj such that τ ( si, sj ) is undefined then report that A is not synchronizable and exit; foreach si, sj, s k ∈ S do compute δ(s k , τ ( si, sj )); C = S; // C will keep track of the current set of states Γ = ε; // Γ is the synchronizing sequence to be constructed while |C| > 1 do // we have two or more states yet to be merged (line 8). The synchronizing sequence (initialized to the empty sequence at line 6) is extended with τ ( s i , s j ) at line 9. Finally, τ ( s i , s j ) is applied to C to update the current set of states. When |C| = 1, this means that Γ accumulated at that point is a synchronizing sequence.

find a pair si, sj ∈ C 2 with minimum |τ ( si, sj )| among all pairs in C 2 ; Γ = Γ τ ( si, sj ); 10 C = δ(C, τ ( si, sj )); n = 1000 n = 2000 n = 4000 p tALL tP M F t P M F t ALL tALL tP M F t P M F t ALL tALL tP M F t P M F t ALL 2 0,
The following results are shown in [START_REF] Eppstein | Reset sequences for monotonic automata[END_REF]Theorem 5]. For an automaton A with n states and p inputs, Phase 1 of Greedy (lines 1-3) can be implemented to run in time O(pn 2 ) and Phase 2 of Greedy (lines 4-10) can be implemented to run in time O(n 3 ). Hence the overall time for Greedy is O(n 3 + pn 2 ).

We performed an experimental analysis to see how much Phase 1 (which we will call as the PMF construction phase5 ) and Phase 2 (the synchronizing sequence construction phase) of the algorithm contribute to the running time in practice for a sequential implementation. Based on these experiments, we observed that PMF construction actually dominates the running time of the algorithm (see Table 1). Hence, in order to improve the performance of Greedy, we developed approaches for parallel implementation of PMF construction, together with some algorithmic modifications, which we explain in Section 4.

Parallelization Approach and Improvements

Algorithm 1 necessarily performs a BFS on the pair automaton A, and a BFS forest rooted at singleton states of A is implicitly obtained. At the roots of the forest (i.e. in the first frontier set F ) we have singleton states of A, which corresponds to the nodes at level 0 of the BFS forest. At each iteration of the algorithm, the current frontier F has all the nodes at level k in the BFS forest. These nodes are processed by Algorithm 2 to compute the next frontier F which are the nodes at level k + 1 in the BFS forest. The processing of the state pairs in F are the tasks to be performed at the current level. To process a state pair, Algorithm 2 considers incoming transitions of the pair (i.e., inverse transitions) based on the δ -1 function (line 4). Hence, the cost of each task can be different. Furthermore, the total number of edges of the tasks in F , i.e., frontier edges, determines the cost of the corresponding level's BFS step F2R execution and this also varies for each level. We used OpenMP for parallel implementation and employed the dynamic scheduling policy (with batches of 512-pairs) since the task costs are not uniform.

Computing a PMF in parallel

When Algorithm 1 is implemented sequentially, handling two consecutive iterations is seamless: using a single queue to enque and deque the frontier pairs suffices to process them in the correct order (i.e. a pair at level k +1 is only found after all level k pairs are found). However, with multiple threads, a barrier (a global synchronization technique) is required after each iteration. Otherwise, a pair from the next frontier can be processed before another pair in the current frontier and an incorrect PMF function τ can be computed. Here we present Algorithm 1 iteratively, and isolate the BFS step F2R from the main flow of the algorithm since it will be our main target for efficiency.

Algorithm 4: BFS step F2R (in parallel)

input : An automaton A = (S, Σ, δ), the frontier F , the remaining set R, τ output: The new frontier F , the new remaining set R , and updated function τ

1 foreach thread t do F t ←-∅ ; 2 foreach si, sj ∈ F in parallel do 3 foreach x ∈ Σ do 4 foreach s i , s j where s i ∈ δ -1 (si, x) and s j ∈ δ -1 (sj, x) do 5 if τ ( s i , s j ) is undefined then // s i , s j ∈ R 6 τ ( s i , s j ) ←-xτ ( si, sj ); 7 F t = F t ∪ { s i , s j }; 8 F ←-∅; 9 foreach thread t do F = F ∪ F t ; 10 let R be R \ F ;
To parallelize BFS step F2R, we partition the current frontier F among multiple threads where only a single thread processes a frontier pair as shown in Algorithm 4 (line 2). Since there is no task-dependency among the pairs, all the threads can simultaneously work. However, a race condition occurs since the next frontier set F is a shared object in the sequential implementation. To break dependency with a lock-free approach, in our parallel implementation, each thread t uses a local frontier array F t and when a new pair from the next frontier is found by thread t, it is immediately added to F t . When two threads find the same pair s i , s j at the same time, both threads insert it to their local frontiers (lines 5-7). Hence, when the local frontiers are combined at the end of each iteration (lines 8-9), the same pair can occur multiple times if no duplicate pair check is applied. In our preliminary experiments, we observed that at most one in a thousands extra pairs are inserted to F when they are allowed. Hence, we let the threads process them since the total extra pair cost is negligible compared to the cost of checking and resolving duplicates.

Another approach for BFS steps

Algorithm 2 and Algorithhm 4 follow a natural and possibly the most common technique to construct the next frontier set F from the current frontier set F by considering the incoming transitions. Another approach to construct the next frontier F function, which we call "Remaining to Frontier (R2F)", is processing the remaining state pairs' edges instead of those in the frontier. As mentioned above, a state pair s i , s j stays in R, i.e., in the remaining pair set, as long as τ ( s i , s j ) stays undefined. In the parallel R2F approach described by Algorithm 5, the threads process the transitions of the remaining state pairs instead of the ones in the frontier. Hence, instead of δ -1 , the original transition function δ is used and the pair found is checked to be in the frontier (lines 5-6). If a pair s i , s j has a transition to a pair s i , s j ∈ F (i.e., if s i , s j is in the next frontier), τ ( s i , s j ) is set and the process ends (lines 7-9). Otherwise, s i , s j is kept in the remaining set (lines [START_REF] Jourdan | Reduced checking sequences using unreliable reset[END_REF][START_REF] Berlinkov | On the probability of being synchronizable[END_REF]. Similar to parallel F2R, we use a local remaining pair array R t for each thread t in the lock-free parallelization of R2F.

Algorithm 5: BFS step R2F (in parallel)

input : An automaton A = (S, Σ, δ), the frontier F , the remaining set R, τ output: The new frontier F , the new remaining set R , and updated function τ

1 foreach thread t do R t ←-∅; 2 foreach si, sj ∈ R in parallel do 3 connected ←-false; 4 foreach x ∈ Σ do 5 s i , s j ←-δ(si, x), δ(sj, x) ; 6 if τ ( s i , s j ) is defined then // s i , s j ∈ F 7 τ ( si, sj ) ←-xτ ( s i , s j ); 8 connected ←-true; 9 break; 10 if not connected then 11 R t = R t ∪ { si, sj }; 12 R ←-∅; 13 foreach thread t do R = R ∪ R t ; 14 let F be R \ R ;

A hybrid approach to construct the next frontier

Since the size of R decreases at each iteration, R2F becomes faster at each step. On the other hand, F2R is expected to be faster than R2F during the earlier iterations. Therefore it makes sense to use a hybrid approach, where either an F2R or an R2F BFS step is used depending on their respective cost for the current iteration. These observations have been used by Beamer et al. to implement a direction-optimized BFS [START_REF] Beamer | Direction-optimizing breadth-first search[END_REF]. Since the cost of each F2R/R2F Algorithm 6: Computing a function τ : S 2 → Σ (Hybrid) input : An automaton A = (S, Σ, δ) output: A function τ : S 2 → Σ 1 foreach singleton s, s ∈ S 2 do τ ( s, s ) = ε; 2 foreach pair si, sj ∈ S 2 do τ ( si, sj ) = undefined; 3 F ←-{ s, s |s ∈ S}; // all singleton states of A 4 R ←-{ si, sj |si, sj ∈ S ∧ si = sj}; // all pair states of A 5 while F is not empty do

6 if |F | < |R| then 7 F, R, τ ←-BFS step F2R(A, F, R, τ ); 8 else 9 F, R, τ ←-BFS step R2F(A, F, R, τ );
iteration depends on the number of edges processed, it is reasonable to compare the number of frontier/remaining pairs' edges to choose the cheaper approach at each iteration as in [START_REF] Beamer | Direction-optimizing breadth-first search[END_REF]. When the BFS is executed on a simple graph, this strategy is easy to apply. However, by only using δ -1 , it takes O(p) time to count a new frontier pair's edges. Overall, the counting process takes O(pn 2 ) time which is expensive considering that the overall sequential complexity is also O(pn 2 ). In this work, we compared the size of R and F instead of the edges to be processed. The total additional complexity due to counting is O(n 2 ) since each pair will be counted only once.

To analyze the validity of our counting heuristic and the potential improvement due to the Hybrid approach described in Algorithm 6, we compared the size of R and F , and the corresponding execution time of each F2R/R2F execution in Figure 2. As the figure shows, counting the pairs instead of transitions can be a good heuristic to guess the cheaper approach in our case. Furthermore, the performance difference of F2R and R2F at the each iteration shows that the proposed Hybrid approach can yield a much better performance.

Implementation details

To store and utilize the δ -1 (s, x) for all x ∈ Σ and s ∈ S, we employ the data structures in Fig. 1 (right). For each symbol x ∈ Σ, we used two arrays ptrs x and js x where the former is of size n+ 1 and the latter is of size n. For each state s ∈ S, ptrs x [s] and ptrs x [s + 1] are the start (inclusive) and end (exclusive) pointers to two js x entries. The array js x stores the ids of the states δ -1 (s, x) in between js x [ptrs x [s]] and js x [ptrs x [s + 1] -1]. This representation has a low memory footprint. Furthermore, we access the entries in the order of their array placement in our implementation hence, it is also good for spatial locality.

The memory complexity of the algorithms investigated in this study is O(n 2 ). For each pair of states, we need to employ an array to store the length of the shortest merging sequence. To do that one can allocate an array of size n 2 , Fig. 3 (left), and given the array index = (i -1) × n + j for a state pair {s i , s j } where 1 ≤ i ≤ j ≤ n, she can obtain the state ids by i = n and j = -((i -1) × n). This simple approach effectively uses only the half of the array since for a state pair {s i , s j }, a redundant entry for {s j , s i } is also stored. In our implementation, Fig. 3 (right), we do not use redundant locations. For an index = i×(i+1) 2 + j the state ids can be obtained by i = √ 1 + 2 -0.5 and j = -i×(i+1)

2

. Preliminary experiments show that this approach, which does not suffer from the redundancy, also have a positive impact on the execution time. That being said, all the algorithms in the paper uses it and this improvement will not have change their relative performance. Fig. 3. Indexing and placement of the state pair arrays. A simple placement of the pairs (on the left) uses redundant places for state pairs {si, sj}, i = j, e.g., {s1, s2} and {s2, s1} in the figure. On the right, the indexing mechanism we used is shown.

Experimental Results

All the experiments in the paper are performed on a single machine running on 64 bit CentOS 6.5 equipped with 64GB RAM and a dual-socket Intel Xeon E7-4870 v2 clocked at 2.30 GHz where each socket has 15 cores (30 in total). For the multicore implementations, we used OpenMP and all the codes are compiled with gcc 4.9.2 with the -O3 optimization flag enabled.

To measure the efficiency of the proposed algorithms, we used randomly generated automatons 6 with n ∈ {1000, 2000, 4000} states and p ∈ {2, 8, 32, 128} inputs. For each (n, p) pair, we randomly generated 20 different automatons and executed each algorithm on these automatons. The values in the figures and the tables are the averages of these 20 executions for each configuration, i.e., algorithm, n and p.

Multicore parallelization of PMF construction

Figure 4 shows the speedups of our parallel F2R implementation over the sequential baseline (that has no parallelism). Since F2R uses the same frontier extension mechanism with the sequential baseline, and R2F employs a completely different one, here we only present the speedup values of F2R. As the figure shows, when p is large, the parallel F2R presents good speedups, e.g., for p = 128, the average speedup is 14.1 with 16 threads. Furthermore, when compared to the single-thread F2R, the average speedup is 15.2 with 16 threads. A performance difference between sequential baseline and single-threaded F2R exists because of the parallelization overhead during the local queue management. Overall, we observed 10% parallelization penalty for F2R on the average over the sequential baseline for all (n, p) pairs. For p values smaller than 128, i.e., 2, 8, and 32, the average speedups are 5.4, 9.1, and 12.8, respectively, with 16 threads. The impact of the parallelization overhead is more for such cases since the amount of the local-queue overhead is proportional to the number of states but not to the number of edges. Consequently, when p decreases the amount of total work decreases and hence, the impact of the overhead increases. Furthermore, since the number of iterations for PMF construction increases with decreasing p, the local queues are merged more for smaller p values. Therefore, one can expect more overhead, and hence, less efficiency for smaller p values as the experiments confirm. Figure 5 compares the execution times of F2R, R2F and Hybrid algorithm for n = 1000 (top) and n = 4000 (bottom) states, p ∈ {2, 8, 32} and {1, 2, 4, 8, 16} threads (the results for n = 2000 are similar but omitted due to space limitations). For better figure scaling, the results for p = 128 is given in Figure 6. An interesting observation is that F2R is consistently faster than R2F for p = 2, however, it is slower otherwise. This can be explained by the difference in the number of required iterations to construct PMF: when p is large, the frontier expands very quickly and the PMF is constructed in less iterations, e.g., for n = 2000, the PMF is generated in 16 iterations for p = 2, whereas only 7 iterations are required for p = 8. Since each edge will be processed once, the runtime of F2R always increases with p, i.e., with the number of edges. However, since the frontier expands much faster, the total number of remaining (R-)pairs processed by the R2F throughout the process will probably decrease. Furthermore, since when the frontier is large, while traversing the edge list of an R-pair, it is more probable to early terminate the traversal and add the R-pair to the next frontier earlier. Surprisingly, when p increases, these may yield a decrease in the R2F runtime (observe the change from p = 2 to p = 8 in Fig. 5). However, once the performance benefits of early termination are fully exploited, an increase on the R2F runtime with increasing p is more probable since the overall BFS work, i.e., the total number of edges, also increases with p (observe the change from p = 8 to p = 32 in Fig. 5).

Observing such performance differences for R2F and F2R on automatons with different characteristics, the potential benefit of a Hybrid algorithm in practice is more clear. As Figure 5 and Figure 6 show, the hybrid approach, which is just a combination of F2R and R2F, is almost always faster than employing a pure F2R or a pure R2F BFS-level expansion. Furthermore, we do not need parallelism to observe these performance benefits: the Hybrid approach works better even when a single thread is used at runtime. For example, when n = 4000 and p = 128, the Hybrid algorithm is 23 and 6 times faster than F2R and R2F, respectively. For the same automaton set, the speedups due to hybridization of the process become 14 and 4 with 16 threads on average.

When the Hybrid algorithm is used, the speedups on the PMF generation phase are given in Figure 7. As the figure shows, thanks to parallelism and good scaling of Hybrid (for large p values), the speedups increase when the number of threads increases. The PMF generation process becomes 95, 165, and 199 times faster when 16 threads used for 1000, 2000, and 4000 state automatons, respectively. Even with single thread, i.e., no parallelization, the Hybrid heuristic is 8, 14, and 21 times faster than the sequential algorithm. Since we generate the PMF to find a synchronizing sequence, a more practical evaluation metric would be the performance improvement over the sequential reset sequence construction process. As Table 1 shows, for Eppstein's Greedy heuristic (also for some other heuristics such as Cycle [START_REF] Trahtman | Some results of implemented algorithms of synchronization[END_REF]), the PMF generation phase dominates the overall runtime. For this reason, we simply conducted an experiment where the Hybrid approach is used to construct the PMF and no further parallelization is applied during the synchronizing sequence construction phase. Table 2 shows the speedups for this experiment for single thread and 16 thread Hybrid executions. As the results show, even when the sequence construction phase is not parallelized, more than 50x and more than 100x improvement is possible for p = 32 and p = 128, respectively. As noted before, F2R based PMF construction has O(pn 2 ) time complexity. R2F based PMF construction, on the other hand, has O(dpn 2 ) time complexity (where d is the diameter of the pair automaton A), since states of A in the remaining set R will be processed at most d times. In practice, however, R2F based construction (and Hybrid computation which also has O(dpn 2 ) time complexity since it performs R2F steps) can beat F2R based construction.

We did not perform an extensive study on automata with larger state numbers, since it takes too long with the sequential baseline implementation. For example, sequential PMF generation takes around 75 minutes for an automaton with 32000 states and 128 letters, whereas our Hybrid implementation completes in 3 minutes.

Conclusion and Future Work

We investigated the efficient implementation and use of modern multicore CPUs to scale the performance of synchronizing sequence generation heuristics. We parallelized one of the well-known heuristics Greedy. We mainly focused on the PMF generation phase (which is employed by almost all the heuristics in the literature), since it is the most time consuming part of Greedy. Even with no parallelization, our algorithmic improvements yielded a 20x speedup on Greedy for automatons with 4000 states and 128 inputs. Furthermore, around 150x speedup has been obtained with 16 threads for the same automata class.

To eliminate threats to validity, we checked and confirmed that the sequence constructed by each algorithm is indeed a synchronizing sequence. We also compared the length of the sequences constructed by the original implementation of Greedy and different versions of Greedy algorithms suggested in this paper. We observed that regardless of the PMF construction approach used, for each pair s i , s j , we obtain the same length |τ ( s i , s j )| for the shortest merging sequences, but the actual shortest merging sequence τ ( s i , s j ) can differ, which causes around ±1% difference in the length of the synchronizing sequences.

As a future work, we will apply our techniques to other heuristics in the literature that are relatively slower than Greedy but can produce shorter synchronizing sequences. For these heuristics, parallelizing only the PMF generation phase may not be sufficient since the synchronizing sequence construction part of these heuristics are much more expensive compared to Greedy. Hence, we aim to parallelize the whole sequence generation process. Another problem we want to study is the use of cutting-edge manycore architectures such as GPUs and FP-GAs to make such heuristics faster and more practical for large scale automatons.

Fig. 1 .

 1 Fig.1. A synchronizable automaton A (left), and the data structures we used to store and process the transition function δ -1 in memory (see Section 4.4 for the details). A synchronizing sequence for A is abbbabbba.

Fig. 2 .

 2 Fig. 2. The number of frontier and remaining vertices at each BFS level and the corresponding execution times of F2R and R2F while constructing the PMF τ for n = 2000 and p = 8 (top) and p = 128 (bottom).

Fig. 4 .

 4 Fig. 4. The speedup of our parallel F2R PMF construction over the sequential PMF construction baseline.

Fig. 5 .

 5 Fig. 5. Comparison of the parallel execution times of the three PMF construction algorithms: (1) F2R, (2) R2F, and (3) hybrid. The figures show the times for n = 1000 (top) and n = 4000 (bottom), p ∈ {2, 8, 32}, with {1, 2, 4, 8, 16} threads (x-axis). For a better readability and figure scaling, the single-thread F2R bars with 32 inputs are allowed to exceed the max value on the y-axis.

Fig. 6 .

 6 Fig. 6. Comparison of the parallel execution times of the three PMF construction algorithms: (1) F2R, (2) R2F, and (3) hybrid. The figures show the times for n = 1000 (left), n = 2000 (middle), and n = 4000 (bottom), p = 128, with {1, 2, 4, 8, 16} threads (x-axis).

Fig. 7 .

 7 Fig. 7. The speedups of the Hybrid PMF construction algorithm with n = 1000 (top), 2000 (middle), 4000 (bottom) and p ∈ {2, 8, 32, 128}. The x-axis shows the number of threads used for the Hybrid execution. The values are computed based on the average sequential PMF construction time over 20 different automatons for each (n, p) pair.

Table 1 .

 1 Sequential PMF construction time (tP M F ), and overall time (tALL) for automata with n ∈ {1000, 2000, 4000} states and p ∈ {2, 8, 32, 128} inputs.

		045	0,042	0,928	0,188	0,175	0,929	1,214	1,158	0,954
	8	0,125	0,122	0,974	0,526	0,513	0,975	2,757	2,698	0,979
	32	0,483	0,480	0,993	2,151	2,138	0,994	9,980	9,919	0,994
	128	2,202	2,199	0,999	9,243	9,229	0,999	39,810	39,749	0,998

Table 2 .

 2 The speedups obtained on Eppstein's Greedy algorithm when the Hybrid PMF construction algorithm is used.

			p (single thread)			p (16 threads)
	n	2	8	32	128	2	8	32	128
	1000	1,2	1,8	13,4	7,5	4,6	10,8	58,2	83,7
	2000	1,2	2,7	2,2	14,0	4,8	13,1	24,3	133,9
	4000	1,1	2,9	3,3	20,7	5,5	14,8	31,7	154,0

Synchronizing sequences are also known as reset sequences, or reset words.

Lines 2-3 of Phase 1 is easily handled as a part of PMF construction by checking if R is empty or not at the end of PMF construction.

For each state s and input x, δ(s, x) is randomly assigned to a state s ∈ S.

Acknowledgements

This work is supported by T ÜB İTAK Grants #114E569 and #115C018.


Eduard P Enoiu

Daniel Sundmark

Adnan Čaušević

Robert Feldt

Paul Pettersson

Mutation-Based Test Generation for PLC Embedded Software using Model Checking

Testing is an important activity in engineering of industrial embedded software. In certain application domains (e.g., railway industry) engineering software is certified according to safety standards that require extensive software testing procedures to be applied for the development of reliable systems. Mutation analysis is a technique for creating faulty versions of a software for the purpose of examining the fault detection ability of a test suite. Mutation analysis has been used for evaluating existing test suites, but also for generating test suites that detect injected faults (i.e., mutation testing). To support developers in software testing, we propose a technique for producing test cases using an automated test generation approach that operates using mutation testing for software written in IEC 61131-3 language, a programming standard for safety-critical embedded software, commonly used for Programmable Logic Controllers (PLCs). This approach uses the Uppaal model checker and is based on a combined model that contains all the mutants and the original program. We applied this approach in a tool for testing industrial PLC programs and evaluated it in terms of cost and fault detection. For realistic validation we collected industrial experimental evidence on how mutation testing compares with manual testing as well as automated decision-coverage adequate test generation. In the evaluation, we used manually seeded faults provided by four industrial engineers. The results show that even if mutation-based test generation achieves better fault detection than automated decision coverage-based test generation, these mutation-adequate test suites are not better at detecting faults than manual test suites. However, the mutation-based test suites are significantly less costly to create, in terms of testing time, than manually created test suites. Our results suggest that the fault detection scores could be improved by considering some new and improved mutation operators (e.g., Feedback Loop Insertion Operator (FIO)) for PLC programs as well as higher-order mutations.

Introduction

Software testing is an important verification and validation activity used to reveal software faults and make sure that actual software behavior matches its expected behavior [START_REF] Ammann | Introduction to Software Testing[END_REF]. Safety-critical and real-time software systems implemented in Programmable Logic Controllers (PLCs) are used in many real-world industrial application domains. One of the programming languages defined by the International Electrotechnical Commission (IEC) for PLCs is the Function Block Diagram (FBD) language. In testing IEC 61131-3 FBD programs in the railway domain, the engineering processes of software development are performed according to safety standards and regulations [START_REF]Railway Application: Communications, Signaling and Processing Systems, Software For Railway Control and Protection Systems[END_REF]. As an alternative to manually testing software, a few techniques for automated test generation have been proposed [START_REF] Fraser | Evosuite: Automatic Test Suite Generation for Objectoriented Software[END_REF][START_REF] Cadar | KLEE: Unassisted and Automatic Generation of High-Coverage Tests for Complex Systems Programs[END_REF]. While high code coverage has historically been used as a proxy for the ability of a test suite to detect faults, recent results (e.g., [START_REF] Inozemtseva | Coverage is Not Strongly Correlated with Test Suite Effectiveness[END_REF]) indicate that code coverage may not be a good measure of fault detection effectiveness. As an alternative to coverage-based test generation, mutation testing has been proposed [START_REF] Demillo | Constraint-based automatic test data generation[END_REF][START_REF] Fraser | Mutation-driven generation of unit tests and oracles[END_REF]. In mutation testing, test cases are generated based on the concept of mutants-small syntactic modifications in the program, intended to imitate real faults. A set of test cases that can distinguish a certain program from its mutants is sensitive to faults, and it thus hypothesized to be good at detecting real faults (a hypothesis that has strong empirical support [START_REF] Just | Are mutants a valid substitute for real faults in software testing?[END_REF]). However, for domain specific languages used in embedded software development (i.e., IEC 61131-3), there is a lack of mature approaches and tools for performing mutation test generation.

In this paper, we describe and evaluate an automated mutation-based test generation approach for IEC 61131-3 embedded software. The main contributions of the paper are:

-An approach for mutation test generation of IEC 61131-3 programs using a model checker by combining all the mutants and the original program into a single combined model that is monitored dynamically. -An evaluation of the approach in an industrial case study. The results show that mutation-adequate test suites are worse at detecting faults than manual test suites with the cost of performing mutation testing being consistently lower than the cost of manually testing IEC 61131-3 software. -The identification of new mutation operators for mutation testing of IEC 61131-3 software. The reduction in fault detection between manual and mutation testing was attributed based on our analysis to an incomplete list of mutation operators for IEC 61131-3 software. We propose new operators simulating this kind of faults (e.g., Feedback loop Insertion Operator (FIO)).

The rest of the paper is organized as follows. Section 2 introduces PLC embedded software, automated test generation and mutation testing. Section 3 describes the approach for mutation test generation for IEC 61131-3 programs using a model checker. Section 4 explains the experimental method, while the results are provided and discussed in Section 5. Finally, Section 6 concludes the paper.

Background and Related Work

This paper describes a method for mutation testing for PLC embedded programs implemented in the IEC 61131-3 FBD language. In this section, we provide a background on PLC embedded software, automated test suite generation and mutation testing.

PLC Embedded Software

Safety-critical embedded systems implemented using Programmable Logic Controllers (PLCs) are used in many industrial application domains such as electric, transportation, chemical, pharmaceutical, etc.. One of the programming languages defined by the International Electrotechnical Commission (IEC) for PLCs is the Function Block Diagram (FBD) language [START_REF]IEC: International Standard on 61131-3 Programming Languages[END_REF]. Programs developed in FBD are compiled into program code, which in turn is compiled into machine code by using specific engineering tools provided by PLC vendors. The motivation for using FBD as the target language in this study comes from the fact that it is the de facto standard in many industrial systems [START_REF] Schwartz | Control System Devices: Architectures and Supply Channels Overview[END_REF], such as the ones in the railway transportation domain. Programs running on a PLC execute in a loop, in which the iteration follows the "read-execute-write" semantics. FBD is popular because of its graphical notations and its usefulness in applications with a high degree of data flow between control components. As shown in Figure 1, predefined logical and/or stateful blocks (i.e., SR, XOR, TOF, LT and TON in Figure 1) and signals (i.e., connections) between blocks represent the behavior of an FBD program. The blocks are supplied by the hardware manufacturer or defined by a developer. PLCs contain particular types of blocks called timers (e.g., TON and TOF) that provide the same functions as timing relays in electrical circuits and are used to activate or deactivate a device after a preset interval of time. For more details on this programming language we refer the reader to the work of John et al. [START_REF] John | IEC 61131-3: Programming Industrial Automation Systems[END_REF].

Automated Test Generation for PLC Embedded Software

In general, automated test generation has been explored in a considerable amount of work [START_REF] Orso | Software testing: a research travelogue[END_REF] in the last couple of years. Numerous techniques for automated test generation using code coverage criteria (e.g., [START_REF] Fraser | Evosuite: Automatic Test Suite Generation for Objectoriented Software[END_REF][START_REF] Cadar | KLEE: Unassisted and Automatic Generation of High-Coverage Tests for Complex Systems Programs[END_REF]) have been proposed in the last decade, since test suites can be created and executed with reduced human effort and cost. However, for domain specific languages used in embedded software development, contributions have been more sparse. For IEC 61131-3 software, a few automated test generation approaches [START_REF] Wu | Automatic Test Case Generation for Structural Testing of FBD[END_REF][START_REF] Jamro | POU-Oriented Unit Testing of IEC 61131-3 Control Software[END_REF][START_REF] Simon | Hauck-Stattelmann: Automatic Test Case Generation for PLC Programs Using Coverage Metrics[END_REF] have been proposed in the last couple of years, but currently there is a lack of tool support. In our previous work, we developed an automated test input generation approach and tool named CompleteTest [START_REF] Enoiu | Automated Test Generation using Model Checking: an Industrial Evaluation[END_REF], which automatically produces test suites for a given coverage criterion and an IEC 61131-3 program written using the FBD language. CompleteTest supports different code coverage criteria with the default criterion being decision coverage.

Mutation Testing

Recent work [START_REF] Gay | The Risks of Coverage-Directed Test Case Generation[END_REF][START_REF] Inozemtseva | Coverage is Not Strongly Correlated with Test Suite Effectiveness[END_REF] suggests that coverage criteria alone can be a poor indication of fault detection in testing. To tackle this issue, researchers have proposed approaches for improving fault detection by using mutation analysis as a test criterion. Mutation analysis is the technique of automatically generating faulty implementations of a program for the purpose of examining the fault detection ability of a test suite [START_REF] Demillo | Hints on Test Data Selection: Help for the Practicing Programmer[END_REF]. A mutant is a new version of a program created by making a small change to the original program. The execution of a test case on the resulting mutant may produce a different output as the original program, in which case we say that the test case kills that mutant. The mutation score is calculated using either an output-only oracle (i.e., strong mutation [START_REF] Woodward | From weak to strong, dead or alive? an analysis of some mutation testing issues[END_REF]) or a state change oracle (i.e., weak mutation [START_REF] Howden | Weak mutation testing and completeness of test sets[END_REF]) against the set of mutants. For all programs, one needs to assess the fault-finding effectiveness of each test suite by calculating the ratio of mutants killed to total number of mutants. When this technique is used to generate test suites rather than evaluating existing ones, it is commonly referred to as mutation testing or mutation-based test generation. Despite its effectiveness [START_REF] Just | Are mutants a valid substitute for real faults in software testing?[END_REF], to the best of our knowledge, no attempt has been made to propose and evaluate mutation testing for PLC embedded software written in the IEC 61131-3 FBD programming language. This motivated us to develop an automated test generation approach based on mutation testing targeting this type of software.

Mutation Test Generation for PLC Embedded Software

Within the last decade model-checking has turned out to be a useful technique for generation of test cases from models [START_REF] Fraser | Testing with Model Checkers: a Survey[END_REF]. In this paper, we describe an approach to automatically generate test suites using a model checker based on mutation testing for PLC embedded software. Overall, the approach is composed of the following steps, mirrored in Figure 2: 1. Mutant Generation. This first step (described in detail in Section 3. The output of this step is a model containing the original structure and behavior of the program together with all inserted mutants. 3. Mutant Annotation. The third step (described in Section 3.3) involves the annotation of the combined model with instrumentation instructions for the detection of each mutant. This means that the mutation detection monitor is used to record the mutant execution and detection, thus for all mutants a property is created for checking the detection of mutants. 4. Test Suite Generation. The fourth step (described in Section 3.4) requires the use of the Uppaal model checker [START_REF] Larsen | Uppaal in a nutshell[END_REF] to generate a set of test cases satisfying the detection of mutants by using the model checker's ability to export abstract traces witnessing a submitted property.

Mutation Generation

To facilitate mutation testing, we begin by generating mutated versions of the original program. The mutation generator parses a given program and processes the structural elements for performing mutations. In particular, for each mutation operator, the program is traversed invoking the corresponding mutation function at all possible locations, each mutation resulting in a separate mutant version of the program. For the creation of mutants, we rely on previous studies that looked at commonly occurring faults in IEC 61131-3 software [START_REF] Oh | Software Safety Analysis of FBD using Fault Trees[END_REF][START_REF] Shin | Empirical Evaluation on FBD Model-based Test Coverage Criteria using Mutation Analysis[END_REF]. We used these common faults in this study for establishing the following mutation operators:

- These mutation operators are systematically applied to the entire program (i.e., blocks, variables, constants, connections) and thus resulting in a set of mutants, each simulating one syntactic change.

Model Aggregation

We start the model aggregation step with the translation of a program and its set of mutants to a timed automata representation. We have shown in a previous study [START_REF] Enoiu | Automated Test Generation using Model Checking: an Industrial Evaluation[END_REF] how the mapping of an IEC 61131-3 program to timed automata is implemented. Timed automata, introduced by Alur and Dill [START_REF] Alur | Model-checking for real-time systems[END_REF], were chosen because there is an already existing formal semantics and tool support for simulation and model-checking using Uppaal [START_REF] Larsen | Uppaal in a nutshell[END_REF] and automated test generation using CompleteTest [START_REF] Enoiu | Automated Test Generation using Model Checking: an Industrial Evaluation[END_REF]. A timed automaton is a standard finite-state automaton extended with time (i.e., real-valued clocks are used for measuring time progress). A model in Uppaal consists of a network of processes that are composed of locations. Transitions between these locations define how the model behaves. The semantics of a timed automaton A is defined in terms of a state transition system, where the state of A is defined as a pair (l, u), where l is a location (i.e. node) and u is a clock assignment. A state of A depends on its current location and on the current values of its clocks. A network of timed automata B 0 ... B n-1 is a parallel composition of n timed automata over synchronization functions (i.e., a! is correlative with a?). Further information on timed automata can be found in [START_REF] Alur | Model-checking for real-time systems[END_REF]. In our previous work [START_REF] Enoiu | Automated Test Generation using Model Checking: an Industrial Evaluation[END_REF] we showed that an IEC 61131-3 FBD program can be transformed to a formal representation containing both its functional and timing behavior. In this study, the model aggregation is using this already developed translation for obtaining the model needed for running mutation-based test generation. Let M be a finite set of of mutants, each of which contains one syntactic change in the original program P . The model aggregation step is applied as follows:

-Create a timed automaton P corresponding to the original FBD -Create a supervision automaton that executes each block and its mutants according to the order of execution. The execution order N is automatically defined according to the general rules included in the IEC 61131-3 standard [START_REF]IEC: International Standard on 61131-3 Programming Languages[END_REF]. This predetermined order directly dictates the data dependency in a program. Basically, each mutated entity executes in parallel with its original counterpart.

As a result of the model aggregation step we consider that the combined model is a closed network of timed automata. This model, briefly shown in Figure 2, contains four processes, two modeling the program and its mutants and the other two supervising the overall execution and monitoring the mutant detection. To show an example of an aggregated model cycle scan, different actions are executed: read(IN) for reading input variables, write(OUT) for updating the output variables, and write(OUT(m i )) for updating the duplicated output variables corresponding to each mutant m i . When the execution order holds, the input variables are updated and the execution continues to the next block.

Mutant Annotation

Informally, our approach is based on the idea that in order to kill all mutants of a specific program, it would be sufficient to (i) annotate the mutants in an FBD program by adding a mutation detection monitor, (ii) formulate a reachability property for the mutation score (i.e., what portion of the existing mutants have been killed), and (iii) find a path from the initial state to some state where the mutation score is 100%. Thus, using auxiliary variables, we annotate the aggregated model such that a condition describing whether a single mutant is killed or not can be expressed.

For annotation, it should be noted that there are different interpretations of how to implement mutation analysis. The most common implementation, called strong mutation deals with the comparison of the original and mutated program outputs at the end of the execution cycle. Another way is weak mutation, which compares the state of the program immediately after the execution of the mutated part of the program. As these implementations can be useful in their different interpretation of mutation analysis, our approach employed both approaches.

Weak Mutation. A mutant is weakly killed in an FBD program if it leads to a block output change (i.e., block infection) compared to the original program behavior. For each mutation operator we define a detection monitor that precisely describes the decision that leads to a change in block output. In model checking we require a reachability property and a mutant detection monitor that guides the search towards detection. We define this weak mutation monitor for individual mutation operators. For each mutant m i in M , where M is the entire set of mutants, there is a weak mutation monitor wm i (M ) that looks at the block output change; if wm i (M ) is 1 then m i is detected. Using a model checker, the aim of weak mutation testing is to achieve a state where all mutants are killed with respect to the block output change. For generating tests for weak mutation we represent the test obligations over a set of variables monitoring the original behavior and its mutants as a reachability property.

Strong Mutation. Weak mutation testing for an FBD program results in a test suite where an internal block is infected; however, a change in block output does not necessarily propagate to an observable program output. Using a model checker, we propose to propagate the mutated behaviors to the output of the program using additional data variables and signals and monitor the change in output using a strong mutation monitor. For each mutant m i in M , where M is the entire set of mutants, the output of each mutant is propagated to the depended blocks until it reaches the program output. There is a strong mutation monitor sm i (M ) that looks at the program output change; if sm i (M ) is 1 then m i is detected. Using a model checker, the aim of strong mutation testing is to achieve a state where all mutants are killed with respect to the program output change. In our scenario, a mutant is killed if there exists a path in the model such that a test input shows that the mutated program output differs from the output of the original program.

Test Generation

In order to generate a test suite for mutation testing of FBD programs using Uppaal, we make use of Uppaal's ability to generate traces witnessing a submitted reachability property. A trace produced by the model checker for a given reachability property defines the set of actions executed on an FBD program which in our case is considered the system model f bd. An example of a diagnostic trace has the following form (f bd 0 )

t1 -→ (f bd 1 ) t2 -→ ... tn -→ (f bd n ),
where (f bd k ) are states of the combined model and a k are either internal synchronization actions, time-delays or read!, execute!, and write! global synchronizations. Test cases are obtained by extracting from the test path the observable actions read! and write! as these actions contain updates on input and output variables. In summary, the output of this step is a set of ordered test cases containing inputs, actual outputs and timing information (i.e., the time parameter in the test suite is expressing timing constraints within one program).

Experimental Evaluation

In order to evaluate the proposed mutation test generation technique, we designed an industrial case study. In particular, we aimed to answer the following research questions:

-RQ1: Does mutation adequate test suites detect more faults than tests suites manually created by industrial engineers or automatically created test suites based on decision coverage? -RQ2: Are mutation adequate test suites less costly than tests suites manually created by industrial engineers or automatically created test suites based on decision coverage?

The case study setup is shown in Figure 3. From a high level view we started the case study by collecting: (i) a set of real industrial programs from a recently developed train control management system (TCMS), and (ii) manual test suites created for the above programs by industrial engineers. The studied programs were already thoroughly tested and are currently used in a set of operational trains. For all programs, test suites were also generated for weak mutation, strong mutation and decision coverage (as detailed below).

In order to measure fault detection, realistic faulty versions of the programs under test are required. However, the data set did not contain any information about what faults occurred during development, as Bombardier Transportation AB does not keep any such data in a format that could be directly collected post-mortem at this level of testing. To overcome this issue, several engineers from Bombardier were asked to manually create a number of faults for the programs considered in this study. We obtained faults from engineers at Bombardier Transportation manually introducing relevant faults in some of the programs considered in this study. Since mutation-based test generation is using an existing program implementation to guide the search, we automatically generate all tests suites using the seeded faults instead of the original programs. This corresponds to the realistic situation where an engineer has made a fault located in the program to be tested. In summary, we used a TCMS system containing 61 programs provided by Bombardier Transportation AB. These programs contained on average per program: 828 lines of IEC 61131-3 FBD code, 22 decisions (i.e., branches), 11 input variables and 5 output variables.

Manually Seeding Faults. For the TCMS programs, we provided four engineers working at Bombardier Transportation AB, who were not involved with the study with a document on doing fault seeding together with all the 61 programs. We asked each engineer to seed faults into the set of programs; we followed a specific fault seeding procedure using the IEC 61131-3 programming tools the engineers are using for developing the programs and instructed them to insert faults that were as realistic as possible. In particular, we instructed the engineers to insert any number of relevant faults, based on their experience, in the set of programs we provided as a TCMS project. We specifically instructed them to try to insert multiple faults in the same program one at the time and seed faults in at least ten programs from the total of 61. To avoid any misunderstanding, the fault seeding procedure document included information about the type of faults we were interested in: any fault that they might have encountered in their experience, as long as the interface (i.e., inputs and outputs) remained the same. This includes, but is not limited to, faults associated with variables, blocks, connections and constants. The fault seeding procedure resulted in 77 faults, versions of 33 (out of 61 in total) original programs containing a single fault (i.e., each fault contained one or more changes in the program). Each of the collected and generated test suites was executed on each of the faulty versions and its original counterpart so that a fault detection score could be calculated. Practically, each faulty variant contained one fault that had been manually seeded. A fault was considered to be detected by a test suite if the output from the faulty program differed from that of the original program.

Test Generation For each faulty program, we ran mutation and decisioncoverage test generation ten times using a random-depth-first search (RDFS) strategy with random seed (i.e., test suites are varying from run to run), each test generation run with a stopping time limit for the search of 10 minutes. The stopping criteria for the search is three-fold: achieving 100% mutation score, reaching the time limit of 10 minutes, or getting a memory exception. We chose a time limit of 10 minutes for the sake of this experiment. In addition, we used manual test suites created by industrial engineers in Bombardier Transportation from a TCMS project delivered already to customers. Manual test suites were collected by using a post-mortem analysis of the test data available. The test suites collected in this study were based on functional specifications expressed in a natural language. Practically, we considered the original TCMS programs and for each faulty program, we executed the test suites produced by manual testing for the original program. Finally, for all test suites we collected the following measures: generation time, execution time, number of test cases and fault detection score. In order to calculate the fault detection score, each test suite was executed on both the original program and its faulty counterpart. In case the results differed between the executions, the fault was considered to be detected.

Measuring Cost. We measured the cost of performing testing focusing on the unit testing process as it is implemented in Bombardier Transportation for testing the programs selected in this case study. For the TCMS system, the creation and execution of test cases is performed by the implementer of the IEC 61131-3 software. In the cost measure, we use the creation cost, the execution cost, and the result check cost. The cost does not include the required tool preparation, the reporting and the maintenance of the test suite. We consider that all cost components related to human effort are depended to the number of test cases. The higher the number of tests cases, the higher are the respective costs. We assume this relationship to be linear with a constant factor representing the average time spend by an engineer in each cost component for a test case. Practically, we measured the costs of these activities directly as an average of the time taken by three industrial engineers (working at Bombardier Transportation implementing some of the IEC 61131-3 programs used in our case study) to perform manual testing.

Experimental Results and Discussion

The case study presented us with a fault detection score and a cost measurements for each of the collected test suites (i.e., manually created test suites by industrial engineers (MAN), mutation-adequate test suites (i.e., weak-mutation testing (WM), strong-mutation testing (SM)) and automatically generated test suites based on decision coverage (DC)). The overall results of this study are summarized in the form of boxplots3 in Figure 4 and5. Fault Detection. To answer RQ1 regarding the fault detection, in terms of detection of manually seeded faults, we focused on comparing all DC, WM, SM and MAN test suites. For all programs, as shown in Figure 4, the fault detection scores obtained by manual written test suites are higher in average with 9% and 6% than those achieved by weak mutation and strong mutation respectively. The difference in fault detection is slightly greater between strong-mutation testing and decision coverage-adequate testing (i.e., a difference of almost 12% on average). To understand how manual test suites achieve better fault detection than mutation-adequate test suites, we examined if the test suites are particularly weak or strong in detecting certain type of faults. We concern this analysis to what kind of faults were detected by manual testing and not by strong-mutation test generation. From a total of 77 faults, we identified eight faults (i.e., for exemplification purposes these faults are named Fault 1-8) that were not detected by any strong-mutation test suite while being detected by manual test suites. To produce meaningful results the remaining 69 faults are not included in this fault detection analysis because there is no consistent difference between manual and strong-mutation test suites. There are some broad trends for eight faults that can be used for explaining at least the difference in fault detection between manual and strong-mutation testing. Test suites written using manual testing are able to detect all of these eight faults. Mutation test suites are achieving a poor selection of test inputs produced for detecting certain faulty behaviors; for six faults, strong mutation testing generated test suites achieving 100% mutation score while for the remaining two faults, the model checker was unable to find a test suite detecting all mutants, given the 10 minutes time limit. It seems that manual testing has a stronger ability to detect these faults than mutation testing because of its inherent advantage of relying also on the specification of the program under test. For four of the faults, multiple changes in the program have been seeded (e.g, two or more blocks and variables have been replaced, deleted or inserted). For example, Fault 1 contains three changes combining three simpler faults corresponding to the application of CRO and VRO mutation operators. Fault 2 contains a combination of seeded changes corresponding to the creation of mutants using LRO and NIO mutation operators. In addition, Faults 3 and 4 contain multiple changes that were not captured by previously defined mutation operators. On the other hand, four of the faults are first order faults containing only one change in the program. A feedback loop signal connecting one of the outputs of the programs with one of the blocks was seeded in Faults 5 and 6. On the other hand, Fault 7 contains an extra logical block that was added to the original program while in Fault 8 a constant variable has been replaced to a non-boundary value. As a direct result, we discuss in Section 5.1 the improvement of mutation-based test generation for PLC software by considering q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SM WM DC MAN (a) Overall Cost Results q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q SM(g+e) WM additional mutation operators not considered before in the literature [START_REF] Oh | Software Safety Analysis of FBD using Fault Trees[END_REF][START_REF] Shin | Empirical Evaluation on FBD Model-based Test Coverage Criteria using Mutation Analysis[END_REF] to model possible faults.

Cost. We interviewed three engineers working on developing and manually testing TCMS software and asked them to estimate the time (in minutes) needed to create, execute and check the result of a test suite. All engineers independently provided similar cost estimations. We averaged the estimated time given by these three engineers and we calculated each individual cost using the following constants: 6.6 minutes for the creation of a test case, 3.3 minutes for the execution of a test case and 2.5 minutes for the checking of the result of a test case. Practically, for answering RQ2, we used these constants and the number of test cases in each test suites to represent the average time spend by an engineer to manually test each program. In addition, for mutation-based and decision coverage-directed test generation the total cost involves both machine and human resources. We calculated the cost of generating and executing a test suite by directly measuring the time required by the tool to run the test generation and the time required to execute each test case. For the cost of checking the test result we used the same average time as for manual testing (i.e. 2.5 minutes for the checking of the result of a test case). The resulting cost measures are reflected in Figure 5. The cost of performing testing using mutation testing either weak or strong is consistently significantly lower than for manually created test suites; automatically generated test suites have a smaller testing cost (110 and 115 minutes shorter testing time on average for WM and SM respectively) than the cost of using manual test suites. A more detailed cost measurement would be needed to obtain more confidence in the cost results obtained in this study.

Discussion

To explore the results of our study we consider the implications for future work and the extent to which mutation testing for PLC programs can be improved.

Improving Mutation Testing for PLC Programs. The results of this study indicate that fault detection scores obtained by manual test suites are better than the ones achieved by mutation testing. While comparing just strong-mutation testing with manual testing, we discovered that some of these faults are not reflected in the mutation operator list used for generating mutation adequate test suites, as described in Section 3.1. From our results, we highlight the need for improving the list of mutation operators used for mutation testing of PLC software by the addition of the following new mutation operators:

-Feedback loop Insertion Operator (FIO) is inserting a signal connecting an output variable to any block that is connected with the input variables. -Logical Block Insertion Operator (LIO) is inserting a logical block between any other two logical blocks in the program. -Logical Block Deletion Operator (LDO) is deleting a logical block and connecting the inputs of this block to the next logical block in the program.

In addition there are couple of already implemented mutation operators (shown in Section 3.1) that can be improved by considering the following operators:

-Value Replacement Operator-Improved (VRO-I) is replacing a value of a constant variable value connected to a block not only with its boundary values but also with a selection of non-boundary values including 0, 1, -1. -Logical Block Replacement Operator-Improved (LRO-I) is replacing a logical block not only with logical blocks from the same category but also with other blocks with Boolean inputs (e.g., replace an AND block with an SR block).

By generating test suites that detect faults created based on these mutation operators, one could improve the goals of mutation testing for PLC programs. In addition, we recommend the use of higher-order mutation [START_REF] Jia | Higher order mutation testing[END_REF] for PLC software in order to find more complex faults.

Mutation Testing using Model Checking. Our study is the first to consider mutation testing using model checking for PLC programs written in IEC 61131-3 FBD language. Model checking is a formal technique based on state exploration that has been applied to mutation testing by either using a process named reflection [START_REF] Black | Modeling and marshaling: Making tests from model checker counterexamples[END_REF], by state machine duplication [START_REF] Okun | Testing with model checker: Insuring fault visibility[END_REF], or by explicitly evaluating the fault coverage over multiple mutants [START_REF] Gargantini | Using model checking to generate fault detecting tests[END_REF][START_REF] Godskesen | Connectivity testing through modelchecking[END_REF] thus creating test cases for manifesting fault propagation. The performance of this kind of approaches is depended not only on the model size but also the time spent on checking each and every mutated model or property against its original counterpart. This way of using the model checker for mutation testing can introduce unnecessary runs of the model checker and can considerably affect the feasibility of these approaches in practice. The method proposed in this study for the IEC 61131-3 FBD language is using a rather different approach for mutation testing, by combining all the mutants and the original model into a single combined model that is monitored dynamically using a model checking approach. By considering this way of utilizing the model checker one could potentially improve the cost of using mutation testing for other languages and models; the detection can be verified in a single run of the model checker for all mutated models rather than considering each individual case and thus removing the unnecessary model checking runs needed for detecting trivial mutants. This needs to be carefully considered in future studies and compared with other approaches on mutation testing using model checking.

Conclusions

In this paper we introduced mutation testing for PLC programs written in IEC 61131-3 programming language using a model checker. We implemented our approach in a tool and used this implementation to evaluate mutation testing on industrial programs and manually seeded faults. Our results show that mutation testing achieves lower fault detection compared to manual testing but with a significant lower cost in terms of testing time. We found out that these fault detection scores can be improved by considering some new and improved mutation operators for PLC programs as well as higher-order mutation.

Fig. 1 :

 1 Fig. 1: An FBD program with six inputs and two outputs.

Fig. 2 :

 2 Fig. 2: Overview of mutation testing for IEC 61131-3 FBD programs.

Fig. 3 :

 3 Fig.3: Overview of the experimental setup used to perform the case study.

  Average Fault Detection Score.

Fig. 4 :

 4 Fig. 4: Fault detection results for manual testing (MAN), decision coveragedirected test generation (DC), weak mutation testing (WM) and strong mutation testing (SM).

Fig. 5 :

 5 Fig. 5: Cost measurement results for manual testing (MAN), decision coveragedirected test generation (DC), weak mutation testing (WM) and strong mutation testing (SM).

  Logic Block Replacement Operator (LRO) replaces a logical block with another block from the same function category (e.g., replacing an XOR block with an OR block), -Comparison Block Replacement Operator (CRO) replaces a comparison block with another block from the same function category (e.g., replacing a Less-Than (LT) block with a Less-or-Equal (LE) block), -Arithmetic Block Replacement Operator (ARO) replaces an arithmetic block with another block from the same function category (e.g., replacing a maximum (MAX) block with a subtraction (ADD) block), -Negation Insertion Operator (NIO) negates an input or output connection (e.g., an input variable IN1 becomes not(IN1)), -Value Replacement Operator (VRO) replaces a value of a constant variable connected to a block (e.g., replacing a constant value (const = 20s) with its boundary values (e.g., const = 19s and const = 21s)), and -Timer Block Replacement Operator (TRO) replaces a timer block with an-

other block from the same function category (e.g., replacing a Timer-On (TON) block with a Timer-Off (TOF) block).

  program, and construct the structure of the program representing the set of blocks b n , set of signals s m and set of variables v p in P: b 1 ... b n , s 1 ... s m and v 1 ... v p . -For each mutant m i in M , created by changing a block, signal or variable in P , create a duplicate version of it (e.g., b 11 is a duplicate of b 1 ) having a different identifier and output than the original. This duplicate version has an interface, consisting of a name identifier. In addition, this duplicate version contains the same inputs as the original behavior, but different output variables and internal parameters in case of a mutated block. The interface is used to access both the block behavior and its duplicated version.

boxes spans from 1st to 3rd quartile, black middle lines mark the median and the whiskers extend up to 1.5x the inter-quartile range and the circle symbols represent outliers.

Acknowledgments

This research was supported by The Knowledge Foundation (KKS) through the following projects: (20130085) Testing of Critical System Characteristics (TOCSYC), Automated Generation of Tests for Simulated Software Systems (AGENTS), and the ITS-EASY industrial research school.


Masaki Suzuki
email: masaki-suzuki@kddilabs.jp

Quentin Plessis
email: qu-plessis@kddilabs.jp

Takeshi Kitahara
email: kitahara@kddilabs.jp

Masato Tsuru
email: tsuru@cse.kyutech.ac.jp

Distribution Visualization for User Behavior Analysis on LTE Network

Keywords: LTE, C-Plane analysis, Self-Organizing Map (SOM)

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

Introduction

Mobile network operators (MNOs) are responsible for providing high quality of communication services. It is very important for them to monitor the communication quality. For this purpose, MNOs tackle to immediately detect the degradation of the communication quality when any incidents occur. The existing approaches are generally either the log-based or the conformance-based. In the log-based approaches, a system monitors messages and system logs of equipment in the LTE network [START_REF] Agrawal | On the design of element management system for node Bs in a 3G wireless network[END_REF]. The system detects hardware errors and link errors. On the other hand, in the conformance-based approaches, a system detects unfamiliar sequence of messages referring the specifications of 3GPP standard. However, there exists a degradation of the communication quality occurring without any error messages or nonconformity. For instance, ping-pong handover is a general phenomenon in mobile networks, which causes inefficient network performance and communication quality [START_REF] Li | Adaptive Frame Switching for UMTS UL-EDCH -Ping-Pong Avoidance[END_REF][START_REF] Kim | A Mobility Management Technique with Simple Handover Prediction for 3G LTE Systems[END_REF]. When one User Equipment (UE) which is moving close to the fringe between multiple evolve node Bs (eNBs) and connecting to one of them, it handovers from one eNB to another eNB, then it often immediately connects back to the former eNB. In the case where the UE stays around the fringe, it sometimes repeatedly handovers from/to these eNBs. In this situation, there exists no evolved packet core (EPC) equipment errors. However, this phenomenon still causes unnecessary control messages in EPC and degrades the communication quality. In such a case, MNOs hardly detect the degradation unless customers report the problem to them.

For the first step to study detecting the degradation without any errors, MNOs have to know how users behave in EPC. In order not to lose the generality, the users' behavior analysis should be exhaustive and comprehensive. However, since EPC signals through various interfaces between function nodes in EPC are mixture of different protocols and IDs, it is difficult to trace users' behavior sequentially.

In this paper, we report a preliminarily implemented system that captures and analyzes C-Plane signals in EPC, quantifies users' behavior, and visualizes the distribution of users' behavior. Then we introduce a case study with the actual C-Plane signals and a typical example for cluster of degraded situation of users' behavior.

Related Works

There exist several studies about users' behavior analysis in mobile networks. In [START_REF] Gorbil | Modeling and Analysis of RRC-Based Signalling Storms in 3G Networks[END_REF], the authors analyze signaling storms based on radio resource control (RRC) protocol. In order to detect anomaly and malicious users' behavior causing signaling storms, [START_REF] Gorbil | Modeling and Analysis of RRC-Based Signalling Storms in 3G Networks[END_REF] firstly models and analyzes the patterns of signals in RRC protocol. Then, it identifies the specific patterns. In [START_REF] Wang | A Control-Plane Traffic Analysis Tool for LTE Network[END_REF], the authors focus on retrieving radio access information from S1-MME and S11 interfaces. As an example, the authors summarize the time transition of the duration of radio access bearer establishment.

To the best of our knowledge, there exists no study analyzing or visualizing users' behavior in EPC. Therefore, for the first step of research, we tackle to visualize users' behavior based on C-Plane signals in EPC.

Implementation

Figure 1 LTE architecture. Figure 2 System architecture. In order to analyze users' behavior in EPC, we implement a distribution visualization system in an actual LTE network which is standardized by 3GPP [START_REF]3rd Generation Partnership Project, General Packet Radio Service (GPRS) enhancements for Evolved Universal Terrestrial Radio Access Network (E-UTRAN) access[END_REF]. Figure 1 briefly depicts the LTE architecture regarding to C-Plane signals. In our implementa-tion, we focus on the signals through S1-MME, S10 and S11 interfaces. They are a mixture of S1 application protocol (S1AP) and Evolved general packet radio service tunneling protocol for control plane (GTPv2-C). Figure 2 shows the architecture of our implementation. Firstly, the capture server captures signals. Secondly, the signal analyzer extracts users' state transition from capture files. Thirdly, the statistics monitor quantify users' behavior based on users' state transition. Finally, the distribution visualizer draws users' behavior distribution using self-organizing map (SOM).

Capture of C-Plane signals and signal analysis Table 1 States in the system

Table 2 Examples of state transition

The process in our implementation starts with the capture of the signals. The implemented system groups the signals by user, then constructs signal sequences by user. After that, it extracts specific patterns of signal sequences. Note that, the implemented system does not identify the specific user. It can only distinguish users using a temporary assigned identifier. Since the temporary identifier is valid for a certain duration of time, the implemented system can trace users' behavior for a short time. Thus the implemented system cannot follow any specific user for a long time, e.g. several hours or longer.

Based on the signal sequences, the system constructs a state transition graph. The state transition graph consists of 5 elements as follows. The input is specific patterns of signaling messages extracted from S1AP and GTPv2-C signals. The states are defined according to 3GPP standard, and determined by the combination of current state and input. In the system, the initial and final state of the transition are ignored since, in actual LTE networks, the initial state should always be the same state and the final state should not be naturally defined. Tables 1 and2 show the lists of states and examples of state transition respectively.

Statistics monitor

After that, the implemented system calculates statistics values in order to quantify the behavior of each user based on his/her state transition. In the implementation, in order to characterize the continuous-time state transition of a user, we adopt the state transition probability matrix (p(n,m)) as well as the average and the variation coefficient of the dwell time , at state n in a transition from state n to state m. The probability p(n,m) from state n to state m is calculated in Equation 1,

, . (1)
The dwell time , is calculated by state transition as in Equation 2,

, , (2)
where, and is the arrival time at state m and n in i-th state transition from state m to n respectively. To gather these values, we describe users' behavior with a multi-dimensional vector. As respecting the definition of states, number of possible state transition is 600. Since we adopt 3 different statistics values, users' behavior described in a 1,800-dimensional space in our implementation.

Distribution visualizer

In order to visualize the distribution in a multi-dimensional space, the distribution visualizer uses the self-organizing map (SOM) [START_REF] Kohonen | The self-organizing map[END_REF]. SOM is an artificial neural network using unsupervised learning to construct a two-dimensional space representing a multi-dimensional space. We can intuitively see the distribution of the users' behavior by mapping the distribution in a multi-dimension into a two-dimension.

According to the SOM algorithm, the distribution visualizer firstly define the vector space based on the entire input data. Secondly, the distribution visualizer plots the quantified user's behavior in an n-dimensional space one by one. Then it transforms the distribution into a two-dimensional space. In the process of the transformation, it draws regular grid of circles (namely, units) in the two-dimensional space. Each unit represents principal components and each plot is located in the closest circle so that the more the behaviors are similar, the closer they are located. The visualizer highlight the specific condition of users in the case where they are labeled in advance and we can compare different conditions of users intuitively.

4

Case study

In order to validate the result of the implemented system and assess its usefulness, we visualize users' behavior based on 24 hours of the actual anonymized C-Plane signals in a large urban area in Japan. In this case study, we intuitively identify the fundamental characteristics of specific users who had experienced ping-pong handovers and labeled in advance. Firstly, the signal analyzer parses the captured signals and constructs signal sequences by users. Then, it extracts state transitions. Figure 3 depicts the state transition diagram. In the figure, the indexes of the nodes are the indexes of the states in Table 1 and the width of edges are the probability of state transition (p(n,m)). For the readability, we ignore the edges which p(n,m) is less than 0.10 in the figure. According to the state transitions, the statistics monitor quantifies the users' behavior in terms of p(n,m), the mean value and variation coefficient of dwell time (t (i) (n,m)). Figure 4 Indexes of units. We define the input space using the entire 24 hours of input data. We prepare 100 units to describe the input space in the 2-dimensional space and the indexes of units are numbered in a left-to-right and bottom-to-top fashion as described in Fig. 4. Figure 5 depicts distribution maps of 0 am, 6 am, 12 pm and 6 pm in the day. In each unit, we plot users' behavior of each hour in gray color. Then we highlight users' who had experienced ping-pong handovers in the period of time in red color.

According to the figures, the number of visualized users are varied by time and the distribution of users are different especially between 0 am and 6 am. The highlighted users, however, located in similar units. Focusing on those unit 2, 3, 13, 14, 15, 23 and 24, the components of them commonly include variation coefficient of ∈ 11,24 , p(11,24), p(11,23), p [START_REF] Li | Adaptive Frame Switching for UMTS UL-EDCH -Ping-Pong Avoidance[END_REF]24). Since ping-pong handovers mean frequent handovers, it is quite understandable those users are likely to belong to those units which include transition from 11 to 23 or 24. However, the visualized results indicate a phenomenon that a number of users who experience ping-pong handovers also start X2 handover right after the initial context setup, which is unexpected by MNOs. Our system enables to highlight the unknown characteristics of ping-pong handovers.

Conclusion and future works

In order to analyze the users' behavior in EPC, we implemented the visualization system for user's behavior distribution. We draw the distribution maps using implemented system with the actual C-Plane data. As future works, we will deeply analyze users' behavior based on multi-hop transitions of states.

Figure 3

 3 Figure 3 State transition diagram.Figure4Indexes of units. We define the input space using the entire 24 hours of input data. We prepare 100 units to describe the input space in the 2-dimensional space and the indexes of units are numbered in a left-to-right and bottom-to-top fashion as described in Fig.4. Figure5depicts distribution maps of 0 am, 6 am, 12 pm and 6 pm in the day. In each unit, we plot users' behavior of each hour in gray color. Then we highlight users' who had experienced ping-pong handovers in the period of time in red color.According to the figures, the number of visualized users are varied by time and the distribution of users are different especially between 0 am and 6 am. The highlighted users, however, located in similar units. Focusing on those unit 2, 3, 13, 14, 15, 23 and 24, the components of them commonly include variation coefficient of

Figure 5

 5 Figure 5 Distribution maps of (a) 0 am (b) 6 am (c) 12 pm and (d) 6pm.