The evaluation script for the TSP-GLS problem is provided below.

```python
# Evaluate TSP-GLS
""" 
GLS is a metaheuristic that enhances local search by dynamically penalizing features of poor-quality solutions. For TSP:
Local Search: Uses standard operators (2-opt, relocate) to improve tours
- Penalty Mechanism: When stuck in local optima, identifies "bad" edges (high utility) and penalizes them
- Utility Function: util = guide[u,v] / (1 + penalty[u,v]) - guides which edges to penalize
- Perturbation: Temporarily modifies objective function by adding penalties, forcing search to explore different regions
- Iterative Process: Alternates between local search and perturbation until convergence
- The guide matrix (generated by `heuristics` function) typically contains edge desirability scores (e.g., edge lengths or other heuristics)
"""
import os
import sys
import numpy as np
import numpy.typing as npt
import numba as nb
import logging
from typing import Dict, List, Tuple, Any
from scipy.spatial import distance_matrix
import argparse
import concurrent.futures
import inspect
import seed_solution as solution_module  # Note: solution module script is generated and saved on the fly


# =====Load function to evolve=====
problem = "tsp_gls"
function_to_evolve = "heuristics"
# Robust way to get the name of the function to evolve
def get_function_name(module, possible_names: list[str]):
    """ 
    Args:   
        module: module to search for heuristic function
        possible_names: list of possible names for function to evolve
    
    Used in eval script.
    
    Usage:
        possible_func_names = ["heuristics", "heuristics_v1", "heuristics_v2", "heuristics_v3"]
        heuristic_name = get_function_name(gpt, possible_func_names)
    """
    for func_name in possible_names:
        if hasattr(module, func_name):
            if inspect.isfunction(getattr(module, func_name)):
                return func_name

possible_func_names = [function_to_evolve, function_to_evolve + "_v1", function_to_evolve + "_v2", function_to_evolve + "_v3"]
heuristic_name = get_function_name(solution_module, possible_func_names)
heuristics = getattr(solution_module, heuristic_name)  # Here use the function name


# =====GLS algorithm=====
FloatArray = npt.NDArray[np.float_]
IntArray = npt.NDArray[np.int_]
usecache = True

class TSPInstance:
    def __init__(self, positions: npt.NDArray[np.float_]) -> None:
        self.positions = positions
        self.n = positions.shape[0]
        self.distmat = distance_matrix(positions, positions) + np.eye(self.n)*1e-5

@nb.njit(nb.float32(nb.float32[:,:], nb.uint16[:], nb.uint16), nogil=True, cache = usecache)
def _two_opt_once(distmat, tour, fixed_i = 0):
    '''Perform a single 2-opt move (swap two edges) to improve the tour.

    Args:
        distmat: Distance matrix between cities
        tour: Current tour (modified in-place)
        fixed_i: If non-zero, only consider moves involving this position

    Returns:
        delta: Cost improvement (negative if improvement found)
    '''
    n = tour.shape[0]
    p = q = 0
    delta = 0
    # Search for best 2-opt move
    for i in range(1, n - 1) if fixed_i==0 else range(fixed_i, fixed_i+1):
        for j in range(i + 1, n):
            node_i, node_j = tour[i], tour[j]
            node_prev, node_next = tour[i-1], tour[(j+1) % n]
            if node_prev == node_j or node_next == node_i:
                continue
            # Calculate cost change: remove edges (i-1,i) and (j,j+1), add edges (i-1,j) and (i,j+1)
            change = (  distmat[node_prev, node_j]
                        + distmat[node_i, node_next]
                        - distmat[node_prev, node_i]
                        - distmat[node_j, node_next])
            if change < delta:
                p, q, delta = i, j, change
    if delta < -1e-6:
        # Apply the 2-opt move: reverse segment between positions p and q
        tour[p: q+1] = np.flip(tour[p: q+1])
        return delta
    else:
        return 0.0

@nb.njit(nb.float32(nb.float32[:,:], nb.uint16[:], nb.uint16), nogil=True, cache = usecache)
def _relocate_once(distmat, tour, fixed_i = 0):
    '''Perform a single relocate move (move a city to different position) to improve the tour.

    Args:
        distmat: Distance matrix between cities
        tour: Current tour (modified in-place)
        fixed_i: If non-zero, only consider moves involving this position

    Returns:
        delta: Cost improvement (negative if improvement found)
    '''
    n = distmat.shape[0]
    delta = p = q = 0
    # Search for best relocate move
    for i in range(1, n) if fixed_i==0 else range(fixed_i, fixed_i+1):
        node = tour[i]
        prev_node = tour[i-1]
        next_node = tour[(i+1)%n]
        for j in range(n):
            if j == i or j == i-1:
                continue
            prev_insert = tour[j]
            next_insert = tour[(j+1)%n]
            # Calculate cost change: remove node from position i, insert between positions j and j+1
            cost = ( - distmat[prev_node, node]
                     - distmat[node, next_node]
                     - distmat[prev_insert, next_insert]
                     + distmat[prev_insert, node]
                     + distmat[node, next_insert]
                     + distmat[prev_node, next_node] )
            if cost < delta:
                delta, p, q = cost, i, j
    if delta >= 0:
        return 0.0
    # Apply the relocate move
    if p<q:
        tour[p:q+1] = np.roll(tour[p:q+1], -1)  # Move node p to position q
    else:
        tour[q:p+1] = np.roll(tour[q:p+1], 1)   # Move node p to position q
    return delta

@nb.njit(nb.float32(nb.float32[:,:], nb.uint16[:]), nogil=True, cache = usecache)
def _calculate_cost(distmat, tour):
    '''Calculate total distance of a TSP tour.

    Args:
        distmat: Distance matrix between cities
        tour: Tour as sequence of city indices

    Returns:
        Total tour distance
    '''
    cost = distmat[tour[-1], tour[0]]  # Closing edge from last to first city
    for i in range(len(tour) - 1):
        cost += distmat[tour[i], tour[i+1]]
    return cost

@nb.njit(nb.float32(nb.float32[:,:], nb.uint16[:], nb.uint16, nb.uint16), nogil=True, cache = usecache)
def _local_search(distmat, cur_tour, fixed_i = 0, count = 1000):
    '''Perform local search using 2-opt and relocate moves until no improvement.

    Args:
        distmat: Distance matrix between cities
        cur_tour: Current tour (modified in-place)
        fixed_i: If non-zero, only consider moves involving this position
        count: Maximum number of iterations

    Returns:
        Total cost improvement achieved
    '''
    sum_delta = 0.0
    delta = -1
    # Continue while improvements are found and iteration limit not reached
    while delta < 0 and count > 0:
        delta = 0
        delta += _two_opt_once(distmat, cur_tour, fixed_i)
        delta += _relocate_once(distmat, cur_tour, fixed_i)
        count -= 1
        sum_delta += delta
    return sum_delta

@nb.njit(nb.void(nb.float32[:,:], nb.float32[:,:], nb.float32[:,:], nb.uint16[:], nb.float32, nb.uint32), nogil=True, cache = usecache)
def _perturbation(distmat, guide, penalty, cur_tour, k, perturbation_moves = 30):
    '''Perform perturbation phase of GLS: penalize "bad" edges and escape local optima.

    Args:
        distmat: Distance matrix between cities
        guide: Guide matrix (e.g., edge lengths) used to compute utility
        penalty: Current penalty matrix (modified in-place)
        cur_tour: Current tour (modified in-place)
        k: Penalty scaling factor
        perturbation_moves: Number of perturbation moves to perform
    '''
    moves = 0
    n = distmat.shape[0]
    while moves < perturbation_moves:
        # Find edge with maximum utility (guide/(1+penalty)) in current tour
        max_util = 0
        max_util_idx = 0
        for i in range(n-1):
            j = i+1
            u, v = cur_tour[i], cur_tour[j]
            # Utility function: guide value divided by (1 + current penalty)
            # Higher utility means edge is more "undesirable" and should be penalized
            util = guide[u, v] / (1.0 + penalty[u, v])
            if util > max_util:
                max_util_idx, max_util = i, util

        # Penalize the edge with maximum utility
        penalty[cur_tour[max_util_idx], cur_tour[max_util_idx+1]] += 1.0
        # Create modified distance matrix with penalties
        edge_weight_guided = distmat + k * penalty

        # Try local search around the penalized edge to escape local optimum
        for fixed_i in (max_util_idx, max_util_idx+1):
            if fixed_i == 0 or fixed_i + 1 == n:
                continue
            delta = _local_search(edge_weight_guided, cur_tour, fixed_i, 1)
            if delta < 0:
                moves += 1  # Count successful perturbation moves

@nb.njit(nb.uint16[:](nb.float32[:,:], nb.uint16), nogil=True, cache = usecache)
def _init_nearest_neighbor(distmat, start):
    '''Construct initial tour using nearest neighbor heuristic.

    Args:
        distmat: Distance matrix between cities
        start: Starting city index

    Returns:
        Initial tour starting from given city
    '''
    n = distmat.shape[0]
    tour = np.zeros(n, dtype=np.uint16)
    visited = np.zeros(n, dtype=np.bool_)
    visited[start] = True
    tour[0] = start
    for i in range(1, n):
        min_dist = np.inf
        min_idx = -1
        # Find nearest unvisited city to the last city in tour
        for j in range(n):
            if not visited[j] and distmat[tour[i-1], j] < min_dist:
                min_dist = distmat[tour[i-1], j]
                min_idx = j
        tour[i] = min_idx
        visited[min_idx] = True
    return tour

@nb.njit(nb.uint16[:](nb.float32[:,:], nb.float32[:,:], nb.uint16, nb.int32, nb.uint16), nogil = True, cache = usecache)
def _guided_local_search(
    distmat, guide, start, perturbation_moves = 30, iter_limit = 1000
) -> npt.NDArray[np.uint16]:
    '''Core GLS algorithm: iteratively apply local search and perturbation.

    Args:
        distmat: Distance matrix between cities
        guide: Guide matrix for utility calculation
        start: Starting city index
        perturbation_moves: Number of perturbation moves per iteration
        iter_limit: Maximum number of GLS iterations

    Returns:
        Best tour found
    '''
    penalty = np.zeros_like(distmat)  # Initialize penalty matrix

    # Initial solution using nearest neighbor + local search
    best_tour = _init_nearest_neighbor(distmat, start)
    _local_search(distmat, best_tour, 0, 1000)
    best_cost = _calculate_cost(distmat, best_tour)
    # Compute penalty scaling factor: 0.1 * (average edge cost)
    k = 0.1 * best_cost / distmat.shape[0]
    cur_tour = best_tour.copy()

    # Main GLS loop
    for _ in range(iter_limit):
        # Perturbation phase: penalize edges and escape local optima
        _perturbation(distmat, guide, penalty, cur_tour, k, perturbation_moves)
        # Local search phase: improve solution with modified penalties
        _local_search(distmat, cur_tour, 0, 1000)
        cur_cost = _calculate_cost(distmat, cur_tour)
        # Update best solution if improved
        if cur_cost < best_cost:
            best_tour, best_cost = cur_tour.copy(), cur_cost
    return best_tour

def guided_local_search(
    distmat: FloatArray,
    guide: FloatArray,
    perturbation_moves: int = 30,
    iter_limit: int = 1000
) -> npt.NDArray[np.uint16]:
    '''Python wrapper for GLS algorithm with type conversion.

    Args:
        distmat: Distance matrix between cities
        guide: Guide matrix for utility calculation
        perturbation_moves: Number of perturbation moves per iteration
        iter_limit: Maximum number of GLS iterations

    Returns:
        Best tour found
    '''
    return _guided_local_search(
        distmat = distmat.astype(np.float32),
        guide = guide.astype(np.float32),
        start = 0,
        perturbation_moves = perturbation_moves,
        iter_limit = iter_limit,
    )

def multi_start_guided_local_search(
    dist: FloatArray,
    guide: FloatArray,
    n_starts: int = 10,
    perturbation_moves = 30,
    iter_limit = 1000
):
    '''Run GLS from multiple starting points in parallel and return best solution.

    Args:
        dist: Distance matrix between cities
        guide: Guide matrix for utility calculation
        n_starts: Number of different starting cities
        perturbation_moves: Number of perturbation moves per iteration
        iter_limit: Maximum number of GLS iterations

    Returns:
        Best tour found across all starting points
    '''
    dist = dist.astype(np.float32)
    guide = guide.astype(np.float32)
    start_nodes = np.arange(n_starts).astype(np.uint16)

    # Run GLS from different starting points in parallel
    with concurrent.futures.ThreadPoolExecutor() as executor:
        futures = []
        for start in start_nodes:
            future = executor.submit(_guided_local_search, dist, guide, start, perturbation_moves = perturbation_moves, iter_limit = iter_limit)
            futures.append(future)
        tours = [f.result() for f in futures]
        # Calculate costs and return the best tour
        costs = np.array([_calculate_cost(dist, tour) for tour in tours])
        best_tour = tours[np.argmin(costs)]
        return best_tour


# =====Evaluation function=====
perturbation_moves = 30
iter_limit = 1200

def calculate_cost(inst: TSPInstance, path: np.ndarray) -> float:
    return inst.distmat[path, np.roll(path, 1)].sum().item()

def eval_heuristic(inst: TSPInstance) -> float:
    heu = heuristics(inst.distmat.copy())
    assert tuple(heu.shape) == (inst.n, inst.n)
    result = guided_local_search(inst.distmat, heu, perturbation_moves, iter_limit)
    # print(result)
    return calculate_cost(inst, result)


# =====Helper functions=====
def get_feature(metrics: Dict[int, float]) -> Tuple[int, ...]:
    """
    Convert the metrics dict to a feature vector

    Args:
        metrics (dict): A mapping of test problem size (int) to a score (float).

    Returns:
        (tuple): a tuple of discretized scores sorted by problem size
    """
    scores = metrics.values()
    features = tuple([int(x) for x in scores])
    return features
    
def get_score(metrics: Dict[int, float]) -> float:
    """
    Convert the metrics dict to a score

    Args:
        metrics (dict): A mapping of test problem size (int) to a score (float).

    Returns:
        (float): a score
    """
    return sum(metrics.values()) / len(metrics)


# =====Main Function=====
if __name__ == '__main__':
    print("TSP-constructive evaluation script running...")
    # -----Parse command line arguments (same for all problems)-----
    parser = argparse.ArgumentParser(description='Evaluation script.')
    parser.add_argument(
        '--root_dir',
        type=str,
        default=os.getcwd(),
        help='Project root directory for loading data (default: current working directory)'
    )
    parser.add_argument(
        '--file_output_prefix',
        type=str,
        default='',
        help='Output file prefix for saving evaluation results. '
             'Absolute path recommended. Files saved as {prefix}filename '
             '(default: empty string, saves to current directory)')
    parser.add_argument(
        '--mode',
        type=str,
        default='val',
        choices=['train', 'val'],
        help='Execution mode: train or val (default: val)'
    )
    parser.add_argument(
        '--problem_size',
        type=int,
        default=50,  # Customize this to your needs
        help='Problem size parameter'
    )
    # Parse arguments
    args = parser.parse_args()
    root_dir = args.root_dir
    file_output_prefix = args.file_output_prefix
    mode = args.mode
    problem_size = args.problem_size
    # Print parsed arguments for verification
    print(f"root_dir: {root_dir}")
    print(f"file_output_prefix: {file_output_prefix}")
    print(f"mode: {mode}")
    #print(f"problem_size: {problem_size}")
    
    # -----Run the evaluation-----
    try:
        basepath = os.path.join(root_dir, "problems", problem , "dataset")
        # Initialize performance and feature variables
        metrics = {}
        # ---Train mode---
        if mode == 'train':
            # Load dataset
            dataset_path = os.path.join(basepath, f"train{problem_size}_dataset.npy")
            node_positions = np.load(dataset_path)
            n_instances = node_positions.shape[0]  # data shape: (n_instances, problem_size, 2)
            print(f"[*] Dataset loaded: {dataset_path} with {n_instances} instances.")
            
            objs = []
            for i in range(n_instances):
                # Invoke evaluation function
                inst = TSPInstance(node_positions[i])
                obj = eval_heuristic(inst)
                print(f"[*] Instance {i}: {obj}")
                objs.append(obj)
            
            print("[*] Average:")
            print(np.mean(objs))
            metrics[problem_size] = np.mean(objs)
        # ---Val mode---
        else:
            for problem_size in [20, 50, 100]:
                dataset_path = os.path.join(basepath, f"val{problem_size}_dataset.npy")
                logging.info(f"[*] Evaluating {dataset_path}")
                node_positions = np.load(dataset_path)
                n_instances = node_positions.shape[0]
                objs = []
                for i in range(n_instances):
                    inst = TSPInstance(node_positions[i])
                    obj = eval_heuristic(inst)
                    objs.append(obj)
                print(f"[*] Average for {problem_size} cities: {np.mean(objs)}")
                metrics[problem_size] = np.mean(objs)
                
        if metrics:
            features = get_feature(metrics)
            score = get_score(metrics)
        else:
            features = None
            score = None

        # -----Print results to stdout (same for all problems)-----
        print('__SANDBOX_RESULT__')        
        print('__METRICS_START__')
        print(repr(metrics))
        print('__METRICS_END__')
        
        print('__FEATURES_START__')
        print(repr(features))
        print('__FEATURES_END__')
        
        print('__SCORE_START__')
        print(repr(score))
        print('__SCORE_END__')
        
        print('__SANDBOX_SUCCESS__')
        
    except Exception as e:
        import traceback
        print('__SANDBOX_ERROR__:')
        print(f'Error type: {type(e).__name__}')
        print(f'Error message: {str(e)}')
        print('Full traceback:')
        traceback.print_exc()
```