Mathieu Pont
email: mathieu.pont@sorbonne-universite.fr.

Jules Vidal
email: jules.vidal@sorbonne-universite.fr.

Julie Delon
email: julie.delon@u-paris.fr

Julien Tierny
email: julien.tierny@sorbonne-universite.fr.

• M Pont

Wasserstein Distances, Geodesics and Barycenters of Merge Trees

come   L'archive ouverte pluridisciplinaire

Wasserstein Distances, Geodesics and Barycenters of Merge Trees Mathieu Pont, Jules Vidal, Julie Delon and Julien Tierny Fig. 1. The merge trees of three members (a-c) of the Isabel ensemble (wind velocity) concisely and visually encode the number and salience of the features of interest found in the data (eyewall and region of high speed wind, blue and cyan). They also describe how these features are globally connected in the data. In these trees, branches with a low persistence (less than 20% of the function range) are shown with small white arcs. The pointwise mean for the three members (d) exhibits 5 salient maxima (due to distinct eyewall locations, blue, cyan and black) and its merge tree is not representative of the input trees (containing at most 3 large features). In contrast, the Wasserstein barycenter (e) is representative of the input trees, with a number and persistence of large branches that better match the input trees (a-c). Our framework for distances, geodesics and barycenters enables a variety of merge tree based applications, including (f) feature tracking, (g) temporal reduction -key frames are automatically identified (white insets) and deleted merge trees (blue insets) are accurately reconstructed with geodesics -and (h) ensemble clustering and summarization -the clusters and centroids automatically computed by our approach provide a visual summary of the main trends of features found in the ensemble.

Abstract-This paper presents a unified computational framework for the estimation of distances, geodesics and barycenters of merge trees. We extend recent work on the edit distance [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] and introduce a new metric, called the Wasserstein distance between merge trees, which is purposely designed to enable efficient computations of geodesics and barycenters. Specifically, our new distance is strictly equivalent to the L 2 -Wasserstein distance between extremum persistence diagrams, but it is restricted to a smaller solution space, namely, the space of rooted partial isomorphisms between branch decomposition trees. This enables a simple extension of existing optimization frameworks [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF] for geodesics and barycenters from persistence diagrams to merge trees. We introduce a task-based algorithm which can be generically applied to distance, geodesic, barycenter or cluster computation. The task-based nature of our approach enables further accelerations with shared-memory parallelism. Extensive experiments on public ensembles and SciVis contest benchmarks demonstrate the efficiency of our approach -with barycenter computations in the orders of minutes for the largest examples -as well as its qualitative ability to generate representative barycenter merge trees, visually summarizing the features of interest found in the ensemble. We show the utility of our contributions with dedicated visualization applications: feature tracking, temporal reduction and ensemble clustering. We provide a lightweight C++ implementation that can be used to reproduce our results.

Index Terms-Topological data analysis, merge trees, scalar data, ensemble data

INTRODUCTION

Modern datasets, acquired or simulated, are continuously gaining in geometrical complexity, thanks to the ever-increasing accuracy of acquisition devices or computing power of high performance systems. This geometrical complexity makes interactive exploration and analysis difficult, which challenges the interpretation of the data by the end users. This motivates the definition of expressive data abstractions, capable of capturing the main features of the data into concise representations, which visually convey the most important information to the users.

In that context, Topological Data Analysis (TDA) [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] forms a family of generic, robust, and efficient techniques whose utility has been demonstrated in a number of visualization tasks [START_REF] Heine | A survey of topology-based methods in visualization[END_REF] for revealing the implicit structural patterns present in complex datasets. Examples of popular application fields include turbulent combustion [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topologybased data segmentation[END_REF][START_REF] Gyulassy | Stability of dissipation elements: A case study in combustion[END_REF][START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF], material sciences [START_REF] Favelier | Visualizing ensembles of viscous fingers[END_REF][START_REF] Gyulassy | Topologically clean distance fields[END_REF][START_REF] Gyulassy | Interstitial and interlayer ion diffusion geometry extraction in graphitic nanosphere battery materials[END_REF], nuclear energy [START_REF] Maljovec | Topology-inspired partition-based sensitivity analysis and visualization of nuclear simulations[END_REF], fluid dynamics [START_REF] Kasten | Two-dimensional timedependent vortex regions based on the acceleration magnitude[END_REF], bioimaging [START_REF] Anderson | Topological data analysis of functional MRI connectivity in time and space domains[END_REF][START_REF] Bock | TopoAngler: Interactive Topology-Based Extraction of Fishes[END_REF][START_REF] Carr | Simplifying Flexible Isosurfaces Using Local Geometric Measures[END_REF], quantum chemistry [START_REF] Bhatia | Topoms: Comprehensive topological exploration for molecular and condensed-matter systems[END_REF][START_REF] Guenther | Characterizing molecular interactions in chemical systems[END_REF][START_REF] Olejniczak | A Topological Data Analysis Perspective on Non-Covalent Interactions in Relativistic Calculations[END_REF] or astrophysics [START_REF] Shivashankar | Felix: A topology based framework for visual exploration of cosmic filaments[END_REF][START_REF] Sousbie | The persistent cosmic web and its filamentary structure: Theory and implementations[END_REF]. Among the data abstractions developed in TDA (see Sec. 1.1), the merge tree [START_REF] Carr | Computing contour trees in all dimensions[END_REF], which describes the global structure of the connected components of the sub-level sets of scalar datasets (Fig. 2), is a prominent example in the visualization literature [START_REF] Bock | TopoAngler: Interactive Topology-Based Extraction of Fishes[END_REF][START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topologybased data segmentation[END_REF][START_REF] Carr | Simplifying Flexible Isosurfaces Using Local Geometric Measures[END_REF].

In practice, in addition to the increasing geometrical complexity of datasets, users are also confronted to the emergence of ensemble datasets, where a given phenomenon is not described with only one dataset, but with a collection of datasets, called ensemble members. Regarding topological features, a topological data abstraction such as the merge tree can be computed for each ensemble member (possibly insitu [START_REF] Ayachit | Paraview catalyst: Enabling in situ data analysis and visualization[END_REF][START_REF] Bauer | In-situ methods, infrastructures, and applications on high performance computing platforms[END_REF]). Then, a major challenge for end users is the interpretation of the resulting ensemble of merge trees. To address this, a statistical analysis framework for merge trees is needed, requiring several key building blocks, such as: distances (to compare merge trees), geodesics (to visualize optimum transitions between them), and barycenters (to visualize one merge tree representative of a set). These building blocks have been well studied for persistence diagrams [START_REF] Lacombe | Large Scale computation of Means and Clusters for Persistence Diagrams using Optimal Transport[END_REF][START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF][START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF]. However, persistence diagrams suffer from a lack of specificity (Fig. 3), which can prevent the identification of distinct feature trends within the ensemble.

This paper addresses this problem by introducing a unified computational framework for the automatic computation of distances, geodesics, barycenters and clusters of merge trees. In particular, we extend recent work on the edit distance [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] and introduce a new metric, called the Wasserstein distance between merge trees, which is purposely designed to enable efficient computations of geodesics (i.e. length minimizing morphings) and barycenters. In that regard, our work can be interpreted as an extension of previous work on the edit distance [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF], to adapt it to the optimization strategy previously developed for the computation of barycenters of persistence diagrams [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF]. We present efficient, taskbased algorithms using shared-memory parallelism, resulting in the computation of distances, geodesics and barycenters in practical times for real-life datasets. We illustrate the utility of each of our contributions in dedicated visualization tasks. First, we show that our distance computation algorithm can be used for a merge-tree based tracking of features through time. Second, we show that our framework for computing geodesics between merge trees can be used for the reliable sub-sampling of temporal sequences of merge trees. Third, we illustrate the utility of our barycenters for clustering ensemble members based on their merge trees, while providing cluster centroids which visually summarize the main features of interest present in each cluster.

Related work

The literature related to our work can be classified into three main groups, reviewed in the following: (i) uncertainty visualization, (ii) ensemble visualization, and (iii) topological methods for ensembles. (i) Uncertainty visualization: Variability in data can be modeled and encoded in several ways. In particular, uncertain datasets capture variability by modeling each point of the domain as a random variable, whose variability is explicitly modeled by an estimator of an a priori probability density function (PDF). The analysis of uncertain data is a notoriously challenging problem in visualization, described in several surveys [START_REF]2008 uncertainty of measurement -part 3: Guide to the expression of uncertainty in measurement (GUM)[END_REF][START_REF] Bonneau | Overview and state-of-the-art of uncertainty visualization[END_REF][START_REF] Johnson | A next step: Visualizing errors and uncertainty[END_REF][START_REF] Maceachren | Visualizing geospatial information uncertainty: What we know and what we need to know[END_REF][START_REF] Pang | Approaches to uncertainty visualization[END_REF][START_REF] Potter | From quantification to visualization: A taxonomy of uncertainty visualization approaches[END_REF]. Early techniques focused on estimating the entropy of the random variables [START_REF] Potter | Visualization of uncertainty without a mean[END_REF], their correlations [START_REF] Pfaffelmoser | Visualization of global correlation structures in uncertain 2D scalar fields[END_REF] or their gradient variations [START_REF] Pfaffelmoser | Visualizing the variability of gradients in uncertain 2D scalar fields[END_REF]. The positional uncertainty of level sets has been studied for several interpolation schemes and PDF models [5-7, 87, 90-93, 100]. Similarly, the positional uncertainty of critical points has been studied for Gaussian [START_REF] Liebmann | Critical points of gaussian-distributed scalar fields on simplicial grids[END_REF][START_REF] Otto | Uncertain 2D vector field topology[END_REF][START_REF] Otto | Uncertain topology of 3D vector fields[END_REF][START_REF] Petz | Probabilistic local features in uncertain vector fields with spatial correlation[END_REF] or uniform distributions [START_REF] Bhatia | Flow visualization with quantified spatial and temporal errors using edge maps[END_REF][START_REF] Günther | Mandatory critical points of 2D uncertain scalar fields[END_REF][START_REF] Szymczak | Hierarchy of stable Morse decompositions[END_REF]. A general limitation of existing methods for uncertain data is their dependence on the specific PDF model for which they have been designed. This reduces their usability for ensemble data, where the PDF estimated from the ensemble members can follow an arbitrary, unknown model. Also, most existing techniques for uncertain data do not consider multi-modal PDF models, which is however necessary when several, distinct trends are present in the ensemble data. (ii) Ensemble visualization: Another way to model and encode variability in data consists in considering ensemble datasets. In this setting, the variability is directly encoded by an ensemble of empirical observations (i.e. the members of the ensemble). Current approaches to ensemble visualization typically compute some geometrical objects describing the features of interest (level sets, streamlines, etc), for each member of the ensemble. Then, an aggregation phase estimates a representative object for the resulting ensemble of geometrical objects. For Fig. 3. The persistence diagram, D( f i ), and the merge tree, T ( f i ), both visually summarize the number, data range and salience of the features of interest present in the data. However, the persistence diagram represents each individual feature independently, while the merge tree additionally describes how they connect together. This results in a lack of specificity for the persistence diagram which can yield identical data representations for significantly distinct datasets (from left to right, the gaussians with white and cyan spheres have been swapped). In contrast, the merge tree captures this nuance and produces two distinct data representations.

instance, spaghetti plots [START_REF] Diggle | The Analysis of Longitudinal Data[END_REF] are a typical example for studying level-set variability, especially for weather data [START_REF] Potter | Ensemble-vis: A framework for the statistical visualization of ensemble data[END_REF][START_REF] Sanyal | Noodles: A tool for visualization of numerical weather model ensemble uncertainty[END_REF]. More specifically, boxplots [START_REF] Whitaker | Contour boxplots: A method for characterizing uncertainty in feature sets from simulation ensembles[END_REF] describe the variability of contours and curves [START_REF] Mirzargar | Curve boxplot: Generalization of boxplot for ensembles of curves[END_REF]. For flow ensembles, Hummel et al. [START_REF] Hummel | Comparative visual analysis of lagrangian transport in CFD ensembles[END_REF] introduce a Lagrangian framework for classification purposes. Clustering techniques have been investigated, to identify the main trends, and their variability, in ensembles of streamlines [START_REF] Ferstl | Streamline variability plots for characterizing the uncertainty in vector field ensembles[END_REF] and isocontours [START_REF] Ferstl | Visual analysis of spatial variability and global correlations in ensembles of iso-contours[END_REF]. However, only few approaches have applied this overall aggregation strategy to topological objects. Favelier et al. [START_REF] Favelier | Persistence Atlas for Critical Point Variability in Ensembles[END_REF] and Athawale et al. [START_REF] Athawale | Uncertainty visualization of 2d morse complex ensembles using statistical summary maps[END_REF] introduced approaches for analyzing the variability of critical points and gradient separatrices respectively. Several techniques attempted to generate an aggregated contour tree from an ensemble based on overlap-driven heuristics [START_REF] Kraus | Visualization of uncertain contour trees[END_REF][START_REF] Wu | A contour tree based visualization for exploring data with uncertainty[END_REF]. Recently, Lohfink et al. [START_REF] Lohfink | Fuzzy contour trees: Alignment and joint layout of multiple contour trees[END_REF] introduced an approach for the consistent layout of multiple contour trees, to support effective visual comparisons between the contour trees of the distinct members of an ensemble. Although the above techniques addressed the visualization of ensembles of topological objects, they did not focus explicitly on the computation of a representative of multiple topological objects, such as barycenters. (iii) Topological methods: Concepts and algorithms from computational topology [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] have been investigated, adapted and extended by the visualization community for more than twenty years [START_REF] Heine | A survey of topology-based methods in visualization[END_REF][START_REF] Yan | Scalar field comparison with topological descriptors: Properties and applications for scientific visualization[END_REF]. Popular topological representations include the persistence diagram [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF][START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] (Sec. 2.2), which represents the population of features of interest in function of their salience, and which can be computed via matrix reduction [START_REF] Bauer | Distributed computation of persistent homology[END_REF][START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF]. The Reeb graph [START_REF] Biasotti | Reeb graphs for shape analysis and applications[END_REF], which describes the connectivity evolution of level sets, has also been widely studied and several efficient algorithms have been documented [START_REF] Doraiswamy | Computing reeb graphs as a union of contour trees[END_REF][START_REF] Parsa | A deterministic o(m log m) time algorithm for the reeb graph[END_REF][START_REF] Pascucci | Robust on-line computation of Reeb graphs: simplicity and speed[END_REF][START_REF] Tierny | Loop surgery for volumetric meshes: Reeb graphs reduced to contour trees[END_REF], including parallel algorithms [START_REF] Gueunet | Task-based Augmented Reeb Graphs with Dynamic ST-Trees[END_REF]. Efficient algorithms have also been documented for its variants, the merge and contour trees [START_REF] Carr | Computing contour trees in all dimensions[END_REF][START_REF] Tarasov | Construction of contour trees in 3d in o(n log n) steps[END_REF] (Sec. 2.3), and parallel algorithms have also been described [START_REF] Acharya | A parallel and memory efficient algorithm for constructing the contour tree[END_REF][START_REF] Carr | Parallel peak pruning for scalable SMP contour tree computation[END_REF][START_REF] Gueunet | Task-Based Augmented Contour Trees with Fibonacci Heaps[END_REF][START_REF] Maadasamy | A hybrid parallel algorithm for computing and tracking level set topology[END_REF]. The Morse-Smale complex [START_REF] Bremer | A Multi-Resolution Data Structure for 2-Dimensional Morse Functions[END_REF][START_REF] Edelsbrunner | Morse-smale complexes for piecewise linear 3-manifolds[END_REF][START_REF] Edelsbrunner | Hierarchical morse complexes for piecewise linear 2-manifolds[END_REF], which depicts the global behaviour of integral lines, is another popular topological data abstraction in visualization [START_REF] Floriani | Morse complexes for shape segmentation and homological analysis: discrete models and algorithms[END_REF]. Robust and efficient algorithms have been introduced for its computation [START_REF] Gyulassy | Shared-Memory Parallel Computation of Morse-Smale Complexes with Improved Accuracy[END_REF][START_REF] Robins | Theory and Algorithms for Constructing Discrete Morse Complexes from Grayscale Digital Images[END_REF][START_REF] Shivashankar | Parallel Computation of 3D Morse-Smale Complexes[END_REF] based on Discrete Morse Theory [START_REF] Forman | A User's Guide to Discrete Morse Theory[END_REF].

Distance metrics, which are necessary ingredients for the computation of barycenters, have been studied for most of the above objects. Inspired by the literature in optimal transport [START_REF] Kantorovich | On the translocation of masses[END_REF][START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF], the Wasserstein distance between persistence diagrams [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] (Sec. 2.2) and its variant the Bottleneck distance [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] have been extensively studied. They are based on a bipartite assignment problem, for which exact [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] and approximate [START_REF] Bertsekas | A new algorithm for the assignment problem[END_REF][START_REF] Kerber | Geometry helps to compare persistence diagrams[END_REF] implementations are publicly available [START_REF] Tierny | The Topology ToolKit[END_REF]. Several similarity measures have been introduced for Reeb graphs [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3d shapes[END_REF] and their variants [START_REF] Saikia | Extended branch decomposition graphs: Structural comparison of scalar data[END_REF]. However, since these measures are not distance metrics (the preservation of the triangle inequality is not specifically enforced), they do not seem conducive to barycenter computation. Stable distance metrics between Reeb graphs [START_REF] Bauer | Measuring distance between Reeb graphs[END_REF] and merge trees [START_REF] Morozov | Interleaving distance between merge trees[END_REF] have been studied from a theoretical point of view but their computation, following an exponential time complexity, is not tractable for practical datasets in general, except if reliable correspondence labels between the nodes of the trees are provided on the input [START_REF] Gasparovic | Intrinsic interleaving distance for merge trees[END_REF], which is not practical either for large ensembles. Distances with polynomial time computation algorithms have also been investigated. Similarly to our overall strategy, Beketayev et al. [START_REF] Beketayev | Measuring the distance between merge trees[END_REF] focus on a dual representation, Fig. 4. In this example with two datasets f i and f j (a), the optimal matchings (gray) with regard to the edit distance D E [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] (b) map a maximum to a saddle (red spheres). The resulting linear interpolation T α ( f i → f j ) (b) does not describe a shortest path between the input trees (c, top):

D E T ( f i ), T α ( f i → f j ) + D E T α ( f i → f j ), T ( f j ) > D E T ( f i ), T ( f j ) .
In contrast, our new metric W T 2 enables linear interpolations (d) which exactly coincide with shortest paths (c, bottom). The numbers included in (c) are the actual values for D E (top) and W T 2 (bottom) for this example.

the branch decomposition tree (BDT, Sec. 2.3), but in contrast to our approach, they estimate their distances by iteratively reducing a target mismatch term, in particular, over a search space significantly larger than ours. Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] specialize efficient algorithms for computing constrained edit distances between trees [121] to the special case of merge trees (see Appendix 1), resulting in a distance which is computable for real-life datasets and with acceptable practical stability. However, it is not conducive to simple barycenter computations. Indeed, the linear interpolation of the optimal node assignments induced by this metric (Fig. 4) does not result in a shortest path, and hence generates inaccurate midpoints (i.e. inaccurate barycenters given two trees). This further implies that there is no clear or simple strategy for the general computation of barycenters according to that metric.

Regarding the estimation of a representative object from a set of topological representations, several approaches emerged recently. A recent line of work [START_REF] Gasparovic | Intrinsic interleaving distance for merge trees[END_REF][START_REF] Yan | A structural[END_REF] introduced a framework for computing a 1-center of a set of merge trees (i.e. minimizing its maximum distance to the set), according to an interleaving distance. However, as documented by its authors, this approach requires pre-existing, reliable correspondence labels between the nodes of all the input trees, which is not practical with real-life datasets (heuristics need to be considered). Also, since they minimize their maximum distance to a set, 1-centers are typically sensitive to outliers, which prevents their usage for estimating trends or supporting clustering tasks (which typically focus on densities rather than maximum distances). This is further evaluated in Sec. 7.2. In contrast, our approach focuses on the estimation of barycenters (instead of 1-centers) and computes a tree which minimizes its average distance to an ensemble of merge trees (instead of its maximum distance), which is less sensitive to outliers, which better captures trends and which supports clustering tasks. Moreover, the node correspondences between the barycenter and the input trees are automatically estimated by our approach via an assignment optimization present at the core of our distance estimation. Thus our method does not require input correspondences, which makes it readily applicable to real-life ensembles. Several methods [START_REF] Lacombe | Large Scale computation of Means and Clusters for Persistence Diagrams using Optimal Transport[END_REF][START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF][START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] have been introduced for the automatic estimation of barycenters (Sec. 2.2) of persistence diagrams (or vectorized variants [START_REF] Adams | Persistence Images: A Stable Vector Representation of Persistent Homology[END_REF][START_REF] Bubenik | Statistical topological data analysis using persistence landscapes[END_REF]). However, the persistence diagram can lack specificity in its data characterization (Fig. 3). This limitation is addressed by our work which focuses instead on merge trees.

Contributions

This paper makes the following new contributions:

1. A practical distance metric between merge trees: We extend recent work on the edit distance [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] and introduce a new distance between merge trees, which, in contrast to previous work, is purposely designed to enable efficient computations of geodesics and barycenters. It can be computed efficiently, it has acceptable practical stability and it has a strong connection to established metrics, which eases its interpretation. Specifically, it can be interpreted as a variant of the L 2 -Wasserstein distance for persistence diagrams, for which we constrain the underlying search space to account for the additional structural information provided by the merge tree.

2. A simple approach for computing geodesics between merge trees: Given our new metric, we present a simple approach for computing geodesics between merge trees. It uses a simple linear interpolation of the assignments resulting from our new metric, enabling the exact computation of geodesics in linear time. This follows from previous work on persistence diagram geodesics [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF] and it is made possible thanks to a new, local normalization procedure, guaranteeing the topological consistency of the interpolated trees. 3. An approach for computing barycenters between merge trees: Our method for geodesics between merge trees enables a straightforward adaptation of previous optimization strategies for persistence diagram barycenters [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF], resulting, to our knowledge, in the first approach for the computation of barycenters of merge trees. 4. Unified computational framework: We present a unified computational framework for the estimation of distances, geodesics, barycenters, and clusters of merge trees. In particular, we introduce an efficient, task-based algorithm adapted from previous work on edit distances [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF]121], which is generically applicable to any of the above tasks. Our algorithm supports shared-memory parallelism, allowing for further accelerations in practice. 5. Applications: We illustrate the utility of each of our contributions with dedicated visualization tasks, including feature tracking, temporal reduction and ensemble clustering and summarization. 6. Implementation: We provide a lightweight C++ implementation of our algorithms that can be used for reproduction purposes.

PRELIMINARIES

This section presents the theoretical background of our work. It contains definitions adapted from the Topology ToolKit [START_REF] Tierny | The Topology ToolKit[END_REF]. We refer the reader to textbooks [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] for an introduction to computational topology.

Input data

The input data is an ensemble of N piecewise linear (PL) scalar fields f i : M → R, with i ∈ {1, . . . , N}, defined on a PL d-manifold M , with d ≤ 3 in our applications. The sub-level set of f i , noted f i -1

-∞ (w) = {p ∈ M | f i (p) < w}, is defined as the pre-image of (-∞, w) by f i . The super-level set of f i is defined symmetrically:

f i -1 +∞ (w) = {p ∈ M | f i (p) > w}.
As w continuously increases, the topology of f i -1 -∞ (w) changes at specific vertices of M , called the critical points of f i [START_REF] Banchoff | Critical points and curvature for embedded polyhedral surfaces[END_REF]. In practice, f i is enforced to contain only isolated, non-degenerate critical points [START_REF] Edelsbrunner | Hierarchical Morse-Smale complexes for piecewise linear 2-manfiolds[END_REF][START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF]. Critical points are classified by their index I i : 0 for minima, 1 for 1-saddles, d -1 for (d -1)-saddles and d for maxima.

Persistence diagrams

The persistence diagram is a visual summary of the topological features (i.e. connected components, independent cycles, voids) of f i -1 -∞ (w). Specifically, each topological feature of f i -1 -∞ (w) can be associated with a unique pair of critical points (c, c ), corresponding to its birth and death. The Elder rule [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF] states that critical points can be arranged according to this observation in a set of pairs, such that each critical point appears in only one pair (c, c ), with f i (c) < f i (c ) and

I i (c) = I i (c ) -1. For instance, if two connected components of f i -1 -∞ (w)
meet at a critical point c , the younger component (created last, in c) dies, in favor of the older one (created first). Then the persistence diagram, noted D( f i ), embeds each pair to a single point in 2D at coordinates f i (c), f i (c ) . The persistence of a pair is given by its height f i (c )f i (c). Then, the persistence diagram provides a visual overview of the features of interest of a dataset (Fig. 2), where salient features stand out from the diagonal while pairs corresponding to noise are located in the vicinity of the diagonal. Note that, in addition to its interest as a visual summary, the persistence diagram captures all the information about the persistent homology groups of the data [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF].

Given two diagrams D( f i ) and D( f j ), a pointwise distance, noted d q (with q > 0), can be introduced in the 2D birth/death space between two points p i = (x i , y i ) ∈ D( f i ) and p j = (x j , y j ) ∈ D( f j ):

d q (p i , p j ) = (|x j -x i | q + |y j -y i | q ) 1/q = p i -p j q .
(1) Fig. 5. Distance and geodesic computation between two persistence diagrams D( f i ) and D( f j ), given the metric W D 2 . The sets P i and P j are shown in transparent (a). The matching induced by the optimal partial assignment φ is shown with dashed lines (c). For two diagrams, the Wasserstein barycenter (b) is given, thanks to the d 2 distance in the birth/death space, by the arithmetic mean of the matched points. Then, the linear interpolation of the matchings (c) describes a geodesic [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF].

By convention, d q (p i , p j ) is set to zero if both p i and p j exactly lie on the diagonal (x i = y i and x j = y j ). Let P i be a subset of the off-diagonal points of D( f i ) and P i its complement (i.e. the other off-diagonal points of D( f i ) not in P i ). Let (φ , P i , P j ) be a partial assignment between D( f i ) and D( f j ), i.e. a bijective map between P i and a subset of offdiagonal points P j of D( f j ), with complement P j (Fig. 5a). Then, the L q -Wasserstein distance, noted W D q , can be introduced as:

W D q D( f i ), D( f j ) = min (φ ,P i ,P j )∈Φ ∑ p i ∈P i d q p i , φ (p i ) q (2) + ∑ p i ∈P i d q p i , ∆(p i ) q (3) + ∑ p j ∈P j d q ∆(p j ), p j q 1/q (4)
where Φ is the set of all possible partial assignments mapping each point

p i ∈ D( f i ) to a point φ (p i ) = p j ∈ D( f j ) (line 2), or to its di- agonal projection, ∆(p i ) = ( x i +y i 2 , x i +y i 2 )
, denoting the removal of the corresponding feature from D( f i ) or D( f j ) (lines 3 and 4). Intuitively, the Wasserstein metric optimizes a matching between the two diagrams, and evaluates their distance given the resulting mismatch. In practice, D( f i ) and D( f j ) are augmented into D ( f i ) and D ( f j ) [START_REF] Kerber | Geometry helps to compare persistence diagrams[END_REF], by injecting the diagonal projections of one diagram into the other (Fig. 5a):

D ( f i ) = D( f i ) ∪ {∆(p j ) | p j ∈ P j } D ( f j ) = D( f j ) ∪ {∆(p i ) | p i ∈ P i }.
This augmentation (Fig. 5a) preserves the distance, while making the assignment problem balanced, and thus easily solvable with traditional algorithms [START_REF] Bertsekas | A new algorithm for the assignment problem[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] (with P i = D ( f i ), P j = D ( f j ) and P i = P j = / 0). Given a set S D = {D( f 1 ), . . . , D( f N )} of persistence diagrams, let F(D, α) be the Fréchet energy of the set, under the metric

W D 2 , with the coefficients α = {α 1 , α 2 , . . . , α N }, such that α i ∈ [0, 1] and ∑ i α i = 1: F(D, α) = ∑ D( f i )∈S D α i W D 2 D, D( f i ) 2 . ( 5
)
Then the diagram D * ∈ D (where D is the space of persistence diagrams) which minimizes F(D, α) is called the Wasserstein barycenter of the set S D (or its Fréchet mean under the metric W D

2 ). In practice, the coefficients α i are all set to the same value (α i = 1/N, ∀i ∈ {1, . . . , N}). When N = 2 and α 1 = α 2 = 0.5 (Fig. 5b), D * becomes a midpoint between D( f i ) and D( f j ) and the set of possible values for α 1 and α 2 (Fig. 5c) describes a geodesic in D (i.e. length minimizing path) with regard to the L 2 -Wasserstein metric [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF].

Merge trees

The join tree, noted T -( f i ), is a visual summary of the connected components of f i -1 -∞ (w) [START_REF] Carr | Computing contour trees in all dimensions[END_REF]. It is a 1-dimensional simplicial complex defined as the quotient space T -( f i ) = M / ∼ by the equivalence relation ∼ which states that p 1 and p 2 are equivalent if

f i (p 1 ) = f i (p 2 ) and if p 1 and p 2 belong to the same connected component of f i -1 -∞ f i (p 1 ) .
The split tree (Fig. 2), noted T + ( f i ), is defined symmetrically and describes the connected components of the super-level set f i -1

+∞ (w).

Each of these two directed trees is called a merge tree, noted generically T ( f i ) in the following. Intuitively, these trees track the creation of connected components of the sub (or super) level sets at their leaves, and merge events at their interior nodes. These trees are often visualized according to a persistence-driven branch decomposition [START_REF] Pascucci | Multi-resolution computation and presentation of contour trees[END_REF], to make the persistence pairs captured by the tree stand out. In this context, a persistent branch is a monotone path on the tree connecting the nodes corresponding to the creation and destruction (according to the Elder rule, Sec. 2.2) of a connected component of sub (or super) level set. Then, the branch decomposition provides a planar layout of the merge tree, where each persistent branch is represented as a vertical segment. The branch decomposition tree (BDT), noted B( f i ), is a directed tree whose nodes are the persistent branches captured by the branch decomposition and whose arcs denote adjacency relations between them in the merge tree. In Fig. 2, the BDTs (top right insets) can be interpreted as the dual of the branch decompositions (bottom right insets, with matching colors): each vertical segment in the branch decomposition (bottom) corresponds to a node in the BDT (top) and each horizontal segment (bottom, denoting an adjacency relation between branches) corresponds to an arc in the BDT. Intuitively, the BDT, like the persistence diagram, describes the population of (extremum-saddle) persistence pairs present in the data. However, unlike the persistence diagram, it additionally captures adjacency relations between them.

WASSERSTEIN DISTANCES BETWEEN MERGE TREES

This section introduces our new distance metric between merge trees, which is specifically designed for the subsequent computation of geodesics (Sec. 4) and barycenters (Sec. 5). For this, we bridge the gap between the edit distance between merge trees [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] and existing work addressing the computation of geodesics and barycenters for persistence diagrams according to the L 2 -Wasserstein distance [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF].

Overview

The end goal of our work is the computation of barycenters of merge trees. For this, we extend the edit distance D E [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] (formalized in Appendix 1, additional material), to make it fit the optimization strategy used for barycenters of persistence diagrams [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF]. Our key idea consists in transforming D E such that it becomes strictly equivalent to the L 2 -Wasserstein distance of persistence diagrams, but given a restricted space of possible assignments, constrained by the structure of the input trees T ( f i ) and T ( f j ), hence its name Wasserstein distance between merge trees. Then, thanks to this compatibility with the L 2 -Wasserstein distance, the assignments resulting from our metric can be directly used for interpolation-based geodesic and barycenter computations (Secs. 4 and 5). Overall, our strategy involves four major modifications to the edit distance D E [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF], detailed in the remainder of this section:

1. To consider assignments between persistence pairs instead of merge tree nodes, we consider an edit distance between the BDTs B( f i ) and B( f j ) (Sec. 2.3) instead of the input merges trees T ( f i ) and T ( f j ) (as done with D E ). This is described in Sec. 3.2. 2. We constrain the assignment search space to the space of rooted partial isomorphisms. Specifically, similarly to D E , we enforce the assignment of disjoint subtrees of B( f i ) to disjoint subtrees of B( f j ). Moreover, in contrast to D E , we additionally extend this constraint by enforcing the destruction of entire subtrees upon the destruction of their root. These two constraints together enforce assignments describing isomorphisms between rooted subtrees of B( f i ) and B( f j ). Such isomorphisms pave the way for interpolationbased geodesics. This is described in Secs. 3.2, 3.3 and 4.1. 3. We introduce a cost model based on the Euclidean distance d 2 to enable geodesic computation by linear interpolation of the assignments in the 2D birth/death space. This is described in Secs. 3.2 and 4.1. 4. We finally extend our metric with a local normalization term, which enforces nested birth-death values, along the interpolation of the assignments, for nested branches. This is described in Sec. 4.2.

Definition and properties

Given two input merge trees, T ( f i ) and T ( f j ), we first consider their BDTs B( f i ) and B( f j ) (Sec. 2.3). Let B i be a subset of the nodes of B( f i ) and B i its complement. Note that each node in B i corresponds to a persistence pair of D( f i ). Let (φ , B i , B j ) be a partial assignment between B i and a subset B j of the nodes of B( f j ) (with complement B j ). Then we introduce the L 2 -Wasserstein distance W T 2 between the BDTs B( f i ) and B( f j ) of the merge trees T ( f i ) and T ( f j ) as:

W T 2 B( f i ), B( f j ) = min (φ ,B i ,B j )∈Φ ∑ b i ∈B i γ b i → φ (b i ) 2 (6) + ∑ b i ∈B i γ(b i → / 0) 2 (7) + ∑ b j ∈B j γ( / 0 → b j ) 2 1/2 ( 8
)
where Φ is the space of constrained partial assignments mapping disjoints subtrees of B( f i ) to disjoint subtrees of B( f j ), and mapping entire subtrees to the empty tree / 0 if their root is itself mapped to / 0. Then, given the

k th direct child of b i , noted b k i , it follows that b k i either maps through φ to a direct child of φ (b i ) ∈ B( f j ) (then b i , b k i ∈ B i
) or to the empty tree / 0 (then the subtree rooted in b k i , noted B( f i , b k i ), belongs to B i ). This further implies that the rooted subtrees B i ⊆ B( f i ) and B j = φ (B i ) ⊆ B( f j ) are isomorphic and we call (φ , B i , B j ) a rooted partial isomorphism. Unlike D E (see Appendix 1) but similarly to W D 2 (Eq. 2), the cost of each operation (mapping, line 18, destruction, line 18, and creation, line 19) is squared, and the square root of the sum of the squared costs is considered as the overall distance.

Next, we define the edit costs as follows (we recall that each branch b i ∈ B( f i ) exactly coincides with a persistence pair p i ∈ D( f i )):

γ b i → φ (b i ) = d 2 b i , φ (b i ) γ(b i → / 0) = d 2 b i , ∆(b i ) (9)
γ( / 0 → b j ) = d 2 ∆(b j ), b j .
Note that the expression of the L 2 -Wasserstein distance W T 2 between merge trees (Eq. 6) is therefore identical to the expression of the Wasserstein distance between persistence diagrams (Eq. 2) for q = 2, at the notable exception of the search space of the partial assignments Φ ⊂ Φ, which is constrained to rooted partial isomorphisms. W T 2 is indeed a distance metric (the proof is included in Appendix 2, supplemental material): it is non-negative and symmetric, it preserves the identity of indiscernibles as well as the triangle inequality. Moreover, since

Φ ⊂ Φ, it follows that W T 2 B( f i ), B( f j ) ≥ W D 2 D( f i ), D( f j )
, which was one of the main motivations of our work (i.e. to exploit the merge tree to define a more discriminative metric, Fig. 3). Similarly to Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF], we mitigate saddle swap instabilities in a preprocessing step, by merging adjacent saddles in the input trees if their difference in scalar value is smaller than a threshold ε 1 ∈ [0, 1] (relative to the largest difference between adjacent saddles). Then, when

ε 1 = 1, it follows that W T 2 B( f i ), B( f j ) = W D 2 D( f i ), D( f j )
. This simple merging strategy significantly improves the practical stability of W T 2 , as empirically studied in Sec. 7.2 (Fig. 14).

Computation

This section describes our algorithm for the recursive exploration of the search space Φ (Eq. 6). It is based on the same recursive traversal as Zhang's algorithm [121], which we simplify as our search space is significantly more constrained. Specifically, as detailed in Appendix 3, our distance evaluation between subtrees (Eq. 11) involves fewer solutions and it is restricted to subtrees rooted at identical depth only.

Given the subtree

B( f i , b) of B( f i ) (rooted in b) and b k the k th direct child of b in B( f i , b), the distance between the subtree B( f i , b)
and the empty tree / 0 is then obtained recursively by:

W T 2 B( f i , b), / 0 = γ(b → / 0) 2 + ∑ k W T 2 B( f i , b k ), / 0 2 1/2
. (10) The first step of our algorithm consists in evaluating W T 2 B( f i , b), / 0 with Eq. 10 for all branches b ∈ B( f i ) (and similarly for B( f j )).

Next, let

F ( f i , b) be the forest of b in B( f i ): F ( f i , b)
is the set of all the subtrees rooted at the k direct children of b:

F ( f i , b) = {B( f i , b 1 ), B( f i , b 2 ), . . . , B( f i , b k )}.
Then the distance between two subtrees B( f i , b i ) and B( f j , b j ) is set to +∞ when b i and b j have distinct depths (gray crosshatching lines in Fig. 6, left). Otherwise (spheres in Fig. 6, left), it is obtained recursively by:

W T 2 B( f i , b i ), B( f j , b j ) = γ(b i → b j ) 2 (11)
+ W T 2 F ( f i , b i ), F ( f j , b j ) 2 1/2
.

Let F i be a subset of the forest F ( f i , b i ) and F i its complement. The distance between two forests is then given recursively by:

W T 2 F ( f i , b i ), F ( f j , b j ) = min (φ ,F i ,F j )∈Φ ∑ f i ∈F i W T 2 f i , φ ( f i ) 2 + ∑ f i ∈F i W T 2 f i , / 0 2 + ∑ f j ∈F j W T 2 f j , / 0 2 1/2
where (φ , F i , F j ) becomes the solution of a local, partial assignment problem between forests, mapping F i to a subset F j ∈ F ( f j , b j ) (with complement F j ) or to the empty tree / 0, and which can be solved with traditional assignment algorithms [START_REF] Bertsekas | A new algorithm for the assignment problem[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] (see Sec. 2.2). Then, the overall distance W T 2 between the two input merge trees is obtained by estimating Eq. 11 at the roots of B( f i ) and B( f j ), and solving recursively the local assignment problems between forests (the recursion returns are illustrated with arrows in Fig. 6). Note that if ε 1 = 1 (Sec. 3.2), all the branches of B( f i ) and B( f j ) get attached to the roots and the recursive local assignment problems between forests (above) become only one, large, assignment problem between all branches. Thus, when ε 1 = 1, this algorithm indeed becomes strictly equivalent to the resolution of the assignment problem involved in W D 2 (Sec. 2.2).

Parallelism

Similarly to Zhang [121], Eqs. 10 and 11 can be estimated recursively.

To avoid redundant computations, the distances between the forests

F ( f i , b i ) and F ( f j , b j ) are stored at the entry (b i , b j ) of a matrix F (of size |B( f i )| × |B( f j )|)
, while the distances between the subtrees B( f i , b i ) and B( f j , b j ) (used within the assignment problems between higher forests) are stored in a matrix T (of the same size, see Fig. 6).

In our work, we additionally express this computation in terms of tasks, to leverage task-based shared memory parallelism. In particular, we initiate a task for each independent term of Eqs. 10 and 11, which is ready for computation (see Fig. 6), as further detailed in Appendix 4. Then, the number of parallel tasks is initially bounded by the number of leaves in the input BDTs (which is typically much larger than the number of cores) and progressively decreases during the computation.

WASSERSTEIN GEODESICS BETWEEN MERGE TREES

This section introduces our approach for the efficient computation of geodesics between merge trees, according to the metric W T 2 (Sec. 2). For this, we leverage the rooted partial isomorphism resulting from the distance computation, as well as linear interpolations of the matchings, as introduced for persistence diagrams [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF].

Definition and properties

Given two input merge trees T ( f i ) and T ( f j ), our approach to geodesic computation (Fig. 7) simply consists in linearly interpolating the rooted partial isomorphism (φ , B i , B j ) resulting from the optimization involved in the computation of W T 2 B( f i ), B( f j ) (Eq. 6). In particular, given the two BDTs B( f i ) and B( f j ), the interpolated BDT, noted B α (

f i → f j ) with α ∈ [0, 1] such that B 0 ( f i → f j ) = B( f i ) and B 1 ( f i → f j ) = B( f j ), is obtained by considering the union of: 1. the linear interpolation B α ⊆ B α ( f i → f j ), between the nodes B i ⊆ B( f i
) and these of B j ⊆ B( f j ), given the isomorphism φ (the trees B i , B j and B α are then isomorphic, Fig. 7):

b(α) = (1 -α)b + αφ (b) ∀b ∈ B i (12)
2. the linear interpolation of the destruction of the subtrees B i , noted

B i α ⊆ B α ( f i → f j ) (B i and B i α are also isomorphic): b(α) = (1 -α)b + α∆(b) ∀b ∈ B i (13)
3. the linear interpolation of the creation of the subtrees B j , noted B j α ⊆ B α ( f i → f j ) (B j and B j α are also isomorphic):

b(α) = (1 -α)∆(b) + αb ∀b ∈ B j . ( 14
)
Similarly to the distance W D 2 between persistence diagrams, since the edit costs involved in the distance W T 2 are Euclidean distances in the birth/death space (Eq. 9), the interpolated branches b(α) of B α ( f i → f j ) can be efficiently computed with the simple linear interpolations described above. As detailed in Appendix 5, the resulting interpolated tree B α ( f i → f j ) is indeed on a geodesic given W T 2 .

From branch decomposition trees to merge trees

The previous section described the computation of geodesics between BDTs, given W T 2 . In this section, given an interpolated BDT B α ( f i → f j ), we describe how to retrieve the corresponding merge tree T α ( f i → f j ) (i.e. a merge tree whose BDT is indeed equal to B α ( f i → f j )).

A requirement for an arbitrary BDT B to be the valid BDT of a merge tree T is that subtrees of B need to respect a nesting condition on their birth/death (i.e. x, y) values (to respect the Elder rule, Sec. 2.2). In particular, given a direct child b k α of a branch b α ∈ B α (

f i → f j ), we need to guarantee that [x k b α , y k b α ] ⊆ [x b α , y b α ].
While this is guaranteed by construction for the subset B α ⊆ B α ( f i → f j ) (B α is isomorphic to B i and B j ), this is not necessarily the case for the subsets of B α ( f i → f j ) involved in subtree creation or destruction (B i α and B j α , Sec. 4.1). In particular, since the branches involved in destructions map independently to the diagonal (Eq. 13), it is possible that the above nesting condition is not respected along their interpolation. This is shown in Fig. 8c (red interpolation), where the resulting merge tree, T α ( f i → f j ), becomes disconnected and hence invalid (i.e. B α ( f i → f j ), Fig. 8d, is connected and not equal to the BDT of T α ( f i → f j ) from Fig. 8c).

In the following, we introduce a pre-processing step for the trees B( f i ) and B( f j ) (together with its inverse post-processing step), which we call local normalization, which addresses this issue and guarantees the above nesting condition, even in case of destruction/creation. Given a direct child b i k of a branch b i ∈ B( f i ), we consider the following local, birth/death normalization

N (b i k ) = N x (b i k ), N y (b i k ) : N x (b i k ) = (x b i k -x b i )/(y b i -x b i ) N y (b i k ) = (y b i k -x b i )/(y b i -x b i ).
Once this pre-process is recursively completed, the Wasserstein distance W T 2 between the locally normalized BDTs, noted N B( f i ) and N B( f j ) is computed as described in Sec. 3.3. Then, the interpolation of the locally normalized BDTs, noted N B α ( f i → f j ) , is evaluated as described in Sec. 4.1. Next, the local normalization is recursively reverted to turn N B α (

f i → f j ) back into B α ( f i → f j ), by explicitly evaluating N -1 (b N α ) for each branch b N α ∈ N B α ( f i → f j )
. Now, even in case of branch destruction, by construction, the birth/death interval of each interpolated branch

N (b α ), noted [N x (b α ), N y (b α )], is included in [0, 1] (since ∆ N (b i ) ⊆ [N x (b i ), N y (b i )] ⊆ [0, 1]
). Therefore, after reverting the local normalization, we have the guarantee that

[x k b α , y k b α ] ⊆ [x b α , y b α ]
for all the branches b α of B α ( f i → f j ). At this stage, B α ( f i → f j ) indeed respects the nesting condition on the birth/death values of all its subtrees. Then, given the dual relation between merge trees and BDTs, the merge tree T α ( f i → f j ) can be simply obtained by creating a vertical branch for each node b α of B α ( f i → f j ) and connecting them according the arcs of B α ( f i → f j ), as illustrated in Fig. 8 (right). The distance W T 2 between N B( f i ) and N B( f j ) then still describes a metric between B( f i ) and B( f j ), such that B α ( f i → f j ) is indeed on a geodesic (see Appendix 6).

Note that the local normalization shrinks all the input branches to the interval [0, 1], irrespective of their original persistence. To mitigate this effect, we introduce a pre-processing step on the input BDTs, Fig. 9. Visual comparison between our barycenter (Sec. 5) and the 1center of Yan et al. [START_REF] Yan | Source Code for a Structural Average of Labeled Merge Trees for Uncertainty Visualization[END_REF][START_REF] Yan | A structural[END_REF]. Left: an ensemble is created with an outlier member f j (red, 7 persistent branches) and 10 noisy versions of a field f i (4 persistent branches). Right: planar view of the ensemble computed by multi-dimensional scaling of W T 2 . The barycenter computed with our approach (cyan) is more similar to the merge trees of f i (same number and persistence of large branches) and hence better captures the overall trend of the ensemble, despite the presence of the outlier f j (red sphere). 2 ). This iterated removal of trees highlights key frames in the sequence (d-g) corresponding to key phases of an asteroid impact simulation [START_REF] Patchett | The IEEE SciVis Contest[END_REF]: initial state (black, time steps 1-5), approach (light green, 6-10), impact (blue, 11-15), aftermath (light blue, [START_REF] Bhatia | Topoms: Comprehensive topological exploration for molecular and condensed-matter systems[END_REF][START_REF] Bhatia | Flow visualization with quantified spatial and temporal errors using edge maps[END_REF][START_REF] Biasotti | Reeb graphs for shape analysis and applications[END_REF][START_REF] Bin Masood | An Overview of the Topology ToolKit[END_REF][START_REF] Bock | TopoAngler: Interactive Topology-Based Extraction of Fishes[END_REF]. In contrast, a similar greedy optimization based on the distance W D 2 between persistence diagrams (red curve) fails at capturing the impact phase (blue) of this sequence. which moves, up the trees, subtrees rooted at branches with a relative persistence smaller than ε 3 , until their persistence relative to their parent becomes smaller than a threshold ε 2 . This has the practical effect of reducing the normalized persistence of small branches corresponding to small features. Overall, ε 1 , ε 2 and ε 3 are the only parameters of our approach and we use a unique, default set of values (ε 1 = 0.05, ε 2 = 0.95 and ε 3 = 0.9) in our experiments (Sec. 7). In the remainder, we will consider that all the input BDTs are normalized this way.

WASSERSTEIN BARYCENTERS OF MERGE TREES

This section introduces our approach for the computation of barycenters of merge trees, for the metric W T 2 (Sec. 2). The resulting barycenters will serve as core tools for clustering ensembles of merge trees (Sec. 6).

Definition

Let S B = {B( f 1 ), B( f 2 ), . . . , B( f N )} be a set of N BDTs. Similarly to Eq. 5, The Fréchet energy, under the metric W T 2 , is given by:

F(B) = ∑ B( f i )∈S B W T 2 B, B( f i ) 2 .
We call a Wasserstein barycenter of S B , a BDT B * ∈ B (where B is the space of BDTs) which minimizes F(B). It is a centroid of the set, i.e. a tree which minimizes the sum of its distances to the set.

Computation

Our distance W T 2 (Sec. 2) is identical to W D 2 , but with a smaller search space, restricted to rooted partial isomorphisms. This enabled an extension of interpolation-based geodesics from persistence diagrams to merge trees (Sec. 4). Given these two components, the strategy presented by Turner et al. [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF] for minimizing the Fréchet energy over the space of persistence diagrams can be directly extended to our framework. For this, we consider an algorithm that resembles a Lloyd relaxation [START_REF] Lloyd | Least squares quantization in PCM[END_REF], and which alternates an (i) assignment and an (ii) update procedure. First, the candidate B is initialized at an arbitrary tree of S B . Then the assignment step (i) computes an optimal assignment (φ i , B i , B B ) between B and each tree B( f i ) ∈ S B . Next, the update step (ii) updates the candidate B to a position in B which minimizes F(B) under the current set of assignments (φ i , B i , B B ) i=1,...,N . This is achieved by moving each branch b ∈ B (in the birth/death space) to the arithmetic mean of the assignments (by generalizing the interpolation defined in Eqs. 12, 13, and 14, to more than two trees):

b ← 1 N ∑ i=1,...,N      φ i (b) if b ∈ B i ∆(b) if b ∈ B i b if b ∈ B B .
This overall assignment/update sequence is then iterated (as discussed in Appendix 7, each iteration of this sequence decreases the Fréchet energy constructively). In our implementation, the algorithm stops and returns the barycenter estimation B * when the Fréchet energy decreased by less than 1% between two consecutive iterations. Given B * , we obtain its dual merge tree T * as described in Sec. 4.2. Fig. 9 illustrates a barycenter computed with this strategy for a toy example.

Parallelism

The N assignment problems (between the candidate B and the trees of the set S B , Sec. 5.2) are independent and can be computed in parallel. However, this naive strategy is subject to load imbalance, as the input trees can have different sizes. Hence, each iteration would be bounded by the sequential execution of the largest of the N assignment problems.

We address this issue by leveraging the task-based parallelization of our distance computation algorithm (Sec. 3.4). In particular, we use a single task pool for all of the N assignment problems. Then, the task environment picks up at runtime the tasks to compute irrespective of their tree of origin, and place them on different threads. This fine scheduling granularity has the beneficial effect of triggering the execution of the tasks of a new assignment problem while a first problem is reaching completion (and thus exploiting less threads, Sec. 3.4). This improves thread load imbalance and thus increases the overall parallel efficiency.

APPLICATIONS

The section illustrates the utility of our contributions (distances, geodesics, and barycenters) in concrete visualization tasks (Fig. 1).

Branch matching for feature tracking

Our distance (Sec. 3.3) relies on the optimization of a partial isomorphism between the input BDTs. Then, the resulting matchings can be used to track features in time-varying data, as studied for persistence diagrams [START_REF] Soler | Lifted Wasserstein matcher for fast and robust topology tracking[END_REF]. Fig. 11 illustrates this on a temporal sequence (SciVis contest 2008 [START_REF] Taylor | The IEEE SciVis Contest[END_REF]). Since W D 2 considers persistence pairs individually, it can generate inconsistent matchings with a typical incorrect crossing in the feature tracking (already visible on synthetic data, Fig. 7). Our distance W T 2 improves this aspect by better preserving the global structure of the data, thanks to our more constrained, merge-tree driven, Fig. 12. Three members (a) of an acquired ensemble, corresponding to distinct volcanic eruptions [START_REF] Feng | The IEEE SciVis Contest[END_REF]. Our clustering approach correctly assigns the members to each cluster (b, distinct colors in the planar view, generated in a post-process by multi-dimensional scaling of W T 2 ). Our centroids (larger spheres in the planar view) provide a visual summary of the features of interest (matching colors) for each cluster. assignment search space. Overall, our matchings provide visual hints to the users, to help them relate features from distinct time steps.

Geodesics for temporal reduction

The topological analysis of time-varying data typically requires the computation of a topological representation, for instance a merge tree, for each time step. Although merge trees are usually orders of magnitude smaller than the original data, the resulting sequence of merge trees can still represent considerable amounts of data. To address this issue, we exploit our geodesic computation (Sec. 4) for the reduction of temporal sequences of merge trees. In particular, we greedily remove from the sequence, one by one, the merge trees which can be accurately reconstructed by simple geodesic computation, until the sequence only contains a target number of merge trees (see the detailed algorithm in Appendix 8). This enables the reliable visualization of time-varying sequences of merge trees at greatly reduced storage costs. The remaining merge trees (white background, Fig. 10) correspond to key frames of the sequence, i.e. time steps of particular significance in terms of the features of interest. In contrast, a similar strategy based on persistence diagram interpolation (red curve) fails at identifying a key frame in one of the key phases of the sequence (impact, in blue). Also, note that the reduced merge trees (reconstructed with geodesics, blue background) are visually highly similar to the trees from the input sequence.

Barycenters for topological clustering

To understand the main trends within an ensemble, in terms of features of interest, it may be desirable to cluster the ensemble by grouping members with a similar topological profile. For this, we adapt the k-means algorithm [START_REF] Celebi | A comparative study of efficient initialization methods for the k-means clustering algorithm[END_REF][START_REF] Elkan | Using the triangle inequality to accelerate k-means[END_REF] to the problem of clustering merge trees. In particular, this can be easily achieved by using our merge tree barycenter computation algorithm (Sec. 5) as the centroid estimation routine of k-means, and by using W T 2 (Sec. 2) to measure the distance between merge trees. Note that in practice, our entire computational framework is implemented in this single clustering algorithm (with a unique task pool), as the above clustering generalizes the barycenter problem (k = 1) as well as the geodesic and distance problems (N = 2).

Figs 12 and13 present clustering examples obtained with this strategy on an acquired [START_REF] Feng | The IEEE SciVis Contest[END_REF] and cosmology ensemble [START_REF] Hentschel | The IEEE SciVis Contest[END_REF]. In both cases, our approach correctly assigns the members to each cluster. Moreover, the centroids computed by our algorithm provide a visual summary of the features of interest found in each cluster, enabling global overviews (Figs. 12, right, and 13, bottom) summarizing the topological profile of each of the main trends found in the ensemble. In both figures, the tree branches of the ensemble members are automatically colored with the color of their matched centroid branch. This matching visualization enables users to visually relate the centroid to concrete features in the data (Fig. 12) and to compare matching features across multiple members (i.e. which have been matched to the same centroid branch, Fig. 13). Then, the centroid, in addition to being a visual summary, also acts as a reference point for the visual comparison of ensemble members.

RESULTS

This section presents experimental results obtained on a computer with two Xeon CPUs (3.2 GHz, 2x10 cores, 96GB of RAM). The input merge trees were computed with FTM [START_REF] Gueunet | Task-Based Augmented Contour Trees with Fibonacci Heaps[END_REF] and pre-processed to discard noisy features (persistence simplification threshold: 0.25% Fig. 13. Eight members (top) of a cosmology ensemble [START_REF] Hentschel | The IEEE SciVis Contest[END_REF], and their merge trees (persistent maxima are displayed with matching colors in the data). Our clustering approach correctly assigns the members to each cluster (distinct colors in the bottom planar view, generated in post-process by multi-dimensional scaling of W T 2 ). Our centroids (large spheres, bottom) provide a visual summary which is well representative of the trees in the cluster (same number and persistence of large branches, automatically color-coded based on their matching to their centroid).

of the data range). We implemented our approach in C++ (with the OpenMP task runtime), as modules for TTK [START_REF] Bin Masood | An Overview of the Topology ToolKit[END_REF][START_REF] Tierny | The Topology ToolKit[END_REF].

Our experiments were performed on a variety of simulated and acquired 2D and 3D ensembles used in previous work [START_REF] Favelier | Persistence Atlas for Critical Point Variability in Ensembles[END_REF] (vorticity and sea surface height) or extracted from past SciVis contests: 2004 (wind velocity magnitude [START_REF] Wang | The IEEE SciVis Contest[END_REF]), 2006 (wavefront velocity magnitude [START_REF] Olsen | The IEEE SciVis Contest[END_REF]), 2008 (ion concentration [START_REF] Taylor | The IEEE SciVis Contest[END_REF]), 2014 (sulfur dioxide concentration [START_REF] Feng | The IEEE SciVis Contest[END_REF]), 2015 (dark matter density [START_REF] Hentschel | The IEEE SciVis Contest[END_REF]), 2016 (salt concentration [START_REF] Garth | The IEEE SciVis Contest[END_REF]), 2017 (pressure [START_REF] Wischgoll | The IEEE SciVis Contest[END_REF] ), 2018 (matter density [START_REF] Patchett | The IEEE SciVis Contest[END_REF]). A detailed specification of these ensembles is provided in Appendix 9.

Time performance

The time complexity of our algorithm for exploring the search space of W T 2 (Sec. 3.3) is similar to that of the edit distance [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF]121]. It takes O(|B| 2 ) steps in practice, with |B| the number of nodes in the input BDTs (in our implementation, each local forest assignment problem is solved with the efficient Auction approximation [START_REF] Bertsekas | A new algorithm for the assignment problem[END_REF] with default parameters). Once W T 2 B( f i ), B( f j ) is computed, the computation of a point on the geodesic (Sec. 4) between B( f i ) and B( f j ) is obtained in O(|B|) steps. Regarding our barycenter computation algorithm (Sec. 5), each of its iterations takes O(N|B| 2 ) steps. Tab. 1 evaluates the practical time performance of our computational framework for the barycenter computation (which includes itself distance and geodesic computations). In sequential mode, we observe that the running time is indeed a function of the number of ensemble members (N) and the average size of the trees (B). It is slightly slower for W T 2 than for W D 2 , but runtimes remain comparable overall. In parallel, speedups are the most important for the largest examples. However, the iterative nature of our barycenter optimization algorithm seems to limit parallel efficiency globally (the end of each iteration still constitutes a strong synchronization). For the smaller examples, the cost of the task runtime seems to become non-negligible in comparison to the actual computation, resulting in moderate speedups. Still, our parallelization significantly reduces runtimes overall, with less than 3 minutes of computation on average and at most 15 minutes for the largest examples.

Framework quality

W T 2 is indeed a distance metric (Appendix 2). It is more discriminative than W D 2 (i.e. W T 2 ≥ W D 2 ,
Sec. 2, Fig. 3). Fig. 14 evaluates empirically its stability. For this, given a scalar field f i , a noisy version f j is created such that f if j ∞ ≤ ε, for increasing values of ε. Then, we observe the evolution of W T 2 B( f i ), B( f j ) , as a function of ε (Fig. 14, right), to estimate how W T 2 varies under input perturbations. For

ε 1 = 1, we have W T 2 = W D 2 (Sec.
2) and the curve evolves nearly linearly (W D 2 is stable [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF]). For other ε 1 values, the curves indicate clear transition points (colored dots) before which W T 2 evolves nearly linearly too. This indicates that for reasonable noise levels (smaller than the ε value of each transition point, vertical lines), W T 2 is also stable and that only mild increases of ε 1 result in fast shifts of these transition points (to an accepted noise level of 64% at ε 1 = 0.15). This illustrates overall that the stability of W T 2 can indeed be controlled with ε 1 and that small values already lead to stable results for reasonable noise levels. A detailed empirical analysis of the other two parameters of our approach (ε 2 , ε 3 , Sec. 4.2) is provided in Appendix 10 (supplemental material).

Next, we study the practical relevance of W T 2 by evaluating our clustering performance. For this, each ensemble of Tab. 1 is associated with a ground truth classification (distinct phases of a time-varying phenomenon, distinct input parameters, etc), by following the companion specifications [START_REF] Feng | The IEEE SciVis Contest[END_REF][START_REF] Garth | The IEEE SciVis Contest[END_REF][START_REF] Hentschel | The IEEE SciVis Contest[END_REF][START_REF] Olsen | The IEEE SciVis Contest[END_REF][START_REF] Patchett | The IEEE SciVis Contest[END_REF][START_REF] Taylor | The IEEE SciVis Contest[END_REF][START_REF] Wang | The IEEE SciVis Contest[END_REF][START_REF] Wischgoll | The IEEE SciVis Contest[END_REF]. Clustering performance is evaluated with accepted scores, namely the normalized mutual information and adjusted rand index (NMI, ARI). When using our barycenters (Sec. 5), our clustering approach (Sec. 6.3) achieves a perfect classification for all ensembles (NMI = ARI = 1). These scores decrease to NMI = 0.78 and ARI = 0.69 on average when using, within k-means, a barycenter of persistence diagrams [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF] (ε 1 = 1), and to NMI = 0.73 and ARI = 0.56 when using the 1-center of Yan et al. [START_REF] Yan | A structural[END_REF] (obtained with the authors' implementation [START_REF] Yan | Source Code for a Structural Average of Labeled Merge Trees for Uncertainty Visualization[END_REF], using leaf labels generated by our distance computation, Sec. 2). This simply confirms experimentally that 1-centers in general are not suited for clustering tasks. A standard clustering approach (multi-dimensional scaling to kD followed by k-means) using the distance D E [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] achieves lower average scores than our approach, with NMI = 0.89 and ARI = 0.85 on average. Overall, this confirms that D E induces more discriminative classifiers than W D 2 , and that our metric W T 2 further improves that. Fig. 15 shows the evolution of the Fréchet energy for our barycenter algorithm (Sec. 5) for various ε 1 values. In practice, the algorithm stops when the Fréchet energy decreases by less than 1% between consecutive iterations, which occurs early in the process.

Fig. 9 provides a visual comparison between our barycenter and the 1-center of Yan et al. [START_REF] Yan | A structural[END_REF] (obtained with the authors' implementation [START_REF] Yan | Source Code for a Structural Average of Labeled Merge Trees for Uncertainty Visualization[END_REF], using leaf labels generated by our distance computation, Sec. 2). This figure confirms the general sensitivity in practice of 1-centers to outliers, and the ability of barycenters to better capture the main trends in the ensemble. From a qualitative perspective, our framework enables the computation of faithful interpolations of merge trees: the reconstructed trees, blue background (Fig. 10), are visually very similar to the input trees. Moreover, our framework produces barycenters (Figs. 1, 12, and 13) which capture well the main features of the input ensemble: for each cluster, the resulting centroid is visually similar to the input trees of the cluster (same number and persistence of large branches). Then, our clustering framework, coupled with our centroids, provides a faithful visual summary of the features of interest, for each of the main trends (i.e. for each cluster) found in the ensemble.

Limitations

The search space associated with our metric W T 2 is constrained to rooted partial isomorphisms. Then, if a matching exists between two BDTs (i.e. if they are not both destructed when optimizing W T

2 ), it has to match their roots together. In other words, W T 2 nearly always matches the most persistent branch of the two trees together, which might be too restrictive (in particular for feature tracking applications). Note however, that W D 2 behaves equivalently: the most persistent branch of B( f i ) corresponds to the component of f i -1 -∞ (w) created in the global minimum of f i , which in principle has infinite persistence and which is typically treated separately when evaluating W D 2 . Similarly to Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF], saddle swap instabilities are handled in our approach by a pre-processing step which merges adjacent saddles (controlled by ε 1 ). An alternative would consist in exploring the space of all possible branch decompositions (not necessarily persistence-driven), as studied by Beketayev et al. [START_REF] Beketayev | Measuring the distance between merge trees[END_REF]. However, the search space would then become significantly larger. Moreover, the nesting of birth/death values within the BDTs would no longer be guaranteed, which is however a key property which we exploit in our framework (Sec. 4). When computing barycenters of persistence diagrams, Vidal et al. [START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] showed that the optimization could be drastically accelerated by introducing persistence pairs progressively along the iterations, while implicitly maintaining previous assignments at each initialization. We leave the study of such a progressive strategy to future work, although the fact that W T 2 handles many small assignment problems (unlike W D

2 ) indicates that such a strategy may result in only modest gains for merge trees. Fig. 14 provides an empirical evaluation of the stability of W T 2 . Similarly to Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF], we believe that the theoretical investigation of the stability of W T 2 goes beyond the scope of this paper and we leave it for future work.

CONCLUSION

In this paper, we presented a computational framework for the estimation of distances, geodesics and barycenters of merge trees, with applications to feature tracking, temporal reduction and ensemble clustering and summarization. Our approach filled the gap between the edit distance [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] and existing optimization frameworks for persistence diagrams [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF]. Our work enables faithful interpolations of merge trees (Fig. 10) and the generation of merge trees representative of a set (Figs. 1, 12, and 13). Moreover, our task-based algorithm enables automatic barycenter computations within minutes for real-life ensembles.

A natural direction for future work is the extension of our framework to other topological data representations, such as Reeb graphs or Morse-Smale complexes. However, the question of defining relevant and computable metrics for these objects is still an active research debate. Moreover, as illustrated by this paper, extending existing metrics to make them conducive to efficient geodesic computation further requires additional efforts. We believe our work is an important practical step towards the definition of a larger statistical framework on the space of merge trees. In the future, based on our framework, we will study the definition of more sophisticated statistical indexes (for instance by investigating a notion of covariance matrix for merge trees), to support even more advanced analyses of large-scale ensemble data.

APPENDIX 1 THE EDIT DISTANCE BETWEEN MERGE TREES [106]

This section formalizes the edit distance introduced by Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] (Sec. 1.1) and discusses some of its technical aspects which make it not conducive to interpolationbased geodesics (Sec. 1.2).

Definition

The edit distance between two merge trees T ( f i ) and T ( f j ), noted D E T ( f i ), T ( f j ) , is defined as follows. Let N i be a subset of the nodes of T ( f i ) and N i its complement. Let φ be a partial assignment between N i and a subset N j of the nodes of T ( f j ) (with complement N j ). Then D E T ( f i ), T ( f j ) is given by:

D E T ( f i ), T ( f j ) = min (φ ,N i ,N j )∈Φ ∑ n i ∈N i γ n i → φ (n i ) (15)
+ ∑ n i ∈N i γ(n i → / 0) (16)
+ ∑ n j ∈N j γ( / 0 → n j ) ( 17
)
where Φ is the space of constrained partial assignments (i.e. φ maps disjoint subtrees of T ( f i ) to disjoint subtrees of T ( f j )) and where γ refers to the cost for: (i) mapping a node

n i ∈ T ( f i ) to a node φ (n i ) = n j ∈ T ( f j ) (line 15), (ii) deleting a node n i ∈ T ( f i ) (line 16
) and (iii) creating a node n j ∈ T ( f j ) (line 17), / 0 being the empty tree. Zhang [121] introduced a polynomial time algorithm for computing a constrained sequence of edit operations with minimal edit distance (Eq. 15), and showed that the resulting distance is indeed a metric if each cost γ for the above three edit operations is itself a metric (non-negativity, identity, symmetry, triangle inequality). Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF] exploited this property to introduce their metric, by defining the following distance-based cost model, where p i and p j stand for the persistence pairs containing the nodes n i ∈ T ( f i ) and n j ∈ T ( f j ):

γ(n i → n j ) = min d ∞ (p i , p j ), γ(n i → / 0) + γ( / 0 → n j ) γ(n i → / 0) = d ∞ p i , ∆(p i ) γ( / 0 → n j ) = d ∞ ∆(p j ), p j .
In our work, we introduce an alternative edit distance which further adheres to the L 2 -Wasserstein distance between persistence diagrams.

Interpolation

As shown in the main manuscript (Fig. 4), the linear interpolation of D E 's matchings does not describe a shortest path (i.e. it generates inaccurate midpoints). A key technical reason for this is that D E involves assignments between nodes (of the input merge trees) and not persistence pairs. This has several consequences. First, given two input trees T ( f i ) and T ( f j ), D E 's matchings may assign a saddle node in T ( f i ) to an extremum node in T ( f j ), resulting in inconsistent interpolations in the data (from a valley to a peak). Second, D E 's matchings can possibly assign two nodes in T ( f i ) belonging to a single persistence pair of f i to nodes in T ( f j ) belonging to distinct persistence pairs in f j . This second phenomenon further challenges interpolation-based geodesics.

W T

2 IS A METRIC As further described in the main manuscript, given two merge trees T ( f i ) and T ( f j ) and their branch decomposition trees (BDTs) B( f i ) and B( f j ), the dissimilarity measure W T 2 B( f i ), B( f j ) is given by:

W T 2 B( f i ), B( f j ) = min (φ ,B i ,B j )∈Φ ∑ b i ∈B i γ b i → φ (b i ) 2 + ∑ b i ∈B i γ(b i → / 0) 2 (18) + ∑ b j ∈B j γ( / 0 → b j ) 2 1/2 ( 19
)
where φ is an isomorphism between B i ⊆ B( f i ) and B j ⊆ B( f j ).

In this section, we argue that W T 2 is a metric.

W T 2 B( f i ), B( f j ) is always non-negative (the costs γ are squared). W T 2 B( f i ), B( f j
) is symmetric (destruction and creation costs are symmetric, lines 18 and 19).

W T 2 B( f i ), B( f j ) = 0 if and only if all costs γ = 0, which only happens if B( f i ) = B( f j ) (the identity is included in Φ ).
We now argue that W T 2 preserves the triangle inequality, given three trees B( f i ), B( f j ) and B( f k ). For this, we follow a classical approach which we detail here for the sake of completeness. First, we argue that a composition of (optimal) partial rooted isomorphisms (from B( f i ) to B( f j ), then from B( f j ) to B( f k )) is itself a valid partial rooted isomorphism (and hence belong to our solution space Φ ) and that its associated cost consequently bounds by above [START_REF] Bock | TopoAngler: Interactive Topology-Based Extraction of Fishes[END_REF]. Second, we argue that this associated cost is itself bounded by above by

W T 2 B( f i ), B( f k ) (Eq.
W T 2 B( f i ), B( f j ) +W T 2 B( f j ), B( f k ) . Let (φ , B i , B j
) be the optimal solution of the partial assignment problem between B( f i ) and B( f j ). φ is a rooted isomorphism (i.e. an isomorphism between rooted subtrees) between a subtree B i of B( f i ) and a subtree B j of B( f j ) (blue, Fig. 16). Equivalently, φ can also be interpreted as a bijection between the arcs of B i and those of B j .

Let (ψ , B j , B k ) be the optimal solution of the partial assignment problem between B( f j ) and B( f k ). ψ is a rooted isomorphism between a subtree B j of B( f j ) and a subtree B k of B( f k ) (green, Fig. 16).

Let B j be the set of nodes of B( f j ) involved in both φ and ψ (B j = B j ∩ B j , in red in Fig. 16, center). Let B i be their pre-image by φ (B i = φ -1 (B j ), in red in Fig. 16, left) and B k their image by ψ (B k = ψ (B j ), in red in Fig. 16, right).

Since both φ and ψ are rooted isomorphisms, their composition ψ • φ is also a (rooted) isomorphism between the subtrees B i of B( f i ) and B k of B( f k ) (equivalently, it is a bijection between the arcs of B i and the arcs of B k ). Then (ψ • φ , B i , B j ) is itself a rooted partial isomorphism and belongs to Φ .

Then, it follows that:

W T 2 B( f i ), B( f k ) ≤ ∑ b∈B i γ b → ψ • φ (b) 2 (20)
+ ∑ b∈B i γ(b → / 0) 2 + ∑ b∈B k γ( / 0 → b) 2 1/2 .
Now, let U, V , W be scalar functions on the nodes of the set 16) such that:

B ik = B i ∪ B i ∪ (B k \ B k ) ∪ B k (cyan subset, Fig.
U(b) =                γ b → ψ • φ (b) for b ∈ B i = φ -1 (B j ) γ(b → / 0) for b ∈ B i \ B i γ(b → / 0) for b ∈ B i γ( / 0 → b) for b ∈ B k γ( / 0 → b) for b ∈ B k \ B k = ψ (B j ∩ B j ). (21)
U describes all the possible individual costs involved in the composition ψ • φ . In particular, we can re-write Eq. 20 as:

W T 2 B( f i ), B( f k ) ≤ U 2 = ∑ b∈B ik U(b) 2 1/2 . ( 22
)
V (b) =                γ b → φ (b) for b ∈ B i = φ -1 (B j ) γ b → φ (b) for b ∈ B i \ B i γ(b → / 0) for b ∈ B i 0 for b ∈ B k γ ψ -1 (b) → / 0 for b ∈ B k \ B k = ψ (B j ∩ B j ) (23)
V describes a subset of the individual costs involved in the optimal rooted partial isomorphism φ . In particular, only the costs involving B j ∩ B j (orange square, Fig. 16, middle) are excluded. Thus, we have:

W T 2 B( f i ), B( f j ) ≥ V 2 = ∑ b∈B ik V (b) 2 1/2 . ( 24
)
W (b) =                γ φ (b) → ψ • φ (b) for b ∈ B i = φ -1 (B j ) γ φ (b) → / 0 for b ∈ B i \ B i 0 for b ∈ B i γ(b → / 0) for b ∈ B k γ ψ -1 (b) → b for b ∈ B k \ B k = ψ (B j ∩ B j ) (25)
Similarly to V , W describes a subset of the individual costs involved in the optimal rooted partial isomorphism ψ . In particular, only the costs involving B k (red square, Fig. 16, right) are excluded. Thus:

W T 2 B( f j ), B( f k ) ≥ W 2 = ∑ b∈B ik W (b) 2 1/2 . ( 26
)
Now, since γ is defined by the Euclidean distance (Equation 9 of the main manuscript), we have for each node b ∈ B ik :

0 ≤ U(b) ≤ V (b) +W (b).
This can be verified by comparing the i th line of Eq. 21 to the sum of the i th lines of Eq. 23 and Eq. 25. Then, we have:

U 2 ≤ V +W 2 . (27)
Now, since the L 2 norm between vectors respects itself the triangle inequality, we have the following inequality:

V +W 2 ≤ V 2 + W 2 . (28)
Then, by combining equations 22, 27, 28, 24, and 26, it follows that:

W T 2 B( f i ), B( f k ) ≤ U 2 ≤ V +W 2 ≤ V 2 + W 2 ≤ W T 2 B( f i ), B( f j ) +W T 2 B( f j ), B( f k )
which concludes the proof.

COMPARISON TO THE EDIT DISTANCE ALGORITHM [121]

In addition to considering squared costs in our edit distance (equations of the section 3.3 of the main manuscript), our algorithm for the exploration of the search space indeed simplifies the approach by Zhang [121] (used by Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF]), as our search space is significantly more constrained. First, since our solution space only considers partial isomorphisms between rooted subtrees, this implies that the destruction of a node (a branch) b j ∈ B( f j ) necessarily implies the destruction of its subtrees, i.e. of its forest F ( f j , b j ). Thus, the admissible solutions in [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF]121] consisting in deleting b j and mapping a subtree B( f i , b i ) to one of the subtrees of b j in the forest F ( f j , b j ) are no longer admissible given our overall solution space Φ . The removal of such solutions drastically simplifies the evaluation of the distance between subtrees (being the minimum of three solutions in [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF], Eq. 12) to the Equation 11 of our main manuscript (containing only one expression to evaluate).

Second, our solution space (rooted partial isomorphisms) also implies that the nodes of B( f i ) can only be assigned to nodes with the same depth in B( f j ). This further implies that the distance between subtrees (Equation 11 of the main manuscript) only needs to be evaluated for subtrees rooted at nodes of identical depth (see Fig. 6 of the main manuscript).

Together, these two simplifications ((i) simpler subtree distance and (ii) distance evaluation restricted to subtrees of identical depth from the root) are the key adaptations of Zhang's algorithm [121] that are required for the exploration of our (more constrained) solution space.

PARALLEL COMPUTATION OF W T 2

In our work, we express the computation of W T 2 in terms of tasks, to leverage task-based shared memory parallelism. First, the Equation 10 of the main manuscript is evaluated. For this, we initiate a task at each leaf of B( f i ). If a task is the last one to compute among all the direct children of a node b ∈ B( f i ), it is then authorized to continue and estimate Equation 10 in b. Atomic counters in b are implemented (and atomically incremented by the task of each child) to determine which child task is the last one to complete, which enables an efficient lightweight synchronization (Fig. 6 of the main manuscript). Overall, Equation 10 is completely estimated with this strategy in a bottom-up fashion. Second, Equation 11(main manuscript) is evaluated similarly, by initiating a task at each leaf b i of B( f i ). In particular, this task will evaluate Equation 11given b i against all subtrees of B( f j ) of identical depth (again using independent tasks initiated at the leaves of B( f j ), see Fig. 6). Similarly to Equation 10, we employ the same lightweight synchronization mechanism based on atomic counters to continue a task over to its parent only when it is the last child task reaching it. Thus, in both cases (Eqs. 10 and 11), the number of parallel tasks is initially bounded by the number of leaves in B( f i ) and B( f j ) (which is typically much larger than the number of cores) and progressively decreases during the computation.

B

α ( f i → f j ) IS A GEODESIC FOR 0 ≤ α ≤ 1
We now argue that W T 2 defines a geodesic space. For this, for any two BDTs B( f i ) and B( f j ), we describe the existence of a path between them whose length is equal to W T 2 B( f i ), B( f j ) (and thus minimal). Let P = B t t∈[0,1] be a path of BDTs parameterized by t. We recall that the length L (P) of P is given by:

L (P) = sup n;0=t 0 ≤t 1 ≤•••≤tn =1 n-1 ∑ k=0 W T 2 B t k , B t k +1 .
Now, let P α be the path corresponding to the interpolation between B( f i ) and B( f j ), as defined in section 4.1 of the main manuscript. We now argue that L (P α ) = W T 2 B( f i ), B( f j ) . Let (φ , B i , B j ) be the optimal rooted partial isomorphism between B( f i ) and B( f j ). Moreover, let B s ( f i → f j ) and B t ( f i → f j ) be two interpolated trees obtained respectively with α = s and α = t, given 0 ≤ s ≤ t ≤ 1. We will note φ s the application on B i ∪ B i ∪ B j defined by interpolation (section 4.1 of the main manuscript):

φ s (b) =      (1 -s)b + sφ (b) for b ∈ B i (1 -s)b + s∆(b) for b ∈ B i sb + (1 -s)∆(b) for b ∈ B j . ( 29
)
φ t is defined similarly for t. Then, as discussed in Sec. 2 of this appendix, since the composition of partial rooted isomorphisms is itself a partial rooted isomorphism, the composition φ t • φ s -1 (which goes from B s ( f i → f j ) to B( f i ) and then from B( f i ) to

B t ( f i → f j ))
does define a valid partial rooted isomorphism between B s ( f i → f j ) and B t ( f i → f j ) and we have:

W T 2 B s ( f i → f j ), B t ( f i → f j ) ≤ ∑ b∈B i ∪B i ∪B j γ φ s (b) → φ t (b) 2 1/2
and we also have by definition of φ s and φ t (Eq. 29):

∑ b∈B i ∪B i ∪B j γ φ s (b) → φ t (b) 2 1/2 = (t -s)W T 2 B( f i ), B( f j ) .
Now, given the triangle inequality on the path P α , we have:

W T 2 B( f i ), B( f j ) ≤ W T 2 B 0 ( f i → f j ), B s ( f i → f j ) + W T 2 B s ( f i → f j ), B t ( f i → f j ) + W T 2 B t ( f i → f j ), B 1 ( f i → f j ) ≤ s + (t -s) + (1 -t) W T 2 B( f i ), B( f j ) = W T 2 B( f i ), B( f j ) .
If follows that the above inequalities are in fact equalities and we have:

W T 2 B s ( f i → f j ), B t ( f i → f j ) = (t -s)W T 2 B( f i ), B( f j ) .
Then, for any subdivision 0 = t 0 ≤ t 1 ≤ • • • ≤ t n = 1 of P α , we have:

n-1 ∑ k=0 W T 2 B t k ( f i → f j ), B t k +1 ( f i → f j ) = n-1 ∑ k=0 (t k+1 -t k )W T 2 B( f i ), B( f j ) = (t n -t 0 )W T 2 B( f i ), B( f j ) = W T 2 B( f i ), B( f j ) .
Thus L (P α ) = W T 2 B( f i ), B( f j ) . Hence the space of merge trees equipped with W T 2 is a geodesic space, and B α ( f i → f j ) constructs paths of minimal length on it.

Asteroid impact

This ensemble is composed of 7 members, given as 3D regular grids (sampled at 300 × 300 × 300, implicitly triangulated by TTK). It has been made available in the context of the SciVis contest 2018 [START_REF] Patchett | The IEEE SciVis Contest[END_REF]. Each member corresponds to the last time step of the simulation of the impact of an asteroid with the sea at the surface of the Earth, for two configurations of asteroid diameter. The considered scalar field is the matter density, which is one of the variables of the simulation which discriminates well the asteroid from the water and the ambient air. This ensemble corresponds to a parameter study (in this case, studying the effect of the asteroid's diameter on the resulting wave), which is a typical task in numerical simulation. In this application, salient maxima capture well the asteroid and large water splashes. Thus, each member is represented by the split tree (capturing maxima). The associated ground-truth classification assigns members computed with similar asteroid diameters to the same class. Thus, the corresponding classification task consists in identifying, for a given member, its correct asteroid diameter class. The ground-truth classification is as follows:

• Class 1 (3 members): yA11, yB11, yC11.

• Class 2 (4 members): yA31, yA32, yB31, yC31

Another selection of the original data has been used for the evaluation of our temporal reduction framework (Fig. 10 of the main manuscript). For this experiment, we used the asteroid diameter configuration "yA31" and considered the following time steps, organized in 4 phases (according to the SciVis contest companion documentation [START_REF] Patchett | The IEEE SciVis Contest[END_REF]):

• Phase 1, initial state (5 time steps): 01141, 03429, 05700, 07920, 09782

Cloud processes

This ensemble is composed of 12 members, given as 2D regular grids (sampled at 1430 × 1557, implicitly triangulated by TTK). Each member corresponds to a time step of the simulation of cloud formations [START_REF] Wischgoll | The IEEE SciVis Contest[END_REF]. For this application, large clouds are well captured by the maxima of the pressure variable (pre-processed with 10 iterations of smoothing). Thus, split trees (capturing maxima) are considered for this ensemble. The associated ground-truth classification assigns each time step to one of the three key phases of the simulation. The corresponding classification task therefore consists in identifying, for each time step, to which phase it belongs. The ground-truth classification is as follows:

• Class 1 (4 members): 0, 5, 10, 15

Viscous fingering

This ensemble is composed of 15 members, given as 3-dimensional point clouds (representing a particle-based flow simulation). Each point cloud is turned into a Eulerian representation of the variables by using the "Gaussian Resampling" filter of ParaView, effectively transforming, via interpolation [START_REF] Shepard | A two-dimensioanl interpolation function for irregularlyspaced data[END_REF], each ensemble member into a 3D regular grid (sampled at 50 × 50 × 50, implicitly triangulated by TTK). The original data has been made available in the context of the SciVis contest 2016 [START_REF] Garth | The IEEE SciVis Contest[END_REF]. Each member corresponds to the last time step of the simulation of a viscous fingering phenomenon, occurring when dissolving salt in water. The considered scalar field is the salt concentration, whose salient maxima capture well the most prominent fingers. Thus, each member is represented by the split tree (capturing maxima). Given the studied physical phenomenon, the simulation approach is not deterministic, resulting in distinct outputs for identical initial configurations. In this application, three distinct solver resolutions have been considered, corresponding to three distinct numbers of particles (resolution code 20: 194k particles, resolution code 30: 544k particles, resolution code 44: 1.7M particles). Thus, this ensemble corresponds to a parameter study (in this case, studying the effect of the input resolution on the output fingering), which is a typical task in numerical simulation. The associated ground-truth classification assigns members with the same input resolution to the same class. Thus, the corresponding classification task consists in identifying, for a given ensemble member, its corresponding particle count. The ground-truth classification is as follows:

• Class 1, resolution 20 (5 members): 20run1, 20run3, 20run4, 20run5, 20run6.

• Class 2, resolution 30 (5 members): 30run1, 30run2, 30run3, 30run4, 30run5

• Class 3, resolution 44 (5 members): 30run1, 30run2, 30run3, 30run4, 30run5

Dark matter

This ensemble is composed of 40 members, given as 3-dimensional point clouds (representing a particle-based simulation). Each point cloud is turned into a Eulerian representation of the variables by using the "Gaussian Resampling" filter of ParaView, effectively transforming, via interpolation [START_REF] Shepard | A two-dimensioanl interpolation function for irregularlyspaced data[END_REF], each ensemble member into a 3D regular grid (sampled at 100 × 100 × 100, implicitly triangulated by TTK). The original data has been made available in the context of the SciVis contest 2015 [START_REF] Hentschel | The IEEE SciVis Contest[END_REF]. Each member corresponds to a time step of a simulation of the universe formation, where regions of high concentration of dark matter form a filament structure known as the cosmic web. The considered scalar field is therefore dark matter density, whose salient maxima capture well large clusters of galaxies. Thus, each member is represented by the split tree (capturing maxima). The associated ground-truth classification assigns each time step to one of the four key phases of the simulation. The corresponding classification task therefore consists in identifying, for each time step, to which phase it belongs. The ground-truth classification is as follows:

• Class 1 (10 members): 0.0200, 0.0300, 0.0400, 0.0500, 0.0600, 0.0700, 0.0800, 0.0900, 0.1000, 0.1100

• Class 2 (10 members): 0.2700, 0.2800, 0.2900, 0.3000, 0.3100, 0.3200, 0.3300, 0.3400, 0.3500, 0.3600

• Class 3 (10 members): 0.5900, 0.6000, 0.6100, 0.6200, 0.6300, 0.6400, 0.6500, 0.6600, 0.6700, 0.6800

• Class 4 (10 members): 0.9100, 0.9200, 0.9300, 0.9400, 0.9500, 0.9600, 0.9700, 0.9800, 0.9900, 1.0000

Volcanic eruptions

This ensemble is composed of 12 members, given as 2D regular grids (sampled at 500 × 500, implicitly triangulated by TTK). Each member corresponds to an observation of a volcanic eruption, obtained by satellite imaging (as this data exhibits a bit of noise, it has been pre-simplified by removing all saddle-maxima pairs with a persistence lower than 0.5% of the data range). The original data has been made available in the context of the SciVis contest 2014 [START_REF] Feng | The IEEE SciVis Contest[END_REF]. The considered scalar field is the sulfur dioxide concentration, for which salient maxima correspond to volcanic eruptions. Thus, each observation is represented by the split tree (capturing maxima). Each member corresponds to a specific acquisition period, itself corresponding to the eruption of one particular volcano at the surface of the Earth. The associated ground-truth classification assigns observations acquired in the same period of time to the same class. The corresponding classification task therefore consists in identifying, for each observation (taken at a specified date), the erupting volcano it corresponds to. The ground-truth classification is as follows:

•

Ionization front (3D)

This ensemble is composed of 16 members, given as 3D regular grids (sampled at 300 × 124 × 124, implicitly triangulated by TTK). Each member corresponds to a time step of a simulation of ionization front propagation [START_REF] Taylor | The IEEE SciVis Contest[END_REF]. For this application, large ionization flares are well captured by salient maxima of the ion concentration. Thus, split trees (capturing maxima) are considered for this ensemble. The associated ground-truth classification assigns each time step to one of the four key phases of the simulation. The corresponding classification task therefore consists in identifying, for each time step, to which phase it belongs. The ground-truth classification is as follows:

• Class 1 (4 members): 0025, 0026, 0027, 0028

• Class 2 (4 members): 0075, 0076, 0077, 0078,

• Class 3 (4 members): 0125, 0126, 0127, 0128

• Class 4 (4 members): 0175, 0176, 0177, 0178,

Ionization front (2D)

This ensemble is a 2D version of the above ensemble, where the dataset providers have selected a 2D slice in the center of the volume (sampled at 600 × 248). The associated classification task is therefore identical.

Earthquake

This ensemble is composed of 12 members, given as 3D regular grids (sampled at 375 × 188 × 50, implicitly triangulated by TTK). Each member corresponds to a time step of the simulation of an earthquake at the San Andreas fault [START_REF] Olsen | The IEEE SciVis Contest[END_REF]. For this application, the shock wave can be tracked with the local maxima of the wave front velocity magnitude (this scalar field is pre-processed to pre-simplify all saddle-maxima pairs with a persistence smaller than 0.05% of the data range). Thus, split trees (capturing maxima) are considered for this ensemble. The associated ground-truth classification assigns each time step to one of the three key phases of the simulation. The corresponding classification task therefore consists in identifying, for each time step, to which phase it belongs. The ground-truth classification is as follows:

• Class 1 (4 members): 002700, 002900, 003100, 003300

• Class 2 (4 members): 007700, 007900, 008100, 008300

• Class 3 (4 members): 011700, 011900, 012100, 012300

Isabel

This ensemble is composed of 12 members, given as 3D regular grids (sampled at 250 × 250 × 50, implicitly triangulated by TTK). Each member corresponds to a time step of the simulation of the Isabel hurricane [START_REF] Wang | The IEEE SciVis Contest[END_REF]. This ensemble has been used in previous work [START_REF] Favelier | Persistence Atlas for Critical Point Variability in Ensembles[END_REF][START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] and the corresponding VTK files are available at the following address: https: //github.com/julesvidal/wasserstein-pd-barycenter. In this application, the eyewall of the hurricane is typically characterized by high wind velocities, well captured by the the maxima of the flow velocity. Thus, split trees (capturing maxima) are considered for this ensemble. The associated ground-truth classification assigns each time step to one of the three key phases (formation, drift, landfall) of the hurricane simulation. The corresponding classification task therefore consists in identifying, for each member, to which key phase it belongs. The ground-truth classification is as follows:

• Class 1 (4 members): 2, 3, 4, 5

Starting vortex

This ensemble is composed of 12 members, given as 2D regular grids (sampled at 1500 × 1000, implicitly triangulated by TTK). It has been generated with the Gerris flow solver [START_REF] Popinet | Gerris Flow Solver[END_REF] and was provided in previous work [START_REF] Favelier | Persistence Atlas for Critical Point Variability in Ensembles[END_REF][START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF]. It is available at the following address: https: //github.com/julesvidal/wasserstein-pd-barycenter. The data models flow turbulence behind a wing, for two ranges of wing inclination angles. The considered scalar field is the orthogonal component of the curl of the flow velocity. This ensemble corresponds to a parameter study (in this case, studying the effect of the wing configuration on turbulence), which is a typical task in numerical simulation. In this application, salient extrema are typically considered as reliable estimations of the center of vortices. Thus, each run is represented by two merge trees (the join tree -capturing minima, and the split tree, capturing maxima), which are processed independently by our algorithms. The associated ground-truth classification assigns members computed with similar inclination angles to the same class. The corresponding classification task therefore consists in identifying, for a given ensemble member, its correct wing configuration class. The ground-truth classification is as follows:

• Class 1 (6 members): Angle=2, Angle=3, Angle=4, Angle=5, Angle=6, Angle=8

• Class 2 (6 members): Angle=38, Angle=39, Angle=40, Angle=41, Angle=42, Angle=43

Sea surface height

This ensemble is composed of 48 members, given as 2D regular grids (sampled at 1440 × 720, implicitly triangulated by TTK). Each member corresponds to an observation of the sea surface height at the surface of the Earth, taken in January, April, July and October 2012. The original data can be found at the following address: https://ecco.jpl.nasa. gov/products/all/. This ensemble has been used in previous work [START_REF] Favelier | Persistence Atlas for Critical Point Variability in Ensembles[END_REF][START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF] and the corresponding VTK files are available at the following address: https://github.com/ julesvidal/wasserstein-pd-barycenter. In this application, the features of interest are the center of eddies, which can be reliably estimated with height extrema. Thus, each observation is represented by two merge trees (the join tree -capturing minima, and the split tree, capturing maxima), which are processed independently by our algorithms. The associated ground-truth classification assigns observations acquired in the same month to the same class. The corresponding classification task therefore consists in identifying, for each observation (taken at a specified date), the season in which it has been acquired. The ground-truth classification is as follows:

•

Vortex street

This ensemble is composed of 45 members, given as 2D regular grids (sampled at 300 × 100, implicitly triangulated by TTK). It has been generated with the Gerris flow solver [START_REF] Popinet | Gerris Flow Solver[END_REF] and was provided in previous work [START_REF] Favelier | Persistence Atlas for Critical Point Variability in Ensembles[END_REF][START_REF] Vidal | Progressive Wasserstein Barycenters of Persistence Diagrams[END_REF]. It is available at the following address: https://github.com/julesvidal/wasserstein-pd-barycenter. The data models flow turbulence behind an obstacle. The considered scalar field is the orthogonal component of the curl of the flow velocity, for 5 fluids of different viscosity. This ensemble corresponds to a parameter study (in this case, studying the effect of viscosity on turbulence), which is a typical task in numerical simulation. In this application, salient extrema are typically considered as reliable estimations of the center of vortices. Thus, each run is represented by two merge trees (the join tree -capturing minima, and the split tree, capturing maxima), which are processed independently by our algorithms. The associated ground-truth classification assigns members computed with similar viscosities to the same class. The corresponding classification task therefore consists in identifying, for a given ensemble member, its correct viscosity class. The ground-truth classification is as follows:

• Class 1 (9 members): Viscosity=100.0, Viscosity=100.1, Viscosity=100.2, Viscos-ity=100.3, Viscosity=100.4, Viscosity=100.5, Viscosity=100.6, Viscosity=100.7, Viscosity=100.9

• Class 2 (9 members): Viscosity=160.0, Viscosity=160.1, Viscosity=160.2, Viscos-ity=160.3, Viscosity=160.4, Viscosity=160.5, Viscosity=160.6, Viscosity=160.7, Viscosity=160.8

• Class 3 (9 members): Viscosity=200.0, Viscosity=200.1, Viscosity=200.2, Viscos-ity=200.3, Viscosity=200.4, Viscosity=200.5, Viscosity=200.6, Viscosity=200.7, Viscosity=200.8

• Class 4 (9 members): Viscosity=50.0, Viscosity=50.1, Viscosity=50.2, Viscos-ity=50.3, Viscosity=50.5, Viscosity=50.6, Viscosity=50.7, Viscosity=50.8, Viscos-ity=50.9

• Class 5 (9 members): Viscosity=60.1, Viscosity=60.2, Viscosity=60.3, Viscos-ity=60.4, Viscosity=60.5, Viscosity=60.6, Viscosity=60.7, Viscosity=60.8, Viscos-ity=60.9

PARAMETER ANALYSIS

In this section, we study the practical effect of the parameters of our approach. In particular, we extend the empirical stability evaluation of our metric with regard to all the parameters of our approach and we illustrate their effect on geodesic computation.

Fig. 2 .

 2 Fig. 2. Critical points (spheres, white: minima, blue: maxima, other: saddles), persistence diagram (bottom left), merge tree (bottom right) and branch decomposition tree (top right) of a clean (a) and noisy (b) 2D scalar field. In both cases, four main hills are clearly represented with salient features in the persistence diagram and the merge tree. Branches with low persistence (less than 10% of the function range) are shown with small white arcs. They correspond to noisy features in the data (b).

Fig. 6 .

 6 Fig. 6. The exploration of the space Φ (a candidate is highlighted in red, right) relies on the evaluation of the sparse matrix T of subtree distances (left). Our task-based algorithm optimizes the parallel computation of independent terms. Spheres of equal radius in T denote independent terms and arrows between the lines of T indicate task dependence (equivalently illustrated with arrows in the input BDTs, right).

Fig. 7 .

 7 Fig. 7. Our geodesic computation extends interpolation-based geodesics from persistence diagrams (a) to merge trees (b). The interpolated BDT B α ( f i → f j ) is obtained by linear interpolation (with local normalization) of the partial isomorphism φ in the birth/death space. In the data, the feature matching (dashed lines) induced by φ with W T 2 (b) better preserves the global structure of the data than φ with W D 2 (a, red crossing).

Fig. 8 .

 8 Fig. 8. Given two scalar fields f i (a) and f j (b), a simple interpolation of the birth/death values of the branches of their BDTs may result in inconsistencies upon branch destruction (red): the interpolated merge tree (black) in (c) is disconnected, unlike the interpolated BDT (d). Our local normalization (Sec. 4.2) addresses this issue by enforcing nested birth/death values for nested branches. This results in a valid interpolated merge tree (e) whose BDT is indeed equal to the interpolated BDT (e).

Fig. 10 .

 10 Fig. 10. Geodesic computation for the reduction (b) of temporal sequences of merge trees (a). Our algorithm greedily removes from the sequence the trees that it can accurately estimate by geodesic computation (trees with blue background (b)). This reduction is also visualized with the three curves in (c), plotting the distance W T 2 to the empty tree / 0 over time (multiple colors: original sequence, cyan: reduction by W T 2 , red: reduction by W D2 ). This iterated removal of trees highlights key frames in the sequence (d-g) corresponding to key phases of an asteroid impact simulation[START_REF] Patchett | The IEEE SciVis Contest[END_REF]: initial state (black, time steps 1-5), approach (light green, 6-10), impact (blue, 11-15), aftermath (light blue,[START_REF] Bhatia | Topoms: Comprehensive topological exploration for molecular and condensed-matter systems[END_REF][START_REF] Bhatia | Flow visualization with quantified spatial and temporal errors using edge maps[END_REF][START_REF] Biasotti | Reeb graphs for shape analysis and applications[END_REF][START_REF] Bin Masood | An Overview of the Topology ToolKit[END_REF][START_REF] Bock | TopoAngler: Interactive Topology-Based Extraction of Fishes[END_REF]. In contrast, a similar greedy optimization based on the distance W D 2 between persistence diagrams (red curve) fails at capturing the impact phase (blue) of this sequence.

Fig. 11 .

 11 Fig. 11. Tracking features (the five most persistent maxima, spheres) in time-varying 2D data (ion density during universe formation [109]): optimal assignment φ (a) of W D 2 (Sec. 2.2), optimal isomorphism φ (b) of W T 2 (Sec. 3.2). Since W D 2 considers persistence pairs individually, it can generate incorrect matchings resulting in a characteristic crossing (a, red). Our distance improves this aspect (b) thanks to its more constrained search space, which better preserves the global structure of the data.

Fig. 14 .

 14 Fig. 14. Empirical stability evaluation. Given an input scalar field f i , a noisy version f j is created by inserting a random noise of increasing amplitude ε (left). The evolution of W T 2 B( f i ), B( f j ) with ε (right), for varying values of ε 1 (Sec. 3.2), indicates clear transition points (colored dots) before which W T2 evolves nearly linearly. Before these transition points (i.e. before these noise levels, vertical lines), W T 2 is stable.

Fig. 15 .

 15 Fig. 15. Evolution of the Fr échet energy for estimating the barycenter (b) for an ensemble of 100 noisy variants of four fields (a). The energy (c) is shown for several ε 1 values (Sec. 3.2). In practice, we stop the algorithm when the energy decreases by less than 1% (vertical lines).

Fig. 16 .

 16 Fig.[START_REF] Bhatia | Topoms: Comprehensive topological exploration for molecular and condensed-matter systems[END_REF]. The composition ψ • φ of two rooted isomorphisms φ (blue) and ψ (green) is itself a rooted isomorphism (red). In this schematic view, the involved subtrees are represented as squares.

Table 1 .

 1 Running times (in seconds, 10 run average) of our approach for the barycenter computation, with respect to W D 2 (ε 1 = 1, Sec. 3.3, sequential) and to our new metric W T 2 (sequential, then with 20 cores).

	Dataset	N	|B|	W D 2 (1 c.)	W T 2 (1 c.) W T 2 (20 c.) Speedup
	Asteroid Impact [84] (3D)	7 1,295	514.71	450.91	93.11	4.84
	Cloud processes [116] (2D)	12 1,209	54.90	124.99	35.14	3.55
	Viscous fingering [45] (3D)	15	118	5.68	5.12	3.89	1.31
	Dark matter [56] (3D)	40 2,592	3,172.37	3,083.24	471.45	6.53
	Volcanic eruptions [41] (2D)	12	811	171.13	140.02	48.52	2.88
	Ionization front [109] (2D)	16	135	10.40	12.10	8.20	1.47
	Ionization front [109] (3D)	16	763	682.76	1,277.72	219.61	5.81
	Earthquake [77] (3D)	12 1,203	191.54	509.59	117.31	4.34
	Isabel [114] (3D)	12 1,338	330.88	284.19	62.70	4.53
	Starting Vortex [39] (2D)	12	124	7.72	5.58	6.11	0.91
	Sea Surface Height [39] (2D)	48 1,787	4,509.78	10,557.07	881.49	11.97
	Vortex Street [39] (2D)	45	23	1.71	1.90	1.44	1.31

•

  Class 2 (4 members): 30, 31, 32, 33 • Class 3 (4 members): 45, 46, 47, 48

  Class 1 (12 members): 20120111, 20120115, 20120116, 20120117, 20120118, 20120119, 20120120, 20120121, 20120123, 20120128, 20120129 • Class 2 (12 members): 20120419, 20120420, 20120421, 20120422, 20120423, 20120424, 20120425, 20120426, 20120427, 20120428, 20120429, 20120430 • Class 3 (12 members): 20120711, 20120712, 20120713, 20120714, 20120715, 20120716, 20120717, 20120718, 20120719, 20120720, 20120721, 20120722 • Class 4 (12 members): 20121008, 20121009, 20121010, 20121011, 20121012, 20121016, 20121017, 20121018, 20121019, 20121020, 20121022, 20121023

ACKNOWLEDGMENTS

This work is partially supported by the European Commission grant ERC-2019-COG "TORI" (ref. 863464, https://erc-tori.github.io/).

WITH NORMALIZED COSTS DEFINES A GEODESIC SPACE

Let W N 2 B( f i ), B( f j ) be a similarity measure between B( f i ) and B( f j ), defined as:

where N is the local normalization described in Section 4.2 of the main manuscript. Since N is invertible, W N 2 inherits all the properties of W T 2 and is also a distance metric. The normalized interpolation B s ( f i → f j ), s ∈ [0, 1] between B( f i ) and B( f j ) is defined as the image by N -1 of the interpolation between the normalized trees N (B( f i )) and N (B( f j )). Then, given s and t such that 0 ≤ s ≤ t ≤ 1, it follows that:

which proves that the space of merge trees equipped with W N 2 is a geodesic space, and that the above normalized interpolation constructs paths of minimal length on it.

MINIMIZING THE FR ÉCHET ENERGY

The optimization algorithm described in Section 5.2 of the main manuscript constructively decreases the Fréchet energy at each iteration. In particular, once a local minimizer of the Fréchet energy is obtained for a fixed assignment with the update step (ii), the subsequent assignment step (i) does further improve the assignments hence iteratively decreasing the Fréchet energy constructively.

Let F be a function of an arbitrary BDT B and of an arbitrary (i.e. not necessarily optimal) set of N rooted partial isomorphisms (φ i , B i , B B ) i=1,...,N between B and the N BDTs of S B :

. Now, let B k be the candidate barycenter at the iteration k of the algorithm and let

be the optimal rooted partial isomorphism between B k and B( f i ), computed by the assignment step of the iteration. Then, we have:

Next, the update step of the iteration k consists in moving B k to B k+1 by placing (in the 2D birth/death space) each branch b ∈ B k at the arithmetic mean of the assignments. Since the arithmetic mean generally minimizes sums of Euclidean distances, we have:

Now, observe that since the previous rooted partial isomorphisms are not optimal anymore for B k+1 , we also have:

Once B k+1 is fixed, all the rooted partial isomorphisms are then optimized again with the assignment step of the iteration k + 1 to attain:

The result of these two steps is that:

Then, each iteration of our algorithm indeed decreases the Fréchet energy. Since there is a finite number of combinations of rooted partial isomorphisms between the barycenter and the N input trees B( f i ), it follows that the algorithm converges, in a finite number of steps, to a local minimum B * of the Fréchet energy (if multiple, equally valued, optimal sets of assignments exist between B * and S B , each one needs to be explored with the update step of our algorithm). In practice, as described in the manuscript, we stop our algorithm when the Fréchet energy has decreased by less than 1% between consecutive iterations. (red, bottom). By construction, the reduction algorithm identifies as key frames the first and last time steps, irrespective of the employed metric.

TEMPORAL REDUCTION ALGORITHM

Let S = {B( f 1 ), B( f 2 ), . . . , B( f N )} be the input temporal sequence of BDTs (we assume a regular temporal sampling). Let K ⊆ S be a set of key frames. Let S = {B ( f 1 ), B ( f 2 ), . . . , B ( f N )} be a reduced temporal sequence, where:

where B( f j ) and B( f k ) are two consecutive trees in K , such that j ≤ i ≤ k and α i = (ij)/(kj). B ( f i ) is then on a geodesic between B( f j ) and B( f k ). We introduce the following distance between the temporal sequences S and S :

d S is indeed a metric since it is a composition of metrics (being the L 2 norm between vectors of BDTs under the metric W T 2 ). Our algorithm for temporal reduction consists in initializing K with the entire input sequence (K ← S ) and then removing greedily, at each iteration, the tree B * from K (K ← K -{B * }) which minimizes d S (S , S ), and which, hence, better preserves the input sequence, until K reaches a target size.

Fig. 17 shows the temporal reduction performed by this algorithm on the Asteroid impact sequence (see Section 9.1). This figure illustrates key frames, which correspond to time steps for which B( f i ) = B ( f i ): these are the time steps which have not been removed from the sequence through the reduction. In particular, this figure compares the usage of two metrics in the reduction algorithm: W T 2 (blue, top) and W D 2 (red, bottom, obtained with ε 1 = 1). By construction, since our reduction algorithm is based on interpolation only, the first and last BDTs in the sequence S are always kept in the reduced sequence S . In other words, the first (leftmost, Fig. 17) and last (rightmost, Fig. 17) time steps are always identified as key frames, irrespective of the employed metric. For this specific example, the second key frame (second from left, Fig. 17) also happens to be identical for both metrics. In contrast to the sequence extremities, the common identification of this time step as a key frame by W T 2 and W D 2 is not obtained by construction: the reduction algorithm did select this key frame in both configurations. Then, only the third key frame (third from left, Fig. 17) is different in this example. In particular, when using W T 2 , the reduction algorithm identifies one key frame per key phase of the simulation (see Section 9.1, each key phase is represented in Fig. 17 with a frame of distinct color). In contrast, when using W D 2 , the third key frame belongs to the same key phase as the last key frame ("Aftermath", light blue frame). Then the reduction driven by W D 2 fails at identifying a key frame for the third key phase ("Impact", dark blue frame). This is confirmed visually in Fig. 17, as the third key frame identified with W T 2 (in blue) seems to represent an intermediate step in the simulation between the second and fourth key frames. In contrast, the third key frame identified by the reduction with W D 2 (in red) is more visually similar to the fourth key frame, and hence possibly more redundant.

DATA SPECIFICATION

This section provides a complete specification of the ensemble datasets used in the paper. In particular, we document the data provenance, its representation, its pre-processing when applicable, and we specify the associated ground-truth classification.

All of these ensemble datasets were extracted from public repositories. We additionally provide a set of scripts which automatically download all of these datasets (at the exception of Asteroid impact and Cloud processes, for which the dataset providers need to be contacted personally), pre-process them with TTK and output them in VTK file format, with the ground-truth classification attached to the files as metadata (i.e. "Field Data" in the VTK terminology). For convenience, we also provide an archive containing the entire curated ensemble datasets (in VTK file format). All of this new material (scripts and curated data) is located at the following address: https://github.com/MatPont/WassersteinMergeTreesData.

Moreover, we also provide in additional material all the ensembles of merge trees computed from these datasets (in the code archive containing the implementation of our method).

Interpretation

The first parameter of our approach is ε 1 ∈ [0, 1]. It dictates the merge of saddles in the input trees, to mitigate saddle swap instabilities, as previously documented by Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF]. Adjacent saddles in the input trees are merged if their relative difference in scalar value (relative to the largest function difference between adjacent saddles) is smaller than ε 1 . For ε 1 = 0, no saddle merge is performed whereas for ε 1 = 1, all saddles are merged and W T 2 becomes equivalent to the L 2 Wasserstein distance between persistence diagrams, noted W D 2 . The local normalization step of our framework (Section 4.2 of the main manuscript) guarantees the topological consistency of the interpolated merge trees (Fig. 8 of the main manuscript). However, this normalization shrinks the birth/death values of all the input branches to the interval [0, 1], irrespective of their original persistence. To mitigate this effect, the input BDTs are pre-processed, so that branches with small initial persistence (i.e. small branches) are not given too much importance in the metric. In particular, small branches are moved up the input BDT if their persistence relative to their parent is larger than ε 2 ∈ [0, 1]. When ε 2 = 0, all branches are moved up to the root of the BDT and again, W T 2 becomes equivalent to W D 2 . When ε 2 = 1, no branch is moved up the BDT and ε 2 has no effect on the outcome (i.e. the input BDT is left unchanged). In practice, we recommend the default value ε 2 = 0.95: if a branch b has a nearly identical persistence to that of its parent b , it is moved higher in the BDT, so that its normalized persistence becomes nearly identical to that of its parent b (instead of being artificially larger due to the local normalization).

The parameter ε 3 ∈ [0, 1] further restricts the application of the above BDT preprocessing, by only considering (for displacement up the BDT) the branches with a relative persistence (with respect to the overall data range) smaller than ε 3 . When ε 3 = 1, all branches are subject to the above pre-processing and ε 2 fully dictates the BDT preprocessing. When ε 3 = 0, no branch is moved up the BDT and the two parameters ε 2 and ε 3 have no effect on the outcome. In practice, we recommend the default value ε 3 = 0.9, which prevents the most persistent branches from moving up the BDT.

Overall, when the parameters (ε 1 , ε 2 , ε 3 ) are set to the values (0, 1, 1), the input trees are not pre-processed by the above procedures (i.e. they are left unchanged) and their structure has a strong impact on W T 2 . When ε 1 = 1 or when ε 2 = 0, W T 2 becomes equivalent to W D 2 and the structure of the input trees has no impact anymore on the metric. In-between values balance the importance of the structure of the trees on the metric. We recommend the default values (0.05, 0.95, 0.9), which provides an acceptable stability with regard to saddle swaps (mitigated by ε 1 ) and which gives a reasonable importance to small branches in the metric (controlled by ε 2 and ε 3 , which are dependent parameters). Fig. 18. Empirical stability evaluation with regard to ε 2 . Given an input scalar field f i , a noisy version f j is created by inserting a random noise of increasing amplitude ε (cf. Figure 14 of the main manuscript). The evolution of W T 2 B( f i ), B( f j ) with ε is reported for varying values of ε 2 .

Fig. 19. Empirical stability evaluation with regard to ε 3 . Given an input scalar field f i , a noisy version f j is created by inserting a random noise of increasing amplitude ε (cf. Figure 14 of the main manuscript). The evolution of W T 2 B( f i ), B( f j ) with ε is reported for varying values of ε 3 .

Metric stability

Figure 14 of the main manuscript provides an empirical stability evaluation of our new metric W T 2 , as a function of an input perturbation, modeled by a random noise of amplitude ε. In particular, this experiment is achieved for several values of ε 1 . The conclusion of this experiment is that W T 2 is not stable when ε 1 = 0 (sudden increase in W T 2 for small values of ε) and that it is stable when ε 1 = 1 (as anticipated [START_REF] Turner | Fréchet Means for Distributions of Persistence Diagrams[END_REF]). For in-between values, W T 2 is stable until a transition point (colored dots in Fig. 14 of the main manuscript), located at increasing noise levels (ε) for increasing values of ε 1 . In particular, for the recommended default value ε 1 = 0.05, W T 2 is stable up to a perturbation noise of amplitude 16% (of the overall data range).

In the following, we perform the same study for the other parameters of our approach, ε 2 and ε 3 . Fig. 18 studies the practical stability of W T 2 , for several values of ε 2 . For this experiment, ε 3 has been set to 1 (then, only ε 2 has an impact on the BDT pre-processing described in the previous section). Moreover, ε 1 has been set to its recommended value, 0.05. Several curves are reported, one per ε 2 values. For ε 2 = 0, all branches are moved up the BDT (irrespective of ε 1 ) and W T 2 becomes equivalent to W D 2 and the corresponding curve (red) exactly coincides with the light blue curve of the Figure 14 of the main manuscript (ε 1 = 1). For ε 2 = 1, the input BDT is not pre-processed at all and the corresponding curve (pink) exactly coincides with the cyan curve of Figure 14 of the main manuscript (obtained for the default value ε 1 = 0.05). In-between values of ε 2 result in continuous transitions between these two extreme cases (blue, green and cyan curves).

Fig. 19 studies the practical stability of W T 2 , for several values of ε 3 . For this experiment, we set ε 1 = 0 and ε 2 = 0, to better isolate the effect of ε 3 . When ε 3 = 1, all the branches of the input BDTs are subject to the BDT pre-processing. Since ε 2 = 0, all branches are moved up to the root and W T 2 becomes equivalent to W D 2 and the corresponding curve (red) exactly coincides with the light blue curve of the Figure 14 of the main manuscript (ε 1 = 1). When ε 3 = 0, no branch is moved up in the input BDTs and the corresponding curve (pink) exactly coincides with the grey curve of the Figure 14 of the main manuscript (ε 1 = 0). In-between values of ε 3 result in transitions between these two extreme cases (blue, green and cyan curves), with transition points (similar to the Figure 14 of the main manuscript), before which W T 2 is stable. Note however, that since it is dependent on ε 2 (default value: 0.95), ε 3 has only a very mild practical impact on the metric.

Geodesic analysis

Figures 20, 21 and 22 respectively illustrate the effect of the parameters ε 1 , ε 2 and ε 3 on the geodesics between merge trees. In particular, each figure shows, on the left, the geodesic obtained with a disabling value of the parameter (no effect on the computation). In contrast, the right side of each figure shows the geodesic obtained with the recommended default value of the parameter, to clearly visualize its impact.

Overall, as discussed in the detailed captions, these three parameters have the effect of moving branches up the input BDTs, hence reducing the structural impact of the trees on the metric, but also improving its stability (as discussed in Section 10.2). In the data, moving a branch up the BDT corresponds to only slight modifications, which consist in reconnecting maxima to distinct saddles. For each parameter, the resulting pre-processing addresses cases where nearby saddles have very close function values, which impacts the stability of the metric. Similarly to Sridharamurthy et al. [START_REF] Sridharamurthy | Edit distance between merge trees[END_REF], we mitigate this effect with ε 1 , but we also introduce ε 2 and ε 3 to specifically limit the importance in the metric of branches with persistence close to that of their parent. Fig. 20. Impact of the parameter ε 1 on geodesic computation (left: ε 1 = 0, right: ε 1 = 0.05). In this example (left), the white branch in T ( f i ) is not matched to the white branch in T ( f j ) as they have distinct depths in the corresponding BDTs (2 versus 1). However, these features are visually similar in the data (Gaussians with the white maximum in f i and f j , bottom left corner of the domain). With ε 1 = 0.05 (right), the saddle of the white branch in T ( f i ) gets merged with its ancestor saddle (whose f i value was less than ε 1 away). Consequently, the white branch gets moved up the BDT (the white branch is attached to the main light blue branch in T ( f i ), right). Since they now have identical depths in the corresponding BDTs, the white branches of T ( f i ) and T ( f j ) can now be matched together (right), which results in an overall matching (and geodesic) between these two trees which better conveys the resemblance between the two scalar fields f i and f j . Equivalently, one can interpret this procedure of saddle merge in the input trees as a modification of the input scalar field, turning f i into f j . In particular, this field modification disconnects the Gaussian with the white maximum from the Gaussian with the dark blue maximum ( f i ) and reconnects it to the Gaussian with the light blue maximum ( f j ). Fig. 21. Impact of the parameter ε 2 on geodesic computation (left: ε 2 = 1, right: ε 2 = 0.95). In this example (left), the white branch in T ( f i ) is not matched to the white branch in T ( f j ) as they have distinct depths in the corresponding BDTs (2 versus 1). Moreover, given the function difference between the white branch's saddle and its ancestor, that branch cannot be moved up the BDT under the effect of the ε 1 procedure (above). The white branch in T ( f i ) has a persistence nearly identical to its parent (cyan). Thus, after local normalization (necessary to guarantee the topological consistency of the interpolated trees), its normalized persistence would become artificially high, which can have an undesirable effect on the metric. The BDT pre-processing addresses this issue and moves up the BDT branches with a relative persistence to their parent larger than ε 2 (recommended default value: 0.95). In this example (right), the white branch in T ( f i ) moves up the BDT and becomes adjacent to the main light blue branch in T ( f i ). Since they now have identical depths in the corresponding BDTs, the white branches of T ( f i ) and T ( f j ) can now be matched together (right), which better conveys the resemblance between the two scalar fields f i and f j . Equivalently, one can interpret this procedure of BDT pre-processing as a modification of the input scalar field, turning f i into f j . In particular, this field modification disconnects the Gaussian with the white maximum from the Gaussian with the cyan maximum ( f i ) and reconnects it to the Gaussian with the light blue maximum ( f j ). Fig. 22. Effect of the parameter ε 3 on geodesic computation (left: ε 3 = 0, right: ε 3 = 0.9). In this example (left), the white branch in T ( f i ) is not matched to the white branch in T ( f j ) as they have distinct depths in the corresponding BDTs (3 versus 1). Applying the above BDT pre-processing (ε 2 ) to all branches would move the cyan branch in T ( f i ) up the BDT, which would prevent it to match to the cyan branch in T ( f j ). The parameter ε 3 restricts the application of the above BDT pre-processing and prevents the movement of the most persistent branches (relative persistence larger than ε 3 , default: 0.9). In this example (right), the white branch in T ( f i ) moves up the BDT and becomes adjacent to the main light blue branch in T ( f i ). Since they now have identical depths in the corresponding BDTs, the white branches of T ( f i ) and T ( f j ) can now be matched together (right), which results in an overall matching (and geodesic) between these two trees which better conveys the resemblance between the two scalar fields f i and f j . Equivalently, one can interpret this procedure on the BDTs as a modification of the input scalar field, turning f i into f j . In particular, this field modification disconnects the Gaussian with the white maximum from the Gaussian with the dark green maximum ( f i ) and reconnects it to the Gaussian with the light blue maximum ( f j ).