VRPLIB/vrplib/parse/parse_distances.py at 7cb9794c2a3534df29253644cd5caa1678b5a9de · PyVRP/VRPLIB · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
import numpy as np


def parse_distances(
    data: list[float],
    edge_weight_type: str,
    edge_weight_format: str | None = None,
    node_coord: np.ndarray | None = None,
    comment: str | None = None,
    **kwargs: float | str | np.ndarray,
) -> np.ndarray:
    """
    Parses the distances. The specification "edge_weight_type" describes how
    the distances should be parsed. The two main ways are to calculate the
    Euclidean distances using the the node coordinates or by parsing an
    explicit distance matrix.

    Parameters
    ----------
    data
        The edge weight data.
    edge_weight_type
        The type of the edge weight data.
    edge_weight_format, optional
        The format of the edge weight data.
    node_coord, optional
        The customer location coordinates.
    comment, optional
        The comment specification in the instance.
    **kwargs, optional
        Optional keyword arguments.

    Returns
    -------
    np.ndarray
        An n-by-n distances matrix.
    """
    if "2D" in edge_weight_type:  # Euclidean distance on node coordinates
        if node_coord is None:
            msg = (
                "Cannot compute Euclidean distances because node coordinates "
                "are not provided."
            )
            raise ValueError(msg)

        distance = pairwise_euclidean(node_coord)

        if edge_weight_type == "EUC_2D":
            return distance

        if edge_weight_type == "FLOOR_2D":
            return np.floor(distance)

        if edge_weight_type == "EXACT_2D":
            return np.round(distance * 1000)

        if edge_weight_type == "CEIL_2D":
            return np.ceil(distance)

    if edge_weight_type == "EXPLICIT":
        if edge_weight_format == "LOWER_ROW":
            return from_lower_row(data)

        if edge_weight_format == "FULL_MATRIX":
            return np.array(data)

    raise ValueError("Edge weight type or format unknown.")


def pairwise_euclidean(coords: np.ndarray) -> np.ndarray:
    """
    Computes the pairwise Euclidean distance between the passed-in coordinates.

    Parameters
    ----------
    coords
        An n-by-2 array of location coordinates.

    Returns
    -------
    np.ndarray
        An n-by-n Euclidean distances matrix.

    """
    coords = np.atleast_2d(coords)

    sq_sum = (coords**2).sum(axis=1)
    sq_dist = np.add.outer(sq_sum, sq_sum) - 2 * (coords @ coords.T)
    np.fill_diagonal(sq_dist, 0)  # avoids minor numerical issues

    return np.sqrt(sq_dist)


def from_lower_row(data: np.ndarray) -> np.ndarray:
    """
    Computes a full distances matrix from a LOWER_ROW edge weight section.

    The input is treated as a continuous 1D stream of values (as specified
    by TSPLIB95), regardless of how the values are wrapped across lines.

    Parameters
    ----------
    data
        Edge weight data, possibly as a ragged array of rows.

    Returns
    -------
    np.ndarray
        An n-by-n distances matrix.
    """
    flattened = np.concatenate(data).astype(float)

    # The flattened data represents the lower triangle of a symmetric matrix.
    # See https://en.wikipedia.org/wiki/Triangular_number.
    # m = n * (n - 1) / 2 => n = (1 + sqrt(1 + 8m)) / 2
    n = (1 + int((1 + 8 * flattened.size) ** 0.5)) // 2

    distances = np.zeros((n, n))
    distances[np.tril_indices(n, k=-1)] = flattened
    distances += distances.T

    return distances