Stéphane Zampelli
email: stephane.zampelli@emn.fr

Yves Deville
email: yves.deville@uclouvain.be

Christine Solnon
email: christine.solnon@liris.cnrs.fr

Solving subgraph isomorphism problems with constraint programming

Keywords: Subgraph isomorphism, Global constraint, Filtering algorithm

come

Introduction

Graphs are widely used in real-life applications to represent structured objects, e.g., molecules, images, or biological networks. In many of these applications, one looks for a copy of a pattern graph into a target graph [START_REF] Conte | Thirty years of graph matching in pattern recognition[END_REF]. This problem, known as subgraph isomorphism, is NP-complete [START_REF] Garey | Computers and Intractability[END_REF] in the general case.

There exist dedicated algorithms for solving subgraph isomorphism problems, such as [START_REF] Ullmann | An algorithm for subgraph isomorphism[END_REF][START_REF] Cordella | An improved algorithm for matching large graphs[END_REF]. However, such dedicated algorithms can hardly be used to solve more general problems, with additional constraints, or approximate subgraph isomorphism problems, such as the one introduced in [START_REF] Zampelli | Approximate constrained subgraph matching[END_REF].

An attractive alternative to these dedicated algorithms is Constraint Programming (CP), which provides a generic framework for solving constraint satisfaction problems (CSP). Indeed, subgraph isomorphism problems may be formulated as CSPs in a straightforward way [START_REF] Régin | Développement d'outils algorithmiques pour l'intelligence artificielle. application à la chimie organique[END_REF][START_REF] Rudolf | Utilizing constraint satisfaction techniques for efficient graph pattern matching[END_REF]. To make CP competitive with dedicated approaches for these problems, [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF] has introduced a filtering algorithm. [START_REF] Zampelli | Approximate constrained subgraph matching[END_REF] has extended this work to approximate subgraph isomorphism, and has shown that CP is competitive with dedicated approaches.

Contribution and outline of the paper

In this paper, we introduce a new global constraint and an associated filtering algorithm for the subgraph isomorphism problem. The filtering exploits the global structure of graphs to achieve a stronger partial consistency than when propagating edge constraints separately. This work takes inspiration from the partition refinement procedure used in [START_REF] Mckay | Practical graph isomorphism[END_REF][START_REF] Fowler | Efficient graph automorphism by vertex partitioning[END_REF][START_REF] Darga | Exploiting structure in symmetry detection for cnf[END_REF] for finding graph automorphisms: the idea is to label every node by some invariant property, such as node degrees, and to iteratively extend labels by considering labels of adjacent nodes. Similar labelings are used in [START_REF] Sorlin | A global constraint for graph isomorphism problems[END_REF][START_REF] Sorlin | A parametric filtering algorithm for the graph isomorphism problem[END_REF] to define filtering algorithms for the graph isomorphism problem: the idea is to remove from the domain of a variable associated with a node v every node the label of which is different from the label of v. The extension of such a label-based filtering to subgraph isomorphism problems mainly requires to define a partial order on labels in order to express compatibility of labels for subgraph isomorphism: this partial order is used to remove from the domain of a variable associated with a node v every node the label of which is not compatible with the label of v.

This paper extends preliminary results presented in [START_REF] Zampelli | Filtering for subgraph isomorphism[END_REF][START_REF] Zampelli | A constraint programming approcah to subgraph isomorphism[END_REF]. In Section 2, we introduce the subgraph isomorphism problem and we give an overview of existing approaches for solving this problem. The basic idea of a labelingbased filtering is described in Section 3: we first introduce the concept of labeling, and show how labelings can be used for filtering; then we show that labelings may be iteratively strengthened by adding information from labels of neighbors. The global constraint and a first filtering algorithm are described in Section 4. We introduce in Section 5 another filtering algorithm which ensures a lower consistency at a lower cost. Experimental results are described in Section 6.

Subgraph Isomorphism

Definitions and notations

A graph G = (N, E) consists of a node set N and an edge set E ⊆ N × N , where an edge (u, v) is a couple of nodes. In this paper, we implicitly consider undirected graphs, such that (u, v) ∈ E ⇒ (v, u) ∈ E. The extension of our work to directed graphs, which is straightforward, is discussed in 3.4.

A subgraph isomorphism problem (SIP) between a pattern graph G p = (N p , E p ) and a target graph G t = (N t , E t ) consists in deciding whether G p is isomorphic to some subgraph of G t . More precisely, one should find an injective function f : N p → N t that associates a different node of the target graph with every node of the pattern graph and that matches edges of the pattern graph, i.e., ∀(u, v) ∈ N p × N p , (u, v) ∈ E p ⇒ (f (u), f (v)) ∈ E t . Note that the subgraph is not necessarily induced so that two pattern nodes that are not linked by an edge may be matched to target nodes which are linked by an edge. The problem is also called subgraph monomorphism problem or subgraph matching in the literature. The function f is called a subgraph isomorphism function.

In the following, we assume G p = (N p , E p ) and G t = (N t , E t ) to be the underlying instance of subgraph isomorphism problem and we assume that N p ∩N t = ∅. We denote #S the cardinality of a set S. We also define N ode = N p ∪ N t , Edge = E p ∪ E t , n p = #N p , n t = #N t , n = #N ode, e = #Edge, d p and d t the maximal degrees of the graphs G p and G t , and d = max(d p , d t ). The set of neighbors of a node i is denoted adj(i) and is defined by adj(i) = {j | (i, j) ∈ Edge}.

Dedicated approaches

The subgraph isomorphism problem is NP-complete. Ullmann's algorithm [START_REF] Ullmann | An algorithm for subgraph isomorphism[END_REF] is one of the first algorithm dedicated to the graph and subgraph isomorphism problem itself. It includes a preprocessing step based on node degrees and basically performs a depth-first search algorithm. The search space is pruned by maintaining arc consistency on edge constraints. Ullmann's algorithm has been the baseline for more advanced algorithms and it is still cited in recent applications (see for example [START_REF] Grochow | Network motif discovery using subgraph enumeration and symmetry-breaking[END_REF]).

Another well known algorithm is vflib [START_REF] Cordella | Performance evaluation of the vf graph matching algorithm[END_REF], which is considered as the state-of-theart for subgraph isomorphism. The key points of vflib are the incremental building of connected graphs with cheap pruning rules and a clever choice of data structures leading to a lower memory footprint. vflib builds incrementally the assignment through a backtracking procedure, stopping whenever its extension is impossible. Assignments are extended by growing two connected graphs inside the pattern and the target graph respectively, and by checking conditions in the neighborhood of the current subgraphs.

[26] also solves subgraph isomorphism with this approach. The popularity of vflib is due to several papers [START_REF] Cordella | Performance evaluation of the vf graph matching algorithm[END_REF][START_REF] Cordella | Graph matching: A fast algorithm and its evaluation[END_REF][START_REF] Cordella | An improved algorithm for matching large graphs[END_REF] demonstrating the superiority of vflib against Ullman's algorithm, the public availability of its source code, and the systematic experiments over a synthetic database for graph matching problems [START_REF] Foggia | A database of graphs for isomorphism and sub-graph isomorphism benchmarking[END_REF].

CP Models for Subgraph Isomorphism

A subgraph isomorphism problem may be formulated as a CSP in a straightforward way [START_REF] Régin | Développement d'outils algorithmiques pour l'intelligence artificielle. application à la chimie organique[END_REF][START_REF] Rudolf | Utilizing constraint satisfaction techniques for efficient graph pattern matching[END_REF][START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF]. A variable x u is associated with every node u of the pattern graph, and its domain D(x u ) is the set of target nodes. A global AllDiff constraint [START_REF] Regin | A filtering algorithm for constraints of difference in CSPs[END_REF] ensures that the matching function is injective. Edge matching is ensured by a set of #E p binary constraints denoted by c2 and defined as follows:

∀ (u, v) ∈ E p , c2(x u , x v ) ≡ ((x u , x v ) ∈ E t ) .
Example 1 Let us consider the subgraph isomorphism problem displayed in Fig. 1.

The corresponding CSP is Different levels of consistency may be considered for these constraints, i.e., forward checking (FC) or arc consistency (AC). For AllDiff, FC is performed on the binary decomposition of the constraint while AC is hyperarc consistency on the n-ary global constraint.

X = {x 1 , x 2 , x 3 } D(x i ) = {a, b, c, d}, ∀x i ∈ X C = (x 1 , x 2 ) ∈ {(a, b), (b, a), (a, d), (d, a), (d, c), (c, d), (b, c), c, b)} ∧ (x 2 , x 3 ) ∈ {(a, b), (b, a), (a, d), (d, a), (d, c), (c, d), (b, c), c, b)} ∧ (x 3 , x 1 ) ∈ {(a, b), (b, a), (a, d), (d, a), (d, c), (c, d), (b, c), c, b)} ∧ AllDiff ({x 1 , x 2 , x 3 })
Larrosa and Valiente have proposed in [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF] a filtering algorithm (called nRF+) dedicated to the subgraph isomorphism problem. The idea is to check, for every pattern node i and every target node a ∈ D(x i ), that all variables associated with adjacent nodes of i can be assigned to different target nodes that are all adjacent to a. More precisely, they define the set

F(i, a) = ∅ if ∃j ∈ adj(i), D(x j ) ∩ adj(a) = ∅ F(i, a) = ∪ j∈adj(i) (D(x j ) ∩ adj(a)) otherwise
and remove a from D(x i ) whenever #F (i, a) < #adj(i). Larrosa and Valiente have shown that this filtering algorithm combined with binary c2 constraints and a global AllDiff constraint allows CP to efficiently solve difficult instances.

Example 2 Let us consider again the subgraph isomorphism problem displayed in Fig. 1, and let us suppose that the domains of x 2 and x 3 have been reduced to {a, b, c}. In this case, we have

F(1, a) = (D(x 2 ) ∩ adj(a)) ∪ (D(x 3 ) ∩ adj(a)) = {b}
As #adj(1) = 2, a is removed from D(x 1 ). For the same reason, c is removed from D(x 1 ). Then, b is removed from D(x 2 ) and D(x 3 ) as the neighbors of b no longer belong to the domains of neighbors of 2 and 3.

Labeling for subgraph isomorphism

This section introduces a labeling process which is used in section 4 to define a filtering algorithm for a global constraint dedicated to the SIP.

Basically, this labeling process aims at filtering, for every pattern node u, the set of target nodes that may be matched with u by a subgraph isomorphism function. To do that, the idea is to associate a label with every node and to define a partial order over these labels. This partial order defines node compatibility, i.e., a target node u may be matched to a pattern node v only if (u, v) belongs to the partial order. Of course, this partial order must be consistent with respect to subgraph isomorphism, i.e., two nodes that may be matched by a subgraph isomorphism function must be compatible. We define such subgraph isomorphism consistent labelings in 3.1. We then show in 3.2 how a labeling may be strengthened by adding information from labels of neighbors. We show in 3.3 that this strengthening step may be iterated until some fixpoint is reached, thus defining a sequence of labelings the strength of which increases at each iteration. In sections 3.1 to 3.3, we always consider undirected graphs, such that

(u, v) ∈ E ⇒ (v, u) ∈ E.
The extension of our approach to directed graphs, which is straightforward, is discussed in 3.4.

Subgraph Isomorphism Consistent Labelings

Definition 1 A labeling l is defined by a triple (L, , α) such that -L is a set of labels that may be associated with nodes; -⊆ L × L is a partial order on L; α : Node → L is a total function assigning a label α(v) to every node v.

A labeling induces a compatibility relation between nodes of the pattern graph and the target graph.

Definition 2

The set of compatible couples of nodes induced by a labeling l = (L, , α) is defined by

CC l = {(u, v) ∈ N p × N t | α(u) α(v)}
This compatibility relation may be used to filter the domain of a variable x u associated with a node u of the pattern graph by removing from it every node v of the target graph such that (u, v) ∈ CC l .

Example 3 Let us consider the subgraph isomorphism problem displayed in Fig. 2. Note that this instance has no solution as G p cannot be mapped into a subgraph of G t . Let us define the labeling l deg = (N, ≤, deg) such that deg is the function which returns node degree, i.e., ∀u ∈ N odes, deg(u) = #adj(u). This labeling assigns the following labels to nodes.

deg(A) = deg(B) = deg(D) = deg(2) = deg(4) = 4 deg(C) = deg(E) = deg(F ) = deg(G) = deg(1) = deg(3) = 3 deg(5) = deg(6) = 2
Hence, the set of compatible couples induced by this labeling is

CC l deg = {(u, v) | u ∈ {2, 4}, v ∈ {A, B, D}} ∪ {(u, v) | u ∈ {1, 3, 5, 6}, v ∈ {A, B, C, D, E, F, G}}
f (v)) ∈ CC l .
In the context of graph isomorphism, such as in [START_REF] Mckay | Practical graph isomorphism[END_REF], as opposed to subgraph isomorphism studied here, an SIC labeling is often called an invariant. In this case, the partial ordering is replaced by an equality condition: two nodes are compatible if they have the same label. Many graph properties, that are "invariant" to subgraph isomorphism, may be used to define SIC labelings such as, e.g., the three following SIC labelings:

l deg = (N, ≤, deg) where deg is the function that returns node degree; l distance k = (N, ≤, distance k ) where distance k is the function that returns the number of nodes that are reachable by a path of length smaller than k; l clique k = (N, ≤, clique k ) where clique k is the function that returns the number of cliques of size k that contains the node.

Strengthening a Labeling

We propose to start from an elementary SIC labeling which is easy to compute such as the l deg labeling defined in Ex. 3, and to iteratively strengthen this labeling. The strength of labelings is defined with respect to the induced compatible couples as follows.

Definition 4 Let l and l be two labelings. l is strictly stronger than l iff CC l ⊂ CC l . l is equivalent to l iff CC l = CC l .

A stronger labeling yields a better filtering, as it contains fewer compatible couples.

To strengthen a labeling, the idea is to extend the label of a node by adding information from the label of its neighbors. This information is a multiset (as several neighbors may have the same label). We shall use the following notations for multisets.

Definition 5 Given an underlying set A, a multiset is a function m : A → N, such that m(a) is the multiplicity (i.e., the number of occurrences) of a in m. The multiset m can also be represented by the bag {{a 0 , . . . , a 0 , a 1 , . . .}} where elements are repeated according to their multiplicity.

Example 4

The multiset m which contains exactly 2 occurrences of a, 3 occurrences of b, and 1 occurrence of c is defined by m(a)=2, m(b)=3, m(c)=1, and ∀x ∈ {a, b, c}, m(x)=0. This multiset may also be represented by {{a, a, b, b, b, c}}.

Given a partial order on a set A, we extend it to a partial order on multisets over A as follows.

Definition 6 Given two multisets m and m over a set A, and a partial order ⊆ A × A, we define m m iff there exists a total injective mapping t : m → m such that ∀a i ∈ m, a i t(a i ).

In other words, m m iff for every element of m there exists a different element of m which is greater or equal.

Example 5 If we consider the classical ordering on N, we have {{3, 3, 4}} {{2, 3, 5, 5}}, but {{3, 3, 4}} is not comparable with {{2, 5, 6}}.

Note that comparing two multisets is not trivial if the order relation on the underlying set A is not total. This point will be handled in the next section.

We now define the labeling extension procedure.

Definition 7 Given a labeling l = (L, , α), the neighborhood extension of l is the labeling l = (L , , α ) such that:

every label of L is composed of a label of L and a multiset of labels of L, i.e., L = L • (L → N); the labeling function α extends every label α(v) by the multiset of the labels of the neighbors of v, i.e., α

(v) = α(v) • m where ∀l i ∈ L, m(l i ) = #{u | (u, v) ∈ Edge ∧ α(u) = l i }; -the partial order on the extended labels of L is defined by l 1 • m 1 l 2 • m 2 iff l 1 l 2 and m 1 m 2 .
Example 6 Let us consider again the subgraph isomorphism problem displayed in Fig. 2, and the labeling l deg = (N, ≤, deg) defined in Ex. 3. The neighborhood extension of l deg is the labeling l = (L , , α ) such that labels of L are composed of numbers from N followed by multisets over N; the labeling function α is defined as follows

α (A) = 4 • {{3, 3, 4, 4}} renamed m 1 α (B) = α (D) = 4 • {{3, 3, 3, 4}} renamed m 2 α (2) = α (4) = 4 • {{2, 2, 3, 3}} renamed m 3 α (C) = 3 • {{4, 4, 4}} renamed m 4 α (E) = α (F ) = α (1) = α (3) = 3 • {{3, 4, 4}} renamed m 5 α (G) = 3 • {{3, 3, 4}} renamed m 6 α (5) = α (6) = 2 • {{4, 4}} renamed m 7 -the partial order is such that m 3 m 1 m 5 m 1 m 7 m 1 m 3 m 2 m 5 m 4 m 7 m 4 m 5 m 5 m 7 m 5
Note that we only display compatibility relationships m i m j such that m i is the label of a node of the pattern graph and m j is the label of a node of the target graph as other relations are useless for filtering purposes.

The set of compatible couples induced by this extended labeling is

CC l = {(u, v) | u ∈ {2, 4}, v ∈ {A, B, D}} ∪ {(u, v) | u ∈ {1, 3, 5, 6}, v ∈ {A, C, E, F }}
As compared to the initial labeling l deg , this set of compatible couples allows one to further remove values B, D and G from the domains of the variables associated with nodes 1, 3, 5 and 6.

The next theorem states that the neighborhood extension of a SIC labeling is a stronger (or equivalent) SIC labeling.

Theorem 1 Let l = (L, , α) be a labeling, and l = (L , , α ) be its neighborhood extension. If l is an SIC labeling, then (i) l is also SIC, and (ii) l is stronger than or equivalent to l.

Proof (i): Let f be a subgraph isomorphism function and v ∈ N p . We show that α (v) α (f (v)), that is α(v) α(f (v)) and m m , where m (resp. m ) is the multiset of the labels of the neighbors of v in G p (resp. of f (v) in G t ). Indeed, -α(v) α(f (v)) because l is SIC;
m m because f induces an injection from m to m that respects α.

(ii) : This is a direct consequence of the partial order on the extended labels in L (Definition 7) : α(u) α(v) is one of the conditions to have α (u) α (v).

Iterative Labeling Strengthening

The strengthening of a labeling described in the previous subsection can be repeated by relabeling nodes iteratively, starting from a given SIC labeling l.

Definition 8 Let l = (L, , α) be an initial SIC labeling. We define the sequence of SIC labelings l i = (L i , i , α i ) such that l 0 = l and l i+1 = neighborhood extension of l i (i ≥ 0).

A theoretical filter can be built on this sequence. Starting from an initial SIC labeling function l = l 0 , we iteratively compute l i+1 from l i and filter domains with respect to the set of compatible couples induced by l i+1 until either a domain becomes empty (thus indicating that the problem has no solution) or reaching some termination condition.

A termination condition is to stop iterating when the sequence reaches a fixpoint, i.e., a step where any further relabeling cannot change the strength of the labeling.

Theorem 2 shows that a fixpoint is reached when both the set of compatible couples and the number of different labels are not changed between two consecutive steps.

Definition 9 Given a labeling l = (L, , α), we note labels(l) the set of labels associated with nodes, i.e., labels(l) = image(α).

Theorem 2 Let l = (L, , α) be a SIC labeling. The following properties hold:

1. each iteration can only remove compatible couples and increase the number of different labels, i.e., ∀k ≥ 0, CC l k+1 ⊆ CC l k and #labels(l k+1 ) ≥ #labels(l k )

2. if at a given iteration k the set of compatible couples and the number of labels are not changed, then they will not change at any further iteration, i.e.,

∀k ≥ 0, if CC l k+1 = CC l k and #labels(l k+1 ) = #labels(l k ), then ∀j > k, CC l j = CC l k and #labels(l j ) = #labels(l k )
3. the fixpoint is reached in at most n p + n t + n p .n t steps, i.e.,

∃k ≤ n p + n t + n p .n t , CC l k+1 = CC l k and #labels(l k+1 ) = #labels(l k ) Proof (1)
The inclusion is a direct consequence of property (ii) in Theorem 1. For the cardinality of the labels, by Definition 7, we have

α k+1 (v) = α k (v).m, with m some multiset. Hence α k+1 (u) = α k+1 (v) only if α k (u) = α k (v). ( 2
)
The set of labels at step k defines a partition P k of N p ∪ N t (such that two nodes belong to a same part if they have the same label). The relation between a set S 1 of P k and a set S 2 of P k+1 must be either S 1 ∩ S 2 = ∅ or S 2 ⊆ S 1 . This follows from Definition 7: α k+1 (v) = α k (v).m implies that a node v must be in a partition α k+1 (v) finer than the partition α k (v). Hence, when #labels(l k ) = #labels(l k+1 ), the underlying partitions are the same for all subsequent steps. In particular the number of labels will never change. It remains to show that

∀ (u, v) ∈ N p × N t if (α k (u) k α k (v)) ⇔ (α k+1 (u) k+1 α k+1 (v)) and (α k (u) = α k (v)) ⇔ (α k+1 (u) = α k+1 (v)) then (α k+1 (u) k+1 α k+1 (v)) ⇔ (α k+2 (u) k+2 α k+2 (v)).
We know that the underlying partitions and the set of compatible couples are the same for step k and step k + 1. This means that there exists an isomorphism preserving between labels(l k ) and labels(l k+1 ).

(3) is a consequence of property (1) of Theorem 2; the fixpoint is reached in at most #labels(l) + #CC l steps, that is (n p + n t ) + n p .n t steps.

Example 7 Let us consider again the subgraph isomorphism problem displayed in Fig. 2, and let us suppose that the sequence of SIC labelings is started from l 0 = l deg = (N, ≤, deg) as defined in Ex. 3. After the first iteration, the neighborhood extension l 1 of l 0 is the labeling displayed in Ex. 6. From labeling l 1 , we compute the following extended labels and partial order:

α 2 (A) = m 1 • {{m 2 , m 2 , m 4 , m 5 }} renamed n 1 α 2 (B) = m 2 • {{m 1 , m 4 , m 5 , m 6 }} renamed n 2 α 2 (C) = m 4 • {{m 1 , m 2 , m 2 }} renamed n 3 α 2 (D) = m 2 • {{m 1 , m 4 , m 5 , m 5 }} renamed n 4 α 2 (E) = m 5 • {{m 1 , m 2 , m 6 }} renamed n 5 α 2 (F ) = m 5 • {{m 2 , m 2 , m 6 }} renamed n 6 α 2 (G) = m 6 • {{m 2 , m 5 , m 5 }} renamed n 7 α 2 (1) = α 2 (3) = m 5 • {{m 3 , m 3 , m 5 }} renamed n 8 2 {n 1 , n 3 } α 2 (2) = α 2 (4) = m 3 • {{m 5 , m 5 , m 7 , m 7 }} renamed n 9 2 {n 4 } α 2 (5) = α 2 (6) = m 7 • {{m 3 , m 3 }} renamed n 10 2 {n 1 , n 3 , n 5 , n 6 }
The set of compatible couples induced by this labeling is

CC l 2 = {(1, A), (1, C), (2, D), (3, A), (3, C), (4, D)} ∪ {(u, v) | u ∈ {5, 6}, v ∈ {A, C, E, F }}
From labeling l 2 , the following extended labels and partial order are computed:

α 3 (A) = n 1 • {{n 2 , n 3 , n 4 , n 5 }} renamed o 1 α 3 (B) = n 2 • {{n 1 , n 3 , n 6 , n 7 }} renamed o 2 α 3 (C) = n 3 • {{n 1 , n 2 , n 4 }} renamed o 3 α 3 (D) = n 4 • {{n 1 , n 3 , n 5 , n 6 }} renamed o 4 α 3 (E) = n 5 • {{n 1 , n 4 , n 7 }} renamed o 5 α 3 (F ) = n 6 • {{n 2 , n 4 , n 7 }} renamed o 6 α 3 (G) = n 7 • {{n 2 , n 5 , n 6 }} renamed o 7 α 3 (1) = α 3 (3) = n 8 • {{n 8 , n 9 , n 9 }} renamed o 8 3 ∅ α 3 (2) = α 3 (4) = n 9 • {{n 8 , n 8 , n 10 , n 10 }} renamed o 9 3 {o 4 } α 3 (5) = α 3 (6) = n 10 • {{n 9 , n 9 }} renamed o 10 3 ∅
The set of compatible couples induced by this labeling is

CC l 3 = {(2, D), (4, D)}
As no node is compatible with nodes 1, 3, 5 and 6, the domains of the associated variables become empty and an inconsistency is detected.

Extension to directed graphs

In directed graphs, the set of neighbors of a node u is partitioned into successors and predecessors. The iterative labeling strengthening procedure may be extended to directed graphs in a very straightforward way, by computing two different multisets that respectively contain the labels of successors and predecessors, instead of one multiset that contains the labels of neighbors. More precisely, in the case of directed graphs, we define the neighborhood extension of a labeling as follows.

Definition 10 Given a labeling l = (L, , α), the neighborhood extension of l is the labeling l = (L , , α ) such that:

every label of L is composed of a label of L and two multisets of labels of L, i.e., L = L • (L → N) • (L → N); the labeling function α extends every label α(v) by the two multisets of the labels of the predecessors and successors of v respectively, i.e., α

(v) = α(v)•m pred •m succ where ∀l i ∈ L, m pred (l i ) = #{u | (u, v) ∈ Edge ∧ α(u) = l i }, and m succ (l i ) = #{u | (v, u) ∈ Edge ∧ α(u) = l i }; -the partial order on the extended labels of L is defined by l 1 • m pred 1 • m succ 1 l 2 • m pred 2 • m succ 2 iff l 1 l 2 , m pred 1 m pred 2 , and m succ 1 m succ 2 .
This neighborhood extension procedure may be iterated until a fixpoint is reached exactly the same way as for undirected graphs.

ILF filtering algorithm

The labeling process introduced in the previous section can be used to define a filtering algorithm for a global constraint dedicated to subgraph isomorphism. Hence, we introduce in 4.1 a global constraint for the subgraph isomorphism problem. Then, we describe in 4.2 an associated filtering algorithm called ILF. We compare in 4.3 the achieved level of consistency with the ones obtained by forward checking and arc consistency on the standard CP model based on AllDiff and c2 constraints and the filtering procedure nRF+ described in 2.3.

Definition of the SIP global constraint

The SIP global constraint is defined by the relation sip(G p , G t , Lvar) where G p = (N p , E p ) and G t = (N t , E t ) are two graphs and Lvar is a set of couples which associates a different variable of the CSP with each different node of N p , i.e., Lvar is a set of #N p couples of the form (x u , u) where x u is a variable of the CSP and u is a node of N p .

Semantically, the global constraint sip(G p , G t , Lvar) is satisfied by an assignment of the variables of Lvar if there exists a subgraph isomorphism function f :

N p → N t such that, ∀(x u , u) ∈ Lvar, x u = f (u).
This global constraint is not semantically global according to [START_REF] Bessière | To be or not to be . . . a global constraint[END_REF] as it may be represented by a semantically equivalent set of AllDiff and c2 constraints as described in 2.3. However, the sip global constraint allows us to exploit the global semantic of SIPs to filter variable domains more efficiently.

Description of the ILF filtering procedure

The filtering procedure associated with the global sip constraint is called Iterative Labeling Filtering (ILF) and is described in Algorithm 1. This algorithm takes in input a global SIP constraint and initial variable domains; it is also parameterized by an initial labeling l 0 and the number k of labeling extension iterations. The different steps of this algorithm are described in sections 4.2.1 to 4.2.4. Its correctness and complexity are discussed in 4.2.5 and 4.2.6.

Filtering of D w.r.t. a labeling l i (line 3)

The successively computed labelings l i are used to filter domains by removing nodes which have incompatible labels. More precisely, for every couple (x u , u) ∈ Lvar we remove from D(x u ) every value v such that α i (u) i α i (v).

This step is done in O(n p • n t ).

Reduction of the target graph w.r.t. domains (line 5)

If a target node u does not belong to any domain, then this node and its incident edges are discarded from the target graph as no pattern node may be matched with u.

More precisely, we define N t ↓ D = N t ∩ (∪ (x u ,u)∈Lvar D(x u )). Each time domains */ for every node u ∈ Np ∪ Nt do

7 mu ← {{α i (v) | (u, v) ∈ Ep ∪ Et}} 8 α i+1 (u) ← α i (u) • mu 9 L i+1 ← {α i+1 (u)|u ∈ Np ∪ Nt} 10 rename labels in L i+1 and α i+1 11 i+1 ← {(α i+1 (u), α i+1 (v)) | (xu, u) ∈ Lvar ∧ v ∈ D(xu) ∧ test(mu, mv, i )} 12 i ← i + 1 13 return D 14 have been reduced, the target graph G t is reduced to the subgraph of G t induced by N t ↓ D.
Note that nodes and edges that are discarded during filtering must be restored when some backtracking occurs. This step is done in O(n t • d) where d is the maximal node degree.

Strengthening of a labeling l i w.r.t. singleton domains (line 6)

If a domain D(x v ) is reduced to a singleton {u}, then nodes u and v are labeled with a new label l uv which is not compatible with any other label, except itself, thus preventing other pattern nodes from being matched with u. More precisely, we add a new label l uv to L i , we assign α i (u) and α i (v) to l uv , and we add (l uv , l uv ) to i . This step is done in O(n p ).

4.2.4

Computation of the neighborhood extension l i+1 of l i (lines 7 to 12)

The neighborhood extension step is composed of the following steps:

lines 7-9: α i+1 is computed from α i . This step is done in O(e) as it basically involves collecting the labels of the neighbours of every node.

lines 10-11: Labels of L i+1 are collected and renamed in order to keep their size in O(d).

This step is done in O(d • n) as there are at most n labels the length of which is bounded by d. We use a hashtable to associate new names to labels: each time a label is computed, if it already belongs to the hashtable then the label is renamed consequently, otherwise a new name is created and the label is added to the hashtable.

Let mu = {{a, a, b, c}} and mv = {{x, y, y, z, z}} be two multisets to be compared, and let the partial order be = {(a, x), (a, y), (b, x), (c, z)}. To decide if test(mu, mv, ) is true or false, we build the following bipartite graph and use Hopcroft algorithm to find a matching which covers mu. On this example, {(a, y), (a, y), (b, x), (c, z)} is a matching that covers mu so that test(mu, mv, ) returns true.

Fig. 3 Modeling the compatibility test problem as a matching problem.

line 12: The partial order i+1 is computed.

For every couple of nodes (u, v) ∈ N p × N t such that v was compatible with u at step i, we call the function test(m u , m v , i ) to determine if m u i m v , i.e., if there exists, for each label occurrence in m u , a distinct label occurrence in m v which is greater or equal according to i . As different nodes may have the same label, we call the test function only when encountering a couple (m u , m v ) for the first time, and we memorize the result in a table so that we can reuse it if the same couple is encountered later.

To implement the test function, we have used the matching algorithm of Hopcroft [START_REF] Hopcroft | An n 5/2 algorithm for maximum matchings in bipartite graphs[END_REF]. Indeed, m u i m v iff there exists a matching that covers m u in the bipartite graph G = ((m u , m v ), E) such that E = {(x, y) ∈ m u × m v | x i y}, i.e., a label occurrence x of m u is linked by an edge with a label occurrence y of m v iff x i y (see example in Fig. 3). The complexity of this algorithm is in

O(#m u • #m v • √ #m u + #m v ).
As the sizes of m u and m v are bounded by the maximal degree d, the test function can be executed in O(d 5/2 ). Hence, the computation of the partial order is done in O(n p • n t • d 5/2 ).

Correctness of the filtering procedure

Property 1 Let l 0 be an SIC labeling and k ∈ N. The ILF(sip(G p , G t , Lvar), D, l 0 , k) procedure is correct, i.e., for every subgraph isomorphism function f from N p to N t such that ∀(x u , u) ∈ Lvar, f (u) ∈ D(x u ), ILF does not remove f (u) from D(x u ).

Proof. Indeed, we have shown in section 3 that the iterative labeling strengthening procedure is correct, i.e., if it is started from an SIC labeling l 0 , then every successive l i labeling is also an SIC labeling that can be used to filter domains without removing solutions. Our filtering procedure is based on this iterative labeling strengthening procedure, but it also integrates information brought by variable domains to filter more values:

the target graph G t is reduced to the subgraph induced by variable domains; labelings are strengthened with respect to singleton domains.

It is clear that this extra filtering does not remove values that may belong to a solution.

Complexity of the filtering procedure

At each iteration, the most expensive step is the computation of the partial order during the neighborhood extension procedure, which is done in O(n p • n t • d 5/2 ). The procedure stops iterating either if the bound k on the number of iterations is reached, or if a domain becomes empty, or if the fixpoint is reached. We have shown in 3.3 that the fixpoint is reached in O(n p .n t ) iterations.

Property 2

The ILF procedure has an overall time complexity of

O(min(k, n p • n t ) • n p • n t • d 5/2 )
We propose in Section 5 an approximated filtering procedure which ensures a lower partial consistency at a lower cost.

Level of consistency of ILF

The level of consistency achieved by the ILF filtering procedure described in the previous section depends on the initial labeling l 0 which is a parameter of this procedure.

We first introduce and compare in 4.3.1 different SIC initial labelings that may be considered. We then compare the level of consistency ensured by ILF with forward checking (FC) and arc consistency (AC) on the standard SIP model and the filtering algorithm nRF+ introduced in section 2.3.

Initial SIC labeling l 0

The iterated neighborhood extension process is started from an initial SIC labeling l 0 which is a parameter. Different initial labelings may be considered, that may have different strength.

The weakest initial labeling that may be used is the labeling l ∅ which associates the same label ∅ to every node, i.e., l ∅ = (L ∅ , α ∅ , ∅ ) such that -L ∅ = {∅}, -∀u ∈ N odes, α ∅ (u) = ∅, and -∅ = {(∅, ∅)}.

A stronger initial labeling is the l deg labeling defined in Ex. 3. Note however that, when l 0 = l ∅ , the labeling l 1 computed after the first iteration of ILF exactly corresponds to l deg as the label α 1 (u) of every node u is ∅.m u with m u (∅) = deg(u). Hence, ILF(sip, D, l deg , k) achieves the same level of consistency as ILF(sip, D, l ∅ , k + 1).

One may define stronger initial labelings by integrating information brought by the domains. Let us define for example the labeling l dom = (L dom , α dom , dom ) such that -L dom = N, α dom assigns a different unique label to every node, and

-∀(u, v) ∈ N p × N t , α dom (u) dom α dom (v) iff v ∈ D(x u ) and deg(u) ≤ deg(v).
When every variable domain contains all the target nodes, l dom is equivalent to l deg . However, if some target nodes have been removed from some variable domains, then l dom is strictly stronger than l deg so that ILF(sip, D, l dom , k) will achieve a stronger consistency than ILF(sip, D, l deg , k).

Note however that l dom introduces many different labels and that α dom is not a total order so that the neighborhood extension procedure may become expensive as most of the time is spent for comparing multisets and the complexity of this step depends on the number of different labels (see 4.2.4).

A good compromise to limit the number of different initial labels while strengthening the initial labeling with respect to information brought by the domains may be obtained with the labeling l domCC = (L domCC , α domCC , domCC ) such that -L domCC = N.N, -∀u ∈ N odes, α domCC (u) = deg(u).id(u) where id(u) is defined as follows. Let G = ((N p , N t ), E) be the bipartite graph such that

E = {(u, v) ∈ N p × N t | v ∈ D(x u )}.
We compute the different connected components of G, and we associate a different id with every different connected component of G. We note id(u) the id of the connected component which contains u. -domCC is the set of couples (deg

1 .id 1 , deg 2 .id 2 ) such that deg 1 ≤ deg 2 and id 1 = id 2 .
l domCC is weaker than l dom , but it introduces fewer labels. When the bipartite graph defined by domains is connected, l domCC is equivalent to l deg . However, if it is not connected, then l domCC is strictly stronger than l deg so that ILF(sip, D, l domCC , k) will achieve a stronger consistency than ILF(sip, D, l deg , k).

Comparison with forward checking on c2 and AllDiff

Our filtering procedure ILF achieves a level of consistency at least as strong as forward checking (FC) on the standard SIP model which combines a global AllDiff constraint with c2 constraints, provided that the number of strengthening iterations (defined by the parameter k) is greater or equal to 1, and whatever the initial labeling l 0 is. In this case, our filtering procedure can be used as a unique filter for SIP, thus offering a global constraint for SIP.

Property 3 ILF( , , ,1) is as strong as F C(AllDiff + c2).

Proof Each time a variable x u associated with a pattern node u is assigned to a target node a, FC filters the domain of every variable x v with respect to c2 and AllDiff constraints as follows:

-For every pattern node v ∈ adj(u), FC of c2(x u , x v ) removes any value b ∈ D(x v ) such that (a, b) is not an edge.

ILF( , , ,1) also removes these values. Indeed, during the first iteration, the same unique label l ua is assigned to nodes u and a (line 6) before computing the extended labeling l 1 . Hence, after the first iteration, we have α 1 (v) 1 α 1 (b) for every node b such that (a, b) is not an edge. Indeed, α

1 (v) = α 0 (v).m v with l ua ∈ m v (as u ∈ adj(v)) whereas α 1 (b) = α 0 (b).m b with l ua ∈ m b (as a ∈ adj(b)) so that m v 0 m b and therefore α 1 (v) 1 α 1 (b).
-FC of AllDiff removes a from the domains of every variable associated with a pattern node v = u. ILF( , , ,1) also removes these values. Indeed, because of the unique label l ua assigned to nodes u and a, the label of every pattern node v = u is not compatible with the label of a so that a is removed from D(x v ) by ILF during the first iteration.

Comparison with arc consistency on c2 and AllDiff

Our filtering procedure ILF is not comparable with arc consistency (AC) on the standard SIP model which combines a global AllDiff constraint with c2 constraints.

Indeed, ILF is able to prove the inconsistency of the SIP instance of Figure 2, as shown in Ex. 7. This is still true even if the initial labeling is l ∅ . As a comparison, AC(AllDiff + c2) is not able to prove the inconsistency and does not reduce any domain.

However, AC(AllDiff +c2) may be stronger than ILF in some cases. Let Example 8 Let us consider two pattern nodes u and v that are linked by an edge, and let us consider the case where the domains of the two variables x u and x v respectively associated with u and v is such that none of the target nodes of D(x u ) and D(x v ) are linked by an edge, i.e., (D(x u ) × D(x v )) ∩ E t = ∅. In this case, the propagation of AC(c2) removes all values from D(x u ) and D(x v ) since each value of D(x u ) (resp. D(x v )) has no support in D(x v ) (resp. D(x u )).

Note that ILF may also be able to reduce the domains of x u and x v to empty sets if the initial labeling l 0 takes into account information brought by domains. This is the case, e.g., if the initial labeling is the l dom labeling introduced in 4.3.1. Indeed, in this case, a different label l u is associated with every different node u, and we have l u 0 l a iff a ∈ D(x u ). Hence, in our example, l u will be compatible only with the set of labels {l a | a ∈ D(x u )} whereas l v will be compatible only with the set of labels {l b | b ∈ D(x v )}. After the first labeling extension, we will have α 1 (u) = l u .m u such that l v ∈ m u as v ∈ adj(u), and for every target node a ∈ D(x u ), α 1 (a) = l a .m a such that ∀b ∈ D(x v ), l b ∈ m a as b ∈ adj(a).

Hence, ∀a ∈ D(x u ), m u 0 m a as m a only contains labels which are not compatible with the label l v which belongs to m u . Therefore, ∀a ∈ D(x u ), α 1 (u) 1 α 1 (a), so that the domain of x u will be reduced to the empty set after one iteration of ILF.

Example 9 Let us consider the case where a subset of n variables associated with pattern nodes all have the same domain D such that #D < n. In this case, the instance has no solution as the subgraph isomorphism function must be injective. This will be detected by AC(AllDiff ) as it trivially detects inconsistency when the number of values is lower than the number of variables.

Note that this kind of inconsistency cannot be detected by ILF, whatever the initial labeling is.

Comparison with nRF+

When the initial labeling is l deg , our filtering procedure ILF is not comparable with the filtering procedure nRF+ proposed by Larrosa and Valiente in [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF]. Indeed, ILF is able to prove the inconsistency of the SIP instance of Fig. 2, as shown in Ex. 7, whereas nRF+ combined with AC(AllDiff + c2) is not able to prove the inconsistency: it is only able to reduce the domains of the variables associated with nodes 2 and 4 to {A, B, D} whereas the domains of the other variables are not reduced. However, ILF is not able to reduce any domain for the SIP instance of Fig. 1 whereas nRF+ is able to reduce some domains as shown in Example 1.

When the initial labeling is l dom , our filtering procedure ILF is at least as strong as nRF+, provided that the number of iterations k is greater or equal to 1.

Property 4 ILF( , ,l dom ,1) is as strong as nRF+.

Proof. nRF+ removes from the domain of a variable associated with a pattern node i every target node a such that #F (i, a) < #adj(i). This value is also removed after the first iteration of ILF. Indeed, let us consider the multisets m i and m a which contain the labels of the neighbours of i and a. When the initial labeling is l dom , the label in m i of a node j adjacent to i is compatible with the label in m a of a node b adjacent to a only if b ∈ D(x j ). Therefore, if #F (i, a) < #adj(i), m i will not be compatible with m a and a will be removed from D(x i ).

Weaker filtering algorithm ILF*

At each iteration of ILF, the most expensive step is the computation of the partial order i+1 during the neighborhood extension procedure, which is done in O(n p • n t • d 5/2 ). This comes from the fact that each call to the test function is in O(d 5/2 ).

However, if i is a total order (instead of a partial order), the function test(m u , m v , i ) may be implemented much more efficiently, by sorting each multiset and matching every label of m u with the smallest compatible label of m v . In this case, the complexity of test is O(d • log d), and the complexity of the neighborhood extension becomes O(n p • n t • d) as the sorting of multisets may be done once. When i is not a total order, one may extend it into a total order ≤. This total order may then be used in the test function to determine if m u ≤ m v . However, the total order introduces new label compatibilities so that test(m u , m v , ≤) may return true while test(m u , m v , i ) returns false. As a consequence, using this approximated order may induce a weaker filtering.

In this section, we first introduce in 5.1 the theoretical framework that defines a new neighborhood labeling extension based on a total order and proves its validity; then we show in 5.2 how to compute the total order; finally, we define in 5.3 the new filtering algorithm called ILF * and we discuss its complexity and the level of consistency it achieves.

Labeling strengthening based on total orders

We define a new neighborhood extension procedure where old labels are extended by multisets of neighbor labels like in Definition 7, but where the partial order over the new extended labels is defined with respect to a total order over the old labels, thus reducing the complexity of its computation.

The next definition gives a simple condition on the total order to ensure its consistency with respect to the partial order, i.e., to ensure that if test(m u , m v , i ) returns true, then test(m u , m v , ≤) also returns true.

Definition 11 Let l = (L, , α) be a labeling. A consistent total order for l is a total order ≤ on L such that ∀u ∈

N p , ∀v ∈ N t , α(u) α(v) ⇒ α(u) ≤ α(v)
We extend the order ≤ on multisets like for partial orders in Definition 6, i.e., m ≤ m iff there exists an injective function t : m → m such that ∀a i ∈ m, a i ≤ t(a i ). Hence, m m ⇒ m ≤ m . Let us note however that this extension of ≤ to multisets only induces a partial order on multisets as some multisets may not be comparable.

We can then define a new neighborhood extension procedure, based on a consistent total order. Definition 12 Let l = (L, , α) be a labeling, and let ≤ be a consistent total order for l. The neighborhood extension of l based on ≤ is the labeling l = (L , , α ) where L and α are defined like in Definition 7, and the order relation ⊆ L × L is defined by

l 1 • m 1 l 2 • m 2 iff l 1 l 2 ∧ m 1 ≤ m 2
The next theorem shows that a neighborhood extension based on a consistent total order may be used in our iterative labeling process, and that it is stronger or equal to the original labeling. However, it may be weaker than the neighborhood extension based on the partial order. Indeed, the total order induces more compatible couples of labels than the partial order.

Theorem 3 Let l = (L, , α), l = (L , , α ), and l = (L , , α ), be three labelings such that l is the neighborhood extension of l and l is the neighborhood extension of l based on a total consistent order.

If l is an SIC labeling, then (i) l is SIC, (ii) l is stronger than (or equal to) l, and (iii) l is stronger than (or equal to) l . Proof (ii) and (iii): For labeling l , we have l

1 • m 1 l 2 • m 2 iff l 1 l 2 ∧ m 1 m 2 .
As ≤ is an extension of , we have m m ⇒ m ≤ m . Hence, CC l ⊆ CC l ⊆ CC l . (i) is a direct consequence of (iii), as l is SIC (Theorem 1).

Computation of a total order

Given a partial order ⊆ L × L, one may easily compute a consistent total order by performing a topological sort of the graph G = (L, ) [START_REF] Cormen | Introduction to Algorithms[END_REF]. Topological sorts are not unique so that different consistent total orders may be derived from a given partial order, and these different total orders may induce prunings of different strength: the less new couples of compatible nodes are introduced by the total order, the better the filtering. However, we conjecture that finding the best consistent total order is NP-hard (see appendix A).

We propose a heuristic algorithm that aims at computing a total order that introduces few new compatible couples without guarantee of optimality. This algorithm is based on the following greedy principle: starting from an empty sequence, one iteratively adds some labels at the end of the sequence, until all labels have been sequenced.

Algorithm 2: Computation of a consistent total order

Input:

the set Lp of labels of pattern nodes the set Lt of labels of target nodes a partial order ⊆ Lp × Lt Output: an increasing sequence S defining a consistent total order for l for every label e ∈ Lp ∩ Lt do The total order is derived from the sequence by stating that a label is smaller than every label occuring after it in the sequence. This algorithm is described in algorithm 2. It first renames the labels of the target graph that are also used in the pattern graph (lines 1-4). This does not modify the set of compatible couples of nodes, induced by the labeling. This renaming step is done in

O(n p • n t ).
Then, at each iteration of our greedy algorithm, we choose a label e t of the target graph which has not yet been sequenced in S (lines 7-9) and then we introduce in the sequence the selected label e t , preceded by every label e p ∈ L p such that e p e t (lines 10-12) in order to preserve the compatibility relation.

Our goal is to sequence the labels of L p as late as possible (in order to minimize the number of target nodes which are compatible with pattern nodes). To this aim, we first compute (line 7) the set of labels e t ∈ L t for which the number of labels e p ∈ L p , e p e t is minimal. To break ties, we then select (line 8) the labels which are already compatible with the maximum number of labels in L t .

The time complexity of this heuristic algorithm is in O(n t • log n t • d p • d t ), which for simplicity will be approximated by Õ(n p • n t • d).

Complexity of ILF *

We denote by ILF * the filtering procedure obtained when calling algorithm 2 before computing the new partial ordering i+1 , i.e., before line 12 of algorithm 1, and replacing the current partial order by the computed total order in the call to the test function.

Property 5

The ILF * procedure has an overall time complexity of

Õ(min(k, n p • n t ) • n p • n t • d)

Considered instances

We evaluate our approach on scale-free networks and on graphs of the GraphBase database used in [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF].

Scale-free networks (classes A to F) Scale-free networks model a wide range of real networks, such as social, Internet, or neural networks [START_REF] Barabasi | Linked: How Everything Is Connected to Everything Else and What It Means[END_REF]. We have randomly generated graphs using a power law distribution of degrees P (d = k) = k -λ . We have made experiments with different values of λ, ranging between 1 and 5, and obtained similar results. Hence, we only report experiments on graphs generated with the standard value λ = 2.5.

We have considered 6 classes of instances, each class containing 20 different instances. These classes are described in Table 1. For each instance, we first generate a connected target graph which node degrees are bounded between d min and d max .

Then, a connected pattern graph is extracted from the target graph by randomly selecting a percentage p n (resp. p e ) of nodes (resp. edges). We have performed experiments with different percentages p n and p e and noticed that when p n and p e are smaller than 0.7 the number of solutions increases so much that instances become much easier to solve. Hence, we only report results obtained on harder instances with p n = 0.9.

All instances of classes A to D are sparse feasible instances such that node degrees range between 5 and 8. Target graphs in classes A to C are undirected and their number of nodes is increased from 200 to 1000 to study scale-up properties. Target graphs in class D are directed, such that their edges have been randomly directed.

All instances of classes E and F are denser undirected instances that have been generated with d min = 20, d max = 300. Instances of class E are feasible ones that have been generated with p e = 90%. Instances of class F are non feasible ones: for these instances, pattern graphs are extracted from target graphs by randomly selecting 90% of nodes and 90% of edges, but after this extraction, 10% of new edges have been randomly added.

For all experiments reported below, each run searches for all solutions of an instance.

GraphBase database (classes G and H)

This database contains graphs with different properties, i.e., simple, acyclic, connected, biconnected, triconnected, bipartite and planar (see [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF]). It contains 113 undirected graphs and 59 directed graphs. We have considered the first 30 directed graphs. This set contains graphs ranging from 10 nodes to 462 nodes. Using these graphs, we have generated 406 directed instances (class and CP(Map) [START_REF] Dooms | Cp(graph): Introducing a graph computation domain in constraint programming[END_REF][START_REF] Deville | Cp(graph+map) for approximate graph matching, 1st International Workshop on Constraint Programming Beyond Finite Integer Domains[END_REF] which provide graph and function domain variables. These new structured domains offer abstractions that simplified the implementation of the filtering algorithms. In all experiments, the filtering algorithms ILF and ILF* have been combined with AC(c2) and AC(AllDiff). Table 2 compares ILF(k) and ILF*(k), respectively based on an exact partial order or an approximated total order (as described in 4.2 and 5), for different limits k on the number of iterations. In all variants, the initial labeling l 0 is the labeling l deg defined in 2.3. Note that the order of l deg is a total order so that in this case the exact and approximated variants are equivalent for k = 1.

Let us first compare ILF and ILF*. The number of failed nodes with ILF*(2) is greater than with ILF(2), but it is smaller than with ILF(1). This shows us that the total order computed by our heuristic algorithm is a quite good approximation of the partial order. When considering CPU-times, we note that ILF*(2) is significantly quicker than ILF [START_REF] Bessière | To be or not to be . . . a global constraint[END_REF].

Table 2 also shows that the best performing variant differs when considering different classes of instances. The directed instances of classes D and G are best solved when k = 0, i.e., with the simple l deg labeling. Instances of classes A, B, C, E, F and H are undirected; they are best solved when k = 1 than when k = 0, i.e., after one iteration of the labeling extension. Instances of class F, which have significantly higher node degrees and are infeasible, are very difficult ones. For these instances, iterative labeling extensions actually improve the solution process and the best results are obtained with ILF * (4).

As a conclusion, these experiments show us that (1) the approximated variant offers a good compromise between filtering's strength and time, and (2) the optimal limit k on the number of iterations depends on the difficulty of instances. Good average results are obtained with ILF*(2). Table 3 Comparison of state-of-the-art approaches. Each line successively reports the percentage of instances that have been solved within a CPU time limit of 600s on an Intel Xeon 3,06 Ghz with 2Go of RAM; the average run time for the solved instances; and the average number of failed nodes in the search tree for the solved instances.

Comparison with state-of-the-art approaches

We now compare results obtained with our filtering procedure ILF*(2), which offers a good compromise between search space reduction and time, with other CP models. We consider four different CP models which combine different kinds of constraints and different levels of consistency as introduced in Section 2.3:

-CP1 performs forward checking on c2 and AllDiff constraints; -CP2 maintains arc consistency on c2 and AllDiff constraints; -CP3 combines forward checking on c2 and AllDiff constraints with the filtering procedure nRF+ of [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF];

-CP4 combines arc consistency on c2 and AllDiff constraints with the filtering procedure nRF+ of [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF].

All CP models have been implemented in Gecode using CP(Graph) and CP(Map). We compare these different CP models with a state-of-the-art algorithm coming from a C++ library vflib.

Conclusion

We introduced a new filtering algorithm for the subgraph isomorphism problem that exploits the global structure of the graph in order to achieve a stronger partial consistency. This work extends a filtering algorithm for graph isomorphism [START_REF] Sorlin | A new filtering algorithm for the graph isomorphism problem[END_REF] where a total order defined on some graph property labeling is strengthened until a fixpoint is reached. The extension to subgraph isomorphism has been theoretically funded. The order is partial and can also be iterated until a fixpoint is reached. However, using such a partial order is very time consuming. Instead, one can map this partial order to a total order. Performing such a mapping is hard, and can be efficiently approxi-mated through a heuristic algorithm. Experimental results show that our propagators are efficient against state-of-the-art propagators and algorithms. One of the advantages of our CP approach is to easily accomodate extra constraints. Further work includes the comparisons of the different CP approaches, including our ILF and ILF* algorithms on different classes of SI problems with side constraints. We will also extend our filters to better take into account the domain of the variables. Although Gecode has a copy strategy, it is worth introducing incrementality in our filtering algorithms, reducing computation between different executions of the propagator. Future work will also include the development of dynamic termination criteria for the iterative labeling, the experimental study of other degree distributions, and the analysis of alternative initial labelings.

Fig. 1

 1 Fig. 1 Instance of subgraph isomorphism problem.

Fig. 2

 2 Fig. 2 Instance of subgraph isomorphism problem.

Algorithm 1 : 1 loop 2 filter 4 reduce

 1124 ILF(sip(G p , G t , Lvar), D, l 0 , k)Input:a global SIP constraint sip(Gp, Gt, Lvar) the initial domain D(xu) of every variable xu such that (xu, u) ∈ Lvar an initial SIC labeling l 0 = (L 0 , 0 , α 0 ) a limit k on the number of iterations Output: filtered domains i ← 0 D w.r.t. l i 3 exit when ∃(xu, u) ∈ Lvar, D(xu) = ∅ or i = k or fixpoint reached Gt w.r.t. D 5 strengthen l i w.r.t. singleton domains 6 /* Computation of the neighborhood extension l i+1 of l i

1 3 add (e , e ) to 4 S ←<> 5 while Lt = ∅ do 6 Cand 7 Cand 8 randomly choose et in Cand 9 for

 3456789 replace e by e in Lt (where e is a new label such that e ∈ Lp ∪ Lt) 2 for every label e ∈ Lp such that e e do ← {et ∈ Lt | #{ep ∈ Lp, ep et} is minimal} ← {et ∈ Cand | ep∈Lp∧ep e t #{e ∈ Lt, ep e} is maximal} every label ep ∈ Lp such that ep et do 10 add ep at the end of S and remove ep from Lp 11 add et at the end of S and remove et from Lt 12 add the remaining elements of Lp at the end of S 13 return S 14

  Definition 3 A labeling l is subgraph isomorphism consistent (SIC) iff for any subgraph isomorphism function f , we have ∀v ∈ N p , (v,

	This set of compatible couples allows one to remove values C, E, F and G from the
	domains of the variables associated with nodes 2 and 4.
	The goal of this work is to find a labeling that filters domains as strongly as possible
	without removing solutions to the subgraph isomorphism problem, i.e., if a node u
	of the pattern graph may be matched to a node v of the target graph by a subgraph
	isomorphism function, then the label of u must be compatible with the label of v. This
	property is called subgraph isomorphism consistency.

  us consider for example the case where both the pattern and the target graphs are regular graphs such that every pattern (resp. target) node has the same degree d p (resp. d t ), with d p ≤ d t . If the initial labeling is l deg , all pattern (resp. target) nodes have the same label d p (resp. d t ) and the fixpoint is reached at the first iteration without reducing any domain as d p ≤ d t . However, AC may be able to reduce some domains in some cases as shown in Ex. 8 for AC(c2) and Ex. 9 for AC(AllDiff ).

Table 1

 1 Description of the considered classes of instances

		Target graph	Pattern graph	directed	feasible
	Class #Nt	d min	dmax	pn	pe	graphs ?	instances ?
	A	200	5	8	90%	90%	no	yes
	B	600	5	8	90%	90%	no	yes
	C	1000	5	8	90%	90%	no	yes
	D	600	5	8	90%	90%	yes	yes
	E	300	20	300 90%	90%	no	yes
	F	300	20	300 90% 90% + 10%	no	no
	6 Experimental Results

Table 2

 2 Comparison of different variants of ILF(k) and ILF*(k) where k gives the maximum number of iterations. Each line successively reports the percentage of instances that have been solved within a CPU time limit of 600s on an Intel Xeon 3,06 Ghz with 2Gb of RAM; the average run time for the solved instances; and the average number of failed nodes in the search tree for the solved instances.G) and 1276 undirected instances (class H) of the subgraph isomorphism problem by considering all couples of graphs (G p , G t ) such that N p < N t .

	9 14.8	-16.7 18.6 29.5	144.3 42.0	-29.1 19.9 17.6

6.2 Comparison of ILF(k) and ILF*(k)

Algorithm 1 has been implemented in Gecode (http://www.gecode.org), using CP(Graph)

Table 3

 3 compares all these approaches and shows us that CP approaches clearly outperform vflib, except for directed instances of class D. When comparing the different CP approaches, we note that the number of failed nodes is reduced when considering higher consistencies, and ILF*(2) always explores much less nodes than other CP models. This allows ILF*(2) to solve more instances for all classes, except for instances of class H which are best solved by CP2.

Acknowledgments

The authors want to thank the anonymous reviewers for their helpful comments. Christine Solnon acknowledges an Anr grant Blanc 07-1 184534: this work was done in the context of project Sattic. This research is also partially supported by the Interuniversity Attraction Poles Programme (Belgian State, Belgian Science Policy).

Appendix A: Complexity of the best total order problem Given a partial order ⊆ N p × N t , our best total order problem involves finding the consistent total order ≤ which introduces the smallest number of new compatibility relationships from N p to N t , i.e., such that #{(u, v) ∈ N p × N t , u ≤ v} is minimal.

Our problem may be formalized as follows. Let us define the directed bipartite graph between N p and N t which contains an arc from u ∈ N p to v ∈ N t if α(u) α(v), and an arc from v ∈ N t to u ∈ N p if ¬(α(u) α(v)). The underlying undirected bipartite graph is complete so that it corresponds a total relation. However, it may not correspond to a total order as it may contain cycles (so that transitivity is not verified). The best total order problem is equivalent to finding the smallest number of arcs from N t to N p that should be reversed so that the graph becomes acyclic.

We conjecture that this problem is NP-hard. Indeed, the feedback arc set problem, which involves finding the smallest number of arcs that should be deleted to make a directed graph acyclic has been shown to be NP-hard, even in the case of complete bipartite graphs [START_REF] Guo | Feedback arc set in bipartite tournaments is np-complete[END_REF]. The only difference between the feedback arc set problem and our best total order problem is that we must not reverse (or delete) arcs from N p to N t (otherwise the total order would not be consistent).


Frederik Kurzrock

Sylvain Cros²

Nicolas Sébastien

Fabrice Chane-Ming

Laurent Linguet

Roland Potthast

Gilles Lajoie

The use of satellite data assimilation methods in regional NWP for solar irradiance forecasting

In common state-of-the-art data assimilation procedures, primarily applied to global NWP models, all types of available observations are assimilated using the established methods (e.g. 3D-Var, 4D-Var, ENKF). The ultimate goal is to find the optimal initial conditions regarding the model's state variables. The typically assimilated satellite observations are radiances and atmospheric motion vectors (AMVs). Assimilating such observations using variational data assimilation or Ensemble Kalman Filter methods affects the model state variables in the first place and therefore influences the representation and evolution of clouds indirectly.

Advantages of satellite observations

Recent developments and their shortcomings

Compared to other sources of observations, geostationary satellites provide many advantages regarding continuous solar irradiance forecasting:

• large spatial coverage ; high spatial and temporal resolution • lack of measured data around islands • satellite radiances are reliable radiometric signatures of cloud presence, properties and evolution • constant accuracy in space and time • comprehensive observations for large parts of the world • future-oriented

In the recent years, different approaches aimed at deriving atmospheric analyses for limited-area models with most realistic cloud features, using geostationary satellite observations. The focus is mostly set on the assimilation of cloud-top information (e.g. cloud top temperature), since optical and thermal sensors are not able to capture information inside clouds (Figure 2). Table 1 sums up recent advances and their main limitations. Although the listed methods generally provide improved intraday (and even day-ahead) cloud cover forecasts in mid-latitudes, a major limitation for all methods is that only clear-sky or completely cloudy cases can be considered. This is because fractional clouds cause a measured signal mixing cold clouds and the warmer Earth surface. 1. Motivation

Indirectly influencing cloud evolution

Photovoltaic energy offers an intermittent production. PV production forecasts are necessary to permit a safe integration into the electricity grids. The correct simulation of cloud evolution is essential to forecast accurately irradiance reaching the ground level. Besides solar power forecasting, accurate short-term forecasts of precipitating and non-precipitating clouds are important for various meteorological applications. Regional NWP is the state-of-the-art tool to provide highly resolved cloud forecasts for lead times beyond 6h. Clouds are currently rather poorly represented in NWP models. The precise simulation of cloud evolution requires accurate atmospheric analyses which are obtained by exploiting atmospheric observations using data assimilation methods. In this work we focus on the assimilation of geostationary satellite data into regional NWP models with the goal of improving short-term cloudiness forecasts at high spatial resolutions. Progress being made regarding the assimilation of cloud-affected radiances into limited-area NWP models (Stengel et al., 2013). However, the indirect assimilation of quantities like radiances or AMVs has not been exploited so far in regional NWP regarding solar irradiance forecasting. The main reason for this is a number of problems that have to be overcome for each individual model domain: the choice of satellite channels, optimal thinning, bias correction, observation and background error estimation. Future investigations address:

Satellite-derived parameters

Influenced

• the combination of both the classical "indirect" assimilation methods with methods that explicitly assimilate cloud physical properties (e.g. cloud top height, cloud emissivity, …).

• the impact of such methods on different lead times, since the information from the driving global model is expected to predominate at a certain point.

• the influence of domain size and grid spacing in connection with satellite data assimilation in regional NWP.

• the use of different limited-area models and satellites, and respectively the application to different geographical regions (e.g. mid-latitudes and tropics). The ultimate goal is an assimilation strategy that can be used independently of the model and satellite and which can

Outline

Figure 1 :

 1 Figure 1: Mashup of a MSG-2 RGB composite (source: EUMETSAT) and a GFS GHI forecast (visualised using Panoply).

Figure 2 :

 2 Figure 2: Sketches of arbitrary satellite IR weighting functions for different cases. (After lecture notes of Tony McNally (ECMWF))


Frederik Kurzrock
email: frederik.kurzrock@univ-reunion.fr

Sylvain Cros

Fabrice Chane-Ming

Roland Potthast

Laurent Linguet

Gilles Lajoie

SATELLITE DATA ASSIMILATION IN REGIONAL NUMERICAL WEATHER PREDICTION AS A KEY FOR BETTER CLOUD COVER FORECASTS IN TROPICAL ENVIRONMENTS

Keywords: Forecasting, Solar Radiation, Grid-Connected, Plant Control, PV System

Although the high amount of solar irradiance in the tropics is an advantage for a profitable PV production, the local meteorological conditions induce a very high variability which is problematic for a safe and gainful injection into the power grid. This issue is even more critical in non-interconnected territories where network stability is an absolute necessity and the injection of PV power has to be limited. The basis for precise cloud evolution and subsequent irradiance forecasts are high quality atmospheric analyses for NWP. Geostationary meteorological satellites provide valuable observations of cloud properties with high spatiotemporal resolutions and allow a pertinent data assimilation. The shortcoming is that optical and thermal channels of satellite sensors do not provide cloud properties from inside clouds. Different existing data assimilation approaches aim at deriving atmospheric analyses with most realistic cloud features, utilising geostationary satellite observations. The potential of assimilating satellite-derived cloud information in regional NWP with focus on irradiance forecasts in tropical regions has not been evaluated so far. Hence, the present work aims at evaluating the potential of geostationary satellite data assimilation in limited-area models applied to the French tropical oversea territories Reunion Island and French Guiana.

INTRODUCTION

The globally installed amount of PV systems is expected to grow from 178 GW in 2014 to 540 GW by 2020 ( [START_REF] Rekinger | Global Market Outlook for Solar Power 2015-2019[END_REF]). This implies that solar power is becoming a major component of the future global energy mix and opens the way towards independence of fossil fuels and thus a reduction of global warming. Nevertheless, the exploitation of solar power is a technological challenge. As the photovoltaic power production reduces in the presence of clouds, its intermittent character constitutes an obstacle concerning grid injection. In France, the injection of intermittent renewable energies is legally capped at 30 %. This threshold has already been exceeded in Corsica (40 %) and in the French tropical oversea territories Guadeloupe (35 %), Reunion Island (39 %). In French Guiana it already reached 25 % [START_REF][END_REF].

Reliable forecasts of the solar power production facilitate a safer injection of PV power into the grids and will in the long term permit the elevation of legal injection limits. The high amount of irradiance in the tropics, combined with the fact that non-interconnected territories are obliged to be energy independent, leads to the importance of accurate irradiance forecasts in such regions if one wants to increase the injection of PV into the grids [3] [4].

Regarding the requirements for temporal (a few minutes) and spatial (a few kilometers) resolutions of day-ahead solar power forecasts, numerical weather prediction (NWP) is the only appropriate forecasting solution. Solar power production forecasts almost entirely depend upon the forecast of the Global Horizontal Irradiance (GHI), directly driven by cloud conditions. Being able to predict cloud cover and its evolution thus means being able to predict solar power production. Since many physical processes which have large impacts on cloud and radiation processes have to be parameterised in NWP models, irradiance predictions by NWP models exhibit a considerable lack of precision regarding solar energy applications [START_REF] Kleissl | Solar Energy Forecasting and Resource Assessment[END_REF] [START_REF] Sengupta | Best Practices Handbook for the Collection and Use of Solar Resource Data for Solar Energy Applications[END_REF].

The temporal and spatial resolutions of global weather models (around 0.25° grid spacing and 3h temporal resolution) do not permit high-precision cloud evolution forecasts. It is therefore necessary to apply dynamical downscaling with regional NWP models in order to increase the spatio-temporal forecast resolution and permit a more realistic simulation of cloud processes.

Irradiance prediction in the tropics is particularly difficult due to pronounced convection, homogenous air masses and large temperature contrasts between land and sea [START_REF] Beucher | Météorologie Tropicale: des alizés au cyclone[END_REF].

CURRENT STATUS OF RESEARCH

Efforts are being made to improve irradiance forecasts by improving the internal model processes (e.g. WRF-Solar [START_REF] Jimenez | WRF-Solar: An augmented NWP model for solar power prediction. Model description and clear sky assessment[END_REF] or optimal model configuration [START_REF] López-Coto | Comparison between several parameterization schemes in WRF for solar forecasting in coastal zones[END_REF]). Nevertheless, such efforts do not necessarily improve the forecast results if the initial conditions of the NWP simulation already lack precision.

Geostationary meteorological satellite data assimilation

The quality of NWP forecasts largely depends upon the initial conditions which are used to launch the NWP model. Data assimilation is the procedure which statistically combines observations and short-range NWP forecasts (first guess) in order to derive the best possible analysis of the atmospheric state on the entire NWP model grid [START_REF] Kalnay | Atmospheric Modeling, Data Assimilation, and Predictability[END_REF]. Various meteorological observations (synoptic stations, buoys, radiosondes, aircraft measurements, satellites, radars, etc.) can potentially be used to derive the initial atmospheric state for NWP. The various types of observations have to be treated differently in the data assimilation process.

Besides the advantages of high temporal and spatial resolution and coverage, geostationary meteorological satellite observations provide pertinent information about cloud presence and evolution. Leveraging the remote sensing technology together with data assimilation techniques is a promising solution for improving solar irradiance forecasts, especially in data-sparse regions like islands.

Different existing approaches aim at deriving atmospheric analyses with most realistic cloud features, utilising satellite observations. Two approaches exist: the direct assimilation of observed radiances and the assimilation of cloud physical properties.

Direct radiance assimilation

This technique is massively performed in global NWP models by weather centers around the world.

A radiative transfer model is used to simulate a radiance observation from the model variables. Within the data assimilation procedure the model's initial state is subsequently modified by statistically combining the synthetic observations and the actual radiance observations measured by the satellite [START_REF] Weng | Advances in Radiative Transfer Modeling in Support of Satellite Data Assimilation[END_REF].

Since the optical and thermal sensors on geostationary meteorological satellites do not always capture information inside clouds, the data assimilation procedure has to adapt to different cloudiness conditions. Therefore, in many current assimilation systems, most cloud-affected radiance observations are discarded [START_REF] Bauer | Assimilating Satellite Observations of Clouds and Precipitation into NWP Models[END_REF].

Since the direct assimilation of radiances brings along a number of critical issues to solve, the method has not yet been exploited in regional models regarding the tropics. Accordingly, an appropriate radiance data assimilation system for GHI forecasts must address the following issues:

 observation and background error estimation  an accurate estimation of atmospheric profiles and cloud microphysical properties  the choice of appropriate satellite channels  optimal thinning of radiance observations  bias correction.

Physical cloud property assimilation

Besides the direct assimilation of radiances, another approach is to assimilate geostationary satellite observations into limited-area models is the assimilation of derived cloud properties. This is an active field of research since approximately the beginning of this century [START_REF] Bayler | NWP Cloud Initialization Using GOES Sounder Data and Improved Modeling of Nonprecipitating Clouds[END_REF] [START_REF] Yucel | Short-Term Performance of MM5 with Cloud-Cover Assimilation from Satellite Observations[END_REF].

Recent research focusses on the assimilation of cloud-top information like cloud mask, cloud top-height, cloud-top temperature, cloud type classification, cloud liquid and ice water path and cloud albedo. In one existing approach the initial water vapour profile of the WRF model has been directly modified using cloud-top properties derived by GOES imagery [START_REF] Mathiesen | A highresolution, cloud-assimilating numerical weather prediction model for solar irradiance forecasting[END_REF]. Another initiative uses MSG imagery, an ensemble Kalman filter and the COSMO NWP model in order to assimilate cloud-top height and relative humidity [START_REF] Schomburg | A concept for the assimilation of satellite cloud information in an Ensemble Kalman Filter: singleobservation experiments[END_REF]). Cloud top property assimilation has also been experimentally and successfully performed in the Rapid Update Cycle of the HIRLAM model, operated by the Royal Netherlands Meteorological Institute [START_REF] De Haan | Cloud Initialization in the Rapid Update Cycle of HIRLAM[END_REF]. All recently published approaches that address cloud property assimilation show improved cloud cover forecasts or solar irradiance forecasts, respectively. This shows the general potential of such techniques to improve cloud cover and solar irradiance forecasts. Nevertheless, special attention to the tropics has not been drawn so far.

Although all approaches from recent years provide improved cloud cover forecasts or irradiance forecasts in mid-latitudes, the major limitation is that only clear-sky or completely cloud cases can be considered. Fractional clouds cause a measured signal mixing cold clouds and warmer Earth surface. As a consequence, clouds are usually considered as single-layer clouds.

AREAS OF STUDY

So far there has not been any research focus on the assimilation of geostationary satellite observations in regional NWP with the goal of improving cloud cover forecasts in the tropics. The aim of this work is to evaluate the benefits of this practice in the noninterconnected French oversea territories Reunion Island and French Guiana (Figure 1). Going along with the global increase of installed PV power, these zones likewise experience a rapid increase, e.g. 8.3 % increase of installed solar power on Reunion Island in 2013 [START_REF]Le bilan énérgétique de la Réunion[END_REF].

While Reunion Island is surrounded by ocean for hundreds of kilometres, French Guiana has both continental and maritime influence. Reunion Island has a distinct orography with up to 3000 m high mountains which facilitate local thermal circulations. There are no major elevations in French Guiana which makes it a good site for studies of coastal clouds which are influenced by sea breeze effects.

OUTLOOK AND GOALS OF THIS WORK

For the first time, we evaluate the potential of assimilating satellite-derived cloud information in tropical regions. In a first study on Reunion Island, Meteosat-7 (IODC service) observations are assimilated in the COSMO model using its Ensemble Kalman filter (EnKF) data assimilation system [START_REF] Schraff | Kilometre-scale ensemble data assimilation for the COSMO model (KENDA)[END_REF]. This serves as an initial evaluation of remotely sensed cloud information assimilation in a state-of-the-art regional NWP modelling system, applied to the tropics. Surface irradiance measurements are used in order to validate the simulation results.

In the next step, a satellite data assimilation strategy for the WRF model is chosen. In the course of this, the results of the experiments using COSMO help approaching diverse issues: the interplay of nesting and satellite data assimilation, satellite channel selection, bias correction and observation error treatment. Applying two models to the same region allows to compare the influences of changes in the data assimilation procedure on the simulation results.

The chosen WRF satellite data assimilation strategy is then being applied and tested on Reunion Island (using Meteosat IODC data) and in French Guiana (using GOES-East data).

CONCLUSION

In conclusion, this work is an innovative research approach, aiming at increasing the injection of solar power into the electricity grids of non-interconnected tropical territories.

The assimilation of radiances and cloud properties into limited-area models is not yet exploited and constitutes a promising approach for improving cloud cover and solar irradiance forecasts in the tropics.

The raised solution could be a major contribution to facilitate the energy autonomy of such territories by an optimal use of the abundant solar energy yield which is largely under-exploited so far.

Reunion Island and French Guiana constitute ideal tropical, yet European experimental sites for testing solutions for solar power injection. Nevertheless, the outcome of this work will be applicable to all noninterconnected territories in tropical areas. It might permit the numerous developing and very populated countries in the tropics to increase their green economic growth by an optimal use of solar energy.

Figure 1 :

 1 Figure 1: Geographic locations of French Guiana and Reunion Island.