#!/usr/bin/env python3
"""Benchmark: Finsler-Randers directed routing — QAP (MIP) vs QUBO (SA/MIP).

Compares four formulation strategies on the same Finsler-Randers test
problems at n=8, 12, 24 directions:

  1. QUBO-SA:        Simulated annealing on symmetric QUBO (existing baseline)
  2. QUBO-MIP:       Exact QUBO optimum via HiGHS MIP linearization
  3. QAP-LP:         Assignment LP relaxation of TSP (no subtour elimination)
  4. QAP-MIP:        TSP with MTZ subtour elimination via HiGHS MIP

Each result reports both symmetric energy (QUBO objective) and asymmetric
path cost (raw Q_ij sum over directed edges) to quantify what is lost
when the Finsler drift β is symmetrized away.
"""

import sys
import os
sys.path.insert(0, os.path.dirname(__file__))

import json
import math
import time
import numpy as np
import traceback

from qaoa_adapter import (
    FinslerMetric, finsler_metric_to_qubo, stochastic_abuse_qubo,
)

try:
    import highspy
    import numpy as np
    HAS_HIGHS = True
except ImportError:
    HAS_HIGHS = False


# ── Test problem generator (mirrors benchmark_finsler_qaoa.py) ─────────

def generate_finsler_problem(n_dirs: int, dimension: int = 3, seed: int = 42) -> tuple:
    rng = np.random.RandomState(seed)
    alpha_mass = rng.uniform(0.5, 2.0, size=dimension).tolist()
    beta_wind = rng.uniform(-0.3, 0.3, size=dimension).tolist()
    metric = FinslerMetric(alpha_mass=alpha_mass, beta_wind=beta_wind, dimension=dimension)
    directions = []
    for k in range(n_dirs):
        vec = rng.randn(dimension)
        norm = math.sqrt(sum(v*v for v in vec))
        directions.append([v / norm for v in vec])

    # Raw asymmetric cost matrix Q_ij = crossing_cost(i→j)
    raw = np.zeros((n_dirs, n_dirs))
    for i in range(n_dirs):
        for j in range(n_dirs):
            if i != j:
                raw[i, j] = metric.crossing_cost(directions[i], directions[j])

    # Symmetrized QUBO matrix (what finsler_metric_to_qubo produces)
    sym = np.zeros((n_dirs, n_dirs))
    for i in range(n_dirs):
        for j in range(i + 1, n_dirs):
            sym[i, j] = raw[i, j] + raw[j, i]

    max_aniso = float(np.max(np.abs(raw - raw.T)))

    return metric, directions, raw, sym, max_aniso


# ── Solvers ─────────────────────────────────────────────────────────────

def solve_qubo_sa(sym_matrix: np.ndarray, time_limit: float = 10.0, seed: int = 42) -> dict:
    n = len(sym_matrix)
    Q_dict = {}
    for i in range(n):
        for j in range(i + 1, n):
            if abs(sym_matrix[i, j]) > 1e-15:
                Q_dict[(i, j)] = sym_matrix[i, j]
    qubo = type("QUBO", (), {"n": n, "matrix": Q_dict, "energy": lambda self, x: sum(
        Q_dict.get((i, j), 0.0) * x[i] * x[j] for i, j in Q_dict)})()
    result = stochastic_abuse_qubo(qubo, method="sa", time_limit=time_limit, seed=seed)
    return result


def _build_binary_mip(num_vars: int, obj: np.ndarray,
                      col_entries: dict[int, list[tuple[int, float]]],
                      row_lower: np.ndarray, row_upper: np.ndarray,
                      time_limit: float) -> highspy.Highs:
    """Build and solve a binary MIP using direct HighsLp construction."""
    num_rows = len(row_lower)
    starts, indices_list, values_list, nnz = [], [], [], 0
    for col in range(num_vars):
        starts.append(nnz)
        for r, v in sorted(col_entries.get(col, []), key=lambda x: x[0]):
            indices_list.append(r)
            values_list.append(v)
            nnz += 1
    starts.append(nnz)

    lp = highspy.HighsLp()
    lp.num_col_ = num_vars
    lp.num_row_ = num_rows
    lp.col_cost_ = obj
    lp.col_lower_ = np.zeros(num_vars)
    lp.col_upper_ = np.ones(num_vars)
    lp.integrality_ = [highspy.HighsVarType.kInteger] * num_vars
    lp.a_matrix_ = highspy.HighsSparseMatrix()
    lp.a_matrix_.format_ = highspy.MatrixFormat.kColwise
    lp.a_matrix_.start_ = np.array(starts, dtype=np.int32)
    lp.a_matrix_.index_ = np.array(indices_list, dtype=np.int32)
    lp.a_matrix_.value_ = np.array(values_list)
    lp.row_lower_ = row_lower
    lp.row_upper_ = row_upper

    h = highspy.Highs()
    h.setOptionValue("time_limit", time_limit)
    h.setOptionValue("output_flag", False)
    h.passModel(lp)
    h.run()
    return h


def _mip_status(h: highspy.Highs) -> str:
    info = h.getInfo()
    val = info.primal_solution_status if hasattr(info, "primal_solution_status") else 0
    return {0: 'unknown', 1: 'infeasible', 2: 'feasible', 3: 'optimal'}.get(val, f'status_{val}')


def solve_qubo_mip(sym_matrix: np.ndarray, time_limit: float = 30.0,
                   diag_penalty: float = 0.0, card_k: int = None) -> dict:
    """Solve QUBO exactly via HiGHS MIP linearization.

    Args:
        sym_matrix: N×N symmetrized matrix (upper triangle stored in Q_dict)
        time_limit: solver time limit
        diag_penalty: negative diagonal term −λ added to each x_i (incentivizes selection)
        card_k: if set, enforce Σ x_i = card_k (cardinality constraint)

    Returns dict with solution, objective, status.
    """
    n = len(sym_matrix)
    if not HAS_HIGHS:
        return {"error": "highspy not installed"}

    Q_dict = {}
    for i in range(n):
        for j in range(i + 1, n):
            if abs(sym_matrix[i, j]) > 1e-15:
                Q_dict[(i, j)] = sym_matrix[i, j]

    pairs = list(Q_dict.keys())
    num_x = n
    num_y = len(pairs)
    num_vars = num_x + num_y
    y_offset = num_x
    pair_to_idx = {p: k for k, p in enumerate(pairs)}

    extra_row = 1 if card_k is not None else 0
    num_rows = 3 * num_y + extra_row

    # Objective
    obj = np.zeros(num_vars)
    for (i, j), coeff in Q_dict.items():
        obj[y_offset + pair_to_idx[(i, j)]] = coeff
    for i in range(n):
        obj[i] = diag_penalty

    # Constraint matrix
    col_entries: dict[int, list[tuple[int, float]]] = {}
    row_idx = 0
    for (i, j), yi in pair_to_idx.items():
        yv = y_offset + yi
        for col, val in [(yv, 1.0), (i, -1.0)]:
            col_entries.setdefault(col, []).append((row_idx, val))
        row_idx += 1
        for col, val in [(yv, 1.0), (j, -1.0)]:
            col_entries.setdefault(col, []).append((row_idx, val))
        row_idx += 1
        for col, val in [(yv, -1.0), (i, 1.0), (j, 1.0)]:
            col_entries.setdefault(col, []).append((row_idx, val))
        row_idx += 1

    if card_k is not None:
        for i in range(n):
            col_entries.setdefault(i, []).append((row_idx, 1.0))
        card_row = row_idx
        row_idx += 1

    # Row bounds:
    #   Row 3*k + 0: y - x_i <= 0     →  upper = 0
    #   Row 3*k + 1: y - x_j <= 0     →  upper = 0
    #   Row 3*k + 2: -y + x_i + x_j <= 1  →  upper = 1
    #   (for k = 0..num_y-1)
    row_lower = np.full(num_rows, -1e30)
    row_upper = np.full(num_rows, 0.0)
    row_upper[2:3 * num_y:3] = 1.0  # every 3rd row from offset 2
    if card_k is not None:
        row_lower[card_row] = float(card_k)
        row_upper[card_row] = float(card_k)

    h = _build_binary_mip(num_vars, obj, col_entries, row_lower, row_upper, time_limit)
    sol = h.getSolution()
    x_vals = list(sol.col_value)
    solution = [int(round(x_vals[i])) for i in range(n)]

    obj_val = sum(Q_dict.get((i, j), 0.0) * solution[i] * solution[j]
                  for i in range(n) for j in range(i + 1, n)) + diag_penalty * sum(solution)

    return {"solution": solution, "objective": obj_val, "status": _mip_status(h)}


def solve_tsp_lp(cost_matrix: np.ndarray, time_limit: float = 30.0) -> dict:
    """Assignment LP relaxation of TSP (no subtour elimination)."""
    n = len(cost_matrix)
    if n <= 2:
        return {"route": list(range(n)), "cost": float(np.sum(cost_matrix)), "status": "trivial"}
    if not HAS_HIGHS:
        return {"error": "highspy not installed"}

    import sys as _sys
    _sys.path.insert(0, os.path.dirname(__file__))
    from qubo_highs import solve_route_lp
    result = solve_route_lp(cost_matrix.tolist(), time_limit=time_limit)
    return result


def solve_tsp_mtz(cost_matrix: np.ndarray, time_limit: float = 60.0) -> dict:
    """TSP with Miller-Tucker-Zemlin subtour elimination via HiGHS MIP.

    Formulation:
      Variables:
        x_ij ∈ {0,1} for all i≠j (edge i→j is in tour)
        u_i  ∈ [0, n-1] for i=1..n-1 (position of node i in tour)

      Minimize  Σ_i Σ_j Q_ij * x_ij

      Subject to:
        Σ_j x_ij = 1  for all i          (flow out)
        Σ_i x_ij = 1  for all j          (flow in)
        u_i - u_j + n*x_ij ≤ n-1         for all i≠j, i,j ≥ 1  (MTZ)
    """
    n = len(cost_matrix)
    if n <= 2:
        return {"route": list(range(n)), "cost": float(np.sum(cost_matrix)), "status": "trivial"}
    if not HAS_HIGHS:
        return {"error": "highspy not installed"}

    # Variable layout:
    #   x_ij for i≠j: binary edge variables (n*(n-1) vars)
    #   u_i for i=1..n-1: continuous position variables (n-1 vars)
    edges = [(i, j) for i in range(n) for j in range(n) if i != j]
    num_edge_vars = len(edges)
    num_cont_vars = n - 1
    num_vars = num_edge_vars + num_cont_vars
    edge_idx = {e: k for k, e in enumerate(edges)}
    u_offset = num_edge_vars
    cont_nodes = list(range(1, n))  # nodes with u_i (skip 0)

    # Constraints:
    #  2n flow + (n-1)*(n-2) MTZ
    num_rows = 2 * n + (n - 1) * (n - 2)

    # Objective: Σ Q_ij * x_ij
    obj = np.zeros(num_vars)
    for (i, j), vi in edge_idx.items():
        obj[vi] = cost_matrix[i, j]

    # Bounds
    col_upper = np.ones(num_vars)
    for k in range(num_cont_vars):
        col_upper[u_offset + k] = float(n - 1)

    # Build constraint matrix (CSC)
    col_entries: dict[int, list[tuple[int, float]]] = {}
    row_ptr = 0

    for i in range(n):
        for j in range(n):
            if i != j:
                col_entries.setdefault(edge_idx[(i, j)], []).append((row_ptr, 1.0))
        row_ptr += 1

    for j in range(n):
        for i in range(n):
            if i != j:
                col_entries.setdefault(edge_idx[(i, j)], []).append((row_ptr, 1.0))
        row_ptr += 1

    for i in cont_nodes:
        for j in cont_nodes:
            if i == j:
                continue
            ui = u_offset + cont_nodes.index(i)
            uj = u_offset + cont_nodes.index(j)
            xij = edge_idx[(i, j)]
            col_entries.setdefault(ui, []).append((row_ptr, 1.0))
            col_entries.setdefault(uj, []).append((row_ptr, -1.0))
            col_entries.setdefault(xij, []).append((row_ptr, float(n)))
            row_ptr += 1

    # Row bounds
    row_lower = np.full(num_rows, -1e30)
    row_upper = np.full(num_rows, 1e30)
    row_lower[:2 * n] = 1.0
    row_upper[:2 * n] = 1.0
    row_upper[2 * n:] = float(n - 1)

    # Build LP directly
    starts, indices_list, values_list, nnz = [], [], [], 0
    for col in range(num_vars):
        starts.append(nnz)
        for r, v in sorted(col_entries.get(col, []), key=lambda x: x[0]):
            indices_list.append(r)
            values_list.append(v)
            nnz += 1
    starts.append(nnz)

    lp = highspy.HighsLp()
    lp.num_col_ = num_vars
    lp.num_row_ = num_rows
    lp.col_cost_ = obj
    lp.col_lower_ = np.zeros(num_vars)
    lp.col_upper_ = col_upper
    lp.a_matrix_ = highspy.HighsSparseMatrix()
    lp.a_matrix_.format_ = highspy.MatrixFormat.kColwise
    lp.a_matrix_.start_ = np.array(starts, dtype=np.int32)
    lp.a_matrix_.index_ = np.array(indices_list, dtype=np.int32)
    lp.a_matrix_.value_ = np.array(values_list)
    lp.row_lower_ = row_lower
    lp.row_upper_ = row_upper
    lp.integrality_ = [highspy.HighsVarType.kInteger] * num_edge_vars + [highspy.HighsVarType.kContinuous] * num_cont_vars

    h = highspy.Highs()
    h.setOptionValue("time_limit", time_limit)
    h.setOptionValue("output_flag", False)
    h.passModel(lp)
    h.run()

    sol = h.getSolution()
    x_vals = list(sol.col_value)

    # Extract tour
    next_node = {}
    for (i, j), vi in edge_idx.items():
        if x_vals[vi] > 0.5:
            next_node[i] = j

    route = [0]
    seen = {0}
    for _ in range(n - 1):
        nxt = next_node.get(route[-1])
        if nxt is None or nxt in seen:
            break
        route.append(nxt)
        seen.add(nxt)

    if len(route) < n:
        route = list(range(n))

    total_cost = sum(cost_matrix[route[i]][route[(i + 1) % n]] for i in range(n))
    obj_val = sum(cost_matrix[i, j] * x_vals[vi] for (i, j), vi in edge_idx.items())

    status_raw = h.getInfoValue("primal_solution_status")[1]
    status_str = {0: 'unknown', 1: 'infeasible', 2: 'feasible', 3: 'optimal'}.get(status_raw, f'status_{status_raw}')

    return {"route": route, "cost": float(total_cost), "objective": float(obj_val), "status": status_str}


# ── Evaluation helpers ──────────────────────────────────────────────────

def asymmetric_subset_cost(solution: list[int], raw: np.ndarray) -> float:
    """Σ_i Σ_j Q_ij * x_i * x_j using raw (asymmetric) Q_ij."""
    cost = 0.0
    for i in range(len(solution)):
        if solution[i]:
            for j in range(len(solution)):
                if solution[j] and i != j:
                    cost += raw[i, j]
    return cost


def symmetric_subset_cost(solution: list[int], sym: np.ndarray) -> float:
    cost = 0.0
    for i in range(len(solution)):
        if solution[i]:
            for j in range(i + 1, len(solution)):
                if solution[j]:
                    cost += sym[i, j]
    return cost


# ── Benchmark engine ────────────────────────────────────────────────────

def run_benchmark(sizes: list[int] = None, time_limits: dict[int, float] = None, seed: int = 42) -> dict:
    if sizes is None:
        sizes = [8, 12, 24]
    if time_limits is None:
        time_limits = {8: 30.0, 12: 60.0, 24: 120.0, 48: 120.0}

    results = {}
    for n in sizes:
        print(f"\n{'='*60}")
        print(f"  Finsler-Randers Routing  n={n}")
        print(f"{'='*60}")

        metric, directions, raw, sym, max_aniso = generate_finsler_problem(n, seed=seed)
        tl = time_limits.get(n, 60.0)
        print(f"  max |Q_ij - Q_ji| = {max_aniso:.6f}  (asymmetry from β wind field)")

        row = {"n": n, "max_anisotropy": max_aniso}

        # ── QUBO-SA (degenerate baseline) ──
        print(f"\n  ── Solver: QUBO-SA (degenerate, all-positive Q_ij) ──")
        try:
            t0 = time.time()
            sa = solve_qubo_sa(sym, time_limit=min(tl, 10.0))
            sa_time = time.time() - t0
            x_sa = sa.get("solution", [0] * n)
            row["qubo_sa"] = {
                "symmetric_energy": round(symmetric_subset_cost(x_sa, sym), 6),
                "asymmetric_cost": round(asymmetric_subset_cost(x_sa, raw), 6),
                "x": x_sa,
                "selected": sum(x_sa),
                "time_s": round(sa_time, 3),
                "status": "ok",
            }
            print(f"    sym_energy={row['qubo_sa']['symmetric_energy']}  asym_cost={row['qubo_sa']['asymmetric_cost']}  selected={sum(x_sa)}/{n}  time={sa_time:.3f}s")
        except Exception as e:
            row["qubo_sa"] = {"status": "error", "detail": str(e)}
            print(f"    FAILED: {e}")

        # ── QUBO-MIP (degenerate, no constraint) ──
        print(f"  ── Solver: QUBO-MIP (unconstrained, trivially selects 0 vars) ──")
        try:
            t0 = time.time()
            mip = solve_qubo_mip(sym, time_limit=min(tl, 10.0))
            mip_time = time.time() - t0
            x_mip = mip.get("solution", [0] * n)
            row["qubo_mip"] = {
                "symmetric_energy": round(symmetric_subset_cost(x_mip, sym), 6),
                "asymmetric_cost": round(asymmetric_subset_cost(x_mip, raw), 6),
                "x": x_mip,
                "selected": sum(x_mip),
                "time_s": round(mip_time, 3),
                "status": mip.get("status", "?"),
            }
            print(f"    sym_energy={row['qubo_mip']['symmetric_energy']}  asym_cost={row['qubo_mip']['asymmetric_cost']}  selected={sum(x_mip)}/{n}  time={mip_time:.3f}s  status={mip.get('status','?')}")
        except Exception as e:
            row["qubo_mip"] = {"status": "error", "detail": str(e)}
            print(f"    FAILED: {e}")

        # ── QUBO-MIP-card: cardinality-constrained (select K = n//2) ──
        K = max(1, n // 2)
        print(f"  ── Solver: QUBO-MIP-card (cardinality-constrained, K={K}) ──")
        try:
            t0 = time.time()
            mip_k = solve_qubo_mip(sym, time_limit=tl, card_k=K)
            mip_k_time = time.time() - t0
            x_mip_k = mip_k.get("solution", [0] * n)
            row["qubo_mip_card"] = {
                "symmetric_energy": round(symmetric_subset_cost(x_mip_k, sym), 6),
                "asymmetric_cost": round(asymmetric_subset_cost(x_mip_k, raw), 6),
                "x": x_mip_k,
                "cardinality": K,
                "selected": sum(x_mip_k),
                "time_s": round(mip_k_time, 3),
                "status": mip_k.get("status", "?"),
            }
            print(f"    sym_energy={row['qubo_mip_card']['symmetric_energy']}  asym_cost={row['qubo_mip_card']['asymmetric_cost']}  selected={sum(x_mip_k)}/{n}  time={mip_k_time:.3f}s  status={mip_k.get('status','?')}")
        except Exception as e:
            row["qubo_mip_card"] = {"status": "error", "detail": str(e)}
            print(f"    FAILED: {e}")

        # ── QAP-LP (assignment relaxation) ──
        print(f"  ── Solver: QAP-LP (assignment relaxation) ──")
        try:
            t0 = time.time()
            lp_res = solve_tsp_lp(raw, time_limit=tl)
            lp_time = time.time() - t0
            route_lp = lp_res.get("route", list(range(n)))
            tour_lp = [(route_lp[i], route_lp[(i + 1) % n]) for i in range(n)]
            row["qap_lp"] = {
                "path_cost": round(lp_res.get("cost", 0.0), 6),
                "time_s": round(lp_time, 3),
                "status": lp_res.get("status", "?"),
                "route": route_lp,
            }
            print(f"    path_cost={row['qap_lp']['path_cost']}  time={lp_time:.3f}s  status={lp_res.get('status','?')}")
        except Exception as e:
            row["qap_lp"] = {"status": "error", "detail": str(e)}
            print(f"    FAILED: {e}")

        # ── QAP-MIP (TSP with MTZ) ──
        print(f"  ── Solver: QAP-MIP (TSP MTZ, asymmetric preserved) ──")
        try:
            t0 = time.time()
            mtz_res = solve_tsp_mtz(raw, time_limit=tl)
            mtz_time = time.time() - t0
            route_mtz = mtz_res.get("route", list(range(n)))
            row["qap_mip"] = {
                "path_cost": round(mtz_res.get("cost", 0.0), 6),
                "time_s": round(mtz_time, 3),
                "status": mtz_res.get("status", "?"),
                "route": route_mtz,
            }
            print(f"    path_cost={row['qap_mip']['path_cost']}  time={mtz_time:.3f}s  status={mtz_res.get('status','?')}")
        except Exception as e:
            row["qap_mip"] = {"status": "error", "detail": str(e)[:200]}
            print(f"    FAILED: {e}")

        # ── Cross-evaluation ──
        print(f"\n  ── Cross-evaluation ──")

        def _nearest_neighbor_tour(nodes: list[int], cost_mat: np.ndarray) -> tuple:
            if len(nodes) <= 1:
                return nodes, 0.0
            unvisited = set(nodes)
            current = nodes[0]
            unvisited.remove(current)
            tour = [current]
            total = 0.0
            while unvisited:
                best_n = -1
                best_c = float("inf")
                for nxt in unvisited:
                    if cost_mat[current, nxt] < best_c:
                        best_c = cost_mat[current, nxt]
                        best_n = nxt
                tour.append(best_n)
                unvisited.remove(best_n)
                total += best_c
                current = best_n
            return tour, total

        for key in ("qubo_sa", "qubo_mip", "qubo_mip_card"):
            entry = row.get(key, {})
            if entry.get("status") in ("ok", "feasible", "optimal"):
                sol = entry.get("x", [])
                selected_nodes = [i for i, v in enumerate(sol) if v]
                if len(selected_nodes) >= 2:
                    tour, tc = _nearest_neighbor_tour(selected_nodes, raw)
                    entry["nn_tour"] = tour
                    entry["nn_tour_asymmetric_cost"] = round(tc, 6)

                    # Also solve exact TSP on just these K nodes using the asymmetric cost
                    K = len(selected_nodes)
                    idx_map = {j: i for i, j in enumerate(selected_nodes)}
                    sub_raw = np.array([[raw[i, j] for j in selected_nodes] for i in selected_nodes])
                    sub_tsp = solve_tsp_mtz(sub_raw, time_limit=min(30.0, tl))
                    if "error" not in sub_tsp and sub_tsp.get("route"):
                        remapped_route = [selected_nodes[i] for i in sub_tsp["route"]]
                        entry["tsp_on_subset_route"] = remapped_route
                        entry["tsp_on_subset_cost"] = round(sub_tsp["cost"], 6)
                        print(f"  {key}  NN-tour={tc:.4f}  TSP-on-subset={sub_tsp['cost']:.4f}  (K={K})")
                    else:
                        print(f"  {key}  NN-tour={tc:.4f}  (K={K}, TSP skipped)")

        # Evaluate TSP tour's subset energy
        qap_entry = row.get("qap_mip", {})
        if qap_entry.get("route"):
            x_tsp = [1] * n
            qap_entry["as_subset_energy"] = round(symmetric_subset_cost(x_tsp, sym), 6)
            qap_entry["as_subset_asymmetric_cost"] = round(asymmetric_subset_cost(x_tsp, raw), 6)
            print(f"  qap_mip (TSP all {n})  path_cost={qap_entry['path_cost']}")

        print(f"\n  ── Summary ──")
        row["max_anisotropy"] = max_aniso

        results[f"n_{n}"] = row

    return results


if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description="Finsler-Randers QAP vs QUBO benchmark")
    parser.add_argument("--sizes", type=int, nargs="+", default=[8, 12, 24])
    parser.add_argument("--time-limit", type=float, default=None,
                        help="Solver time limit (default: per-size limits)")
    parser.add_argument("--seed", type=int, default=42)
    parser.add_argument("--output", "-o", help="JSON output path")
    args = parser.parse_args()

    if args.time_limit is not None:
        tl = {s: args.time_limit for s in args.sizes}
    else:
        tl = None

    results = run_benchmark(sizes=args.sizes, time_limits=tl, seed=args.seed)

    output = json.dumps(results, indent=2)
    if args.output:
        with open(args.output, "w") as f:
            f.write(output)
        print(f"\nWrote results to {args.output}")
    else:
        print(f"\n{'-'*60}")
        print(output)