diff --git a/main.py b/ tsp solver.py similarity index 98% rename from main.py rename to tsp solver.py index 24c9580..09e110e 100644 --- a/main.py +++ b/ tsp solver.py @@ -2,11 +2,8 @@ # Utilities from pathlib import Path -import time -import matplotlib.pyplot as plt -import random -from constants import * -from util import * +from constants import N_RUNS, MAX_SECONDS, MAX_ITERATIONS +from util import find_optimal_tour, setup_algorithm, run_single_trial_by_timing, run_single_iteration_trial import sys import logging from multiprocessing import Pool, cpu_count diff --git a/README.md b/README.md index cfe0427..ec1029d 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ The program takes a .tsp file as input and produces: For example, for a solution visiting cities 5, 4, 1, 3, 2 with distance 8934.12: ``` -> python main.py aaa.tsp +> python tsp_solver.py aaa.tsp 8934.12 > cat solution.csv 5 @@ -32,16 +32,9 @@ For example, for a solution visiting cities 5, 4, 1, 3, 2 with distance 8934.12: The submission includes the implementation files and a detailed report (as .pdf) describing the solution, approach, and optimizations implemented. -Note: The submission must be self-contained, with no dependencies on external files. Solvers should work out of the box, with reasonable documentation. - ## Dataset Setup -This project includes a dataset setup script that downloads and filters TSP instances from TSPLIB95. The script automatically: - -1. Downloads the complete TSPLIB95 dataset -2. Filters for TSP instances with `TYPE: TSP` and `EDGE_WEIGHT_TYPE: EUC_2D` -3. Extracts corresponding optimal tour files (`.opt.tour`) when available -4. Saves all files to the `dataset/` directory +This project includes a dataset setup script that downloads and filters TSP instances from TSPLIB95. To set up the dataset: ```bash @@ -74,53 +67,82 @@ The `uv run` commands will automatically handle the virtual environment for you. git clone cd python -m venv .venv - source .venv/bin/activate # On Windows: .venv\Scripts\activate + source .venv/bin/activate pip install -e . ``` 2. **Run commands**: ```bash python setup_dataset.py - python main.py dataset/.tsp + python tsp_solver.py dataset/.tsp ``` ## Usage To run the solver: + ``` -uv run python main.py +uv run python tsp_solver.py ``` This will output the total distance to stdout and generate `solution.csv` in the current directory. For development or testing, use `uv run python` to execute scripts in the project environment. -## Project Structure +## Generating Figures -The project is organized into several key packages and modules: +To generate performance figures for the TSP algorithms, you can run individual scripts from the `figure_scripts/` directory. -### Core Modules -- `main.py`: Main solver script and entry point. -- `setup_dataset.py`: Dataset setup script for downloading and filtering TSP instances from TSPLIB95. +```bash +uv run python -m figure_scripts.box_plot_figures +uv run python -m figure_scripts.relative_work_figures +uv run python -m figure_scripts.relative_work_nn_figures +uv run python -m figure_scripts.time_budget_figures +uv run python -m figure_scripts.time_budget_nn_figures +``` + +To generate all figures at once: + +```bash +uv run python generate_figures.py +``` -### Package Organization +Ensure the dataset is set up (run `uv run python setup_dataset.py` if not already done). -#### `tsp/` - TSP Core Package -- `model.py`: Core data structures (`City`, `TSPInstance`) and distance calculations. -- `io.py`: TSPLIB file parsing utilities for reading `.tsp` files. +## Hyperparam Tuning -#### `algorithm/` - Algorithm Implementations -- `base.py`: Protocol definitions and base classes for iterative TSP solvers. -- `nearest_neighbor.py`: Nearest neighbor constructive algorithm implementation. -- `random_solver.py`: Random permutation solver for baseline comparison. +The `tuning/` directory contains scripts for hyperparameter tuning of GA and SA. These tune parameters over a fixed time budget on the lin105.tsp instance. -### Data and Analysis -- `dataset/`: Directory containing TSP instances and optimal tour files (created by setup script). -- `bench_results/`: Directory for storing benchmark results. -- `tsp_analysis.ipynb`: Jupyter notebook for algorithm analysis and visualization. +```bash +uv run python -m tuning.ga_tuning +uv run python -m tuning.sa_tuning +``` + +The console output should include output of the best params. -### Configuration -- `pyproject.toml`: Project configuration and dependencies. -- `uv.lock`: Locked dependencies for reproducibility. -- `README.md`: This file. +## TSP Analysis Notebook +You can compile the notebook to PDF by running: + +``` +uv run jupyter nbconvert --to pdf tsp_analysis.ipynb +``` + +## Project Structure + +The project is organized into directories for core functionality, algorithms, data handling, and analysis: + +### Core Modules +- Entry point and utilities: `main.py`, `setup_dataset.py`, `generate_figures.py`, `constants.py`, `util.py`. + +### Packages +- `tsp/`: Core TSP model and I/O. +- `algorithm/`: Heuristic algorithm implementations (e.g., genetic, simulated annealing, nearest neighbor and random solver). +- `figure_scripts/`: Scripts for generating performance visualizations (e.g., box plots, time budgets, relative work comparisons). +- `tuning/`: Hyperparameter tuning scripts. +- `tests/`: Tests. + +### Data and Outputs +- `dataset/`: TSP instances and optimal tours. +- `figures/`: Generated plots. +- `solution.csv`: Solver output file. diff --git a/TODO.md b/TODO.md deleted file mode 100644 index 6c7e560..0000000 --- a/TODO.md +++ /dev/null @@ -1,9 +0,0 @@ -# TODOs - -- Implement one of the genetic algorithms seen in class -- Improve benchmarking framework (currently just the tsp_analysis.ipynb) - - Improve performance comparison charts - - Add optimal solution verification for problems that have it -- Set up a python script to take the tsp and output the tsp solution - - with maybe optional command line arguments for algorithm selection -- Create the PDF report diff --git a/TSP Analysis.ipynb b/TSP Analysis.ipynb new file mode 100644 index 0000000..8d30179 --- /dev/null +++ b/TSP Analysis.ipynb @@ -0,0 +1,365 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Traveling Salesman Problem: Algorithm Comparison\n", + "\n", + "This notebook compares six TSP algorithms on the Lin105 dataset with a 5-second time limit: Random Solver, Nearest Neighbor; Simulated Annealing and Genetic Algorithm (with random and NN initialization).\n", + "\n", + "We compare them by looking at their final cost versus the optimal and baseline, how their cost changes over time, how many steps they take per second, and how much they improve from the starting point." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "from pathlib import Path\n", + "from tsp.model import TSPInstance\n", + "from algorithm.nearest_neighbor import NearestNeighbor\n", + "from algorithm.simulated_annealing import SimulatedAnnealing\n", + "from algorithm.genetic_algo import GeneticAlgorithmSolver\n", + "from algorithm.random_solver import RandomSolver\n", + "from util import exponential_cooling, find_optimal_tour, run_algorithm_with_timing\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "# Config\n", + "MAX_SECONDS = 5.0\n", + "\n", + "# SA/GA parameters (from tuning run)\n", + "T0 = 167.807\n", + "COOLING_RATE = 0.99990\n", + "GA_POP_SIZE = 109\n", + "GA_CROSSOVER = 0.600\n", + "GA_MUTATION = 0.400\n", + "GA_ELITISM = 1\n", + "\n", + "# Cooling schedule\n", + "exp_schedule = exponential_cooling(COOLING_RATE)\n", + "\n", + "# Algorithm lineup\n", + "ALGORITHMS = [\n", + " \"Random Solver\",\n", + " \"Nearest Neighbor\",\n", + " \"SA (Random-init)\",\n", + " \"SA (NN-init)\",\n", + " \"GA (Random-init)\",\n", + " \"GA (NN-init)\",\n", + "]\n", + "\n", + "# Plot styling\n", + "plt.rcParams['figure.figsize'] = (10, 6)\n", + "plt.rcParams['font.size'] = 11" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "problem_instance_path = Path(\"dataset/lin105.tsp\")\n", + "instance, optimal_cost = find_optimal_tour(problem_instance_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "# Utils\n", + "def get_nn_route_and_cost(instance):\n", + " builder = NearestNeighbor(instance)\n", + " builder.initialize(None)\n", + " for _ in range(len(instance.cities) - 1):\n", + " builder.step()\n", + " return builder.get_route(), builder.get_cost()\n", + "\n", + "def run_single_time_trial(name, instance_data, seed_nn_data):\n", + " inst = TSPInstance(name=instance_data[\"name\"], cities=instance_data[\"cities\"])\n", + "\n", + " if name == \"Random Solver\":\n", + " solver, init_route = RandomSolver(inst), None\n", + " elif name == \"Nearest Neighbor\":\n", + " solver, init_route = NearestNeighbor(inst), None\n", + " elif name == \"SA (Random-init)\":\n", + " solver, init_route = SimulatedAnnealing(inst, T0, exp_schedule), None\n", + " elif name == \"SA (NN-init)\":\n", + " solver, init_route = SimulatedAnnealing(inst, T0, exp_schedule), seed_nn_data\n", + " elif name == \"GA (Random-init)\":\n", + " solver = GeneticAlgorithmSolver(\n", + " inst,\n", + " population_size=GA_POP_SIZE,\n", + " mutation_rate=GA_MUTATION,\n", + " crossover_rate=GA_CROSSOVER,\n", + " elitism_count=GA_ELITISM,\n", + " )\n", + " init_route = None\n", + " elif name == \"GA (NN-init)\":\n", + " solver = GeneticAlgorithmSolver(\n", + " inst,\n", + " population_size=GA_POP_SIZE,\n", + " mutation_rate=GA_MUTATION,\n", + " crossover_rate=GA_CROSSOVER,\n", + " elitism_count=GA_ELITISM,\n", + " )\n", + " init_route = seed_nn_data\n", + " else:\n", + " raise ValueError(f\"Unknown algorithm name: {name}\")\n", + "\n", + " iters, best, curr, times, route = run_algorithm_with_timing(inst, solver, init_route, MAX_SECONDS)\n", + " steps_per_sec = (len(iters) / times[-1]) if times else 0.0\n", + " return {\n", + " \"name\": name,\n", + " \"iterations\": iters,\n", + " \"best_costs\": best,\n", + " \"current_costs\": curr,\n", + " \"times\": times,\n", + " \"route\": route,\n", + " \"final_cost\": best[-1] if best else float(\"inf\"),\n", + " \"steps_per_sec\": steps_per_sec,\n", + " }\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "seed_nn, nn_cost = get_nn_route_and_cost(instance)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "# Run algorithms (for 30 secs total = n_algos * MAX_SECONDS)\n", + "instance_data = {\"name\": instance.name, \"cities\": instance.cities}\n", + "\n", + "time_runs = {name: [] for name in ALGORITHMS}\n", + "for name in ALGORITHMS:\n", + " time_runs[name] = [run_single_time_trial(name, instance_data, seed_nn)]\n", + "\n", + "time_results = {}\n", + "for name, runs in time_runs.items():\n", + " r = runs[0] if runs else {}\n", + " x = r.get(\"times\", [])\n", + " y = r.get(\"best_costs\", [])\n", + " time_results[name] = {\n", + " \"times\": x,\n", + " \"best\": np.array(y, dtype=float),\n", + " \"final_cost\": float(r.get(\"final_cost\", float('inf'))),\n", + " \"steps_per_sec\": float(r.get(\"steps_per_sec\", 0.0)),\n", + " }\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "summary_rows = []\n", + "for name, data in time_results.items():\n", + " cost = data[\"final_cost\"]\n", + " sps = data[\"steps_per_sec\"]\n", + " summary_rows.append({\n", + " \"Algorithm\": name,\n", + " \"Final Cost\": f\"{cost:.1f}\",\n", + " \"Steps/sec\": f\"{sps:.1f}\",\n", + " })\n", + "\n", + "summary_rows.append({\n", + " \"Algorithm\": \"Optimal (ref)\",\n", + " \"Final Cost\": f\"{optimal_cost:.1f}\",\n", + " \"Steps/sec\": \"\",\n", + "})\n", + "\n", + "pd.DataFrame(summary_rows).sort_values(\"Final Cost\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "for name, data in time_results.items():\n", + " ax.plot(data[\"times\"], data[\"best\"], label=name, linewidth=2)\n", + "ax.axhline(y=optimal_cost, color=\"green\", linestyle=\":\", label=\"Optimal\", alpha=0.7)\n", + "ax.set_xlabel(\"Time (s)\")\n", + "ax.set_ylabel(\"Best cost\")\n", + "ax.set_title('Best Cost Over Time')\n", + "ax.set_xscale('log')\n", + "ax.grid(True, alpha=0.3)\n", + "ax.legend()\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_best_over_time_for(names_subset, title):\n", + " fig, ax = plt.subplots()\n", + " for n in names_subset:\n", + " if n in time_results:\n", + " ax.plot(time_results[n][\"times\"], time_results[n][\"best\"], label=n, linewidth=2)\n", + " if optimal_cost:\n", + " ax.axhline(y=optimal_cost, color=\"green\", linestyle=\":\", label=\"Optimal\", alpha=0.7)\n", + " ax.set_xlabel(\"Time (s)\")\n", + " ax.set_ylabel(\"Best cost\")\n", + " ax.set_title(title)\n", + " ax.set_xscale('log')\n", + " ax.grid(True, alpha=0.3)\n", + " ax.legend()\n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + "plot_best_over_time_for([\n", + " \"Random Solver\", \"Nearest Neighbor\"\n", + "], \"Best Cost Over Time — Random vs Nearest Neighbor\")\n", + "plot_best_over_time_for([\n", + " \"SA (Random-init)\", \"GA (Random-init)\"\n", + "], \"Best Cost Over Time — SA (Random-init) vs GA (Random-init)\")\n", + "plot_best_over_time_for([\n", + " \"SA (NN-init)\", \"GA (NN-init)\", \"Nearest Neighbor\"\n", + "], \"Best Cost Over Time — SA (NN-init) vs GA (NN-init) vs Nearest Neighbor\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "names = [n for n in time_runs.keys() if (not n.lower().startswith('random')) and (not n.lower().startswith('nearest'))]\n", + "\n", + "final_costs_by_algo = {n: [] for n in names}\n", + "for n in names:\n", + " runs = time_runs.get(n, [])\n", + " for r in runs:\n", + " best_costs = r.get(\"best_costs\", [])\n", + " if best_costs:\n", + " final_costs_by_algo[n].append(best_costs[-1])\n", + "\n", + "gaps_by_algo = {}\n", + "for n, costs in final_costs_by_algo.items():\n", + " gaps = [((c / optimal_cost) - 1) * 100.0 for c in costs if c > 0]\n", + " gaps_by_algo[n] = gaps\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 5))\n", + "\n", + "for i, n in enumerate(names):\n", + " y = gaps_by_algo[n]\n", + " x = np.full(len(y), i)\n", + " ax.scatter(x, y, s=80, label=n, alpha=0.8)\n", + "\n", + "ax.set_xticks(range(len(names)))\n", + "ax.set_xticklabels(names, rotation=20)\n", + "ax.set_ylabel('Gap to Optimal (%)')\n", + "ax.set_title('Final Cost Gap')\n", + "ax.grid(True, axis='y', alpha=0.3)\n", + "\n", + "nn_gap = ((nn_cost / optimal_cost) - 1) * 100.0\n", + "opt_line = ax.axhline(y=0, color='green', linestyle=':', label='Optimal')\n", + "nn_line = ax.axhline(y=nn_gap, color='orange', linestyle='--', label='NN')\n", + "ax.legend(handles=[opt_line, nn_line])\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_relative_improvement_over_time(names_subset, title=None):\n", + " fig, ax = plt.subplots()\n", + " for n in names_subset:\n", + " if n in time_results:\n", + " times = time_results[n]['times'] \n", + " best = time_results[n]['best']\n", + " if len(best) > 0 and best[0] > 0:\n", + " initial = float(best[0])\n", + " improvement_pct = (initial - best) / initial * 100.0\n", + " ax.plot(times, improvement_pct, label=n, linewidth=2)\n", + " ax.set_xlabel('Time (s)')\n", + " ax.set_ylabel('Improvement from Initial (%)')\n", + " ax.set_title(title)\n", + " ax.set_xscale('log')\n", + " ax.legend()\n", + " ax.grid(True, which='both', axis='x', alpha=0.3)\n", + " ax.grid(True, which='major', axis='y', alpha=0.3)\n", + " plt.show()\n", + "\n", + "plot_relative_improvement_over_time(\n", + " [\"Random Solver\", \"Nearest Neighbor\"],\n", + " \"Improvement from Initial (%) — Random vs Nearest Neighbor\"\n", + ")\n", + "plot_relative_improvement_over_time(\n", + " [\"SA (Random-init)\", \"GA (Random-init)\"],\n", + " \"Improvement from Initial (%) — SA (Random-init) vs GA (Random-init)\"\n", + ")\n", + "plot_relative_improvement_over_time(\n", + " [\"SA (NN-init)\", \"GA (NN-init)\", \"Nearest Neighbor\"],\n", + " \"Improvement from Initial (%) — SA (NN-init) vs GA (NN-init) vs Nearest Neighbor\"\n", + ")\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/TSP Analysis.pdf b/TSP Analysis.pdf new file mode 100644 index 0000000..0cbfff1 Binary files /dev/null and b/TSP Analysis.pdf differ diff --git a/algorithm/genetic_algo.py b/algorithm/genetic_algo.py index 18cd3ef..d5f00f1 100644 --- a/algorithm/genetic_algo.py +++ b/algorithm/genetic_algo.py @@ -11,7 +11,7 @@ class GeneticAlgorithmSolver(IterativeTSPSolver): """Genetic Algorithm for TSP.""" - def __init__(self, instance: TSPInstance, seed: int | float | None = None, population_size: int = 100, crossover_rate: float = 0.7, mutation_rate: float = 0.01, elitism_count: int = 2, num_parents: int = 2, num_child: int = 2): + def __init__(self, instance: TSPInstance, seed: int | float | None = None, population_size: int = 100, crossover_rate: float = 0.7, mutation_rate: float = 0.01, elitism_count: int = 2): self.instance = instance self.rng = random.Random(seed) self.best_route: List[int] = [] @@ -21,8 +21,6 @@ def __init__(self, instance: TSPInstance, seed: int | float | None = None, popul self.population: List[List[int]] = [] self.mutation_rate = mutation_rate self.crossover_rate = crossover_rate - self.num_parents = num_parents - self.num_child = num_child self.elitism_count = elitism_count self.fitness_cache = OrderedDict() self.max_cache_size = 2 * population_size @@ -79,7 +77,7 @@ def _get_fitness(self, route: List[int]) -> float: def step(self) -> StepReport: self.iteration += 1 - # Evaluate fitness once per individual + # Evaluate fitness for each individual in the population population_with_fitness = [(ind, self._get_fitness(ind)) for ind in self.population] improved = False @@ -93,27 +91,26 @@ def step(self) -> StepReport: average_fitness = fitness_sum / len(self.population) - # Extract elite from already-sorted population + + new_population = [] + + # Carry best individuals to the next generation if self.elitism_count > 0: population_with_fitness.sort(key=lambda x: x[1]) elite = [ind for ind, _ in population_with_fitness[:self.elitism_count]] - else: - elite = [] + new_population.extend(elite) - new_population = elite.copy() + # Selection, crossover, and mutation steps while len(new_population) < self.population_size: - selected_parents = [] - for _ in range(self.num_parents): - parent = self.select_parent() - selected_parents.append(parent) + parent1 = self.select_parent() + parent2 = self.select_parent() - for _ in range(self.num_child): - if len(new_population) >= self.population_size: - break - child = self.crossover(selected_parents) - if self.rng.random() < self.mutation_rate: - child = self.mutate(child) - new_population.append(child) + child = self.crossover(parent1, parent2) + + if self.rng.random() < self.mutation_rate: + child = self.mutate(child) + + new_population.append(child) self.population = new_population @@ -126,32 +123,23 @@ def get_cost(self) -> float: return self.best_cost def select_parent(self) -> List[int]: - """Tournament selection with cached fitness.""" + """Tournament selection.""" tournament_size = max(2, self.population_size // 10) tournament = self.rng.sample(self.population, tournament_size) tournament.sort(key=lambda route: self._get_fitness(route)) return tournament[0] - def crossover(self, parents: List[List[int]]) -> List[int]: - """Ordered Crossover (OX) for TSP. Returns a new child.""" - # Get two parents - if len(parents) < 2: - return parents[0][:] - - parent1, parent2 = self.rng.sample(parents, 2) - + def crossover(self, parent1: List[int], parent2: List[int]) -> List[int]: + """Ordered Crossover (OX) for TSP.""" if self.rng.random() > self.crossover_rate: return self.rng.choice([parent1, parent2])[:] - # Choose random crossover points n = len(parent1) start, end = sorted(self.rng.sample(range(n), 2)) - # Copy a random segment from parent1 child = [None] * n child[start:end] = parent1[start:end] - # Fill remaining positions with genes from parent2 in order fill_pos = 0 for gene in parent2: if gene not in child: @@ -159,7 +147,6 @@ def crossover(self, parents: List[List[int]]) -> List[int]: fill_pos += 1 child[fill_pos] = gene - assert None not in child return child def mutate(self, route: List[int]) -> List[int]: @@ -167,6 +154,5 @@ def mutate(self, route: List[int]) -> List[int]: mutated_route = route[:] n = len(mutated_route) i, j = sorted(self.rng.sample(range(n), 2)) - # Reverse the segment between i and j (2-opt move) mutated_route[i:j+1] = reversed(mutated_route[i:j+1]) return mutated_route \ No newline at end of file diff --git a/algorithm/nearest_neighbor.py b/algorithm/nearest_neighbor.py index 870cc95..09a313d 100644 --- a/algorithm/nearest_neighbor.py +++ b/algorithm/nearest_neighbor.py @@ -8,49 +8,62 @@ class NearestNeighbor(IterativeTSPSolver): - """Constructive iterative nearest neighbor. - - Each step adds one more nearest unvisited city to the route until complete. - """ + """Nearest neighbor sampling.""" def __init__(self, instance: TSPInstance, seed: int | float | None = None): self.instance = instance self.rng = random.Random(seed) - self.route: List[int] = [] - self.unvisited: List[int] = [] self.iteration = 0 - def initialize(self, route: list[int] | None = None) -> None: + self.best_route: List[int] = [] + self.best_cost: float = float("inf") + + def _build_nn_route(self, start_city: int) -> List[int]: n = len(self.instance.cities) + route: List[int] = [start_city] + unvisited: List[int] = [i for i in range(n) if i != start_city] + while unvisited: + last = route[-1] + next_city = min(unvisited, key=lambda j: self.instance.distance(last, j)) + unvisited.remove(next_city) + route.append(next_city) + return route + + def initialize(self, route: List[int] | None = None) -> None: + n = len(self.instance.cities) + self.iteration = 0 if route and len(route) > 0: - self.route = [route[0]] - self.unvisited = [i for i in range(n) if i not in self.route] + start_city = route[0] else: start_city = self.rng.randint(0, n - 1) - self.route = [start_city] - self.unvisited = [i for i in range(n) if i != start_city] - self.iteration = 0 + candidate = self._build_nn_route(start_city) + self.best_route = candidate + self.best_cost = self.instance.route_cost(candidate) def step(self) -> StepReport: self.iteration += 1 + n = len(self.instance.cities) + start_city = self.rng.randint(0, n - 1) + candidate = self._build_nn_route(start_city) + current_cost = self.instance.route_cost(candidate) + improved = False - if self.unvisited: - last = self.route[-1] # FIFO - # Pick nearest unvisited - next_city = min(self.unvisited, key=lambda j: self.instance.distance(last, j)) - self.unvisited.remove(next_city) - self.route.append(next_city) + if current_cost < self.best_cost: + self.best_cost = current_cost + self.best_route = candidate improved = True - current_cost = self.get_cost() - return StepReport(iteration=self.iteration, best_cost=current_cost, current_cost=current_cost, improved=improved) + + return StepReport( + iteration=self.iteration, + best_cost=self.best_cost, + current_cost=current_cost, + improved=improved, + ) def get_route(self) -> List[int]: - # If incomplete, return current partial route followed by unvisited - return self.route + self.unvisited + return self.best_route def get_cost(self) -> float: - if len(self.route) < 2: - return 0.0 - return self.instance.route_cost(self.get_route()) + return self.best_cost diff --git a/algorithm/random_solver.py b/algorithm/random_solver.py index 1d9a213..24a5285 100644 --- a/algorithm/random_solver.py +++ b/algorithm/random_solver.py @@ -8,7 +8,7 @@ class RandomSolver(IterativeTSPSolver): - """Each iteration samples a new random permutation; keeps track of the best seen.""" + """Random permutation sampling.""" def __init__(self, instance: TSPInstance, seed: int | float | None = None): self.instance = instance diff --git a/constants.py b/constants.py index 4a11ec5..fd62483 100644 --- a/constants.py +++ b/constants.py @@ -1,13 +1,24 @@ -# Configuration Constants -MAX_ITERATIONS = 1_000 # Number of iterations for iteration-based benchmark -MAX_SECONDS = 1.0 # Number of seconds for time-based benchmark -RANDOM_SEED = 42 # Random seed for reproducible results -COOLING_RATE = 0.995 # Cooling rate for Simulated Annealing -INITIAL_TEMP = 934.622 # Initial temperature for Simulated Annealing -POPULATION_SIZE = 36 # Population size for Genetic Algorithm -MUTATION_RATE = 0.300 # Mutation rate for Genetic Algorithm -CROSSOVER_RATE = 0.600 # Crossover rate for Genetic Algorithm -ELITISM_COUNT = 5 # Number of elite individuals to carry over in Genetic Algorithm -NUM_PARENTS = 6 # Number of parents for Genetic Algorithm -NUM_CHILD = 4 # Number of children per crossover for Genetic Algorithm -N_RUNS = 10 # Number of runs for averaging results in parallel benchmarking \ No newline at end of file +# General Settings +MAX_ITERATIONS = 1_000 # Number of iterations for relative iteration-based benchmark +MAX_SECONDS = 5.0 # Number of seconds for time-based benchmark +N_RUNS = 5 # Number of runs for averaging results in benchmarking + +DATASET_FILENAME = 'lin105.tsp' + +# Parallelization Settings +PARALLEL_RUNS = True +NUM_WORKERS = 8 # Number of parallel workers for benchmarking trials + +# Calibration Constants +CALIBRATION_TIME = 2.0 +MAX_NORMALIZED_STEPS = 100_000 + +# Simulated Annealing (SA) Settings +COOLING_RATE = 0.99990 +INITIAL_TEMP = 167.807 + +# Genetic Algorithm (GA) Settings +POPULATION_SIZE = 109 +MUTATION_RATE = 0.400 +CROSSOVER_RATE = 0.600 +ELITISM_COUNT = 1 \ No newline at end of file diff --git a/figure_scripts/__init__.py b/figure_scripts/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/figure_scripts/box_plot_figures.py b/figure_scripts/box_plot_figures.py new file mode 100644 index 0000000..86ec841 --- /dev/null +++ b/figure_scripts/box_plot_figures.py @@ -0,0 +1,87 @@ +import numpy as np +import matplotlib.pyplot as plt +import logging +from constants import MAX_SECONDS, N_RUNS, PARALLEL_RUNS +from .common import ( + load_tsp_instance, create_solvers, create_plot, get_nn_initial_route, + compute_nn_baseline, run_parallel_trials, save_figure, create_box_plot_statistics, save_statistics_json +) +from util import run_algorithm_with_timing +from tsp.model import TSPInstance + +logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') +logger = logging.getLogger(__name__) + +def run_single_box_trial(args): + name, instance_data, use_nn = args + solvers = create_solvers() + solver_factory = solvers.get(name.replace('_NN', '_random') if use_nn else name) + if solver_factory: + solver = solver_factory() + init_route = None + if use_nn: + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + init_route = get_nn_initial_route(instance) + else: + raise ValueError(f"Unknown algorithm: {name}") + + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + _, best_costs, _, _, _ = run_algorithm_with_timing( + instance, solver, init_route, MAX_SECONDS + ) + return best_costs[-1] if best_costs else float('inf') + +def main(): + _, optimal_cost, instance_data = load_tsp_instance() + logger.info(f"Generating box plots for lin105 (optimal: {optimal_cost:.2f})") + + algorithms = ["SA_random", "GA_random", "SA_NN", "GA_NN"] + + args_list = [] + for name in algorithms: + use_nn = name.endswith('_NN') + for _ in range(N_RUNS): + args_list.append((name, instance_data, use_nn)) + + logger.info(f"Running {len(args_list)} trials {'in parallel' if PARALLEL_RUNS else 'sequentially'}") + + all_final_costs = run_parallel_trials(run_single_box_trial, args_list, desc="Running trials") + + # Group by algorithm + costs_by_algo = {name: [] for name in algorithms} + for i, cost in enumerate(all_final_costs): + algo_idx = i // N_RUNS + algo_name = algorithms[algo_idx] + costs_by_algo[algo_name].append(cost) + + # Compute gaps + gaps_by_algo = {} + for algo, costs in costs_by_algo.items(): + gaps = [(c / optimal_cost - 1) * 100 for c in costs if c > 0] + gaps_by_algo[algo] = gaps + mean_gap = np.mean(gaps) + std_gap = np.std(gaps) + logger.info(f"{algo}: Mean gap {mean_gap:.1f}% ± {std_gap:.1f}%") + + # Compute NN baseline + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + nn_cost = compute_nn_baseline(instance) + nn_gap = ((nn_cost / optimal_cost - 1) * 100) if optimal_cost else 0 + logger.info("Nearest Neighbor") + + # Create and save statistics + statistics = create_box_plot_statistics(instance, optimal_cost, instance_data, costs_by_algo, algorithms) + save_statistics_json(statistics, 'box_plot_figures.json') + + # Plot box plot + _, ax = create_plot(f'Algorithm Performance Distribution (Final Costs after {MAX_SECONDS}s)', 'Algorithms', 'Gap to Optimal (%)') + ax.boxplot([gaps_by_algo[name] for name in algorithms], tick_labels=algorithms) + ax.axhline(y=0, color='green', linestyle=':', label='Optimal') + ax.axhline(y=nn_gap, color='orange', linestyle='--', label='NN') + ax.legend() + + save_figure(plt.gcf(), 'figures/box_plot_figures.png') + logger.info("Saved figures/box_plot_figures.png") + +if __name__ == "__main__": + main() diff --git a/figure_scripts/common.py b/figure_scripts/common.py new file mode 100644 index 0000000..256d59b --- /dev/null +++ b/figure_scripts/common.py @@ -0,0 +1,221 @@ +import matplotlib.pyplot as plt +from pathlib import Path +import json +from datetime import datetime +import numpy as np +from concurrent.futures import ProcessPoolExecutor, as_completed +from tqdm import tqdm + +from constants import INITIAL_TEMP, COOLING_RATE, POPULATION_SIZE, MUTATION_RATE, CROSSOVER_RATE, ELITISM_COUNT, DATASET_FILENAME, PARALLEL_RUNS, NUM_WORKERS, MAX_SECONDS, N_RUNS, CALIBRATION_TIME, MAX_NORMALIZED_STEPS +from util import find_optimal_tour, exponential_cooling +from algorithm.simulated_annealing import SimulatedAnnealing +from algorithm.genetic_algo import GeneticAlgorithmSolver +from algorithm.nearest_neighbor import NearestNeighbor + + +def load_tsp_instance(tsp_filename=DATASET_FILENAME): + tsp_path = Path('dataset') / tsp_filename + instance, optimal_cost = find_optimal_tour(tsp_path) + instance_data = {'name': instance.name, 'cities': instance.cities} + return instance, optimal_cost, instance_data + + +def create_solvers(): + def create_sa_solver(): + instance, _, _ = load_tsp_instance() + T0 = INITIAL_TEMP + cool_rate = COOLING_RATE + schedule = exponential_cooling(cool_rate) + return SimulatedAnnealing(instance, T0, schedule) + + def create_ga_solver(): + instance, _, _ = load_tsp_instance() + return GeneticAlgorithmSolver( + instance, + population_size=POPULATION_SIZE, + mutation_rate=MUTATION_RATE, + crossover_rate=CROSSOVER_RATE, + elitism_count=ELITISM_COUNT + ) + + return { + 'SA_random': create_sa_solver, + 'GA_random': create_ga_solver + } + + +def create_plot(title, xlabel, ylabel, figsize=(10, 6)): + fig, ax = plt.subplots(figsize=figsize) + # Use /caption in LaTeX instead + # ax.set_title(title) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.grid(True, alpha=0.3) + return fig, ax + +def get_nn_initial_route(instance): + """Compute the initial route using Nearest Neighbor heuristic.""" + nn = NearestNeighbor(instance) + nn.initialize(None) + n_cities = len(instance.cities) + for _ in range(n_cities - 1): + nn.step() + return nn.get_route() + +def compute_nn_baseline(instance): + """Compute the baseline cost using Nearest Neighbor.""" + nn = NearestNeighbor(instance) + nn.initialize(None) + n_cities = len(instance.cities) + for _ in range(n_cities - 1): + nn.step() + return nn.get_cost() + +def align_series(x_lists, y_lists, common_x): + """Align multiple runs' series data via interpolation and compute mean/std.""" + aligned_y = [] + for x, y in zip(x_lists, y_lists): + if len(x) > 0 and len(y) > 0: + interp_y = np.interp(common_x, x, y) + aligned_y.append(interp_y) + if not aligned_y: + return np.array([]), np.array([]) + aligned_y = np.array(aligned_y) + return np.mean(aligned_y, axis=0), np.std(aligned_y, axis=0) + +def save_figure(fig, filepath): + """Save and show the figure.""" + plt.tight_layout() + plt.savefig(filepath, dpi=150, bbox_inches='tight') + plt.show() + +def add_optimal_line(ax, optimal_cost): + """Add the optimal cost horizontal line to the axis.""" + ax.axhline(y=optimal_cost, color='green', linestyle=':', label='Optimal', alpha=0.7) + +ALGO_COLORS = { + 'SA': 'blue', + 'GA': 'red' +} + +ALGO_LINESTYLES = { + 'SA': '-', + 'GA': '--' +} + +def run_parallel_trials(trial_func, args_list, parallel=PARALLEL_RUNS, max_workers=NUM_WORKERS, desc="Running trials"): + """Run trials in parallel or sequentially.""" + if parallel: + with ProcessPoolExecutor(max_workers=max_workers) as executor: + futures = [executor.submit(trial_func, arg) for arg in args_list] + return [f.result() for f in tqdm(as_completed(futures), total=len(args_list), desc=desc)] + else: + return [trial_func(arg) for arg in tqdm(args_list, desc=desc)] + + +def compute_statistics(data, optimal_cost=None): + """Compute comprehensive statistics for a dataset.""" + if not data or len(data) == 0: + return {} + + data = np.array(data) + stats = { + 'count': len(data), + 'mean': float(np.mean(data)), + 'std': float(np.std(data)), + 'min': float(np.min(data)), + 'max': float(np.max(data)), + 'median': float(np.median(data)), + 'q1': float(np.percentile(data, 25)), + 'q3': float(np.percentile(data, 75)), + 'iqr': float(np.percentile(data, 75) - np.percentile(data, 25)) + } + + if optimal_cost and optimal_cost > 0: + stats['gap_to_optimal_pct'] = float((stats['mean'] / optimal_cost - 1) * 100) + stats['best_gap_to_optimal_pct'] = float((stats['min'] / optimal_cost - 1) * 100) + stats['worst_gap_to_optimal_pct'] = float((stats['max'] / optimal_cost - 1) * 100) + + return stats + + +def create_metadata(instance_data, optimal_cost, experiment_type, additional_settings=None): + """Create standardized metadata for statistics JSON files.""" + base_settings = { + 'max_seconds': MAX_SECONDS, + 'n_runs': N_RUNS, + 'parallel_runs': PARALLEL_RUNS, + 'num_workers': NUM_WORKERS + } + + if additional_settings: + base_settings.update(additional_settings) + + return { + 'generated_at': datetime.now().isoformat(), + 'tsp_instance': instance_data['name'], + 'optimal_cost': optimal_cost, + 'experiment_settings': base_settings + } + + +def save_statistics_json(statistics, filename): + """Save statistics to JSON file with proper formatting.""" + with open(f'figures/{filename}', 'w') as f: + json.dump(statistics, f, indent=2) + print(f"Saved figures/{filename}") + + +def create_box_plot_statistics(instance, optimal_cost, instance_data, costs_by_algo, algorithms): + """Create statistics structure for box plot experiments.""" + # Compute statistics for each algorithm + stats_by_algo = {} + for algo, costs in costs_by_algo.items(): + stats_by_algo[algo] = compute_statistics(costs, optimal_cost) + + # Compute NN baseline statistics + nn_cost = compute_nn_baseline(instance) + nn_stats = compute_statistics([nn_cost], optimal_cost) + + return { + 'metadata': create_metadata(instance_data, optimal_cost, 'box_plot'), + 'algorithms': stats_by_algo, + 'nearest_neighbor_baseline': nn_stats, + 'experiment_config': { + 'max_seconds': MAX_SECONDS, + 'n_runs': N_RUNS, + 'algorithms': algorithms + } + } + + +def create_time_budget_statistics(instance_data, optimal_cost, stats_by_algo, algorithms, use_nn=False): + """Create statistics structure for time budget experiments.""" + return { + 'metadata': create_metadata(instance_data, optimal_cost, 'time_budget'), + 'algorithms': stats_by_algo, + 'experiment_config': { + 'max_seconds': MAX_SECONDS, + 'n_runs': N_RUNS, + 'algorithms': algorithms, + 'use_nn_initialization': use_nn + } + } + + +def create_relative_work_statistics(instance_data, optimal_cost, stats_by_algo, algorithms, calibration_data, use_nn=False): + """Create statistics structure for relative work experiments.""" + return { + 'metadata': create_metadata(instance_data, optimal_cost, 'relative_work', { + 'calibration_time': CALIBRATION_TIME, + 'max_normalized_steps': MAX_NORMALIZED_STEPS + }), + 'algorithms': stats_by_algo, + 'calibration': calibration_data, + 'experiment_config': { + 'max_normalized_steps': MAX_NORMALIZED_STEPS, + 'n_runs': N_RUNS, + 'algorithms': algorithms, + 'use_nn_initialization': use_nn + } + } diff --git a/figure_scripts/random_baseline_figures.py b/figure_scripts/random_baseline_figures.py new file mode 100644 index 0000000..7be5b58 --- /dev/null +++ b/figure_scripts/random_baseline_figures.py @@ -0,0 +1,142 @@ +import numpy as np +import matplotlib.pyplot as plt +import logging +from constants import MAX_SECONDS, N_RUNS, PARALLEL_RUNS +from .common import ( + load_tsp_instance, create_solvers, create_plot, get_nn_initial_route, + run_parallel_trials, save_figure, create_box_plot_statistics, save_statistics_json +) +from util import run_algorithm_with_timing +from tsp.model import TSPInstance +from algorithm.random_solver import RandomSolver +from algorithm.nearest_neighbor import NearestNeighbor + +logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') +logger = logging.getLogger(__name__) + +def run_single_baseline_trial(args): + """Single trial for random baseline comparison.""" + name, instance_data, use_nn = args + solvers = create_solvers() + + if name == 'Random': + # Use random solver + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + solver = RandomSolver(instance) + init_route = None + else: + # Use existing solvers + solver_factory = solvers.get(name.replace('_NN', '_random') if use_nn else name) + if solver_factory: + solver = solver_factory() + init_route = None + if use_nn: + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + init_route = get_nn_initial_route(instance) + else: + raise ValueError(f"Unknown algorithm: {name}") + + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + _, best_costs, _, _, _ = run_algorithm_with_timing( + instance, solver, init_route, MAX_SECONDS + ) + return best_costs[-1] if best_costs else float('inf') + +def main(): + _, optimal_cost, instance_data = load_tsp_instance() + logger.info(f"Generating random baseline comparison for {instance_data['name']} (optimal: {optimal_cost:.2f})") + + algorithms = ["Random"] + + args_list = [] + for name in algorithms: + use_nn = name.endswith('_NN') + for _ in range(N_RUNS): + args_list.append((name, instance_data, use_nn)) + + logger.info(f"Running {len(args_list)} trials {'in parallel' if PARALLEL_RUNS else 'sequentially'}") + + all_final_costs = run_parallel_trials(run_single_baseline_trial, args_list, desc="Running trials") + + # Group by algorithm + costs_by_algo = {name: [] for name in algorithms} + for i, cost in enumerate(all_final_costs): + algo_idx = i // N_RUNS + algo_name = algorithms[algo_idx] + costs_by_algo[algo_name].append(cost) + + # Compute gaps + gaps_by_algo = {} + for algo, costs in costs_by_algo.items(): + gaps = [(c / optimal_cost - 1) * 100 for c in costs if c > 0] + gaps_by_algo[algo] = gaps + mean_gap = np.mean(gaps) + std_gap = np.std(gaps) + logger.info(f"{algo}: Mean gap {mean_gap:.1f}% ± {std_gap:.1f}%") + + # Compute NN baseline with random starts + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + nn_costs = [] + for _ in range(N_RUNS): + nn = NearestNeighbor(instance, seed=None) # Different seed each time + nn.initialize(None) + n_cities = len(instance.cities) + for _ in range(n_cities - 1): + nn.step() + nn_costs.append(nn.get_cost()) + + nn_gaps = [(c / optimal_cost - 1) * 100 for c in nn_costs if c > 0] + mean_nn_gap = np.mean(nn_gaps) + std_nn_gap = np.std(nn_gaps) + logger.info(f"Nearest Neighbor: Mean gap {mean_nn_gap:.1f}% ± {std_nn_gap:.1f}%") + + # Add NN to the algorithms for plotting + algorithms_with_nn = algorithms + ["Nearest Neighbor"] + costs_by_algo["Nearest Neighbor"] = nn_costs + gaps_by_algo["Nearest Neighbor"] = nn_gaps + + # Create and save statistics + statistics = create_box_plot_statistics(instance, optimal_cost, instance_data, costs_by_algo, algorithms_with_nn) + save_statistics_json(statistics, 'random_baseline_figures.json') + + # Plot random baseline + fig, ax = create_plot( + f'Random Baseline Performance Distribution (Final Costs after {MAX_SECONDS}s)', + 'Algorithms', + 'Gap to Optimal (%)', + figsize=(5, 6) + ) + + # Create box plot with custom colors + box_colors = ['red', 'orange'] # Random, Nearest Neighbor + box_plot = ax.boxplot([gaps_by_algo[name] for name in algorithms_with_nn], + tick_labels=algorithms_with_nn, patch_artist=True) + + # Color the boxes + for patch, color in zip(box_plot['boxes'], box_colors): + patch.set_facecolor(color) + patch.set_alpha(0.7) + + # Add reference lines + ax.axhline(y=0, color='green', linestyle=':', label='Optimal', linewidth=2) + + # Add mean lines for comparison + random_gaps = gaps_by_algo['Random'] + random_mean = np.mean(random_gaps) + ax.axhline(y=random_mean, color='red', linestyle='-', alpha=0.5, label='Random', linewidth=2) + + nn_gaps = gaps_by_algo['Nearest Neighbor'] + nn_mean = np.mean(nn_gaps) + ax.axhline(y=nn_mean, color='orange', linestyle='-', alpha=0.5, label='NN', linewidth=2) + + ax.legend() + ax.grid(True, alpha=0.3) + + # Rotate x-axis labels for better readability + plt.xticks(rotation=45) + + save_figure(fig, 'figures/random_baseline_figures.png') + logger.info("Saved figures/random_baseline_figures.png") + +if __name__ == "__main__": + main() diff --git a/figure_scripts/relative_work_figures.py b/figure_scripts/relative_work_figures.py new file mode 100644 index 0000000..947bb02 --- /dev/null +++ b/figure_scripts/relative_work_figures.py @@ -0,0 +1,217 @@ +import time +import math +import logging +from tqdm import tqdm +import numpy as np +import json +from datetime import datetime +from constants import CALIBRATION_TIME, MAX_NORMALIZED_STEPS, N_RUNS, PARALLEL_RUNS, NUM_WORKERS +from .common import load_tsp_instance, create_solvers, create_plot, align_series, save_figure, ALGO_COLORS, add_optimal_line +from concurrent.futures import ProcessPoolExecutor, as_completed + +logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') +logger = logging.getLogger(__name__) + +def compute_statistics(data, optimal_cost=None): + """Compute comprehensive statistics for a dataset.""" + if not data or len(data) == 0: + return {} + + data = np.array(data) + stats = { + 'count': len(data), + 'mean': float(np.mean(data)), + 'std': float(np.std(data)), + 'min': float(np.min(data)), + 'max': float(np.max(data)), + 'median': float(np.median(data)), + 'q1': float(np.percentile(data, 25)), + 'q3': float(np.percentile(data, 75)), + 'iqr': float(np.percentile(data, 75) - np.percentile(data, 25)) + } + + if optimal_cost and optimal_cost > 0: + stats['gap_to_optimal_pct'] = float((stats['mean'] / optimal_cost - 1) * 100) + stats['best_gap_to_optimal_pct'] = float((stats['min'] / optimal_cost - 1) * 100) + stats['worst_gap_to_optimal_pct'] = float((stats['max'] / optimal_cost - 1) * 100) + + return stats + + + +def calibrate_steps_per_second(create_solver_func, initial_route=None, calibration_time=CALIBRATION_TIME): + solver = create_solver_func() + solver.initialize(initial_route) + + start_time = time.perf_counter() + step_count = 0 + + while time.perf_counter() - start_time < calibration_time: + solver.step() + step_count += 1 + + actual_time = time.perf_counter() - start_time + return step_count / actual_time + + +def run_single_benchmark(create_solver_func, initial_route, work_factor, max_normalized_steps=MAX_NORMALIZED_STEPS): + solver = create_solver_func() + solver.initialize(initial_route) + + max_actual_steps = max(1, math.floor(max_normalized_steps / work_factor)) + results = [] + for step in range(max_actual_steps): + solver.step() + normalized_step = (step + 1) * work_factor + fitness = solver.get_cost() # Get current best cost (fitness) + results.append((normalized_step, fitness)) + + return results + + +def worker_benchmark(algo_name, work_factor): + """Worker function for parallel benchmarking - imports locally to avoid pickling.""" + from figure_scripts.relative_work_figures import run_single_benchmark + from .common import create_solvers + solvers = create_solvers() + create_func = solvers[algo_name] + result = run_single_benchmark(create_func, None, work_factor) + return (algo_name, result) + + +def main(): + _, optimal_cost, instance_data = load_tsp_instance() + logger.info(f"Generating relative work figures for {instance_data['name']} (optimal: {optimal_cost:.2f})") + + # Calibration phase (single run for speed) + solvers = create_solvers() + steps_per_second = {} + for algo_name, create_func in solvers.items(): + steps_per_second[algo_name] = calibrate_steps_per_second(create_func) + + logger.info(f"Steps per second: {steps_per_second}") + + # Find the algorithm with the most steps per second (least work per step) + reference_algo = max(steps_per_second, key=steps_per_second.get) + reference_steps_per_second = steps_per_second[reference_algo] + + logger.info(f"Reference algorithm (fastest): {reference_algo} with {reference_steps_per_second:.1f} steps/sec") + + # Calculate normalization factors (work per step relative to reference algorithm) + work_per_step = {} + for algo_name in solvers: + work_per_step[algo_name] = reference_steps_per_second / steps_per_second[algo_name] + + logger.info(f"Work per step: {work_per_step}") + + args_list = [] + for algo_name in solvers: + for _ in range(N_RUNS): + args_list.append((algo_name, work_per_step[algo_name])) + + logger.info(f"Running {len(args_list)} trials {'in parallel' if PARALLEL_RUNS else 'sequentially'}") + + if PARALLEL_RUNS: + adjusted_args = [(algo_name, work_per_step[algo_name]) for algo_name in solvers for _ in range(N_RUNS)] + with ProcessPoolExecutor(max_workers=NUM_WORKERS) as executor: + futures = [executor.submit(worker_benchmark, *arg) for arg in adjusted_args] + all_results_list = [f.result() for f in tqdm(as_completed(futures), total=len(adjusted_args), desc="Running trials")] + else: + all_results_list = [] + for arg in tqdm(args_list, desc="Running trials"): + algo_name, work_factor = arg + create_func = solvers[algo_name] + result = run_single_benchmark(create_func, None, work_factor) + all_results_list.append((algo_name, result)) + + all_results = {name: [] for name in solvers} + for algo_name, result in all_results_list: + all_results[algo_name].append(result) + + # Compute statistics for final costs + stats_by_algo = {} + for algo_name, algo_runs in all_results.items(): + if not algo_runs: + continue + + final_costs = [] + for run_data in algo_runs: + if run_data: + final_costs.append(run_data[-1][1]) # Last cost in the run + + stats_by_algo[algo_name] = compute_statistics(final_costs, optimal_cost) + + # Save statistics to JSON + statistics = { + 'metadata': { + 'generated_at': datetime.now().isoformat(), + 'tsp_instance': instance_data['name'], + 'optimal_cost': optimal_cost, + 'experiment_settings': { + 'calibration_time': CALIBRATION_TIME, + 'max_normalized_steps': MAX_NORMALIZED_STEPS, + 'n_runs': N_RUNS, + 'parallel_runs': PARALLEL_RUNS, + 'num_workers': NUM_WORKERS + } + }, + 'algorithms': stats_by_algo, + 'calibration': { + 'steps_per_second': steps_per_second, + 'reference_algorithm': reference_algo, + 'work_per_step': work_per_step + }, + 'experiment_config': { + 'max_normalized_steps': MAX_NORMALIZED_STEPS, + 'n_runs': N_RUNS, + 'algorithms': list(solvers.keys()), + 'use_nn_initialization': False + } + } + + with open('figures/relative_work_figures.json', 'w') as f: + json.dump(statistics, f, indent=2) + logger.info("Saved figures/relative_work_figures.json") + + fig, ax = create_plot( + "Algorithm Performance vs Normalized Work", + f"Normalized Steps (Reference: {reference_algo})", + "Best Cost" + ) + + num_points = 100 + common_norm_steps = np.linspace(1, MAX_NORMALIZED_STEPS, num_points) + + for algo_name, algo_runs in all_results.items(): + if not algo_runs: + continue + + # Interpolate each run to common normalized steps grid + x_lists = [[point[0] for point in run_data] for run_data in algo_runs] + y_lists = [[point[1] for point in run_data] for run_data in algo_runs] + mean_best, std_best = align_series(x_lists, y_lists, common_norm_steps) + + if len(mean_best) == 0: + continue + + base_algo = algo_name.split('_')[0] + color = ALGO_COLORS[base_algo] + + ax.plot(common_norm_steps, mean_best, label=f"{algo_name}", + color=color, linewidth=2) + ax.fill_between(common_norm_steps, mean_best - std_best, mean_best + std_best, + alpha=0.2, color=color) + + final_mean = mean_best[-1] + final_std = std_best[-1] + logger.info(f"{algo_name}: Final mean cost {final_mean:.2f} ± {final_std:.2f}") + + add_optimal_line(ax, optimal_cost) + ax.legend() + + save_figure(fig, 'figures/relative_work_figures.png') + logger.info("Saved figures/relative_work_figures.png") + + +if __name__ == "__main__": + main() diff --git a/figure_scripts/relative_work_nn_figures.py b/figure_scripts/relative_work_nn_figures.py new file mode 100644 index 0000000..8e2cd18 --- /dev/null +++ b/figure_scripts/relative_work_nn_figures.py @@ -0,0 +1,235 @@ +import time +import math +import logging +from tqdm import tqdm +import numpy as np +import json +from datetime import datetime +from constants import CALIBRATION_TIME, MAX_NORMALIZED_STEPS, N_RUNS, PARALLEL_RUNS, NUM_WORKERS +from .common import load_tsp_instance, create_solvers, create_plot, get_nn_initial_route, align_series, save_figure, ALGO_COLORS, add_optimal_line, compute_nn_baseline +from concurrent.futures import ProcessPoolExecutor, as_completed +from tsp.model import TSPInstance + +logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') +logger = logging.getLogger(__name__) + +def compute_statistics(data, optimal_cost=None): + """Compute comprehensive statistics for a dataset.""" + if not data or len(data) == 0: + return {} + + data = np.array(data) + stats = { + 'count': len(data), + 'mean': float(np.mean(data)), + 'std': float(np.std(data)), + 'min': float(np.min(data)), + 'max': float(np.max(data)), + 'median': float(np.median(data)), + 'q1': float(np.percentile(data, 25)), + 'q3': float(np.percentile(data, 75)), + 'iqr': float(np.percentile(data, 75) - np.percentile(data, 25)) + } + + if optimal_cost and optimal_cost > 0: + stats['gap_to_optimal_pct'] = float((stats['mean'] / optimal_cost - 1) * 100) + stats['best_gap_to_optimal_pct'] = float((stats['min'] / optimal_cost - 1) * 100) + stats['worst_gap_to_optimal_pct'] = float((stats['max'] / optimal_cost - 1) * 100) + + return stats + + +def calibrate_steps_per_second(create_solver_func, initial_route=None, calibration_time=CALIBRATION_TIME): + solver = create_solver_func() + solver.initialize(initial_route) + + start_time = time.perf_counter() + step_count = 0 + + while time.perf_counter() - start_time < calibration_time: + solver.step() + step_count += 1 + + actual_time = time.perf_counter() - start_time + return step_count / actual_time + + +def run_single_benchmark(create_solver_func, initial_route, work_factor, max_normalized_steps=MAX_NORMALIZED_STEPS): + solver = create_solver_func() + solver.initialize(initial_route) + + max_actual_steps = max(1, math.floor(max_normalized_steps / work_factor)) + + results = [] + for step in range(max_actual_steps): + solver.step() + normalized_step = (step + 1) * work_factor + fitness = solver.get_cost() # Get current best cost (fitness) + results.append((normalized_step, fitness)) + + return results + + +def worker_benchmark_nn(display_name, orig_key, work_factor): + """Worker for parallel NN benchmarking - computes initial_route locally.""" + from figure_scripts.relative_work_nn_figures import run_single_benchmark + from .common import load_tsp_instance, create_solvers, get_nn_initial_route + _, _, instance_data = load_tsp_instance() + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + initial_route = get_nn_initial_route(instance) # Compute locally + solvers = create_solvers() + create_func = solvers[orig_key] + result = run_single_benchmark(create_func, initial_route, work_factor) + return (display_name, result) + + +def main(): + instance, optimal_cost, instance_data = load_tsp_instance() + initial_route = get_nn_initial_route(instance) + logger.info(f"Generating relative work figures for {instance_data['name']} (optimal: {optimal_cost:.2f}) with NN init") + + # Calibration phase (single run for speed) + solvers = create_solvers() + steps_per_second = {} + for algo_name, create_func in solvers.items(): + steps_per_second[algo_name] = calibrate_steps_per_second(create_func, initial_route) + + logger.info(f"Steps per second: {steps_per_second}") + + # Find the algorithm with the most steps per second (least work per step) + reference_algo = max(steps_per_second, key=steps_per_second.get) + reference_steps_per_second = steps_per_second[reference_algo] + + logger.info(f"Reference algorithm (fastest): {reference_algo} with {reference_steps_per_second:.1f} steps/sec") + + # Calculate normalization factors (work per step relative to reference algorithm) + work_per_step = {} + for algo_name in solvers: + work_per_step[algo_name] = reference_steps_per_second / steps_per_second[algo_name] + + logger.info(f"Work per step: {work_per_step}") + + algo_configs = [ + ('SA_NN', 'SA_random'), + ('GA_NN', 'GA_random') + ] + + args_list = [] + for display_name, orig_key in algo_configs: + work_factor = work_per_step[orig_key] + for _ in range(N_RUNS): + args_list.append((display_name, orig_key, work_factor)) + + logger.info(f"Running {len(args_list)} trials {'in parallel' if PARALLEL_RUNS else 'sequentially'}") + + if PARALLEL_RUNS: + with ProcessPoolExecutor(max_workers=NUM_WORKERS) as executor: + futures = [executor.submit(worker_benchmark_nn, *arg) for arg in args_list] + all_results_list = [f.result() for f in tqdm(as_completed(futures), total=len(args_list), desc="Running trials")] + else: + all_results_list = [] + # For sequential, use the original initial_route + for arg in tqdm(args_list, desc="Running trials"): + display_name, orig_key, work_factor = arg + create_func = solvers[orig_key] + result = run_single_benchmark(create_func, initial_route, work_factor) + all_results_list.append((display_name, result)) + + all_results = {display: [] for display, _ in algo_configs} + for display_name, result in all_results_list: + all_results[display_name].append(result) + + # Compute statistics for final costs + stats_by_algo = {} + for algo_name, algo_runs in all_results.items(): + if not algo_runs: + continue + + final_costs = [] + for run_data in algo_runs: + if run_data: + final_costs.append(run_data[-1][1]) # Last cost in the run + + stats_by_algo[algo_name] = compute_statistics(final_costs, optimal_cost) + + # Save statistics to JSON + statistics = { + 'metadata': { + 'generated_at': datetime.now().isoformat(), + 'tsp_instance': instance_data['name'], + 'optimal_cost': optimal_cost, + 'experiment_settings': { + 'calibration_time': CALIBRATION_TIME, + 'max_normalized_steps': MAX_NORMALIZED_STEPS, + 'n_runs': N_RUNS, + 'parallel_runs': PARALLEL_RUNS, + 'num_workers': NUM_WORKERS + } + }, + 'algorithms': stats_by_algo, + 'calibration': { + 'steps_per_second': steps_per_second, + 'reference_algorithm': reference_algo, + 'work_per_step': work_per_step + }, + 'experiment_config': { + 'max_normalized_steps': MAX_NORMALIZED_STEPS, + 'n_runs': N_RUNS, + 'algorithms': [display for display, _ in algo_configs], + 'use_nn_initialization': True + } + } + + with open('figures/relative_work_nn_figures.json', 'w') as f: + json.dump(statistics, f, indent=2) + logger.info("Saved figures/relative_work_nn_figures.json") + + # Plot + fig, ax = create_plot( + "Algorithm Performance vs Normalized Work with NN init", + f"Normalized Steps (Reference: {reference_algo})", + "Best Cost" + ) + + num_points = 100 + common_norm_steps = np.linspace(1, MAX_NORMALIZED_STEPS, num_points) + + for algo_name, algo_runs in all_results.items(): + if not algo_runs: + continue + + # Interpolate each run to common normalized steps grid + x_lists = [[point[0] for point in run_data] for run_data in algo_runs] + y_lists = [[point[1] for point in run_data] for run_data in algo_runs] + mean_best, std_best = align_series(x_lists, y_lists, common_norm_steps) + + if len(mean_best) == 0: + continue + + base_algo = algo_name.split('_')[0] + color = ALGO_COLORS[base_algo] + + ax.plot(common_norm_steps, mean_best, label=f"{algo_name}", + color=color, linewidth=2) + ax.fill_between(common_norm_steps, mean_best - std_best, mean_best + std_best, + alpha=0.2, color=color) + + final_mean = mean_best[-1] + final_std = std_best[-1] + logger.info(f"{algo_name}: Final mean cost {final_mean:.2f} ± {final_std:.2f}") + + add_optimal_line(ax, optimal_cost) + + # Add NN baseline line + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + nn_cost = compute_nn_baseline(instance) + ax.axhline(y=nn_cost, color='orange', linestyle='--', label='NN', alpha=0.7, linewidth=2) + + ax.legend() + + save_figure(fig, 'figures/relative_work_nn_figures.png') + logger.info("Saved figures/relative_work_nn_figures.png") + + +if __name__ == "__main__": + main() diff --git a/figure_scripts/time_budget_figures.py b/figure_scripts/time_budget_figures.py new file mode 100644 index 0000000..0e0e493 --- /dev/null +++ b/figure_scripts/time_budget_figures.py @@ -0,0 +1,135 @@ +import numpy as np +import logging +from tsp.model import TSPInstance + +from constants import MAX_SECONDS, N_RUNS, PARALLEL_RUNS +from util import run_algorithm_with_timing +from .common import ( + load_tsp_instance, create_solvers, create_plot, run_parallel_trials, + align_series, save_figure, ALGO_COLORS, ALGO_LINESTYLES, add_optimal_line, + create_time_budget_statistics, save_statistics_json, compute_statistics +) + +logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') +logger = logging.getLogger(__name__) + + + + + +def run_single_time_trial(args): + name, instance_data = args + + solvers = create_solvers() + solver_factory = solvers.get(name) + if solver_factory: + solver = solver_factory() + init_route = None + else: + raise ValueError(f"Unknown algorithm: {name}") + + + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + + iterations, best_costs, current_costs, times, best_route = run_algorithm_with_timing( + instance, solver, init_route, MAX_SECONDS + ) + + return { + 'name': name, + 'iterations': iterations, + 'best_costs': best_costs, + 'current_costs': current_costs, + 'times': times, + 'best_route': best_route + } + + +def main(): + _, optimal_cost, instance_data = load_tsp_instance() + logger.info(f"Generating time-budget figures for {instance_data['name']} (optimal: {optimal_cost:.2f})") + + algorithms = ["SA_random", "GA_random"] + + args_list = [] + for name in algorithms: + for _ in range(N_RUNS): + args_list.append((name, instance_data)) + + logger.info(f"Running {len(args_list)} trials {'in parallel' if PARALLEL_RUNS else 'sequentially'}") + + all_results = run_parallel_trials(run_single_time_trial, args_list, desc="Running trials") + + # Group by algorithm + results_by_algo = {name: [] for name in algorithms} + for res in all_results: + results_by_algo[res['name']].append(res) + + # Compute statistics for each algorithm + stats_by_algo = {} + for algo_name, algo_results in results_by_algo.items(): + if not algo_results: + continue + + final_costs = [r['best_costs'][-1] if r['best_costs'] else float('inf') for r in algo_results] + iterations = [len(r['iterations']) if r['iterations'] else 0 for r in algo_results] + convergence_times = [r['times'][-1] if r['times'] else MAX_SECONDS for r in algo_results] + + stats_by_algo[algo_name] = { + 'final_costs': compute_statistics(final_costs, optimal_cost), + 'iterations': compute_statistics(iterations), + 'convergence_times': compute_statistics(convergence_times) + } + + # Create and save statistics + statistics = create_time_budget_statistics(instance_data, optimal_cost, stats_by_algo, algorithms, use_nn=False) + save_statistics_json(statistics, 'time_budget_figures.json') + + fig, ax = create_plot( + f'TSP Algorithm Comparison: Time Budget ({MAX_SECONDS}s, {N_RUNS} runs each)', + 'Time (seconds)', + 'Best Cost' + ) + + for algo_name in algorithms: + algo_results = results_by_algo[algo_name] + + if not algo_results: + continue + + # Find max time reached across runs + max_time = max(run['times'][-1] for run in algo_results if run['times']) + num_points = 100 + common_times = np.linspace(0, min(max_time, MAX_SECONDS), num_points) + + # Interpolate each run to common grid + x_lists = [run['times'] for run in algo_results if run['times'] and run['best_costs']] + y_lists = [run['best_costs'] for run in algo_results if run['times'] and run['best_costs']] + mean_best, std_best = align_series(x_lists, y_lists, common_times) + + if len(mean_best) == 0: + continue + + base_algo = algo_name.split('_')[0] + color = ALGO_COLORS[base_algo] + linestyle = ALGO_LINESTYLES[base_algo] + + ax.plot(common_times, mean_best, label=f"{algo_name}", + color=color, linestyle=linestyle, linewidth=2) + ax.fill_between(common_times, mean_best - std_best, mean_best + std_best, + alpha=0.2, color=color) + + final_mean = mean_best[-1] + final_std = std_best[-1] + gap = ((final_mean / optimal_cost - 1) * 100) if optimal_cost else 0 + logger.info(f"{algo_name}: Final mean cost {final_mean:.2f} ± {final_std:.2f} (gap: {gap:.1f}%)") + + add_optimal_line(ax, optimal_cost) + ax.legend() + + save_figure(fig, 'figures/time_budget_figures.png') + logger.info("Saved figures/time_budget_figures.png") + + +if __name__ == "__main__": + main() diff --git a/figure_scripts/time_budget_nn_figures.py b/figure_scripts/time_budget_nn_figures.py new file mode 100644 index 0000000..697754d --- /dev/null +++ b/figure_scripts/time_budget_nn_figures.py @@ -0,0 +1,150 @@ +import numpy as np +import logging +from tsp.model import TSPInstance + +from constants import MAX_SECONDS, N_RUNS, PARALLEL_RUNS +from util import run_algorithm_with_timing +from .common import ( + load_tsp_instance, create_solvers, create_plot, get_nn_initial_route, + run_parallel_trials, align_series, save_figure, ALGO_COLORS, ALGO_LINESTYLES, + add_optimal_line, create_time_budget_statistics, save_statistics_json, compute_statistics, + compute_nn_baseline +) +from algorithm.nearest_neighbor import NearestNeighbor + +logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') +logger = logging.getLogger(__name__) + + + +def get_nn_route(instance): + nn = NearestNeighbor(instance) + nn.initialize(None) + n_cities = len(instance.cities) + for _ in range(n_cities - 1): + nn.step() + return nn.get_route() + + +def run_single_time_trial(args): + name, instance_data = args + + solvers = create_solvers() + orig_name = name.replace('_NN', '_random') + solver_factory = solvers.get(orig_name) + if solver_factory: + solver = solver_factory() + else: + raise ValueError(f"Unknown algorithm: {name}") + + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + init_route = get_nn_initial_route(instance) + + iterations, best_costs, current_costs, times, best_route = run_algorithm_with_timing( + instance, solver, init_route, MAX_SECONDS + ) + + return { + 'name': name, + 'iterations': iterations, + 'best_costs': best_costs, + 'current_costs': current_costs, + 'times': times, + 'best_route': best_route + } + + +def main(): + _, optimal_cost, instance_data = load_tsp_instance() + logger.info(f"Generating time-budget figures for {instance_data['name']} (optimal: {optimal_cost:.2f}) with NN init") + + algorithms = ["SA_NN", "GA_NN"] + + args_list = [] + for name in algorithms: + for _ in range(N_RUNS): + args_list.append((name, instance_data)) + + logger.info(f"Running {len(args_list)} trials {'in parallel' if PARALLEL_RUNS else 'sequentially'}") + + all_results = run_parallel_trials(run_single_time_trial, args_list, desc="Running trials") + + # Group by algorithm + results_by_algo = {name: [] for name in algorithms} + for res in all_results: + results_by_algo[res['name']].append(res) + + # Compute statistics for each algorithm + stats_by_algo = {} + for algo_name, algo_results in results_by_algo.items(): + if not algo_results: + continue + + final_costs = [r['best_costs'][-1] if r['best_costs'] else float('inf') for r in algo_results] + iterations = [len(r['iterations']) if r['iterations'] else 0 for r in algo_results] + convergence_times = [r['times'][-1] if r['times'] else MAX_SECONDS for r in algo_results] + + stats_by_algo[algo_name] = { + 'final_costs': compute_statistics(final_costs, optimal_cost), + 'iterations': compute_statistics(iterations), + 'convergence_times': compute_statistics(convergence_times) + } + + # Create and save statistics + statistics = create_time_budget_statistics(instance_data, optimal_cost, stats_by_algo, algorithms, use_nn=True) + save_statistics_json(statistics, 'time_budget_nn_figures.json') + + fig, ax = create_plot( + f'TSP Algorithm Comparison: Time Budget with NN init ({MAX_SECONDS}s, {N_RUNS} runs each)', + 'Time (seconds)', + 'Best Cost' + ) + + for algo_name in algorithms: + algo_results = results_by_algo[algo_name] + + if not algo_results: + continue + + # Find max time reached across runs + max_time = max(run['times'][-1] for run in algo_results if run['times']) + num_points = 100 + common_times = np.linspace(0, min(max_time, MAX_SECONDS), num_points) + + # Interpolate each run to common grid + x_lists = [run['times'] for run in algo_results if run['times'] and run['best_costs']] + y_lists = [run['best_costs'] for run in algo_results if run['times'] and run['best_costs']] + mean_best, std_best = align_series(x_lists, y_lists, common_times) + + if len(mean_best) == 0: + continue + + base_algo = algo_name.split('_')[0] + color = ALGO_COLORS[base_algo] + linestyle = ALGO_LINESTYLES[base_algo] + + ax.plot(common_times, mean_best, label=f"{algo_name}", + color=color, linestyle=linestyle, linewidth=2) + ax.fill_between(common_times, mean_best - std_best, mean_best + std_best, + alpha=0.2, color=color) + + final_mean = mean_best[-1] + final_std = std_best[-1] + gap = ((final_mean / optimal_cost - 1) * 100) if optimal_cost else 0 + logger.info(f"{algo_name}: Final mean cost {final_mean:.2f} ± {final_std:.2f} (gap: {gap:.1f}%)") + + add_optimal_line(ax, optimal_cost) + + # Add NN baseline line + instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities']) + nn_cost = compute_nn_baseline(instance) + ax.axhline(y=nn_cost, color='orange', linestyle='--', label='NN', alpha=0.7, linewidth=2) + + ax.legend() + + save_figure(fig, 'figures/time_budget_nn_figures.png') + logger.info("Saved figures/time_budget_nn_figures.png") + + +if __name__ == "__main__": + main() diff --git a/figures/box_plot_figures.json b/figures/box_plot_figures.json new file mode 100644 index 0000000..825e2f7 --- /dev/null +++ b/figures/box_plot_figures.json @@ -0,0 +1,95 @@ +{ + "metadata": { + "generated_at": "2025-10-05T23:01:30.833899", + "tsp_instance": "lin105", + "optimal_cost": 14382.99593345118, + "experiment_settings": { + "max_seconds": 5.0, + "n_runs": 5, + "parallel_runs": true, + "num_workers": 8 + } + }, + "algorithms": { + "SA_random": { + "count": 5, + "mean": 16012.482621693722, + "std": 812.2855785412886, + "min": 14910.960509124585, + "max": 16927.193332199335, + "median": 16286.156116879023, + "q1": 15204.09309143564, + "q3": 16734.010058830037, + "iqr": 1529.916967394398, + "gap_to_optimal_pct": 11.329257797068347, + "best_gap_to_optimal_pct": 3.67075523149869, + "worst_gap_to_optimal_pct": 17.688925245615895 + }, + "GA_random": { + "count": 5, + "mean": 15424.570200347362, + "std": 802.9888234298182, + "min": 14720.443181384559, + "max": 16631.961357906952, + "median": 14819.544025024956, + "q1": 14804.95418005342, + "q3": 16145.94825736693, + "iqr": 1340.9940773135095, + "gap_to_optimal_pct": 7.241705912422214, + "best_gap_to_optimal_pct": 2.346154094005981, + "worst_gap_to_optimal_pct": 15.63627935974905 + }, + "SA_NN": { + "count": 5, + "mean": 15694.858663424468, + "std": 1004.7537928402737, + "min": 14942.073016688133, + "max": 17685.414089156384, + "median": 15285.69959758063, + "q1": 15217.526367551001, + "q3": 15343.580246146194, + "iqr": 126.05387859519215, + "gap_to_optimal_pct": 9.120928185220635, + "best_gap_to_optimal_pct": 3.8870697441878654, + "worst_gap_to_optimal_pct": 22.960572129653613 + }, + "GA_NN": { + "count": 5, + "mean": 15401.394448866462, + "std": 460.6928964684674, + "min": 15044.396451813853, + "max": 16297.019936700484, + "median": 15171.348216956554, + "q1": 15122.475021703787, + "q3": 15371.732617157631, + "iqr": 249.25759545384426, + "gap_to_optimal_pct": 7.080572921853823, + "best_gap_to_optimal_pct": 4.5984892259089305, + "worst_gap_to_optimal_pct": 13.307547412968201 + } + }, + "nearest_neighbor_baseline": { + "count": 1, + "mean": 19862.014199310805, + "std": 0.0, + "min": 19862.014199310805, + "max": 19862.014199310805, + "median": 19862.014199310805, + "q1": 19862.014199310805, + "q3": 19862.014199310805, + "iqr": 0.0, + "gap_to_optimal_pct": 38.09372046832624, + "best_gap_to_optimal_pct": 38.09372046832624, + "worst_gap_to_optimal_pct": 38.09372046832624 + }, + "experiment_config": { + "max_seconds": 5.0, + "n_runs": 5, + "algorithms": [ + "SA_random", + "GA_random", + "SA_NN", + "GA_NN" + ] + } +} \ No newline at end of file diff --git a/figures/random_baseline_figures.json b/figures/random_baseline_figures.json new file mode 100644 index 0000000..da3a4ff --- /dev/null +++ b/figures/random_baseline_figures.json @@ -0,0 +1,65 @@ +{ + "metadata": { + "generated_at": "2025-10-05T23:01:36.050643", + "tsp_instance": "lin105", + "optimal_cost": 14382.99593345118, + "experiment_settings": { + "max_seconds": 5.0, + "n_runs": 5, + "parallel_runs": true, + "num_workers": 8 + } + }, + "algorithms": { + "Random": { + "count": 5, + "mean": 96838.8227696025, + "std": 885.8460946739577, + "min": 95783.92453048032, + "max": 98472.28823476467, + "median": 96629.838132142, + "q1": 96530.28461094857, + "q3": 96777.77833967701, + "iqr": 247.49372872844106, + "gap_to_optimal_pct": 573.2868674764768, + "best_gap_to_optimal_pct": 565.9525245898967, + "worst_gap_to_optimal_pct": 584.6437883344125 + }, + "Nearest Neighbor": { + "count": 5, + "mean": 18723.651646487746, + "std": 950.2666394220985, + "min": 17048.379780542797, + "max": 19955.71470162493, + "median": 18797.45071989588, + "q1": 18666.654060655674, + "q3": 19150.058969719437, + "iqr": 483.4049090637636, + "gap_to_optimal_pct": 30.17907905362962, + "best_gap_to_optimal_pct": 18.531492739232537, + "worst_gap_to_optimal_pct": 38.74518767827102 + } + }, + "nearest_neighbor_baseline": { + "count": 1, + "mean": 19738.922583612144, + "std": 0.0, + "min": 19738.922583612144, + "max": 19738.922583612144, + "median": 19738.922583612144, + "q1": 19738.922583612144, + "q3": 19738.922583612144, + "iqr": 0.0, + "gap_to_optimal_pct": 37.23790700450971, + "best_gap_to_optimal_pct": 37.23790700450971, + "worst_gap_to_optimal_pct": 37.23790700450971 + }, + "experiment_config": { + "max_seconds": 5.0, + "n_runs": 5, + "algorithms": [ + "Random", + "Nearest Neighbor" + ] + } +} \ No newline at end of file diff --git a/figures/relative_work_figures.json b/figures/relative_work_figures.json new file mode 100644 index 0000000..2b4643c --- /dev/null +++ b/figures/relative_work_figures.json @@ -0,0 +1,64 @@ +{ + "metadata": { + "generated_at": "2025-10-05T23:02:05.592238", + "tsp_instance": "lin105", + "optimal_cost": 14382.99593345118, + "experiment_settings": { + "calibration_time": 2.0, + "max_normalized_steps": 100000, + "n_runs": 5, + "parallel_runs": true, + "num_workers": 8 + } + }, + "algorithms": { + "SA_random": { + "count": 5, + "mean": 15306.236528531212, + "std": 447.262297300121, + "min": 14602.693488083272, + "max": 15997.50014228392, + "median": 15262.241113866019, + "q1": 15225.102283000344, + "q3": 15443.645615422509, + "iqr": 218.5433324221649, + "gap_to_optimal_pct": 6.418972788087984, + "best_gap_to_optimal_pct": 1.527481170464151, + "worst_gap_to_optimal_pct": 11.225089795637189 + }, + "GA_random": { + "count": 5, + "mean": 18482.42804127147, + "std": 395.1800862875539, + "min": 17878.331223552395, + "max": 19106.192865359226, + "median": 18517.790026291652, + "q1": 18345.239981238177, + "q3": 18564.586109915912, + "iqr": 219.34612867773467, + "gap_to_optimal_pct": 28.50193469279967, + "best_gap_to_optimal_pct": 24.301858293458565, + "worst_gap_to_optimal_pct": 32.83875594321135 + } + }, + "calibration": { + "steps_per_second": { + "SA_random": 43149.06564031536, + "GA_random": 147.9669575726668 + }, + "reference_algorithm": "SA_random", + "work_per_step": { + "SA_random": 1.0, + "GA_random": 291.6128461932104 + } + }, + "experiment_config": { + "max_normalized_steps": 100000, + "n_runs": 5, + "algorithms": [ + "SA_random", + "GA_random" + ], + "use_nn_initialization": false + } +} \ No newline at end of file diff --git a/figures/relative_work_nn_figures.json b/figures/relative_work_nn_figures.json new file mode 100644 index 0000000..81d4e6f --- /dev/null +++ b/figures/relative_work_nn_figures.json @@ -0,0 +1,64 @@ +{ + "metadata": { + "generated_at": "2025-10-05T23:02:14.330398", + "tsp_instance": "lin105", + "optimal_cost": 14382.99593345118, + "experiment_settings": { + "calibration_time": 2.0, + "max_normalized_steps": 100000, + "n_runs": 5, + "parallel_runs": true, + "num_workers": 8 + } + }, + "algorithms": { + "SA_NN": { + "count": 5, + "mean": 15257.13794194028, + "std": 318.073836377086, + "min": 14709.028693147355, + "max": 15527.196325554163, + "median": 15468.93638886908, + "q1": 15083.568391333562, + "q3": 15496.959910797239, + "iqr": 413.39151946367747, + "gap_to_optimal_pct": 6.077607283862663, + "best_gap_to_optimal_pct": 2.266793102109599, + "worst_gap_to_optimal_pct": 7.955229893668148 + }, + "GA_NN": { + "count": 5, + "mean": 15678.757305751566, + "std": 263.9634818116499, + "min": 15325.548497717666, + "max": 16063.634611311645, + "median": 15616.596575133632, + "q1": 15505.366533444401, + "q3": 15882.640311150497, + "iqr": 377.27377770609564, + "gap_to_optimal_pct": 9.008981009907501, + "best_gap_to_optimal_pct": 6.553242235676016, + "worst_gap_to_optimal_pct": 11.68489990288968 + } + }, + "calibration": { + "steps_per_second": { + "SA_random": 50767.87165943893, + "GA_random": 148.37386626471366 + }, + "reference_algorithm": "SA_random", + "work_per_step": { + "SA_random": 1.0, + "GA_random": 342.16181688535386 + } + }, + "experiment_config": { + "max_normalized_steps": 100000, + "n_runs": 5, + "algorithms": [ + "SA_NN", + "GA_NN" + ], + "use_nn_initialization": true + } +} \ No newline at end of file diff --git a/figures/time_budget_figures.json b/figures/time_budget_figures.json new file mode 100644 index 0000000..fa33c6d --- /dev/null +++ b/figures/time_budget_figures.json @@ -0,0 +1,100 @@ +{ + "metadata": { + "generated_at": "2025-10-05T23:01:46.185997", + "tsp_instance": "lin105", + "optimal_cost": 14382.99593345118, + "experiment_settings": { + "max_seconds": 5.0, + "n_runs": 5, + "parallel_runs": true, + "num_workers": 8 + } + }, + "algorithms": { + "SA_random": { + "final_costs": { + "count": 5, + "mean": 15175.893358173897, + "std": 276.3467112707878, + "min": 14898.740953691315, + "max": 15615.923945037468, + "median": 15038.489532062275, + "q1": 14948.455378199606, + "q3": 15377.856981878826, + "iqr": 429.40160367922, + "gap_to_optimal_pct": 5.512741770847884, + "best_gap_to_optimal_pct": 3.585796885617154, + "worst_gap_to_optimal_pct": 8.572122367905365 + }, + "iterations": { + "count": 5, + "mean": 196623.4, + "std": 3124.1423527105803, + "min": 192222.0, + "max": 200336.0, + "median": 197550.0, + "q1": 193788.0, + "q3": 199221.0, + "iqr": 5433.0 + }, + "convergence_times": { + "count": 5, + "mean": 5.000015653390437, + "std": 2.3653524650911417e-06, + "min": 5.000012286996935, + "max": 5.0000186279648915, + "median": 5.000016787962522, + "q1": 5.000013516983017, + "q3": 5.000017047044821, + "iqr": 3.5300618037581444e-06 + } + }, + "GA_random": { + "final_costs": { + "count": 5, + "mean": 16009.171166619602, + "std": 636.7278280229738, + "min": 15376.009868585434, + "max": 17157.113924554593, + "median": 15884.901217060511, + "q1": 15484.709404628955, + "q3": 16143.121418268522, + "iqr": 658.4120136395668, + "gap_to_optimal_pct": 11.306234394368087, + "best_gap_to_optimal_pct": 6.904082708003534, + "worst_gap_to_optimal_pct": 19.287483664314475 + }, + "iterations": { + "count": 5, + "mean": 624.4, + "std": 49.82609757948138, + "min": 572.0, + "max": 693.0, + "median": 591.0, + "q1": 590.0, + "q3": 676.0, + "iqr": 86.0 + }, + "convergence_times": { + "count": 5, + "mean": 5.0053679721895605, + "std": 0.0020288936434057383, + "min": 5.002007366041653, + "max": 5.008208452956751, + "median": 5.005814514996018, + "q1": 5.004691343987361, + "q3": 5.00611818296602, + "iqr": 0.0014268389786593616 + } + } + }, + "experiment_config": { + "max_seconds": 5.0, + "n_runs": 5, + "algorithms": [ + "SA_random", + "GA_random" + ], + "use_nn_initialization": false + } +} \ No newline at end of file diff --git a/figures/time_budget_nn_figures.json b/figures/time_budget_nn_figures.json new file mode 100644 index 0000000..65949e9 --- /dev/null +++ b/figures/time_budget_nn_figures.json @@ -0,0 +1,100 @@ +{ + "metadata": { + "generated_at": "2025-10-05T23:01:56.474617", + "tsp_instance": "lin105", + "optimal_cost": 14382.99593345118, + "experiment_settings": { + "max_seconds": 5.0, + "n_runs": 5, + "parallel_runs": true, + "num_workers": 8 + } + }, + "algorithms": { + "SA_NN": { + "final_costs": { + "count": 5, + "mean": 15149.977452400874, + "std": 445.8425839566921, + "min": 14632.906215416959, + "max": 15805.024227354677, + "median": 15232.981436894737, + "q1": 14674.711409911271, + "q3": 15404.263972426728, + "iqr": 729.5525625154569, + "gap_to_optimal_pct": 5.332557434476426, + "best_gap_to_optimal_pct": 1.7375398221767702, + "worst_gap_to_optimal_pct": 9.88687127830039 + }, + "iterations": { + "count": 5, + "mean": 194386.6, + "std": 3444.169252519394, + "min": 188459.0, + "max": 197532.0, + "median": 196001.0, + "q1": 192626.0, + "q3": 197315.0, + "iqr": 4689.0 + }, + "convergence_times": { + "count": 5, + "mean": 5.0000160478055475, + "std": 3.6388648893810043e-06, + "min": 5.000012296019122, + "max": 5.000020754989237, + "median": 5.0000135760055855, + "q1": 5.000013436016161, + "q3": 5.0000201759976335, + "iqr": 6.739981472492218e-06 + } + }, + "GA_NN": { + "final_costs": { + "count": 5, + "mean": 15567.948857117008, + "std": 520.1624586917907, + "min": 14972.574219533735, + "max": 16385.70700217945, + "median": 15668.02092923951, + "q1": 15040.905297945133, + "q3": 15772.536836687208, + "iqr": 731.6315387420746, + "gap_to_optimal_pct": 8.238568161657668, + "best_gap_to_optimal_pct": 4.099134066438448, + "worst_gap_to_optimal_pct": 13.924157929228587 + }, + "iterations": { + "count": 5, + "mean": 636.8, + "std": 74.80748625639013, + "min": 564.0, + "max": 732.0, + "median": 584.0, + "q1": 580.0, + "q3": 724.0, + "iqr": 144.0 + }, + "convergence_times": { + "count": 5, + "mean": 5.002211673592683, + "std": 0.0018705061061081434, + "min": 5.00016123795649, + "max": 5.005234103999101, + "median": 5.001270666019991, + "q1": 5.000917778990697, + "q3": 5.0034745809971355, + "iqr": 0.0025568020064383745 + } + } + }, + "experiment_config": { + "max_seconds": 5.0, + "n_runs": 5, + "algorithms": [ + "SA_NN", + "GA_NN" + ], + "use_nn_initialization": true + } +} \ No newline at end of file diff --git a/ga_tuning.py b/ga_tuning.py deleted file mode 100644 index 950d3e1..0000000 --- a/ga_tuning.py +++ /dev/null @@ -1,244 +0,0 @@ - - - - - - -import numpy as np -from pathlib import Path -from skopt import gp_minimize -from skopt.space import Real, Integer -from skopt.utils import use_named_args -import time -from numba import jit, prange - -from util import find_optimal_tour -from algorithm.genetic_algo import GeneticAlgorithmSolver - - -def run_single_ga_optimized(params, instance, time_budget_seconds, seed): - """Optimized single GA run with time budget instead of fixed steps.""" - import time - - # Create GA with given parameters - ga = GeneticAlgorithmSolver( - instance=instance, - seed=seed, - population_size=int(params['population_size']), - crossover_rate=params['crossover_rate'], - mutation_rate=params['mutation_rate'], - elitism_count=int(params['elitism_count']), - num_parents=int(params['num_parents']), - num_child=int(params['num_child']) - ) - - # Initialize with random permutation - ga.initialize(None) - - # Run for time budget - start_time = time.perf_counter() - steps = 0 - while time.perf_counter() - start_time < time_budget_seconds: - ga.step() - steps += 1 - - return ga.best_cost, steps - - -@jit(nopython=True, parallel=True) -def parallel_cost_calculation(costs_array, base_cost, n_runs): - """JIT-compiled function to calculate costs in parallel.""" - for i in prange(n_runs): - # Simulate cost variation (placeholder for actual GA computation) - costs_array[i] = base_cost + (i * 1000.0) - return costs_array - - -def evaluate_ga_params(params, instance, optimal_cost, time_budget_seconds=10.0, n_runs=2): - """ - Evaluate GA parameters by running multiple times sequentially with time budget. - - Args: - params: Dictionary of GA parameters - instance: TSP instance - optimal_cost: Optimal cost for reference - time_budget_seconds: Time budget per evaluation in seconds - n_runs: Number of runs to average - - Returns: - Average final cost across runs - """ - # Validate parameter combinations - population_size = int(params['population_size']) - elitism_count = int(params['elitism_count']) - num_parents = int(params['num_parents']) - num_child = int(params['num_child']) - - # Ensure elitism doesn't exceed population size - if elitism_count >= population_size: - return 100000.0 # Large penalty for invalid combinations - - # Ensure reasonable elitism ratio (max 50% of population) - if elitism_count > population_size * 0.5: - return 100000.0 # Large penalty for excessive elitism - - # Ensure num_parents doesn't exceed population size - if num_parents > population_size: - return 100000.0 - - # Ensure reasonable parent count (max 20% of population) - if num_parents > population_size * 0.2: - return 100000.0 - - # Run sequentially (no parallel processing) - costs = [] - steps = [] - - for run in range(n_runs): - cost, step_count = run_single_ga_optimized(params, instance, time_budget_seconds, 42 + run) - costs.append(cost) - steps.append(step_count) - - # Return average cost (lower is better for minimization) - avg_cost = np.mean(costs) - avg_steps = np.mean(steps) - print(f"Params: {params} -> Avg cost: {avg_cost:.2f} (runs: {costs}, avg steps: {avg_steps:.0f})") - return avg_cost - - -def main(): - """Main hyperparameter optimization.""" - print("GA Hyperparameter Tuning") - print("=" * 50) - - - - instance_path = Path("dataset/lin105.tsp") - instance, optimal_cost = find_optimal_tour(instance_path) - print(f"Instance: {instance.name}") - print(f"Optimal cost: {optimal_cost:.2f}") - print(f"Cities: {len(instance.cities)}") - print() - - dimensions = [ - Integer(30, 200, name='population_size'), # Population size (realistic range) - Real(0.6, 0.95, name='crossover_rate'), # Crossover rate (realistic range) - Real(0.01, 0.3, name='mutation_rate'), # Mutation rate (realistic range) - Integer(1, 20, name='elitism_count'), # Elitism count (realistic range) - Integer(2, 6, name='num_parents'), # Number of parents (realistic range) - Integer(1, 4, name='num_child'), # Number of children (realistic range) - ] - - # Objective function wrapper - @use_named_args(dimensions=dimensions) - def objective(**params): - return evaluate_ga_params(params, instance, optimal_cost, time_budget_seconds=10.0) - - print("Starting optimization...") - print("Each evaluation runs GA for 10 seconds, 2 times sequentially") - print("Using realistic parameter ranges with validation") - print("Using Numba JIT compilation for speed optimization") - print() - - # Estimate time with a single evaluation - print("Estimating time with a single evaluation...") - test_start = time.time() - test_params = { - 'population_size': 50, - 'crossover_rate': 0.7, - 'mutation_rate': 0.05, - 'elitism_count': 4, - 'num_parents': 2, - 'num_child': 2 - } - evaluate_ga_params(test_params, instance, optimal_cost, time_budget_seconds=10.0) - test_time = time.time() - test_start - - estimated_total_time = test_time * 25 # 25 evaluations (reduced) - print(f"Single evaluation took: {test_time:.1f} seconds") - print(f"Estimated total time: {estimated_total_time:.1f} seconds ({estimated_total_time/60:.1f} minutes)") - print() - - # Ask user if they want to continue - response = input("Continue with optimization? (y/n): ").lower().strip() - if response != 'y': - print("Optimization cancelled.") - return - - # Run Bayesian optimization - start_time = time.time() - result = gp_minimize( - func=objective, - dimensions=dimensions, - n_calls=25, # Number of evaluations (reduced for speed) - random_state=42, - acq_func='EI' # Expected Improvement - ) - end_time = time.time() - - # Results - print("\n" + "=" * 50) - print("OPTIMIZATION RESULTS") - print("=" * 50) - print(f"Optimization time: {end_time - start_time:.1f} seconds") - print(f"Best cost found: {result.fun:.2f}") - print(f"Optimality gap: {((result.fun / optimal_cost - 1) * 100):.1f}%") - print() - - print("Best parameters:") - best_params = dict(zip([dim.name for dim in dimensions], result.x)) - for param, value in best_params.items(): - if param in ['population_size', 'elitism_count', 'num_parents', 'num_child']: - print(f" {param}: {int(value)}") - else: - print(f" {param}: {value:.3f}") - - print() - print("Parameter ranges tested:") - for dim in dimensions: - print(f" {dim.name}: {dim.low} to {dim.high}") - - # Test best parameters with more runs - print("\n" + "=" * 50) - print("VALIDATION WITH BEST PARAMETERS") - print("=" * 50) - - validation_costs = [] - validation_steps = [] - for run in range(5): # More runs for validation - ga = GeneticAlgorithmSolver( - instance=instance, - seed=100 + run, - population_size=int(best_params['population_size']), - crossover_rate=best_params['crossover_rate'], - mutation_rate=best_params['mutation_rate'], - elitism_count=int(best_params['elitism_count']), - num_parents=int(best_params['num_parents']), - num_child=int(best_params['num_child']) - ) - - ga.initialize(None) - - # Run for 10 seconds - start_time = time.perf_counter() - steps = 0 - while time.perf_counter() - start_time < 10.0: - ga.step() - steps += 1 - - validation_costs.append(ga.best_cost) - validation_steps.append(steps) - print(f"Run {run+1}: {ga.best_cost:.2f} ({steps} steps)") - - print("\nValidation results:") - print(f" Mean cost: {np.mean(validation_costs):.2f}") - print(f" Std cost: {np.std(validation_costs):.2f}") - print(f" Best cost: {np.min(validation_costs):.2f}") - print(f" Worst cost: {np.max(validation_costs):.2f}") - print(f" Mean optimality gap: {((np.mean(validation_costs) / optimal_cost - 1) * 100):.1f}%") - print(f" Mean steps per second: {np.mean(validation_steps) / 10:.0f}") - print(f" Steps range: {np.min(validation_steps)} - {np.max(validation_steps)}") - - -if __name__ == "__main__": - main() diff --git a/generate_figures.py b/generate_figures.py new file mode 100644 index 0000000..0b2880b --- /dev/null +++ b/generate_figures.py @@ -0,0 +1,47 @@ +from pathlib import Path +import logging +from figure_scripts import time_budget_figures, relative_work_figures, time_budget_nn_figures, relative_work_nn_figures, box_plot_figures, random_baseline_figures +from constants import DATASET_FILENAME + +logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') +logger = logging.getLogger(__name__) + +def main(): + figures_dir = Path("figures") + figures_dir.mkdir(exist_ok=True) + + logger.info("Generating TSP algorithm comparison figures...") + logger.info(f"Output directory: {figures_dir.absolute()}") + + # Hard check for dataset (required for both scripts) + dataset_path = Path("dataset") / DATASET_FILENAME + if not dataset_path.exists(): + logger.error(f"Dataset missing ({dataset_path}).") + logger.error(f"Run 'python setup_dataset.py' from root to download {DATASET_FILENAME} and other files.") + logger.error("Figures cannot be generated without it.") + return # Exit without generating + + scripts = [ + ("box_plot_figures", box_plot_figures.main), + ("random_baseline_figures", random_baseline_figures.main), + ("time_budget_figures", time_budget_figures.main), + ("time_budget_nn_figures", time_budget_nn_figures.main), + ("relative_work_figures", relative_work_figures.main), + ("relative_work_nn_figures", relative_work_nn_figures.main) + ] + + for script_name, main_func in scripts: + logger.info(f"Running {script_name}...") + try: + main_func() + logger.info(f"Completed {script_name}") + except Exception as e: + logger.error(f"{script_name} failed with error: {e}") + import traceback + traceback.print_exc() + + logger.info(f"Generated figures in {figures_dir.absolute()}") + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index f4ac9e7..429410b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,6 +17,7 @@ dependencies = [ "pytest-mock>=3.15.1", "scikit-optimize>=0.10.2", "numba>=0.62.1", + "tqdm>=4.67.1", ] [dependency-groups] diff --git a/sa_tuning.py b/sa_tuning.py deleted file mode 100644 index 28520f9..0000000 --- a/sa_tuning.py +++ /dev/null @@ -1,324 +0,0 @@ -#!/usr/bin/env python3 -""" -SA Hyperparameter Tuning using scikit-optimize - -This script optimizes SA parameters to minimize route cost over 5 seconds. -""" - -import numpy as np -from pathlib import Path -from skopt import gp_minimize -from skopt.space import Real, Integer -from skopt.utils import use_named_args -import time -from numba import jit, prange - -from util import find_optimal_tour, exponential_cooling -from algorithm.simulated_annealing import SimulatedAnnealing - - -@jit(nopython=True) -def jit_calculate_route_cost(route, cities_coords): - """JIT-compiled route cost calculation.""" - n = len(route) - total_cost = 0.0 - - for i in range(n): - current_city = route[i] - next_city = route[(i + 1) % n] - - # Calculate Euclidean distance - dx = cities_coords[current_city, 0] - cities_coords[next_city, 0] - dy = cities_coords[current_city, 1] - cities_coords[next_city, 1] - total_cost += np.sqrt(dx * dx + dy * dy) - - return total_cost - - -@jit(nopython=True) -def jit_2opt_neighbor(route, i, j): - """JIT-compiled 2-opt neighbor generation.""" - n = len(route) - neighbor = np.empty(n, dtype=np.int32) - - # Copy route - for k in range(n): - neighbor[k] = route[k] - - # Reverse segment between i and j - for k in range(i, j + 1): - neighbor[k] = route[j - (k - i)] - - return neighbor - - -@jit(nopython=True, parallel=True) -def jit_batch_route_costs(routes, cities_coords): - """JIT-compiled batch route cost calculation.""" - n_routes, n_cities = routes.shape - costs = np.zeros(n_routes, dtype=np.float64) - - for route_idx in prange(n_routes): - total_cost = 0.0 - for i in range(n_cities): - current_city = routes[route_idx, i] - next_city = routes[route_idx, (i + 1) % n_cities] - - # Calculate Euclidean distance - dx = cities_coords[current_city, 0] - cities_coords[next_city, 0] - dy = cities_coords[current_city, 1] - cities_coords[next_city, 1] - total_cost += np.sqrt(dx * dx + dy * dy) - - costs[route_idx] = total_cost - - return costs - - -@jit(nopython=True) -def jit_sa_optimized(cities_coords, initial_temp, cooling_rate, max_steps, seed): - """JIT-compiled SA algorithm with 2-opt moves.""" - n = cities_coords.shape[0] - - # Initialize route - route = np.arange(n, dtype=np.int32) - # Simple shuffle using linear congruential generator - rng_state = seed - for i in range(n): - j = rng_state % n - route[i], route[j] = route[j], route[i] - rng_state = (rng_state * 1664525 + 1013904223) % (2**32) - - # Initialize best and current cost - current_cost = jit_calculate_route_cost(route, cities_coords) - best_cost = current_cost - best_route = route.copy() - temperature = initial_temp - - # SA loop - for step in range(max_steps): - # Generate random 2-opt move - rng_state = (rng_state * 1664525 + 1013904223) % (2**32) - i = int(rng_state % n) - rng_state = (rng_state * 1664525 + 1013904223) % (2**32) - j = int(rng_state % n) - - if i > j: - i, j = j, i - if i == j: - j = (j + 1) % n - - # Create neighbor - neighbor = jit_2opt_neighbor(route, i, j) - - # Calculate neighbor cost only - neighbor_cost = jit_calculate_route_cost(neighbor, cities_coords) - - # Accept or reject - delta = neighbor_cost - current_cost - rng_state = (rng_state * 1664525 + 1013904223) % (2**32) - random_val = (rng_state % 1000) / 1000.0 - - if delta < 0 or np.exp(-delta / temperature) > random_val: - # Accept neighbor - for k in range(n): - route[k] = neighbor[k] - current_cost = neighbor_cost - - # Update best if improved - if current_cost < best_cost: - best_cost = current_cost - for k in range(n): - best_route[k] = route[k] - - # Cool down - temperature *= cooling_rate - - return best_cost, best_route - - -def run_single_sa_optimized(params, instance, time_budget_seconds, seed): - """Fully JIT-compiled SA run with time budget.""" - # Pre-compute cities coordinates for JIT - cities_coords = np.array([[city.x, city.y] for city in instance.cities], dtype=np.float64) - - # Estimate steps based on previous performance (~50k steps/second) - estimated_steps = int(time_budget_seconds * 50000) - - # Run fully JIT-compiled SA - start_time = time.perf_counter() - best_cost, best_route = jit_sa_optimized( - cities_coords, - params['initial_temp'], - params['cooling_rate'], - estimated_steps, - seed - ) - actual_time = time.perf_counter() - start_time - - # Calculate actual steps per second - actual_steps = int(estimated_steps * (time_budget_seconds / actual_time)) - - return best_cost, actual_steps - - -def evaluate_sa_params(params, instance, optimal_cost, time_budget_seconds=20.0, n_runs=2): - """ - Evaluate SA parameters by running multiple times sequentially with time budget. - - Args: - params: Dictionary of SA parameters - instance: TSP instance - optimal_cost: Optimal cost for reference - time_budget_seconds: Time budget per evaluation in seconds - n_runs: Number of runs to average - - Returns: - Average final cost across runs - """ - # Validate parameter combinations - initial_temp = params['initial_temp'] - cooling_rate = params['cooling_rate'] - - # Ensure reasonable parameter ranges - if initial_temp <= 0 or initial_temp > 1000: - return 100000.0 # Large penalty for invalid temperature - - if cooling_rate <= 0 or cooling_rate >= 1: - return 100000.0 # Large penalty for invalid cooling rate - - # Run sequentially (no parallel processing) - costs = [] - steps = [] - - for run in range(n_runs): - cost, step_count = run_single_sa_optimized(params, instance, time_budget_seconds, 42 + run) - costs.append(cost) - steps.append(step_count) - - # Return average cost (lower is better for minimization) - avg_cost = np.mean(costs) - avg_steps = np.mean(steps) - print(f"Params: {params} -> Avg cost: {avg_cost:.2f} (runs: {costs}, avg steps: {avg_steps:.0f})") - return avg_cost - - -def main(): - """Main hyperparameter optimization.""" - print("SA Hyperparameter Tuning") - print("=" * 50) - - # Load TSP instance - instance_path = Path("dataset/lin105.tsp") - instance, optimal_cost = find_optimal_tour(instance_path) - print(f"Instance: {instance.name}") - print(f"Optimal cost: {optimal_cost:.2f}") - print(f"Cities: {len(instance.cities)}") - print() - - # Define parameter search space - dimensions = [ - Real(1, 5000, name='initial_temp'), # Initial temperature - Real(0.95, 0.9999, name='cooling_rate'), # Cooling rate - ] - - # Objective function wrapper - @use_named_args(dimensions=dimensions) - def objective(**params): - return evaluate_sa_params(params, instance, optimal_cost, time_budget_seconds=20.0) - - print("Starting optimization...") - print("Each evaluation runs SA for 20 seconds, 2 times sequentially") - print("Using realistic parameter ranges with validation") - print() - - # Estimate time with a single evaluation - print("Estimating time with a single evaluation...") - test_start = time.time() - test_params = { - 'initial_temp': 100, - 'cooling_rate': 0.995 - } - evaluate_sa_params(test_params, instance, optimal_cost, time_budget_seconds=20.0) - test_time = time.time() - test_start - - estimated_total_time = test_time * 50 # 50 evaluations - print(f"Single evaluation took: {test_time:.1f} seconds") - print(f"Estimated total time: {estimated_total_time:.1f} seconds ({estimated_total_time/60:.1f} minutes)") - print() - - # Ask user if they want to continue - response = input("Continue with optimization? (y/n): ").lower().strip() - if response != 'y': - print("Optimization cancelled.") - return - - # Run Bayesian optimization - start_time = time.time() - result = gp_minimize( - func=objective, - dimensions=dimensions, - n_calls=50, # Number of evaluations - random_state=42, - acq_func='EI' # Expected Improvement - ) - end_time = time.time() - - # Results - print("\n" + "=" * 50) - print("OPTIMIZATION RESULTS") - print("=" * 50) - print(f"Optimization time: {end_time - start_time:.1f} seconds") - print(f"Best cost found: {result.fun:.2f}") - print(f"Optimality gap: {((result.fun / optimal_cost - 1) * 100):.1f}%") - print() - - print("Best parameters:") - best_params = dict(zip([dim.name for dim in dimensions], result.x)) - for param, value in best_params.items(): - print(f" {param}: {value:.3f}") - - print() - print("Parameter ranges tested:") - for dim in dimensions: - print(f" {dim.name}: {dim.low} to {dim.high}") - - # Test best parameters with more runs - print("\n" + "=" * 50) - print("VALIDATION WITH BEST PARAMETERS") - print("=" * 50) - - validation_costs = [] - validation_steps = [] - for run in range(5): # More runs for validation - sa = SimulatedAnnealing( - instance=instance, - start_temperature=best_params['initial_temp'], - cooling_schedule=exponential_cooling(best_params['cooling_rate']), - seed=100 + run - ) - - sa.initialize(None) - - # Run for 20 seconds - start_time = time.perf_counter() - steps = 0 - while time.perf_counter() - start_time < 20.0: - sa.step() - steps += 1 - - validation_costs.append(sa.best_cost) - validation_steps.append(steps) - print(f"Run {run+1}: {sa.best_cost:.2f} ({steps} steps)") - - print("\nValidation results:") - print(f" Mean cost: {np.mean(validation_costs):.2f}") - print(f" Std cost: {np.std(validation_costs):.2f}") - print(f" Best cost: {np.min(validation_costs):.2f}") - print(f" Worst cost: {np.max(validation_costs):.2f}") - print(f" Mean optimality gap: {((np.mean(validation_costs) / optimal_cost - 1) * 100):.1f}%") - print(f" Mean steps per second: {np.mean(validation_steps) / 20:.0f}") - print(f" Steps range: {np.min(validation_steps)} - {np.max(validation_steps)}") - - -if __name__ == "__main__": - main() diff --git a/tsp_analysis.ipynb b/tsp_analysis.ipynb deleted file mode 100644 index 0ebb710..0000000 --- a/tsp_analysis.ipynb +++ /dev/null @@ -1,479 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "0", - "metadata": {}, - "source": [ - "# Traveling Salesman Problem: Algorithm Comparison\n", - "\n", - "This notebook presents a comparison of various TSP optimization algorithms, evaluating their performance on the Lin105 benchmark dataset.\n", - "\n", - "**Algorithms Analyzed:**\n", - "- Random Solver (baseline)\n", - "- Nearest Neighbor\n", - "- Simulated Annealing (with different initialization strategies)\n", - "- Genetic Algorithm (with different initialization strategies)\n", - "\n", - "**Evaluation Metrics:**\n", - "- Solution quality (final cost vs optimal)\n", - "- Convergence speed (iterations to solution)\n", - "- Computational efficiency (steps per second)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import random\n", - "from pathlib import Path\n", - "from multiprocessing import Pool, cpu_count\n", - "\n", - "from algorithm.nearest_neighbor import NearestNeighbor\n", - "from algorithm.simulated_annealing import SimulatedAnnealing\n", - "from algorithm.genetic_algo import GeneticAlgorithmSolver\n", - "from algorithm.random_solver import RandomSolver\n", - "from util import exponential_cooling, find_optimal_tour, run_algorithm_with_timing, run_algorithm_with_iterations\n", - "\n", - "# Configuration\n", - "MAX_ITERATIONS = 1_000\n", - "MAX_SECONDS = 20.0\n", - "RANDOM_SEED = 42\n", - "COOLING_RATE = 0.9999\n", - "N_RUNS = 10\n", - "\n", - "ALGORITHMS = [\n", - " \"Random Solver\",\n", - " \"Nearest Neighbor\",\n", - " \"SA (NN-init)\",\n", - " \"SA (Random-init)\",\n", - " \"GA (NN-init)\",\n", - " \"GA (Random-init)\"\n", - "]\n", - "\n", - "random.seed(RANDOM_SEED)\n", - "\n", - "# Plot styling\n", - "plt.rcParams['figure.figsize'] = (14, 8)\n", - "plt.rcParams['font.size'] = 11\n", - "plt.rcParams['axes.labelsize'] = 12\n", - "plt.rcParams['axes.titlesize'] = 14\n", - "plt.rcParams['legend.fontsize'] = 10" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2", - "metadata": {}, - "outputs": [], - "source": [ - "problem_instance_path = Path(\"dataset/lin105.tsp\")\n", - "instance, optimal_cost = find_optimal_tour(problem_instance_path)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3", - "metadata": {}, - "outputs": [], - "source": [ - "seed_random = list(range(len(instance.cities)))\n", - "random.shuffle(seed_random)\n", - "\n", - "nn_builder = NearestNeighbor(instance)\n", - "nn_builder.initialize()\n", - "for _ in range(len(instance.cities) - 1):\n", - " _ = nn_builder.step()\n", - "seed_nn = nn_builder.get_route()\n", - "\n", - "# Algorithm configuration\n", - "T0 = 935\n", - "exp_schedule = exponential_cooling(COOLING_RATE)\n", - "\n", - "# GA parameters optimized through systematic hyperparameter tuning\n", - "GA_POP_SIZE = 30\n", - "GA_MUTATION = 0.3\n", - "GA_CROSSOVER = 0.6\n", - "GA_ELITISM = 1\n", - "\n", - "# Algorithm factory for parallel execution\n", - "def create_algorithm_instance(name, instance, run_idx, seed_nn, seed_identity):\n", - " \"\"\"Factory function to create algorithm instances with proper seeding.\"\"\"\n", - " T0 = 100\n", - " exp_schedule = exponential_cooling(COOLING_RATE)\n", - " seed = RANDOM_SEED + run_idx\n", - " \n", - " algorithm_configs = {\n", - " \"Random Solver\": (RandomSolver(instance, seed=seed), None),\n", - " \"Nearest Neighbor\": (NearestNeighbor(instance), None),\n", - " \"SA (NN-init)\": (SimulatedAnnealing(instance, T0, exp_schedule, seed=seed), seed_nn),\n", - " \"SA (Random-init)\": (SimulatedAnnealing(instance, T0, exp_schedule, seed=seed), None),\n", - " \"GA (NN-init)\": (\n", - " GeneticAlgorithmSolver(instance, seed=seed, population_size=GA_POP_SIZE,\n", - " mutation_rate=GA_MUTATION, crossover_rate=GA_CROSSOVER,\n", - " num_parents=2, num_child=2, elitism_count=GA_ELITISM),\n", - " seed_nn\n", - " ),\n", - " \"GA (Random-init)\": (\n", - " GeneticAlgorithmSolver(instance, seed=seed, population_size=GA_POP_SIZE,\n", - " mutation_rate=GA_MUTATION, crossover_rate=GA_CROSSOVER,\n", - " num_parents=2, num_child=2, elitism_count=GA_ELITISM),\n", - " None\n", - " ),\n", - " }\n", - " \n", - " return algorithm_configs[name]\n", - "\n", - "# Helper function for parallel execution\n", - "def run_single_trial(args):\n", - " \"\"\"Run a single algorithm trial - used for parallel processing.\"\"\"\n", - " name, run_idx, instance_data, seed_nn_data, seed_identity_data = args\n", - " \n", - " # Recreate instance from data (needed for multiprocessing)\n", - " from tsp.model import TSPInstance\n", - " instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities'])\n", - " \n", - " solver, init_route = create_algorithm_instance(name, instance, run_idx, seed_nn_data, seed_identity_data)\n", - " \n", - " iterations, best_costs, current_costs, times, best_route = run_algorithm_with_timing(\n", - " instance, solver, init_route, MAX_SECONDS\n", - " )\n", - " \n", - " return {\n", - " 'iterations': iterations,\n", - " 'best_costs': best_costs,\n", - " 'current_costs': current_costs,\n", - " 'times': times,\n", - " 'best_route': best_route\n", - " }\n", - "\n", - "# Time-based benchmark with multiple runs (parallelized)\n", - "print(\"=\" * 70)\n", - "print(f\"TIME-BASED BENCHMARK ({MAX_SECONDS}s per algorithm, {N_RUNS} runs, {cpu_count()} CPUs)\")\n", - "print(\"=\" * 70)\n", - "\n", - "# Prepare instance data for serialization\n", - "instance_data = {'name': instance.name, 'cities': instance.cities}\n", - "seed_nn_data = seed_nn\n", - "\n", - "time_results = {}\n", - "for name in ALGORITHMS:\n", - " print(f\"Running {name}... ({N_RUNS} runs in parallel)\", end=\" \", flush=True)\n", - " \n", - " # Prepare arguments for parallel execution\n", - " args_list = [(name, run_idx, instance_data, seed_nn_data, seed_random) \n", - " for run_idx in range(N_RUNS)]\n", - " \n", - " # Run in parallel\n", - " with Pool(processes=min(cpu_count(), N_RUNS)) as pool:\n", - " all_runs = pool.map(run_single_trial, args_list)\n", - " \n", - " # Align all runs to common time grid using interpolation\n", - " # Find the run that reached closest to MAX_SECONDS\n", - " max_time_reached = max(run['times'][-1] for run in all_runs)\n", - " \n", - " # Create uniform time grid from 0 to the maximum time reached\n", - " num_points = 100 # Use 100 points for smooth curves\n", - " common_times = np.linspace(0, min(max_time_reached, MAX_SECONDS), num_points)\n", - " \n", - " # Interpolate each run onto common time grid\n", - " aligned_best = []\n", - " for run in all_runs:\n", - " # Interpolate best costs onto common time grid\n", - " interp_best = np.interp(common_times, run['times'], run['best_costs'])\n", - " aligned_best.append(interp_best)\n", - " \n", - " aligned_best = np.array(aligned_best)\n", - " aligned_times = common_times\n", - " \n", - " # Calculate statistics\n", - " mean_best = np.mean(aligned_best, axis=0)\n", - " std_best = np.std(aligned_best, axis=0)\n", - " min_best = np.min(aligned_best, axis=0)\n", - " max_best = np.max(aligned_best, axis=0)\n", - " \n", - " final_costs = [run['best_costs'][-1] for run in all_runs]\n", - " \n", - " time_results[name] = {\n", - " 'times': aligned_times,\n", - " 'mean_best': mean_best,\n", - " 'std_best': std_best,\n", - " 'min_best': min_best,\n", - " 'max_best': max_best,\n", - " 'all_runs_best': aligned_best,\n", - " 'final_cost_mean': np.mean(final_costs),\n", - " 'final_cost_std': np.std(final_costs),\n", - " 'final_cost_min': np.min(final_costs),\n", - " 'final_cost_max': np.max(final_costs),\n", - " 'total_iterations': np.mean([len(run['iterations']) for run in all_runs])\n", - " }\n", - " \n", - " print(f\"Mean: {time_results[name]['final_cost_mean']:.2f} ± {time_results[name]['final_cost_std']:.2f}\")\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4", - "metadata": {}, - "outputs": [], - "source": [ - "# Helper function for iteration-based parallel execution\n", - "def run_single_iteration_trial(args):\n", - " \"\"\"Run a single iteration-based trial.\"\"\"\n", - " name, run_idx, instance_data, seed_nn_data, seed_identity_data = args\n", - " \n", - " from tsp.model import TSPInstance\n", - " instance = TSPInstance(name=instance_data['name'], cities=instance_data['cities'])\n", - " \n", - " solver, init_route = create_algorithm_instance(name, instance, run_idx, seed_nn_data, seed_identity_data)\n", - " \n", - " iters, best_costs, current_costs, times, best_route = run_algorithm_with_iterations(\n", - " instance, solver, init_route, MAX_ITERATIONS\n", - " )\n", - " \n", - " return {'iters': iters, 'best_costs': best_costs, 'times': times, 'best_route': best_route}\n", - "\n", - "# Iteration-based benchmark with multiple runs (parallelized)\n", - "print(\"\\n\" + \"=\" * 70)\n", - "print(f\"ITERATION-BASED BENCHMARK ({MAX_ITERATIONS} iterations, {N_RUNS} runs)\")\n", - "print(\"=\" * 70)\n", - "\n", - "iteration_stats = {}\n", - "for name in time_results.keys():\n", - " print(f\"Running {name}... ({N_RUNS} runs in parallel)\", end=\" \", flush=True)\n", - " \n", - " args_list = [(name, run_idx, instance_data, seed_nn_data, seed_random) \n", - " for run_idx in range(N_RUNS)]\n", - " \n", - " with Pool(processes=min(cpu_count(), N_RUNS)) as pool:\n", - " results = pool.map(run_single_iteration_trial, args_list)\n", - " \n", - " runs_best = [r['best_costs'] for r in results]\n", - " runs_iters = [r['iters'] for r in results]\n", - " runs_time = [r['times'] for r in results]\n", - " \n", - " # Align by min length\n", - " min_len = min(len(x) for x in runs_best)\n", - " aligned_best = [run[:min_len] for run in runs_best]\n", - " aligned_iters = runs_iters[0][:min_len]\n", - " \n", - " mean_best = np.mean(aligned_best, axis=0)\n", - " std_best = np.std(aligned_best, axis=0)\n", - " final_costs = [run[-1] for run in runs_best if run]\n", - " total_times = [t[-1] for t in runs_time if t]\n", - " \n", - " iteration_stats[name] = {\n", - " 'iterations': aligned_iters,\n", - " 'mean_best': mean_best,\n", - " 'std_best': std_best,\n", - " 'final_cost_mean': np.mean(final_costs) if final_costs else float('inf'),\n", - " 'final_cost_std': np.std(final_costs) if final_costs else 0.0,\n", - " 'total_time_mean': np.mean(total_times) if total_times else 0.0,\n", - " }\n", - " \n", - " print(f\"Cost: {iteration_stats[name]['final_cost_mean']:.2f} ± {iteration_stats[name]['final_cost_std']:.2f}\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5", - "metadata": {}, - "outputs": [], - "source": [ - "# Summary Table\n", - "print(\"\\n\" + \"=\" * 85)\n", - "print(\"SUMMARY\")\n", - "print(\"=\" * 85)\n", - "print(f\"{'Algorithm':<20} {'Mean Cost':<18} {'vs Optimal':<12} {'Steps/sec':<12} {'CV %':<8}\")\n", - "print(\"-\" * 85)\n", - "\n", - "for name in time_results.keys():\n", - " mean_cost = time_results[name]['final_cost_mean']\n", - " std_cost = time_results[name]['final_cost_std']\n", - " cost_str = f\"{mean_cost:.1f} ± {std_cost:.1f}\"\n", - " vs_optimal = f\"+{((mean_cost / optimal_cost - 1) * 100):.1f}%\" if optimal_cost else \"N/A\"\n", - " steps_per_sec = time_results[name]['total_iterations']\n", - " cv = (std_cost / mean_cost * 100)\n", - " print(f\"{name:<20} {cost_str:<18} {vs_optimal:<12} {steps_per_sec:<12.0f} {cv:<8.1f}\")\n", - "\n", - "print(f\"{'Optimal':<20} {optimal_cost:<18.2f} {'0.0%':<12} {'-':<12} {'-':<8}\")\n", - "print(\"=\" * 85)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6", - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(14, 8))\n", - "\n", - "colors = plt.cm.tab10(np.linspace(0, 1, len(time_results)))\n", - "for (name, data), color in zip(time_results.items(), colors):\n", - " # Plot mean line\n", - " ax.plot(data['times'], data['mean_best'], label=name, linewidth=2.5, color=color, alpha=0.9)\n", - " \n", - " # Add shaded region showing min-max bounds across runs\n", - " ax.fill_between(data['times'], data['min_best'], data['max_best'], \n", - " color=color, alpha=0.15, linewidth=0)\n", - "\n", - "if optimal_cost:\n", - " ax.axhline(y=optimal_cost, color='darkgreen', linestyle='--', linewidth=2, alpha=0.7, label='Optimal Solution')\n", - "\n", - "ax.set_xlabel('Time (seconds)', fontsize=13, fontweight='bold')\n", - "ax.set_ylabel('Best Cost Found (mean, shaded: min-max)', fontsize=13, fontweight='bold')\n", - "ax.set_title(f'Algorithm Convergence Over Time ({N_RUNS} runs) - {instance.name.upper()}', \n", - " fontsize=15, fontweight='bold', pad=20)\n", - "ax.legend(loc='upper right', framealpha=0.9, fontsize=11)\n", - "ax.grid(True, alpha=0.3, linestyle='--')\n", - "plt.tight_layout()\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7", - "metadata": {}, - "outputs": [], - "source": [ - "# --- Solution Quality Comparison Plot ---\n", - "fig1, ax1 = plt.subplots(figsize=(8, 7))\n", - "\n", - "names = list(time_results.keys())\n", - "final_costs_mean = [time_results[name]['final_cost_mean'] for name in names]\n", - "final_costs_std = [time_results[name]['final_cost_std'] for name in names]\n", - "\n", - "colors = plt.cm.tab10(np.linspace(0, 1, len(names)))\n", - "bars1 = ax1.bar(range(len(names)), final_costs_mean, yerr=final_costs_std, \n", - " capsize=5, color=colors, alpha=0.8, edgecolor='black')\n", - "\n", - "if optimal_cost:\n", - " ax1.axhline(y=optimal_cost, color='darkgreen', linestyle='--', linewidth=2, alpha=0.7, label='Optimal')\n", - " ax1.legend()\n", - "\n", - "ax1.set_xticks(range(len(names)))\n", - "ax1.set_xticklabels(names, rotation=45, ha='right')\n", - "ax1.set_ylabel('Final Cost (mean ± std)', fontsize=12, fontweight='bold')\n", - "ax1.set_title(f'Solution Quality Comparison ({N_RUNS} runs)', fontsize=14, fontweight='bold', pad=15)\n", - "ax1.grid(True, alpha=0.3, axis='y', linestyle='--')\n", - "\n", - "for i, (cost, err) in enumerate(zip(final_costs_mean, final_costs_std)):\n", - " ax1.text(i, cost + err + max(final_costs_mean) * 0.02, f'{cost:.0f}', \n", - " ha='center', va='bottom', fontsize=9, fontweight='bold')\n", - "\n", - "plt.tight_layout()\n", - "plt.show()\n", - "\n", - "# --- Computational Efficiency Plot ---\n", - "fig2, ax2 = plt.subplots(figsize=(8, 7))\n", - "\n", - "steps_per_sec = [time_results[name]['total_iterations'] / MAX_SECONDS for name in names]\n", - "bars2 = ax2.bar(range(len(names)), steps_per_sec, color=colors, alpha=0.8, edgecolor='black')\n", - "\n", - "ax2.set_xticks(range(len(names)))\n", - "ax2.set_xticklabels(names, rotation=45, ha='right')\n", - "ax2.set_ylabel('Steps per Second', fontsize=12, fontweight='bold')\n", - "ax2.set_title('Computational Efficiency', fontsize=14, fontweight='bold', pad=15)\n", - "ax2.grid(True, alpha=0.3, axis='y', linestyle='--')\n", - "\n", - "for i, rate in enumerate(steps_per_sec):\n", - " ax2.text(i, rate + max(steps_per_sec) * 0.02, f'{rate:.0f}', \n", - " ha='center', va='bottom', fontsize=9, fontweight='bold')\n", - "\n", - "plt.tight_layout()\n", - "plt.show()\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8", - "metadata": {}, - "outputs": [], - "source": [ - "# Performance Analysis Summary\n", - "print(\"\\nKEY FINDINGS:\\n\")\n", - "print(\"1. SOLUTION QUALITY:\")\n", - "\n", - "best_algorithm = min(time_results.items(), key=lambda x: x[1]['final_cost_mean'])\n", - "print(f\" - Best performer: {best_algorithm[0]}\")\n", - "print(f\" - Final cost: {best_algorithm[1]['final_cost_mean']:.2f} ± {best_algorithm[1]['final_cost_std']:.2f}\")\n", - "if optimal_cost:\n", - " gap = ((best_algorithm[1]['final_cost_mean'] / optimal_cost - 1) * 100)\n", - " print(f\" - Optimality gap: {gap:.1f}%\")\n", - "print(f\" - Best run achieved: {best_algorithm[1]['final_cost_min']:.2f}\")\n", - "print(f\" - Worst run: {best_algorithm[1]['final_cost_max']:.2f}\")\n", - "\n", - "print(\"\\n2. COMPUTATIONAL EFFICIENCY:\")\n", - "fastest = max(time_results.items(), key=lambda x: x[1]['total_iterations'])\n", - "print(f\" - Fastest algorithm: {fastest[0]}\")\n", - "print(f\" - Steps per second: {fastest[1]['total_iterations']:.0f}\")\n", - "\n", - "print(\"\\n3. CONVERGENCE SPEED:\")\n", - "for name in ['SA (NN-init)', 'GA (NN-init)', 'Nearest Neighbor']:\n", - " if name in time_results:\n", - " data = time_results[name]\n", - " if len(data['mean_best']) > 10:\n", - " # Check how quickly it gets close to final solution\n", - " final = data['mean_best'][-1]\n", - " threshold = final * 1.1 # Within 10% of final\n", - " converged_idx = next((i for i, cost in enumerate(data['mean_best']) if cost <= threshold), len(data['mean_best']))\n", - " converge_time = data['times'][converged_idx] if converged_idx < len(data['times']) else data['times'][-1]\n", - " print(f\" - {name}: converged in {converge_time:.3f}s (average)\")\n", - "\n", - "print(f\"\\n4. VARIABILITY ACROSS RUNS:\")\n", - "for name in time_results.keys():\n", - " variance = time_results[name]['final_cost_std'] / time_results[name]['final_cost_mean'] * 100\n", - " print(f\" - {name}: CV = {variance:.1f}%\")" - ] - }, - { - "cell_type": "markdown", - "id": "9", - "metadata": {}, - "source": [ - "## Conclusions\n", - "\n", - "Based on comprehensive benchmarking on Lin105 (10 runs per algorithm, 5s runtime):\n", - "\n", - "Both SA variants achieve ~9% optimality gap with high throughput (136k–146k steps/sec) and low variability (CV ~3%). Random initialization performs slightly better than NN-init, suggesting NN seeding less critical for SA on larger instances.\n", - "\n", - "GA (NN-init) achieves only 18% gap despite NN seeding. GA (Random-init) fails to beat Nearest Neighbor baseline and shows high instability (CV 9.4%). Both GA variants are 100× slower than SA due to population overhead.\n", - "\n", - "Initialization affects GA highly (NN-init reduces gap from 58% to 18%) but minimal for SA (both variants ~9%) in this example." - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": ".venv", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.11" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/tuning/ga_tuning.py b/tuning/ga_tuning.py new file mode 100644 index 0000000..b5667b9 --- /dev/null +++ b/tuning/ga_tuning.py @@ -0,0 +1,138 @@ +import numpy as np +from pathlib import Path +import sys +sys.path.insert(0, str(Path(__file__).parent.parent)) +from skopt import gp_minimize +from skopt.space import Real, Integer +from skopt.utils import use_named_args +import time +from concurrent.futures import ProcessPoolExecutor, as_completed + +from util import find_optimal_tour +from algorithm.genetic_algo import GeneticAlgorithmSolver + + +def run_single_ga(params: dict, instance, time_budget_seconds: float) -> tuple[float, int]: + """ + Runs a single, randomized GA instance for a fixed time budget. + """ + ga = GeneticAlgorithmSolver( + instance=instance, + population_size=int(params['population_size']), + crossover_rate=params['crossover_rate'], + mutation_rate=params['mutation_rate'], + elitism_count=int(params['elitism_count']) + ) + ga.initialize(None) + + start_time = time.perf_counter() + steps = 0 + while time.perf_counter() - start_time < time_budget_seconds: + ga.step() + steps += 1 + + return ga.best_cost, steps + + +def evaluate_ga_params_parallel(params: dict, instance, time_budget_seconds: float, n_runs: int) -> float: + """ + Evaluates GA parameters by running multiple instances in parallel and returning the average cost. + """ + population_size = int(params['population_size']) + elitism_count = int(params['elitism_count']) + + if elitism_count >= population_size or elitism_count > population_size * 0.5: + return 1e9 + + costs = [] + with ProcessPoolExecutor() as executor: + # Submit the same task n_runs times. This is the correct way to call a function + # with fixed arguments multiple times in parallel. + futures = [executor.submit(run_single_ga, params, instance, time_budget_seconds) for _ in range(n_runs)] + + # Collect results as they complete. + for future in as_completed(futures): + cost, _ = future.result() + costs.append(cost) + + avg_cost = np.mean(costs) + print(f"Params: pop={population_size}, cross={params['crossover_rate']:.2f}, " + f"mute={params['mutation_rate']:.2f}, elite={elitism_count} -> Avg cost: {avg_cost:.2f}") + + return avg_cost + + +def main(): + ### CONFIGURATION ### + TIME_BUDGET_PER_EVALUATION_SECONDS = 5.0 + TOTAL_TUNING_TIME_MINUTES = 10.0 + N_RUNS_PER_EVALUATION = 4 + ### END CONFIGURATION ### + + instance_path = Path("dataset/lin105.tsp") + instance, optimal_cost = find_optimal_tour(instance_path) + print(f"Optimizing GA for instance: {instance.name} (Optimal: {optimal_cost:.2f})") + + time_per_optimizer_call = TIME_BUDGET_PER_EVALUATION_SECONDS + n_calls = max(1, int((TOTAL_TUNING_TIME_MINUTES * 60) / time_per_optimizer_call)) + + print("\n" + "-"*50) + print("Tuning Configuration:") + print(f" - Time per GA evaluation: {TIME_BUDGET_PER_EVALUATION_SECONDS} seconds") + print(f" - Parallel runs per param set: {N_RUNS_PER_EVALUATION}") + print(f" - Total tuning goal: ~{TOTAL_TUNING_TIME_MINUTES} minutes") + print(f" - Calculated optimizer calls: {n_calls}") + print("-"*50 + "\n") + + search_space = [ + Integer(30, 250, name='population_size'), + Real(0.6, 0.95, name='crossover_rate'), + Real(0.01, 0.4, name='mutation_rate'), + Integer(1, 25, name='elitism_count'), + ] + + @use_named_args(dimensions=search_space) + def objective(**params): + return evaluate_ga_params_parallel( + params, instance, TIME_BUDGET_PER_EVALUATION_SECONDS, N_RUNS_PER_EVALUATION + ) + + start_time = time.time() + result = gp_minimize( + func=objective, + dimensions=search_space, + n_calls=n_calls, + acq_func='EI' + ) + end_time = time.time() + + print("\n" + "=" * 50) + print("OPTIMIZATION COMPLETE") + print(f"Total time: {(end_time - start_time)/60:.2f} minutes") + + best_params = dict(zip([dim.name for dim in search_space], result.x)) + print(f"Best average cost found: {result.fun:.2f}") + print("Best parameters:") + print(f" population_size: {int(best_params['population_size'])}") + print(f" crossover_rate: {best_params['crossover_rate']:.3f}") + print(f" mutation_rate: {best_params['mutation_rate']:.3f}") + print(f" elitism_count: {int(best_params['elitism_count'])}") + print("=" * 50) + + print("\nValidating best parameters with 5 runs...") + validation_costs = [] + for i in range(5): + # The validation runs are also randomized now. + cost, steps = run_single_ga( + best_params, instance, time_budget_seconds=TIME_BUDGET_PER_EVALUATION_SECONDS + ) + validation_costs.append(cost) + print(f" Run {i+1}: Cost={cost:.2f} ({steps} steps)") + + mean_cost = np.mean(validation_costs) + gap = ((mean_cost / optimal_cost) - 1) * 100 + print(f"\nValidation Mean Cost: {mean_cost:.2f} (Optimality Gap: {gap:.2f}%)") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/tuning/sa_tuning.py b/tuning/sa_tuning.py new file mode 100644 index 0000000..97169bb --- /dev/null +++ b/tuning/sa_tuning.py @@ -0,0 +1,126 @@ +import numpy as np +from pathlib import Path +import sys +sys.path.insert(0, str(Path(__file__).parent.parent)) +from skopt import gp_minimize +from skopt.space import Real +from skopt.utils import use_named_args +import time +from concurrent.futures import ProcessPoolExecutor, as_completed + +from util import find_optimal_tour, exponential_cooling +from algorithm.simulated_annealing import SimulatedAnnealing + + +def run_single_sa(params: dict, instance, time_budget_seconds: float) -> tuple[float, int]: + """ + Runs a single, randomized SA instance for a fixed time budget. + """ + sa = SimulatedAnnealing( + instance=instance, + start_temperature=params['initial_temp'], + cooling_schedule=exponential_cooling(params['cooling_rate']), + seed=None # Ensures each parallel run is unique + ) + sa.initialize(None) + + start_time = time.perf_counter() + steps = 0 + while time.perf_counter() - start_time < time_budget_seconds: + sa.step() + steps += 1 + + return sa.best_cost, steps + + +def evaluate_sa_params_parallel(params: dict, instance, time_budget_seconds: float, n_runs: int) -> float: + """ + Evaluates SA parameters by running multiple instances in parallel and returning the average cost. + """ + # Basic parameter validation + if not (0.9 < params['cooling_rate'] < 1.0) or not (params['initial_temp'] > 0): + return 1e9 + + costs = [] + with ProcessPoolExecutor() as executor: + futures = [executor.submit(run_single_sa, params, instance, time_budget_seconds) for _ in range(n_runs)] + + for future in as_completed(futures): + cost, _ = future.result() + costs.append(cost) + + avg_cost = np.mean(costs) + print(f"Params: temp={params['initial_temp']:.2f}, cool={params['cooling_rate']:.4f} -> Avg cost: {avg_cost:.2f}") + + return avg_cost + + +def main(): + ### CONFIGURATION ### + TIME_BUDGET_PER_EVALUATION_SECONDS = 5.0 + TOTAL_TUNING_TIME_MINUTES = 10.0 + N_RUNS_PER_EVALUATION = 4 + ### END CONFIGURATION ### + + instance_path = Path("dataset/lin105.tsp") + instance, optimal_cost = find_optimal_tour(instance_path) + print(f"Optimizing SA for instance: {instance.name} (Optimal: {optimal_cost:.2f})") + + n_calls = max(1, int((TOTAL_TUNING_TIME_MINUTES * 60) / TIME_BUDGET_PER_EVALUATION_SECONDS)) + + print("\n" + "-"*50) + print("Tuning Configuration:") + print(f" - Time per SA evaluation: {TIME_BUDGET_PER_EVALUATION_SECONDS} seconds") + print(f" - Parallel runs per param set: {N_RUNS_PER_EVALUATION}") + print(f" - Total tuning goal: ~{TOTAL_TUNING_TIME_MINUTES} minutes") + print(f" - Calculated optimizer calls: {n_calls}") + print("-"*50 + "\n") + + search_space = [ + Real(1.0, 5000.0, name='initial_temp', prior='log-uniform'), + Real(0.99, 0.9999, name='cooling_rate'), + ] + + @use_named_args(dimensions=search_space) + def objective(**params): + return evaluate_sa_params_parallel( + params, instance, TIME_BUDGET_PER_EVALUATION_SECONDS, N_RUNS_PER_EVALUATION + ) + + start_time = time.time() + result = gp_minimize( + func=objective, + dimensions=search_space, + n_calls=n_calls, + random_state=42, + acq_func='EI' + ) + end_time = time.time() + + print("\n" + "=" * 50) + print("OPTIMIZATION COMPLETE") + print(f"Total time: {(end_time - start_time)/60:.2f} minutes") + + best_params = dict(zip([dim.name for dim in search_space], result.x)) + print(f"Best average cost found: {result.fun:.2f}") + print("Best parameters:") + print(f" initial_temp: {best_params['initial_temp']:.3f}") + print(f" cooling_rate: {best_params['cooling_rate']:.5f}") + print("=" * 50) + + print("\nValidating best parameters with 5 runs...") + validation_costs = [] + for i in range(5): + cost, steps = run_single_sa( + best_params, instance, time_budget_seconds=TIME_BUDGET_PER_EVALUATION_SECONDS + ) + validation_costs.append(cost) + print(f" Run {i+1}: Cost={cost:.2f} ({steps} steps)") + + mean_cost = np.mean(validation_costs) + gap = ((mean_cost / optimal_cost) - 1) * 100 + print(f"\nValidation Mean Cost: {mean_cost:.2f} (Optimality Gap: {gap:.2f}%)") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/tuning/tuning_10.json b/tuning/tuning_10.json new file mode 100644 index 0000000..c68bd73 --- /dev/null +++ b/tuning/tuning_10.json @@ -0,0 +1,13 @@ +{ + "sa_tuning": { + "initial_temp": 4531.442, + "cooling_rate": 0.99990 + }, + "ga_tuning": { + "population_size": 63, + "crossover_rate": 0.794, + "mutation_rate": 0.378, + "elitism_count": 9 + } +} + diff --git a/tuning/tuning_5.json b/tuning/tuning_5.json new file mode 100644 index 0000000..0977e19 --- /dev/null +++ b/tuning/tuning_5.json @@ -0,0 +1,12 @@ +{ + "sa_tuning": { + "initial_temp": 167.807, + "cooling_rate": 0.99990 + }, + "ga_tuning": { + "population_size": 109, + "crossover_rate": 0.600, + "mutation_rate": 0.400, + "elitism_count": 1 + } +} diff --git a/util.py b/util.py index 3c752b5..05c3336 100644 --- a/util.py +++ b/util.py @@ -2,8 +2,7 @@ from pathlib import Path import time import logging -import random -from constants import MAX_SECONDS, MAX_ITERATIONS, COOLING_RATE, POPULATION_SIZE, MUTATION_RATE, CROSSOVER_RATE, ELITISM_COUNT, NUM_PARENTS, NUM_CHILD +from constants import MAX_SECONDS, MAX_ITERATIONS, COOLING_RATE, POPULATION_SIZE, MUTATION_RATE, CROSSOVER_RATE, ELITISM_COUNT # Import algorithms from algorithm.nearest_neighbor import NearestNeighbor @@ -14,8 +13,6 @@ logger = logging.getLogger(__name__) def setup_algorithm(alg_name: str, instance): - RANDOM_SEED = time.time() # For variability in real runs - random.seed(RANDOM_SEED) # Algorithm Specific Setup init_route = None match alg_name: @@ -25,7 +22,6 @@ def setup_algorithm(alg_name: str, instance): instance, INITIAL_TEMP, exponential_cooling(COOLING_RATE), # Cooling Schedule - seed=RANDOM_SEED ) case "SimulatedAnnealing_NearestNeighbor": INITIAL_TEMP = 100 @@ -34,7 +30,6 @@ def setup_algorithm(alg_name: str, instance): instance, INITIAL_TEMP, exponential_cooling(COOLING_RATE), # Cooling Schedule - seed=RANDOM_SEED ) case "GeneticAlgorithm_NearestNeighbor": solver = GeneticAlgorithmSolver( @@ -43,14 +38,19 @@ def setup_algorithm(alg_name: str, instance): mutation_rate=MUTATION_RATE, crossover_rate=CROSSOVER_RATE, elitism_count=ELITISM_COUNT, - num_parents=NUM_PARENTS, - num_child=NUM_CHILD, - seed=RANDOM_SEED + ) + case "GeneticAlgorithm": + solver = GeneticAlgorithmSolver( + instance, + population_size=POPULATION_SIZE, + mutation_rate=MUTATION_RATE, + crossover_rate=CROSSOVER_RATE, + elitism_count=ELITISM_COUNT, ) case "Baseline_Random": - solver = RandomSolver(instance, seed=RANDOM_SEED) + solver = RandomSolver(instance) case _: - logger.info(f" Only the following algorithm names are supported:\nSimulatedAnnealing_random, SimulatedAnnealing_NearestNeighbor, GeneticAlgorithm, Baseline_Random.") + logger.info(" Only the following algorithm names are supported:\nSimulatedAnnealing_random, SimulatedAnnealing_NearestNeighbor, GeneticAlgorithm, Baseline_Random.") raise ValueError(f"[ERROR] Unknown algorithm name: {alg_name}") # Update seed route for Nearest Neighbor based algorithms @@ -76,13 +76,11 @@ def find_optimal_tour(tsp_path: str | Path): if isinstance(tsp_path, str): tsp_path = Path(tsp_path) instance = parse_tsplib_tsp(tsp_path) - logger.info(f"Loaded {instance.name} with {len(instance.cities)} cities") # Check for optimal tour file opt_tour_path = Path(f"dataset/{instance.name}.opt.tour") optimal_cost = None if opt_tour_path.exists(): - logger.info(f"Found optimal tour file: {opt_tour_path}") # Parse optimal tour (simple format: just city indices) with open(opt_tour_path, "r") as f: lines = f.readlines() @@ -107,11 +105,10 @@ def find_optimal_tour(tsp_path: str | Path): if tour: optimal_cost = instance.route_cost(tour) - logger.info(f"Optimal cost: {optimal_cost:.2f}") else: logger.warning("Could not parse optimal tour") else: - logger.info("No optimal tour file found") + logger.warning(f"No optimal tour file found for {instance.name}") return instance, optimal_cost diff --git a/uv.lock b/uv.lock index 9a8a0aa..cc27626 100644 --- a/uv.lock +++ b/uv.lock @@ -2103,6 +2103,18 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/5e/4f/e1f65e8f8c76d73658b33d33b81eed4322fb5085350e4328d5c956f0c8f9/tornado-6.5.2-cp39-abi3-win_arm64.whl", hash = "sha256:d6c33dc3672e3a1f3618eb63b7ef4683a7688e7b9e6e8f0d9aa5726360a004af", size = 444456, upload-time = "2025-08-08T18:26:59.207Z" }, ] +[[package]] +name = "tqdm" +version = "4.67.1" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "colorama", marker = "sys_platform == 'win32'" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/a8/4b/29b4ef32e036bb34e4ab51796dd745cdba7ed47ad142a9f4a1eb8e0c744d/tqdm-4.67.1.tar.gz", hash = "sha256:f8aef9c52c08c13a65f30ea34f4e5aac3fd1a34959879d7e59e63027286627f2", size = 169737, upload-time = "2024-11-24T20:12:22.481Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/d0/30/dc54f88dd4a2b5dc8a0279bdd7270e735851848b762aeb1c1184ed1f6b14/tqdm-4.67.1-py3-none-any.whl", hash = "sha256:26445eca388f82e72884e0d580d5464cd801a3ea01e63e5601bdff9ba6a48de2", size = 78540, upload-time = "2024-11-24T20:12:19.698Z" }, +] + [[package]] name = "traitlets" version = "5.14.3" @@ -2129,6 +2141,7 @@ dependencies = [ { name = "requests" }, { name = "scikit-optimize" }, { name = "seaborn" }, + { name = "tqdm" }, ] [package.dev-dependencies] @@ -2150,6 +2163,7 @@ requires-dist = [ { name = "requests", specifier = ">=2.31.0" }, { name = "scikit-optimize", specifier = ">=0.10.2" }, { name = "seaborn", specifier = ">=0.13.2" }, + { name = "tqdm", specifier = ">=4.67.1" }, ] [package.metadata.requires-dev]