Bug with `LOWER_ROW` edge weight

[WGMRouting]

This is related to #40. The from_eilon workaround introduced there replaced a line-wrapping bug with an indexing bug.
Both code paths for parsing EDGE_WEIGHT_FORMAT: LOWER_ROW produce incorrect distance matrices. This is demonstrable against the published BKS for E-n13-k4 (optimal value 247).

The Two Bugs

1. from_lower_row incorrectly depends on physical line wrapping

The standard function executes distances[i + 1, : i + 1] = triangular[i], expecting each inner list to correspond to exactly one triangular row. TSPLIB95 specifies that EDGE_WEIGHT_SECTION is a continuous 1D stream where line breaks carry no structural meaning. If a physical line contains values spanning multiple triangular rows (or if all values are on a single line), the parser scrambles the matrix or crashes.

2. from_eilon uses the wrong (i,j) mapping (UPPER_ROW indexing)

The Eilon workaround correctly flattens the data to 1D but maps values using itertools.combinations(range(n), r=2). This produces pairs in (0,1), (0,2), (0,3) order—which is the UPPER_ROW convention. The LOWER_ROW order per TSPLIB95 is (1,0), (2,0), (2,1), (3,0). These diverge from the third value onward, silently assigning distances to the wrong node pairs. (Note: If this was implemented to fix a previous issue, it replaced a line-wrapping bug with an indexing bug).

Proof with E-n13-k4

The CVRPLIB BKS for E-n13-k4 is 247. This instance uses EDGE_WEIGHT_FORMAT: LOWER_ROW.

• vrplib v2.1.0 (from_eilon path): d(0,3) = 21, d(1,2) = 5, solver optimal = 290 (Fails)
• Correct LOWER_ROW (flat stream, strictly lower-triangular): d(0,3) = 23, d(1,2) = 21, solver optimal = 247 (Matches BKS)

Minimal Reproducible Examples

Here are two minimal .vrp files proving both paths fail. Both represent the exact same LOWER_ROW sequence: 10, 20, 30, 40, 50, 60.

MRE 1: Testing from_eilon indexing bug

NAME : test_eilon
COMMENT : (Eilon et al, test)
TYPE : CVRP
DIMENSION : 4
EDGE_WEIGHT_TYPE : EXPLICIT
EDGE_WEIGHT_FORMAT : LOWER_ROW
EDGE_WEIGHT_SECTION
10
20 30
40 50 60
DEMAND_SECTION
1 0
2 10
3 10
4 10
DEPOT_SECTION
1
-1
EOF

Expected: d(0,3) = 40, d(1,2) = 30
Actual (vrplib v2.1.0): d(0,3) = 30, d(1,2) = 40 (Wrong values due to combinations mapping)

MRE 2: Testing from_lower_row physical line bug

NAME : test_standard
COMMENT : standard test
TYPE : CVRP
DIMENSION : 4
EDGE_WEIGHT_TYPE : EXPLICIT
EDGE_WEIGHT_FORMAT : LOWER_ROW
EDGE_WEIGHT_SECTION
10 20 30 40 50 60
DEMAND_SECTION
1 0
2 10
3 10
4 10
DEPOT_SECTION
1
-1
EOF

Expected: The matrix parses successfully with d(0,3) = 40, d(1,2) = 30.
Actual (vrplib v2.1.0): ValueError: could not broadcast input array from shape (6,) into shape (1,) (Fails to parse entirely)

Suggested Fix

Replace both from_lower_row and the from_eilon special case with a single, vectorised NumPy function that flattens incoming data natively (handling both ragged nested arrays and flat arrays) and mathematically enforces the TSPLIB95 bounds.

import math
import itertools
import numpy as np

def from_lower_row(data):
    # 1. Robust flattening: Handle both ragged nested sequences and flat arrays
    if isinstance(data, np.ndarray):
        flattened = np.asarray(data).ravel()
    elif isinstance(data, (list, tuple)):
        # If the first element is also a sequence, flatten the nested structure
        if len(data) > 0 and isinstance(data[0], (list, tuple, np.ndarray)):
            flattened = np.fromiter(itertools.chain.from_iterable(data), dtype=float)
        else:
            flattened = np.asarray(data, dtype=float).ravel()
    else:
        flattened = np.asarray(data, dtype=float).ravel()

    m = flattened.size

    # 2. Exact mathematical inverse using the quadratic formula
    # n^2 - n - 2m = 0  =>  n = (1 + sqrt(1 + 8m)) / 2
    disc = 1 + 8 * m
    root = math.isqrt(disc)
    
    if root * root != disc:
        raise ValueError(f"Invalid LOWER_ROW size: {m} elements is not a valid triangular sequence.")

    n = (1 + root) // 2

    # 3. Vectorized assignment using true row-major lower triangular indexing
    distances = np.zeros((n, n), dtype=flattened.dtype)
    distances[np.tril_indices(n, k=-1)] = flattened
    distances += distances.T  # Mirror across diagonal for symmetry

    return distances

Implementation & Cleanup Notes

I have tested this locally by patching parse_distances.py. To implement this cleanly, the following legacy elements should also be removed from that file:

1. Remove Eilon Routing: Delete the if comment is not None and "Eilon" in comment: logic from parse_distances(). The new function handles Eilon instances natively.
2. Delete Dead Functions: The from_eilon() and is_triangular_number() functions are now entirely obsolete and should be removed.
3. Clean Imports: Remove from itertools import combinations, but ensure import math and import itertools are added at the top.

This eliminates the fragile comment-string detection, making the parser robust and mathematically correct across all compliant instances.
Tested against vrplib v2.1.0.

[leonlan]

@WGMRouting thanks for reporting this issue. I'll have a look at this later next week.

[leonlan]

Fixed and released in v2.2.0.

[WGMRouting]

Happy to help! Thanks for the rapid response and fixing/implementing quickly