diff --git a/.gitignore b/.gitignore index 3a9118a..12dd477 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ /results/ -/data/ \ No newline at end of file +**/data/ +**/__pycache__/ diff --git a/docs/source/index.rst b/docs/source/index.rst index aa5ddc6..ad05af8 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -13,6 +13,8 @@ RL4Ising overview/introduction overview/key_concepts overview/content + overview/hardware_reqs + .. toctree:: :maxdepth: 2 diff --git a/docs/source/overview/hardware_reqs.rst b/docs/source/overview/hardware_reqs.rst new file mode 100644 index 0000000..c40635e --- /dev/null +++ b/docs/source/overview/hardware_reqs.rst @@ -0,0 +1,10 @@ +# Software Requirements for FOCI Cluster + +A100 - CUDA 11.4 or later supported +tensorflow-2.12.0 (minimum version for CUDA 11.8) +python 3.8-3.11 +Minor compatibility is available for all versions within 11 meaning new code and old hardware versions should be supported with the right packages. + +Sources: + + diff --git a/src/baseline/.ipynb_checkpoints/gurobi_qubo-checkpoint.ipynb b/src/baseline/.ipynb_checkpoints/gurobi_qubo-checkpoint.ipynb new file mode 100644 index 0000000..c5a3060 --- /dev/null +++ b/src/baseline/.ipynb_checkpoints/gurobi_qubo-checkpoint.ipynb @@ -0,0 +1,177 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Prerequisites\n", + "Need to have Gurobi Optimizer and obtain license.\n", + "\n", + "1. Install Gurobi\n", + "- https://www.gurobi.com/downloads/gurobi-software/\n", + "- https://support.gurobi.com/hc/en-us/articles/4534161999889-How-do-I-install-Gurobi-Optimizer\n", + "2. Obtain License\n", + "- https://portal.gurobi.com/iam/licenses/list/\n", + "\n", + "\n", + "%pip install gurobipy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Imports\n", + "import os\n", + "import json\n", + "import numpy as np\n", + "import networkx as nx\n", + "from natsort import natsorted\n", + "import gurobipy as gp\n", + "from math import sqrt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def read_nxgraph(filename: str) -> nx.Graph(): # type: ignore\n", + " graph = nx.Graph()\n", + " with open(filename, 'r') as file:\n", + " # lines = []\n", + " line = file.readline()\n", + " is_first_line = True\n", + " while line is not None and line != '':\n", + " if '//' not in line:\n", + " if is_first_line:\n", + " strings = line.split(\" \")\n", + " num_nodes = int(strings[0])\n", + " num_edges = int(strings[1])\n", + " nodes = list(range(num_nodes))\n", + " graph.add_nodes_from(nodes)\n", + " is_first_line = False\n", + " else:\n", + " node1, node2, weight = line.split()\n", + " graph.add_edge(int(node1), int(node2), weight=weight)\n", + " line = file.readline()\n", + " return graph" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def gurobi_maxcut(graph):\n", + "\n", + " # Create QUBO matrix\n", + " nodes = len(list(graph.nodes))\n", + " J = nx.to_numpy_array(graph)\n", + "\n", + " # Construct gurobi model\n", + " model = gp.Model(\"maxcut_qubo\")\n", + "\n", + " # Set time limit 1 hr\n", + " model.Params.LogFile = \"../logs/gurobi.log\"\n", + " model.setParam('TimeLimit', 3600)\n", + " model.Params.LogToConsole = 1\n", + " model.setParam(\"MIPGap\", 0.0)\n", + "\n", + " # Create variable for each vertex\n", + " x = model.addVars(nodes, vtype=gp.GRB.BINARY)\n", + "\n", + " objective = gp.quicksum(\n", + " -J[i, j] * (2 * x[i] - 1) * (2 * x[j] - 1) * (1/sqrt(nodes))\n", + " for i in range(nodes) for j in range(i+1, nodes) if J[i, j] != 0.0\n", + " )\n", + " model.setObjective(objective, gp.GRB.MINIMIZE)\n", + " \n", + " # Solve\n", + " model.optimize()\n", + " obj_val = model.ObjVal\n", + " obj_bnd = model.ObjBound\n", + " solution = \"\"\n", + " for i in range(nodes):\n", + " solution += f\"{int(x[i].X)}\"\n", + "\n", + " return obj_val, solution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# graph = read_nxgraph(\"../data/vna/ER/100_SK_seed40.txt\")\n", + "# obj_val, sol = gurobi_maxcut(graph)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 4081189 3127413 -8.84763 56 71 -7.15263 -8.84763 23.7% 29.7 1190s\n" + ] + } + ], + "source": [ + "# Parse all files in data_dir \n", + "data_dir = f\"../data/vna/SK\"\n", + "input_files = [ f for f in os.listdir(data_dir) ]\n", + "input_files = natsorted(input_files)\n", + "\n", + "# Write directories & files\n", + "results_dir = f\"../results/{'/'.join(data_dir.split('/')[2:])}\"\n", + "solutions_dir = f\"../solutions/{'/'.join(data_dir.split('/')[2:])}\"\n", + "os.makedirs(results_dir, exist_ok=True)\n", + "os.makedirs(solutions_dir, exist_ok=True)\n", + "\n", + "i = 0\n", + "for file in input_files:\n", + " i += 1\n", + " if i > 25: break\n", + "\n", + " graph = read_nxgraph(f'{data_dir}/{file}') \n", + " obj_val, sol = gurobi_maxcut(graph)\n", + "\n", + " with open(f\"{results_dir}/GUROBI.txt\", \"a\") as f:\n", + " f.write(f\"{obj_val}\\n\")\n", + "\n", + " with open(f\"{solutions_dir}/GUROBI.txt\", \"a\") as f:\n", + " f.write(f\"{sol}\\n\")\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/src/baseline/README.md b/src/baseline/README.md index bee3a2b..274bb45 100644 --- a/src/baseline/README.md +++ b/src/baseline/README.md @@ -1,11 +1,95 @@ # Baseline MIP Solvers - - -| Solvers | +| Solvers | |-----------| | Gurobi | | ILOG CPLEX| | COPT | +## gurobi.py + +### Single file + +python script.py path/to/graph1.txt -t 3600 + +### Multiple specific files + +python script.py path/to/graph1.txt path/to/graph2.txt -t 3600 + +### Directory Batch Mode (Recommended) + +Automatically finds and processes all .txt files recursively within a given directory, outputting the results to a CSV. + +python script.py -i src/data/graphs -o results/gurobi/newest --csv_out batch_results.csv + +--- + +### Command Line Arguments + +| Flag | Long Name | Default | Description | +| :--- | :--- | :--- | :--- | +| paths | N/A | [specific test file] | Unflagged arguments at the end of the command are treated as manual file paths. | +| -t | --time_limit | 3600 | Time limit for the Gurobi solver per graph, in seconds. | +| -i | --in_dir | None | Input directory. Processes all .txt files found inside. | +| -o | --out_dir | results/gurobi/newest | Output directory. All .log files, individual .txt solution files, and the CSV are saved here. Creates the folder if it does not exist. | +| N/A | --csv_out | gurobi_batch_results.csv | The name of the aggregate CSV file. It is saved directly inside the --out_dir. | +| N/A | --retry_license | False | Retry entries that previously failed due to license limits. | +| N/A | --retry_timeouts | False | Retry entries that previously hit the time limit. | + +### Output Structure + +For every successfully solved graph, the script generates three pieces of data in the --out_dir: + +1. [graph_name].log: The raw solver output straight from Gurobi, containing simplex iterations, heuristic steps, and final bounds. +2. [graph_name].txt: A summarized output file containing the Objective Value, Objective Bound, Duration, MIP Gap, and the raw binary sequence of the best solution. +3. CSV Entry: A single row in the aggregate .csv file containing the Absolute Path to the original graph and its Objective Value (or state flag like TIMEOUT / LICENSE_LIMIT). + +--- + +### State Management + +When running in Directory Batch Mode (-i), the script executes a Synchronization Phase before solving any graphs. It cross-references the existing output CSV with the .log files in the output directory to map the state of the data. + +### Smart Resumption & Interruption Handling + +If a batch is canceled midway using Ctrl+C, or if the system crashes: + +* Clean Exits: Pressing Ctrl+C immediately halts the batch. The script will not attempt to write incomplete data to the CSV. +* Orphan Cleanup: On the next run, the script detects .log files containing "Solve interrupted". It automatically deletes these junk logs and incomplete .txt files, forcing a clean rerun of that specific graph. +* Resume Capability: Fully completed graphs recorded in the CSV (or validated by a successful .log file) are permanently skipped, allowing large batches to resume exactly where they left off. + +### Timeout Management + +If Gurobi reaches a defined time limit (e.g., 3600 seconds) before finding the optimal solution: + +* Gurobi gracefully stops and saves the best objective value found up to that point. +* The script records this partial solution to the .txt and .log files. +* In the CSV, the value is appended with a (TIMEOUT) flag (e.g., -70.5 (TIMEOUT)) to allow for easy filtering of sub-optimal runs during data analysis. + +### Directory Timeout Skipping + +When processing large nested datasets, graphs within the same subfolder often share similar complexity and execution times. To optimize batch processing, the script tracks timeouts at the directory level: + +* Subfolder Monitoring: If any single graph triggers a time limit (TIMEOUT) during execution, the script flags the parent subfolder. +* Automatic Bypassing: Once a subfolder is flagged, the script instantly bypasses all remaining unprocessed .txt files within that specific subfolder and moves forward to the next directory. +* Prevention of Bottlenecks: This mechanism prevents the solver from stalling on a directory filled with massive graphs, allowing the batch to continue clearing out simpler subfolders elsewhere in the dataset. + +### License Limit Protection + +If a graph is too large for the current Gurobi license, the solver will throw a size-limit exception. + +* Instead of crashing the entire batch or entering an infinite retry loop, the script catches this specific error. +* It writes LICENSE_LIMIT to the CSV for that specific graph. +* On future runs, the script recognizes LICENSE_LIMIT as a "completed" state and skips the graph permanently, saving processing time. + +### Duplicate Conflict Handling + +If the script detects that a file has been processed twice with conflicting results (e.g., folders were moved and the script was run twice), it triggers an interactive terminal prompt: + +[!] DUPLICATE CONFLICT: 100_SK_seed33.txt + 1: Keep First Value -> -70.5 + 2: Keep Second Value -> -71.2 + 3: Delete Both (forces rerun) + +If the duplicate values are identical, the script silently resolves the conflict and keeps one entry, requiring no manual input. diff --git a/src/baseline/format_gurobi.py b/src/baseline/format_gurobi.py new file mode 100644 index 0000000..7683b46 --- /dev/null +++ b/src/baseline/format_gurobi.py @@ -0,0 +1,150 @@ +import csv +import re +import os +import argparse +from collections import defaultdict + + +ACRONYMS = { + "EA": "Edwards-Anderson", + "SK": "Sherrington-Kirkpatrick", + # "FM": "Ferromagnetic", + # "AFM": "Antiferromagnetic", + # "RRG": "Random Regular Graph", + # "VNA": "Vertex-Normalized Adjacency", + # "BA": "Barabási-Albert", + # "ER": "Erdős-Rényi" +} + + +def natural_sort_key(text): + """ + Splits a string into string and integer components so that lists sort naturally. + Example: 10x10 comes before 20x20. 2 comes before 10. + """ + return [int(c) if c.isdigit() else c.lower() for c in re.split(r"(\d+)", str(text))] + + +def format_header_name(text): + """ + Cleans up folder names by replacing underscores with spaces and + expanding known acronyms. + """ + text = text.replace("_", " ") + words = text.split() + expanded_words = [ACRONYMS.get(w, w) for w in words] + return " ".join(expanded_words) + + +def format_csv(input_csv, output_csv): + if not os.path.exists(input_csv): + print(f"[!] Input file '{input_csv}' not found.") + return + + # Data Structure: grouped_data[hierarchy_tuple][size_label] = [(seed, energy)] + grouped_data = defaultdict(lambda: defaultdict(list)) + + with open(input_csv, mode="r", newline="") as infile: + reader = csv.reader(infile) + header = next(reader, None) + + for row in reader: + if len(row) < 2: + continue + + raw_path = row[0].replace("\\", "/") + obj_val = row[1] + + clean_path = re.sub(r"^.*?/data/[^/]+/", "", raw_path) + + parts = clean_path.split("/") + if len(parts) < 2: + continue + + filename = parts[-1] + original_size_folder = parts[-2] + hierarchy_folders = parts[:-2] + + size_match = re.search(r"(\d+x\d+)", filename) + if size_match: + size_label = size_match.group(1) + else: + size_label = original_size_folder + + # Always grab the last number in the filename + numbers = re.findall(r"\d+", filename) + seed = int(numbers[-1]) if numbers else 0 + + hierarchy_tuple = tuple(hierarchy_folders) + grouped_data[hierarchy_tuple][size_label].append((seed, obj_val)) + + with open(output_csv, mode="w", newline="") as outfile: + writer = csv.writer(outfile) + + sorted_hierarchies = sorted( + grouped_data.keys(), key=lambda x: [natural_sort_key(part) for part in x] + ) + + for hierarchy in sorted_hierarchies: + sizes_dict = grouped_data[hierarchy] + + for folder_name in hierarchy: + writer.writerow([format_header_name(folder_name)]) + + sorted_sizes = sorted(sizes_dict.keys(), key=natural_sort_key) + + size_row = [] + for size in sorted_sizes: + size_row.extend([size, "", ""]) + writer.writerow(size_row) + + header_row = [] + for _ in sorted_sizes: + header_row.extend(["Seed", "Min Energy", ""]) + writer.writerow(header_row) + + for size in sorted_sizes: + sizes_dict[size].sort(key=lambda x: x[0]) + + max_rows = max(len(sizes_dict[size]) for size in sorted_sizes) + + for i in range(max_rows): + data_row = [] + for size in sorted_sizes: + lst = sizes_dict[size] + if i < len(lst): + seed_val, energy_val = lst[i] + data_row.extend([seed_val, energy_val, ""]) + else: + data_row.extend(["", "", ""]) + + writer.writerow(data_row) + + writer.writerow([]) + + print(f"Spreadsheet formatting complete. Saved to '{output_csv}'.") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Format Gurobi batch CSV into a readable spreadsheet grid." + ) + + parser.add_argument( + "input_csv", + type=str, + default="gurobi_batch_results.csv", + nargs="?", + help="The raw CSV output from the Gurobi batch solver.", + ) + + parser.add_argument( + "-o", + "--output_csv", + type=str, + default="formatted_results.csv", + help="The name of the new formatted output file.", + ) + + args = parser.parse_args() + format_csv(args.input_csv, args.output_csv) diff --git a/src/baseline/gurobi.lic b/src/baseline/gurobi.lic new file mode 100644 index 0000000..f442554 --- /dev/null +++ b/src/baseline/gurobi.lic @@ -0,0 +1,20 @@ +# DO NOT EDIT THIS FILE except as noted +# +# +# Place this file in your home directory that takes precedence +# or in one of the following shared locations: +# * C:\gurobi\ on Windows +# * /opt/gurobi/ on Linux +# * /Library/gurobi/ on Mac OS X +# Or set the environment variable GRB_LICENSE_FILE to point to this file, +# it will override the default locations +# +LICENSEID=2777688 +TYPE=ACADEMIC +VERSION=13 +HOSTNAME=balltoucher +HOSTID=5942879a +USERNAME=garrett +EXPIRATION=2027-02-10 +KEY=MKEFLWWA +CKEY=MKEFLWWA diff --git a/src/baseline/gurobi.py b/src/baseline/gurobi.py index 41cd04c..a9a3269 100644 --- a/src/baseline/gurobi.py +++ b/src/baseline/gurobi.py @@ -1,25 +1,27 @@ import networkx as nx -import numpy as np import os -import math import gurobipy as gp +import argparse +import csv +import re from utils import read_nxgraph from utils import float_to_binary from utils import base64_encode -def save_to_file(model, file_path, time_limit, print_terminal=True): - - file_dir = f"results/gurobi/{"/".join(file_path.split("/")[:-1])}" - file_name = file_path.split("/")[-1] +def save_to_file(model, file_name, time_limit, out_dir, print_terminal=True): obj_val = model.ObjVal obj_bnd = model.ObjBound duration = model.Runtime mip_gap = model.MIPGap - best_solution = "".join([ float_to_binary(x.X) for x in model.getVars() ]) + best_solution = "".join([float_to_binary(x.X) for x in model.getVars()]) best_encoded = base64_encode(best_solution) - with open(f"{file_dir}/{file_name}", "w") as f: + if not os.path.exists(out_dir): + os.makedirs(out_dir) + + base_name = os.path.basename(file_name) + with open(os.path.join(out_dir, base_name), "w") as f: f.write(f"Objective Value: {obj_val}\n") f.write(f"Objective Bound: {obj_bnd}\n") f.write(f"Duration: {duration}\n") @@ -37,9 +39,10 @@ def save_to_file(model, file_path, time_limit, print_terminal=True): print(f"Best Solution Encoded: {best_encoded}") print(f"Best Solution Raw: {best_solution}\n") + return obj_val -def gurobi_solve(graph, file_path, time_limit=3600): +def gurobi_solve(graph, file_name, time_limit, out_dir): nodes = len(list(graph.nodes)) J = nx.to_numpy_array(graph) @@ -48,35 +51,370 @@ def gurobi_solve(graph, file_path, time_limit=3600): x = model.addVars(nodes, vtype=gp.GRB.BINARY) objective = gp.quicksum( -J[i, j] * (2 * x[i] - 1) * (2 * x[j] - 1) - for i in range(nodes) for j in range(i+1, nodes) if J[i, j] != 0.0 + for i in range(nodes) + for j in range(i + 1, nodes) + if J[i, j] != 0.0 ) model.setObjective(objective, gp.GRB.MINIMIZE) - model.setParam('TimeLimit', time_limit) + model.setParam("TimeLimit", time_limit) model.setParam("MIPGap", 0.0) - file_dir = f"results/gurobi/{"/".join(file_path.split("/")[:-1])}" - file_name = file_path.split("/")[-1] - if not os.path.exists(file_dir): - os.makedirs(file_dir) - model.Params.LogFile = f"{file_dir}/{file_name.split(".")[0]}.log" + if not os.path.exists(out_dir): + os.makedirs(out_dir) + + base_name = os.path.basename(file_name) + model.Params.LogFile = os.path.join(out_dir, f"{base_name.split('.')[0]}.log") model.optimize() - save_to_file(model, file_path, time_limit) -def single_shot_instance(file_name, time_limit): + if model.Status == gp.GRB.INTERRUPTED: + return "INTERRUPTED" + + obj_val = save_to_file(model, file_name, time_limit, out_dir) + + if model.Status == gp.GRB.TIME_LIMIT: + return f"{obj_val} (TIMEOUT)" + + return obj_val + + +def single_shot_instance(file_name, time_limit, out_dir): graph = read_nxgraph(file_name) - file_path = file_name[5:] - gurobi_solve(graph, file_path, time_limit) + return gurobi_solve(graph, file_name, time_limit, out_dir) + -def multiple_shot_instance(file_names, time_limit): +def multiple_shot_instance(file_names, time_limit, out_dir): for file_name in file_names: - single_shot_instance(file_name, time_limit) + single_shot_instance(file_name, time_limit, out_dir) + + +def directory_shot_instance( + in_dir, time_limit, out_dir, output_csv, retry_license=False, retry_timeouts=False +): + if not os.path.isdir(in_dir): + print(f"\n[!] ERROR: The input directory '{in_dir}' does not exist.") + return + + if not os.path.exists(out_dir): + os.makedirs(out_dir) + + csv_path = os.path.join(out_dir, output_csv) + csv_data = {} + synced_something = False + + if os.path.exists(csv_path): + with open(csv_path, mode="r", newline="") as f: + reader = csv.reader(f) + next(reader, None) + for row in reader: + if len(row) >= 2: + abs_path = row[0] + new_val = str(row[1]).strip() + + if abs_path in csv_data: + old_val = csv_data[abs_path] + if old_val == new_val: + synced_something = True + else: + print( + f"\n[!] DUPLICATE CONFLICT: {os.path.basename(abs_path)}" + ) + print(f" 1: Keep First Value -> {old_val}") + print(f" 2: Keep Second Value -> {new_val}") + print(f" 3: Delete Both (forces rerun)") + while True: + choice = input("Choose an option (1/2/3): ").strip() + if choice == "1": + synced_something = True + break + elif choice == "2": + csv_data[abs_path] = new_val + synced_something = True + break + elif choice == "3": + del csv_data[abs_path] + synced_something = True + break + else: + print("Invalid input. Please enter 1, 2, or 3.") + else: + csv_data[abs_path] = new_val + + for root, dirs, files in os.walk(in_dir): + rel_path = os.path.relpath(root, in_dir) + current_out_dir = os.path.join(out_dir, rel_path) + + for file in files: + if not file.endswith(".txt") or file == output_csv: + continue + + abs_path = os.path.abspath(os.path.join(root, file)) + log_file_name = f"{file.split('.')[0]}.log" + log_path = os.path.join(current_out_dir, log_file_name) + + if abs_path in csv_data: + val = csv_data[abs_path] + is_binary = ( + all(c in "01" for c in val.replace(" (TIMEOUT)", "")) + and len(val.replace(" (TIMEOUT)", "")) >= 8 + ) + + if "LICENSE_LIMIT" in val or "size-limited" in val.lower(): + if retry_license: + del csv_data[abs_path] + synced_something = True + else: + csv_data[abs_path] = "LICENSE_LIMIT" + + elif "(TIMEOUT)" in val and retry_timeouts: + del csv_data[abs_path] + synced_something = True + + elif "INTERRUPTED" in val or "ERROR" in val: + del csv_data[abs_path] + synced_something = True + + elif is_binary: + if os.path.exists(log_path): + with open(log_path, "r", errors="ignore") as log_f: + log_content = log_f.read() + if ( + "size-limited" in log_content.lower() + or "too large" in log_content.lower() + ): + if retry_license: + del csv_data[abs_path] + synced_something = True + else: + csv_data[abs_path] = "LICENSE_LIMIT" + synced_something = True + elif "Solve interrupted" in log_content: + del csv_data[abs_path] + synced_something = True + elif "Time limit reached" in log_content and retry_timeouts: + del csv_data[abs_path] + synced_something = True + else: + match = re.search( + r"Objective Value:\s*([-\d\.]+)", log_content + ) + if match: + obj_val = match.group(1) + if "Time limit reached" in log_content: + obj_val += " (TIMEOUT)" + csv_data[abs_path] = obj_val + synced_something = True + else: + del csv_data[abs_path] + synced_something = True + else: + del csv_data[abs_path] + synced_something = True + else: + if os.path.exists(log_path): + with open(log_path, "r", errors="ignore") as log_f: + log_content = log_f.read() + if "Solve interrupted" in log_content: + del csv_data[abs_path] + synced_something = True + elif "Time limit reached" in log_content: + if retry_timeouts: + del csv_data[abs_path] + synced_something = True + elif "(TIMEOUT)" not in val: + csv_data[abs_path] = f"{val} (TIMEOUT)" + synced_something = True + + else: + if os.path.exists(log_path): + with open(log_path, "r", errors="ignore") as log_f: + log_content = log_f.read() + + if ( + "size-limited" in log_content.lower() + or "too large" in log_content.lower() + ): + if retry_license: + os.remove(log_path) + txt_out_path = os.path.join(current_out_dir, file) + if os.path.exists(txt_out_path): + os.remove(txt_out_path) + else: + csv_data[abs_path] = "LICENSE_LIMIT" + synced_something = True + elif "Solve interrupted" in log_content: + os.remove(log_path) + txt_out_path = os.path.join(current_out_dir, file) + if os.path.exists(txt_out_path): + os.remove(txt_out_path) + elif "Time limit reached" in log_content and retry_timeouts: + os.remove(log_path) + txt_out_path = os.path.join(current_out_dir, file) + if os.path.exists(txt_out_path): + os.remove(txt_out_path) + else: + match = re.search( + r"Objective Value:\s*([-\d\.]+)", log_content + ) + if match: + obj_val = match.group(1) + if "Time limit reached" in log_content: + obj_val += " (TIMEOUT)" + csv_data[abs_path] = obj_val + synced_something = True + + if synced_something or not os.path.exists(csv_path): + with open(csv_path, mode="w", newline="") as f: + writer = csv.writer(f) + writer.writerow(["Absolute Path", "Objective Value"]) + for p, v in csv_data.items(): + writer.writerow([p, v]) + + timed_out_folders = set() + for p, v in csv_data.items(): + if "(TIMEOUT)" in v: + timed_out_folders.add(os.path.dirname(p)) + + with open(csv_path, mode="a", newline="") as f: + writer = csv.writer(f) + + for root, dirs, files in os.walk(in_dir): + abs_root = os.path.abspath(root) + rel_path = os.path.relpath(root, in_dir) + current_out_dir = os.path.join(out_dir, rel_path) + + if not os.path.exists(current_out_dir): + os.makedirs(current_out_dir) + + if abs_root in timed_out_folders: + unprocessed_count = sum( + 1 + for f_name in files + if f_name.endswith(".txt") + and f_name != output_csv + and os.path.abspath(os.path.join(root, f_name)) not in csv_data + ) + if unprocessed_count > 0: + print( + f"\n[*] Skipping {unprocessed_count} unprocessed file(s) in {os.path.basename(root)} due to previous timeout." + ) + continue + + for file in files: + if not file.endswith(".txt") or file == output_csv: + continue + + abs_path = os.path.abspath(os.path.join(root, file)) + + if abs_path in csv_data: + continue + + print(f"\nProcessing: {file}") + try: + obj_val = single_shot_instance( + file_name=os.path.join(root, file), + time_limit=time_limit, + out_dir=current_out_dir, + ) + + if obj_val == "INTERRUPTED": + print("\n[!] Run stopped by user. Halting batch safely.") + return + + writer.writerow([abs_path, obj_val]) + csv_data[abs_path] = obj_val + + if "(TIMEOUT)" in str(obj_val): + timed_out_folders.add(abs_root) + print( + f"[*] Timeout reached! Skipping any remaining files in '{os.path.basename(root)}'." + ) + break + + except KeyboardInterrupt: + print("\n[!] Run stopped by user. Halting batch safely.") + return + except Exception as e: + error_str = str(e).lower() + if "size-limited" in error_str or "too large" in error_str: + writer.writerow([abs_path, "LICENSE_LIMIT"]) + csv_data[abs_path] = "LICENSE_LIMIT" + else: + writer.writerow([abs_path, f"ERROR: {e}"]) + + f.flush() + + print(f"\nBatch processing complete. Results saved to '{csv_path}'.") + if __name__ == "__main__": - file_name = "../../data/DRL/2D/10/DRL_10_ID0.txt" - time_limit = 3600 - single_shot_instance(file_name, time_limit) + parser = argparse.ArgumentParser( + description="Run Gurobi MaxCut QUBO solver on graphs." + ) + + parser.add_argument( + "paths", + nargs="*", + default=["../../src/data/1D/VNA/Chain/32/ising_chain_32_seed1.txt"], + help="One or more file paths to process. Defaults to a specific test file if omitted.", + ) + + parser.add_argument( + "-t", + "--time_limit", + type=int, + default=3600, + help="Time limit for the solver in seconds (default: 3600).", + ) + + parser.add_argument( + "-i", + "--in_dir", + type=str, + help="Input directory. Process all .txt files in this folder.", + ) + + parser.add_argument( + "-o", + "--out_dir", + type=str, + default="results/gurobi/newest", + help="Output directory for logs, individual text files, and the CSV.", + ) + + parser.add_argument( + "--csv_out", + type=str, + default="gurobi_batch_results.csv", + help="Name of the output CSV file (saved inside out_dir).", + ) + + parser.add_argument( + "--retry_license", + action="store_true", + help="Retry entries that previously failed due to license limits.", + ) + + parser.add_argument( + "--retry_timeouts", + action="store_true", + help="Retry entries that previously hit the time limit.", + ) + + args = parser.parse_args() - # file_names = [ "data/vna/ER/40x40_uniform_seed1.txt", "data/vna/ER/40x40_uniform_seed15.txt", "data/vna/ER/40x40_uniform_seed25.txt" ] - # time_limit = 3600 - # multiple_shot_instance(file_names, time_limit) \ No newline at end of file + if args.in_dir: + directory_shot_instance( + args.in_dir, + args.time_limit, + args.out_dir, + args.csv_out, + args.retry_license, + args.retry_timeouts, + ) + elif len(args.paths) == 1: + single_shot_instance(args.paths[0], args.time_limit, args.out_dir) + elif len(args.paths) > 1: + multiple_shot_instance(args.paths, args.time_limit, args.out_dir) + else: + print("No paths or directories were provided.") diff --git a/src/tutorials/VCA/.ipynb_checkpoints/VNA_Tutorial-checkpoint.ipynb b/src/tutorials/VCA/.ipynb_checkpoints/VNA_Tutorial-checkpoint.ipynb new file mode 100644 index 0000000..a6fb725 --- /dev/null +++ b/src/tutorials/VCA/.ipynb_checkpoints/VNA_Tutorial-checkpoint.ipynb @@ -0,0 +1,608 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "99676ca3", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/garrett/.conda/envs/vca/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:526: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n", + " _np_qint8 = np.dtype([(\"qint8\", np.int8, 1)])\n", + "/home/garrett/.conda/envs/vca/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:527: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n", + " _np_quint8 = np.dtype([(\"quint8\", np.uint8, 1)])\n", + "/home/garrett/.conda/envs/vca/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:528: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n", + " _np_qint16 = np.dtype([(\"qint16\", np.int16, 1)])\n", + "/home/garrett/.conda/envs/vca/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:529: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n", + " _np_quint16 = np.dtype([(\"quint16\", np.uint16, 1)])\n", + "/home/garrett/.conda/envs/vca/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:530: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n", + " _np_qint32 = np.dtype([(\"qint32\", np.int32, 1)])\n", + "/home/garrett/.conda/envs/vca/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:535: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n", + " np_resource = np.dtype([(\"resource\", np.ubyte, 1)])\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Successfully imported local modules.\n" + ] + } + ], + "source": [ + "import tensorflow as tf\n", + "import numpy as np\n", + "import random\n", + "import time\n", + "import argparse\n", + "from math import sqrt\n", + "\n", + "try:\n", + " from config import config\n", + " from DilatedRNN import DilatedRNNWavefunction\n", + " from utils import Fullyconnected_localenergies, Fullyconnected_diagonal_matrixelements\n", + " from vca import vca_solver\n", + " print(\"Successfully imported local modules.\")\n", + "except ImportError:\n", + " print(\"Local .py files not found. Please ensure they are in the same directory or paste the classes below.\")\n", + "\n", + "tf.compat.v1.disable_eager_execution()" + ] + }, + { + "cell_type": "markdown", + "id": "b676f9fb", + "metadata": {}, + "source": [ + "# Section 1: Ising Model as Energy Minimization" + ] + }, + { + "cell_type": "markdown", + "id": "29c29caa", + "metadata": {}, + "source": [ + "## 1.1 What is the Ising Model?\n", + "The Ising model was originally developed in statistical physics to describe the **ferromagnetism** of solid materials. It assumes that a system consists of a collection of **spins** arranged on a lattice, where each spin can only take one of two discrete values:\n", + "* $s_i = +1$ (Spin Up)\n", + "* $s_i = -1$ (Spin Down)\n", + "\n", + "In this implementation, the spins are mapped to binary values $s \\in \\{0, 1\\}$ to be compatible with neural network processing." + ] + }, + { + "cell_type": "markdown", + "id": "e6c98273", + "metadata": {}, + "source": [ + "## 1.2 Mathematical Nature & The Hamiltonian\n", + "The **Hamiltonian** is a function in physics that represents the total energy of a system. For a fully connected Ising model, the energy function $E(\\mathbf{s})$ is defined as:\n", + "\n", + "$$E(\\mathbf{s}) = -\\sum_{i < j} J_{ij} \\sigma_i \\sigma_j - \\sum_i h_i \\sigma_i$$\n", + "\n", + "* **$J_{ij}$ (Interaction Term)**: Defines the coupling strength between spins $i$ and $j$.\n", + " * If $J_{ij} > 0$: Spins tend to align in the same direction (Ferromagnetic).\n", + " * If $J_{ij} < 0$: Spins tend to align in opposite directions (Anti-ferromagnetic).\n", + "* **$h_i$ (External Field)**: Represents the influence of an external magnetic field or bias on an individual spin. In this code, a transverse field `Bx` is introduced to add non-diagonal perturbations.\n", + "* **Code Example**: In `utils.py`, the `Jz` matrix stores all interaction weights, and the function `Fullyconnected_diagonal_matrixelements` calculates the classical energy based on these weights.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "b4a03bde", + "metadata": {}, + "outputs": [], + "source": [ + "def Fullyconnected_diagonal_matrixelements(Jz, samples):\n", + " numsamples = samples.shape[0]\n", + " N = samples.shape[1]\n", + " energies = np.zeros((numsamples), dtype = np.float64)\n", + "\n", + " for i in range(N-1):\n", + " values = np.expand_dims(samples[:,i], axis = -1)+samples[:,i+1:]\n", + " valuesT = np.copy(values)\n", + " valuesT[values==2] = +1 #If both spins are up\n", + " valuesT[values==0] = +1 #If both spins are down\n", + " valuesT[values==1] = -1 #If they are opposite\n", + "\n", + " energies += np.sum(valuesT*(-Jz[i,i+1:]), axis = 1)\n", + "\n", + " return energies" + ] + }, + { + "cell_type": "markdown", + "id": "50879217", + "metadata": {}, + "source": [ + "## 1.3 Problem Type: Combinatorial Optimization\n", + "From a computer science perspective, the Ising model is a classic **Combinatorial Optimization Problem**:\n", + "* **Goal**: To find a specific arrangement among $2^N$ possible discrete configurations that minimizes the total energy of the system.\n", + "* **Applications**: Many famous NP-hard problems, such as the **Max-Cut Problem**, **Traveling Salesperson Problem (TSP)**, or **Graph Coloring**, can be exactly mapped onto the energy minimization of an Ising model." + ] + }, + { + "cell_type": "markdown", + "id": "ad844f1b", + "metadata": {}, + "source": [ + "## 1.4 The Essence of NP-hardness\n", + "Finding the lowest energy state (the **Ground State**) of an Ising model is extremely difficult for the following reasons:\n", + "\n", + "1. **Exponential State Space**: For $N$ spins, there are $2^N$ possible states. Brute-force search becomes impossible as $N$ grows.\n", + "2. **Frustration**: When the distribution of $J_{ij}$ is complex (as in Spin Glass models), the system cannot satisfy all coupling terms simultaneously leading towards numerous equivalent ground states.\n", + "3. **Complex Energy Landscape**: Frustration leads to a highly non-convex energy landscape filled with numerous **Local Minima**.\n", + "4. **Mathematical Bottleneck**: Due to these high energy barriers, traditional optimization algorithms easily get stuck in local optima, failing to find the global optimum. This is the core of its NP-hardness." + ] + }, + { + "cell_type": "markdown", + "id": "304e4c2d", + "metadata": {}, + "source": [ + "# Section 2:Simulated Annealing as a Classical Heuristic" + ] + }, + { + "cell_type": "markdown", + "id": "0435161b", + "metadata": {}, + "source": [ + "## 2.1 What is Simulated Annealing\n", + "Simulated Annealing (SA) is a classical **probabilistic optimization algorithm** designed to find low-energy (or high-quality) solutions in large, discrete, and non-convex search spaces. It was originally inspired by the physical annealing process in metallurgy, but in optimization it should be understood primarily as a **randomized search heuristic** rather than a physical simulation.\n", + "\n", + "Formally, simulated annealing defines a **discrete-time, inhomogeneous Markov chain** over the solution space. Starting from an initial solution, the algorithm repeatedly proposes local moves within a predefined neighborhood structure. If a proposed move improves the objective value, it is always accepted. Otherwise, the move may still be accepted with a probability that depends on both the **increase in cost** and a **temperature parameter** that decreases over time.\n", + "\n", + "The key feature of simulated annealing is its ability to **accept worse solutions during the early stages of the search**, which allows the algorithm to escape local minima that commonly trap greedy or purely descent-based methods. As the temperature is gradually lowered according to a cooling schedule, the algorithm becomes increasingly conservative and focuses on refining high-quality solutions.\n", + "\n", + "From a theoretical perspective, the acceptance rule of simulated annealing is closely related to the **Boltzmann (Gibbs) distribution**. When the temperature is held fixed, the induced Markov chain admits a Boltzmann distribution as its invariant distribution. As the temperature approaches zero, this distribution concentrates on the set of global minima of the objective function. A rigorous analysis of this viewpoint and the convergence properties of simulated annealing can be found in the classical work of Bertsimas and Tsitsiklis (1993)." + ] + }, + { + "cell_type": "markdown", + "id": "1e86dff1", + "metadata": {}, + "source": [ + "## 2.2 How Simulated Annealing Works\n", + "\n", + "This section reconstructs the simulated annealing procedure by following the basic elements defined in the classical formulation of Bertsimas and Tsitsiklis (1993).\n", + "\n", + "### Solution Space and Objective\n", + "\n", + "Simulated annealing operates on a **finite set of feasible solutions**, denoted by $S$. \n", + "Each solution represents a valid configuration of the problem (e.g., a spin configuration in the Ising model).\n", + "\n", + "An objective (or energy) function $J$ is defined on $S$, assigning a real-valued cost to each solution. \n", + "The goal of the algorithm is to find solutions with **minimal objective value**, ideally belonging to the set of global minima.\n", + "\n", + "### Neighborhood Structure\n", + "\n", + "For each solution $i \\in S$, a **neighborhood set** $S(i)$ is defined. \n", + "This set specifies which solutions can be reached from $i$ through a single local modification.\n", + "\n", + "The neighborhood structure determines how the algorithm explores the solution space and encodes the notion of a *local move*.\n", + "\n", + "### Randomized Move Proposal\n", + "\n", + "At each iteration, the algorithm selects a candidate solution from the neighborhood of the current solution. \n", + "This selection is randomized according to a fixed probability rule, ensuring that the search does not follow a deterministic path.\n", + "\n", + "This randomness is essential for exploration and allows the algorithm to visit different regions of the solution space.\n", + "\n", + "### Temperature and Acceptance Behavior\n", + "\n", + "The algorithm maintains a **temperature parameter** that changes over time according to a predefined cooling schedule. \n", + "The temperature controls how willing the algorithm is to accept worse solutions.\n", + "\n", + "- At high temperature, the algorithm is more exploratory and may accept worse solutions frequently.\n", + "- At low temperature, the algorithm becomes more conservative and favors improvements.\n", + "\n", + "This mechanism enables simulated annealing to escape local minima early in the search and gradually refine solutions as the temperature decreases.\n", + "\n", + "### State Evolution\n", + "\n", + "Starting from an initial solution, the algorithm repeatedly:\n", + "1. Proposes a neighboring solution at random,\n", + "2. Evaluates its objective value,\n", + "3. Decides whether to accept or reject the move based on the current temperature.\n", + "\n", + "Through this process, the algorithm generates a sequence of solutions that gradually concentrates on low-energy regions of the solution space.\n", + "\n", + "---\n", + "\n", + "**Remark.** \n", + "From a theoretical perspective, this iterative process defines a time-inhomogeneous Markov chain over the solution space. When the temperature is held constant, the long-run behavior of the algorithm favors low-cost solutions, providing the foundation for simulated annealing as a global optimization heuristic." + ] + }, + { + "cell_type": "markdown", + "id": "1339f983", + "metadata": {}, + "source": [ + "## 2.2.1 Boltzmann Distribution (Intuition)\n", + "\n", + "This section provides additional intuition for the **acceptance behavior** described in the previous section. \n", + "It explains *why* simulated annealing uses a temperature-dependent probabilistic acceptance rule.\n", + "\n", + "### Preference for Low-Energy Solutions\n", + "\n", + "The Boltzmann distribution expresses a simple and intuitive idea:\n", + "\n", + "**solutions with lower energy (or lower cost) are more likely than solutions with higher energy**.\n", + "\n", + "This preference is controlled by a temperature parameter:\n", + "- At **high temperature**, differences in energy matter less, and many solutions are treated similarly.\n", + "- At **low temperature**, low-energy solutions are strongly preferred over higher-energy ones.\n", + "\n", + "This temperature-dependent preference mirrors how physical systems behave under thermal fluctuations.\n", + "\n", + "### Boltzmann Form (Interpretation)\n", + "\n", + "When the temperature is held fixed, the search process favors solutions according to a Boltzmann-like preference:\n", + "\n", + "$$\n", + "\\pi(i) \\propto \\exp\\!\\left(-\\frac{J(i)}{T}\\right),\n", + "$$\n", + "\n", + "where $J(i)$ denotes the objective (or energy) of solution $i$, and $T$ is the temperature.\n", + "\n", + "This distribution is **not computed explicitly** during the algorithm. \n", + "Rather, the local acceptance rule of simulated annealing is designed so that the overall search behavior is consistent with this preference.\n", + "\n", + "### Key Takeaway\n", + "\n", + "**The Boltzmann distribution explains why simulated annealing behaves the way it does.** \n", + "It provides a theoretical lens for understanding how temperature controls exploration and exploitation, while the algorithm itself remains a simple, local, randomized search procedure.\n" + ] + }, + { + "cell_type": "markdown", + "id": "89c5f294", + "metadata": {}, + "source": [ + "# Section 3: From Simulated Annealing to Learing Methods" + ] + }, + { + "cell_type": "markdown", + "id": "bf3454a3", + "metadata": {}, + "source": [ + "## 3.1 Limitations of Simulated Annealing\n", + "Simulated annealing is a well-established classical heuristic, but it exhibits several limitations when applied to large-scale and highly structured optimization problems.\n", + "\n", + "The search behavior of simulated annealing is **manually designed**.\n", + "Key components such as neighborhood moves, acceptance rules, and temperature schedules are specified a priori and remain fixed, limiting the algorithm’s ability to adapt to instance-specific structure.\n", + "\n", + "Moreover, simulated annealing does not **accumulate or reuse experience**.\n", + "Each run is independent, and information obtained during previous searches is not leveraged to improve future performance.\n", + "\n", + "In rugged energy landscapes such as spin glass systems, the resulting random walk may spend a substantial amount of time exploring suboptimal regions before reaching low-energy configurations.\n", + "\n", + "These limitations motivate learning-based approaches, which aim not to eliminate randomness, but to **learn heuristics that more effectively guide the search process**." + ] + }, + { + "cell_type": "markdown", + "id": "e20e7ab7", + "metadata": {}, + "source": [ + "## 3.2 Why learning methods make sense\n", + "\n", + "Compared to classical simulated annealing, learning-based approaches do not aim to replace stochastic search mechanisms, but rather to augment and adapt them. Classical simulated annealing relies on hand-designed components—such as neighborhood moves and temperature schedules—that remain fixed across problem instances. Learning-based methods retain the stochastic nature of annealing, while introducing learnable components that can adapt these mechanisms based on data or accumulated search experience.\n", + "\n", + "In this sense, learning serves as a way to parameterize and guide the annealing process rather than eliminating it. While each simulated annealing run still performs randomized exploration in the solution space, learning-based approaches can amortize experience across multiple instances, enabling the search to be biased toward more promising regions of the energy landscape. This experience-driven guidance is particularly valuable in large-scale NP-hard problems, where purely random exploration may suffer from slow convergence and high computational cost.\n", + "\n", + "Importantly, learning-based methods do not circumvent the NP-hard nature of the problem or remove randomness from the search process. Instead, they improve practical efficiency by shaping the stochastic dynamics of simulated annealing, allowing randomized search to operate more effectively within complex energy landscapes." + ] + }, + { + "cell_type": "markdown", + "id": "382bf9a2", + "metadata": {}, + "source": [ + "# Section 4: VNA as Example" + ] + }, + { + "cell_type": "markdown", + "id": "841dd0cf", + "metadata": {}, + "source": [ + "Here is the link to VNA Github: https://github.com/VectorInstitute/VariationalNeuralAnnealing/tree/main?tab=readme-ov-file" + ] + }, + { + "cell_type": "markdown", + "id": "1cfaf27c", + "metadata": {}, + "source": [ + "## 4.1 Motivation\n", + "\n", + "Simulated annealing solves optimization problems by approximately sampling a sequence of Boltzmann distributions, but its effectiveness is often limited by slow sampling dynamics in rugged or glassy energy landscapes. In practical, finite-time runs, the system deviates from equilibrium, leading to suboptimal solutions.\n", + "\n", + "Variational Neural Annealing addresses this limitation by replacing the exact Boltzmann distribution with a flexible, parameterized model that can be efficiently sampled. Instead of annealing the physical system itself, VNA anneals a learnable probability distribution whose parameters are optimized using variational principles and energy feedback. Modern autoregressive neural networks provide an ideal parameterization, as they allow exact and fast sampling even when representing complex, highly structured landscapes. This shift from trajectory-based sampling to distribution-based learning enables more efficient exploration of low-energy configurations under limited computational budgets." + ] + }, + { + "cell_type": "markdown", + "id": "76e94747", + "metadata": {}, + "source": [ + "## 4.2 Algorithm Logic\n", + "\n", + "This section provides a step-by-step walkthrough of a learning-based annealing method for Ising optimization. \n", + "Each conceptual step is explicitly linked to a corresponding function in the implementation, allowing readers to connect algorithmic ideas with concrete code components.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "02fbd61b", + "metadata": {}, + "source": [ + "### 1. Load an Ising Instance\n", + "\n", + "**Relevant functions:** \n", + "- `config.read_graph(graph_path)` \n", + "- `config.__init__()`\n", + "\n", + "We begin by loading an Ising problem instance from file. The graph structure and coupling strengths are encoded in a dense interaction matrix , which defines the Ising Hamiltonian. This step specifies the optimization problem and is independent of the solver." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fac15adb", + "metadata": {}, + "outputs": [], + "source": [ + "# Initialize the configuration with a specific Ising instance\n", + "instance_path = \"dataset/EA/EA_10x10/10x10_uniform_seed1.txt\" \n", + "cfg = config(instance_path, seed=42)\n", + "\n", + "# Verify the system scale and problem dimensions\n", + "print(f\"System Size (N): {cfg.N}\") # Number of spins (sites)\n", + "print(f\"J-Matrix Shape: {cfg.Jz.shape}\") # Should be an (N x N) matrix\n", + "\n", + "# Physical Consistency Check: The interaction J_ij must equal J_ji (Symmetry)\n", + "# In physics, the coupling between spin i and spin j is mutual.\n", + "is_symmetric = (cfg.Jz == cfg.Jz.T).all()\n", + "print(f\"Is the J-matrix symmetric? {is_symmetric}\")\n", + "\n", + "# Inspect a sub-section of the interaction matrix\n", + "# This represents a local 'snapshot' of the couplings between the first 5 spins.\n", + "print(\"\\nLocal Sampling of J-Matrix (Top 5x5):\")\n", + "print(cfg.Jz[:5, :5])" + ] + }, + { + "cell_type": "markdown", + "id": "334658f9", + "metadata": {}, + "source": [ + "### 2. Define the Neural Sampler\n", + "\n", + "**Relevant class:** \n", + "- `DilatedRNNWavefunction`\n", + "\n", + "The neural sampler parameterizes a probability distribution over spin configurations. We use an autoregressive recurrent neural network to model the joint distribution of spins, enabling exact and efficient sampling without Markov chain dynamics. The network serves as a learnable heuristic that guides the randomized search process." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e84ee6cc", + "metadata": {}, + "outputs": [], + "source": [ + "# Each RNN layer uses the same number of hidden units.\n", + "# cfg.num_layers controls how many dilated layers are used.\n", + "# cfg.num_units controls the capacity of each RNN cell.\n", + "units = [cfg.num_units] * cfg.num_layers\n", + "# Instantiate the autoregressive neural sampler\n", + "rnn_sampler = DilatedRNNWavefunction(\n", + " systemsize=cfg.N, # Number of spins in the Ising model\n", + " units=units, # Hidden units per RNN layer\n", + " layers=cfg.num_layers, # Number of dilated RNN layers\n", + " activation=cfg.activation_function,\n", + " seed=cfg.seed\n", + ")\n", + "\n", + "print(\"Number of spins (N):\", cfg.N)\n", + "print(\"Number of RNN layers:\", cfg.num_layers)\n", + "print(\"Hidden units per layer:\", cfg.num_units)\n", + "print(\"RNN cells per layer:\", len(rnn_sampler.rnn[0]))\n", + "print(\"Total number of RNN cells:\", cfg.num_layers * cfg.N)\n", + "#Each layer contains one RNN cell per spin position. \n", + "#This design allows the model to capture both short-range and long-range dependencies between spins through dilated connections." + ] + }, + { + "cell_type": "markdown", + "id": "4e8ac4fd", + "metadata": {}, + "source": [ + "### 3. Sampling Spin Configurations\n", + "\n", + "**Relevant method:** \n", + "- `DilatedRNNWavefunction.sample(numsamples, inputdim)`\n", + "\n", + "At each iteration, the model generates a batch of complete spin configurations by sampling from the learned distribution. Each sample corresponds to a candidate solution to the Ising optimization problem. Randomness is preserved, but the sampling is biased by the learned model parameters.\n", + "\n", + "---" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18f1722f", + "metadata": {}, + "outputs": [], + "source": [ + "# Number of configurations to sample at each iteration\n", + "num_samples = cfg.numsamples\n", + "\n", + "# IMPORTANT:\n", + "# The following call builds the sampling operations in the\n", + "# TensorFlow 1.x computation graph.\n", + "# `samples` and `log_probs` are symbolic Tensors, not concrete values yet.\n", + "samples, log_probs = rnn_sampler.sample(\n", + " numsamples=num_samples,\n", + " inputdim=2\n", + ")\n", + "\n", + "print(\"Symbolic sample tensor shape:\", samples.shape)\n", + "print(\"Symbolic log-probability tensor shape:\", log_probs.shape)\n", + "\n", + "# In TensorFlow 1.x, tensors must be explicitly evaluated\n", + "# inside a Session to obtain NumPy arrays.\n", + "with tf.compat.v1.Session(graph=rnn_sampler.graph) as sess:\n", + " sess.run(tf.compat.v1.global_variables_initializer())\n", + " \n", + " # Run the sampling operation\n", + " samples_val, log_probs_val = sess.run([samples, log_probs])\n", + "\n", + "# At this point, `samples_val` and `log_probs_val` are NumPy arrays\n", + "print(\"Sample shape (NumPy):\", samples_val.shape)\n", + "print(\"First 5 sampled configurations:\")\n", + "print(samples_val[:5])\n" + ] + }, + { + "cell_type": "markdown", + "id": "8b6e8646", + "metadata": {}, + "source": [ + "\n", + "### 4. Energy Evaluation\n", + "\n", + "**Relevant functions:** \n", + "- `Fullyconnected_diagonal_matrixelements(Jz, samples)` \n", + "- `Fullyconnected_localenergies(...)` (optional)\n", + "\n", + "The quality of sampled configurations is evaluated using the Ising Hamiltonian. Energy values provide the only feedback signal used to guide learning. Lower-energy configurations correspond to better approximate solutions.\n", + "\n", + "---" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a30414a", + "metadata": {}, + "outputs": [], + "source": [ + "# Jz encodes the pairwise couplings of the Ising Hamiltonian\n", + "Jz = cfg.Jz\n", + "\n", + "# Compute classical Ising energies for each sampled configuration\n", + "energies = Fullyconnected_diagonal_matrixelements(Jz, samples_val)\n", + "\n", + "print(\"Energy array shape:\", energies.shape)\n", + "print(\"Mean energy:\", energies.mean())\n", + "print(\"Minimum energy:\", energies.min())" + ] + }, + { + "cell_type": "markdown", + "id": "069dbc2e", + "metadata": {}, + "source": [ + "### 5. Training with Annealing\n", + "\n", + "**Relevant function:** \n", + "- `vca_solver(config)`\n", + "\n", + "Model parameters are updated using energy feedback under a temperature-controlled objective. A warmup phase stabilizes learning, followed by an annealing schedule that gradually shifts the focus from exploration to exploitation. Unlike classical simulated annealing, annealing here regulates the learning objective rather than the sampling dynamics.\n", + "\n", + "Training begins with a **warmup phase**, during which the model stabilizes before annealing is applied. \n", + "This is followed by an **annealing phase**, where the temperature is gradually reduced. At higher temperatures, the learning objective encourages exploration across diverse configurations. As the temperature decreases, updates increasingly focus on concentrating probability mass around low-energy regions.\n", + "\n", + "A key distinction from classical simulated annealing is that annealing in VCA does not control the sampling dynamics. Instead, it regulates the learning objective, acting as a curriculum that guides how the sampling distribution evolves over time." + ] + }, + { + "cell_type": "markdown", + "id": "e0c22d43", + "metadata": {}, + "source": [ + "### 6. Final Sampling for Benchmarking\n", + "\n", + "**Relevant code block:** \n", + "- Final sampling section in `vca_solver`\n", + "\n", + "After the annealing process completes, the trained model is used to generate a large number of spin configurations efficiently. The minimum and average energies of these samples are reported as benchmark metrics, reflecting the practical performance of the learned heuristic." + ] + }, + { + "cell_type": "markdown", + "id": "191dcfe9", + "metadata": {}, + "source": [ + "## 4.3 Execution Cell\n", + "To run the whole solver in this Notebook, use the following code block. Ensure you have the problem instance file (e.g., `10x10_uniform_seed1.txt`) in your working directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6e8f749", + "metadata": {}, + "outputs": [], + "source": [ + "# 1. Path to your problem instance\n", + "instance_path = \"dataset/EA/EA_10x10/10x10_uniform_seed1.txt\" \n", + "seed = 0\n", + "\n", + "# 2. Initialize configuration\n", + "vca_config = config(instance_path, seed)\n", + "\n", + "# 3. Run the solver\n", + "# This will output the annealing progress, energy (E), and free energy (F)\n", + "mean_energies, min_energies = vca_solver(vca_config)\n", + "\n", + "print(f\"\\nFinal Results:\")\n", + "print(f\"Minimum Energy Found: {min_energies}\")\n", + "print(f\"Mean Energy: {mean_energies}\")" + ] + }, + { + "cell_type": "markdown", + "id": "f395513c", + "metadata": {}, + "source": [ + "**Reference**\n", + "\n", + "D. Bertsimas and J. N. Tsitsiklis, *Simulated Annealing*, Statistical Science, vol. 8, no. 1, pp. 10–15, 1993. \n", + "Available at: https://www.mit.edu/~dbertsim/papers/Optimization/Simulated%20annealing.pdf" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (rl4ising)", + "language": "python", + "name": "rl4ising" + }, + "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.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}