Antoine Houdard

Statistical Modeling of the Patches DC Component for Low-Frequency Noise Reduction

This work is devoted to patch-based image denoising. Assuming an additive white Gaussian noise (AWGN) on patches, we derive corresponding models on centered patches and on their DC components. Then we propose to treat the two components separately. Finally, we provide experiments with the recent denoising method HDMI [1] and we show that our approach improves the denoising quality, particularly for residual low frequency noise.

Introduction

Context In this paper, we focus on image denoising, which aims to estimate an image p u from its noisy observation v " u `e P R n ,

where e " N p0, σ 2 I n q is an additive white Gaussian noise (AWGN) and u the underlying clean image. Some of the recent denoising methods are based on a statistical modeling of the image patches [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Afonso M Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF][START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (hdmi)[END_REF]. For each patch i P t1, . . . , nu seen as vector of size p " s ˆs, model (1) yields

Y i " x i `Ni , (2)
where Y i P R p is the observed random vector modeling the i-th patch, x i P R p is the underlying clean patch and N i " N p0, σ 2 I p q is a Gaussian white noise. The idea is then to set a prior model on the clean patch x i seen as a realization of a random vector X i . The model therefore rewrites

Y i " X i `Ni , (3)
and Bayes' theorem yields the posterior X i |Y i . Finally, each clean patch can be estimated with the conditional expectation

p x i " ErX i |Y i " y i s. (4) 1
Convenient priors for computing this estimator (4) are Gaussian priors [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] or Gaussian mixture models (GMM) [START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF][START_REF] Afonso M Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF][START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (hdmi)[END_REF]. The use of these priors has been widely studied and it appears that the covariance matrix of these models can encode local structures up to some contrast change [START_REF] Delon | Gaussian Priors for Image Denoising[END_REF]. This permits to regroup more patches under the same Gaussian model and then allows for a better estimate of its parameters. There is however a drawback: grouping patches in this way makes the mean of the model less informative. This yields an estimate for each patch that has some bias. This produces the low frequency residual noise that appears in the result of model-based patch-based denoising methods. Figure 1 (b) illustrates this phenomenon in the case of the HDMI method [START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (hdmi)[END_REF], with strong noise and small patches. A large part of this low frequency noise seems to come from a poor estimation of the mean of each patch. Indeed, the image (d) from figure 1 made up of the mean of each patch of the noisy image (a) has the same patterns than the image (c) which is made up of the mean of each patch of the denoised image (b). Moreover, replacing the mean of each denoised patches from (b) with the true mean of the oracle (f), yields a denoised image (c) that is way better, both in terms of Peak Signal to Noise Ratio (PSNR) and visual quality, than the image (b). In addition, some methods from the literature [START_REF] Afonso M Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] also seem to suggest that removing the mean -also called the DC component -of the patches may improve the denoising quality.

Proposed work

In this work, we propose to study the decomposition of the patches into the DC component and the centered component for denoising purposes. To do so, we define the centered observed random variable Y c i " Y i ´s Y i 1 p , where

s Y i " 1 p p ÿ j"1 Y i pjq, (5)
is the mean of Y i and 1 p " p1, . . . , 1q P R p . The model (3) can then be divided into the two following problems

s Y i " s X i `s N i P R, (6)
and

Y c i " X c i `N c i P R p . (7)
We propose to model the noise components N c i (section 2) and s N i (section 3) of these two problems. Then we suggest in section 4 various solutions that can enhance the denoising results of the patch-based denoising method HDMI [START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (hdmi)[END_REF], and we provide numerical experiments for these solutions.

Modeling the centered noise

The centered noise is defined by N c i " N i ´Ď N i 1 p and then the j-th entry of N c i is since N i is a Gaussian random vector, N c i is also Gaussian. Then, we can compute its mean and its covariance matrix coordinate by coordinate. That gives for the mean

N c i pjq " N i pjq ´1 p p ÿ k"1 N i pkq, (8)
E rN c i pjqs " E rN i pjqs ´1 p p ÿ k"1 E rN j pkqs " 0. (9)
And for the covariance matrix, we have @k, l P t1, . . . , pu

E rN c i pkqN c i plqs " E rN i pkqN i plqs ´1 p p ÿ m"1 E rN i pmq pN c i plq `N c i pkqqs `1 p 2 p ÿ m"1 p ÿ n"1 E rN i pmqN c i pnqs " σ 2 ˆδkl ´2 p `p p 2 " σ 2 ˆδkl ´1 p ˙.
Finally, we have

N c i " N p0, Σ N c i q with Σ N c i " σ 2 p ¨p ´1 ´1 ¨¨¨´1 ´1 p ´1 . . . . . . . . . ´1 ´1 ¨¨¨´1 p ´1 ‹ ‹ ‹ ‹ ' . (10)
Since Σ N c i is a real symmetric matrix, there exists an orthonormal basis that diagonalizes it. Given that its eigenvalues are p (of multiplicity p ´1) and 0, we can build an orthogonal matrix Q such that

Σ N c i " Q ˆσ2 I p´1 0 0 0 ˙QT . (11)
It is worth noticing that the unit eigenvector corresponding to the eigenvalue 0 is 1 ? p p1, . . . , 1q. The change of basis Q T applied to the centered noise then yields Q T N c i " N p0, diagpσ 2 I p´1 , 0qq. Finally, centering the noise implies a dimension reduction. The total variance of the centered noise defined as

TVarpN c i q " Er}N c i } 2 2 s (12) " TrpΣ N c i q " pp ´1qσ 2 (13) " p ´1 p TVarpN i q (14)
is also reduced by a factor pp ´1q{p.

Modeling of the DC component

Since the initial noise model on a patch N i is a Gaussian vector, its mean is a Gaussian random variable s N i " N p0, σ 2 p q. Then, reshaping problem (6) as an image yields the new image denoising problem with additive Gaussian noise

s Y " s X `s N , (15)
where s Y , s X and s N are the images of R n whose values at pixel i are s Y i , s X i and s N i . The major difference is that the noise is now colored. Indeed, if we consider two random variables s N i and s N j within a same area of s ˆs pixels, they are issued from two overlapping patches and thus are not independent.

However, we can still perform patch-based image denoising on this problem: let us consider patches of the same size p " s ˆs from this new image. We define the patches Z i " π i p s Y q, W i " π i p s Xq and M i " π i p s N q, where π i is the i-th patch extraction operator. We consider the patch noise model

Z i " W i `Mi , (16)
with M i modeling the noise. Since M i " p s N i1 , . . . s N ip q, all its entries are linear combinations of the noise components of the problem (1) that are i.i.d following N p0, σ 2 q. Therefore, all linear combinations of entries of M i are also Gaussian. This shows that M i is a Gaussian vector. We can now compute its mean and covariance matrix.

The mean of M i is obviously 0 p the coefficients of the covariance matrix Σ Mi are given by

pΣ Mi q kl " E " s N i k , s N i l ‰ " σ 2 p 2 C kl , (17)
where C kl is the number of common pixels between the two patches of the original image from which s N i k and s N i l are derived. This yields after counting

Σ Mi " σ 2 p 2 B b B (
‹ ‹ ‹ ‹ ' , ( 19
)
and b is the Kronecker product. In order to use a denoising method that has been designed for AWGN, we want to perform a change of basis for the data.

To do so, we study the structure of Σ Mi . First, we show that B is symmetric positive-definite. Using the Sylvester's criterion, it is sufficient to show that all of the leading principal minors of B are positive.

Finally, developing the determinant with respect to the first column and repeating the two last steps yields

m d " p2s ´d `1q2 d´2 ą 0. (23)
This shows that B is a positive-definite matrix. Then B b B is also symmetric positive-definite as the Kronecker product of two positive-definite matrices and the Cholesky decomposition yields an invertible matrix L such that B b B " LL T . Therefore, the problem

L ´1Z i " L ´1W i `L´1 M i , (24)
is an AWGN problem with noise variance σ 2 {p 2 I p and a denoising method such as HDMI can be used to find an estimate { L ´1W i that gives an estimate of

W i by x W i " L { L ´1W i .

Experiments

In this section, we take advantage of the previous modeling in order to propose various strategies for improving the denoising result of model-based patch-based denoising methods. For this testing part, we propose to use the HDMI method [START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (hdmi)[END_REF] which has the advantage of using only statistical tools. The principle of this method is rather simple:

1. It uses a GMM with intrinsic lower dimensions to model clean patches;

2. This model is inferred on the noisy patches with an expectation-maximization (EM) algorithm;

Denoising the DC component

In section 3, we have proposed to reshape the DC component of all patches into an image and to extract patches of this image in order to perform patch-based denoising on it. We showed that the noise model on these patches is an AWGN in a given basis L. We can therefore apply the HDMI method directly on (24).

Figure 2 shows the result of the denoising of the DC images from the images simpson and lena with a noise of standard deviation σ " 30{255. Note that the results are quite good since the problem ( 24) is an easier problem than the original one (3). Indeed, the dynamic of the DC image is quite the same as the one of the original image whereas the dynamic of the noise is reduced by a factor p. Therefore, the signal-to-noise ratio of the problem ( 16) is about p times larger than the one of the original problem (3). Finally, with this step, we obtain for each patch i an accurate estimate x Ď X i of its DC component.

Denoising the patches

In order to perform the final patch denoising, we expose hereafter two strategies.

For each of these strategies, we propose numerical examples and we discuss the improvements they bring for the denoising task. Images denoised with method 4.2.1 (lena 31.19dB, simpson 32.58dB).

Correction of the DC component afterwards

The first idea is to run the HDMI method on the original patches y i , yielding for each patch i an estimate x X i , and then to correct this estimate with the estimated DC component x Ď X i . That is, to consider the final estimate

Est 1 pX i q " x X i ´Ď x X i 1 p `x Ď X i 1 p . (25)
This strategy improves the denoising quality compared to the original denoising, as illustrated in Figure 3. The low frequency noise is visually reduced and in term of PSNR, the results show an improvement of « 0.1dB for both images. However, this does not use the modeling of the centered noise from section 2, and the grouping of the patches from the HDMI model remains the same.

Learn a model on the centered patches

The second strategy is to learn the model on the centered patches X c in order to group together patches with various contrasts. To do so, we can apply the change of basis Q T from section 2 on the problem (7). This yields

Q T Y c i " Q T X c i `QT N c i , (26)
with the last line of this equality being Ď Y c i " Ď X c i which is always true since Ď Y c i " Ď X c i " 0 by definition. Removing this last dimension then yields an AWGN problem in dimension p´1 that can be denoised with the HDMI algorithm. This gives us an estimate of the centered patch x X c i and therefore the final estimate Figure 4 shows the result of this strategy. The advantage of learning the model on the centered patches is that the patches with same structure but different contrasts are grouped together (figure 4). This can make the estimate of the covariance matrix easier. However, this centered problem has a lower signalto-noise ratio than the original problem since the centering of the patches has reduced the dynamic. Therefore, the algorithm is more likely to missclassify some fine textures. This strategy still leads to a visual and PSNR improvement but seems not to be better than the previous solution.

Est 2 pX i q " x X c i `x Ď X i 1 p . (27)

Discussion

Finally, we propose in table 1 the results of the two presented strategies for the three images simpson, lena and barbara with different level of noise. This shows that the first strategy always improve the denoising quality, with a larger improvement for images that have less textures like simpson. The second strategy seems also to work well in this precise case but fail to improve the result in the case of barbara which has a lot of fine texture. The other trend that appears is the better improvement for high variance noise than for low variance noise. This is due to the use of the same patch size for all experiments. Indeed, for a given patch size, the residual low frequency noise increases with the variance of the initial noise.

Conclusion

In this work, we have studied the effect of patch centering for model-based patchbased denoising methods. For this purpose, we have proposed a modeling of the centered noisy patches and a modeling of the DC component of the noise. These modeling have led us to strategies for improving the result of existing denoising methods, especially to reduce low-frequency residual noise. The results obtained with the HDMI method on grayscale images show improvements, both visually and in term of PSNR.

In future work, we would like to study the links between this approach and multiscale frameworks. Then, in a second step, we would like to propose a discussion on the strategy to be adopted for color images.

Figure 1 :

 1 Figure 1: (a) noisy image (standard deviation 30{255). (b) image denoised with HDMI [1] (patches 7 ˆ7), PSNR = 31.92 dB. (c) image from the denoised patches of (b) with DC component corrected with oracle value (f), PSNR = 33.78 dB. (d) DC component of the patches of (a). (e) DC component of the patches of (b). (f ) DC component of the oracle image (ground truth).

Figure 2 :

 2 Figure 2: Left: noisy DC image. Middle: denoised DC image ('simpson' PSNR = 40.64 dB , 'lena' PSNR = 40.89 dB). Right: oracle DC image. Top line is the from the simpson image and bottom line from the lena image both with noise of standard deviation σ " 30{255.

Figure 3 :

 3 Figure 3: From left to right. Noisy images (standard deviation 30{255). HDMI resulting groups. Images denoised by HDMI (lena 31.09dB, simpson 32.46dB). Images denoised with method 4.2.1 (lena 31.19dB, simpson 32.58dB).

Figure 4 :

 4 Figure 4: From left to right. Noisy images (standard deviation 30{255). Resulting groups of method 4.2.2. Denoising with method 4.2.2 (lena 31.14dB, simpson 32.58dB).

  These minors of size d P t1, . . . , su are given by

	Adding the first column in the last one yields
		ˇˇˇˇˇˇˇˇˇs		ps ´1q	¨¨¨1	ˇˇˇˇˇˇˇˇˇ,
	m d " p2s ´d `1q	ps ´1q . . .		s . . .	. . . . . .	. . . 1	(21)
			ps ´d `1q	¨¨¨ps ´1q 1
	then subtracting the second and the last column to the first one gives
				ˇˇˇˇˇˇˇˇˇ2 ps ´1q	¨¨¨1	ˇˇˇˇˇˇˇˇˇ,
	m d " p2s ´d `1q	0 . . .	s . . .		. . . . . .	. . . 1
				0	¨¨¨ps ´1q 1
		ˇˇˇˇˇˇˇˇˇs	ps ´1q		¨¨¨ps ´d `1q	ˇˇˇˇˇˇˇˇˇ,
	m d "	ps ´1q . . .		s . . .		. . . . . .	. . . ps ´1q	(20)
		ps ´d `1q		¨¨¨ps ´1q	s

Table 1 :

 1 Results in PSNR (dB) of the enhancement methods from section 4.2.1 and 4.2.2 for the images simpson, lena and barbara with different noise levels.

	Image	σ HDMI [1] sec. 4.2.1	diff.	sec. 4.2.2	diff.
		10	38.89	38.93	+0.04	38.93	+0.04
	Simpson	20	34.88	34.97	+0.09	34.96	+0.08
		30	32.46	32.58	+0.12	32.58	+0.12
		10	35.81	35.82	+0.01	35.81	+0.00
	Lena	20	32.86	32.91	+0.05	32.89	+0.03
		30	31.09	31.19	+0.10	31.14	+0.05
		10	34.81	34.82	+0.01	34.79	-0.02
	Barbara	20	31.42	31.45	+0.03	31.42	+0.00
		30	29.38	29.42	+0.04	29.36	-0.02

The clean patches are estimated with the conditional expectation (4), which has a closed-form and is numerically stable (proposition 1 of[START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (hdmi)[END_REF]).


Dmytro Shytyi

High availability of VNF Orchestrator

Keywords: NFV, Tacker, resiliency, nodes recovery, placement. !

Network Function Virtualization allows transformation from physical proprietary hardware to general purpose servers with virtual services. Depending on the available resources in the network, virtual services (functions) could be placed accordingly. The special entity is required too place functions and manage resources in the network. This entity is an Orchestrator. Usage of Orchestrator raises the problem of which action to undertake in case of orchestrator failure. The problem transforms into questions like: where to place these instances? How to manage fail overs? How many instances should be placed? The Orchestrator has risk of failure thus, to provide the high availability, this paper consider the fail over and placement problems.

INTRODUCTION

I NCREASING demand on different network services makes network operators and big companies to increase the number and variety of complex proprietary hardware. The combination of increasing operational costs, complex integration and maintenance of such hardware is strong limitation of increasing services. Network function vitalization introduces the different way to architect networks. This technology aims to convert different network hardware with specific functions as a virtual services that are running on the general purpose servers (GPS), network nodes such that those services could be moved to a different locations in the network instead of needing the installation of new equipment. To make available the transition from the various complex proprietary hardware to easy managed software (virtual functions/services) on the GPS, Network Function Virtualization platforms as OpenStack [START_REF] Sefraoui | OpenStack: Toward an Open-Source Solution for Cloud Computing[END_REF] appeared.

While services require to be properly managed, coordinated, arranged and resources require to be aggregated, problems related to automation of network, storage, performance and provisioning appear. Frequently automation closely relates to an orchestration (coordination of the resources and networks needed to set up cloud-based services, applications). Correspondingly the entity responsible for orchestration tasks without requiring direct human interaction is defined as orchestrator.

To be able to perform orchestration tasks OpenStack includes the Tacker (Manegement and Orchestration) service [14] based on MANO framework [START_REF] Mijumbi | Management and Orchestration Challenges in Network Functions Virtualization[END_REF]. Normally Tacker server could be responsible (supervision) for multiple Open-Stack clusters. As any component it has risk of failure and in case of Tacker instance failure further operations on NFV will be impossible to perform. Thus our goal is to provide high availability for Tacker, that we are going to achieve with appropriate placement of Tacker servers and recovery method in case of failure.

RELATED WORK

The latest related studies on controller placement performed in Software Defined Networks. The most common objective for the problem is to minimize controller-switch latency studied by Heller et al. [START_REF] Heller | The Controller Placement Problem[END_REF], where they presented the correlation between number of SDN controlers, their locations, and controller-to-switch latency's. Jimenez et al. [START_REF] Jimenez | On the controller placement for designing a distributed SDN control layer[END_REF] performed study on controller placement in SDN with k controllers by taking into account metrics such as latency. Yao et al. [START_REF] Yao | On the Capacitated Controller Placement Problem in Software Defined Networks[END_REF] presents placement algorithm by taking into account capacity of SDN controller. Particularly, it has been widely studied to locate a controller close to switches. Finally, mechanisms to perform failover were studied by Obadia et al. [START_REF] Obadia | Failover Mechanisms for Distributed SDN Controllers[END_REF]. They described and tested a greedy failover strategy for neighbor active SDN controllers to take over the control of orphan switches in SDN netowrks.

Topology Model

In terms of methodology that includes fail tolerance approach and controller placement, also creation of network topology model considered in the section 5.1 to perform simulations. Related work in this area includes a list of models that were studied and helped to define our model.

Balanced tree model gave an inspiration to use structure inside domains as balanced tree. However it did not satisfy desired parameters of model as presented tree has one root node.

Dorogovtsev-Goltsev-Mendes graph [16] that starts by creating three nodes and tree edges, making a triangle, and then add one node at a time. In this model each time a node is added, an edge is chosen randomly and the node is connected via two new edges, however even with a fact that this model distribution follows Power Law, it does not consider the case when Tacker servers could reach each other through intermediate nodes (OpenStack clusters).

Margulis-Gabber-Galil model [START_REF] Goldreich | Basic Facts about Expander Graphs[END_REF] satisfied clustering structure and not directly connected Tacker instances, however it does not present tree-type structure of OpenStack clusters for each Tacker Server.

Caveman graph model [START_REF] Schaeffer | Graph clustering[END_REF] satisfied multiple domains and Tacker servers were not connected directly, however it still did not allow to obtain tree structure. There were more models analyzed that could satisfy some requirements but not all of them. Those models are Waxman graph model [START_REF] Thomas | Generation and Analysis of Random graphs to model internetworks[END_REF], Random Powerlaw Tree model [START_REF] Carvalhoa | A generative power-law search tree model[END_REF], Watts-Strogatz Graph model [START_REF] Watts | Collective dynamics of small-world networks[END_REF].

Therefore we propose our tree relaxed model that inspired from the tree model and has a feature similar to Caveman relaxed graph model [START_REF] Judd | Behavioral Conflict and Fairness in Social Networks[END_REF].

Placement algorithm

The part of Tacker resiliency methodology is a placement algorithm. Popular algorithms, similar to the algorithm proposed in this paper, such as k-means, k-medioids, kcenter are specializing on specific metric or properties thus have some differences. For example k-means generally is used on the euclidean space and taking into account special unknown data. These clustering algorithms are from the set of algorithms that are splitting the dataset S into several parts (1, ..N ∈ S) and try to minimize the distances d(u, c) from units u ∈ Si of the part S i to the unit c ∈ S i Main advantages of presented in this paper algorithm is the simplicity and execution time. However as a heuristic algorithm it does not always give an optimal solution. Also the number of domains to split the dataset should be predetermined. Finally the algorithm spends redundant processing time to calculate distances between arch of the k domains and N Since most nodes after several iterations stay in the same domains, there are unnecessary redundant computation of shortest paths. Thus the presented algorithm could be improved during further work.

To cope with unknown number of domains there exist approaches that could be applied. Such approaches are silhouette method, the elbow method [START_REF] Kodinariya | Review on determining number of Cluster in K-Means Clustering[END_REF].The Elbow method could be described as a function of the number of clusters: the final number of clusters will be the one, that does not give a much better result if another cluster is added. Another way is the Silhouette method, that provides a graphical representation of how well each object lies within its cluster.

One of the most recent active field where related work appears -is the Software Defined Networks (SDN) controller resiliency. However, some parameters such as time constrains and frequency of operation execution/message exchange in SDN highly differs from Tacker and OpenStack clusters message exchange density. For example the SDN controller placement highly depends on optimal placement of the controller. For example OpenFlow switches ask controller with high frequency how to operate on new packets. Density of data exchange between OpenStack clusters -Tacker server (to orchestrate, to manage vitalized services, resources) and density of packets exchanges in SDN between controller and switches (while new "unknown" packet appears on switch) is significantly differs. Thus Tacker-OpenStack cluster communication model is less sensitive for close to optimal placement than SDN communication model.

There exist different approaches to determine domains such as Girvan-Newman method and Louivan method. Both of them try to distinguish domains by using different techniques. However they still require modifications to provide the center of domains. The Girvan-Newman method [START_REF] Girvan | Community structure in social and biological networks[END_REF] focuses on the edges that are less central, those that are between domains. Instead of adding the strongest edges to empty vertex set, they are constructed by progressively removing edges from the original graph. Where edge e ∈ E betweenness is defined as the number of shortest path between pairs of the vertices that run through the e. It is a measure of influence of a node over the flow of information between the nodes. However it does not give proper results as it determines some of nodes from different physical locations (different domains) as those that should be carried by the same Tacker server.

Another popular algorithm is based on Louvain method [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. It is divided in two phases that are repeated iteratively. First, a different community is assigned to each node of the network. So, in this initial partition there are as many communities as there are nodes. Then, for each node considered the neighbours j of i and evaluated the gain of modularity that would take place by removing i from its community and by placing it in the community of j. The node i is then placed in the community for which this gain is maximum (in case of a tie we use a breaking rule), but only if this gain is positive. If no positive gain is possible, i stays in its original community. This process is applied repeatedly and sequentially for all nodes until no further improvement can be achieved and the first phase is then complete. However it also could not provide the proper domains determination, in term of current presented approach (network topology)

Contrary to the partitioning algorithm type (Our algorithm), hierarchical clustering is another basic type of domain determinition.. Hierarchical clustering create a decomposition of D such that the decomposition is represented by a tree that iteratively splits D into smaller subsets until each subset consists of only one object. In such a tree, each node represents a cluster of D.

Hierarchical algorithms do not require number of domains as an input. However as shown in paper M. Kaur et al. [START_REF] Kaur | Comparison Between K-Mean and Hierarchical Algorithm Using Query Redirection[END_REF] hierarchical algorithms are slower than "K-Means"type algorithms, thus slower than our algorithm. As a hierarchical algorithms there exists Density-based spatial clustering of applications with noise (DBSCAN) [START_REF] Ester | A Density-Based Algorithm for Discovering Clusters inLarge Spartial Databases with Noise[END_REF] that also does not require to specify the number of clusters in the data as Hierarchical algorithms. Another algorithm is "Ordering points to identify the clustering structure" (OPTICS). The basic approach of OPTICS is similar to DB-SCAN, but instead of maintaining a set of known, but so far unprocessed cluster members, a priority queue (e.g. using an indexed heap) is used.

Contribution

Our contribution is a methodology that provides high availability for Tacker service, more precisely the combination of different techniques that together give desired result specifically to the Tacker service.

The methodology consists of combination of 3 parts:

1) first part answers the question: "How to handle the failure of Tacker server(s)."; 2) second part consider the appropriate location of Tacker Service to perform NVF management and orchestration. 3) last part includes creation of model where network topology represented by a graph to estimate the first two parts.

Moreover as in the studied case the Tacker service has specific requirements, parts of methodology have differences from general cases. Specifically to Tacker as many OpenStack clusters as possible should be connected directly to the Tacker servers, however the case where OpenStack clusters are connected to the Tacker server through another OpenStack cluster also should be considered. Therefore we propose the modification in placement algorithm: a central node of each domain on each step in our algorithm we chose according to the Relative Point Centrality (Algorithm 2) of the node in the network, for example another metric like "degree centrality" does not satisfy the Tacker service requirement, as in the last case the central node will be chose according to the biggest degree, but in our case we also have to consider 2 and more hops connections between OpenStack clusters and Tacker servers.

WHAT IS TACKER

Tacker -is an implementation of ETSI MANO (Management and Orchestration) architectural framework whose mission is to Orchestrate network services [START_REF]Network Functions Virtualisation (NFV); Management and Orchestration[END_REF]. Tacker identifies the next blocks: Network Function Virtualization (NFV) Orchestrator,Virtual Network Function (VNF) Manager, NFV Catalogue (NFVC).

The NFV Orchestator (NFVO) manages the lifecycle of network services and performs the orchestration of resources across multiple VIMs. More precisely, the NFVO functions are: network service installation and lifecycle management (update, scale, collection of performance statistics and events); management of network services deployment templates; validation of virtual resource requests from VNF Managers; create, update, delete VNF Forwarding Graphs (network services topology); vitalized resources management based on geolocation, resource usage and allocation policies; collect the information about resource usage by VNF instances.

The VNF Manager (VNFM) manages the lifecycle of VNF instances. Each VNF instance has an association with VNFM and VNFM can be assigned to a multiple instances. The functions of VNFM on assigned VNFs are: item update or upgrade the VNF ( i.e. VNF's software update and configuration changes); create VNF; change the VNF's configuration/capacity; release associated to VNF resources; automated or assisted healing of VNF instance.

The NFV Catalogue (NFVC) is representation of the repository with all VNF packages. It supports the management of the VNFD by the NFVO operations.

Another module Virtualised Infrastructure Manager (VIM) is used for managing and controlling the network, computing and storage resources. The VIM is not a Tacker block, however it is important to mention that OpenStack acts as a VIM. The VIM directly interconnects with Tacker participates in orchestration of the upgrade, release, allocation and usage optimization of resources. This block manages the the association of physical resources such as storage, network and computational facilities to the virtualised resources.

ORCHESTRATOR RESILIENCY/PROBLEM STATE-

MENT

As mentioned before Tacker servers supervise the Openstack clusters and in case of Tacker instance failure further management and orchestration on NFV in supervised OpenStack clusters will be impossible. Resilience to Tacker orchestator failure could be achieved by deploying additional Tacker instances on the network. When several Tacker instances should be deployed among domains (set of Open-Stack Clusters), we consider "resource allocation" problem. Where Tacker instances are "resources" that should be "allocated" among competing domains. To allocate Tacker instances appropriately we need to consider the next questions:

• Where to place Tacker orchestrators among several competing domains? Also taking into account affordable expenses on service provisioning; appropriate latency bounds between Tacker instances and members of domains (that are carried by Tacker instances).

• How to cope with orphan OpenStack clusters after Tacker server failure? Appropriate REassociation of orphane clusters between the Tacker instances in case of one particular instance failure?.

TACKER RESILIENCY

We believe that a methodology that is used to achieve the goal (to provide high availability for Tacker) should consists of several elements.

Firstly we have to consider the development of the model that represents the network environment. Further we have to consider the simulation of recovery of orphan switches by alive Tacker instances that requires existing model of network environment.

Secondly we should consider an approach to assign the orphan switches to the Tacker servers. The approach is presented in subsection 5.

Thirdly the position of Tacker nodes should be taken into account to place 1, N Tacker servers in the simulated network environment such that they could benefit from their location (for example low latency's from Tacker severs to OpenStack clusters).

Model of network environment

Normally in the network environment we have different distances (thus latency's) between OpenStack clusters, Tacker servers.. Therefore it is more desirable to connect new OpenStack clusters to existing OpenStack clusters such that delays (thus distances) to the Tacker server would be minimal. Also it was important to have Tacker servers not connected directly between each other, but through the intermediate nodes, as domains (set of OpenStack Clusters) could be located in different places. Consequently the last requirement to the model is to organize Tacker Servers and OpenStack clusters into domains.

Edges (links) that connect different domains have higher latency than edges which connect OpenStack clusters and Tacker inside domains. In current model for the sake of simplicity all off the links between domains have weight/latency = 20. If create the model without links with weight = 20, each domain is represented by tree [START_REF] Comer | The Ubiquitous B-Tree[END_REF], where the head of the tree is an orchestrator.

Due to the model features OpenStack clusters linked to the node of the domain d ∈ D which has the minimum distance between such OpenStack cluster and node from domain d ∈ D. As each node in domain has some limited degree, i.e. while the number of neighbours of each node is less than limit, the OpenStack Clusters want to be connected to the head of the tree to have shortest distance.

After analysis of existing models in section 2 and network topology with Tacker servers, we discovered that it is required to create a different model. Our model is inspired from tree model as it could represent the separated domains of the network. (Above we described that each domain itself without links with weight = 20 is a tree).

The difference of our model from tree model that we propose tree relaxed model, which has a feature similar to Caveman relaxed graph model [START_REF] Judd | Behavioral Conflict and Fairness in Social Networks[END_REF]. When graph is formed from each community randomly chose one edge. The next step is delete this edge that connects nodes inside of domain and add it in the way such that it would connect 2 different domains. In the presented model, we skip the deletion of edge that connects the vertices inside of domain and add the new edge with weight = 20 between domains.

An example of the model that presents random topology (nodes and weights) with some constrains (parameters) is presented in the figure 1. Where initial parameters are:

• N all -total number of nodes in the network topology (32 elements in the figure);

• N t -number of Tacker servers (number of domains is 3);

• N degree -number of neighbours per node (limit is 11).

Finally our model allows to run Tacker services placement algorithms and assign orphan switches algorithms on the generated graphs.

Recovery method

This subsection presents an approach to restore control on orphan OpenStack clusters by acting Tacker instances in case of one particular instance failure. The way to perform reassignment between the Tacker instances is presented below:

1) Each Tacker instances carries a set of clusters.

2) When particular Tacker instance become offline: a) Tackers are informed that there exists orphan switches. b) Tacker instances compare to which Tacker server from particular OpenStack cluster the distance is smaller. c) Orphan OpenStack cluster is assigned to the Tacker with shortest distance.

Simulation of algorithm to assign the orphan switches to the Tacker Domains was performed in Python language and presented in the figure 1. For example when we have 3 Tacker in the network. 1 Tacker went down. 2 other servers, each after another assign orphaned nodes to their domain.

Tacker placement

Firstly we want to introduce latency metric in networks. For a network represented by a graph G(V, E). Latency represented by distance from source s ∈ V to the Tacker t ∈ Vd(s, t).

We distinguish several types of distance bounds from OpenStack clusters to the Tacker server:

• maximum distance bound;

• average distance bound.

Maximum distance bound

The goal of current problem is to select suitable placement p ∈ P lacements in such way to minimize the maximal distance from OpenStack cluster s to the Tacker server t.

d max (p) = min p∈P lacement max s∈V d(s, t)
Average distance bound This problem relates to the difficulty of choosing k ∈ K ∈ V , where k is a set of Tacker servers that should be placed from the pool K of Tacker servers and V is a subset of all vertexes, c ∈ C ∈ V are OpenStack clusters, n -number of vertex in the graph G(E, V ). Where the average distance from Tacker servers to OpenStack clusters represented by the formula:

d average (p) = 1 n s∈V min p∈P lacements d(k, t)
Lower value of d average gives lower latency between OpenStack clusters and Tacker servers.

Placement algorithm

In practice, the most common approach to solve the type of such problems, is described in method that is called "Lloyd's Algorithm" [START_REF] Secord | Weighted Voronoi Stippling[END_REF]. This method could give an optimal solution with NP-hard complexity. On each iteration the the center of mass of the k domains that are represented by k Tackers is a solution. These centers are used for a re-clustering. The iterations repeat until the algorithm is converged.

We combine "Voronoi iteration" solution [START_REF] Secord | Weighted Voronoi Stippling[END_REF] with the Relative Point Centrality (RPC). Where: 1) "Voronoi iteration" -is a partitioning of a plane into regions based on distance to points in a specific subset of the plane.

2) relative point centrality was suggested by

Beauchamp in [START_REF] Linton | Centrality in Social Networks Conceptual Clarification[END_REF].

The combination of "Voronoi iteration" solution and RPC is presented as the pseudo code below. Choose random vertex Z The placement algorithm proposed in this section works as follows: we repeat algorithm until Tacker placement does not change. Firstly we randomly choose initial Tacker placement. Secondly in main loop we assign not assigned vertexes to the domains that are located closer based on the distance. Distance is shortest path that consider edge values. After we calculate the RCP for all nodes in domain, and choose one node per domain with the highest RPC as a Tacker Server. Thirdly the Tacker placements obtained in the end of algorithm we will use in next iteration of main loop as initial Tacker placement, if the placement is different from previous.

Algorithm 1 Placement algorithm

The RPC is used to determine suitable placement for Tacker server based on the distance to OpenStack Clusters. As each subgraph G is connected graph, the relative point centrality of each node is a measure of centrality of the node in the G based on distances to the neighbours in subgraph G. For vertice u the relative point centrality is the sum of the shortest path distances from u to all n -1 other vertices. Thus higher RPC node has, the more central it is and the closer it is to all other nodes. Therefore higher RPC -the better result. The relative point centrality formula is presented below: On each iteration the new tacker instance is choose by algorithm according to RPC. Where RPC monotonically increases (or not changes) after each iteration. When the tacker placement set does not change anymore the algorithm is converged to a local optimum as it has similar behaviour to algo that is presented in K-Means Clustering documentation [START_REF]K-Means Clustering Documentation[END_REF].

Relative point centrality(u) = n -1 n-1 v=1 d(v,
To converge to a global optimum the algorithm should be iteratively executed with estimation of RPC up to infinity number of times. On each turn the average RPC between the domains should be recorded. Within next iteration initial placement of Tacker servers changes or not based on the RPC and with more iterations we have gain a better result that could be expressed by formula below. t 0 , ..tN -tacker placements. P ← t 0 , ..t N inDomains. w -current candidate configuration based on P. w -new candidate configuration based on P .

(ω-> ω ) = (ω) if RPC(ω ) -RPC(ω) <= 0 (ω ) else

DISCUSSION

The simulations performed with using of our network model show that proposed method i.e. the combination of proposed network model, placement and recovery algorithms answer the questions raised in section 4. With the proposed method we may provide Tacker resiliency when the real network topology is similar to our network model. The method includes definition of where to place Tacker orchestrators among several domains and provide the method that defines the behaviour of system in case of Tacker server failure. Also we present in figure 2 the performance of placement algorithm. The figure 2a shows how frequently algorithm gives the result that is expected given different degree and number of nodes. It is supposed that more that few OpenStack Clusters normally connected to a Tacker Servers (Tacker servers have big degree), so we could see that figure 1a shows that with bigger degree algorithm give better ratio of expected/obtained result. The figure 2b shows how many iterations (0-25) of the algorithm should be performed to satisfy the average latency between Tacker servers and OpenStack clusters. We repeatedly execute the placement algorithm and we stop when the latency bound (30, 40, 50, 60 units) is achieved. The figure 2c shows the dependency of number of nodes in network topology, average latency between Tacker servers and OpenStack clusters and execution time of the algorithm. The experiments were performed with different number of node degree and figure 2c shows that bigger degree 5 versus 50 does not significantly affect the performance that means that with using of this algorithm we are not obliged to consider constrains such as "How many OpenStack clusters we may connect to the Tacker" or "How many OpenStack clusters we may connect to other OpenStack clusters".

CONCLUSION

This paper presents a methodology to achieve Tacker resiliency. As an outcome we present a combination of elements that together give a method to achieve Tacker resiliency. More precisely the way to recover orphan switches and the way of proper placement of the Tacker servers were proposed. With our model of network topology the simulations of placement and recovery algorithms were performed. Presented results of proposed placement algorithm are in agreement with expectations and further it is be possible to improve the speed of the algorithm. At present, the placement algorithm is a NP-hard and it's basic principle is focusing on the core of domains and by calculating Relative Point Centrality for each vertex -the sum of all shortest path distances from u to all n-1 other vertices. A further work direction, due to the simplicity of ideas, could be their application in the network analysis not only to solve the problems related to Tacker resiliency.

Fig. 1 :

 1 Fig. 1: Demo of Recovery method is presented, where nodes 0,1,2 are Tacker instances thus there are 3 domains and all other nodes are OpenStack clusters.

foreach Y ∈ d do 11 :←

 11 Find the vertex Y with dist(Z, Y ) shortest 12: if dist(Z, Y ) < shortestP ath then 13: shortestP ath = dist(Z, Y ) foreach d i ∈ Domains do 19:Add nodes from d i to the subgraph of G i Choose vertex with RP Cmax as tacker g

Algorithm 2 6 : 8 :

 268 u)The sum of d(v,u) depends on the number of nodes in the G therefore closeness is normalized by the sum of minimum possible d(v,u). Relative Point Centrality (RPC)1: G ← V, E 2: foreach v ∈ V do 3:sp length ← ShortestP ath(v, other vertexes) 4: sp sum = sp length 5: # normalization to N = all nodes -1 RP C[u] = (len(sp length ) -1.0)/sp sum 7: s = (len(sp length ) -1.0)/len(G) -1 RP C[u] * == s 9: end. Algorithm performance analysis is presented in the figure 2

Fig. 2 :

 2 Fig. 2: Performance analysis of placement algorithm

1 :

 1 N ← number of domains 2: G ← V, E 3: p initial ← t 0 , ..t N inDomains 4: while p curr ! = p prev do

	5:	while Not Assigned Vertexes is not ∅
	6:	ShortestP ath ← ∞
	7:


Jens Gustedt

Maxime Mogé

Memory access classification for vertical task parallelism

Keywords: automatic parallelization, iterative algorithms, parallel control flow, memory access classification

We present a paradigm and implementation of a parallel control flow model for algorithmic patterns of two nested loops; an outer iteration loop and an inner data traversal loop. It is centered around memory access patterns. Other than dataflow programming it emphasizes on upholding the sequential modification order of each data object. As a consequence the visible side effects on any object can be guaranteed to be identical to a sequential execution. Thus the set of optimizations that are performed are compatible with C's abstract state machine and compilers could perform them, in principle, automatically and unobserved. We present two separate implementations of this model. The first in C++ uses overloading of the operator[] to instrument the memory accesses. The second in Modular C uses annotations and code transformations for the two nested loops. Thereby the code inside the loops may stay as close as possible to the original code such that optimization of that code is not impacted unnecessarily. These implementations show promising results for appropriate benchmarks from polybench and rodinia.

I. INTRODUCTION AND OVERVIEW

Race conditions are one of the main challenges of automatic parallelization, and over last decades a lot of effort has been put into understanding and mastering them, in particular into their detection in faulty executions, see e.g [START_REF] Banerjee | A theory of data race detection[END_REF]- [START_REF] Wang | Data race detection based on dependence analysis[END_REF]. Another direction of research has been to attempt to provide a race free task schedule once that a dependency graph is given [START_REF] Sánchez | On efficient distributed deadlock avoidance for realtime and embedded systems[END_REF]. Unfortunately, the need to manually specify a dependency graph can be a burden for the programmer and also defeats automatic parallelization.

In our opinion, to be closer to current programming practice, detection of dependencies and thus possible parallelism should be based on the features that are present in every day's imperative programming languages. In these, dependency between program statements are usually not explicit but they are forced indirectly via data accesses when one statement that reads a data has to be scheduled after another statement that wrote it.

Therefore parallelization of loops with irregular access patterns requires a runtime data dependency analysis that is followed by a tool that ensures that dependencies are enforced. There are currently mainly two types of tools: speculative and scheduling based.

Speculative techniques such as LRPD [START_REF] Dang | The r-lrpd test: speculative parallelization of partially parallel loops[END_REF] identify the dependencies during the execution of the loop. The loop is run in parallel as a doall and followed by a test to check if the computations are correct. If not, a rollback to the previous correct state is performed and a sequential execution of the loop is carried out. For certain types of data accesses, namely affine or nearly affine functions over loop counts, polyhedral models can be used to provide race-free static task parallelism [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF], and these approaches have also recently be used successfully to build runtime systems that are able to rollback execution if access violations occurred [START_REF] Martinez Caamaño | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF]. For programs with low effective parallelism, this can lead to a significant overhead due to rollback followed by sequential execution when many conflicts occur between threads.

Scheduling based methods preprocess the loop to compute a dependency graph at runtime and use this information to do a parallel scheduling. The Inspector-Executor model [START_REF] Saltz | Run-time parallelization and scheduling of loops[END_REF], [START_REF] Venkat | Automating wavefront parallelization for sparse matrix computations[END_REF] consists of three phases: inspection, scheduling and execution. First the program is instrumented to explicitly compute a dependency graph. Then a parallel scheduling of the iterations is derived, and then the iterations are run in parallel in wavefronts, using synchronization between consecutive wavefronts. August et al. [START_REF] August | Automatically exploiting cross-invocation parallelism using runtime information[END_REF] propose a scheduling based method that does not explicitly compute the dependency graph, but instead overlaps the inspection and scheduling. A dedicated scheduling thread dynamically ensures that there is no conflict between threads, and allows cross invocation parallelism. However this method still needs an inspector that is run at the beginning of each iteration. The limitations of these methods are inherent to the Inspector-Executor model: it needs one inspector per iteration. In the case of a cyclic dependency between data and computation of shared addresses the inspector is basically the whole loop body.

In contrast to that, we will not restrict ourselves to a specific fine grained access projection model, but instead focus on certain type of applications, namely those that repeat the same data access pattern (coined data traversal) over a set of iterations, Section II. This is a pattern that is found widely in the field, e.g many iterative algorithms that traverse matrices or geometric objects fall into that category. Using the assumption of a constant data access pattern, we detect dependencies and derive an implicit scheduling at runtime using Ordered Read-Write Locks (ORWL) during some initial iterations, see [START_REF] Clauss | Iterative Computations with Ordered Read-Write Locks[END_REF]. After this instrumentation phase, we do not need to check any condition or recompute dependencies, thus eliminating the drawbacks of the Inspector-Executor model.

As a basis for our argumentation, here, we will use C's model of side effects on data [START_REF]Programming languages -C, cor. 1[END_REF]. C has an abstract state machine that can be used to describe the effects of any valid program and compilers are allowed to perform all optimizations that respect that machine model. Our present work emphasizes on the fact that our parallelization fully respects the computation in that machine model. A fine grained (theoretical) model is presented that guarantees that all computational results are exactly as they would have been produced by the originating sequential program, Section III. Based on [START_REF] Clauss | Iterative Computations with Ordered Read-Write Locks[END_REF] we are able to prove that our parallel execution model is correct, fair and deadlock free.

Being much too fine-grained for practical utilization, we have to coarsen our model by grouping programming steps into "meta-steps" and by classifying objects into "meta-objects", Section IV. We are able to prove that by doing so the good properties of our model are maintained, and present two different strategies that can be used for a memory access classification.

Our approach has been implemented twice, Section V; first by using C++ and its ability to overload the operator[]. By that we are able to dynamically instrument the access pattern of complex code without having to rely on the programmer. A second implementation in Modular C [START_REF] Gustedt | Modular C[END_REF] provides #pragma annotations for a code transformation that shows to be more efficient, but that still needs manual annotations of the data accesses. Both implementations are tested with a set of benchmarks from the polybench and the rodinia collections, Section VI, and show very satisfying speedups that are close to the maximum that we could expect.

II. ITERATIONS OVER DATA TRANSVERSALS

In this work we restrict ourselves to a specific framework, namely programs that are dominated by two levels of loops.

1) An outer loop that we call iteration, e.g, an iteration over a simulated time or an iterative approximation of some numerical quantity. 2) One or several inner loops that we call data traversal tasks or just tasks, e.g loops that visit all elements of a matrix or that iterate over a geometric domain such as the facets of a polyhedral object description. Our main assumption for data accesses is the following:

The data access pattern does not depend on the iteration.

That is, we assume that each outer iteration visits the same data in the same order as the previous ones, e.g a traversal of a matrix or of geometry elements would be done in the same order in each iteration. That does not mean that a task can't use the iteration variable for its computation, only that its data access should not depend on it.

The inner part of a data traversal task t is called a step, denoted by T t (i, p), characterized by the triple (i, t, p), where i is the actual iteration, t is the task and p is the actual position in the traversal. That is, we assume that the program identifies code segments that perform these steps, and that the information about the current iteration, task and position in the traversal is available.

Steps are sequenced, that is there is a linear order in which the steps are executed by the program. This order corresponds to the lexicographic ordering of the indices (i, t, p).

Listing 1 illustrates what we have in mind here. It shows an idealized code excerpt of rodinia's hotspot3D benchmark, where we use a split directive (pragma) and a special task ∶∶ for construct to identify a data traversal task. A step is then the entire code inside the {}. It has access to the iteration variable i and the traversal position z.

Another feature of that implementation is also shown, namely data accesses are instrumented explicitly by means of directives that proclaim that the step may modify (lvalue) or just read (rvalue) a specific data. Evidently, such data access annotations could be avoided with sufficient compiler support that would track access to data. But as this implementation entirely works in source and C does not allow for operator overloading, this must be left to future enhancements of the platform.

A. Execution model: loop exchange and parallelism

In general, we suppose that a sequential program P has the following form:

for (all iterations i) for (all tasks t) for (all positions p in T t data traversal) T t (i, p);

Our goal is to improve the performance of P by parallelizing it. The main idea is to automatically perform a loop exchange. After the exchange the outer loops are in fact over the positions p of the data traversals and (at least conceptually) an independent thread is launched for each such position. Inside such a thread, the iterations are performed over the individual steps that are all associated to the same position p.

Data races could occur when two inner steps would access the same data concurrently. Therefore we will have to develop a mechanism that helps us to detect when one step has to wait for another one because of such a concurrent access. As already implicit in the hotspot3D example, our notion of dependency of steps is not an abstract model of data or control flow, but directly deduced from access to objects of the program. It directly follows C's notion of side effects [START_REF]Programming languages -C, cor. 1[END_REF], that is, of modifications effected to data.

Step S 1 directly depends on step S 0 if

• S 1 is sequenced after S 0 and • S 1 reads or writes a data object o that S 0 has written. By that, direct dependency between two steps can never be induced by only read accesses to a common data.

We transitively close the direct dependency relation to a partial ordering ℘ the dependency relation on the steps that are executed by a program. Two steps S 0 and S 1 are independent if neither of them depends on the other and if so, they can be executed concurrently without introducing a race condition. Observe that two independent steps may be from the same or different tasks, iterations and traversal positions.

The notions of dependence and independence easily extend to sets of steps 0 and S 1 .

B. Parallelizations

Fig. 1 shows an example of the initial section of a direct dependency graph and different scheduling strategies that have been applied. In all four examples, lines correspond to concurrent executions of steps at a given time. Columns visualize all steps that deal with the same position in their respective data traversal. For the sake of the example, we make the simplified assumption that all steps have the same execution time. Sequential execution. A sequential scheduling of the graph is shown in Fig. 1(a). Here, all steps are strictly executed in iteration and traversal order. The rectangle in the figure shows the pattern of the steady state, that is the pattern of the graph that is repeated over and over again until the final iteration. Horizontal parallel execution is similar to what OpenMP would do with a for loop that is prefixed with a parallel for directive. Only here we will guarantee that data dependencies between steps of the same task will be honored. Fig. 1(b) shows that our example only allows for a limited horizontal parallelization that respects dependencies: the small rectangles inside the steady state visualize the steps from the same task that can be parallelized. Vertical parallel execution is illustrated in Fig. 1(c). Here, we start the steps as early as possible under the constraint to still execute the steps of the same task and iteration in order. Vertical and horizontal parallel execution as a combination of both parallel modes, shows the best available parallelism, see Fig. 1(d). Here it not only allows to run steps of the same iteration in parallel (e.g T 1 (0, 0) and T 1 (0, 1)) but also steps of the same task but from different iterations, e.g T 0 (0, 3) and T 0 (1, 0). Or stated differently, a new iteration of a task may already start before the previous has completely finished.

III. A FINE GRAINED EXECUTION MODEL

We aim to develop an execution model that is compatible with C's abstract machine in a way that the parallelizations that we propose do not change the observable state of a program.

Observable changes in the states of that abstract machine occur through so-called side effects. By neglecting IO and similar external events for the moment, the interesting side effects in our context are modifications of the data objects that a program manipulates. Such manipulations occur through operations (such as assignment = or increment ++) or C library functions (such as memcpy).

For each object o in an execution of a program the C standard requires that accesses to o are properly sequenced, that is, that it can be deduced from the program structure if

• a write operation A on o provides the value for a read operation B of o;

• two write operation A and B on o are not separated by another write operation. Read-only operations in turn are not necessarily expected to be properly sequenced. E.g the two evaluations of i in a function call f (i, i) may be executed in any order.

For our context here, we can assume that an access A to object o is characterized by a triple (i, t, p) for step T t (i, p), and, if it is a write access, by the new value that is stored in o at the end of step T t (i, p). In particular, we don't distinguish if a step T t (i, p) does several accesses to the same object o and whether read and write access occur inside the same step.

The modification order M (o) = M 0 , M 1 , . . . is the precise list of changes that are applied to an object o. The access order L(o) is M (o) together with the read accesses to o, that is where we insert the read accesses in sequence order between write accesses:

L(o) = W 0 , R 0 0 , R 1 0 , . . . , W 1 , R 0 0 , R 1 0 , .
. . Since we are assuming that memory accesses do not depend on the iteration, L(o) is of the form

L(o) = L(o, 0), L(o, 1), . . .
where L(o, x) is a predetermined list (x, t 0 , p 0 ), (x, t 1 , p 1 ), . . . where we textually substitute x by i for each iteration i.

By changing the perspective from data objects to steps, for each step T t (i, p) we can now identify the objects o = o 0 , o 1 , . . . that it accesses as well as the list position (i, t, p, o) in L(o, i) where this access occurs.

A. Fifo execution

Such data access lists can now be used to establish a parallel execution model for the program steps, that takes care of the coherence requirements. For each object o we establish a FIFO data structure F(o) that follows L(o) and for which we define time 0 T 0 (0, 0)

5 5 1 T 0 (0, 1) Ø Ø 5 5 2 T 0 (0, 2) 5 5 3 T 0 (0, 3) Ø Ø 4 T 1 (0, 0) 5 T 1 (0, 1) Ø Ø 5 5 6 T 1 (0, 2) 7 T 1 (0, 3) Ø Ø 8 T 0 (1, 0) 5 5 9 T 0 (1, 1) Ø Ø 5 5 10 T 0 (1, 2)
5 5

11 T 0 (1, 3) Ø Ø 12 T 1 (1, 0) 13 T 1 (1, 1) Ù Ù 5 5 14 T 1 (1, 2) 15 T 1 (1, 3) (a) sequential execution time 0 T 0 (0, 0) 5 5 1 T 0 (0, 1) Ø Ø 5 5 2 T 0 (0, 2) 5 5 3 T 0 (0, 3) Õ Õ 4 T 1 (0, 0) T 1 (0, 1) Õ Õ 5 5 5 T 1 (0, 2) T 1 (0, 3) Ø Ø 6 T 0 (1, 0) 5 5 7 T 0 (1, 1) Ø Ø 5 5 8 T 0 (1, 2)
5 5 time 0 T 0 (0, 0)

9 T 0 (1, 3) Õ Õ 10 T 1 (1, 0) T 1 (1, 1) Ö Ö 5 5 11 T 1 (1, 2) T 1 (1, 3) (b
5 5 1 T 0 (0, 1) { { 5 5 2 T 1 (0, 0)
T 0 (0, 2)

5 5 3 T 1 (0, 1) { { 5 5
T 0 (0, 3)

{ { 4 T 0 (1, 0) 5 5
T 1 (0, 2)

5 T 0 (1, 1) { { 5 5 T 1 (0, 3) { { 6 T 1 (1, 0) T 0 (1, 2) 5 5 7 T 1 (1, 1) 1 1
T 0 (1, 3) (c) vertical parallelism: task starts executing as soon conflicts with preceding tasks are resolved. There are 2 parallel threads in steady state. time 0 T 0 (0, 0)

5 5 1 T 0 (0, 1) { { 5 5
2 T 1 (0, 0)

T 1 (0, 1) { { % % T 0 (0, 2)
5 5

3 T 0 (1, 0) 5 5
T 0 (0, 3)

{ { 4 T 0 (1, 1)
{ { 5 5

T 1 (0, 2)

T 1 (0, 3) { { 5 T 1 (1, 0) T 1 (1, 1) 1 1 T 0 (1, 2) 1 1
(d) vertical and horizontal parallelism: steps are scheduled with first-fit. On average, 8 3 = 2.6 parallel threads are in steady state. init, release and acquire operations. They ensure that any write access is only granted when all previous accesses in L(o) have been released, and that all read-only accesses that follow the same write access can be honored concurrently. Now we can augment each step of the program such that it is protected from race conditions by acquiring all FIFOs for all accessed objects and by releasing them thereafter. A protected step Prot⟪T t ⟫(i, p) (1) acquires all FIFOs of accessed objects (2) executes T t (i, p) and (3) releases them thereafter.

Clearly, substituting Prot⟪T t ⟫ for T t for all t in the whole program does not change its visible side effects: an acquire operation would only be issued after all previous positions in L(o) are released, and all operations T t themselves would still appear in sequence order.

We can now proceed to our final goal, the loop exchange. A step iteration StepIt⟪T t ⟫ performs the iteration loop

StepIt⟪T t ⟫(p) { for (i in iteration) Prot⟪T t ⟫(i,p); } Now our parallelized program can be formulated as for (all tasks t) do parallel for (all positions p in T t data traversal) do parallel StepIt⟪T t ⟫(p);

Theorem 1: The parallelized program is well defined and produces the same side effects as the original program.

A proof of this theorem follows from the fact that the dependency relation ℘ is acyclic and from the theory developed for Ordered Read-Write Locks (ORWL), see [START_REF] Clauss | Iterative Computations with Ordered Read-Write Locks[END_REF]. This formalizes a model for autonomous execution of tasks or processes that manage their shared resources by means of cyclic FIFOs that are attached to data resources. Here, the ORWL model can guarantee correctness, equity and liveness of execution.

IV. COARSENING Generally for the framework that we described so far, we are unfortunately not in a coarse grained setting.

• Data traversals can have many steps, exceeding the possible number of threads for current architectures.

• Prot⟪T t ⟫ introduces runtime overhead for each memory access.

• Additionally, we also would have to maintain a FIFO for each data object.

• There is an important overhead to initiate FIFOs for all objects. To be useful, we have to coarsen our model with respect to two aspects. First, we have to ensure that the number of steps (and thus the number of threads) does not exceed the limits of existing platforms. Second, we have to avoid that each data object gives rise to its own FIFO.

A. Groups of steps

To reduce the number of steps, we transform the original program by grouping steps: for (i in iteration) for (t in tasks) for (group g of positions in T t traversal) for (p in g) T t (iteration, p);

The easiest strategy to form such groups of positions is to divide the positions p 0 , p 1 , . . . , p r in the traversal into a fixed amount m of intervals of successive positions. This allows to rewrite a group g of steps into a "meta-step" MetaStep m ⟪T t ⟫:

for (p in g) T t (i, p);

And thus we can reformulate the program as

for (i in iteration) for (t in tasks) for (g of meta-positions in MetaStepm⟪T t ⟫ traversal) MetaStepm⟪T t ⟫(i, g);

As a consequence this rewrite of the program fits exactly the same model as the original program; it consists of an iteration loop that is nested with data traversal loops such that the data access pattern does not change with the iterations. The induced dependency relation ℘ m between groups g 1 and g 2 is given by the union of all (p 1 , p 2 ) ∈ g 1 ×g 2 such p 2 reads data written by p 1 . Since by definition of the groups all p 1 ∈ g 1 come before all p 2 ∈ g 2 such a union never induces a cycle.

Lemma 1:

The induced dependency relation ℘ m between meta-steps is acyclic.

B. Classes of objects

To bundle the management of FIFOs for several data objects at a time we partition the set of data elements D into m ′ meta-objects C = {c k ⊂ D}. Such a classification of objects into meta-objects gives rise to an induced dependency relation ℘ m C between meta-steps in a obvious way. We have Lemma 2: ℘ m C is an partial order extension of ℘ m . Corollary 1: ℘ m C is an acyclic, irreflexive and transitive.

We then attribute one FIFO F(c) to each meta-object c ∈ C and define the protected meta-step analogously as before. Clearly, the resulting scheduling is more restrictive, but, because of Corollary 1 its validity is not impacted.

Theorem 2:

The parallelized program that is coarsened with respect to a step grouping of m and an object classifier C is well defined and produces the same side effects as the original program.

To assemble the meta-objects we investigated two strategies. a) Classify groups of contiguous memory accesses in ranges.: Suppose the main data access is through a vector V[j] for j = 0, ..., N -1 and that the access to that vector is regular such that the progression of the data traversal corresponds to a linear access to the vector. To classify the access to V we can then fix a chunk size M and define the metaobjects Contig(j) as the j th chunk of size M. The advantage of Contig(.) is that the membership for any data object to its class is easily computed. Therefore the overhead that is needed to establish the FIFOs for the meta-steps is relatively small. b) Classify random memory access by ownership.: Vector classification does not work if data access is through lists, e.g if we do a data traversal over a complicated data structure or if our access pattern follows some random order. The class Owner(g) is the set of objects for which g is the first writer in the sequence.

We may prepare the computation of Owner(g) by first running several sequential iterations: first we have to observe the first writer for each object, and then we have to classify all other accesses with respect to that writer. After this initial setup of the meta-objects, the remaining iterations can be done in parallel.

Compared to Contig, Owner classification has a startup that is is more expensive. On the other hand, no assumption is made about access order or about a specific organization of the data.

V. IMPLEMENTATIONS

We have implemented the Owner classifier with two different approaches. The first uses C++ and overloading of the operator[] and the second uses Modular C [START_REF] Gustedt | Modular C[END_REF] and #pragma directives for a source-to-source transformation.

Both build on the EiLck library1 for the FIFO functionality. EiLck is a standalone library written in Modular C for Linux. It implements cyclic access FIFOs on the basis of Linux' futex system call that efficiently combines atomic access with a scheduling queue in the system kernel, [START_REF] Franke | Fuss, futexes and furwocks: Fast userlevel locking in Linux[END_REF]. All scheduling of tasks is done implicitly via these FIFOs; EiLck executes the application without dedicated management threads. EiLck is compatible with different thread libraries (C11, C++11, POSIX, OpenMP) and is interfaced to C and It efficiently avoids overhead for busy waiting if a FIFO access is congested and provides a minimal run time overhead in case it is not.

A. C++: Operator overloading and explicit thread creation

Besides the use of Eilck, our C++ implementation is based on standard tools, mainly operator[] and the explicit creation of threads for the execution of the tasks. operator[] is used to instrument every single memory access to a set of vectors that the programmer designates as being shared between steps. In addition, the programmer also identifies the iteration loop and registers the data traversal tasks as functions that are executed in separate C++11 threads.

Then, during execution, data traversals are grouped and protected automatically. For the grouping of the steps into meta-steps, each task counts its number of calls to operator[] and uses this count as an indication of the traversal position. For the classification of the memory accesses we use a dynamic std::map to register the owners and then classify each access according to that owner. A set of FIFOs is constructed from this, and used in subsequent iterations.

This approach has the advantage that we don't even have to identify steps. Any code that supplies sufficient amount of work can automatically be split by this operator[] approach into steps and meta-steps, respectively. Listing 2 shows a complete parallelization of the hotspot3D benchmark using this implementation. The approach has several disadvantages, though:

• All meta-steps of a task are executed in the same thread.

• Control operations by operator[] interleave the computation.

• The operator[] needs a thread local variable for the position in the traversal.

B. Modular C: pragmas and source-to-source transformation

Listing 3 shows a complete parallelization of the hotspot3D benchmark using our Modular C implementation. Here we use two forms of Modular C directives. The first, insert, is for a construct that precedes some statement or block and potentially changes its interpretation. The second, amend, can be used to rewrite an entire code snippet as necessary.

Here, both directives refer to an external script, mctask, that implements our parallelization. The main features shown in that example are the following • An iterate directive replaces the outer iteration loop.

• A split directive prefixes the data traversal.

• A task ∶∶ for (or task ∶∶ while) loop specifies the data traversal.

• Data access directives lvalue and rvalue indicate all memory accesses that could be subject to race conditions. The iterate directive uses amend to transform the code of the iteration. It creates two modified copies of the code, one for an initial instrumentation phase and one for the steady state iterations. For the first phase, lvalue and rvalue are replaced by calls to access classification functions, for the second they are removed such that the steady-state iterations can use the original code directly. The split directive together with the task ∶∶ for ensures that the loop is interpreted as a single task. It groups the steps into meta-steps and ensures that the protection against race conditions is inserted between them. But in contrast to a steps directive, see below, it still only launches one thread per task, and we only will see vertical parallelism if we only use split.

The transformed parallel program can be quite effective because the inner data traversal is exactly as programmed originally. In particular, the program maintains the same optimization opportunities, e.g for vectorization. Also, compared to our C++ version, there is no need for thread local variables. On the other hand, this implementation here requires that the structure of the program makes the iteration and data traversal loops apparent.

Listing 3 also features two other directives that are be quite useful. With the innermost amend, the alternate directive implements an automatic duplication that alternates the roles of tOut and tIn. The next level of amend with a duplicate directive duplicates its inner part, again, such that at the end we have four copies of the inner part, giving rise to four different tasks.

In total we have 4 × 8 = 32 meta-steps. Here, because of the split directive, we have one thread for each task, so 4 tasks in total that are parallelized vertically. To gain also horizontal parallelism we can replace the split directive by a steps directive (not shown) with 32 threads in total.

VI. BENCHMARKS

We tested our implementations with those benchmarks from the polybench and the rodinia series that fulfill our requirements, and inside a real life application, the SOFA frameworkwww.sofa-framework.org.

The test platform is a Linux machine with an Intel Xeon E5-2650 v3 processor, 2 CPU sockets, 10 real cores (and 20 hyperthreaded cores) each. Although we did not expect much parallelism beyond a speedup factor of 2 or 3, see below, we chose such a platform with relatively many cores to avoid an artificial bound on the parallelism. For the C++ implementation, we adapted the code using C++ vectors.

Fig. 2 shows the speedups that we obtained with these relatively simple tests, once we have reached the steady state. We compare the original sequential implementation seq to our implementations using Modular C with split and steps directives, and C++. Generally, the expected parallelism is bound by the number of tasks that we can identify in the benchmark, which in the examples is some small number. But for some of the tests we unroll the iteration loop (2x or 4x factor) to increase the number of tasks for more parallelism.

With Modular C, we see that all but one of the benchmarks achieve worthy speedup factors that range from 1.2 or 1.3 (srad2D), over 1.5 (fdtd), 1.7 (adi and heat3d), 2.2 (hotspot3D), to up to 4 for seidel_2d. Indeed, the latter parallelizes perfectly with the number of duplicates of the data traversal that we produce.

Most of these parallelizations are straight forward code annotations as we have seen them for the hotspot3D benchmark above. There is no clear winner between the split and steps directives; adding more threads is not always an advantage. The parallelization of srad2D is a little bit more involved than the others, but when we divide that traversal loop into several such that we obtain 4, 5 or 7 tasks in total, we observe a speedup of about 1.3.

Only the nearest neighbor (nn) benchmark from the rodinia suite detaches from this. It is completely dominated by the IO for reading the input data.

With the C++ version, most of the tested benchmarks achieve worthy speedup, too, although lower than with Modular C. The overhead induced by the overloading and the inhibited compiler optimizations can seriously impact the performances (e.g. only 1.1 speedup for heat3d and even a slowdown for adi). However, a simple unrolling of the iteration allows us to create more tasks and gives a much better speedup (1.7 for heat3d with duplicated tasks, and up to 2.9 with 4x unrolling).

The cost of the startup steps is highly dependent on the data objects considered and the application. With default settings, the Modular C version uses fine grained data objects (cache lines) with a cost ranging from 2 to 120 sequential iterations. This granularity can be manually adapted to reduce this cost. The C++ version, uses coarser data objects (rows of 2D matrices) with a cost of 3 sequential iterations for seidel2D, 12 for heat3d, 16 for hotspot3D.

We also tested our C++ implementation on a large scale application, the SOFA framework, with interesting performance: we were able to maintain the runtime overhead induced by our parallelization mechanism low enough to get a significant speedup on a specific test simulation.

VII. CONCLUSION AND OUTLOOK

We have shown that a semi-automatic parallelization tool for a special class of iterative algorithm can be implemented. Provided that the data access pattern does not change in the iteration loop, execution is guaranteed to follow the sequential modification order for each data. Thereby we can e.g guarantee that numerical iterations lead to exactly the same results as sequential execution and that all proofs for convergence translate directly to the parallel version. We showed that our implementations leads to satisfying speedups for a variety of benchmarks.

For the moment our implementations both need intervention from the programmer. Both need an identification of the iterations and data traversals. Also, they are limited in their ability to instrument the memory access. The C++ version has a certain runtime overhead for the overloaded operator[]. The Modular C version needs manual annotation of the accesses and has a high startup cost because of its use of fine grained data objects. We think that all these shortcomings can be overcome by implementing the memory access instrumentation directly in the compiler.

The need of identification could be circumvented by using a framework such as Apollo [START_REF] Martinez Caamaño | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF] that allows to speculate on a particular pattern, and then to rollback if the expected pattern was not respected.

Fig. 1 .

 1 Fig. 1. Schedulings of the direct dependency graphs of steps. Arrows represent data dependencies, numbers on the left represent execution times of the steps in their line. Rectangles illustrate patterns of the steady state.

  Listing 2 PARALLELIZATION OF RODINIA'S hotspot3D BENCHMARK USING OPERATOR OVERLOADING -THE CODE WAS SIMPLIFIED FOR CLARITY (E.G. OMISSION OF TEMPLATE PARAMETERS). // encapsulate computation in a separate function void computeTask(VectorWrapper tOut_w, VectorWrapper tIn_w, ...) { for(size_t z = 0; z < nz; z++) { size_t z0 = ... ; size_t z1 = ... ; for(size_t y = 0; y < ny; y++) { size_t y0 = ... ; size_t y1 = ... ; for(size_t x = 0; x < nx; x++) { size_t x0 = ... ; size_t x1 = ... ; tOut_w[z][y][x] = tIn_w[z][y][x]*cc + tIn_w[z][y0][x]*cn + ... ; } } } } VectorWrapper tOut_w(tOut), tIn_w(tIn); // overload operator[] using a vector wrapper IterationLoopIdentifier loopId = ...; // identify an iteration loop for (size_t i = 0; i < numiter; ++i) { beginLoopIteration(loopId); threadPool->submit(new TraversalTask(computeTask,tOut_w,tIn_w)); // start task 1 threadPool->submit(new TraversalTask(computeTask,tIn_w,tOut_w)); // start task 2 endLoopIteration(loopId); }

Fig. 2 .

 2 Fig.2. Speedups of steady state iterations achieved for benchmarks from polybench and rodinia with the Modular C and the C++ implementations.

  ) horizontal parallelism: if possible steps of the same task are executed in parallel. On average,8 6 = 1.3 parallel threads are in the steady state.

  . alternate AND duplicate UNROLL THE ITERATION LOOP FOUR TIMES SUCH THAT FOUR TASKS CAN BE GENERATED.

	Listing 3
	THE COMPLETE PARALELLIZATION OF RODINIA'S hotspot3D
	#pragma CMOD insert mctask groups = 8
	#pragma CMOD amend mctask iterate = numiter
	/ * unroll for more parallelism * /
	#pragma CMOD amend mctask duplicate = 2
	/ * unroll to alternate buffers * /
	#pragma CMOD amend mctask alternate tOut tIn
	#pragma CMOD insert mctask split
	task ∶∶ for(size_t z = 0; z < nz; z++) {
	size_t z0 = ... ; size_t z1 = ... ;
	/ * insert some miraculous data coherence enforcement * /
	# pragma CMOD insert mctask lvalue = tOut[z]
	# pragma CMOD insert mctask rvalue = \
	tIn[z], tIn[z0], tIn[z1]
	for(size_t y = 0; y < ny; y++) {
	size_t y0 = ... ; size_t y1 = ... ;
	for(size_t x = 0; x < nx; x++) {
	size_t x0 = ... ; size_t x1 = ... ;
	/ * do the computation * /
	tOut[z][y][x] = tIn[z][y][x]*cc
	+ tIn[z][y0][x]*cn + ... ;
	}
	}
	}
	#pragma CMOD done
	#pragma CMOD done
	#pragma CMOD done

BENCHMARK

http://cmod.gforge.inria.fr/eilck.html