# -*- coding: utf-8 -*-

    Implementation of the modified Gromov–Hausdorff (mGH) distance

    between compact metric spaces induced by unweighted graphs. This

    code complements the results from "Efficient estimation of a

    Gromov–Hausdorff distance between unweighted graphs" by V. Oles et

    al. (https://arxiv.org/pdf/1909.09772). The mGH distance was first

    defined in "Some properties of Gromov–Hausdorff distances" by F.

    Mémoli (Discrete & Computational Geometry, 2012).

    Author: Vladyslav Oles

    ===================================================================

    Usage examples:

    1) Estimating the mGH distance between 4-clique and single-vertex

    graph from their adjacency matrices. Note that it suffices to fill

    only the upper triangle of an adjacency matrix.

    >>> AG = [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 1], [0, 0, 0, 0]]

    >>> AH = [[0]]

    >>> lb, ub = gromov_hausdorff(AG, AH)

    >>> lb, ub

    (0.5, 0.5)

    2) Estimating the mGH distance between cycle graphs of length 2 and

    5 from their adjacency matrices. Note that the adjacency matrices

    can be given in both dense and sparse SciPy formats.

    >>> AI = np.array([[0, 1], [0, 0]])

    >>> AJ = sps.csr_matrix(([1] * 5, ([0, 0, 1, 2, 3], [1, 4, 2, 3, 4])), shape=(5, 5))

    >>> lb, ub = gromov_hausdorff(AI, AJ)

    >>> lb, ub

    (0.5, 1.0)

    3) Estimating all pairwise mGH distances between multiple graphs

    from their adjacency matrices as an iterable.

    >>> As = [AG, AH, AI, AJ]

    >>> lbs, ubs = gromov_hausdorff(As)

    >>> lbs

    array([[0. , 0.5, 0.5, 0.5],

           [0.5, 0. , 0.5, 1. ],

           [0.5, 0.5, 0. , 0.5],

           [0.5, 1. , 0.5, 0. ]])

    >>> ubs

    array([[0. , 0.5, 0.5, 0.5],

           [0.5, 0. , 0.5, 1. ],

           [0.5, 0.5, 0. , 1. ],

           [0.5, 1. , 1. , 0. ]])

    ===================================================================

    Notations:

    |X| denotes the number of elements in set X.

    X → Y denotes the set of all mappings of set X into set Y.

    V(G) denotes vertex set of graph G.

    mGH(X, Y) denotes the modified Gromov–Hausdorff distance between

    compact metric spaces X and Y.

    row_i(A) denotes the i-th row of matrix A.

    PSPS^n(A) denotes the set of all permutation similarities of

    all n×n principal submatrices of square matrix A.

    PSPS^n_{i←j}(A) denotes the set of all permutation similarities of

    all n×n principal submatrices of square matrix A whose i-th row is

    comprised of the entries in row_j(A).

    ===================================================================

    Glossary:

    Distance matrix of metric space X is a |X|×|X| matrix whose

    (i, j)-th entry holds the distance between i-th and j-th points of

    X. By the properties of a metric, distance matrices are symmetric

    and non-negative, their diagonal entries are 0 and off-diagonal

    entries are positive.

    Curvature is a generalization of distance matrix that allows

    repetitions in the underlying points of a metric space. Curvature

    of an n-tuple of points from metric space X is an n×n matrix whose

    (i, j)-th entry holds the distance between the points from i-th and

    j-th positions of the tuple. Since these points need not be

    distinct, the off-diagonal entries of a curvature can equal 0.

    n-th curvature set of metric space X is the set of all curvatures

    of X that are of size n×n.

    d-bounded curvature for some d > 0 is a curvature whose

    off-diagonal entries are all ≥ d.

    Positive-bounded curvature is a curvature whose off-diagonal

    entries are all positive, i.e. the points in the underlying tuple

    are distinct. Equivalently, positive-bounded curvatures are

    distance matrices on the subsets of a metric space.

"""

# To sample √|X| * log (|X| + 1) mappings from X → Y by default.

        AG, AH=None, mapping_sample_size_order=DEFAULT_MAPPING_SAMPLE_SIZE_ORDER):

    """

    Estimate the mGH distance between simple unweighted graphs,

    represented as compact metric spaces based on their shortest

    path lengths.

    Parameters

    -----------

    AG: np.array (|V(G)|×|V(G)|)

        (Sparse) adjacency matrix of graph G, or an iterable of

        adjacency matrices if AH=None.

    AH: np.array (|V(H)|×|V(H)|)

        (Sparse) adjacency matrix of graph H, or None.

    mapping_sample_size_order: np.array (2)

        Parameter that regulates the number of mappings to sample when

        tightening upper bound of the mGH distance.

    Returns

    --------

    lb: float

        Lower bound of the mGH distance, or a square matrix holding

        lower bounds of pairwise mGH distances if AH=None.

    ub: float

        Upper bound of the mGH distance, or a square matrix holding

        upper bounds of pairwise mGH distances if AH=None.

    """

    # Form iterable with adjacency matrices.

                             "2 graphs to discriminate")

    else:

    # Find lower and upper bounds of each pairwise mGH distance between

    # the graphs.

            # Transform adjacency matrices of a pair of graphs to

            # distance matrices.

            # Find lower and upper bounds of the mGH distance between

            # the pair of graphs.

                DX, DY, mapping_sample_size_order=mapping_sample_size_order)

        # Symmetrize matrices with lower and upper bounds of pairwise

        # mGH distances between the graphs.

    else:

    """

    Represent simple unweighted graph as compact metric space (with

    integer distances) based on its shortest path lengths.

    Parameters

    -----------

    AG: np.array (|V(G)|×|V(G)|)

        (Sparse) adjacency matrix of simple unweighted graph G.

    Returns

    --------

    DG: np.array (|V(G)|×|V(G)|)

        (Dense) distance matrix of the compact metric space

        representation of G based on its shortest path lengths.

    """

    # Convert adjacency matrix to SciPy format if needed.

    # Compile distance matrix of the graph based on its shortest path

    # lengths.

    # Ensure compactness of metric space, represented by distance

    # matrix.

        # Extract largest connected component of the graph.

    # Cast distance matrix to optimal integer type.

    """

    Given a metric space X induced by simple unweighted graph,

    cast its distance matrix to the smallest signed integer type,

    sufficient to hold all its entries.

    Parameters

    -----------

    DX: np.array (|X|×|X|)

        Distance matrix of X.

    Returns

    --------

    DX: np.array (|X|×|X|)

        Distance matrix of X, cast to optimal type.

    """

    # Type is signed integer to allow subtractions.

    """

    Determine smallest signed integer type sufficient to hold a value.

    Parameters

    -----------

    value: non-negative integer

    Returns

    --------

    optimal_int_type: np.dtype

        Optimal signed integer type to hold the value.

    """

                          if value <= np.iinfo(int_type).max)

    """

    For X, Y metric spaces induced by simple unweighted graphs, find

    lower and upper bounds of mGH(X, Y).

    Parameters

    ----------

    DX: np.array (|X|×|X|)

        (Integer) distance matrix of X.

    DY: np.array (|Y|×|Y|)

        (Integer) distance matrix of Y.

    mapping_sample_size_order: np.array (2)

        Parameter that regulates the number of mappings to sample when

        tightening upper bound of mGH(X, Y).

    Returns

    --------

    lb: float

        Lower bound of mGH(X, Y).

    ub: float

        Upper bound of mGH(X, Y).

    """

    # Ensure distance matrices are of integer type.

    # Cast distance matrices to signed integer type to allow

    # subtractions.

    # Estimate mGH(X, Y).

        DX, DY, mapping_sample_size_order=mapping_sample_size_order, double_lb=double_lb)

    """

    For X, Y metric spaces induced by simple unweighted graphs, find

    lower bound of mGH(X, Y).

    Parameters

    ----------

    DX: np.array (|X|×|X|)

        (Integer) distance matrix of X.

    DY: np.array (|Y|×|Y|)

        (Integer) distance matrix of Y.

    Returns

    --------

    double_lb: float

        Lower bound of 2*mGH(X, Y).

    """

    # Obtain trivial lower bound of 2*mGH(X, Y) from

    # 1) mGH(X, Y) ≥ 0.5*|diam X - diam Y|;

    # 2) if X and Y are not isometric, mGH(X, Y) ≥ 0.5.

    # Initialize lower bound of 2*mGH(X, Y).

    # Try tightening the lower bound.

        # Try proving 2*mGH(X, Y) ≥ d using d-bounded curvatures of

        # X and Y of size 3×3 or larger. 2×2 curvatures are already

        # accounted for in trivial lower bound.

    """

    Find a largest-size d-bounded curvature of metric space X induced

    by simple unweighted graph.

    Parameters

    ----------

    DX: np.array (|X|×|X|)

        (Integer) distance matrix of X.

    diam_X: int

        Largest distance in X.

    Returns

    --------

    K: np.array (n×n)

        d-bounded curvature of X of largest size; n ≤ |X|.

    """

    # Initialize curvature K with the distance matrix of X.

        # Pick a row (and column) with highest number of off-diagonal

        # distances < d, then with smallest sum of off-diagonal

        # distances ≥ d.

                      np.sum(np.ma.masked_less(K, d), axis=0).data

        # Remove the row and column from K.

    """

    For X, Y metric spaces induced by simple unweighted graph, try to

    confirm 2*mGH(X, Y) ≥ d using K, a d-bounded curvature of X.

    Parameters

    ----------

    d: int

        Lower bound candidate for 2*mGH(X, Y).

    K: np.array (n×n)

        d-bounded curvature of X; n ≥ 3.

    DY: np.array (|Y|×|Y|)

        Integer distance matrix of Y.

    max_diam: int

        Largest distance in X and Y.

    Returns

    --------

    lb_is_confirmed: bool

        Whether confirmed that 2*mGH(X, Y) ≥ d.

    """

    # If K exceeds DY in size, the Hausdorff distance between the n-th

    # curvature sets of X and Y is ≥ d, entailing 2*mGH(X, Y) ≥ d (from

    # Theorem A).

                      confirm_lb_using_bounded_curvature_row(d, K, DY, max_diam)

    """

    For X, Y metric spaces induced by simple unweighted graph, and K,

    a d-bounded curvature of X, try to confirm 2*mGH(X, Y) ≥ d using

    some row of K.

    Parameters

    ----------

    d: int

        Lower bound candidate for 2*mGH(X, Y).

    K: np.array (n×n)

        d-bounded curvature of X; n ≥ 3.

    DY: np.array (|Y|×|Y|)

        Integer distance matrix of Y; n ≤ |Y|.

    max_diam: int

        Largest distance in X and Y.

    Returns

    --------

    lb_is_confirmed: bool

        Whether confirmed that 2*mGH(X, Y) ≥ d.

    """

    # Represent row of K as distance distributions, and retain those

    # that are maximal by the entry-wise partial order.

        represent_distance_matrix_rows_as_distributions(K, max_diam))

    # Represent rows of DY as distance distributions.

    # For each i ∈ 1,...,n, check if ||row_i(K) - row_i(L)||_∞ ≥ d

    # ∀L ∈ PSPS^n(DY), which entails that the Hausdorff distance

    # between the n-th curvature sets of X and Y is ≥ d, and therefore

    # 2*mGH(X, Y) ≥ d (from Theorem B).

        # For fixed i, check if ||row_i(K) - row_i(L)||_∞ ≥ d

        # ∀L ∈ PSPS^n_{i←j}(DY)  ∀j ∈ 1,...,|Y|, which is equivalent to

        # ||row_i(K) - row_i(L)||_∞ ≥ d  ∀L ∈ PSPS^n(DY).

            # For fixed i and j, checking ||row_i(K) - row_i(L)||_∞ ≥ d

            # ∀L ∈ PSPS^n_{i←j}(DY) is equivalent to solving a linear

            # (bottleneck) assignment feasibility problem between the

            # entries of row_i(K) and row_j(DY).

                K_max_rows_distance_distributions[i], DY_rows_distance_distributions[j], d)

    """

    Given a metric space X induced by simple unweighted graph,

    represent each row of its distance matrix as the frequency

    distribution of its entries. Entry 0 in each row is omitted.

    Parameters

    ----------

    DX: np.array (n×n)

        (Integer) distance matrix of X.

    max_d: int

        Upper bound of the entries in DX.

    Returns

    --------

    DX_rows_distributons: np.array (n×max_d)

        Each row holds frequencies of each distance from 1 to

        max_d in the corresponding row of DX. Namely, the (i, j)-th

        entry holds the frequency of distance (max_d - j) in row_i(DX).

    """

    # Add imaginary part to distinguish identical distances from

    # different rows of D.

        DX + 1j * np.arange(len(DX))[:, None], return_counts=True)

    # Type is signed integer to allow subtractions.

    # Construct index pairs for distance frequencies, so that the

    # frequencies of larger distances appear on the left.

        (np.imag(unique_distances).astype(optimal_int_type),

         max_d - np.real(unique_distances).astype(max_d.dtype))

    # Fill frequency distributions of the rows of DX.

    # Remove (unit) frequency of distance 0 from each row.

    """

    Given frequency distributions of entries in M positive integer

    vectors of size p, find unique maximal vectors under the following

    (entry- wise) partial order: for v, u vectors, v < u if and only if

    there exists a bijection f: {1,...,p} → {1,...,p} such that

    v_k < u_{f(k)}  ∀k ∈ {1,...,p}.

    Parameters

    ----------

    distributions: np.array (M×max_d)

        Frequency distributions of entries in the m vectors; the

        entries are bounded from above by max_d.

    Returns

    --------

    unique_max_distributions: np.array (m×max_d)

        Unique frequency distributions of the maximal vectors; m ≤ M.

    """

        np.cumsum(distributions - distributions[:, None, :], axis=2)

        np.all(pairwise_distribution_differences >= 0, axis=2),

        np.any(pairwise_distribution_differences > 0, axis=2))

        # `np.unique` is not implemented in NumPy 1.12 (Python 3.4).

            {tuple(distribution) for distribution in distributions[distributions_are_max]})

    """

    For positie integer vectors v of size p and u of size q ≥ p, check

    if there exists injective f: {1,...,p} → {1,...,q}, such that

    |v_k - u_{f(k)}| < d  ∀k ∈ {1,...,p}

    Parameters

    ----------

    v_distribution: np.array (max_d)

        Frequency distribution of entries in v; the entries are bounded

        from above by max_d.

    u_distribution: np.array (max_d)

        Frequency distribution of entries in u; the entries are bounded

        from above by max_d.

    d: int

        d > 0.

    Returns

    --------

    is_assignment_feasible: bool

        Whether such injective f: {1,...,p} → {1,...,q} exists.

    """

        # Find reversed v distribution index of smallest v entries yet

        # to be assigned. Then find index in reversed u distribution of

        # smallest u entries to which the b entries can be assigned to.

                     if reversed_v_distribution[i] > 0)

            # All v entries are assigned.

        else:

        # Find reversed u distribution index of smallest u entries to

        # which v entries, corresponding to a given reversed v

        # distribution index, can be assigned to.

                                                 len(reversed_u_distribution) - 1) + 1)

                     if reversed_u_distribution[j] > 0)

            # No u entries left to assign the particular v entries to.

    # Copy to allow modifications and stay pure; reverse to be

    # compatible with distributions of different size.

    # Injectively assign v entries to u entries if their difference

    # is < d, going from smallest to largest entries in both v and u,

    # until all v entries are assigned or such assignment proves

    # infeasible.

        else:

    # The assignment is feasible if and only if for some injective f,

    # |v_k - u_{f(k)}| < d  ∀k ∈ {1,...,p}.

    """

    For X, Y metric spaces, find an upper bound of mGH(X, Y).

    Parameters

    ----------

    DX: np.array (|X|×|X|)

        Distance matrix of X.

    DY: np.array (|Y|×|Y|)

        Distance matrix of Y.

    mapping_sample_size_order: np.array (2)

        Parameter that regulates the number of mappings to sample when

        tightening the upper bound.

    double_lb: float

        Lower bound of 2*mGH(X, Y).

    Returns

    --------

    double_ub: float

        Upper bound of 2*mGH(X, Y).

    """

    # Find upper bound of smallest distortion of a mapping in X → Y.

        DX, DY, mapping_sample_size_order=mapping_sample_size_order, goal_distortion=double_lb)

    # Find upper bound of smallest distortion of a mapping in Y → X.

        DY, DX, mapping_sample_size_order=mapping_sample_size_order,

        goal_distortion=ub_of_X_to_Y_min_distortion)

                              mapping_sample_size_order=DEFAULT_MAPPING_SAMPLE_SIZE_ORDER,

                              goal_distortion=0):

    """

    For X, Y metric spaces, find an upper bound of smallest distortion

    of a mapping in X → Y by heuristically constructing some mappings

    and choosing the smallest distortion in the sample.

    Parameters

    ----------

    DX: np.array (|X|×|X|)

        Distance matrix of X.

    DY: np.array (|Y|×|Y|)

        Distance matrix of Y.

    mapping_sample_size_order: np.array (2)

        Exponents of |X| and log (|X|+1) in their product that defines

        how many mappings from X → Y to sample.

    goal_distortion: float

        No need to look for distortion smaller than this.

    Returns

    --------

    ub_of_min_distortion: float

        Upper bound of smallest distortion of a mapping in X → Y.

    """

    # Compute the numper of mappings to sample.

        np.array([len(DX), np.log(len(DX) + 1)]) ** mapping_sample_size_order)))

    # Construct each mapping in X → Y in |X| steps by choosing the image

    # of π(i)-th point in X at i-th step, where π is randomly sampled

    # |X|-permutation. Image of each point is chosen to minimize the

    # intermediate distortion at each step.

    """

    # For X, Y metric spaces and |X|-permutation π, construct a mapping

    from X → Y in |X| steps by choosing the image of π(i)-th point in X

    at i-th step. The image of each point is chosen to minimize the

    intermediate distortion at the corresponding step.

    DX: np.array (|X|×|X|)

        Distance matrix of X.

    DY: np.array (|Y|×|Y|)

        Distance matrix of Y.

    pi: np.array (|X|)

        |X|-permutation specifying the order in which the points in X

        are mapped.

    Returns

    --------

    mapped_xs_images: list

        image of the constructed mapping.

    distortion: float

        distortion of the constructed mapping.

    """

    # Map π(1)-th point in X to a random point in Y, due to the

    # lack of better criterion.

    # Map π(i)-th point in X for each i = 2,...,|X|.

        # Choose point in Y that minimizes the distortion after

        # mapping π(i)-th point in X to it.

            np.abs(DX[x, mapped_xs] - DY[:, mapped_xs_images]), axis=1)

        # Map π(i)-th point in X to the chosen point in Y.