Mugurel Ionut
email: mugurel.andreica@cs.pub.ro

Algorithmic Techniques for Several Optimization Problems Regarding Distributed Systems with Tree Topologies

As the development of distributed systems progresses, more and more challenges arise and the need for developing optimized systems and for optimizing existing systems from multiple perspectives becomes more stringent. In this paper I present novel algorithmic techniques for solving several optimization problems regarding distributed systems with tree topologies. I address topics like: reliability improvement, partitioning, coloring, content delivery, optimal matchings, as well as some tree counting aspects. Some of the presented techniques are only of theoretical interest, while others can be used in practical settings.

Introduction

Distributed systems are being increasingly developed and deployed all around the world, because they present efficient solutions to many practical problems. However, as their development progresses, many problems related to scalability, fault tolerance, stability, efficient resource usage and many other topics need to be solved. Developing efficient distributed systems is not an easy task, because many system parameters need to be fine tuned and optimized. Because of this, optimization techniques are required for designing efficient distributed systems or improving the performance of existing, already deployed ones. In this paper I present several novel algorithmic techniques for some optimization problems regarding distributed systems with a tree topology.

Trees are some of the simplest non-trivial topologies which appear in reallife situations. Many of the existing networks have a hierarchical structure (a tree or tree-like graph), with user devices at the edge of the network and router backbones at its core. Some peer-to-peer systems used for content retrieval and indexing have a tree structure. Multicast content is usually delivered using multicast trees. Furthermore, many graph topologies can be reduced to tree topologies, by choosing a spanning tree or by covering the graph's edges with edge disjoint spanning trees [START_REF] Roskind | A Note on Finding Minimum-Cost Edge-Disjoint Spanning Trees[END_REF]. In a tree, there exists a unique path between any two nodes. Thus, the network is quite fragile. The fragility is compensated by the simplicity of the topology, which makes many decisions become easier.

This paper is structured as follows. Section 2 defines the main notations which are used in the rest of the paper. In Section 3 I consider the minimum weight cycle completion problem in trees. In Section 4 I discuss two tree partitioning problems and in Section 5 I consider two content delivery optimization problems. In Section 6 I solve several optimal matching problems in trees and powers of trees and in Section 7 I analyze the first fit online coloring heuristic, applied to trees. In Section 8 I consider three other optimization and tree counting problems. In Section 9 I discuss related work and in Section 10 I conclude and present future work.

Notations

A tree is an undirected, connected, acyclic graph. A tree may be rooted, in which case a special vertex r will be called its root. Even if the tree is unrooted, we may choose to root it at some vertex. In a rooted tree, we define parent(i) as the parent of vertex i and ns(i) as the number of sons of vertex i. For a leaf vertex i, ns(i) = 0 and for the root r, parent(r) is undefined. The sons of a vertex i are denoted by s(i, j) (1 ≤ j ≤ ns(i)). A vertex j is a descendant of vertex i if (parent(j) = i) or parent(j) is also a descendant of vertex i. We denote by T (i) the subtree rooted at vertex i, i.e. the part of the tree composed of vertex i and all of its descendants (together with the edges connecting them). In the paper, the terms node and vertex will be used with the same meaning.

A matching M of a graph G is a set of edges of the graph, such that any two edges in the set have distinct endpoints (vertices). A maximum matching is a matching with maximum cardinality (maximum number of edges).

Minimum Weight Cycle Completion of a Tree

We consider a tree network with n vertices. For m pairs of vertices (i, j) which are not adjacent in the tree, we are given a weight w(i, j) (we can consider w(i, j) = +∞ for the other pairs of vertices). We want to connect some of these m pairs (i.e. add extra edges to the tree), such that, in the end, every vertex of the tree belongs to exactly one cycle. The objective consists of minimizing the total weight of the edges added to the tree. For the unweighted case (w(i, j) = 1) and when we can connect any pair of vertices which is not connected by a tree edge, there exists the following simple greedy algorithm [START_REF] Scortaru | National Olympiad in Informatics[END_REF]. We select an arbitrary root vertex r and then traverse the tree bottom-up (from the leaves towards the root). For each vertex i we will compute a value l(i), representing the largest number of vertices on a path P (i) starting at i and continuing in T (i), such that every vertex j ∈ (T (i) \ P (i)) belongs to exactly one cycle and the vertices in P (i) are the only ones who do not belong to a cycle. We denote by e(i) the second endpoint of the path (the first one being vertex i). For a leaf vertex i, we have l(i) = 1 and e(i) = i. For a non-leaf vertex i, we first remove from its list of sons the sons s(i, j) with l(s(i, j)) = 0, update ns(i) and renumber the other sons starting from 1. If i remains with only one son, we set l(i) = l(s(i, 1)) + 1 and e(i) = e(s(i, 1)). If i remains with ns(i) > 1 sons, we will sort them according to the values l(s(i, j)), such that l(s(i, 1)) ≤ l(s(i, 2)) ≤ . . . ≤ l(s(i, ns(i))). We will connect by an edge the vertices e(s(i, 1)) and e(s(i, 2)). This way, every vertex on the paths P (s(i, 1)) and P (s(i, 2)), plus the vertex i, belong to exactly one cycle. For the other sons s(i, j) (3 ≤ j ≤ ns(i)), we will have to connect s(i, j) to e(s(i, j)). This will only be possible if l(s(i, j)) ≥ 3; otherwise, the tree admits no solution. Afterwards, we set l(i) = 0. If the root r has only one son, then we must have l(r) ≥ 3, such that we can connect r to e(r).

For the general case, I will describe a dynamic programming algorithm (as the greedy algorithm cannot be extended to this case). We will again root the tree at an arbitrary vertex r, thus defining parent-son relationships. For each vertex i, we will compute two values: wA(i)=the minimum total weight of a subset of edges added to the tree such that every vertex in T (i) belongs to exactly one cycle, and wB(i)=the minimum total weight of a subset of edges added to the tree such that every vertex in (T (i) \ {i}) belongs to exactly one cycle (and vertex i belongs to no cycle). We will compute the values from the leaves towards the root. For a leaf vertex i, we have wA(i) = +∞ and wB(i) = 0. For a non-leaf vertex i, we have: wB(i) = ns(i) j=1 wA(s(i, j)). In order to compute wA(i) we will first traverse T (i) and for each vertex j, we will compute wAsum(i, j)=the sum of all the wA(p) values, where p is a son of a vertex q which is located on the path from i to j (P (i . . . j)) and p does not belong to P (i . . . j). We have wAsum(i, i) = wB(i) and for the other vertices j we have wAsum(i, j) = wAsum(i, parent(j)) -wA(j) + wB(j). Now we will try to add an edge, such that it closes a cycle in the tree which contains vertex i. We will first try to add edges of the form (i, j), where j is a descendant of i (but not a son of i, of course) -these will be called type 1 edges. Adding such an edge (i, j) provides a candidate value wcand(i, i, j) for wA(i): wcand(i, i, j) = wAsum(i, j) + w(i, j). We will then consider edges of the form (p, q) (p = i and q = i), where the lowest common ancestor of p and q (LCA(p, q)) is vertex i -these will be called type 2 edges (we consider every pair of distinct sons s(i, a) and s(i, b), and for each such pair we consider every pair of vertices p ∈ T (s(i, a)) and q ∈ T (s(i, b)) and verify if the edge (p, q) can be added to the tree). Adding such an edge (p, q) provides a candidate value wcand(i, p, q) for wA(i): wcand(i, p, q)=wAsum(i, p)+wAsum(i, q)-wB(i)+w(p, q). wA(i) will be equal to the minimum of the candidate values wcand(i, * , * ) (or to +∞ if no candidate value exists). We can implement the algorithm in O(n 2 ) time, which is optimal in a sense, because m ≤ (n

• (n -1)/2 -n + 1), which is O(n 2 ). wA(r)
is the answer to our problem and we can find the actual edges to add to the tree by tracing back the way the wA( * ) and wB( * ) values were computed.

However, when the number m of edges which can be added to the tree is significantly smaller, we can improve the time complexity to O((n + m) • log(n)). We will compute for each of the m edges (i, j) the lowest common ancestor of the vertices i and j (LCA(i, j)) in the rooted tree. This can be achieved by preprocessing the tree in O(n) time and then answering each LCA query in O(1) time [START_REF] Bender | The LCA Problem revisited[END_REF]. If LCA(i, j) = k, then we will add the edge (i, j) to a list Ledge(k). Then, for each non-leaf vertex i, we will traverse the edges in Ledge(k). For each edge (p, q) we can easily determine if it is of type 1 (i = p or i = q) or of type 2 and use the corresponding equation. However, we need the values wAsum(i, p) and wAsum(i, q). Instead of recomputing these values from scratch, we will update them incrementally. It is obvious that wAsum(parent(i), p)=wAsum(i, p)+wB(parent(i))-wA(i). We will preprocess the tree, by assigning to each vertex i its DFS number DF Snum(i) (DF Snum(i)=j if vertex i was the j th distinct vertex visited during a DFS traversal of the tree which started at the root). Then, for each vertex i, we compute DF Smax(i)=the maximum DFS number of a vertex in its subtree. For a leaf node i, we have DF Smax(i) = DF Snum(i). For a non-leaf vertex i, DF Smax(i)=max{DF Snum(i), DF Smax(s(i, 1)), . . . , DF Smax(s(i, ns(i)))}. We will maintain a segment tree, using the algorithmic framework from [START_REF] Andreica | Optimal Offline TCP Sender Buffer Management Strategy[END_REF]. The operations we will use are range addition update and point query. Initially, each leaf i (1 ≤ i ≤ n) has a value v(i) = 0. Before computing wA(i) for a vertex i, we set the value of leaf DF Snum(i) in the segment tree to wB(i). Then, for each son s(i, j), we add the value (wB(i) -wA(s(i, j))

to the interval [DF Snum(s(i, j)), DF Smax(s(i, j))] (range update). We can obtain wAsum(i, p) for any vertex p ∈ T (i) by querying the value of the cell DF Snum(p) in the segment tree: we start from the (current) value of the leaf DF Snum(p) and add the update aggregates uagg stored at every anestor node of the leaf in the segment tree. Queries and updates take O(log(n)) time each.

If the objective is to minimize the largest weight W max of an edge added to the tree, we can binary search W max and perform the following feasibility test on the values W cand chosen by the binary search: we consider only the "extra" edges (i, j) with w(i, j) ≤ W cand and run the algorithm described above for these edges; if wA(r) = +∞, then W cand is feasible.

Tree Partitioning Techniques 4.1 Tree Partitioning with Lower and Upper Size Bounds

Given a tree with n vertices, we want to partition the tree into several parts, such that the number of vertices in each part is at least Q and at most k • Q (k ≥ 1). Each part P must have a representative vertex u, which does not necessarily belong to P . However, (P ∪ {u}) must form a connected subtree. I will present an algorithm which works for k ≥ 3. We root the tree at any vertex r, traverse the tree bottom-up and compute the parts in a greedy manner. For each vertex i we compute w(i)=the size of a connected component C(i) in T (i), such that vertex i ∈ C(i) , |C(i)| < Q, and all the vertices in (T (i) \ C(i)) were split into parts satisfying the specified properties. For a leaf vertex i, w(i) = 1 and C(i) = {i}. For a non-leaf vertex i, we traverse its sons (in any order) and maintain a counter ws(i)=the sum of the w(s(i, j)) values of the sons traversed so far. If ws(i) exceeds Q-1 after considering the son s(i, j), we form a new part from the connected components C(s(i, last son + 1)), . . . , C(s(i, j)) and assign vertex i as its representative. Then, we reset ws(i) to 0. (last son < j) is the previous son where ws(i) was reset to 0 (or 0, if ws(i) was never reset to 0).

After considering every son of vertex i, we set w(i) = ws(i) + 1 and the component C(i) is formed from the components C(s(i, j)) which were not used for forming a new part, plus vertex i. If ws(i) + 1 = Q, then we form a new part from the component C(i) and set w(i) = 0 and C(i) = {}. During the algorithm, the maximum size of any part formed is 2 • Q -2. At the end of the algorithm, we may have that w(r) > 0. In this case, the vertices in C(r) were not assigned to any part. However, at least one vertex from C(r) is adjacent to a vertex assigned to some part P . Then, we can extend that part P in order to contain the vertices in C(r). This way, the maximum size of a part becomes 3 • Q -3. The pseudocode of the first part of the algorithm is presented below. In order to compute the parts, we maintain for each vertex i a value part(i), which is 0, initially (0 means that the vertex was not assigned to any part). In order to assign distinct part numbers, we will maintain a global counter part number, whose initial value is 0. The first part of the algorithm has linear time complexity (O(n)). The second part (adding C(r) to an already existing part) can also be performed in linear time, by searching for an edge (p, q), such that part(p) = 0 and part(q) > 0 (there are only n -1 = O(n) edges in a tree).

LowerUpperBoundTreePartitioning(Q, i): if (ns(i)=0) then w(i)=1 else ws(i)=last son=0 for j=1 to ns(i) do // j=1,2,. . . ,ns(i) LowerUpperBoundTreePartitioning(Q, s(i,j)) ws(i)=ws(i)+w(s(i,j)) if (ws(i) ≥ Q) then part number=part number + 1; last son=j; ws(i)=0
for k=last son+1 to j do AssignPartNumber(s(i,k), part number)

w(i)=ws(i)+1 if (w(i) ≥ Q) then part number=part number + 1; w(i)=0
AssignPartNumber(i, part number)

AssignPartNumber(i, part number): if (part(i) = 0) then return() part(i)=part number for j=1 to ns(i) do AssignPartNumber(s(i,j), part number)

Connected Tree Partitioning

I will now present an efficient algorithm for identifying k connected parts of given sizes in a tree (if possible), subject to minimizing the total cost. Thus, given a tree with n vertices, we want to find k vertex-disjoint components (called parts), such that the

i th part (1 ≤ i ≤ k) has sz(i) vertices (sz(1)+sz(2)+. . .+sz(k) ≤ n and sz(i) ≤ sz(i + 1) for 1 ≤ i ≤ k -1)
. Each tree edge (i, j) has a cost ce(i, j) and each tree vertex i has a cost cv(i). We want to minimize the sum of the costs of the vertices and edges which do not belong to any part. An edge (i, j) belongs to a part p if both vertices i and j belong to part p.

In order to obtain k connected components of the given sizes we need to keep Q -k edges of the tree and remove the others, where Q = sz(1) + . . . + sz(n). We could try all the ((n-1) choose (Q-k)) possibilities of choosing Q-k edges out of the n -1 edges of the tree. For each possibility, we obtain k ′ = n -Q + k connected components with sizes sz ′ (1) ≤ sz ′ (2) ≤ . . . ≤ sz ′ (k ′ ); in case of several components with equal sizes, we sort them in increasing order of the total cost of the vertices in them. Then, we must have sz(j) = sz ′ (k ′ -k + j) and the total cost of the possibility is the sum of the costs of the removed edges plus the sum of the costs of the vertices in the components 1, 2, . . . , k ′ -k (which should have only one vertex each, if the size conditions hold). However, this approach is quite inefficient in most cases. I will present an algorithm with time complexity O(n 3 • 3 k ). We root the tree at an arbitrary vertex r. Then, we compute a table Cmin(i, j, S)=the minimum cost of obtaining from T (i) the parts with indices in the set S and, besides them, we are left with a connected component consisting of j vertices which includes vertex i and, possibly, several vertices which are ignored (if j = 0, then every vertex in T (i) is assigned to one of the parts in S or is ignored). We compute this table bottom-up:

ConnectedTreePartitioning(i):
for each S ⊆ {1, 2, . . . , k} do for j=0 to n do Cmin(i,j,S)=+∞ Cmin(i, 1, {})=0; Cmin(i, 0, {})=cv(i) for x=1 to ns(i) do ConnectedTreePartitioning(s(i,x)) for each S ⊆ {1, 2, . . . , k} do for j=0 to n do Caux(i,j,S)=Cmin(i,j,S); Cmin(i,j,S)=+∞ for each S ⊆ {1, 2, . . . , k} do for j=0 to n do for each W ⊆ S do for q=0 to qlimit(j) do Cmin(i,j,S)=min{Cmin(i,j,S), Caux(i,j-q,S \ W ) + extra cost(i,s(i,x),q) + Cmin(s(i,x),q,W)} for each S ⊆ {1, 2, . . . , k} do for j=0 to n do if (Cmin(i, j, S) < +∞) then for q=1 to k do if ((j=sz(q)) and (q / ∈ S)) then Cmin(i,0,S ∪ {q})=min{Cmin(i,j,S), Cmin(i,0,S ∪ {q})} We define extra cost(i, son x i, q)=if (q > 0) then return(0) else return(ce(i, son x i)) and qlimit(j)=max{j-1,0}. The algorithm computes Cmin(i,*,*) from the values of vertex i's sons, using the principles of tree knapsack. The total amount of computations for each vertex is O(ns(i)•3 k •n 2 ). Summing over all the vertices, we obtain O(n 3 •3 k ). The minimum total cost is Cmin(r, 0, {1, 2, . . . , k}) (if this value is +∞, then we cannot obtain k parts with the given sizes). In order to find the actual parts, we need to trace back the way the Cmin( * , * , * ) values were computed, which is a standard procedure. When the sum of the sizes of the k parts is n, then every vertex belongs to one part.

Content Delivery Optimization Problems

Minimum Number of Unicast Streams

We consider a directed acyclic graph G with n vertices and m edges. Every directed edge (u, v) has a lower bound lbe G (u, v), an upper bound ube G (u, v) and a cost ce G (u, v). Every vertex u has a lower bound lbv G (u), an upper bound ubv G (u) and a cost cv G (u). We need to determine the minimum number of unicast communication streams p and a path for each of the p streams, such that the number of stream paths npe(u, v) containing an edge (u, v) satisfies lbe G (u, v) ≤ npe(u, v) ≤ ube G (u, v) and the number of paths npv(u) containing a vertex u satisfies lbv G (u) ≤ npv(u) ≤ ubv G (u). Each vertex u can be a source node, a destination node, both or none. A stream's path may start at any source node and finish at any destination node. Moreover, for the number of streams p, we want to compute the paths such that the sum S over all the values

(npe(u, v) -lbe G (u, v))•ce G (u, v) and (npv(u) -lbv G (u))•cv G (u) is minimum.
Particular cases of this problem have been studied previously. When lbv G (u) =1 and ubv G (u) = 1 for every vertex u, lbe G (u, v) = 0 and ube G (u, v) = +∞ for every directed edge (u, v), all the costs are 0, and every vertex is a source and destination node, we obtain the minimum path cover problem in directed acyclic graphs, which is solved as follows [START_REF] Cormen | Introduction to Algorithms[END_REF]. Construct a bipartite graph B with n vertices x 1 , . . . , x n on the left side and n vertices y 1 , . . . , y n on the right side. We add an edge (x i , y j ) in B if the directed edge (i, j) appears in G. Then, we compute a maximum matching in B. If the cardinality of this matching is C, then we need p = n -C streams. The paths are computed as follows. Having an edge (x i , y j ) in the maximum matching means that the edge (i, j) in G belongs to some stream's path. If two edges (x i , y j ) and (x j , y k ) in B belong to the matching, then the edges (i, j) and (j, k) in G belong to the path of the same stream. For non-zero costs, we compute a minimum (total) weight matching in B (where every edge (x i , y j ) has a weight equal to ce(i, j)).

In order to solve the problem I mentioned, we will use a standard transformation and construct a new graph G ′ where every vertex u is represented by two vertices u in and u out . For every directed edge (u, v) in G, we add an edge (u out , v in ) in G ′ , with the same cost and lower and upper bounds. We also add a directed edge from u in to u out in G ′ (for every vertex u in G), with cost cv G (u), lower bound lbv G (u) and upper bound ubv G (u). Then we add two special vertices s (source) and t (sink) to G ′ . For every source node u in G, we add a directed edge (s, u in ) in G ′ , with lower bound and cost 0 and upper bound +∞. For every destination node v in G, we add a directed edge (v out , t), with lower bound and cost 0 and upper bound +∞. We also add the edges (s, t) and (t, s) with lower bound and cost 0 and upper bound +∞. The resulting graph G ′ has costs, lower and upper bounds only on its edges and not on its vertices. In order to compute the minimum number of communication streams which satisfy the constraints imposed by G, it is enough to compute a (minimum cost) minimum feasible flow in G ′ , from s to t. Decomposing the flow into unit-flow paths (in order to obtain the path of each communication stream) can then be done easily. We repeatedly perform a graph traversal (DFS or BFS) from s to t in G ′ , considering only directed edges with positive flow on them. From the traversal tree, by following the "parent" pointers, we can find a path P from s to t, containing only edges with positive flow. We compute the minimum flow f P on any edge of P , transform P into f P unit paths and then decrease the flow on the edges in P by f P . If we remove the first and last vertices on any unit path (i.e. s and t), we obtain a path from a vertex u in to a vertex v out , where u is a source node in G and v is a destination node in G. We will use the algorithm presented in [START_REF] Cormen | Introduction to Algorithms[END_REF] for determining a feasible flow (not necessarily minimum) in a flow network with lower and upper bounds on its edges. We will denote this algorithm by A(F, s, t) (F is the flow network given as argument, s is the source vertex and t is the sink vertex). I will describe A(F, s, t) briefly. We construct a new graph F ′ from F , as follows. We maintain all the vertices and edges in F . For every directed edge (u, v) in F , the directed edge (u, v) in F ′ has the same cost, lower bound 0 and upper bound (ube F (u, v) -lbe F (u, v)). We add two extra vertices s ′ and t ′ and the following zero-cost directed edges: (s ′ , u) and (u, t ′ ) for every vertex u in F (including s and t). The lower bound of every edge will be 0. The upper bound of a directed edge (s ′ , u) in F ′ is equal to the sum of the lower bounds of the directed edges ( * , u) in F . The upper bound of every directed edge (u, t ′ ) in F ′ is equal to the sum of the lower bounds of the directed edges (u, * ) in F . The algorithm A(F, s, t) computes a minimum cost maximum flow g in the graph F ′ (which, as stated, only has upper bounds); if all the costs are 0, only a maximum flow is computed. If g is equal to the sum of the upper bounds of the edges (s ′ , * ) (or, equivalently, of the edges ( * , t ′ )), then a feasible flow from s to t exists in F : the flow on every directed edge (u, v) in F will be lbe F (u, v) plus the flow on the edge (u, v) in F ′ .

We will first run the algorithm on G ′ (i.e. call A(G ′ , s, t)) in order to verify if a feasible flow exists). If no feasible flow exists, then the constraints cannot be satisfied by any number of streams. Otherwise, we construct a graph G ′′ from G ′ , by adding a new vertex snew and a zero-cost directed edge (snew, s) with lower bound 0 and upper bound x. snew will be the new source vertex and x is a parameter which is used in order to limit the amount of flow entering the old source vertex s. We will now perform a binary search on x, between 0 and gmax, where gmax is the value of the feasible flow computed by calling A(G ′ , s, t). The feasibility test consists of verifying if there exists a feasible flow in the graph G ′′ (i.e. calling A(G ′′ , snew, t)). The minimum value of x for which a feasible flow exists in G ′′ is the value of the minimum feasible flow in G ′ , from s to t. Obtaining the feasible flow in G ′ from the feasible flow in G ′′ is trivial: for every directed edge (u, v) in G ′ , we set its amount of flow to the flow of the same edge (u, v) in G ′′ . The time complexity of the presented algorithm is O(M F (n, m) • log(gmax)), where gmax is a good upper bound on the value of a feasible flow and M F (n, m) is the best time complexity of a (minimum cost) maximum flow algorithm in a directed graph with n vertices and m edges.

Degree-Constrained Minimum Spanning Tree

In [START_REF] Magnanti | Optimal Trees[END_REF], the following problem was considered: given an undirected graph with n verices and m edges, where each edge (i, j) has a weight w(i, j) > 0, compute a spanning tree M ST of minimum total weight, such that a special vertex r has degree exactly k in M ST . A solution was proposed, based on using a parameter d and setting the cost of each edge (r, j) adjacent to r, c(r, j) = d + w(r, j); the cost of the other edges is equal to their weight. Parameter d can range from -∞ to +∞. We denote by M ST (d)=the minimum spanning tree using the cost functions defined previously. When d = -∞, M ST (d) contains the maximum number of edges adjacent to r. For d = +∞, M ST (d) contains the minimum number of edges adjacent to r. We define the function ne(d)=the number of edges adjacent to r in M ST (d). ne(d) is non-increasing on the interval [-∞, +∞]. We will binary search the smallest value dopt of the parameter d in the interval [-∞, +∞], such that ne(dopt) ≤ k. We will finish the binary search when the length of the search interval is smaller than a small constant ε > 0.

If ne(dopt) = k, then the edges in M ST (dopt) form the required minimum spanning tree. If ne(dopt) < k, then ne(dopt -ε) > k. We define S(d)=the set of edges adjacent to vertex r in M ST (d). It is easy to prove that S(dopt) is included in S(dopt -ε). The required minimum spanning tree is constructed in the following manner. The edges adjacent to vertex r will be the edges in S(dopt), to which we add (k -ne(dopt)) arbitrary edges from the set S(doptε) \ S(dopt). Once these edges are fixed, we construct the following graph G: we set the cost of the chosen edges to 0 and the cost of the other edges (i, j) to w(i, j). We now compute a minimum spanning tree M ST G in G. The edges in M ST G are the edges of the minimum spanning tree of the original graph, in which vertex r has degree exactly k. The time complexity of this approach is O(m • log(m) • log(DM AX)), where DM AX denotes the range over which we search the parameter d. When m is not too large (i.e. m is not of the order O(n 2 )), this represents an improvement over the O(n 2 ) solution given in [START_REF] Magnanti | Optimal Trees[END_REF].

6 Matching Problems

Maximum Weight Matching in an Extended Tree

Let's consider a rooted tree (with vertex r as the root). Each vertex i has a weight w(i). We want to find a matching in the following graph G (extended tree), having the same vertices as T and an edge (x, y) between two vertices x and y, if: (i) x and y are adjacent in the tree; (ii) x and y have the same parent in the tree. The weight of an edge (x, y) in G is |w(x) -w(y)|. The weight of a matching is the sum of the weights of its edges. We are interested in a maximum weight matching in the graph G. For each vertex i, we sort its sons s(i, 1), . . . , s(i, ns(i)) in non-decreasing order of their weights, i.e. w(s(i, 1)) ≤ . . . ≤ w(i, ns(i)). We will compute for each vertex i two values: A(i)=the maximum weight of a matching in T (i) if vertex i is the endpoint of an edge in the matching and B(i)=the maximum weight of a matching in T (i) if vertex i is not the endpoint of any edge in the matching. In order to compute these values, we will compute the following tables for every vertex i: CA(i, j, k)=the maximum weight of a matching in T (i) if vertex i is the endpoint of an edge in the matching and we only consider its sons s(i, j), s(i, j+1), . . . , s(i, k) (and their subtrees). Similarly, we have CB(i, j, k), where vertex i does not belong to any edge in the matching. The maximum weight of a matching is max{A(r), B(r)}. The actual matching can be computed easily, by tracing back the way the A(i), B(i), CA(i, * , * ) and CB(i, * , * ) values were computed. A recursive algorithm (called with r as its argument) is given below. The time complexity is O(ns(i) 2 ) for a vertex i and, thus, O(n 2 ) overall.

MaximumWeightMatching-ExtendedTree(i): if (ns(i)=0) then A(i)=B(i)=0 else for j=1 to ns(i) do MaximumWeightMatching-ExtendedTree(s(i,j)) for j=1 to ns(i) do CA(i, j, j -1)= -∞; CA(i, j, j) = |w(i) -w(s(i, j))| + B(s(i, j)) CB(i, j, j -1)= 0; CB(i, j, j)=max{A(s(i,j)), B(s(i,j))} for count=1 to (ns(i)-1) do for j=1 to (ns(i)-count) do k = j + count CA(i,j,k)=max{|w(s(i, j)) -w(s(i, k))| + B(s(i,j)) + B(s(i,k)) + CA(i, j + 1, k -1), |w(i) -w(s(i, j))| + B(s(i,j)) + CB(i, j+1, k), |w(i) -w(s(i, k))| + B(s(i,k)) + CB(i, j, k-1), max{A(s(i,j)), B(s(i,j))} + CA(i, j+1, k), max{ A(s(i,k)), B(s(i,k))} + CA(i, j, k-1)} CB(i,j,k)=max{|w(s(i, j)) -w(s(i, k))| + B(s(i,j)) + B(s(i,k)) + CB(i, j + 1, k -1), max{A(s(i,j)), B(s(i,j))} + CB(i, j+1, k), max{A(s(i,k)), B(s(i,k))} + CB(i, j, k-1)} A(i)=CA(i,1,ns(i)); B(i)=CB(i,1,ns(i))

Maximum Matching in the Power of a Graph

The k th power G k (k ≥ 2) of a graph G is a graph with the same set of vertices as G, where there exists an edge (x, y) between two vertices x and y if the distance between x and y in G is at most k. The distance between two vertices (x, y) in a graph is the minimum number of edges which need to be traversed in order to reach vertex y, starting from vertex x. A maximum matching in G k of a graph G can be found by restricting our attention to a spanning tree T of G. The following linear algorithm (called with i = r), using observations from [START_REF] Qinglin | Factors and Factor Extensions[END_REF], solves the problem (we consider that, initially, no vertex is matched):

MaximumMatchingGk(i): if (ns(i)=0
) then return() else last son=0 for j=1 to ns(i) do // j=1,2,. . . ,ns(i) MaximumMatchingGk(s(i,j))

if (not matched(s(i,j)) then if (last son = 0) then last son = s(i,j) else add edge (last son, s(i,j)) to the matching matched(last son) = matched(s(i,j)) = true; last son = 0 if (last son > 0) then add edge (i, last son) to the matching matched(i) = matched(last son) = true

First Fit Online Tree Coloring

A very intuitive algorithm for coloring a graph with n vertices is the first-fit online coloring heuristic. We traverse the vertices in some order v(1), v(2), . . . , v(n). We assign color 1 to v(1) and for i = 2, . . . , n, we assign to v(i) the minimum color c(i) ≥ 1 which was not assigned to any of its neighbours v(j) (j < i).

A tree is 2-colorable: we root the tree at any vertex r and then compute for each vertex i its level in the tree (distance from the root); we assign the color 1 to the vertices on even levels and the color 2 to those on odd levels. However, in some situations, we might be forced to process the vertices in a given order. In this case, it would be useful to compute the worst-case coloring that can be obtained by this heuristic, i.e. the largest number of colors that are used, under the worst-case ordering of the tree vertices (Grundy number). I will present an O(n • log(log(n))) algorithm for this problem, similar in nature to the linear algorithm presented in [START_REF] Hedetniemi | A Linear Algorithm for the Grundy (Coloring) Number of a Tree[END_REF]. For each vertex i, we will compute cmax(i)=the largest color the can be assigned to vertex i in the worst-case, if vertex i is the last vertex to be colored. The value max{cmax(i)|1 ≤ i ≤ n} is the largest number of colors that can be assigned by the first fit online coloring heuristic.

We will root the tree at an arbitrary vertex r. The algorithm consists of two stages. In the first stage, the tree is traversed bottom-up and for each vertex i we compute c(1, i)=the largest color that can be assigned to vertex i, considering only the tree T (i). For a leaf vertex i, we have c(1, i) = 1. For a non-leaf vertex i, we will sort its sons s(i, 1), . . . , s(i, ns(i)), such that c(1, s(i, 1)) ≤ c(1, s(i, 2)) ≤ . . . ≤ c(1, s(i, ns(i))). We will initialize c(1, i) to 1 and then consider the sons in the sorted order. When we reach son s(i, j), we compare c(1, s(i, j)) with c(1, i). If c(1, s(i, j)) ≥ c(1, i), then we increment c(1, i) by 1 (otherwise, c(1, i) stays the same). The justification of this algorithm is the following: if a vertex i can be assigned color c(1, i) in some ordering of the vertices in T (i), then there exists an ordering in which it can be assigned any other color c ′ , such that 1 ≤ c ′ ≤ c(1, i). Then, when traversing the sons and reaching a son s(i, j) with c(1, s(i, j)) ≥ c(1, i), we consider an ordering of the vertices in T (s(i, j)), where the color of vertex s(i, j) is c(1, i); thus, we can increase the maximum color that can be assigned to vertex i.

After the bottom-up tree traversal, we have cmax(r) = c(1, r), but we still have to compute the values cmax(i) for the other vertices of the tree. We could do that by rooting the tree at every vertex i and running the previously described algorithm, but this would take O(n 2 • log(log(n))) time. However, we can compute these values faster, by traversing the tree vertices in a top-down manner (considering the tree rooted at r). For each vertex i, we will compute colmax(parent(i), i)=the maximum color that can be assigned to parent(i) if we remove T (i) from the tree and afterwards we consider parent(i) to be the (new) root of the tree. We will use the values c(2, i) as temporary storage variables. c(2, i) is initialized to c(1, i), for every vertex i. When computing cmax(i), we consider that vertex i is the root of the tree. Let's assume that we computed the value cmax(i) of a vertex i and now we want to compute the value cmax(j) of a vertex j which is a son of vertex i. We remove j from the list of sons of vertex i and add parent(i) to this list (parent(i)=vertex i's parent in the tree rooted at the initial vertex r). We now need to lift vertex j above vertex i and make j the new root of the tree. In order to do this, we will recompute the value c(2, i), which is computed similarly to c(1, i), except that we consider the new list of sons for vertex i (and their c(2, * ) values). Afterwards, we add vertex i to the list of sons of vertex j. We will compute the value cmax(j) similarly to the value c(1, j), using the values c(2, * ) of vertex j's sons (instead of the c(1, * ) values of the sons). After computing cmax(j) we restore the lists of sons of vertices i and j to their original states (as if the tree were rooted at the initial vertex r). After computing the values cmax(u) of all the descendants u of a vertex j, we reset the value c(2, j) to c(1, j).

Both traversals take O(n • log(n)) time, if we sort the ns(i) sons of every vertex i in O(ns(i) • log(ns(i))) time. However, it has been proved in [START_REF] Hedetniemi | A Linear Algorithm for the Grundy (Coloring) Number of a Tree[END_REF] that the minimum number of vertices of a tree with the Grundy number q is 2 q-1 , which is the binomial tree B(q -1). The binomial tree B(0) consists of only one vertex. The binomial tree B(k ≥ 1) has a root vertex with k neighbors; the i th of these neighbors (0 ≤ i ≤ k -1) is the root of a B(i) binomial tree. Thus, every value c(1, * ), c(2, * ) and cmax( * ) can be represented using O(log(log(n))) bits. We can use radix-sort and obtain an O(n • log(log(n))) time complexity. The pseudocode of the functions is given below. The main algorithm consists of calling FirstFit-BottomUp(r), initializing the c(2, * ) values to the c(1, * ) values, setting cmax(r) = c(1, r) and then calling FirstFit-TopDown(r) Compute(i, idx): sort the sons of vertex i, such that c(idx,s(i,1))≤ . . . ≤c(idx,s(i,ns(i))) c(idx,i)=1 for j=1 to ns(i) do if (c(idx,s(i,j))≥c(idx,i)) then c(idx,i)=c(idx,i)+1

FirstFit-BottomUp(i): for j=1 to ns(i) do FirstFit-BottomUp(s(i,j)) Compute(i, 1) FirstFit-TopDown(i): if (i = r) then remove vertex i from the list of sons of parent(i) add parent(parent(i)) to the list of sons of parent(i) (if parent(i) = r) Compute(parent(i),2); colmax(parent(i),i)=c(2,parent(i)) add parent(i) to the list of sons of vertex i Compute(i,2); cmax(i)=c(2,i) restore the original lists of sons of the vertices parent(i) and i for j=1 to ns(i) do FirstFit-TopDown(s(i,j)) c(2,i)=c (1,i) 8 Other Optimization and Counting Problems In this subsection I consider the following optimization problem: We are given a sequence of n leaves and each leaf i (1 ≤ i ≤ n) has a height h(i). We want to construct a (strict) binary tree with n -1 internal nodes, such that, in an inorder traversal of the tree, we encounter the n leaves in the given order. The height of an internal node i is h(i) = 1+max{h(lef tson(i), h(rightson(i))} (the height of the leaves is given). We are interested in computing a tree whose root has minimum height. A straight-forward dynamic programming solution is the following: compute Hmin(i, j)=the minimum height of a tree containing the leaves i, i + 1, . . . , j. We have: Hmin(i, j)=1 + min i≤k≤j-1 max{Hmin(i, k), Hmin(k + 1, j)}. Hmin (1, n) is the answer to our problem. However, the time complexity of this algorithm is O(n 3 ), which is unsatisfactory. An optimal, linear-time algorithm was given in [START_REF] Mu | On Building Trees with Minimum Height, Relationally[END_REF]. The main idea of this algorithm is the following. We traverse the leaves from left to right and maintain information about the rightmost path of the optimal tree for the first i leaves. Then, we can add the (i + 1) st leaf by modifying the rightmost path of the optimal tree for the first i leaves. Let's assume that we processed the first i leaves and the optimal tree for these leaves contains, on its rightmost path, the vertices v(1), v(2), . . . , v(nv(i)), in order, from the root to the rightmost leaf (v(1) is the root). Let's assume that the heights of the subtrees rooted at these vertices are hv(1), . . . , hv(nv(i)). It is easy to build this tree for i = 1 and i = 2 (it is unique). When adding the (i + 1) st leaf, we traverse the rightmost path from nv(i) down to 2. Assume that we are considering the vertex v(j). If hv(j -1) < (2 + max{hv(j), h(i + 1)}), then we disconsider the vertex v(j) from the rightmost path and move to the next vertex (v(j -1)). Let's assume that the path now contains the vertices v(1), . . . , v(nv ′ (i)). We replace vertex v(nv ′ (i)) by a new vertex vnew, whose left son will be v(nv ′ (i)) (together with its subtree) and whose right son will be the (i + 1) st leaf. The height of the new vertex will be 1 + max{hv(nv ′ (i)), h(i + 1)}. The rightmost path of the optimal tree behaves like a stack and, thus, the overall time complexity is linear.

I will present a sub-optimal O(n • log(n)) time algorithm which is interesting on its own. The algorithm is similar to Huffman's algorithm for computing optimal prefix-free codes, except that it maintains the order of the leaves. A suggestion that such an approach might work was given to me by C. Gheorghe. At step i (1 ≤ i ≤ n-1) of the algorithm, we will have n-i+1 subtrees of the optimal tree. Each subtree j contains an interval of leaves [lef tleaf (j), rightleaf (j)] and its height is h(j). We will combine the two adjacent subtrees j and j + 1 whose combined height (1+max{height(subtree j),height(subtree j+1)}) is minimum among all the O(n) pairs of adjacent subtrees. At the first step, the n subtrees are represented by the n leaves, whose heights are given. A straightforward implementation of this idea leads to an O(n 2 ) algorithm. However, the processing time can be improved by using two segment trees [START_REF] Andreica | Optimal Offline TCP Sender Buffer Management Strategy[END_REF], A and B, with n and n-1 leaves, respectively. Each node q of a segment tree corresponds to an interval of leaves [lef t(q), right(q)] (leaves are numbered starting from 1). Each leaf node of the segment tree A can be in the active or inactive state. Each node q of A (whether leaf or internal node) maintains a value nactive(q), denoting the number of active leaves in its subtree. Initially, each of the n leaves of A is active and the nactive( * ) values are initialized appropriately, in a bottom-up manner (1, for a leaf node, and nactive(lef tson(q)) + nactive(rightson(q)), for an internal node q). Segment tree B has n -1 leaves and each node of B (leaf or internal node) stores a value hc. If leaf i (1 ≤ i ≤ n -1) is active in A, then hc(leaf i)=1+max{h(i), h(j)}, where j > i is the next active leaf. If leaf i is not active in A or is the last active leaf, then hc(leaf i)=+∞. The value hc of each internal node q of B is the minimum among all the hc values of the leaves in node q's subtree, i.e. hc(node q)=min{hc(leftson(q)), hc(rightson(q))}. Moreover, each node q of B maintains the number lnum of the leaf in its subtree which gives the value hc(node q). We have lnum(leaf i)=i and lnum(internal node q)=if (hc(leftson(q)) ≤ hc(rightson(q))) then lnum(leftson(q)) else lnum(rightson(q)).

At each step i (1 ≤ i ≤ n -1), each active leaf is the leftmost leaf of a subtree of the optimal tree. After every step, the number of active leaves decreases by 1. We can find in O(log(n)) time the pair of adjacent subtrees to combine. The height of the combination of these subtrees is hc(root node of B), the leftmost leaf of the first subtree is i=lnum(root node of B) and that of the second subtree is j = next active(i). We define the function next active by using two other functions: rank(i) and unrank(r). rank(i) returns the number of active leaves before leaf i (0 ≤ rank(i) ≤ nactive(root node of A)-1). unrank(r) returns the index of the leaf whose rank is r. The two functions are inverses of each other: unrank(rank(i)) = i and rank(unrank(r)) = r. We have rank(i)=rank'(i, root node of A), unrank(r)=unrank'(r, root node of A) and next active(i)=unrank(rank(i) + 1). rank'(i, q): if (q is a leaf node) then if (left(q)=right(q)=i) then return(0) else return(-1) else if (i > right(lef tson(q))) then return(nactive(leftson(q))+rank'(i, rightson(q))) else return(rank'(i, leftson(q))) unrank'(r, q): if (q is a leaf node) then if (r > 0) then return(-1) else return(left(q)) else if (nactive(lef tson(q)) ≤ r) then return(unrank'(r-nactive(leftson(q)), rightson(q))) else return(unrank'(r, leftson(q)))

The functions rank, unrank and next active take O(log(n)) time each. After obtaining the indices of the two active leaves i and j whose corresponding subtrees are united (by adding a new internal node whose left son is the root of i's subtree and whose right son is the root of j's subtree), we mark leaf j as inactive. We do this by traversing the segment tree A from leaf j towards the root (from j to parent(j), parent(parent(j)), . . . , root node of A) and decrement by 1 the nactive values of the visited nodes. Then, we change the h values of leaves i and j. We set h(i)=hc(root node of B) and h(j) = +∞. After this, we will also change the hc values associated to the leaves i and j in the segment tree B. The new hc value of leaf j will be +∞. If i is now the last active leaf, then hc(leaf i) becomes +∞, too. Otherwise, let j ′ = next active(i), the next active leaf after i (at this point, leaf j is not active anymore). We will change hc(leaf node i) to (1 + max{h(i), h(j ′ )}). After changing the hc value of a leaf k, we traverse the tree from leaf k towards the root (visiting all of k's ancestors, in order, starting from parent(k) and ending at the root of B). For each ancestor node q, we recompute hc(node q) as min{hc(lef tson(q)), hc(rightson(q))}.

The Number of Trees with a Fixed Number of Leaves

In order to compute the number of labeled trees with n vertices and exactly p leaves, we will compute a table N T (i, j)=the number of trees with i vertices and exactly j leaves (1 ≤ j ≤ i ≤ n). Obviously, we have N T (1, 1) = N T (2, 2) = 1 and N T (i, j) = 0 for i = 1, 2 and j = i. For i > 2, we have N T (i, i) = 0 and for 1 ≤ j ≤ i -1, we will proceed as follows. The j leaves can be chosen in C(i, j) ways (i choose j). After choosing the identifiers of the j leaves, we will conceptually remove the leaves from the tree, thus remaining with a tree having i -j vertices and any number of leaves k (1 ≤ k ≤ j). Each of the j leaves that we conceptually removed is adjacent to one of these k vertices. Furthermore, each of these k vertices is adjacent to at least one of the j leaves from the larger tree. Thus, we need to compute the number of surjective functions f from a domain of size j to a domain of size k. We will denote this value by N F (j, k). This is a "classical" problem, but I will present a simple solution, nevertheless. We have N F (0, 0) = 1 and N F (j, k) = 0, if j < k. In order to compute the values for k ≥ 1 and j ≥ k, we will consider every number g of values x from the set {1, . . . , j} for which f (x) = k. Once g is fixed, we have C(j, g) ways of choosing the g values from the set {1, . . . , j}. For each such possibility we have N F (j -g, k -1) ways of extending it to a surjective function. Thus, N F (j, k) = j g=1 C(j, g) • N F (j -g, k -1). We can tabulate all the N F ( * , * ) values in O(n 3 ) time (after tabulating the combinations C( * , * ) in O(n 2 ) time, first). With the N F ( * , * ) values computed, we have N T (i, j) = C(i, j)• j k=1 (N T (i-j, k)•N F (j, k)). We can easily compute each entry N T (i, j) in O(n) time, obtaining an O(n 3 ) overall time complexity. The technique of performing dynamic programming on successive layers of leaves of a tree is also useful in several other counting problems.

The Number of Trees with Degree Constraints

We want to compute the number of unlabeled, rooted trees with n ≥ 2 vertices, such that the (degree / number of sons) of each vertex belongs to a set S, which is a subset of {0, 1, 2, . . . , n -1}. By (a/b) we mean that a refers to the degreeconstrained problem and b refers to the number-of-sons-constrained problem (everything else being the same). Because every tree with n ≥ 2 vertices must contain at least a leaf (a vertex of degree 1) and at least one vertex with at least 1 son, the set S will always contain the subset ({1}/{0, 1}). We will compute a table N T (i, j, p)=the number of trees with i vertices, such that the root has degree j (j sons) and the maximum number of vertices in the subtree of any son of the root is p; moreover, except perhaps the tree root, the (degrees/numbers of sons) of all the other vertices belong to the set S. Because the trees are unlabeled, we can sort the sons of each vertex in non-decreasing order of the numbers of vertices in their subtrees. Thus, we will compute the table N T in increasing order of p. N T (1, 0, p) = 1 and N T (1, j > 0, p) = N T (i ≥ 2, j, 0) = 0. For p ≥ 1 and i ≥ 2, we have: Grids), which have many practical applications: minimum weight cycle completion (reliability improvement), constrained partitioning (distributed coordination and control), minimum number of streams and degree-constrained minimum spanning trees (efficient content delivery), optimal matchings (data replication and resource allocation), coloring (resource management and frequency allocation) and tree counting aspects. All these problems are variations or extensions of problems which have been previously posed in other research papers. The presented techniques are either better (faster or more general) than the previous solutions or easier to implement.

8. 1

 1 Building a (Constrained) Tree with Minimum Height

N T (i, j, p)=N T (i, j, p -1)+

• CR(T T (p), k) T T (p) is the total number of trees with p vertices, for which the (degree / number of sons) of the root is equal to some ((x -1)/(x)), x ∈ S, and the (degrees / numbers of sons) of the other vertices belong to the set S. By CR(i, j) we denote combinations with repetitions of i elements, out of which we choose j. Because the argument i can be very large, we cannot tabulate CR(i, j). Instead, we will compute it on the fly. We know that CR(i, j) = C(i + j -1, j) and that C(i, j) = i-j+1 j • C(i, j -1). Thus, CR(i, j) can be computed in O(j) time. Before computing any value N T ( * , * , p), we need to compute and store the values T T (p), T T (p) = x∈S N T (p, ((x-1)/(x)), p-1), and CR(T T (p), k),

The memory storage can be reduced from O(n 3 ) to O(n 2 ), by noticing that the values N T ( * , * , p) are computed based only on the values N T ( * , * , p -1). Thus, we can maintain these values only for the most recent two values of p.

A less efficient method is to compute the numbers T ok(i)=the number of trees with i vertices, such that each vertex satisfies the (degree/number of sons) constraints. T ok(1) = T ok(2) = 1. We will make use of the T T (i) values defined previously, except that they will be computed differently. For every i ≥ 2, we consider every possible number x of sons of the tree root and compute N T 2(i, x)=the number of trees with i vertices, such that the tree root has x sons and all the other vertices satisfy the (degree/number of sons) constraints. We will generate all the possibilities (y(1), y(2), ..., y(i -1)), with 0 ≤ y(j) ≤ i-1 j (1 ≤ j ≤ i-1) and y(1) + . . . + y(i-1)=x. y(j) is the number of sons of the tree root which have j vertices in their subtrees. The number of trees "matching" such a partition is equal to i-1 j=1 CR(T T (j), y(j)). N T 2(i, x) is computed by summing the numbers of trees "matching" every partition. Afterwards, if x ∈ S, we add N T 2(i, x) to T ok(i). If x = ((y -1)/(y)) and y ∈ S, then we add N T 2(i, x) to T T (i). N T 2(i, x) may be added to both T ok(i) and T T (i).

Related Work

Reliability analysis and improvement techniques for distributed systems were considered in [START_REF] Chen | Reliability Analysis of Distributed Systems Based on a Fast Reliability Algorithm[END_REF][START_REF] Kumar | Performance and reliability optimization for distributed computing systems[END_REF]. Reliability analysis and optimization for tree networks in particular were considered in [START_REF] Scortaru | National Olympiad in Informatics[END_REF][START_REF] Andreica | Reliability Analysis of Tree Networks Applied to Balanced Content Replication[END_REF][START_REF] Abachi | Reliability analysis of tree, torus and hypercube message passing architectures[END_REF]. Different kinds of tree partitioning algorithms, based on optimizing several objectives, were proposed in [START_REF] Frederickson | Optimal algorithms for tree partitioning[END_REF][START_REF] Cordone | A subexponential algorithm for the coloured tree partition problem[END_REF][START_REF] Wu | On the Uniform Edge-Partition of a Tree[END_REF]. Problems related to tree coloring were studied in [START_REF] Hedetniemi | A Linear Algorithm for the Grundy (Coloring) Number of a Tree[END_REF]. Content delivery in distributed systems is a subject of high practical and theoretical interest and is studied from multiple perspectives. Communication scheduling in tree networks was considered in many papers (e.g. [START_REF] Henzinger | Scheduling multicasts on unit-capacity trees and meshes[END_REF]) and the optimization of content delivery trees (multicast trees) was studied in [START_REF] Cui | Maxmin overlay multicast: rate allocation and tree construction[END_REF].

Conclusions and Future Work

In this paper I considered several optimization problems regarding distributed systems with tree topologies (e.g. peer-to-peer networks, wireless networks,


Pierre Massip

Philippe Blanc

Lucien Wald
email: lucien.wald@mines-paristech.fr

Data fusion: taking into account the modulation transfer function in ARSISbased pansharpening methods

Keywords: Image fusion, pan-sharpening, multispectral, panchromatic, PSF, MTF, ARSIS, remote sensing

Multispectral images provided by satellite have a poor spatial resolution while panchromatic images (PAN) exhibit a spatial resolution two or four times better. Data fusion is a mean to synthesize MS images at higher spatial resolution than original by exploiting the high spatial resolution of the PAN. This process is often called pan-sharpening. The synthesized multispectral images should be as close as possible to those that would have been acquired by the corresponding sensors if they had this high resolution. The methods based on the concept "Amélioration de la Résolution Spatiale par Injection de Structures" (ARSIS) concept are able to deliver synthesized images with good spectral quality but whose geometrical quality can still be enhanced. We propose to consider the characteristics of the sensor to improve the geometrical quality. We take explicitly into account the modulation transfer function (MTF) of the sensor in the fusion process. Though this study is limited in methods and data, we observe a better restitution of the geometry and an improvement in the majority of quality indices classically used in pan-sharpening. The communication also presents a means to assess the respect of the synthesis property from a MTF point of view.

Introduction

Image fusion is a vivid topic in remote sensing, among other domains. Fusion was defined by the EARSeL Working Group [START_REF] Wald | Some terms of reference in data fusion[END_REF] as "a formal framework in which are expressed means and tools for the alliance of data originating from different sources. It aims at obtaining information of a greater quality, although the exact definition of 'greater quality' will depend on the application". Here we specifically focus on the synthesis of multispectral (MS) images at a higher spatial resolution than original by exploiting the high spatial resolution of another image.

Current satellites provide images of two kinds with opposite characteristics. Multispectral (MS) images have a good spectral quality but a poor spatial resolution whereas panchromatic (PAN) images have on the contrary a high spatial resolution but with a poorer spectral quality. These two types of images allow to well identify observed structures through different information by geometrical information on the [Tapez un texte] one hand and by the spectral signature on the other hand. The interest of combining both types of images has been demonstrated by many authors.

Data fusion is a means to synthesize MS images at higher spatial resolution than original by exploiting the high spatial resolution of the PAN. This process is often called pan-sharpening. Here, the aim of fusion is to perform a high-quality transformation of the MS content when increasing the spatial resolution from the original one to that of the PAN image. The problem may be seen as the inference of the information that is missing to the original MS images for the construction of the MS images synthesized at a better resolution [START_REF] Ranchin | Fusion of high spatial and spectral resolution images: The ARSIS concept and its implementation[END_REF][START_REF] Thomas | Assessment of the quality of fused products[END_REF][START_REF] Thomas | Synthesis of multispectral images to high spatial resolution: A critical review of fusion methods based on remote sensing physics[END_REF]. [START_REF] Ranchin | Image fusion-The ARSIS concept and some successful implementation schemes[END_REF][START_REF] Wald | Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images[END_REF][START_REF] Li | Spatial quality evaluation of fusion of different resolution images[END_REF] establish that the fused images should obey several properties, that were synthesized by the EARSeL Working Group [START_REF] Thomas | Assessment of the quality of fused products[END_REF][START_REF] Thomas | Synthesis of multispectral images to high spatial resolution: A critical review of fusion methods based on remote sensing physics[END_REF]. The consistency property expresses that the fused images should offer a strong consistency with the original data set: if a fused product is downsampled to its original low resolution, the original MS image should be retrieved. The synthesis property deals with the quality of the high-quality transformation: any synthetic image should be as identical as possible to the image that the corresponding sensor would observe with the highest spatial resolution, if existent. This holds for each individual modality MS and the multi-modality set.

Recent published works tend to demonstrate that the methods based on the ARSIS concept better synthesize images with respect to the expected spectral properties [START_REF] Ranchin | Image fusion-The ARSIS concept and some successful implementation schemes[END_REF][START_REF] Thomas | Synthesis of multispectral images to high spatial resolution: A critical review of fusion methods based on remote sensing physics[END_REF][START_REF] Yocky | Multiresolution wavelet decomposition image merger of Landsat Thematic Mapper and SPOT panchromatic data[END_REF][START_REF] Aiazzi | Context-driven fusion of high spatial and spectral resolution images based on oversampled multiresolution analysis[END_REF][START_REF] Wang | A comparative analysis of image fusion methods[END_REF][START_REF] Terretaz | Comparison of different methods to merge SPOT P and XS data : Evaluation in an urban area[END_REF][START_REF] Alparone | Comparison of pansharpening algorithms: Outcome of the 2006 GRS-S data fusion contest[END_REF]. However, the need to refine the spatial/geometrical quality of synthesized image is often underlined [START_REF] Ranchin | Image fusion-The ARSIS concept and some successful implementation schemes[END_REF] [7] [START_REF] Alparone | Comparison of pansharpening algorithms: Outcome of the 2006 GRS-S data fusion contest[END_REF][START_REF] Thomas | Fusion d'images de résolutions spatiales différentes (Fusion of images of different spatial resolutions[END_REF]. Geometry can be defined as the spatial organization of the radiometry in the image [START_REF] Peyré | Géométrie multi-échelles pour les images et les textures[END_REF]. Characteristics of the imaging instrument are important with respect to the restitution of the geometry of the landscape. Besides its spectral and spatial resolutions, an important characteristic of optical instruments is the impulse response which is named point spread function (PSF). The PSF has an associate value, the modulation transfer function (MTF). The MTF is the modulus of the Fourier transform of the PSF. According to [START_REF] Grois | Theoretical approach to CMOS APS PSF and MTF modeling -evaluation[END_REF], the MTF determines the upper limits to the image quality, i.e., the image resolution or sharpness and describes the image quality in terms of contrast as a function of spatial frequency. The difference in sampling rates between MS and PAN images implies a difference in the MTF of these images. This difference must be taken into account when synthesizing MS images at the PAN spatial resolution. A few fusion methods already include this difference in MTF through the multiscale model (MSM) in the ARSIS concept [START_REF] Ranchin | Image fusion-The ARSIS concept and some successful implementation schemes[END_REF] [13] [START_REF] Aiazzi | A MTF-based spectral distortion minimizing model for pan-sharpening of very high resolution multispectral images of urban areas[END_REF][START_REF] Garzelli | Interband structure modeling for Pan-sharpening of very high-resolution multispectral images[END_REF][START_REF] Aiazzi | MTF-tailored multiscale fusion of highresolution MS and Pan imagery[END_REF].

Our manner to treat the MTF is different from these previous works. We take explicitly into account the MTF of the sensor in the fusion process for each modality, independently from the MSM.

Objective

We now express the synthesis property in terms of geometry: the restitution of the radiometric organization of the observed landscape by the synthesized image should be as identical as possible to the organization that the corresponding sensor would observe with the highest spatial resolution, if existent. The work presented here aims at including the MTF into an existing fusion method in order to demonstrate that taking into account the difference in MTF leads to a better respect of the synthesis property from the geometrical point of view without degrading the quality of the other aspects of the synthesized images.

[Tapez un texte]

As said before, we selected an existing method belonging to the ARSIS concept. This method is called ATWT-M3, where ATWT denotes the MSM and means "à trous" wavelet transform, and M3 is the inter-band structure model described by [START_REF] Ranchin | Fusion of high spatial and spectral resolution images: The ARSIS concept and its implementation[END_REF], [START_REF] Wang | A comparative analysis of image fusion methods[END_REF]. This method has been shown to provide good results in most cases; it already includes the difference in MTF in its MSM. Therefore, it is considered as an efficient method against which we will compare our own results. Whereas the results will not be the same than those obtained here, our way of taking into account the difference in MTF may easily apply to other methods belonging to the ARSIS concept and likely less easier, to other methods.

By explicitly integrating the difference in MTF, the aim is to give a 'real' MTF to the synthesized images. In figure 1, are drawn the MTFs for the four MS channels of the Pleiades sensor, discussed later, and that for the PAN, as a function of the frequency f normalized with the sampling frequency of the PAN. In this case, the spatial resolution of the PAN image (70 cm) is 4 times better than the MS (280 cm). Note that in this figure, the MTF of the MS are zero for relative frequencies equal to, or greater than, 0.25 (=1/4). It means that one cannot distinguish details less than 280 cm in size in the MS images. On the contrary, the MTF of PAN is equal to 0.4. The difference in MTF between low and high spatial resolution images is evidenced in this figure: the MTF of MS has a sharper decrease with frequency than that of PAN. Of particular interest here is the relative Nyquist frequency for MS, equal to 0.125. It denotes the smallest frequency that can be detected in MS images, i.e., twice the pixel size. When increasing the spatial resolution by injecting high-frequencies into the original MS images, the MTF of the synthesized image should look like that of PAN. However, without modification of the MTF during the fusion process, the resulting MTF exhibits a discontinuity for the relative frequency 0.125 as schematically drawn in Fig. 2. To attenuate the gap between the value of MS and PAN at this frequency, a solution consist in "raising" the MTF of the MS frequencies in the range [0, 0.125] so that it is close to that of a high resolution image. Doing so provides a MTF closer to a "real" MTF and similar contrast in the image for a same frequency without considering its origin: PAN or MS.

Data set

The data set used for this work is a set of images resulting from the simulation of the Pleiades mission of the Centre National d'Etudes Spatiales (CNES). These images are simulated from images at a higher resolution from the airborne sensor Pelican. CNES kindly provides us with MS images at both high (70 cm) and low (280 cm) spatial resolution and a PAN image at high resolution (70 cm). The tests have been done on an excerpt of the city of Toulouse in France.

These images result from simulation; it means that we have the so-called original MS images at 280 cm, the PAN image at 70 cm and the actual MS images at 70 cm which constitute the reference that should be attained. This is not the common situation in which the actual MS images at high resolution are unknown. However, it is a very convenient situation to assess the influence of a parameter such as the MTF because all parameters are known. During the creation of the simulated Pleiades images, the MTF is controlled by CNES for each image and resolution. Consequently, the MTF of MS, PAN and actual MS images are known for this study.

[Tapez un texte]

Modification of MTF for fusion in ARSIS based methods

As said before, we exploit an existing method: ATWT-M3, and we create a new one: ATWT-M3-FTMadapted by taking the MTF explicitly into account in the method. We transform the MTF before the decomposition of the MS image by the undecimated wavelet transform ("à trous" wavelet transform). The application of the undecimated transform implies a resampling of the original MS images to the high resolution. The modification of the MTF is done on this resampled MS image.

The transformation starts with a deconvolution of the original MS images at low resolution. The deconvolution kernel corresponds to the detector part of the MTF -the optical MTF is not taken into account.-Then, we perform the resampling of this image that has now only the optical MTF. The resampling is of spline type. Then we convolute the resampled images with the MTF of a high resolution detector.

After that, the method ATWT-M3 is applied which is based on undecimated wavelet transform [START_REF] Ranchin | Fusion of high spatial and spectral resolution images: The ARSIS concept and its implementation[END_REF]. For each modality, the high frequencies wavelet coefficients from the PAN image are transformed into MS high frequencies wavelet coefficients by applying an affine function. The parameters of this function are determined by a leastsquare fitting between the PAN and MS low-frequencies wavelet coefficients [START_REF] Ranchin | Fusion of high spatial and spectral resolution images: The ARSIS concept and its implementation[END_REF]. Finally, the inverse wavelet transform is applied to obtain the synthesized MS images at high resolution.

Results

The results of the two methods: ATWT-M3 and ATWT_M3_FTMadapted are compared to the reference image. As recommended by several authors [START_REF] Wald | Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images[END_REF][START_REF] Li | Spatial quality evaluation of fusion of different resolution images[END_REF] [6] [START_REF] Alparone | Comparison of pansharpening algorithms: Outcome of the 2006 GRS-S data fusion contest[END_REF], we perform statistical comparisons and a visual analysis of the possible discrepancies.

Statistical results

Quality budgets are obtained by the very same tool than that used for the "data fusion contest 2006" [START_REF] Alparone | Comparison of pansharpening algorithms: Outcome of the 2006 GRS-S data fusion contest[END_REF]. Table 1 exhibits the result of mono-modal indices: relative bias (biasRel), relative difference variance in percent (diffVarRel), relative standard deviation (σRel), correlation coefficient (cc), and correlation coefficient of high frequency (ccHF), as well multimodal indices: ERGAS and mean SAM.

The ideal values are presented for each index. This table shows that the bias for each method is very close to the ideal value; this is characteristic of the methods based on the ARSIS concept [START_REF] Ranchin | Fusion of high spatial and spectral resolution images: The ARSIS concept and its implementation[END_REF] [7] [START_REF] Alparone | Comparison of pansharpening algorithms: Outcome of the 2006 GRS-S data fusion contest[END_REF]. The bias is smaller for ATWT-M3 than for ATWT-M3-FTMadapted but in both cases, it remains very close to 0. For the other parameters, one notes a large increase in quality budget when adapting the MTF. All these indices demonstrate improvements. Nevertheless, they do not inform on the aspect of the MTF of the synthesized images. Consequently, we perform a comparison of MTF between the reference and the synthesized images in order to evaluate if the modification of the MTF meets the expectations. This comparison is based on the cross correlation in the Fourier domain (of the PSD of the images) and is made relative to the MTF of the reference image: [Tapez un texte]

The MTF rel denotes the behaviour of the MTF of the synthesized image compared to that of the reference image, as a function of the relative frequency (Fig. 3). The ideal result is a straight line for a value of one. Fig. 3 displays the results for band B2. corresponding to the wavelength of the red.

In this figure, the MTF rel with adaptation of MTF (ATWT-M3-FTMadapted) is closer from the ideal value 1 than that of the method ATWT-M3, and this for all frequencies. Consequently, the improvement brought by the adaptation of the MTF is demonstrated in this case. One may note in this figure that the MTF rel exhibits a trough between 0.1 and 0.15. It is believed that this is due to the fact that only the low frequency of the resampled image are injected in the synthesized image and also because information between the PAN and each MS band are not the same [START_REF] Wang | A comparative analysis of image fusion methods[END_REF] which make that the PAN information is not fully exploited to create the high frequencies wavelet coefficients. This is under investigation and should be confirmed.

Visual analysis

Conclusion

We have shown how the difference in MTF between low and high resolution can be taken into account into a "classical" fusion method based on the ARSIS concept. In a given case, we have observed a better restitution of the geometry and an improvement in the majority of quality indices classically used in pan-sharpening. A new measure has been proposed to evaluate the respect of the synthesis property from a geometrical point of view.

[Tapez un texte] Take into account the difference of MTF between MS and PAN images allow to have synthesized images of better quality. This work can be apply to others concept. It could be interesting to observe the difference of behavior of this technique depending on the concept it's applied on. The MTF relative that is given by our tool let us observe improvement and weakness of the MTF of the fuse image. Though this study is limited in methods and data, the present results are encouraging and may constitute a new way to improve the restitution of geometrical features by already efficient fusion methods.

Figure 1 .Figure 2 .

 12 Figure 1. Modulation transfer function for MS and PAN images.

Figure 3 .

 3 Figure 3. MTF of ATW-M3 and ATWT-3-FTM adapted relative to the MTF reference.

Figure 4 Figure 4 .

 44 Figure 4 exhibits three images: one without MTF adaptation (a), one with MTF adaptation (b) and the image of reference (c).In this excerpt we observe that the contour line of the white building is sharper in (b) than in (a) and that the image (b) is the closest to (c). This holds also for color composites and here again, we find an improvement when adapting the MTF.

Table 1 .

 1 Comparison of quality budget for ATWT-M3 and ATWT-M3-FTMadapted

	modalities Indexes	Ideal	ATWT-M3	ATWT-M3-FTMadapted	Delta in percent
		biaisRel	0	-0,001	0,006	-670
		diffVarRel	0	14,4	3,3	77
	bleu	σRel	0	8,7	7,2	18
		Cc	1	0,9833	0,9875	0.4
		ccHF	1	0,9203	0,9224	0.2