Pierre-Louis Giscard

Paul Rochet
email: paul.rochet@univ-nantes.fr

Enumerating simple paths from connected induced subgraphs

Keywords: Directed graph, self-avoiding walks, simple cycles, Hamiltonian paths, dominating sets, trace monoid, labeled adjacency matrix, MSC: 05C22, 05C30, 05C38, 05E99

We present an exact formula for the ordinary generating series of the simple paths, also called self-avoiding walks, between any two vertices of a graph. Our formula involves a sum over the connected induced subgraphs and remains valid on weighted and directed graphs. An intermediary result of our approach provides a unifying vision that conciliates several concepts used in the literature for counting simple paths. We obtain a relation linking the Hamiltonian paths and cycles of a graph to its connected dominating sets.

Introduction

Counting simple paths, that is trajectories on a graph that do not visit any vertex more than once, is a problem of fundamental importance in enumerative combinatorics [START_REF] Madras | The self-avoiding walk[END_REF] with numerous applications, e.g. in sociology [START_REF] Ross | On the determination of redundancies in sociometric chains[END_REF][START_REF] Cartwright | Structural balance: a generalization of heider's theory[END_REF]. Several exact formulas for counting simple paths have been discovered over the last 60 years. They make use of a variety of tools such as the inclusion-exclusion principle [START_REF] Bax | Inclusion and exclusion algorithm for the hamiltonian path problem[END_REF][START_REF] Björklund | Determinant sums for undirected hamiltonicity[END_REF][START_REF] Richard | Dynamic programming meets the principle of inclusion and exclusion[END_REF], recursive expressions of the adjacency matrix [START_REF] Khomenko | Identifying certain types of parts of a graph and computing their number[END_REF], sieves [START_REF] Bax | A finite-difference sieve to count paths and cycles by length[END_REF] or immanantal equations [START_REF] Gordon | The number of n-cycles in a graph[END_REF].

Maybe one of the most elegant result of this sort comes from a recent unpublished work by Perepechko and Voropaev [START_REF] Perepechko | The Number of Fixed Length Cycles in an Undirected Graph. Explicit Formulas in Case of Small Lengths[END_REF], who obtained an exact formula for the number of simple cycles of a given length in the graph, involving a sum over all induced subgraphs. 1 While the result only applies to cycles, we show that it can be extended to count simple paths with any specified end points. This leads to a unification of several known formulas for counting simple paths into a single matrix equation. In spite of its higher level of generality, our proof relies on a simple binomial expansion.

The generalization of Perepechko and Voropaev's result leads to an important simplification that makes use of the additivity of adjacency matrices of disjoint graphs. Thanks to this property, the generalized formula can be reduced to connected induced subgraphs. From a computational point of view, this reduction is particularly effective for sparse graphs where connected induced subgraphs. A remarkable consequence is an expression that links the Hamiltonian paths of a graph to its connected dominating sets.

Although the original version of Perepechko and Voropaev's result applies to undirected graphs only, all the results presented in this paper are proven for directed graphs and remain valid on weighted graphs. The proofs are surprisingly simple, owing to the use of a particular trace monoid formalism, that induces some desirable algebraic properties to paths, treated as algebraic objects [START_REF] Giscard | Algebraic Combinatorics on Trace Monoids: Extending Number Theory to Walks on Graphs[END_REF]. In this particular setting, a path is defined as a word whose letters are the oriented edges of the graph, endowed with a partially commutative structure. This formalism, developed by Cartier and Foata in 1969 in their seminal work [START_REF] Cartier | Problemes combinatoires de commutation et réarrangements[END_REF], has managed to remain inexplicably unnoticed in the graph literature, in spite of its generality and inherent simplicity.

The trace monoid setting used to define paths on the graph is described in Section 2. The matrix generating series P(z) of simple paths in a directed graph G is introduced in Section 3, along with our first result, expressing P(z) in function of the adjacency matrix of the induced subgraphs of G. We discuss some implications and how the result stands out in the literature. Finally, we present in Section 4 a simplified expression of P(z) and the resulting relation linking the Hamiltonian paths of the graph to its connected dominating sets. Proofs of the results are postponed to the Appendix.

The trace monoid structure of paths

We start by presenting the algebraic structure on paths that was introduced in [START_REF] Cartier | Problemes combinatoires de commutation et réarrangements[END_REF]. We argue that this setting is very natural as it coincides with the standard views on paths, cycles or similar typical objects commonly used in graph theory. The objective is to offer a rigorous algebraic formalism that is both general and perfectly suited to describe the objects obtained from analytical manipulations of the adjacency matrix.

Let G = (V, E) be a directed graph with vertex set V = {1, ..., N} and edge set E ⊆ V 2 , which may contain self-loops. Each edge (i, j) ∈ E of the digraph is identified with a letter ω ij and endowed with the following commutation rule: two different edges ω ij , ω i ′ j ′ commute if and only if they originate from different vertices, i.e. if i = i ′ . The resulting partially commutative monoid, or trace monoid, is defined as the sets M of finite sequences of letters ω ij , called words, obeying this commutation rule. Formally, M is the trace monoid with alphabet Σ M = {ω ij , (i, j) ∈ E} and independence relation

I M = {(ω ij , ω i ′ j ′ ) : i = i ′ }.
The multiplication on M is defined naturally as the concatenation: m.m ′ = mm ′ . Remark that this operation is not commutative in general.

In this setting, a path on the digraph G can be viewed as a word with a contiguous representation. Precisely, a path p of length ℓ ≥ 1 is a word composed of ℓ edges that can be rearranged (following to the commutation rule) so that each new edge starts where the previous ends, e.g. p = ω ii 1 ω i 1 i 2 ...ω i ℓ-1 j . Paths appear naturally through analytical transformations of the labeled adjacency matrix W, with general term W ij = ω ij if (i, j) ∈ E and W ij = 0 otherwise. In this partially commutative framework, paths of a given length k ≥ 1 are enumerated in the k-th power of W:

(W k ) ij = p: i→j ℓ(p)=k p , i, j = 1, ..., N,
where the sum runs over all paths p of length ℓ(p) = k from i to j on G. Here, the matrix multiplication uses the multiplication of words, i.e. the concatenation, which is not commutative. Thus, the trace monoid formalism provides a natural way to distinguish paths composed of the same edges but in different orders.

The standard terminology applies to paths, defined as elements of the trace monoid

M. A path p = ω ii 1 ω i 1 i 2 ...ω i ℓ-1 j
is open if its end vertices i, j are different and closed otherwise. A closed path is also called a cycle. An edge ω ij is a path of length one from i to j while self-loops ω ii and backtracks ω ij ω ji are cycles of length one and two respectively. By convention, the empty path 1 is considered a cycle of zero length.

A non-empty path is simple if it does not visit the same vertex twice, or equivalently, if all its edges commute and are different. This is another important consequence of the trace monoid structure of paths: their definition do not require specifying end vertices. For an open path, the end vertices are inherent to the contiguous representation, which is unique in this case. Indeed, one shows easily that rearranging the edges of an open path into a different contiguous sequence would require to break the commutation rule. Cycles, on the other hand, may have several contiguous representations, each one with a different starting point. For instance, the length 2 cycle ω ij ω ji from i to itself has another contiguous representation, namely ω ji ω ij , which starts from j. We refer to [START_REF] Espinasse | Relations between connected and self-avoiding hikes in labelled complete digraphs[END_REF][START_REF] Giscard | Algebraic Combinatorics on Trace Monoids: Extending Number Theory to Walks on Graphs[END_REF] for a more comprehensive description of the framework.

The generating series of simple paths

Expressions for counting simple cycles of small lengths on undirected graphs from traces of the adjacency matrix have been discovered by various authors since the 1950's [START_REF] Ross | On the determination of redundancies in sociometric chains[END_REF][START_REF] Harary | On the number of cycles in a graph[END_REF][START_REF] Alon | Finding and counting given length cycles[END_REF]. In 1972, two Ukrainian mathematicians, Khomenko and Golovko, published an exact formula for counting simple cycles of arbitrary length on any graph [START_REF] Khomenko | Identifying certain types of parts of a graph and computing their number[END_REF]. The formula results from an elaborate use of the inclusion-exclusion principle and involves a function defined through a complicated recursion. Possibly owing to its inherent difficulty and lengthy proof, Khomenko and Golovko's result has remained largely unnoticed. In 2009, the Russian mathematicians Sergey Perepechko and Anton Voropaev managed to greatly simplify the formula of Khomenko-Golovko [START_REF] Perepechko | The Number of Fixed Length Cycles in an Undirected Graph. Explicit Formulas in Case of Small Lengths[END_REF], likely by solving the difficult recursion equation involved in it. The formula obtained by Perepechko and Voropaev was only announced in a conference talk and, to the best of our knowledge, never published.

Let S be the set of non-empty subsets of V , including V itself. The cardinality of S ∈ S is denoted by |S|. For a matrix M indexed by the vertices of the graph (typically, the adjacency matrix A or the labeled version W), define the restriction M S of M to S ∈ S by

M S,ij = M ij if i, j ∈ S, 0 otherwise, i, j = 1, ..., N.
Perepechko and Voropaev's result can be stated as follows.

Theorem 3.1 [START_REF] Perepechko | The Number of Fixed Length Cycles in an Undirected Graph. Explicit Formulas in Case of Small Lengths[END_REF]. Let G = (E, V ) be an undirected graph on N vertices with adjacency matrix A. The number π ℓ of simple (unoriented) cycles of length ℓ is given by

π ℓ = 1 2 ℓ ℓ i=2 (-1) ℓ-i N -i N -ℓ S:|S|=i Tr A ℓ S , 3 ≤ ℓ ≤ N (1)
where the inner sum runs over all non-empty subsets S ⊆ V of size i.

Equation ( 1), which we shall now refer to the Perepechko-Voropaev formula, seems to have had little more impact than Khomenko-Golovko's result so far in spite of its deep elegance. This is possibly because it involves a sum over subgraphs, thus producing π ℓ in O(N ℓ log N) operations [START_REF] Perepechko | The Number of Fixed Length Cycles in an Undirected Graph. Explicit Formulas in Case of Small Lengths[END_REF]. Since the formula was devised primarily within and for a numerical computing community, its value may have been gauged from its sole computational cost.

The alternating signs in the Perepechko-Voropaev formula bears similarities with Bax's approach for counting Hamiltonian paths [START_REF] Bax | Inclusion and exclusion algorithm for the hamiltonian path problem[END_REF], based on the inclusion-exclusion principle. In fact, applying the trace operator to the matrix H n defined in [START_REF] Bax | Inclusion and exclusion algorithm for the hamiltonian path problem[END_REF] recovers exactly the Perepechko-Voropaev formula for ℓ = N. In the following, we show how these two results, among others, can be summarized into a single matrix equation. Consider the matrix P(z) whose (i, j)-entry is defined to be the ordinary generating function of simple paths from i to j, i.e.

P ij (z) = p: i→j p simple p z ℓ(p) , for z a formal variable. The information relative to simple paths and simple cycles on the digraph is entirely summarized in P(z), making it a natural object of interest. Our first result gives an expression of P(z) in function of the labeled adjacency matrices of the induced subgraphs of G. For clarity, we distinguish between open and closed paths, by writing P(z) = P cl (z) + P op (z), where -P cl (z) is the matrix generating series of simple cycles. Its i-th diagonal entry is the generating function of the simple cycles visiting i, while the entries outside the diagonal are zero.

-P op (z) is the generating matrix of open simple paths. For i = j, its (i, j)-entry is the generating function of simple paths from i to j, and its diagonal is zero.

For a square matrix M, we denote by Diag(M) the diagonal matrix obtained by setting to zero all non-diagonal entries in M. We are now in position to state our first result.

Theorem 3.2. It holds

i) P op (z) = S∈S (zW S ) |S|-1 (I -zW S ) N -|S| , ii) P cl (z) = S∈S Diag (zW S ) |S| (I -zW S ) N -|S| ,
where |S| denotes the cardinality of S.

We recover Perepechko-Voropaev's formula by attributing the value ω ij = 1 to all directed edges, thus replacing W by A, and taking the trace of P(z) in Theorem 3.2:

Tr P(z) = S∈S Tr (zA S ) |S| (I -zA S ) N -|S| = S∈S N -|S| k=0 (-1) k N -|S| k Tr (zA S ) |S|+k = S∈S N ℓ=|S| z ℓ (-1) ℓ-|S| N -|S| ℓ -|S| Tr A ℓ S = N ℓ=0 z ℓ ℓ i=0 (-1) ℓ-i N -i N -ℓ S:|S|=i Tr A ℓ S .
On the other hand, it follows from the definition of P(z),

Tr P(z) = ℓ≥1 ℓ π ℓ z ℓ ,
and identifying the coefficients of degree ℓ ≥ 3 in Tr P(z) recovers Eq. ( 1), up to a factor 1/21 . This provides a simple proof of the Perepechko-Voropaev formula, which otherwise only follows from Khomenko and Golovko's difficult approach. It shows furthermore that the formula is valid on directed graphs and for cycles of length ℓ ≥ 1.

Counting simple paths from weakly connected sets

A digraph is said to be weakly connected if replacing all its directed edges by undirected edges produces a connected undirected graph. The expression of P(z) can be reduced to a sum over weakly connected subgraphs of G owing to the simple property that the adjacency matrix of a digraph with more than one weakly connected component is block diagonal. In particular, letting G(S) = G(C 1 ) ∪ . . . ∪ G(C k ) be the decomposition of G(S) into weakly connected components, we have for all n ≥ 1,

W n S = W n C 1 + . . . + W n C k ,
Let C ⊆ S denote the non-empty subsets of V for which the resulting induced subgraph is weakly connected. For C ∈ C, the weak neighborhood N(C) of C in G is the set of vertices in V \ C that can reach and/or be reached from C in one step. Formally,

N(C) = {i ∈ V \ C : ∃j ∈ C, (i, j) ∈ E and/or (j, i) ∈ E}.
Of course, this definition recovers the classical definition of neighborhood in undirected graphs.

Theorem 4.1. The matrix generating series of open and closed simple paths verify:

i) P op (z) = C∈C (zW C ) |C|-1 (I -zW C ) |N (C)| , ii) P cl (z) = C∈C Diag (zW C ) |C| (I -zW C ) |N (C)| .
From a computational point of view, the restriction to weakly connected subgraphs provides a clear improvement for counting simple paths, especially for sparse graphs for which connected induced subgraphs are relatively rare.

A direct application concerns Hamiltonian paths, i.e. simple paths of maximal length, which appear in the terms of maximal degrees in the expressions of P op (z) and P cl (z). From the formulas in Theorem 4.1, the terms of maximal degrees correspond to connected subsets C such that C ∪ N(C) = V , i.e. dominating sets. This leads to an expression of the number of Hamiltonian paths that only involves the connected dominating sets of G.

For the next result, we denote by H the Hamiltonian path counting matrix, whose (i, j)-entry gives the number of Hamiltonian paths from i to j. Proposition 4.2. Let D be the set of weakly connected dominating sets in G,

H = D∈D (-1) N -|D| A N -1 D + Tr A N D N I .
The number of (oriented) Hamiltonian cycles appears in any diagonal entry of H,

π N = 1 N D∈D (-1) N -|D| Tr A N D ,
where N -|D| can be replaced by |N(D)| due to the dominating property. This expression improves on Bax's and Khomenko and Golovko's (Theorem 4 in [START_REF] Khomenko | Identifying certain types of parts of a graph and computing their number[END_REF]) results for Hamiltonian cycles in that the sum is restricted to connected dominating sets. Similarly, the total number of open Hamiltonian paths follows by summing all non-diagonal entries of H.

A quick inspection of the proof reveals that, similarly as for P(z), open and closed Hamiltonian paths can be dealt with separately yielding a slightly stronger version of the result, namely

H op = D∈D (-1) N -|D| A N -1 D and H cl = 1 N D∈D (-1) N -|D| Tr A N D I.
From a computational point of view, Proposition 4.2 shows that the problem of counting Hamiltonian paths in a digraph can be polynomially reduced to enumerating the weakly connected dominating sets. For undirected graphs, the later issue has received some attention in the literature [START_REF] Haynes | Fundamentals of domination in graphs[END_REF], although the algorithms developed in this purpose mainly focus on minimal connected dominating sets [START_REF] Mamadou | Enumeration of minimal dominating sets and variants[END_REF]. Nevertheless, connected dominating sets are easily tractable from the minimal ones by adding arbitrary vertices. We believe that this simplification may suffice to reduce the time complexity of some of the algorithms mentioned in [START_REF] Björklund | Determinant sums for undirected hamiltonicity[END_REF] for counting Hamiltonian paths, although further work is necessary at this point.

Proof of Theorem 4.1

We prove the result for P op (z), the argument is similar for P cl (z). We have, from Theorem The proof for P cl (z) is identical.

Proof of Proposition 4.2.

Take W = A in Theorem 4.1 and isolate the term of maximal degree. In P op (z), this term writes

3 . 2 P

 32 op (z) = S∈S (zW S ) |S|-1 (I -zW S ) N -|S| .Let C ⊆ V such that G(C) is weakly connected (i.e. C ∈ C). We regroup the terms of the sum involving W C which correspond to subsets S containing C as a weakly connected component. These subsets are of the form S = C ∪ T for T ⊆ V \ (C ∪ N(C)). We haveT ⊆V \(C∪N (C)) (zW C ) |C∪T |-1 (I -zW C ) N -|C∪T | = (zW C ) |C|-1 (I -zW C ) |N (C)| T ⊆V \(C∪N (C)) (zW C ) |T | (I -zW C ) N -|C|-|N (C)|-|T | . Let k = N -|C| -|N(C)|, remark that T ⊆V \(C∪N (C)) (zW C ) |T | (I -zW C ) k-|T | = k j=0 k j (zW C ) j (I -zW C ) k-j = I.Thus, P op (z) = S∈S (zW S ) |S|-1 (I -zW S ) N -|S| = C∈C (zW C ) |C|-1 (I -zW C ) |N (C)| .

D∈D (- 1 ) 1 D,

 11 |N (D)| A N -where we used that |N(D)| = N -|D| due to the dominating property. This provides the off-diagonal part of H. Its diagonal part follows from noting that simple cycles of length N appear in every term of maximal degree in P cl (z).

The absence of the factor 1/2 is simply due to the modified definition of simple cycles on digraphs which requires to take the orientation in consideration.

Acknowledgements

P.-L. Giscard is grateful for the financial support from the Royal Commission for the Exhibition of 1851.

Appendix

Proof of Theorem 3.2.

For a path p, V (p) denotes the set of vertices visited by p (its support). To prove the first part, write

We exploit the fact that a path p is both simple and open if, and only if,

For a path p from i to j to appear in the (i, j)-entry of W

|S|-1+k S

, one needs that V (p) ⊆ S and ℓ(p) = |S| -1 + k. Thus, the coefficient of p z ℓ(p) in the right-hand side of the above equation is given by

.

It is now only a matter of simplification to show that this coefficient equals


Mourad Benmalek
email: mbenmalek@esi.dz

Yacine Challal
email: ychallal@esi.dz

MK-AMI: efficient Multi-group Key management scheme for secure communications in AMI systems

Keywords: Smart Grid (SG), Advanced Metering Infrastructure (AMI), Cyber Security, Key Management

The Smart Grid (SG) is widely considered to be the informationization of the power grid. Advanced Metering Infrastructure (AMI) has been regarded as a key component of the SG. The critical role of AMI in the SG has made this system a privileged target of cyber attacks. Consequently, AMI security is of very high importance for the security of the SG. For this reason, Key Management has been identified as one of the most challenging topics in AMI development because of the great scale of SG and dynamism of connected clients with respect to tariff programs. This paper proposes a new efficient Multi-group Key management for AMI (MK-AMI) to secure data communications in the smart grid. It is a novel key management scheme that can support unicast, multicast and broadcast communications. An analysis of security and performance, and a comparison of our scheme with recently proposed schemes illustrate that MK-AMI achieves efficient key management and induces low storage and communications overheads compared to existing solutions.

I. INTRODUCTION

Smart Grid (SG), also called intelligent grid, intelligrid or futuregrid refers to the next generation power grid in which the electricity distribution and management is upgraded by incorporating advanced two-way flows of electricity and information and pervasive computing capabilities for improved control, agility, efficiency, reliability, economy and safety [START_REF] Sabbah | A survey of networking challenges and routing protocols in smart grids[END_REF]. A practical example of the benefits of introducing the smart grid includes the greater availability of electricity to homes at a lower cost, and the integration of distributed and renewable power generation such as local solar [START_REF] Fadlullah | Toward Intelligent Machine-to-Machine Communications in Smart Grid[END_REF].

Advanced Metering Infrastructure (AMI) has been regarded as a key component of the smart grid. It inherits the two-way communication capability of smart grid and introduces new opportunities for consumers and suppliers: it is responsible for collecting consumer's time-based data and transmitting them to the AMI host system, and it is also responsible for implementing price signals and control commands to perform necessary control actions [START_REF] Anzalchi | A survey on security assessment of metering infrastructure in Smart Grid systems[END_REF]. AMI is a privileged target for security attacks with potentially great damage against infrastructures and privacy. Consequently, security is one of the most challenging topics in AMI development Some of the key AMI security requirements are described below [START_REF] Anzalchi | A survey on security assessment of metering infrastructure in Smart Grid systems[END_REF]:

• Confidentiality-Requirement that customer's sensitive data is accessible only to authorized systems, customers do not want unauthorized people or marketing firms to know how much energy they are using, what their pattern of energy usage is, or other energy-related information. • Integrity-Requirement that the transmitted data from smart meters to the utility as well as control commands (such as Demand Response DR mechanisms that enable customers to cut down energy usage at peak times for example) are authentic, complete, and without unauthorized deletions, modifications, or additions. • Availability-Requirement that data (informations and control commands) is accessible by the smart meters and the utility whenever they need it. • Accountability (non-repudiation)-Requirement that the entities receiving the data will not subsequently deny receiving and vice versa. To meet these security requirements, cryptographic countermeasures must be deployed. However, cryptographic mechanisms for AMI require also an efficient key management. Inadequate key management can result in possible key disclosure to attackers, and even jeopardizing the entire goal of secure communications in AMI. Therefore, key management is a critical process to ensure the secure operation of AMI.

Several key management schemes (KMS) have been proposed [START_REF] Kamto | Light-weight key distribution and management for Advanced Metering Infrastructure[END_REF]- [START_REF] Nabeel | Scalable end-to-end security for advanced metering infrastructures[END_REF], but none of them can completely satisfy the security requirements mentioned previously. Hence, we propose a new key management scheme based on an efficient and scalable multi-group key graph technique to secure unicast, multicast, and broadcast communications in a SG network while achieving the security requirements of AMI.

The remainder of this paper is organized as follows. We discuss related works in Section II. In Section III we present our multi-group key management scheme. A security and performance analysis is performed in Section IV. Finally, we draw our conclusions in Section V.

II. RELATED WORKS

In recent years, several schemes have been proposed to secure communications for AMI in smart grid. According to [START_REF] Nabeel | Scalable end-to-end security for advanced metering infrastructures[END_REF], key management has been identified as a fundamental security challenge in an AMI.

A key distribution and management scheme for large customer networks to achieve authentication, privacy and data confidentiality in AMI is proposed by Kamto et al. in [START_REF] Kamto | Light-weight key distribution and management for Advanced Metering Infrastructure[END_REF].

Yan et al. [START_REF] Yan | A secure and reliable in-network collaborative communication scheme for advanced metering infrastructure in smart grid[END_REF] proposed an integrated approach in which trust services, integrity and data privacy could be provided by mutual authentications. In [START_REF] Li | Multicast authentication in the smart grid with one time signature[END_REF], Li and Cao proposed a onetime signature scheme to address the problem of preventing message forgery attacks in multicast communications. The proposed scheme presents a significant reduction in the storage and communication overhead, but only focuses on communication integrity and do not address confidentiality.

Nicanfar et al. [START_REF] Nicanfar | Smart grid authentication and key management for unicast and multicast communications[END_REF] developed a key management protocol for data communication between the utility server and customers' smart meters based on the concept of ID-Based public/private key pair model [START_REF] Shamir | Identity-based cryptosystems and signature schemes[END_REF]. Although the proposed key management protocol aims to reduce the computation overheads, the synchronization process still demands considerable computation efforts. Wu and Zhou [START_REF] Wu | Fault-tolerant and scalable key management for smart grid[END_REF] combines symmetric key technique based on the Needham-Schroeder authentication protocol [START_REF] Needham | Using encryption for authentication in large networks of computers[END_REF] and elliptic curve public key technique [START_REF] Miller | Use of Elliptic Curves in Cryptography[END_REF] to provide a novel key management scheme for smart grid assuring strong security, fault-tolerance, efficiency and scalability. In the work of Xia and Wang [START_REF] Xia | Secure Key Distribution for the Smart Grid[END_REF], the authors showed that Wu and Zhou scheme [START_REF] Wu | Fault-tolerant and scalable key management for smart grid[END_REF] is vulnerable to the man-in-the-middle attack and proposed an improvement for this scheme based on a trusted third party. However, these two schemes do not support secure multicast communications.

A key management scheme is proposed by Liu et al. [START_REF] Liu | A Key Management Scheme for Secure Communications of Advanced Metering Infrastructure in Smart Grid[END_REF] to secure unicast, multicast, and broadcast communications in AMI. Authors pretend that this scheme was based on the key graph management approach [START_REF] Wong | Secure group communication using key graphs[END_REF] but it suffers from a lack of scalability due to inefficient key management that results in non-negligible communication overhead for such a large-scale system. Moreover, we found that Liu's et al. scheme is not tolerant to packet loss. Wan et al. [START_REF] Wan | SKM: Scalable Key Management for Advanced Metering Infrastructure in Smart Grids[END_REF] proposed an improvement for Liu's et al. scheme (called SKM) that combines an adapted identity-based cryptosystem [START_REF] Chen | Identity based authenticated key agreement protocols from pairings[END_REF] and One-way Function Trees (OFT) approach [START_REF] Mcgrew | Key establishment in large dynamic groups: Using one-way function trees[END_REF] for multicast key management. The use of an OFT separately for each DR project (DR projects are programs designed to decrease electricity consumption or shift it from on-peak to off-peak periods depending on consumers' preferences) results in nonnegligible overhead for key storage.

Recently, Nabeel et et al. [START_REF] Nabeel | Scalable end-to-end security for advanced metering infrastructures[END_REF] proposed a PUF-based key management scheme for advanced metering infrastructures. Although, their scheme supports decentralized key management, the scheme requires a PUF hardware device.

III. MK-AMI: EFFICIENT MULTI-GROUP KEY

MANAGEMENT FOR AMI We introduce a new scalable and efficient key management scheme that we call efficient Multi-group Key management scheme for secure communications in AMI systems (MK-AMI). It is based on a novel multi-group key graph structure that supports the management of multiple Demand Response projects simultaneously for each customer.

A. Assumptions

1) The AMI complies with the architecture illustrated in Fig. 1 and involves Smart Meters (SMs) and Distributed Energy Resources at the user end, communication networks to connect two ends, Meter Data Management Systems (MDMS) and the means to integrate the collected data into software application systems at the utility end.

2) A specific default DR project is mandatory for all users of the SG. This default DR project will be used by MDMS to broadcast control messages or informations to all customers.

3) Except the mandatory DR project, any user can join or leave any DR project at any time.

B. Initialization of the KMS

Let us consider a set of n smart meters. Initially:

• A specific method of securely exchanging cryptographic keys over a public channel is used to establish individual keys {k 1 , ..., k n } between the MDMS and SMs. Theses individual keys (refreshed periodically) will be used to secure unicast communications, and to generate the multigroup key graph for secure multicast communications. • The MDMS must generate a group key GK 0 (refreshed periodically) for the default DR project. This key will be generated and transmitted through secure channels for each SM, and will be used to secure messages transmitted in broadcast mode. In Table I, we summarize the terminology that we will use throughout the remaining of this paper.

C. Group Key Management

In our solution, we propose a secure, efficient and scalable management of group keys in AMI. To address the scalability issue, key graph techniques can be used. Specifically, we adopt a variation of OFC (One-way Function Chain) [START_REF] Canetti | Multicast security: a taxonomy and some efficient constructions[END_REF]. We consider that the MDMS and all users individually compute all the keys of interior nodes using a pseudo-random function. However, the group keys are always chosen by the MDMS.

Moreover, as the users can subscribe to multiple DR projects at the same time, an intuitive solution is to use a key tree for each DR project. But, this naive application of OFC may be costly and induces a non-negligible key storage overhead. To reduce storage and communication costs in key management, we propose a novel multi-group key graph structure.

1) Multi-group Key Graph Structure: Our multi-group key graph structure can be modeled as shown in Fig. 2: in the lower level, each OFC tree represents a set of users with the same first DR project subscription, the leaf node of the tree is a user's individual key and tree's root is the DR project's group key. The graph in the upper level represents combinations of root keys for users subscribing to multiple DR projects at the same time. Our key structure has the following properties:

• A user only belongs to one OFC tree in the multigroup key graph corresponding to his Home DR project. He holds a copy of his leaf secret key and all keys corresponding to the nodes in the path from his leaf to the root in this tree, • A user has all group keys of the other DR projects to which he is subscribed, • if a user leaves his Home DR project and remains subscribed to one or more DR projects, he will shift to a new OFC tree. These features ensure that a user will not subscribe and pay for the same project multiple times.

2) Rekeying operations: In our solution, when a user subscribes or leaves a DR project, rekeying consists of 3 operations: joining/leaving an OFC tree, shifting among trees.

a) Leave procedure: The leave procedure deals with the case when a user unsubscribes from a DR project (u i leaves DR j ). Let φ j = {u l /u l subscribed to DR j }, Let X jk = {u l /u l ∈ φ j and Home DR

(u l ) = DR k }, Let ω k = {u l /u l ∈ X jk and DR k ∈ set(u i )}, Let δ k = {u l /u l /
∈ X jk and DR k ∈ set(u i )}.

• Case 1: We consider a user who subscribed to one or multiple DR projects and leaves his Home DR project: The MDMS updates and renews keys according to Algorithm 1.

Algorithm 1 : Update keys when user leaves Home DR Function leaveHomeDR (ui, DRj) ; 1: Apply standard OFC approach in DRj tree to update GKj (GK j represents the new group key); 2: If N sub(ui) = 1 : 3:

MDMS → X jk :

Enc(Enc(GK j , GK k ), GKj)

4: Else : 5: MDMS → ω k : Enc(Enc(GK j , k right(GK k ) ), GKj) Enc(Enc(GK j , k lef t(GK k ) ), GKj) 6: MDMS → δ k :
Enc(Enc(GK j , GK k ), GKj)

7:

Shift user ui to OFC tree corresponding to his second subscription DRx using standard OFC apporach (without updating key GKx that ui already has)

• Case 2: We consider a user who is subscribed to multiple DR projects and leaves one DR project which is not his Home DR project (u i leaves DR j ) The MDMS updates and renews keys according to Algorithm 2.

Let

DR x = Home DR(u i ), Let ψ k = {u l /u l ∈ ω k and DR k = DR x }, Let π i = {k l /k l = right(k c ) or k l = lef t(k c ), k c ∈ P ath(u i )}.
Algorithm 2: Update keys when user leaves DR project Function leaveDR (ui, DRj) ; 1: Update GKj (GK j is the new group key); 2: MDMS → Xjj :

Enc(GK j , k right(GK j ) ) Enc(GK j , k lef t(GK j ) ) 3: MDMS → Xjx : kα∈π i Enc(Enc(GK j , kα), GKj) 4: MDMS → ψ k : Enc(Enc(GK j , k right(GK k ) ), GKj) Enc(Enc(GK j , k lef t(GK k ) ), GKj) 5: MDMS → δ k : Enc(Enc(GK j , GK k ), GKj)
Example: Let us consider the key graph in Fig. 2. When u 1 (subscribed to DR 1 , DR 2 and DR 3 ) leaves DR 2 which

MDMS → {u 6 , u 7 } : Enc(Enc(GK 2 , k 6-7 ), GK 2 ) (3) (c) update GK 2 for users in ψ k : MDMS → ψ 3 : Enc(Enc(GK 2 , k 14-15 ), GK 2 ) (4) MDMS → ψ 3 : Enc(Enc(GK 2 , k 16-17 ), GK 2 ) (5) (d) update GK 2 for users in δ k : MDMS → δ 4 : Enc(Enc(GK 2 , GK 4 ), GK 2 ) (6)
b) Join procedure: Algorithm 3 deals with the case when a user subscribes to a new DR project (u i joins j ). Send the new group key GK j to ui ; 4:

Send a notification to all users in φj about the application of the one-way function; 5: Else : 6:

Send a notification to all users in X jk about the application of the one-way function; 7:

Apply standard OFC approach in DRj tree whithout updating GKj.

IV. PERFORMANCE EVALUATION A. Security Analysis

1) Backward and forward Sercery: When a new user joins a DR project, he cannot learn previous group keys because he does not have access to previous group key, and hence backward secrecy is preserved. Moreover, when a user leaves a DR project, all affected keys will be changed and redistributed securely which prevents the departing user from having access to the new keys and hence forward secrecy is preserved.

2) Collusion freedom: Any set of users unsubscribed from a set of DR projects cannot deduce the current used DR projects keys, because all affected keys when any user leaves a DR project will be updated and new keys are independent.

B. Performance Analysis 1) Storage Cost:

We approximate the storage cost with the number of keys stored in the MDMS/SMs (Table II).

2) Communication Cost: The number of keys to be updated varies according to the position of the joining/leaving member.

a) Leave procedure

• Case 1: When u i leaves his Home DR project DR j (user subscribed only to one DR project) :

comCost = (h j + N pr -1)|K| (7)
|K| : the size of the key in bit.

• Case 2: When u i leaves his Home DR project DR j (DR k is the new Home DR project):

comCost = (h j + 2.A + B + h k )|K| + c (8)
A = N sub(u i ) B = N pr -N sub(u i )
The " + c" term is to specify on which group key we must apply the one-way function c = log 2 (N pr ).

• Case 3: When u i leaves one DR project DR j which is not his Home DR project DR l :

comCost = (2 + h l + 2.A + B)|K| (9)
b) Join procedure

• Case 1: When u i joins his Home DR project DR j :

comCost = h j |K| + c (10)
• Case 2: When u i joins a new DR project DR j which is not his Home DR project DR l : (m j -1) + 1

comCost = |K| + c (11)
N sub(u i ) j=1 (log 2 m j +1)+1 MK-AMI 2 Npr j=1 (m j -1) + 1 log 2 (|Home DR(u i )|)+ N sub(u i ) + 1
* n is the number of SMs, Npr is the number of DR projects, mj is the number of j th DR project members, N sub(ui) is the number of DR projects to which subscribes user ui. 3) Simulation:

• Simulation Model: We consider a SG with 1 million users.

The utility provides 15 DR projects to users. We assume that users arrival is modeled as a Poisson process with parameter λ (users/months), and given that there are no statistical studies of DR projects membership behavior for the moment, we assume that membership duration in each DR projects follows an Exponential law with parameter µ.

A typical user session starts by a join event, which can be followed by one or more events. At the end of a membership in a DR project, a user leaves this DR project. We will consider a session of 24 months. Average arrival rate λ is of 10000 users/month, and average membership duration µ is 4 months. We will use a 128b long symmetric keys.

• Simulation Results

Storage Cost

For MDMS, the storage cost can be afforded using special key servers. In contrast, the storage capacity of SMs is limited to 4-12 KB [START_REF] Wang | Cyber security in the Smart Grid: Survey and challenges[END_REF]. Fig. 3 (a), (b) and (c) show the average storage cost induced at SMs with respect to the number of subscribed DR projects. We can see that in our scheme, a SM stores fewer keys than that in SKM+ (reduction reaches 93% while number of subscribers is about 150000 and a user can subscribe to 15 DR projects at the same time) and little more keys than that in Liu's et al. scheme. This can be explained as follows: in Liu's et al. scheme, a SM stores one key for each subscribed DR project. In SKM+, authors used an OFT for each DR project, the number of keys stored will increase significantly when a user subscribes to new DR projects. Whereas, in MK-AMI the number of subscribed DR projects does not affect significantly the storage cost.

Fig. 4 (a), (b) and (c) show the average storage cost in SMs with respect to DR projects' size. In Liu's et al. scheme, the projects's size does not affect significantly the storage cost, SMs store only the group keys. Whereas,in SKM+ and MK-AMI the DR projects' size affects the storage cost, as the number of users increases, the storage cost increases due to the rise of the height of the used key trees, but SMs store much fewer keys in MK-AMI with respect to SKM+.

Communication Cost

Fig. 5 (a) and (b) show a comparison of average communication cost per event (join/leave) with respect to the number of subscribed DR projects at the same time. The bandwidth overhead due for the Liu's et al. scheme is remarkably higher than that of our scheme because of the inefficient multicast key management. In MK-AMI, bandwith overhead reduction reaches 99% with respect to Liu's et al. scheme (Fig. 5 (b)). Note that although SKM+ has less communication overhead than MK-AMI, but the difference is not significant. V. CONCLUSION

In this paper, we proposed a new key management scheme for AMI in SG. MK-AMI is an efficient and scalable key management scheme supporting unicast, multicast, as well as broadcast communications. The proposed scheme uses a novel multi-group key graph technique that supports the management of multiple and dynamic Demand Response projects simultaneously for each customer. In addition, the proposed KMS guarantees both forward and backward secrecy. The detailed security analysis and performance evaluation show that MK-AMI is secure and efficient for AMI systems in smart grid.

Fig. 1 :

 1 Fig. 1: Basic Components of an AMI System

Fig. 2 :

 2 Fig. 2: Example of our multi-group key graph structure

Algorithm 3 :

 3 Update keys when user joins a DR project Function joinDR (ui, DRj) ; 1: GK j = H(GKj); 2: If N sub(ui) ≥ 1 : 3:

Fig. 3 : 15 Fig. 4 :

 3154 Fig. 3: Average storage cost in SMs with respect to number of subscribed DR projects

Fig. 6 (

 6 Fig.5(a) and (b) show a comparison of average communication cost per event (join/leave) with respect to the number of subscribed DR projects at the same time. The bandwidth overhead due for the Liu's et al. scheme is remarkably higher than that of our scheme because of the inefficient multicast key management. In MK-AMI, bandwith overhead reduction reaches 99% with respect to Liu's et al. scheme (Fig.5 (b)). Note that although SKM+ has less communication overhead than MK-AMI, but the difference is not significant.Fig.6(a) and (b) show a comparison of average communication cost per event for the three schemes with respect to the

Fig. 5 :Fig. 6 :

 56 Fig. 5: Average communication cost by event with respect to number of DR projects

TABLE I :

 I Notation Table Number of the i th DR project members h i Height of i th OFC tree h i = log 2 (m i ) All keys corresponding to the nodes in the path from u i 's individual key to Home DR(u i ) Group key right(k i ) Right children of node k i in the tree lef t(k i ) Left children of node k i in the tree

	Notation	Description
	H(.)	A One-way hash function
	n	Number of SMs
	m i
	Npr	Number of DR projects
	N sub(u i )	Number of DR projects to which subscribes user u i
	Home DR(u i ) First DR project to which subscribes user u i
	set(u i )	Set of DR projects to which subscribes user u i
	DR i	The ith DR project
	GK i	Group key of DR i
	P ath(u i )
	Enc(M, k)	Message M encrypted with key k
	A → B : M	A sends a message to B

TABLE II :

 II Storage Cost

	Scheme		MDMS	Storage Overhead	SM i
	Liu's et al., 2013 [10]		n + Npr + 1	N sub(u i ) + 2
	SKM+, 2014 [11]	2	Npr j=1


Liquid nitrogen to room temperature thermometry using niobium nitride thin films Olivier Bourgeois, * Emmanuel André, Cristina Macovei, and Jacques Chaussy Institut Néel, CNRS-UJF, 25 avenue des Martyrs, 38042 Grenoble, France (Dated: February 21, 2007)

Niobium nitride thin film thermometry has been developed for the temperature range of 70K to 300K. The deposition parameters have been optimized in order to get the best performances, i.e. the highest temperature coefficient of resistance (TCR), up to 300K. The TCR is found to be largely higher than 1% as the temperature is lowered from 300K, up to 6% at 77K. These significant performances are compared to the one of regular platinum thermometer as well as to other resistive thermometer: semiconductor type or amorphous metal to insulator transition materials. It is discussed how the properties of the NbN thin films could be due to a high temperature Mott transition.

Materials exhibiting metal to insulator transition have been already widely used as resistive thermometer at low temperature [1,2] and over larger temperature ranges [3]. At low temperature, the interesting properties of niobium nitride (NbN) as well as amorphous Nb-Si have been extensively used to build highly sensitive thin film resistive thermometers due to the significant increase of resistance as the temperature is lowered [2,3,4]. There is nowadays an increasing need for such highly sensitive thermometers for measurement closer to the room temperature, especially in calorimetry for applications in nanomagnetism [5], biophysics [6,7] or infra-red bolometry [8] where specific minute thermal or thermodynamic signatures are being explored.

We report here an original work exploring the performances for thermometry of highly nitrogen doped niobium thin films up to room temperature. The films studied here show characteristics around one order of magnitude better than regular metallic resistive thermometer like platinum thin film in terms of temperature coefficient. These performances are comparable in magnitude to what has been obtained in amorphous Nb-Si [3] or semiconducting like systems as amorphous YBaCuO thin films [8]. By increasing the nitrogen doping we are able to push the metal to insulator transition to much higher temperature than before. The insulating nature of these nitrogen doped niobium materials arises already at room temperature with reasonable thermometer impedance allowing the electrical measurement with a regular low noise measurement chain. Moreover, the deposition process allows the integration of such thin films into any micro or nanofabricated devices.

The NbN thin films were all prepared following the same technique. The films are deposited by reactive sputtering using a pulsed power supply magnetron on a 3.5N purity Nb target in a mixture of nitrogen and argon gases. The base pressure is about 10 -7 mbar and during the deposition the pressure in the chamber is regulated at 2×10 -2 mbar. The composition of the gas mixture is controlled by two mass flow debimeter. The proportion * olivier.bourgeois@grenoble.cnrs.fr of the two different gases is an adjustable parameter used to control the level of nitrogen doping in the film. The exact proportion of argon and nitrogen in the plasma is also controlled via a mass spectrometer connected to the sputtering chamber. The mass spectrometer also works at low pressure to check the possible existence of impurities like organic materials, oxygen and the presence of water; their presence has been evaluated to be less than 10 ppm. The sample holder is maintained at room temperature using circulating water avoiding any overheating of the thin film during the deposition procedure. The plasma is created by a pulsed power supply with a working frequency tuned between 0 and 350 kHz, the higher the frequency the higher the nitrogen doping in the film. The typical thickness of the films measured with a profiler is about 100nm. The thin films are deposited on sapphire substrate on the voltage and current electrodes as it is depicted in the inset of the Fig. 1.

After the deposition, the films are annealed under vacuum at 150 • C for 2 hours. This is a crucial step in the process to stabilize the electrical properties of the thermometer. It prevents the transducer from spurious aging processes and improve its intrinsic temperature coefficient; this has been also reported by other authors [4,9]. Four different samples are presented in this study, three (sample # 1, 2 and 3) adapted to high temperature measurement and one, sample #4, used as a reference sample. The reference sample is also a NbN thin film having a lower nitrogen doping; this kind of thermometer are regularly used for low temperature calorimetric measurements [10]. Rutherford backscattering (RBS) analysis giving the atomic ratio of nitrogen and niobium is used as a posteriori control of the composition of the films. All the deposition parameters are summarized in the Table I; it is worth noticing that as it will be shown in the following the eletrical properties of the films are very sensitive to the atomic proprotion of nitrogen and niobium: x=N/Nb.

The electrical characterizations are performed using the four probe measurement technique: two leads used to applied the current and two leads for the measurement of the voltage. The thin film thermometer is rectangular but the surface between the two voltage leads where the current density is uniform is a square of 1mm by 1mm. The measuring current varies from 1nA to 1µA; IV characteristic were controlled to be linear in that current range. The typical resistance of the films at room temperature lies between 800 Ohm to few kiloOhm. A specific low noise measurement chain has been used based on lock-in amplifier technique and very low noise preamplifier. The resistances have been measured versus the temperature between 70 to 300 K in a 4He cryostat. The variation of the resistances per square of the four samples is presented in the Fig. 1. We choose to plot the resistance per square (R ) because it is then easier to adapt the thermometer for different applications depending on the chosen geometry of the device transducer. The resistance per square is the resistance of a thin film having a square surface in the plane of the circulating current.

A good criterion for a first characterization of the NbN film is the measurement of the resistive ratio (RR), the ratio between the resistance at 77 K and the resistance at 300K. For a low temperature thermometer that ratio has to be around 3 [1], but in our case in order to get a good thermometer in the temperature range of 70K to 300K the RR has been pushed up to 100. All the electrical characterization data are given in the Table II.

It is clear that all the samples exhibit a significant increase of resistance as the temperature is lowered, except for the reference sample. For some sample the increase represents several order of magnitude; it is correlated to a metal to insulator transition occuring early in temperature as it will be discussed further at the end of the paper.

For application in oscillating temperature calorimetry (or ac calorimetry) the figure of merit we favour is I for the details of the deposition parameters) versus the temperature. All samples exhibit an increasing resistance as the temperature is lowered due to a metal to insulator transition. The inset depicts the geometry of the contacts, four electrodes (I+, I-, V+, V-) crossing the thin film. the temperature coefficient of resistance (TCR) given by: α = -1/R × dR/dT ; the minus sign is only here to give a positive value for α. Indeed, the signal to noise ration in a ac calorimetry experiment is directly proportional to α. The TCR of each samples have been calculated by numerical derivation and are given in the Fig. 2. We can first notice that at 70K the TCR is up to 6% which is very high when it is compared to a regular thin platinum film in which the TCR is around 1%, idem at 300K the NbN TCR is around 1% and the TCR of a thin film of platinum is not above 0.25%. Hence the performance of such thin amorphous film can be of great interest for thermometry between 70 and 300K and allows many applications for highly sensitive room temperature calorimetry. It is worth to stress that in the three samples # 1, 2 and 3, the TCR is about the same, which shows that whatever the resistance at room temperature the figure of merit is relatively constant; this is something important in terms of reproducibility. The TCR of the reference sample (#4) has been given as a comparison, it is indeed much smaller at room temperature than that of the other samples by one order of magnitude.

If the nitrogen doping is still increased a saturation in the electrical properties is observed resulting in an increase of resistance but in no change in the temperature coefficient. This behavior is of no interest because the Johnson noise of the transducer increases with no gain in the temperature coefficient leading to a decrease of the global performances of the thermometry. The properties of these films seem comparable to what was reported about Nb-Si in reference [9]; it has been used also for calorimetric measurements up to room temperature [12]. Finaly, the NbN thermometer seems at 300K a little bit less performant than amorphous YBaCuO as reported in ref [8].

In order to better describe the insulating character of the NbN thin films the resistance has been plotted on a logarithmic scale versus the inverse of temperature to the power 1/4; in the Mott law of disordered insulator the resistance is indeed given by R(T ) = R 0 exp(T 0 /T ) 1/4 [13] in the model of variable range hoping in three dimension samples. If the films follow a Mott law the variation of the resistance should be linear in the plot as it is illustrated in the Fig. 3. It is interesting to point out that even at high temperature this materials keeps its Mott insulator features, even if there is no current theoretical extension of this model at room temperature. Nevertheless, we anticipate no major problem prolonging the concept of Mott transition up to 300K. This confirms the origin on the increase of resistance as the temperature is lowered.

In conclusion, we are able, by increasing the number of nitrogen dopant in niobium thin film, to strongly improve the performances of NbN thermometer between liquid nitrogen and room temperature. This has been done by dragging the Mott transition up to much higher temperature. This type of thermometer will be easily implemented as transducer in micro and nanodevice in calorimetric and more generally in thermal and thermodynamic measurement for highly sensitive experiment.

FIG. 1 :

 1 FIG.1: Resistance per square of the four different samples (measured every mK) (see table I for the details of the deposition parameters) versus the temperature. All samples exhibit an increasing resistance as the temperature is lowered due to a metal to insulator transition. The inset depicts the geometry of the contacts, four electrodes (I+, I-, V+, V-) crossing the thin film.

FIG. 2 :

 2 FIG.2: TCR of the four samples versus temperature obtained by numerical derivation of the measured resistance. The value of the TCR increases from about 1% at 300K to 6% at 70K.

FIG. 3 :

 3 FIG.3:The resistance of the NbN thin films are plotted on a logarythmic scale against 1/T 1/4 in order to show the Mott behavior of these amorphous disordered metal-insulator materials. The reference sample #4 has not been plotted for the sake of simplicity. The behavior of ln(R) for the reference sample is almost constant in temperature at 10.2.

TABLE I :

 I Sputtering deposition parameters of the different NbN samples. τAr and τN 2 represent the gas (Ar, N) flow in the plasma in sccm (standard cubic centimeter per minute), f the excitation frequency of the plasma in kHz, x=N/Nb the composition of the film measured by RBS (see text) and t the deposition time in minute. The sample # 4 is considered as a reference sample.

	samples	τAr	τN 2	f	x	t
	#1	1	30	350	1.65	10
	#2	1	30	350	1.65	10
	#3	4	30	350	1.6	10
	#4 (ref.)	4	14	150	1.4	5

TABLE II :

 II Summary of the resistive behavior of the NbN thin films. The resistance given in the table is the resistance per square in Ohm.The sample # 4 is considered as a reference sample.

	sample R300K R77K ρ (mΩ.cm) RR α300K (%) α77K (%)
	#1	2678 186600	26.78	70	0.9	5.3
	#2	1745 61100	17.45	35	0.7	3.5
	#3	643 62200	6.43	96.4	0.83	5.3
	#4 (ref.) 469	750	4.69	1.6	0.1	0.27

We would like to thank M. Aprili, J-L. Bret, P. Brosse-Maron, T. Fournier, J-L. Garden, P. Lachkar, M. Tourniol for usefull help and F. Lalu for Rutherford backscattering measurements.


D Fraix-Burnet

S Girard

Statistics for Astrophysics: Clustering and Classification

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 TO R

Didier Fraix-Burnet 1,2 Abstract. I present a small introduction to R useful for the School. You need to execute and play with the commands.

Introduction

R is a freely available language and environment for statistical computing and graphics which provides a wide variety of statistical and graphical techniques: linear and nonlinear modelling, statistical tests, time series analysis, classification, clustering, etc. It comprises libraries of nearly all the statistical literature. Please consult the R project homepage 1 for further information.

R is an object-oriented interpretated language, made of an ensemble of packages. Packages are installed through a compilation of codes C++ and FORTRAN (gcc compiler). R is command-line driven, making it a wonderful tool for complex scripts and repeated procedures.

R is used by most statisticians in the world mainly for research and development. There are more and more packages for astronomical purposes, with a package to manage FITS format: FITSio.

CRAN 2 (Comprehensive R Archive Network) is a network of ftp and web servers around the world that store identical, up-to-date, versions of code and documentation for R. Please use the CRAN mirror nearest to you to minimize network load.

There are many and many help websites and blogs 3 . Simply asking your question in your favorite internet search engine, you will surely find the corresponding help and even full scripts. To make a search on the web regarding R, type "r-cran" instead of "R"! For MATLAB users, there are tables of command equivalence4 .

You can find a Short Reference Card5 on the School website with the Codes and Data6 .

Installation of R

You can download a binary version for your platform. For Linux users, you may check whether R is in the distribution.

Warning: once R is installed, if your are behind a proxy, it should be indicated with a line like: http_proxy='YOUR PROXY ADDRESS:3128/' in a file ∼/.Renviron or at the end of the file /etc/R/Renviron to be able to install and update packages.

Graphical User Interface (GUI) or Integrated Development Environment

A non-necessary but appreciable complement is to have a GUI interface to R which provides an integrated environment with the command line terminal, the plot window and an editor for the script files among others. This is extremely useful to save commands into a script file since you can execute any selection in this file.

I would strongly recommend RStudio7 . Otherwise there are many others, one being included in the Mac binary.

Packages

Packages are ensembles of programs, commands and scripts. They must be installed, i.e. compiled, once. Then they can be loaded to be used or detached when you are done.

Warning: commands of the newly installed packages hide commands with same name already in memory until the new package is uninstalled ("detached"). These commands are indicated when the package is loaded. Some commands automatically adapt their behaviour to the type of objects (for instance, plot can have several slightly different behaviours). For this reason, and except for some basic set of packages, loaded packages are not kept into memory after session's end.

Packages can be installed directly from a repository (you will be proposed a choice of mirrors) or from a previously loaded file.

To install packages:

i n s t a l l . packages ( " f p c " )

To load packages into memory and access its commands:

l i b r a r y ( f p c )
To remove a package from memory:

detach ( f p c )
Sometimes it is useful to update packages:

update . packages ( r e p o s=" h t t p : // c r a n . univ-l y o n 1 . f r /" , ask=F , c h e c k B u i l t=T)

1.4 Files R uses two important files which are saved under the current directory:

• .Rhistory: the history of the commands

• .RData: the environment file made of the functions and objects in memory. The loaded packages are not kept when you quit R to avoid command confusions in the next session.

It is a good idea to have different directories for different projects. You can move to the right directory before loading R, or change directory within R:

setwd ( " t h e d i r e c t o r y path " )

RStudio encourages you to create projects so that it is very easy to handle several works and switch from one to the other.

Basics

To open R, type R.

When in R , the prompt > will appear. It will not be indicated in the following command lines (to ease copy and paste operations...).

To quit:

q ( )
To obtain help:

help ( plot ) ? plot ?? plot # t

o f i n d a command i n u n l o a d e d p a c k a g e s

Help pages all have the same structure with description, usage, format (options), Details, Sources and Examples at least. Do not hesitate to run the examples.

Commands

Commands are always followed by closed parentheses. Otherwise it prints the content of the object.

l s ( ) # p r i n t s t h e r e s u l t o f t h e command l s # p r i n t s t h e command whi ch i n t h i s c a s e i s a f u n c t i o n

If parentheses are not closed, then the symbol + appears asking for a missing closing ).

Two commands written on the same line have to be separated by ";". The function rm() removes (definitvely) an object from memory:

a <-" H e l l o " ; a rm( a ) ; a

Objects

Objects have values and can be of several data types. To give the value 1 to object a, simply enter equivalently:

a <-1 a = 1
To print values on the console: a a * 2 print ( a ) # a i s an o b j e c t print ( " a " ) # ' ' a ' ' i s a c h a r a c t e r exp ( p i )

Operators

a = 3 ; b = 6 a <= b a != b ( b-3==a ) & ( b>=a ) ( b ==a ) | ( b>=a )

Functions

Functions are defined as follows:

my function <-function ( par=2) {par ˆ2} # d e f i

n e s t h e f u n c t i o n

# " m y f u n c t i o n " r e q u i r i n g one parameter " par " # whi ch i s g i v e n t h e d e f a u l t v a l u e "2" my function

# p r i n t s t h e f u n c t i o n m y f u n c t i o n my function ( ) # e x e c u t e m y f u n c t i o n w i t h t h e d e f a u l t v a l u e s # o f t h e p a r a m e t e r s my function ( 3 ) # e x e c u t e m y f u n c t i o n f o r t h e v a l u e "3"

It is a (very) good idea to write a complex function into a file myfunction.R. It can be loaded by sourcing the file: source ( ' ' my function .R' ' )

3 Data types

Variables

Vectors: the simplest data structure. Lists: lists are a general form of vector in which the various elements need not be of the same type, and are often themselves vectors or lists. Lists provide a convenient way to return the results of a statistical computation.

A = matrix ( 1 : 1 5 , ncol =5) x=l i s t (mat=A, text=" t e s t l i s t " , v ec=y ) x [ [ 2 ] ] ; x$ v ec
Matrices: matrices or more generally arrays are multi-dimensional generalizations of vectors. In fact, they are vectors that can be indexed by two or more indices and will be printed in special ways.

A = matrix ( 1 : 1 5 , ncol =5); A B = matrix ( 1 : 1 5 , nc =5 , byrow=TRUE) ; B A[ 1 , 3 ] ; A[ , 2 ] ; A[ 2 , ] ; A [ 1 : 3 , 1 : 3 ] i s . matrix (A) as . l i s t (A)
Arrays: an array can be considered as a multiply subscripted collection of data entries.

array ( 1 : 1 2 ,c ( 2 ,3 ,2 ) ) array ( c ( 1 : 8 ,rep ( 1 ,8 ) , seq ( 0 , 1 , l e n =8)) , dim = c ( 2 , 4 , 3 ) )

Data frames: data frames are matrix-like structures, in which the columns can be of different types. Think of data frames as data matrices with one row per observational unit but with (possibly) both numerical and categorical variables. Many experiments are best described by data frames: the treatments are categorical but the response is numeric. t = c ( 1 4 7 , 1 3 2 , 1 5 6 , 1 6 7 , 1 5 6 , 1 4 0 ) p = c ( 5 0 , 4 6 , 4 7 , 6 2 , 5 8 , 4 5 ) s = c ( "M" , "F" , "F" , "M" , "M" , "F" ) H = data . frame

Indices

It is important to understand the indices since it can become somewhat tricky for lists or complex data frames. The following examples should help you in this respect.

d = c ( 2 , 3 , 5 , 8 , 4 , 6 ) ; d d [ 2 ] d [ 2 : 3 ] d [ -3] d [ -( 1 : 2 ) ] ; d[-c ( 1 , 2 ) ] ( 1 : length ( d ) ) [ d>5] which ( d>5) c1 <-c ( 1 , 2 , 3 ) c2 <-c ( " a " , "b" ) c <-l i s t ( c1 , c2 , 1 0 0 ) c c [ 2 ] c [ [ 2 ] ] c [ [ 2 ] ] [ 2 ] c [ [ 3 ] ] c [ 3 ] c [ [ 3 ] ] + 1 c [ 3 ] + 1
The entry "NA" (Not Available) indicates a missing value:

d [ 3 ] =NA d summary( d ) i s . na ( d ) any( i s . na ( d ) ) a l l ( i s . na ( d ) ) help (NA)
Some functions can handle these missing values in different ways (ignore or replace).

Objects classes

Objects can be of different classes: integer, double, numeric, character, factor, logical. Factors provide compact ways to handle categorical data. Note that in RStudio, if you get the error "Error in plot.new() : figure margins too large", try to increase the size of the plot window. Here are some useful commands especially for command line sessions:

dev . new( ) # opens a new g r a p h i c a l window x11 ( ) # a l s o works dev . l i s t ( ) # l i s t a l l a v a i l a b l e g r a p h i c a l windows dev . se t ( n ) # t o a c t i v a t e t h e windows "n" dev . cur ( ) # g i v e s t h e number o f t h e c u r r e n t a c t i v e window dev . o f f ( n ) # c l o s e s window "n" plot . new( ) # e r a s e c u r r e n t window ( a l s o frame ( ) )

An illustration with simple classification tools

As an illustration of both the graphical capabilities of R and some basic classification techniques using kmeans and DBSCAN, we use the built-in database iris.

? kmeans kmeans colnames ( i r i s ) i r i s 2 <-i r i s [ , -5 ] colnames ( i r i s 2 ) checkkmeans <-function ( n s t a r t =1000){ i r i s . kmeans <-kmeans ( i r i s 2 , 3 , n s t a r t=n s t a r t ) plot ( i r i s 2 [ c ( " S e p a l . Length " , " S e p a l . Width " ) ] , c o l= i r i s . kmeans$ c l u s t e r ) points ( i r i s . kmeans$ c e n t e r s [ , c ( " S e p a l . Length " , " S e p a l . Width" ) ]

, c o l =1:3 , pch=" * " , cex =5) }

checkkmeans ( 1 )

# r e p e a t t h e command many t i m e s and compare t h e r e s u l t s checkkmeans ( 1 0 ) # idem . I s i t more s t a b l e ? i r i s . kmeans <-kmeans ( i r i s 2 , 3 , n s t a r t =1000) Now let's perform a classification using DBSCAN: l i b r a r y ( f p c ) i r i s . ds <-dbscan ( i r i s 2 , e p s =0.42 , MinPts =5) table ( i r i s $ S p e c i e s , i r i s . ds $ c l u s t e r ) p l o t c l u s t e r ( i r i s 2 , i r i s . ds $ c l u s t e r ) detach ( f p c ) Since there is generally not one best method, it is always a good idea to compare the results obtained using several approaches. Here is the comparison between kmeans and DBSCAN of the above analyses through a convenient graphical display:

table ( i r i s $ S p e c i e s , i r i s . kmeans$ c l u s t e r ) table ( i r i s . kmeans$ c l u s t e r , i r i s . ds $ c l u s t e r ) op <-par ( mfrow=c ( 2 , 1 ) ) plot ( i r i s 2 [ c ( " S e p a l . Length " , " S e p a l . Width " ) ] , c o l= i r i s . kmeans$ c l u s t e r ) points ( i r i s . kmeans$ c e n t e r s [ , c ( " S e p a l . Length " , " S e p a l . Width" ) ]

, c o l =1:3 , pch=" * " , cex =5) plot ( i r i s 2 [ c ( " S e p a l . Length " , " S e p a l . Width " ) ] ,

c o l= i r i s . ds $ c l u s t e r ) par ( op )

Other graphical examples

Here are other examples of graphics: x = runi f ( 5 0 , 0 , 2 ) y = runi f ( 5 0 , 0 , 2 ) plot ( x , y , main = " T i t r e " , x l a b = " a b s c i s s a " , y l a b = " o r d i n a t e " , c o l = " d a r k r e d " ) abl i ne ( h = . 6 , v = . 6 ) text ( . 6 , . 6 , "comment on t h e graph " ) c o l o r s ( ) 5 Import/export read . table ( f i l e . choose ( ) , s t r i n g s A s F a c t o r s = FALSE, h e a d e r=T, quote=" \ " " , na . s t r i n g s=" ? " , s k i p =0 ,comment . c h a r="#") tab <-read . table ( f i l e = ' ' m y f i l e . t x t ' ' ) write . table ( tab , f i l e=" m y f i l e 2 . t x t " , quote=FALSE, s e p=" " , na=" ? " ,row . names=F)

h i s t . norm = function ( n , c o l ) { x = rnorm ( n ) h = h i s t ( x , plot = F) s = sd ( x ) m = mean( x ) y lim = c ( 0 , 1 . 2 * max(max( h$density ) , 1/ ( s * sqrt ( 2 * p i ) ) ) ) x l
Beware:

save ( x , f i l e=" whatev er . RData" ) load ( " whatev er . RData" , e n v i r=n e w e n v i r )

The first command saves the object x in a (R) binary format file named whatever.RData. There is a compression option.

The second command loads the object included in the (R) binary format named whatever.RData. Be careful because this replaces all existing objects with the same names (here object x ) in the current environment. To be safe, if you are not sure, you should load the files in a new environment as indicated. In this case two objects with the same name can exist, each in a different environment.

To define the new environment and use it, this is rather simple:

n e w e n v i r <-new . env ( ) l s ( n e w e n v i r ) n e w e n v i r$x attach ( n e w e n v i r ) x detach ( n e w e n v i r ) rm( x , e n v i r=n e w e n v i r )

R

  c l a s s ( a ) c l a s s ( s [ 3 ] ) c l a s s ( c (TRUE,TRUE, FALSE) ) s t r ( s ) summary( s ) s t r ( as . factor ( s ) ) l e v e l s ( s ) l e v e l s ( as . factor ( s ) ) summary( as . factor ( s ) ) is very powerful for graphics 8 .

  x1 <-1 ; x2 <-3 ; x3 <-5 ; y1 <-2 ; y2 <-4 ; y3 <-3 plot ( x1 , y1 , x lim=range ( x1 , x2 , x3 ) , y lim=range ( y1 , y2 , y3 ) ) points ( x2 , y2 , c o l= ' b l u e ' ) points ( x3 , y3 , c o l= ' r e d ' ) l i n e s ( c ( x1 , y1 ) , c ( x2 , y2 ) ) abl i ne ( 1 , c ( 1 , 2 ) )

  a b = " Histogram and normal a p p r o x i m a t i o n " y l a b = " " main = paste ( " Gaussian sample : n=" , n ) h i s t ( x , f r e q = F , y lim = ylim , x l a b = x lab , y l a b = y lab ,c o l = col , main = main ) curve (dnorm( x , m, s ) , add = T, lwd = 2 ) } op = par ( m f c o l = c ( 1 , 3 ) ) h i s t .norm ( 2 0 0 , c o l = " y e l l o w " ) h i s t . norm ( 8 0 0 , c o l = " d a r k g o l d e n r o d " ) h i s t . norm ( 3 2 0 0 , c o l = " b l u e " ) par ( op )x = rnorm ( 1 0 0 0 ) h = h i s t ( x , f r e q = F) l i n e s ( density ( x ) )x = rnorm ( 1 0 0 ) op <-par ( mfrow=c( 2 , 2 ) , mar=c ( 3 , 3 , 1 , 1 ) , mgp=c ( 1 . 5 , 0 . 5 , 0) , oma=c ( 3 , 0 , 2 , 0 ) ) # t o p l a y w i t h margi ns b e t w e e n p l o t s boxplot ( x ) boxplot ( x , h o r i z o n t a l = T) boxplot ( x , c o l = " r e d " ) boxplot ( x , c o l = " o r a n g e " , b o r d e r = " d a r k b l u e " ) par ( op ) demo( graphics )

http://mathesaurus.sourceforge.net/octave-r.html

https://cran.r-project.org/doc/contrib/Short-refcard.pdf

The R codes for this chapter and the required data files can be found at http://stat4astro2015.sciencesconf.org/

http://www.rstudio.org/

https://en.wikibooks.org/wiki/R_Programming/Graphics https://www.stat.auckland.ac.nz/ ~paul/RG2e/ http://zoonek2.free.fr/UNIX/48_R/03.html http://zoonek2.free.fr/UNIX/48_R/04.html