Merge pull request #12 from Jij-Inc/09-routing

Introduce 09-routing

Author: ksk-jij

ommx_quantum_benchmarks/qoblib/09_routing/models/integer_linear/dat_reader.py

@@ -0,0 +1,151 @@
+from __future__ import annotations
+from math import sqrt
+
+
+def read_vrp_tsplib(
+    path: str,
+    vehicle_limit: int,
+    *,
+    euc2d_round: bool = True,  # TSPLIB EUC_2D convention: round distances to nearest integer
+    depot_policy: str = "first",  # Handling of multiple depots (only "first" supported)
+) -> dict:
+    """Read a TSPLIB-style CVRP instance and return JijModeling instance data.
+
+    Supports the standard TSPLIB format for the Capacitated Vehicle Routing
+    Problem (CVRP), including parsing coordinates, demands, and depot
+    definitions. Distances are computed with the EUC_2D convention.
+
+    Args:
+        path (str): Path to the TSPLIB-format CVRP file.
+        vehicle_limit (int): Maximum number of vehicles allowed.
+        euc2d_round (bool, optional): Whether to round EUC_2D distances to
+            the nearest integer (default True, per TSPLIB convention).
+        depot_policy (str, optional): Policy for selecting the depot if
+            multiple depots are provided. Currently only "first" is supported.
+
+    Returns:
+        Dict: A dictionary formatted for JijModeling instance_data with keys:
+            - "n" (int): Number of nodes (dimension).
+            - "VEHICLE_LIMIT" (int): Vehicle limit.
+            - "CAPACITY" (int): Vehicle capacity.
+            - "DEMAND" (List[int]): Demand at each node (0-based).
+            - "D" (List[List[int|int]]): Distance matrix (symmetric).
+            - "DEPOT" (int): Depot index (0-based).
+    """
+    # --- 1) Read file lines ---
+    with open(path, "r", encoding="utf-8") as f:
+        raw_lines = [ln.strip() for ln in f]
+
+    lines = [ln for ln in raw_lines if ln != ""]
+
+    # --- 2) Parse header ---
+    header: dict[str, str] = {}
+    idx = 0
+    while idx < len(lines):
+        ln = lines[idx]
+        if ln.upper().endswith("SECTION") or ln.upper() == "EOF":
+            break
+        if ":" in ln:
+            key, val = ln.split(":", 1)
+            header[key.strip().upper()] = val.strip()
+        idx += 1
+
+    try:
+        dim = int(header["DIMENSION"])
+        cap = int(header["CAPACITY"])
+        edge_type = header.get("EDGE_WEIGHT_TYPE", "EUC_2D").upper()
+        if edge_type not in ("EUC_2D",):
+            raise ValueError(f"Unsupported EDGE_WEIGHT_TYPE: {edge_type}")
+    except KeyError as e:
+        raise ValueError(f"Missing required header field: {e}")
+
+    # --- 3) Section indices ---
+    def find_section(name: str) -> int:
+        nameU = name.upper()
+        for k in range(idx, len(lines)):
+            if lines[k].upper().startswith(nameU):
+                return k
+        return -1
+
+    sec_coord = find_section("NODE_COORD_SECTION")
+    sec_demand = find_section("DEMAND_SECTION")
+    sec_depot = find_section("DEPOT_SECTION")
+    sec_eof = find_section("EOF")
+    if sec_coord < 0 or sec_demand < 0 or sec_depot < 0:
+        raise ValueError(
+            "Missing one of required sections: NODE_COORD_SECTION / DEMAND_SECTION / DEPOT_SECTION"
+        )
+    if sec_eof < 0:
+        sec_eof = len(lines)
+
+    # --- 4) Node coordinates ---
+    coords: list[tuple[int, int]] = []
+    k = sec_coord + 1
+    while k < len(lines) and k < sec_demand:
+        ln = lines[k]
+        if ln.upper().endswith("SECTION"):
+            break
+        parts = ln.split()
+        if len(parts) >= 3:
+            x = int(parts[1])
+            y = int(parts[2])
+            coords.append((x, y))
+        k += 1
+    if len(coords) != dim:
+        raise ValueError(f"NODE_COORD_SECTION count {len(coords)} != DIMENSION {dim}")
+
+    # --- 5) Demands ---
+    demand: list[int] = [0] * dim
+    k = sec_demand + 1
+    while k < len(lines) and k < sec_depot:
+        ln = lines[k]
+        if ln.upper().endswith("SECTION"):
+            break
+        parts = ln.split()
+        if len(parts) >= 2:
+            idx1 = int(parts[0])
+            dem = int(parts[1])
+            demand[idx1 - 1] = dem
+        k += 1
+
+    # --- 6) Depot section ---
+    depots_1b: list[int] = []
+    k = sec_depot + 1
+    while k < len(lines) and k < sec_eof:
+        ln = lines[k]
+        if ln.upper().endswith("SECTION"):
+            break
+        v = ln.split()[0]
+        if v == "-1":
+            break
+        depots_1b.append(int(v))
+        k += 1
+    if not depots_1b:
+        raise ValueError("DEPOT_SECTION is empty")
+    if depot_policy != "first":
+        raise ValueError(f"Unsupported depot_policy={depot_policy}")
+    depot0 = depots_1b[0] - 1
+
+    # --- 7) Distance matrix ---
+    def euc2d(a: tuple[int, int], b: tuple[int, int]) -> float:
+        return sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2)
+
+    D: list[list[float]] = [[0.0] * dim for _ in range(dim)]
+    for u in range(dim):
+        for v in range(dim):
+            if u == v:
+                D[u][v] = 0.0
+            else:
+                d = euc2d(coords[u], coords[v])
+                D[u][v] = round(d) if euc2d_round else d
+
+    # --- 8) Return instance_data ---
+    instance_data = {
+        "n": dim,
+        "VEHICLE_LIMIT": int(vehicle_limit),
+        "CAPACITY": cap,
+        "DEMAND": demand,
+        "D": D,
+        "DEPOT": depot0,
+    }
+    return instance_data

ommx_quantum_benchmarks/qoblib/09_routing/models/integer_linear/model.py

@@ -0,0 +1,96 @@
+import jijmodeling as jm
+
+
+def build_vrp_ilp() -> jm.Problem:
+    """Create Vehicle Routing Problem (VRP) ILP model.
+
+    Formulates the capacitated vehicle routing problem with a single depot
+    using a standard integer linear programming approach. The objective is
+    to minimize total travel distance subject to vehicle limit and capacity
+    constraints.
+
+    Parameters (placeholders):
+        - n (int): Number of nodes (including depot).
+        - VEHICLE_LIMIT (int): Maximum number of vehicles allowed.
+        - CAPACITY (int): Capacity of each vehicle.
+        - DEMAND (ndarray, shape (n,)): Demand at each node.
+        - D (ndarray, shape (n, n)): Distance matrix.
+        - DEPOT (int): Index of the depot node (0-based).
+
+    Variables:
+        - x[i, j] ∈ {0,1}: 1 if arc (i → j) is used, 0 otherwise.
+        - y[i] ∈ ℤ, [0, CAPACITY]: Load upon arrival at node i.
+
+    Objective:
+        - Minimize total distance: Σ_{i, j} D[i, j] · x[i, j].
+
+    Constraints:
+        - Each customer visited exactly once (excl. depot).
+        - Flow conservation for all non-depot nodes.
+        - Departures from depot ≤ VEHICLE_LIMIT.
+        - Capacity propagation (MTZ-style) to prevent subtours.
+        - Capacity bounds at all customer nodes.
+
+    Returns:
+        jm.Problem: JijModeling problem instance with variables, objective,
+        and constraints for the capacitated VRP.
+    """
+    # ------------ Placeholders ------------
+    n = jm.Placeholder("n")
+    VEHICLE_LIMIT = jm.Placeholder("VEHICLE_LIMIT")
+    CAPACITY = jm.Placeholder("CAPACITY")
+    DEMAND = jm.Placeholder("DEMAND", ndim=1)  # shape (n,)
+    D = jm.Placeholder("D", ndim=2)  # shape (n,n)
+    DEPOT = jm.Placeholder("DEPOT")  # single depot index (0-based)
+
+    # ------------ Indices ------------
+    i = jm.Element("i", belong_to=(0, n))
+    j = jm.Element("j", belong_to=(0, n))
+    h = jm.Element("h", belong_to=(0, n))
+
+    # ------------ Variables ------------
+    x = jm.BinaryVar("x", shape=(n, n), description="arc i->j used")
+    y = jm.IntegerVar(
+        "y",
+        shape=(n,),
+        lower_bound=0,
+        upper_bound=CAPACITY,
+        description="load upon arrival at node",
+    )
+
+    # ------------ Problem & Objective ------------
+    problem = jm.Problem("vrp_ilp", sense=jm.ProblemSense.MINIMIZE)
+    problem += jm.sum([i, j], D[i, j] * x[i, j])  # (10) minimize total distance
+
+    # ------------ Constraints ------------
+    # (11) Each customer visited exactly once (excluding depot)
+    problem += jm.Constraint(
+        "customer_visited_once",
+        jm.sum([(j, j != i)], x[i, j]) == 1,
+        forall=[(i, i != DEPOT)],
+    )
+
+    # (12) Flow conservation for non-depot nodes
+    problem += jm.Constraint(
+        "flow_conservation",
+        jm.sum([(i, i != h)], x[i, h]) - jm.sum([(i, i != h)], x[h, i]) == 0,
+        forall=[(h, h != DEPOT)],
+    )
+
+    # (13) Vehicle limit: departures from depot ≤ VEHICLE_LIMIT
+    problem += jm.Constraint(
+        "vehicle_limit", jm.sum([(j, j != DEPOT)], x[DEPOT, j]) <= VEHICLE_LIMIT
+    )
+
+    # (14) Capacity propagation (MTZ-style), exclude depot and i=j
+    problem += jm.Constraint(
+        "capacity_limit",
+        y[j] >= y[i] + DEMAND[j] * x[i, j] - CAPACITY * (1 - x[i, j]),
+        forall=[i, (j, (j != DEPOT) & (j != i))],
+    )
+
+    # (15) Capacity bounds at nodes
+    problem += jm.Constraint("capacity_limit_node_ub", y[i] <= CAPACITY, forall=[i])
+    problem += jm.Constraint("capacity_limit_node_lb", DEMAND[i] <= y[i], forall=[i])
+
+    return problem

ommx_quantum_benchmarks/qoblib/09_routing/models/integer_linear/ommx_create.py

@@ -0,0 +1,114 @@
+import os
+import glob
+import jijmodeling as jm
+from ommx.artifact import ArtifactBuilder
+from model import build_vrp_ilp
+from dat_reader import read_vrp_tsplib
+from sol_reader import parse_vrp_solution_file
+
+
+def _pick_solution_file(sol_dir: str, base: str) -> str | None:
+    """Pick first existing solution file for a basename."""
+    candidates = [
+        os.path.join(sol_dir, f"{base}.opt.sol"),
+        os.path.join(sol_dir, f"{base}.best.sol"),
+        os.path.join(sol_dir, f"{base}.bst.sol"),
+        os.path.join(sol_dir, f"{base}.sol"),
+    ]
+    for p in candidates:
+        if os.path.exists(p):
+            return p
+    return None
+
+
+def batch_process(
+    inst_dir: str = "../../instances",
+    sol_root: str = "../../solutions",
+    output_directory: str = "./ommx_output",
+):
+    """
+    Process instances from a QBench JSON file and corresponding solution files,
+    convert them into OMMX artifacts, and save them to the output directory.
+
+    Parameters:
+        inst_dir (str): Path to the gph..
+        sol_root (str): Path to the root solutions directory.
+        output_directory (str): Path to save the generated .ommx files.
+    """
+    os.makedirs(output_directory, exist_ok=True)
+
+    problem = build_vrp_ilp()
+
+    vrp_paths = sorted(glob.glob(os.path.join(inst_dir, "*.vrp")))
+    if not vrp_paths:
+        print(f"[MIS] No .vrp files found under: {inst_dir}")
+        return
+
+    processed_count = 0
+    error_count = 0
+
+    for vrp_path in vrp_paths:
+        base = os.path.splitext(os.path.basename(vrp_path))[0]
+        try:
+            print(f"[{base}] Reading: {vrp_path}")
+            instance_data = read_vrp_tsplib(vrp_path, vehicle_limit=4, euc2d_round=True)
+
+            interpreter = jm.Interpreter(instance_data)
+            ommx_instance = interpreter.eval_problem(problem)
+
+            sol_path = _pick_solution_file(sol_root, base)
+            solution = None
+            if sol_path:
+                try:
+                    print(f"  → Evaluating solution: {sol_path}")
+                    objective_value, solution_dict = parse_vrp_solution_file(
+                        sol_path,
+                        depot=0,
+                        demand=instance_data["DEMAND"],
+                        n=instance_data["n"],
+                    )
+                    solution = ommx_instance.evaluate(solution_dict)
+                    if (
+                        solution.feasible
+                        and abs(solution.objective - objective_value) < 1e-6
+                    ):
+                        print(
+                            f"    objective={solution.objective}, feasible={solution.feasible}"
+                        )
+                    else:
+                        print(
+                            "    ! Objective mismatch or infeasible; will save instance only."
+                        )
+                        solution = None
+                except Exception as sol_err:
+                    print(f"    ! Solution evaluation failed: {sol_err}")
+                    solution = None
+
+            out_path = os.path.join(output_directory, f"{base}.ommx")
+            if os.path.exists(out_path):
+                os.remove(out_path)
+
+            builder = ArtifactBuilder.new_archive_unnamed(out_path)
+            builder.add_instance(ommx_instance)
+            if solution is not None:
+                builder.add_solution(solution)
+            builder.build()
+
+            print(f"  ✓ Created: {out_path}")
+            print("-" * 50)
+            processed_count += 1
+
+        except Exception as e:
+            print(f"[{base}] Error: {e}")
+            print("-" * 50)
+            error_count += 1
+
+    print(f"Batch complete — processed: {processed_count}, errors: {error_count}")
+
+
+if __name__ == "__main__":
+    batch_process(
+        inst_dir="../../instances",
+        sol_root="../../solutions",
+        output_directory="./ommx_output",
+    )