runner.py

#
# Copyright 2026 Huawei Technologies Co., Ltd.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#

"""
Main hardware microbenchmarking runner.

This script provides the main entry point for running hardware microbenchmarks
and uses the modular components for parsing, statistics, plotting, and benchmarking.
"""

import os
import subprocess
import math
import numpy as np
import pandas as pd
import datetime
import argparse
import logging
import sys
import shlex
from pathlib import Path
from hw_characterization.scripts.parser import log_info, _require_env, consecutive_bytes_to_elemsize, format_bytes
from hw_characterization.scripts.numa_utils import expand_cpu_ranges
from hw_characterization.scripts.regression_analysis import least_squares_fit
from hw_characterization.scripts.plotters import plot_all_regression_stats, plot_all_regression_stats_gbytes, plot_PaHiS_submodel
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# Global configuration
MAX_BUFF_SIZE = 1 * 1024 * 1024 * 1024  # 1 GB

# Constants for unit conversions
BYTES_TO_GB = 1024**3  # Convert bytes to GB
BYTES_TO_GB_PER_SEC = 1e9  # Convert bytes/second to GB/s

# Map of microbenchmark kernels to their properties:
KERNEL_NAMES = ["chase", "sync", "copy", "daxpy", "dot", "triad"]
KERNEL_NAME_TO_IDX = {name: idx for idx, name in enumerate(KERNEL_NAMES)}
# - Number of IO vectors
KERNEL_VECNUM = [1, 0, 2, 2, 2, 3]
# - Number of reads
KERNEL_READS  = [1, 1, 1, 1, 2, 2]
# - Number of writes
KERNEL_WRITES = [1, 1, 1, 1, 0, 1]

DEFAULT_BW_KERNELS = "copy,daxpy,triad"  # default selection for --bw-kernels

def aggregate_values(values, method='mean'):
    """Aggregate an array/list of values using the specified method."""
    if method == 'mean':
        return float(np.mean(values))
    elif method == 'min':
        return float(np.min(values))
    elif method == 'max':
        return float(np.max(values))
    else:
        raise ValueError(f"Unknown aggregation method: {method}")

# Global variables for logging and file tracking
logger = None
file_tracker = set()
step_summaries = []
CPU_ALLOWED_FROM_TOPO = None  # Optional CPU allow-list from topology file

def setup_logging(verbose=False):
    """Setup dual logging: detailed to runner.log, summary to stdout"""
    global logger
    
    # Get build directory for log file
    build_dir = os.environ.get("HW_BUILD_DIR", "./build")
    log_file = os.path.join(build_dir, "runner.log")
    
    # Create logger
    logger = logging.getLogger('runner')
    logger.setLevel(logging.DEBUG)
    
    # Clear any existing handlers
    for handler in logger.handlers[:]:
        logger.removeHandler(handler)
    
    # File handler for detailed logging
    file_handler = logging.FileHandler(log_file, mode='w')
    file_handler.setLevel(logging.DEBUG)
    file_formatter = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')
    file_handler.setFormatter(file_formatter)
    logger.addHandler(file_handler)
    
    # Console handler for summary output (or debug if verbose)
    console_handler = logging.StreamHandler(sys.stdout)
    console_handler.setLevel(logging.DEBUG if verbose else logging.INFO)
    console_formatter = logging.Formatter('%(message)s')
    console_handler.setFormatter(console_formatter)
    logger.addHandler(console_handler)
    
    return log_file

def log_summary(step_name, status, details=None, files_generated=None):
    """Log a summary step to both file and console"""
    global step_summaries
    
    summary = {
        'step': step_name,
        'status': status,
        'timestamp': datetime.datetime.now().isoformat(),
        'details': details or {},
        'files': files_generated or []
    }
    step_summaries.append(summary)
    
    # Log to file (detailed)
    logger.debug(f"STEP: {step_name}")
    logger.debug(f"STATUS: {status}")
    if details:
        for key, value in details.items():
            logger.debug(f"  {key}: {value}")
    if files_generated:
        logger.debug(f"FILES GENERATED: {len(files_generated)}")
        for file_path in files_generated:
            logger.debug(f"  - {file_path}")
    
    # Log to console (summary)
    status_symbol = "✓" if status == "SUCCESS" else "✗" if status == "ERROR" else "⚠"
    logger.info(f"{status_symbol} {step_name}")
    if details and len(details) <= 3:  # Only show key details in console
        for key, value in details.items():
            logger.info(f"    {key}: {value}")
    if files_generated and len(files_generated) <= 5:  # Only show first few files
        logger.info(f"    Files: {len(files_generated)} generated")
        for file_path in files_generated[:3]:
            logger.info(f"      - {file_path}")
        if len(files_generated) > 3:
            logger.info(f"      ... and {len(files_generated) - 3} more")

def track_file(file_path):
    """Track a generated file"""
    global file_tracker
    if file_path:
        file_tracker.add(str(file_path))

def get_generated_files_summary():
    """Get summary of all generated files organized by type"""
    if not file_tracker:
        return {}
    
    files_by_type = {}
    for file_path in file_tracker:
        path = Path(file_path)
        if path.exists():
            if path.suffix == '.csv':
                files_by_type.setdefault('CSV Data', []).append(file_path)
            elif path.suffix == '.png':
                files_by_type.setdefault('Plots', []).append(file_path)
            elif path.suffix == '.hpp':
                files_by_type.setdefault('C++ Headers', []).append(file_path)
            elif path.suffix == '.log':
                files_by_type.setdefault('Logs', []).append(file_path)
            else:
                files_by_type.setdefault('Other', []).append(file_path)
    
    return files_by_type

# Directory configuration - hardcoded paths + configurable ones
TOPOLOGY_DIR = "hw_characterization/topologies"  # Hardcoded
RESULTS_DIR = "hw_characterization/results"  # Hardcoded  
CPP_MODEL_DIR = f"hw_characterization/results/{os.environ['HW_MODEL']}/output_headers"  # Hardcoded with HW_MODEL
BUILD_DIR = os.environ.get("HW_BUILD_DIR", "./build")  # Configurable

# Check if required environment variables are set
if not os.environ.get("HW_MODEL"):
    print("\n" + "="*80)
    print("ERROR: HW_MODEL environment variable not set!")
    print("="*80)
    print("Please source your system configuration file:")
    print("  source config_${HW_MODEL}.sh")
    print("Or set HW_MODEL manually:")
    print("  export HW_MODEL=your_system_name")
    print("="*80 + "\n")
    exit(1)

# Check if ALP_PATH is set (optional for basic functionality)
ALP_PATH = os.environ.get("ALP_PATH")
if not ALP_PATH:
    print("\n" + "="*80)
    print("WARNING: ALP_PATH environment variable not set!")
    print("="*80)
    print("The runner will work for hardware characterization, but ALP-dependent")
    print("features will be disabled. To enable full functionality:")
    print("  source config_${HW_MODEL}.sh")
    print("Or set ALP_PATH manually:")
    print("  export ALP_PATH=/path/to/ALP")
    print("="*80 + "\n")
else:
    print(f"ALP_PATH set to: {ALP_PATH}")


class BenchmarkConfig:
    """
    Configuration class to encapsulate common benchmark parameters.
    
    This class stores all the parameters that are passed around between functions,
    reducing the number of individual parameters needed in function signatures.
    """
    
    def __init__(self, system_name, consecutive_bytes, num_threads, cacheline_size=64, 
                 thread_offset=0, forced_intercept=-1, level_name=None, 
                 create_sync_level="", alloc_policy="close", 
                 numa_level_reached=0, system_cores=0, bw_kernel_indices=None,
                 kernel_aggregator='mean'):
        """
        Initialize benchmark configuration.
        
        Args:
            system_name: Name of the system being benchmarked
            consecutive_bytes: Consecutive byte size for benchmarks
            num_threads: Number of threads to use
            cacheline_size: Cache line size in bytes (default: 64)
            thread_offset: Thread offset for CPU affinity (default: 0)
            forced_intercept: Forced intercept value for regression (default: -1)
            level_name: Name of the memory level (optional)
            create_sync_level: GLOBAL_SYNC level creation method (default: "")
            alloc_policy: NUMA allocation policy (default: "close")
            numa_level_reached: The level of the NUMA node that we've reached (default: 0)
            system_cores: Total number of cores in the system (default: 0)
            bw_kernel_indices: List of kernel indices to use for bandwidth measurements
                               (default: None → uses copy, daxpy, triad)
            kernel_aggregator: Method to aggregate across kernel measurements ('mean', 'min', 'max')
        """
        self.system_name = system_name
        self.consecutive_bytes = consecutive_bytes
        self.num_threads = num_threads
        self.cacheline_size = cacheline_size
        self.thread_offset = thread_offset
        self.forced_intercept = forced_intercept
        self.level_name = level_name
        self.create_sync_level = create_sync_level
        self.alloc_policy = alloc_policy
        self.numa_level_reached = numa_level_reached
        self.system_cores = system_cores
        self.bw_kernel_indices = bw_kernel_indices if bw_kernel_indices is not None else [
            KERNEL_NAME_TO_IDX[k] for k in DEFAULT_BW_KERNELS.split(',')
        ]
        self.kernel_aggregator = kernel_aggregator
        
        # Derived properties
        self.elemsize = consecutive_bytes_to_elemsize(consecutive_bytes)
        self.chunk_elems = consecutive_bytes // self.elemsize
        self.cacheline_elems = cacheline_size // self.elemsize
        
    def __str__(self):
        """String representation of the configuration."""
        return (f"BenchmarkConfig(system={self.system_name}, "
                f"k={self.consecutive_bytes}B, threads={self.num_threads}, "
                f"cacheline={self.cacheline_size}B, level={self.level_name}, "
                f"alloc_policy={self.alloc_policy}, numa_level={self.numa_level_reached}, "
                f"system_cores={self.system_cores})")
    
    def __repr__(self):
        """Detailed string representation of the configuration."""
        return (f"BenchmarkConfig(system_name='{self.system_name}', "
                f"consecutive_bytes={self.consecutive_bytes}, "
                f"num_threads={self.num_threads}, "
                f"cacheline_size={self.cacheline_size}, "
                f"thread_offset={self.thread_offset}, "
                f"forced_intercept={self.forced_intercept}, "
                f"level_name='{self.level_name}', "
                f"create_sync_level={self.create_sync_level}, "
                f"alloc_policy='{self.alloc_policy}', "
                f"numa_level_reached={self.numa_level_reached}, "
                f"system_cores={self.system_cores})")
    
    def copy_with(self, **kwargs):
        """
        Create a copy of this config with some parameters overridden.
        
        Args:
            **kwargs: Parameters to override in the new config
            
        Returns:
            BenchmarkConfig: New configuration with overridden parameters
        """
        new_config = BenchmarkConfig(
            system_name=self.system_name,
            consecutive_bytes=self.consecutive_bytes,
            num_threads=self.num_threads,
            cacheline_size=self.cacheline_size,
            thread_offset=self.thread_offset,
            forced_intercept=self.forced_intercept,
            level_name=self.level_name,
            create_sync_level=self.create_sync_level,
            alloc_policy=self.alloc_policy,
            numa_level_reached=self.numa_level_reached,
            system_cores=self.system_cores,
            bw_kernel_indices=self.bw_kernel_indices,
            kernel_aggregator=self.kernel_aggregator
        )
        
        for key, value in kwargs.items():
            if hasattr(new_config, key):
                setattr(new_config, key, value)
            else:
                raise ValueError(f"Unknown parameter: {key}")
        
        # Recalculate derived properties
        new_config.elemsize = consecutive_bytes_to_elemsize(new_config.consecutive_bytes)
        new_config.chunk_elems = new_config.consecutive_bytes // new_config.elemsize
        new_config.cacheline_elems = new_config.cacheline_size // new_config.elemsize
        
        return new_config


def process_benchmark_results(h, av_timers, min_timers, max_timers, bench_idx, config):
    """
    Process benchmark results and generate regression plots.
    
    Args:
        h: Array of working set sizes
        av_timers: Array of average timer measurements
        min_timers: Array of minimum timer measurements
        max_timers: Array of maximum timer measurements
        bench_idx: Benchmark index (0=chase, 1=sync, 2=copy, 3=daxpy, 4=dot, 5=triad)
        config: BenchmarkConfig object containing benchmark parameters
        
    Returns:
        tuple: (gi_avg, Li_avg) - bandwidth and latency parameters
    """
    h_datumsize = KERNEL_VECNUM[bench_idx] * (config.cacheline_size if config.consecutive_bytes <= config.cacheline_size else config.consecutive_bytes)
    h_footprint = np.array([h[i] * h_datumsize for i in range(len(h))])
    h_bytes_read = np.array([h[i] * KERNEL_READS[bench_idx] * config.consecutive_bytes for i in range(len(h))])
    h_bytes_written = np.array([h[i] * KERNEL_WRITES[bench_idx] * config.consecutive_bytes for i in range(len(h))])
    h_total = np.array([h_bytes_read[i] + h_bytes_written[i] for i in range(len(h))])

    if config.forced_intercept == -42:
        h_parsed, av_timers_parsed, gi_avg, Li_avg, h_removed, timers_removed =\
            [1], av_timers, 0, av_timers[0], [], []
    else:
        h_parsed, av_timers_parsed, gi_avg, Li_avg, h_removed, timers_removed = least_squares_fit(
            h_total, av_timers, LR_reps=1, forced_intercept=config.forced_intercept) 
        print("Least squares fit: gi_avg = %e (%.2lf Gb/s), Li_avg = %e" % (gi_avg, 1 / gi_avg / BYTES_TO_GB if gi_avg != 0 else 0, Li_avg))

    regression_plot_flag = True
    if regression_plot_flag:
        # Create combined time plot
        import seaborn as sns
        import matplotlib.pyplot as plt
        sns.set_style("whitegrid")
        plt.figure(figsize=(12, 8))

        KB_memfoot_converter = lambda h_val:  1.0 * h_val / ((KERNEL_READS[bench_idx] + KERNEL_WRITES[bench_idx])* config.consecutive_bytes) \
        * h_datumsize / 1024

        # Plot the fitted data
        plot_all_regression_stats(h_parsed, av_timers_parsed, gi_avg, Li_avg, KB_memfoot_converter, label_prefix="Avg", line_color="blue")

        # Plot the removed outliers if any
        if len(h_removed) > 0:
            h_removed_KB = [KB_memfoot_converter(h_val) for h_val in h_removed]
            plt.scatter(h_removed_KB, timers_removed, color='red', s=100, alpha=0.8, 
                    marker='x', linewidths=3, label='Removed outliers (avg)', zorder=5)
        
        plt.title(f"Benchmark: {KERNEL_NAMES[bench_idx]} | Level: {config.level_name} | Threads: {config.num_threads}", fontsize=14, fontweight='bold')
        plt.legend()

        save_path = f"{RESULTS_DIR}/{config.system_name}/{config.alloc_policy}/k_{config.consecutive_bytes}/t_{config.num_threads}/plots/time/regression_{KERNEL_NAMES[bench_idx]}_{config.level_name}.png"
        os.makedirs(os.path.dirname(save_path), exist_ok=True)
        plt.savefig(save_path, dpi=300, bbox_inches='tight')
        plt.close()
        print(f"Time regression plot saved to: {save_path}")
        log_info(f"Saved time regression plot: {save_path}")
        
        # Create combined GBytes plot
        plt.figure(figsize=(12, 8))

        gbytes_converter = lambda h_val, t_val:  h_val / t_val / BYTES_TO_GB

        plot_all_regression_stats_gbytes(h_parsed, av_timers_parsed, gi_avg, Li_avg, KB_memfoot_converter, gbytes_converter, label_prefix="Avg", line_color="blue")

        # Plot the removed outliers in GBytes if any
        if len(h_removed) > 0:
            gbytes_removed = []
            h_to_KB_valid = []

            for h_val, t_val in zip(h_removed, timers_removed):
                if t_val > 0:  # Avoid division by zero
                    h_to_KB_valid.append(KB_memfoot_converter(h_val))
                    gbytes_removed.append(gbytes_converter(h_val, t_val))

            if len(gbytes_removed) > 0:
                plt.scatter(h_to_KB_valid, gbytes_removed, color='red', s=100, alpha=0.8, 
                        marker='x', linewidths=3, label='Removed outliers (avg)', zorder=5)

        plt.title(f"Benchmark: {KERNEL_NAMES[bench_idx]} | Level: {config.level_name} | Threads: {config.num_threads}", fontsize=14, fontweight='bold')
        plt.legend()

        save_path = f"{RESULTS_DIR}/{config.system_name}/{config.alloc_policy}/k_{config.consecutive_bytes}/t_{config.num_threads}/plots/gbytes/regression_{KERNEL_NAMES[bench_idx]}_{config.level_name}_gbytes.png"
        os.makedirs(os.path.dirname(save_path), exist_ok=True)
        plt.savefig(save_path, dpi=300, bbox_inches='tight')
        plt.close()
        log_info(f"Saved bandwidth regression plot: {save_path}")
    
    return gi_avg, Li_avg

def calculate_spread_cpu_affinity(num_threads, config):
    """
    Calculate CPU affinity for spread policy to distribute threads across NUMA nodes.
    Uses the numa_utils module to get actual NUMA topology.
    
    Args:
        num_threads: Number of threads
        config: BenchmarkConfig object containing system information
        
    Returns:
        list: List of CPU IDs for thread affinity
    """
    from hw_characterization.scripts.numa_utils import get_numa_topology, get_cores_per_numa_node
    num_numa_nodes = get_numa_topology()['num_nodes']
    cores_per_numa = get_cores_per_numa_node()
    print(f"Spread policy: {num_numa_nodes} NUMA nodes, {cores_per_numa} cores per node")
    
    # Calculate CPU affinity for spread policy
    cpu_affinity = []
    threads_per_node = num_threads // num_numa_nodes
    residual = num_threads % num_numa_nodes
    
    # Calculate how many threads each NUMA node will have
    # First 'residual' nodes get (threads_per_node + 1) threads, rest get threads_per_node
    threads_per_node_list = [threads_per_node + 1 if i < residual else threads_per_node 
                             for i in range(num_numa_nodes)]
    
    # Build CPU affinity using block distribution (consecutive threads per NUMA node)
    # Extra threads are already accounted for in threads_per_node_list
    thread_counter = 0
    
    for numa_node in range(num_numa_nodes):
        # Assign consecutive threads to this NUMA node
        threads_for_this_node = threads_per_node_list[numa_node]
        for thread_in_node in range(threads_for_this_node):
            cpu_id = config.thread_offset + numa_node * cores_per_numa + thread_in_node
            cpu_affinity.append(cpu_id)
            thread_counter += 1
    
    print(f"Spread CPU affinity: {cpu_affinity}")
    return cpu_affinity


def PaHiS_benchmark_run(bench_idx, buffer_size, sub_buffer_size, config):
    """
    Run PaHiS benchmark and process results.
    
    Args:
        bench_idx: Benchmark index (0=chase, 1=sync, 2=copy, 3=daxpy, 4=dot, 5=triad)
        buffer_size: Buffer size in bytes
        sub_buffer_size: Sub-buffer size in bytes
        config: BenchmarkConfig object containing benchmark parameters
        
    Returns:
        tuple: (gi_avg, Li_avg) - bandwidth and latency parameters
    """
    # Use first num_threads cores for OMP_PLACES, formatted as "{0,1,2,...}"
    threads = config.num_threads

    filename = f"{RESULTS_DIR}/{config.system_name}/{config.alloc_policy}/k_{config.consecutive_bytes}/t_{config.num_threads}/" + \
    f"results/{config.level_name}_{KERNEL_NAMES[bench_idx]}_buf{int(buffer_size)}_sub{int(sub_buffer_size)}_t{config.num_threads}.csv"
    
    # Check if file already exists
    if os.path.exists(filename):
        print(f"Loading existing results from: {filename}")
        try:
            df = pd.read_csv(filename)
            h = df['h'].tolist() # (df['h'] * threads).tolist() #
            av_timers = df['av_timers'].tolist()
            min_timers = df['min_timers'].tolist()
            max_timers = df['max_timers'].tolist()
            
            print(f"Loaded {len(h)} data points from existing file")
            
            # Process the loaded data
            return process_benchmark_results(h, av_timers, min_timers, max_timers, 
                                           bench_idx, config)
        except Exception as e:
            print(f"Error loading existing file: {e}, running new benchmark")

    env = os.environ.copy()
    env["OMP_NUM_THREADS"] = str(threads)
    # Build requested affinity based on policy
    if config.alloc_policy == "spread":
        requested = calculate_spread_cpu_affinity(threads, config)
    else:
        requested = [config.thread_offset + i for i in range(threads)]

    # If topology provided CPU_ALLOWED_LIST, filter/pad from it; else keep requested
    global CPU_ALLOWED_FROM_TOPO
    if CPU_ALLOWED_FROM_TOPO:
        allowed = set(expand_cpu_ranges(CPU_ALLOWED_FROM_TOPO))
        adjusted = [c for c in requested if c in allowed]
        if len(adjusted) < threads:
            adjusted += [c for c in sorted(allowed) if c not in adjusted][:threads - len(adjusted)]
    else:
        adjusted = requested

    env["GOMP_CPU_AFFINITY"] = " ".join(str(cpu_id) for cpu_id in adjusted)
    cmd = [
        f"{BUILD_DIR}/bin/PaHiS_benchmark_numa_ELEM_SIZE-{config.elemsize}_CHUNK_ELEMS-{config.chunk_elems}_CACHELINE_ELEMS-{config.cacheline_elems}",
        str(bench_idx),
        str(int(buffer_size)),
        str(int(sub_buffer_size)),
        str(threads),
        config.alloc_policy,
        str(config.numa_level_reached if config.numa_level_reached > 1 else 0)
    ]
    print("Running:", " ".join(cmd))
    print("OMP_NUM_THREADS:", env["OMP_NUM_THREADS"])
    print("GOMP_CPU_AFFINITY:", env["GOMP_CPU_AFFINITY"])
    print(
        f"NUMA parameters: alloc_policy={config.alloc_policy}, interleave_flag={config.numa_level_reached if config.numa_level_reached > 1 else 0}")
    log_info(
        f"Launching NUMA benchmark kernel_idx={bench_idx} level={config.level_name} threads={threads} k={config.consecutive_bytes}B alloc_policy={config.alloc_policy} interleave_flag={config.numa_level_reached if config.numa_level_reached > 1 else 0}")
    result = subprocess.run(cmd, env=env, capture_output=True, text=True)
    
    # Parse bandwidth from output
    h = []
    av_timers = []
    min_timers = []
    max_timers = []
    new_line_counter = 0
    for line in result.stdout.splitlines():
        if new_line_counter == 1:
            min_timers.append(float(line.split("Min=")[1].split("s")[0].strip()))
            new_line_counter += 1
        elif new_line_counter == 2:
            max_timers.append(float(line.split("Max=")[1].split("s")[0].strip()))
            new_line_counter += 1
        elif new_line_counter == 3:
            av_timers.append(float(line.split("Avg=")[1].split("s")[0].strip()))
            new_line_counter = 0
        elif "\th:" in line:
            h.append(threads*int(line.split("h:")[1].split(",")[0].strip()))
            new_line_counter = 1
    
    # Create DataFrame with the specified columns
    df_data = []
    for i in range(len(h)):
        row = {
            # Cat 1: Configuration parameters (same for each row)
            'level_name': config.level_name,
            'bench_idx': bench_idx,
            'elemsize': config.elemsize,
            'chunk_elems': config.chunk_elems,
            'cacheline_elems': config.cacheline_elems,
            'buffer_size': int(buffer_size),
            'sub_buffer_size': int(sub_buffer_size),
            'threads': threads,
            'interleave_flag': config.numa_level_reached if config.numa_level_reached > 1 else 0,
            # Cat 2: Measurement data (different per row)
            'h': h[i], # PANASTAS NOTE: This depends on how you translate gi - is it per thread or total)
            'av_timers': av_timers[i],
            'min_timers': min_timers[i],
            'max_timers': max_timers[i]
        }
        df_data.append(row)
    
    # Save to CSV
    df = pd.DataFrame(df_data)
    os.makedirs(os.path.dirname(filename), exist_ok=True)
    df.to_csv(filename, index=False)
    print(f"Results saved to: {filename}")
    log_info(f"Saved CSV: {filename}")
    return process_benchmark_results(h, av_timers, min_timers, max_timers,
                                   bench_idx, config)

def benchmark_setup(cur_level, prev_level, config):
    """
    Setup benchmark parameters for a given memory level.
    
    Args:
        cur_level: Current memory level data
        prev_level: Previous memory level data
        config: BenchmarkConfig object containing benchmark parameters
        
    Returns:
        tuple: (pi_mult, mem_sub, buffer_size, sub_buffer_size)
    """
    memi = cur_level['memi']
    pi = cur_level['pi']
    level_ctr = cur_level['level']
    level_name = cur_level['level_name']
    num_threads = config.num_threads

    prev_pi_mult = int(prev_level['pi_mult'])
    prev_mi = int(prev_level['mi']) 

    pi_mult = pi * prev_pi_mult
    bench_pi = 1 if prev_pi_mult > num_threads else min(num_threads, pi_mult)//prev_pi_mult
    if pi == 1: 
        print("Level %d(%s): pi = 1, quasi-level with pi_mult = %d" % (level_ctr, level_name, pi_mult))

    # For last level (main memory), do not use the actual sum of NUMA-node memories bellow.
    # Use the sum of LLC memory sizes instead (NUMA communication cost will be estimated via interleaving)
    if config.numa_level_reached:
        mem_sub = bench_pi * int(prev_level['mem_sub'])
        config.numa_level_reached *= pi
    else:
        mem_sub = bench_pi*(max(int(prev_level['mem_sub']), prev_mi))

    if (level_name == "NUMANode"):
        config.numa_level_reached = 1

    mi = memi - mem_sub

    if mi <= 0 :
        print("\nERROR: Level %d memory (%s) smaller than sum of previous memories (%s)" 
                % (level_ctr, format_bytes(memi*1024), format_bytes(mem_sub*1024)))
        exit(1)
    elif mi/2 <= mem_sub:
        print("\nWARNING: Level %d memory (%s) close to sum of previous memories (%s)" 
                % (level_ctr, format_bytes(memi*1024), format_bytes(mem_sub*1024)))

    # Determine bench_mem_multiplier based on allocation policy
    if config.alloc_policy == 'spread':
        bench_mem_multiplier = min(num_threads, config.system_cores // pi_mult)
    elif config.alloc_policy == 'close':
        # Default/'close' allocation
        bench_mem_multiplier = (num_threads // pi_mult) if (num_threads // pi_mult) else 1
    else:
        print(f"\nERROR: Unsupported allocation policy '{config.alloc_policy}'")
        print(f"Supported policies: 'close', 'spread'")
        exit(1)

    sub_buffer_size = mem_sub * 1024 * bench_mem_multiplier
    buffer_size = memi * 1024 * bench_mem_multiplier
    if not sub_buffer_size:
       # For lowest level cache hierarchy, increase the buffer to try to find BW
       buffer_size *= 2
    
    # Ensure buffer size is a multiple of cacheline size
    CACHELINE_SIZE = 64  # 64 bytes cache line size
    if buffer_size % CACHELINE_SIZE != 0:
        buffer_size = (buffer_size // CACHELINE_SIZE) * CACHELINE_SIZE
    if sub_buffer_size % CACHELINE_SIZE != 0:
        sub_buffer_size = (sub_buffer_size // CACHELINE_SIZE) * CACHELINE_SIZE

    if buffer_size > MAX_BUFF_SIZE:
        print("Warning: Buffer size %s exceeds MAX_BUFF_SIZE (%s), using MAX_BUFF_SIZE instead." 
            % (format_bytes(buffer_size), format_bytes(MAX_BUFF_SIZE)))
        buffer_size = MAX_BUFF_SIZE
        if sub_buffer_size > MAX_BUFF_SIZE/2:
            sub_buffer_size = MAX_BUFF_SIZE/2

    print(f"Running level {level_ctr} ({level_name}) benchmarks.\n"
            f"\tpi={pi}, pi_mult={pi_mult}, memi={format_bytes(memi*1024)}, mem_sub={format_bytes(mem_sub*1024)}\n"
            f"\tthreads={num_threads}, bench_mem_multiplier={bench_mem_multiplier}, bench_pi={bench_pi}, buffer_size={format_bytes(buffer_size)}, sub_buffer_size={format_bytes(sub_buffer_size)}")

    return pi_mult, mem_sub, buffer_size, sub_buffer_size

def create_hw_tree(config):
    """
    Create hardware tree by running benchmarks for all memory levels.
    
    Args:
        config: BenchmarkConfig object containing benchmark parameters
        
    Returns:
        DataFrame: Combined benchmark results for all levels
    """
    from hw_characterization.scripts.parser import parse_pahis_log_to_df, get_topology_file_path
    
    # Use hardcoded paths + configurable BUILD_DIR
    results_dir = RESULTS_DIR
    build_dir = BUILD_DIR
    
    log_info(f"Start analysis: system={config.system_name}, threads={config.num_threads}, k={config.consecutive_bytes}B")
    log_info(f"Dirs: topo={TOPOLOGY_DIR}, results={results_dir}, build={build_dir}")
    global CPU_ALLOWED_FROM_TOPO
    hw_topo, config.system_cores, config.thread_offset, CPU_ALLOWED_FROM_TOPO, config.r_scalar, config.r_vec = parse_pahis_log_to_df(
        get_topology_file_path(config.system_name))
    
    # Update config with thread_offset from topology
    print(hw_topo.to_markdown())
    print(f"system_cores={config.system_cores}, target_threads={config.num_threads}")
    log_info(f"Parsed topology; system_cores={config.system_cores}")
    output_df = pd.DataFrame()
    config.numa_level_reached = 0
    for level_ctr in range(len(hw_topo)):
        pi = int(hw_topo.iloc[level_ctr]['pi'])
        memi = int(hw_topo.iloc[level_ctr]['mi'])
        hw_topo.loc[hw_topo['level'] == level_ctr, 'memi'] = memi
        level_name = str(hw_topo.iloc[level_ctr]['level_name'])
        if level_ctr == 0:
            hw_topo.loc[hw_topo['level'] == level_ctr, 'pi_mult'] = int(hw_topo.iloc[level_ctr]['pi'])
            hw_topo.loc[hw_topo['level'] == level_ctr, 'mem_sub'] = 0
            continue

        curr_level = hw_topo.loc[hw_topo.index[hw_topo['level'] == level_ctr][0]]
        prev_level = hw_topo.loc[hw_topo.index[hw_topo['level'] == level_ctr-1][0]]
        pi_mult, mem_sub, buffer_size, sub_buffer_size = benchmark_setup(
            curr_level, prev_level, config)
        hw_topo.loc[hw_topo.index[hw_topo['level'] == level_ctr], 'pi_mult'] = pi_mult
        hw_topo.loc[hw_topo.index[hw_topo['level'] == level_ctr], 'mem_sub'] = mem_sub

        gi = []
        Li = []
        #Run chase benchmark to get a good latency estimation
        level_config = config.copy_with(level_name=level_name, forced_intercept=0)
        g_est, l_est = PaHiS_benchmark_run(0, buffer_size, sub_buffer_size, level_config)
        x_values = np.linspace(mem_sub+1, memi, 10000)
        x_range = np.array(x_values)
        predicted_times = g_est * x_range + l_est
        mean_time = np.mean(predicted_times/x_range)
        Li = mean_time * config.cacheline_size
        print(f"Access latency (Li) for level {level_ctr}: {Li:e} seconds")
        log_info(f"Li (access) for level {level_ctr} = {Li:e} s")

        idx_row = hw_topo.index[hw_topo['level'] == level_ctr][0]
        hw_topo.at[idx_row, 'Li'] = Li
        row = hw_topo.loc[idx_row]
        for idx in config.bw_kernel_indices:
            kernel_config = config.copy_with(level_name=level_name, forced_intercept=0)
            g, l = PaHiS_benchmark_run(idx, buffer_size, sub_buffer_size, kernel_config)
            if np.isnan(g) or np.isnan(l):
                print(f"(g, l) = ({g}, {l}), Retrying with forced intercept {Li/config.cacheline_size}")
                retry_config = config.copy_with(level_name=level_name, forced_intercept=Li/config.cacheline_size)
                g, l = PaHiS_benchmark_run(idx, buffer_size, sub_buffer_size, retry_config)
            if np.isnan(g) or np.isnan(l):
                hw_topo.loc[hw_topo['level'] == level_ctr, 'memi']*= 2
                print(f"(g, l) = ({g}, {l}), Retrying with double max buffer size {hw_topo.loc[hw_topo['level'] == level_ctr, 'memi']}")
                final_retry_config = config.copy_with(level_name=level_name, forced_intercept=Li/config.cacheline_size)
                g, l = PaHiS_benchmark_run(idx, buffer_size, sub_buffer_size, final_retry_config)
            
            record = row.to_dict()
            record.update({'consecutive_bytes': int(config.consecutive_bytes), 'threads': int(config.num_threads), 'interleave_flag': config.numa_level_reached if config.numa_level_reached > 1 else 0, 'buffer_size': int(buffer_size),
             'sub_buffer_size': int(sub_buffer_size),'idx': int(idx), 'gi': float(g)})
            output_df = pd.concat([output_df, pd.DataFrame([record])], ignore_index=True)
            gi.append(g)
        gi = aggregate_values(gi, config.kernel_aggregator)
        hw_topo.at[idx_row, 'gi'] = gi

    # Add GLOBAL_SYNC level if requested
    if config.create_sync_level:
        sync_config = config.copy_with(
            level_name='GLOBAL_SYNC', forced_intercept=-42)
        _, Li_sync = PaHiS_benchmark_run(1, 0, 0, sync_config)
        print(f"SYNC latency (Li_sync): {Li_sync:e} seconds")
        log_info(f"Li_sync = {Li_sync:e} s")
        
        # Get last level's values
        last_row = hw_topo.iloc[-1]
        sync_level = {
            'level': int(last_row['level']) + 1,
            'level_name': "GLOBAL_SYNC",
            'pi': 1,
            'gi': 0.0,  # GLOBAL_SYNC has no bandwidth (synchronization only)
            'Li': Li_sync,
            'mi': last_row['mi'],
            'memi': last_row['memi'],
            'pi_mult': last_row['pi_mult'],
            'mem_sub': last_row['mem_sub'],
            'consecutive_bytes': config.consecutive_bytes,
            'threads': config.num_threads,
            'buffer_size': last_row['buffer_size'] if 'buffer_size' in last_row else 0,
            'sub_buffer_size': last_row['sub_buffer_size'] if 'sub_buffer_size' in last_row else 0,
            'idx': 0
        }
        # Add to hw_topo and output_df
        hw_topo = pd.concat([hw_topo, pd.DataFrame([sync_level])], ignore_index=True)
        output_df = pd.concat([output_df, pd.DataFrame([sync_level])], ignore_index=True)

    hw_topo = hw_topo[['level', 'level_name', 'pi', 'mi', 'gi', 'Li',
                           'pi_mult', 'memi', 'mem_sub']]
    print(output_df.to_markdown())
    return output_df


def _detect_violators(agg, level_names):
    """
    Detect violator positions in an aggregated value sequence that should be
    monotonically non-decreasing.

    Rules:
    - Edge levels (first/last) are always eligible for flagging.
    - Middle levels are only eligible if their two neighbors are compliant with
      each other (left <= right). Otherwise the middle level is skipped because
      the real problem may lie in a neighbor.
    - Deadlock fallback: if no violators are flagged but violations still exist,
      flag the single level with the largest deviation from avg(left, right).

    Args:
        agg: list of float — aggregated values per level position (GLOBAL_SYNC excluded)
        level_names: list of str — level names corresponding to each position

    Returns:
        list of int — positions of flagged violators (may be empty if monotonic)
    """
    n = len(agg)
    if n <= 1:
        return []

    violators = []

    for i in range(n):
        is_first = (i == 0)
        is_last = (i == n - 1)

        if is_first:
            # Bump at first level: higher than its right neighbor
            if agg[i] > agg[i + 1]:
                violators.append(i)
        elif is_last:
            # Dip at last level: lower than its left neighbor
            if agg[i] < agg[i - 1]:
                violators.append(i)
        else:
            # Middle level: only eligible if neighbors are compliant with each other
            left, right = agg[i - 1], agg[i + 1]
            if left <= right:
                # Bump: higher than both neighbors
                if agg[i] > right and agg[i] > left:
                    violators.append(i)
                # Dip: lower than both neighbors
                elif agg[i] < left and agg[i] < right:
                    violators.append(i)

    # Deadlock fallback: violations exist but no violator was flagged
    if not violators:
        has_violation = any(agg[i] > agg[i + 1] for i in range(n - 1))
        if has_violation:
            # Flag the single level with the largest deviation from avg(left, right)
            best_idx = -1
            best_dev = -1.0
            for i in range(n):
                if i == 0:
                    expected = 0.75 * agg[1]
                elif i == n - 1:
                    expected = 1.5 * agg[n - 2]
                else:
                    expected = (agg[i - 1] + agg[i + 1]) / 2.0
                dev = abs(agg[i] - expected)
                if dev > best_dev:
                    best_dev = dev
                    best_idx = i
            if best_idx >= 0:
                violators.append(best_idx)

    return violators


def _compute_target(agg, i):
    """
    Compute the interpolated target value for a violator at position i.

    - First level:  0.75 * right neighbor
    - Last level:   1.5  * left neighbor
    - Middle level: average of left and right neighbors
    """
    n = len(agg)
    if i == 0:
        return 0.75 * agg[1]
    elif i == n - 1:
        return 1.5 * agg[n - 2]
    else:
        return (agg[i - 1] + agg[i + 1]) / 2.0


def _enforce_monotonicity_column(df_adjusted, group_mask, levels, level_names_map,
                                  column, kernel_aggregator, max_iterations=10):
    """
    Iteratively enforce monotonicity on a single column ('gi' or 'Li') for one
    (threads, alloc_policy, consecutive_bytes) group using violator detection,
    neighbor interpolation, and proportional scaling.

    Args:
        df_adjusted: the mutable DataFrame being adjusted in-place
        group_mask: boolean Series selecting this group's rows
        levels: sorted list of level numbers (GLOBAL_SYNC excluded)
        level_names_map: dict mapping level number to level name
        column: 'gi' or 'Li'
        kernel_aggregator: 'mean', 'min', or 'max'
        max_iterations: safety limit on convergence iterations

    Returns:
        int — total number of individual row adjustments made
    """
    total_adjustments = 0

    # Filter out GLOBAL_SYNC levels from consideration
    active_levels = [l for l in levels if level_names_map.get(l, 'Unknown') != "GLOBAL_SYNC"]
    n = len(active_levels)
    if n <= 1:
        return 0

    for iteration in range(1, max_iterations + 1):
        # Step 1: Compute aggregated values per level
        agg = []
        for level in active_levels:
            mask = group_mask & (df_adjusted['level'] == level)
            agg.append(aggregate_values(df_adjusted.loc[mask, column], kernel_aggregator))

        # Check if already monotonic
        is_monotonic = all(agg[i] <= agg[i + 1] for i in range(n - 1))
        if is_monotonic:
            if iteration > 1:
                logger.debug(f"      {column} converged after {iteration - 1} iteration(s)")
            break

        level_names_list = [level_names_map.get(l, 'Unknown') for l in active_levels]

        # Step 2: Detect violators
        violators = _detect_violators(agg, level_names_list)

        if not violators:
            # Should not happen (not monotonic but no violators), but guard anyway
            logger.debug(f"      {column} iteration {iteration}: non-monotonic but no violators detected, stopping")
            break

        logger.debug(f"      {column} iteration {iteration}: agg={[f'{v:.4e}' for v in agg]}")
        logger.debug(f"      {column} iteration {iteration}: violators at positions {violators} "
                      f"({[level_names_list[v] for v in violators]})")

        # Steps 3 & 4: Compute targets and apply proportional scaling
        for v in violators:
            level = active_levels[v]
            level_name = level_names_map.get(level, 'Unknown')
            mask = group_mask & (df_adjusted['level'] == level)
            current_agg = agg[v]

            target = _compute_target(agg, v)

            if current_agg == 0:
                # Guard against division by zero: set all rows to target directly
                logger.debug(f"        {level_name}: agg=0, setting all rows to {target:.6e}")
                df_adjusted.loc[mask, column] = target
            else:
                scale = target / current_agg
                df_adjusted.loc[mask, column] *= scale
                logger.debug(f"        {level_name}: {current_agg:.6e} -> {target:.6e} "
                              f"(scale={scale:.4f}, {mask.sum()} rows)")

            total_adjustments += mask.sum()
    else:
        # max_iterations reached without convergence
        agg_final = []
        for level in active_levels:
            mask = group_mask & (df_adjusted['level'] == level)
            agg_final.append(aggregate_values(df_adjusted.loc[mask, column], kernel_aggregator))
        logger.warning(f"      {column} did not converge after {max_iterations} iterations. "
                        f"Final agg={[f'{v:.4e}' for v in agg_final]}")

    return total_adjustments


def enforce_parameter_monotonicity(df, kernel_aggregator='mean', keep_levels=None):
    """
    Enforce monotonicity constraints on bandwidth (gi) and latency (Li) values in the dataframe.

    Uses iterative violator detection: levels whose aggregated value breaks the
    expected non-decreasing pattern are identified and adjusted via neighbor
    interpolation with proportional scaling (preserving per-kernel variability).

    Violator detection rules:
    - Edge levels (first/last) are always eligible for flagging.
    - Middle levels are only eligible if their two neighbors are compliant with
      each other (left <= right).
    - Deadlock fallback: when no violator is flagged but violations persist,
      the level with the largest deviation from its expected position is flagged.

    Args:
        df: DataFrame with benchmark results containing columns: 'level', 'level_name', 'gi', 'Li',
            'threads', 'alloc_policy', 'consecutive_bytes', etc.
        kernel_aggregator: Method to aggregate across kernel rows ('mean', 'min', 'max')
        keep_levels: Optional list of level names to restrict enforcement to (None = all levels).
                     When provided, only these levels participate in monotonicity checks,
                     preventing excluded levels (e.g. L1Cache) from influencing kept ones.

    Returns:
        DataFrame: DataFrame with adjusted gi and Li values that satisfy monotonicity
    """
    logger.info(f"Enforcing parameter monotonicity (aggregator={kernel_aggregator}) on raw benchmark data...")
    if keep_levels:
        logger.info(f"  Restricting to keep_levels: {keep_levels}")
    logger.debug(f"  Input dataframe shape: {df.shape}")

    # Create a copy to avoid modifying the original
    df_adjusted = df.copy()

    # Track how many values were adjusted
    total_gi_adjustments = 0
    total_li_adjustments = 0

    # Group by configuration (threads, alloc_policy, consecutive_bytes) to process each separately
    groupby_cols = ['threads', 'alloc_policy', 'consecutive_bytes']

    for group_key, group_df in df.groupby(groupby_cols):
        threads, alloc_policy, consecutive_bytes = group_key
        logger.debug(f"  Processing: threads={threads}, policy={alloc_policy}, k={consecutive_bytes}B")

        # Get mask for this group in the adjusted dataframe
        group_mask = ((df_adjusted['threads'] == threads) &
                     (df_adjusted['alloc_policy'] == alloc_policy) &
                     (df_adjusted['consecutive_bytes'] == consecutive_bytes))

        # Get current state from adjusted dataframe
        group_df_current = df_adjusted[group_mask].sort_values('level').copy()

        # Get unique levels in order (excluding level 0 which doesn't have measurements)
        levels = [l for l in sorted(group_df_current['level'].unique()) if l > 0]
        if not levels:
            continue

        # Build level name mapping (get first occurrence of each level)
        level_names_map = {}
        for level in levels:
            level_data = group_df_current[group_df_current['level'] == level]
            if not level_data.empty:
                level_names_map[level] = level_data.iloc[0]['level_name']

        # Filter to keep_levels if specified (only enforce monotonicity on these levels)
        if keep_levels is not None:
            levels = [l for l in levels if level_names_map.get(l, 'Unknown') in keep_levels]
            if not levels:
                continue

        # Enforce monotonicity for gi (bandwidth)
        logger.debug(f"    Enforcing gi monotonicity across {len(levels)} levels: {[level_names_map.get(l, '?') for l in levels]}")
        total_gi_adjustments += _enforce_monotonicity_column(
            df_adjusted, group_mask, levels, level_names_map,
            'gi', kernel_aggregator)

        # Enforce monotonicity for Li (latency)
        logger.debug(f"    Enforcing Li monotonicity across {len(levels)} levels")
        total_li_adjustments += _enforce_monotonicity_column(
            df_adjusted, group_mask, levels, level_names_map,
            'Li', kernel_aggregator)

    if total_gi_adjustments > 0 or total_li_adjustments > 0:
        logger.info(f"  Monotonicity adjustments: {total_gi_adjustments} gi rows, {total_li_adjustments} Li rows adjusted")
    else:
        logger.info(f"  No monotonicity adjustments needed")

    logger.debug(f"  Output dataframe shape: {df_adjusted.shape}")
    return df_adjusted


def generate_all_subsets(level_names):
    """
    Generate all possible subsets of level names that include NodeMem (excluding empty set).
    
    Args:
        level_names: List of level names
        
    Returns:
        list: List of all possible subsets as lists of level names (all containing NodeMem)
    """
    from itertools import combinations
    
    all_subsets = []
    n = len(level_names)
    
    # Generate all subsets from size 1 to n, but only keep those that include NodeMem
    for r in range(1, n + 1):
        for subset in combinations(level_names, r):
            subset_list = list(subset)
            # Only include subsets that contain NodeMem
            if 'NodeMem' in subset_list:
                all_subsets.append(subset_list)
    
    return all_subsets


def _interpolate_monotonic_ls(thread_count, current_policy, fallback_ls):
    """
    Resolve GLOBAL_SYNC latency for empirical_monotony mode using the
    precomputed monotonic dataset.
    """
    adj_mon_key = f'adjusted_monotonic_{current_policy}'
    adj_data = globals().get(adj_mon_key)
    if not adj_data:
        logger.debug(f"    Warning: Adjusted monotonic data not found for {current_policy}, using empirical value")
        return fallback_ls

    adj_threads = adj_data.get('threads', [])
    adj_ls = adj_data.get('ls', [])
    if not adj_threads or not adj_ls:
        logger.debug(f"    Warning: Adjusted monotonic data empty for {current_policy}, using empirical value")
        return fallback_ls

    if thread_count in adj_threads:
        idx = adj_threads.index(thread_count)
        ls_value = adj_ls[idx]
        logger.debug(f"    Using adjusted monotonic data for GLOBAL_SYNC: threads={thread_count}, ls={ls_value:.6e} (exact match)")
    else:
        adj_threads_arr = np.array(adj_threads)
        adj_ls_arr = np.array(adj_ls)
        min_threads = adj_threads_arr.min()
        max_threads = adj_threads_arr.max()

        if thread_count < min_threads:
            ls_value = adj_ls_arr[np.argmin(adj_threads_arr)]
            logger.debug(f"    Using adjusted monotonic data for GLOBAL_SYNC: threads={thread_count} < min({min_threads}), using min ls={ls_value:.6e}")
        elif thread_count > max_threads:
            ls_value = adj_ls_arr[np.argmax(adj_threads_arr)]
            logger.debug(f"    Using adjusted monotonic data for GLOBAL_SYNC: threads={thread_count} > max({max_threads}), using max ls={ls_value:.6e}")
        else:
            ls_value = float(np.interp(thread_count, adj_threads_arr, adj_ls_arr))
            logger.debug(f"    Using adjusted monotonic data for GLOBAL_SYNC: threads={thread_count}, ls={ls_value:.6e} (interpolated from {adj_threads})")

    if ls_value <= 0:
        logger.debug(f"    ERROR: Adjusted monotonic ls is non-positive ({ls_value:.6e})")
        exit(1)

    return ls_value


def _predict_polynomial_ls(thread_count, current_policy, fallback_ls):
    """
    Resolve GLOBAL_SYNC latency for empirical_poly mode using stored
    polynomial coefficients and enforce boundary equality across policies.
    """
    lr_key = f'lr_coeffs_{current_policy}'
    lr_coeffs = globals().get(lr_key)
    if not lr_coeffs:
        logger.debug(f"    Warning: LR coefficients not found for {current_policy}, using empirical value")
        return fallback_ls

    a = lr_coeffs.get('a', 0.0)
    b = lr_coeffs.get('b', 0.0)
    c = lr_coeffs.get('c', 0.0)
    model_name = lr_coeffs.get('model_name', 'Unknown')

    predicted_ls = a * (thread_count ** 2) + b * thread_count + c
    if predicted_ls <= 0:
        logger.debug(f"    ERROR: Predicted Li_sync is non-positive ({predicted_ls:.6e})")
        logger.debug(f"    This violates the constraint that latency must be positive.")
        exit(1)

    boundaries = globals().get('global_sync_boundaries', {})
    boundary_ls = boundaries.get('ls', {}) if isinstance(boundaries, dict) else {}
    min_threads = boundaries.get('min_threads') if isinstance(boundaries, dict) else None
    max_threads = boundaries.get('max_threads') if isinstance(boundaries, dict) else None

    if boundary_ls:
        if min_threads is not None and thread_count == min_threads and min_threads in boundary_ls:
            predicted_ls = boundary_ls[min_threads]
        elif max_threads is not None and thread_count == max_threads and max_threads in boundary_ls:
            predicted_ls = boundary_ls[max_threads]

    logger.debug(f"    Using {model_name} prediction for GLOBAL_SYNC: Li_sync = {predicted_ls:.6e}")
    return predicted_ls


def compute_global_sync_latency(thread_count, level_name, current_policy, agr_ls, create_sync_level):
    """
    Central helper for determining GLOBAL_SYNC latency adjustments.
    """
    if level_name != "GLOBAL_SYNC":
        return agr_ls

    if create_sync_level == "empirical_monotony":
        return _interpolate_monotonic_ls(thread_count, current_policy, agr_ls)

    if create_sync_level == "empirical_poly":
        return _predict_polynomial_ls(thread_count, current_policy, agr_ls)

    return agr_ls


def enforce_polynomial_boundary_alignment(a, b, c, min_threads, max_threads, boundary_map, policy_name):
    """
    Adjust polynomial coefficients so that p(min_threads) and p(max_threads)
    match the cross-policy boundary averages.
    """
    if not boundary_map:
        return a, b, c, False

    if min_threads not in boundary_map or max_threads not in boundary_map:
        return a, b, c, False

    if max_threads == min_threads:
        return a, b, c, False

    desired_min = boundary_map[min_threads]
    desired_max = boundary_map[max_threads]

    denom = max_threads - min_threads
    adjusted_b = ((desired_max - desired_min) - a * (max_threads ** 2 - min_threads ** 2)) / denom
    adjusted_c = desired_min - a * (min_threads ** 2) - adjusted_b * min_threads

    logger.debug(
        f"    Applied boundary alignment for {policy_name}: "
        f"b {b:.6e}->{adjusted_b:.6e}, c {c:.6e}->{adjusted_c:.6e} "
        f"(targets: ls({min_threads})={desired_min:.6e}, ls({max_threads})={desired_max:.6e})"
    )
    return a, adjusted_b, adjusted_c, True


def adjust_group_with_ramp(
    group_points,
    prev_max_ls=None,
    monotony_factor=1.2,
    boundary_averages=None,
    min_threads=None,
    max_threads=None,
):
    """
    Given a list of points that violate weighted monotonicity, build an
    increasing ramp from min(ls) to max(ls) with the specified step size.
    The ramp is shifted if needed to satisfy the previous maximum value
    scaled by the monotony factor.
    """
    if not group_points:
        return []

    sorted_points = sorted(group_points, key=lambda g: g['threads'])
    ls_values = [p['ls_avg'] for p in sorted_points]
    contains_min = min_threads is not None and any(p['threads'] == min_threads for p in sorted_points)
    contains_max = max_threads is not None and any(p['threads'] == max_threads for p in sorted_points)
    boundary_min_value = boundary_averages.get(min_threads) if boundary_averages and contains_min and min_threads in boundary_averages else None
    boundary_max_value = boundary_averages.get(max_threads) if boundary_averages and contains_max and max_threads in boundary_averages else None

    min_candidates = ls_values.copy()
    max_candidates = ls_values.copy()
    if boundary_min_value is not None:
        min_candidates.append(boundary_min_value)
    if boundary_max_value is not None:
        max_candidates.append(boundary_max_value)

    min_ls = min(min_candidates) if min_candidates else 0.0
    max_ls = max(max_candidates) if max_candidates else min_ls
    span = max_ls - min_ls
    count = len(sorted_points)

    # Ensure the starting point satisfies the previous requirement
    min_target = min_ls
    if boundary_min_value is not None:
        min_target = boundary_min_value
    if prev_max_ls is not None:
        min_target = max(min_target, prev_max_ls * monotony_factor)
    max_target = max_ls
    if boundary_max_value is not None:
        max_target = boundary_max_value
    if max_target < min_target:
        max_target = min_target

    if count == 1:
        return [{
            'threads': sorted_points[0]['threads'],
            'ls_avg': min_target,
            'count': sorted_points[0].get('count', 1)
        }]

    # Step defined as (max_ls - min_ls) / group_member_num per specification
    step = span / count if count > 0 else 0.0
    ramp_values = []
    for idx in range(count):
        base_value = min_target + step * idx
        if idx == count - 1:
            base_value = max_target
        if idx > 0:
            base_value = max(base_value, ramp_values[-1] * monotony_factor)
        ramp_values.append(base_value)

    adjusted_points = []
    for point, value in zip(sorted_points, ramp_values):
        if boundary_min_value is not None and point['threads'] == min_threads:
            value = boundary_min_value
        if boundary_max_value is not None and point['threads'] == max_threads:
            value = boundary_max_value
        if adjusted_points:
            value = max(value, adjusted_points[-1]['ls_avg'] * monotony_factor)
        adjusted_points.append({
            'threads': point['threads'],
            'ls_avg': value,
            'count': point.get('count', 1)
        })
    return adjusted_points


def finalize_weighted_group(
    current_group,
    merged_data,
    monotony_factor=1.2,
    context_label="",
    boundary_averages=None,
    min_threads=None,
    max_threads=None,
):
    """
    Finalize a group of points by enforcing the weighted monotonicity rule.
    Returns True if any adjustments were applied.
    """
    if not current_group:
        return False

    adjusted = False
    group_threads = [g['threads'] for g in current_group]
    min_threads_in_group = min(group_threads)
    prev_candidates = [m['ls_avg'] for m in merged_data if m['threads'] < min_threads_in_group]
    prev_max_ls = max(prev_candidates) if prev_candidates else None

    if len(current_group) == 1:
        finalized_ls = current_group[0]['ls_avg']
        finalized_threads = current_group[0]['threads']
        required_min = prev_max_ls * monotony_factor if prev_max_ls is not None else None
        boundary_override = None
        if boundary_averages:
            if finalized_threads == min_threads and min_threads in boundary_averages:
                boundary_override = boundary_averages[min_threads]
            elif finalized_threads == max_threads and max_threads in boundary_averages:
                boundary_override = boundary_averages[max_threads]

        if boundary_override is not None:
            logger.debug(f"{context_label}Using boundary average for threads={finalized_threads}: {boundary_override:.6e}")
            finalized_ls = boundary_override

        if required_min is not None and finalized_ls < required_min:
            adjusted = True
            logger.debug(f"{context_label}⚠ Adjusted single-point group: threads={finalized_threads}, ls={finalized_ls:.6e} -> {required_min:.6e} (weighted monotonicity)")
            finalized_ls = required_min
        else:
            logger.debug(f"{context_label}✓ No violation, finalized single-point group: threads={finalized_threads}, ls={finalized_ls:.6e}")

        merged_data.append({
            'threads': finalized_threads,
            'ls_avg': finalized_ls,
            'count': int(current_group[0].get('count', 1))
        })
    else:
        adjusted_points = adjust_group_with_ramp(
            current_group,
            prev_max_ls,
            monotony_factor,
            boundary_averages=boundary_averages,
            min_threads=min_threads,
            max_threads=max_threads,
        )
        original_ls = [f"{p['ls_avg']:.6e}" for p in current_group]
        adjusted_ls = [f"{p['ls_avg']:.6e}" for p in adjusted_points]
        logger.debug(f"{context_label}⚠ Weighted adjustment for merged group:")
        logger.debug(f"{context_label}  Threads: {[p['threads'] for p in current_group]}")
        logger.debug(f"{context_label}  Original ls: {original_ls}")
        logger.debug(f"{context_label}  Adjusted ramp: {adjusted_ls}")
        merged_data.extend(adjusted_points)
        adjusted = True

    return adjusted
def parse_hw_parameters(system_name, df, threads_list=None, keep_levels=None, aggregate_method='mean', write_cpp_header=True, alloc_policy=None, all_alloc_policies=None, generate_subsets=False, create_sync_level="", all_combined_df=None, r_scalar=0.0, r_vec=0.0):
    """
    Parse combined benchmark dataframe and output hardware parameters for predict_k_PaHiS_kernel.
    Generates a system-specific header file with ALL possible subset combinations of levels.
    
    Args:
        system_name: name of the system
        df: DataFrame with benchmark results
        threads_list: List of thread counts to include (None for all available threads)
        keep_levels: List of level names to keep (None to keep all levels) - used as base levels for subsets
        aggregate_method: Method to aggregate values ('mean', 'min', 'max')
        write_cpp_header: Whether to write parameters to C++ header file
        alloc_policy: Current NUMA allocation policy being processed (e.g., "close", "spread")
        all_alloc_policies: List of all allocation policies being processed (e.g., ["close", "spread"])
        generate_subsets: If True, generate all subset combinations of levels (always including NodeMem). If False, only use the exact levels given (default: False)
    
    Returns:
        dict: Dictionary of hardware parameters per thread configuration
    """
    # If threads_list is None, use all unique thread counts in the dataset
    if threads_list is None:
        threads_list = sorted(df['threads'].unique())
    
    # Get system name from the first entry
    system_id = system_name.replace('-', '_').replace(' ', '_').lower()
    
    # Get all available levels from the dataset
    all_available_levels = sorted(df['level'].unique())
    
    # Create level name mapping
    level_names_mapping = {}
    for level in all_available_levels:
        if level == 0:  # Skip level 0 as it doesn't have measurement data
            continue
        level_data = df[df['level'] == level]
        if not level_data.empty:
            level_name = level_data['level_name'].iloc[0]
            level_names_mapping[level_name] = level
    
    # Determine base levels for subset generation
    if keep_levels is not None:
        base_levels = [name for name in keep_levels if name in level_names_mapping]
    else:
        base_levels = list(level_names_mapping.keys())
    
    # Always require NodeMem to be in the base levels
    if 'NodeMem' in level_names_mapping and 'NodeMem' not in base_levels:
        raise ValueError("NodeMem level must always be included in keep_levels. Please add 'NodeMem' to your --keep-levels argument.")
    
    print(f"Base levels for subset generation: {base_levels}")
    
    # Generate level subsets based on generate_subsets flag
    if not generate_subsets:
        # Only use the exact levels specified, no subsets
        all_level_subsets = [base_levels]
        print(f"No subsets mode: using only the specific levels: {base_levels}")
    else:
        # Generate all possible subsets of base levels
        all_level_subsets = generate_all_subsets(base_levels)
        logger.debug(f"Generated {len(all_level_subsets)} level subset combinations:")
        for i, subset in enumerate(all_level_subsets):
            logger.debug(f"  {i+1}: {subset}")

    # Create output directory if it doesn't exist
    if write_cpp_header:
        cpp_base_dir = CPP_MODEL_DIR
        os.makedirs(cpp_base_dir, exist_ok=True)
        # Simplified filename - just use system name since all metadata is now in the structure
        header_file = f"{cpp_base_dir}/hw_params_{system_id}.hpp"
        logger.info(f"Header file generation enabled: {header_file}")
        
        # Determine policies for header
        if all_alloc_policies is not None:
            policies_str = ', '.join(all_alloc_policies)
        else:
            policies_str = alloc_policy
            
        # Start writing the header file with includes and defines
        subset_mode_str = f"ALL subset combinations ({len(all_level_subsets)} total)" if generate_subsets else "Exact levels only (no subsets)"
        
        # Reconstruct the command line from sys.argv
        # Properly quote arguments that contain spaces or special characters
        quoted_args = [shlex.quote(arg) for arg in sys.argv[1:]]
        command_line = "python3 runner.py " + " ".join(quoted_args)
        
        header_content = f"""
// Auto-generated hardware parameters for {system_name}
// Allocation policies: {policies_str}
// Base levels: {', '.join(base_levels) if base_levels else 'all'}
// {subset_mode_str}
// Generated on {datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")}
//
// Command used to generate this file:
// {command_line}

#ifndef HW_PARAMS_{system_id.upper()}_HPP
#define HW_PARAMS_{system_id.upper()}_HPP

#include "cost_models.hpp"

// Array of hardware parameters for different thread configurations and level subsets
const std::vector<cost_models::HW_model::HWParameters> hw_models_vector = {{

"""
    else:
        logger.warning("Header file generation is DISABLED (write_cpp_header=False)")
    
    # Process each allocation policy, thread configuration, and level subset combination
    hw_params_dict = {}
    all_hw_models = []
    all_thread_ids = []
    all_policy_ids = []
    all_level_bitmasks = []
    
    model_index = 0
    
    # Determine which policies to process
    if all_alloc_policies is not None:
        policies_to_process = all_alloc_policies
    else:
        policies_to_process = [alloc_policy]
    
    for policy_idx, current_policy in enumerate(policies_to_process):
        logger.info(f"  Processing policy {policy_idx + 1}/{len(policies_to_process)}: {current_policy}")
        logger.debug(f"\nProcessing policy {policy_idx + 1}/{len(policies_to_process)}: {current_policy}")
        
        # Filter data for this policy
        policy_df = df[df.get('alloc_policy', current_policy) == current_policy] if 'alloc_policy' in df.columns else df
        
        # Show thread counts as a clean list
        clean_threads = [int(t) for t in threads_list]
        logger.info(f"    Threads: {clean_threads}")
        
        for thread_count in threads_list:
            # Filter data for this thread count and policy
            thread_df = policy_df[policy_df['threads'] == thread_count]
            
            if thread_df.empty:
                logger.info(f"      Warning: No data for thread count {thread_count} with policy {current_policy}, skipping")
                logger.debug(f"Warning: No data for thread count {thread_count} with policy {current_policy}, skipping")
                continue
            
            # Find max consecutive_bytes in the dataset
            max_consecutive_bytes = thread_df['consecutive_bytes'].max()
            
            logger.debug(f"\nProcessing {thread_count} thread(s) with policy {current_policy}:")
            
            # Process each level subset combination
            for subset_idx, level_subset in enumerate(all_level_subsets):
                logger.debug(f"  Subset {subset_idx + 1}/{len(all_level_subsets)}: {level_subset}")
                
                # Get level numbers for this subset
                subset_levels = [level_names_mapping[name] for name in level_subset]
                subset_levels.sort()
                
                # Collect data for levels in this subset
                subset_level_data = []
                for level in subset_levels:
                    level_data = thread_df[thread_df['level'] == level]
                    if not level_data.empty:
                        level_name = level_data['level_name'].iloc[0]
                        subset_level_data.append({
                            'level': level,
                            'level_name': level_name,
                            'data': level_data
                        })
                
                if not subset_level_data:
                    print(f"    Warning: No data for subset {level_subset}, skipping")
                    continue
                
                # Initialize parameter lists
                g_list = []
                ls_list = []
                m_list = []
                p_list = []
                
                # Calculate pi adjustments for skipped levels
                # We need to account for all pi values from skipped intermediate levels
                
                # Get all available levels from the full dataset for this thread count
                all_levels_df = thread_df.groupby('level').first().reset_index()
                all_levels_df = all_levels_df.sort_values('level')
                
                # Process levels in this subset
                for i, level_info in enumerate(subset_level_data):
                    level = level_info['level']
                    level_name = level_info['level_name']
                    level_data = level_info['data']
                    
                    # Get current level's pi
                    pi = int(level_data['pi'].iloc[0])
                    
                    # Calculate accumulated pi for skipped levels (multiplicative)
                    accumulated_pi = pi  # Start with current level's pi
                    
                    # Find the previous level in the subset
                    if i > 0:
                        prev_level_in_subset = subset_level_data[i-1]['level']
                        
                        # Multiply pi values from all levels between previous subset level and current level
                        for intermediate_level in range(prev_level_in_subset + 1, level):
                            intermediate_level_data = all_levels_df[all_levels_df['level'] == intermediate_level]
                            if not intermediate_level_data.empty:
                                intermediate_pi = int(intermediate_level_data['pi'].iloc[0])
                                accumulated_pi *= intermediate_pi
                                logger.debug(f"    Multiplied by pi={intermediate_pi} from skipped level {intermediate_level} ({intermediate_level_data['level_name'].iloc[0]})")
                    else:
                        # First level in subset - multiply pi values from all levels before it
                        for level_before in range(level):
                            level_before_data = all_levels_df[all_levels_df['level'] == level_before]
                            if not level_before_data.empty:
                                level_before_pi = int(level_before_data['pi'].iloc[0])
                                accumulated_pi *= level_before_pi
                                logger.debug(f"    Multiplied by pi={level_before_pi} from level {level_before} ({level_before_data['level_name'].iloc[0]}) before {level_name}")
                    
                    if accumulated_pi != pi:
                        logger.debug(f"    Accumulated pi for level {level} ({level_name}): {pi} -> {accumulated_pi}")
                    
                    # Get latency (ls) from consecutive_bytes=1
                    latency_data = level_data[level_data['consecutive_bytes'] == 1]
                    if latency_data.empty:
                        # If no data for consecutive_bytes=1, use smallest available
                        min_bytes = level_data['consecutive_bytes'].min()
                        latency_data = level_data[level_data['consecutive_bytes'] == min_bytes]
                        logger.debug(f"    Warning: No data for consecutive_bytes=1 at level {level}, using {min_bytes} instead")
                    
                    # Get bandwidth (g) from max consecutive_bytes
                    bandwidth_data = level_data[level_data['consecutive_bytes'] == max_consecutive_bytes]
                    if bandwidth_data.empty:
                        # If no data for max consecutive_bytes, use largest available
                        avail_bytes = level_data['consecutive_bytes'].max()
                        bandwidth_data = level_data[level_data['consecutive_bytes'] == avail_bytes]
                        logger.debug(f"    Warning: No data for consecutive_bytes={max_consecutive_bytes} at level {level}, using {avail_bytes} instead")
                    
                    # Calculate values based on aggregate_method
                    ls_values = latency_data['gi'].astype(float)
                    g_values = bandwidth_data['gi'].astype(float)
                    
                    agr_ls = aggregate_values(ls_values, aggregate_method)
                    agr_g = aggregate_values(g_values, aggregate_method)
                    
                    ls_value = compute_global_sync_latency(
                        thread_count=thread_count,
                        level_name=level_name,
                        current_policy=current_policy,
                        agr_ls=agr_ls,
                        create_sync_level=create_sync_level
                    )
                    ls_list.append(ls_value)
                    
                    # GLOBAL_SYNC has no bandwidth (synchronization only)
                    if level_name == "GLOBAL_SYNC":
                        g_list.append(0.0)
                        logger.debug(f"    GLOBAL_SYNC bandwidth set to 0.0 (synchronization only)")
                    else:
                        g_list.append(agr_g)
                    
                    # Get memory size (in bytes) - mi is KB, convert to bytes
                    mi = int(level_data['mi'].iloc[0]) * 1024
                    m_list.append(mi)
                    
                    # Use the accumulated pi value
                    p_list.append(accumulated_pi)
                
                # Calculate d_numa: the level (1-indexed) at which NUMA effects start (first level after LLC)
                # d_numa = 0 means no NUMA effects (all levels are cache or no levels after cache)
                d_numa = 0
                last_cache_idx = -1
                for idx, name in enumerate(level_subset):
                    if "Cache" in name:
                        last_cache_idx = idx
                
                if last_cache_idx >= 0 and last_cache_idx + 1 < len(level_subset):
                    # NUMA effects start at the first non-GLOBAL_SYNC level after the last cache
                    for idx in range(last_cache_idx + 1, len(level_subset)):
                        if level_subset[idx] != "GLOBAL_SYNC":
                            d_numa = idx + 1  # Convert to 1-indexed level number
                            break
                
                # Encode policy in d_numa: positive for close, negative for spread
                if current_policy == "spread":
                    d_numa = -d_numa  # Make negative for spread policy
                
                logger.debug(f"    Calculated d_numa={d_numa} (last cache at idx {last_cache_idx}, NUMA starts at level {abs(d_numa) if d_numa != 0 else 'N/A'}, policy={current_policy})")
                
                # Monotonicity for g and ls is already enforced in enforce_parameter_monotonicity()
                # before subset generation, so the values in df are already monotonic
                
                # Adjust p for this thread/policy configuration:
                # Instead of storing raw topology p (same for all configs), adjust so
                # p[l] = number of active sub-components at level l for this thread count
                # and policy, preserving hierarchical nesting.
                P_phys = []
                acc = 1
                for pi_val in p_list:
                    acc *= pi_val
                    P_phys.append(acc)
                total_cores = P_phys[-1]

                if current_policy == "close":
                    # Close: threads pack bottom-up
                    active = [min(thread_count, P_phys[l]) for l in range(len(p_list))]
                else:
                    # Spread: threads distributed evenly from top-down
                    active = [min(P_phys[l], math.ceil(thread_count * P_phys[l] / total_cores))
                              for l in range(len(p_list))]

                p_adj = [active[0]]
                for l in range(1, len(p_list)):
                    p_adj.append(min(p_list[l], max(1, math.ceil(active[l] / active[l - 1]))))
                
                logger.debug(f"    p adjustment: physical={p_list} -> adjusted={p_adj} "
                             f"(threads={thread_count}, policy={current_policy})")
                p_list = p_adj

                # Compute P (cumulative product of adjusted p)
                P_list = []
                acc = 1
                for pi_val in p_list:
                    acc *= pi_val
                    P_list.append(acc)
                
                # Create hardware parameters for this thread count and subset
                hw_params = {
                    'system_string': system_name,
                    'threads': thread_count,
                    'd': len(g_list),
                    'd_numa': d_numa,
                    'r_scalar': r_scalar if r_scalar is not None else 0.0,
                    'r_vec': r_vec if r_vec is not None else 0.0,
                    'g': g_list,
                    'ls': ls_list,
                    'm': m_list,
                    'p': p_list,
                    'P': P_list,
                    'level_names': level_subset
                }
                
                # Store in dictionary with unique key
                key = f"{current_policy}_{thread_count}_{subset_idx}"
                hw_params_dict[key] = hw_params
                
                # Add to header file
                if write_cpp_header:
                    header_content += f"""    // Hardware parameters for {thread_count} thread(s), policy: {current_policy}, subset: {', '.join(level_subset)}
   {{
    /* d            */ {hw_params['d']},
    /* d_numa       */ {hw_params['d_numa']},
    /* r_scalar     */ {hw_params['r_scalar']:e},
    /* r_vec        */ {hw_params['r_vec']:e},
    /* g            */ {{{', '.join([f"{val:e}" for val in hw_params['g']])}}},
    /* ls           */ {{{', '.join([f"{val:e}" for val in hw_params['ls']])}}},
    /* m            */ {{{', '.join([f"{val:d}" for val in hw_params['m']])}}},
    /* p            */ {{{', '.join([f"{val:d}" for val in hw_params['p']])}}},
    /* P            */ {{{', '.join([f"{val:d}" for val in hw_params['P']])}}}}}{',' if model_index < len(policies_to_process) * len(threads_list) * len(all_level_subsets) - 1 else ''}
"""
                
                # Store metadata for configuration structure
                all_hw_models.append(hw_params)
                all_thread_ids.append(threads_list.index(thread_count))
                
                # Calculate bitmask for this subset
                bitmask = 0
                for level_name in level_subset:
                    level_idx = base_levels.index(level_name)
                    bitmask |= (1 << level_idx)
                all_level_bitmasks.append(bitmask)
                
                # Policy ID (index in the policies list)
                all_policy_ids.append(policy_idx)
                
                model_index += 1
                
                # Print summary for this subset
                logger.debug(f"    Levels: {len(level_subset)}, d={hw_params['d']}, d_numa={hw_params['d_numa']}, bitmask={bitmask:0{len(base_levels)}b}, policy_id={policy_idx}")
            
            # Log subset processing summary for this thread/policy
            subset_summary = {
                'Threads': thread_count,
                'Policy': current_policy,
                'Subsets': len(all_level_subsets),
                'Models': len([k for k in hw_params_dict.keys() if k.startswith(f"{current_policy}_{thread_count}_")])
            }
            logger.debug(f"Completed {thread_count} threads, {current_policy} policy: {subset_summary['Models']} models from {subset_summary['Subsets']} subsets")
    
    # Create hardware parameter plots before writing header file
    if write_cpp_header and keep_levels:
        from hw_characterization.scripts.plotters import plot_hw_params_by_level
        
        # Create plots directory in the results structure
        plots_base_dir = f"{RESULTS_DIR}/{system_name}"
        os.makedirs(plots_base_dir, exist_ok=True)
        
        logger.debug(f"\nCreating hardware parameter plots for levels: {keep_levels}")
        # Pass LR coefficients and original data if available
        lr_coeffs_dict = {}
        original_data_dict = {}
        adjusted_monotonic_dict = {}
        if create_sync_level == "empirical_poly":
            for policy in policies_to_process:
                lr_key = f'lr_coeffs_{policy}'
                if lr_key in globals():
                    lr_coeffs_dict[policy] = globals()[lr_key]
                    # Extract original GLOBAL_SYNC data for this policy
                    if all_combined_df is not None:
                        # Debug: check what policies are available
                        available_policies = all_combined_df.get('alloc_policy', 'N/A').unique()
                        logger.debug(f"  Available policies in data: {available_policies}")
                        
                        # Filter for this specific policy and GLOBAL_SYNC level
                        policy_sync_data = all_combined_df[
                            (all_combined_df['alloc_policy'] == policy) & 
                            (all_combined_df['level_name'] == 'GLOBAL_SYNC')
                        ]
                    else:
                        policy_sync_data = df[
                            (df['alloc_policy'] == policy) & 
                            (df['level_name'] == 'GLOBAL_SYNC')
                        ]
                    
                    logger.debug(f"  Found {len(policy_sync_data)} GLOBAL_SYNC data points for policy {policy}")
                    if not policy_sync_data.empty:
                        # Keep all data points (both consecutive_bytes values)
                        logger.debug(f"  Thread counts: {policy_sync_data['threads'].values}")
                        logger.debug(f"  Li values: {policy_sync_data['Li'].values}")
                        logger.debug(f"  Consecutive bytes: {policy_sync_data['consecutive_bytes'].values}")
                        
                        original_data_dict[policy] = {
                            'threads': policy_sync_data['threads'].values,
                            'ls_values': policy_sync_data['Li'].values,
                            'consecutive_bytes': policy_sync_data['consecutive_bytes'].values
                        }
        
        # Collect adjusted monotonic data for empirical_monotony
        if create_sync_level == "empirical_monotony":
            for policy in policies_to_process:
                adj_mon_key = f'adjusted_monotonic_{policy}'
                if adj_mon_key in globals():
                    adjusted_monotonic_dict[policy] = globals()[adj_mon_key]
                    logger.debug(f"  Found adjusted monotonic data for policy {policy}: {len(adjusted_monotonic_dict[policy]['threads'])} points")
                    
                    # Also store original data for plotting
                    if all_combined_df is not None:
                        policy_sync_data = all_combined_df[
                            (all_combined_df.get('alloc_policy', policy) == policy) & 
                            (all_combined_df['level_name'] == 'GLOBAL_SYNC')
                        ]
                        if not policy_sync_data.empty:
                            original_data_dict[policy] = {
                                'threads': policy_sync_data['threads'].values,
                                'ls_values': policy_sync_data['Li'].values,
                                'consecutive_bytes': policy_sync_data.get('consecutive_bytes', [None] * len(policy_sync_data)).values
                            }
                            logger.debug(f"  Also stored original data for policy {policy}: {len(original_data_dict[policy]['threads'])} points")
        
        plot_files = plot_hw_params_by_level(hw_params_dict, system_name, keep_levels, plots_base_dir, lr_coeffs_dict, original_data_dict, adjusted_monotonic_dict)
        for plot_file in plot_files:
            track_file(plot_file)
        
        # Log plotting summary
        logger.debug(f"Generated {len(plot_files)} hardware parameter plots")
        logger.debug(f"Plot directory: {plots_base_dir}")
        for policy in policies_to_process:
            policy_plots = [f for f in plot_files if f'/policy_plots/{policy}/' in f or f'/{policy}/' in f]
            logger.debug(f"  {policy} policy: {len(policy_plots)} plots")
    
    # Finalize and write the header file
    if write_cpp_header:
        if 'header_file' not in locals():
            logger.error("ERROR: write_cpp_header=True but header_file was not initialized. This should not happen!")
            raise RuntimeError("Header file path not initialized")
        # Close the hw_models array and create the HWParameter_configurations struct
        header_content += """
};

// HWParameter_configurations struct
const cost_models::HW_model::HWParameter_configurations dis_system_params = {
    /* threads_options_num */ {""" + ', '.join([str(t) for t in threads_list]) + """},
    /* policy_options_str */ {""" + ', '.join([f'"{policy}"' for policy in policies_to_process]) + """},
    /* level_options_str  */ {""" + ', '.join([f'"{name}"' for name in base_levels]) + """},
    /* hw_model_thread_id */ {""" + ', '.join([str(tid) for tid in all_thread_ids]) + """},
    /* hw_model_policy_id */ {""" + ', '.join([str(pid) for pid in all_policy_ids]) + """},
    /* hw_model_level_bithash */ {""" + ', '.join([str(bm) for bm in all_level_bitmasks]) + """},
    /* hw_models */ hw_models_vector
};

#endif // HW_PARAMS_""" + system_id.upper() + """_HPP
"""
        
        # Write to file
        try:
            with open(header_file, 'w') as f:
                f.write(header_content)
            logger.info(f"C++ hardware parameters written to: {header_file}")
            logger.debug(f"Total models generated: {len(all_hw_models)}")
            logger.debug(f"Allocation policies: {len(policies_to_process)} ({', '.join(policies_to_process)})")
            logger.debug(f"Thread configurations: {len(threads_list)}")
            logger.debug(f"Level subsets: {len(all_level_subsets)}")
            logger.debug(f"Models per policy per thread: {len(all_level_subsets)}")
            logger.debug(f"Models per policy: {len(all_hw_models) // len(policies_to_process) if policies_to_process else 0}")
        except Exception as e:
            logger.error(f"Failed to write header file {header_file}: {e}")
            raise
    
    return hw_params_dict


def generate_default_threads(system_name):
    """
    Generate default thread counts based on the system's topology file.
    
    Produces the sequence: {1, 2, 4, ..., system_max/4, system_max/2, system_max}
    where system_max is the product of all pi values in the topology (i.e. total cores),
    and the leading part doubles (powers of 2) up to system_max/4.
    
    Args:
        system_name: Name of the system (used to locate the topology file)
        
    Returns:
        list: Sorted list of unique thread counts
    """
    from hw_characterization.scripts.parser import parse_pahis_log_to_df, get_topology_file_path

    topology_file = get_topology_file_path(system_name)
    if not os.path.exists(topology_file):
        print(f"Warning: Topology file {topology_file} not found, "
              f"using fallback thread list [1,2,4,8,16,32,64]")
        return [1, 2, 4, 8, 16, 32, 64]

    _, system_max, _, _, _, _ = parse_pahis_log_to_df(topology_file)

    if system_max is None or system_max <= 0:
        print(f"Warning: Could not determine system_max from {topology_file}, "
              f"using fallback thread list [1,2,4,8,16,32,64]")
        return [1, 2, 4, 8, 16, 32, 64]

    threads = set()

    # Powers of 2: 1, 2, 4, ... up to system_max // 4
    t = 1
    while t <= system_max // 4:
        threads.add(t)
        t *= 2

    # Add the three final anchor points
    if system_max // 4 > 0:
        threads.add(system_max // 4)
    if system_max // 2 > 0:
        threads.add(system_max // 2)
    threads.add(system_max)

    return sorted(threads)


def _parse_args():
    """Parse command line arguments for the microbenchmarking script."""
    parser = argparse.ArgumentParser(description="Run PaHiS microbenchmarks and extract models")
    parser.add_argument(
        "--model_coeffs",
        type=str,
        default="header",
        choices=["microbenchmarks", "header", "topology"],
        help="Source of model coefficients: "
             "'header' (default, use a pre-initialized hw_params_{system}.hpp header file from models/ directory, "
             "skip microbenchmark execution and header generation), "
             "'microbenchmarks' (run full microbenchmark suite to generate hardware parameters), "
             "'topology' (NOT IMPLEMENTED - reserved for future loading of model coefficients directly from topology files)",
    )
    parser.add_argument(
        "--threads",
        type=str,
        default="",
        help="Comma-separated list of thread counts to benchmark (e.g., 1,2,4,8). "
             "If not provided, defaults to {1,2,...,system_max/4, system_max/2, system_max} "
             "based on the system topology.",
    )
    parser.add_argument(
        "--consecutive-bytes",
        type=str,
        default="1,64",
        help="Comma-separated list of consecutive byte sizes (e.g., 1,64)",
    )
    parser.add_argument(
        "--keep-levels",
        type=str,
        default="",
        help="Comma-separated list of level names to keep (empty means keep all)",
    )
    parser.add_argument(
        "--create-sync-level",
        type=str,
        default="empirical_poly",
        choices=["", "empirical", "empirical_poly", "empirical_monotony"],
        help="GLOBAL_SYNC level creation method: "
             "'\' (default, no sync level), "
             "'empirical' (empirical benchmarks, use measured Li directly), "
             "'empirical_poly' (empirical benchmarks with polynomial fitting for thread scaling), "
             "'empirical_monotony' (empirical benchmarks with averaging per thread count, monotonicity enforcement, and direct use of adjusted values - no polynomial fitting)",
    )
    parser.add_argument(
        "--core-policy",
        type=str,
        default="close",
        help="NUMA allocation policy(ies): 'close' (threads on consecutive cores), 'spread' (threads distributed across NUMA nodes), or comma-separated list like 'close,spread' (default: close)",
    )
    parser.add_argument(
        "--generate-subsets",
        action="store_true",
        default=False,
        help="Generate all subset combinations of levels (always including NodeMem). "
             "By default, only the exact specified levels are used (default: False).",
    )
    parser.add_argument(
        "--bw-kernels",
        type=str,
        default=DEFAULT_BW_KERNELS,
        help=f"Comma-separated list of kernel names to use for bandwidth (gi) measurements. "
             f"Available: {', '.join(KERNEL_NAMES)}. "
             f"(default: {DEFAULT_BW_KERNELS})",
    )
    parser.add_argument(
        "--kernel-aggregator",
        type=str,
        default="mean",
        choices=["min", "max", "mean"],
        help="Aggregation method used across kernel measurements (copy, daxpy, triad, etc.) "
             "when computing a single gi/ls value per level. Applied consistently in "
             "monotonicity enforcement, hw_topo summary, and final parameter extraction. "
             "(default: mean)",
    )
    parser.add_argument(
        "--clean-dataset",
        type=str,
        default="no",
        choices=["no", "models", "plots", "results", "all"],
        help="Clean benchmark data under hw_characterization/results/{system}/ before running. "
             "'no' (default): do nothing. "
             "'models': delete all models/ subdirs (HW_model.csv files). "
             "'plots': delete all plots/ subdirs. "
             "'results': delete all results/ subdirs (per-kernel CSV logs). "
             "'all': delete everything under the system results directory (full clean).",
    )
    parser.add_argument(
        "--verbose",
        action="store_true",
        help="Enable verbose debug output to console (default: False, debug only in log file)",
    )
    return parser.parse_args()


def clean_dataset(system_name, mode):
    """
    Clean benchmark data under hw_characterization/results/{system}/ before running.
    
    Args:
        system_name: Name of the system (e.g., "AMDEPYC9634")
        mode: Cleaning mode — one of "models", "plots", "results", or "all"
    """
    import shutil
    system_results_dir = os.path.join(RESULTS_DIR, system_name)
    
    if not os.path.isdir(system_results_dir):
        logger.info(f"Clean: nothing to clean — {system_results_dir} does not exist")
        return
    
    logger.info("")
    logger.info("="*60)
    logger.info("DATASET CLEANING")
    logger.info("="*60)
    
    logger.info(f"  Mode: {mode}")
    logger.info(f"  Target: {system_results_dir}")
    
    if mode == "all":
        # Delete everything under the system results directory (preserve output_headers)
        removed_dirs = []
        removed_files = []
        for entry in sorted(os.listdir(system_results_dir)):
            entry_path = os.path.join(system_results_dir, entry)
            if os.path.isdir(entry_path):
                shutil.rmtree(entry_path)
                removed_dirs.append(entry)
            elif os.path.isfile(entry_path):
                os.remove(entry_path)
                removed_files.append(entry)
        if removed_dirs:
            logger.info(f"  Removed directories: {', '.join(removed_dirs)}")
        if removed_files:
            logger.info(f"  Removed files: {', '.join(removed_files)}")
        total = len(removed_dirs) + len(removed_files)
        logger.info(f"  Total: {total} entries removed")
    else:
        # Delete specific subdirectory type (models, plots, or results)
        removed_paths = []
        for root, dirs, _files in os.walk(system_results_dir):
            if os.path.basename(root) == mode:
                rel_path = os.path.relpath(root, system_results_dir)
                # Count files inside before removing
                file_count = sum(len(f) for _, _, f in os.walk(root))
                shutil.rmtree(root)
                removed_paths.append((rel_path, file_count))
        if removed_paths:
            for rel_path, file_count in removed_paths:
                logger.info(f"  Removed: {rel_path}/ ({file_count} files)")
        logger.info(f"  Total: {len(removed_paths)} '{mode}/' subdirectories removed")
    
    logger.info("="*60)


def compile_models_with_hw_params(system_name):
    """
    Compile models directory with generated hardware parameters.
    Uses the same include strategy as ALP exporting.
    """
    import subprocess
    import shutil
    import os
    from pathlib import Path
    
    # Get paths
    repo_root = Path.cwd()
    models_dir = repo_root / "models"
    build_dir = Path(os.environ.get("HW_BUILD_DIR", "./build"))
    models_build_dir = build_dir / "model_build"  # Use model build subdirectory
    
    # Convert system name to valid system_id (same logic as CMakeLists.txt)
    system_id = system_name.lower().replace("-", "_").replace(" ", "_")
    
    # Check if hardware parameter file exists in output_headers
    hw_params_source = Path(CPP_MODEL_DIR) / f"hw_params_{system_id}.hpp"
    hw_params_target = models_dir / f"hw_params_{system_id}.hpp"
    
    if not hw_params_source.exists():
        logger.debug(f"Warning: Hardware parameter file {hw_params_source} not found.")
        logger.debug("Models compilation skipped - run hardware characterization first.")
        return False
    
    # Copy hardware parameter file to models directory (same strategy as ALP export)
    shutil.copy2(hw_params_source, hw_params_target)
    logger.debug(f"Copied hardware parameters: {hw_params_source} -> {hw_params_target}")
    
    logger.debug(f"Compiling models with hardware parameters: {hw_params_target}")
    
    try:
        # Ensure main build directory exists
        models_build_dir.mkdir(exist_ok=True)
        
        # Run CMake with HW_MODEL environment variable (build into main build directory)
        cmake_cmd = [
            "cmake", 
            "-S", str(models_dir), 
            "-B", str(models_build_dir),
            f"-DHW_MODEL={system_name}"
        ]
        
        logger.debug(f"Running: {' '.join(cmake_cmd)}")
        result = subprocess.run(cmake_cmd, capture_output=True, text=True, cwd=repo_root)
        
        if result.returncode != 0:
            logger.debug(f"CMake configuration failed:")
            logger.debug(f"stdout: {result.stdout}")
            logger.debug(f"stderr: {result.stderr}")
            return False
        
        # Compile
        make_cmd = ["make", "-C", str(models_build_dir), "-j"]
        logger.debug(f"Running: {' '.join(make_cmd)}")
        result = subprocess.run(make_cmd, capture_output=True, text=True)
        
        if result.returncode != 0:
            logger.debug(f"Compilation failed:")
            logger.debug(f"stdout: {result.stdout}")
            logger.debug(f"stderr: {result.stderr}")
            return False
        
        logger.debug("✓ Models compiled successfully!")
        
        # List generated executables
        bin_dir = models_build_dir / "bin"
        executables = list(bin_dir.glob("*")) if bin_dir.exists() else []
        if executables:
            logger.debug("Generated executables:")
            for exe in executables:
                if exe.is_file() and exe.stat().st_mode & 0o111:  # Check if executable
                    logger.debug(f"  - {exe.name}")
        
        return True
        
    except Exception as e:
        logger.debug(f"Error compiling models: {e}")
        return False


def run_model_unit_tests(system_name):
    """
    Run unit tests for model compilation and execution.
    This validates that the models compile and run correctly with the generated hardware parameters.
    """
    import subprocess
    import sys
    from pathlib import Path
    
    logger.debug(f"Running unit tests for system: {system_name}")
    
    # Get paths
    repo_root = Path.cwd()
    tests_dir = repo_root / "tests"
    test_file = tests_dir / "test_model_execution.py"
    
    if not test_file.exists():
        logger.debug(f"Warning: Unit test file {test_file} not found.")
        logger.debug("Unit tests skipped.")
        return False
    
    try:
        # Run pytest on the model execution tests
        pytest_cmd = [
            sys.executable, "-m", "pytest", 
            str(test_file),
            "-v",  # verbose output
            "--tb=short",  # shorter traceback format
            "-x",  # stop on first failure
        ]
        
        logger.debug(f"Running: {' '.join(pytest_cmd)}")
        result = subprocess.run(pytest_cmd, capture_output=True, text=True, cwd=repo_root)
        
        if result.returncode == 0:
            logger.debug("✓ All model unit tests passed!")
            # Log a summary of test results
            if "PASSED" in result.stdout:
                logger.debug("Test Summary:")
                for line in result.stdout.split('\n'):
                    if "PASSED" in line or "FAILED" in line or "SKIPPED" in line:
                        logger.debug(f"  {line.strip()}")
            return True
        else:
            logger.debug("✗ Some model unit tests failed:")
            logger.debug("STDOUT:")
            logger.debug(result.stdout)
            logger.debug("STDERR:")
            logger.debug(result.stderr)
            return False
            
    except Exception as e:
        logger.debug(f"Error running unit tests: {e}")
        return False


def _main():
    """Main function for running hardware microbenchmarks."""
    # Parse args first to check for verbose flag
    args = _parse_args()
    
    # Setup logging first
    log_file = setup_logging(verbose=args.verbose)
    logger.info("="*60)
    logger.info("PaHiS - Hardware Characterization")
    logger.info("="*60)
    logger.info(f"Detailed log: {log_file}")
    logger.info("")
    
    try:
        system_string_local = _require_env("HW_MODEL")
        
        # Handle --model_coeffs mode
        if args.model_coeffs == "topology":
            logger.error("ERROR: --model_coeffs=topology is not yet implemented.")
            logger.error("This mode will support loading model coefficients directly from topology files in a future release.")
            raise NotImplementedError(
                "--model_coeffs=topology is not yet implemented. "
                "Use 'microbenchmarks' to run the full benchmark suite, or "
                "'header' to use a pre-initialized hw_params header file."
            )
        
        logger.info(f"Model coefficients source: {args.model_coeffs}")
        
        # args already parsed above for verbose flag
        if args.threads:
            threads_list = [int(x) for x in args.threads.split(',') if x.strip()]
        else:
            threads_list = generate_default_threads(system_string_local)
            logger.info(f"Auto-generated thread list from topology: {threads_list}")
        k_list = [int(x) for x in args.consecutive_bytes.split(',') if x.strip()]
        keep_levels = [s for s in args.keep_levels.split(',') if s.strip()] if args.keep_levels else None
        
        # Parse bandwidth kernel selection
        bw_kernel_names = [s.strip() for s in args.bw_kernels.split(',') if s.strip()]
        for name in bw_kernel_names:
            if name not in KERNEL_NAME_TO_IDX:
                raise ValueError(f"Unknown bandwidth kernel '{name}'. Available: {', '.join(KERNEL_NAMES)}")
        bw_kernel_indices = [KERNEL_NAME_TO_IDX[name] for name in bw_kernel_names]
        
        # Parse allocation policies (support comma-separated list)
        alloc_policies = [s.strip() for s in args.core_policy.split(',') if s.strip()]
        
        # Display input parameters section
        logger.info("")
        logger.info("="*60)
        logger.info("INPUT PARAMETERS")
        logger.info("="*60)
        logger.info(f"System: {system_string_local}")
        logger.info(f"Thread Configurations: {len(threads_list)} ({min(threads_list)}-{max(threads_list)} threads)")
        logger.info(f"Consecutive Bytes: {len(k_list)} configs {k_list}")
        logger.info(f"Allocation Policies: {len(alloc_policies)} ({', '.join(alloc_policies)})")
        logger.info(f"Memory Levels: {keep_levels if keep_levels else 'All available levels'}")
        logger.info(f"BW Kernels: {', '.join(bw_kernel_names)} (indices {bw_kernel_indices})")
        logger.info(f"Sync Level Method: {args.create_sync_level if args.create_sync_level else 'None'}")
        logger.info(f"Subset Generation: {'Enabled (all combinations)' if args.generate_subsets else 'Disabled (exact levels only)'}")
        logger.info(f"Model Coefficients: {args.model_coeffs}")
        logger.info("="*60)
        
        # Log configuration summary
        config_details = {
            'System': system_string_local,
            'Model Coefficients': args.model_coeffs,
            'Threads': f"{len(threads_list)} configs ({min(threads_list)}-{max(threads_list)})",
            'Consecutive Bytes': f"{len(k_list)} configs ({k_list})",
            'Policies': f"{len(alloc_policies)} ({', '.join(alloc_policies)})",
            'Keep Levels': keep_levels if keep_levels else "All levels",
            'BW Kernels': f"{', '.join(bw_kernel_names)}"
        }
        log_summary("Configuration", "SUCCESS", config_details)
        
        # Log detailed config to file
        logger.debug(f"Processing allocation policies: {alloc_policies}")
        logger.debug(f"Threads list: {threads_list}")
        logger.debug(f"Consecutive bytes list: {k_list}")
        logger.debug(f"Keep levels: {keep_levels}")
        logger.debug(f"Create sync level: {args.create_sync_level}")
        logger.debug(f"Generate subsets: {args.generate_subsets}")

        # Use configured directories
        results_dir = RESULTS_DIR  # Use hardcoded results directory
        cpp_model_dir = CPP_MODEL_DIR
        
        # Dataset cleaning (runs before any benchmarks)
        if args.clean_dataset != "no":
            clean_dataset(system_string_local, args.clean_dataset)
        
        if args.model_coeffs == "header":
            # ============================================================
            # HEADER MODE: Use a pre-initialized hw_params header file
            # Skip microbenchmark execution and header generation entirely
            # ============================================================
            import shutil
            
            system_id = system_string_local.replace('-', '_').replace(' ', '_').lower()
            header_source = Path.cwd() / "models" / f"hw_params_{system_id}.hpp"
            header_dest = Path(CPP_MODEL_DIR) / f"hw_params_{system_id}.hpp"
            
            logger.info("")
            logger.info("="*60)
            logger.info("LOADING PRE-INITIALIZED HEADER")
            logger.info("="*60)
            logger.info(f"Source: {header_source}")
            
            if not header_source.exists():
                logger.error(f"ERROR: Pre-initialized header file not found: {header_source}")
                logger.error(f"Expected file: models/hw_params_{system_id}.hpp")
                logger.error("Please ensure the header file exists in the models/ directory,")
                logger.error("or use --model_coeffs=microbenchmarks to generate it.")
                raise FileNotFoundError(
                    f"Pre-initialized header file not found: {header_source}. "
                    f"Expected models/hw_params_{system_id}.hpp in the repository root."
                )
            
            # Copy header to the output_headers directory so compile step can find it
            os.makedirs(os.path.dirname(header_dest), exist_ok=True)
            shutil.copy2(header_source, header_dest)
            track_file(str(header_dest))
            
            logger.info(f"Copied header to: {header_dest}")
            
            # Parse and display header metadata from comment block
            header_meta = {}
            with open(header_source, 'r') as f:
                for line in f:
                    line = line.strip()
                    if not line.startswith('//'):
                        if line and not line.startswith('#'):
                            break  # End of comment block
                        continue
                    comment = line[2:].strip()
                    if comment.startswith('Auto-generated hardware parameters for'):
                        header_meta['System'] = comment.split('for', 1)[1].strip()
                    elif comment.startswith('Allocation policies:'):
                        header_meta['Allocation policies'] = comment.split(':', 1)[1].strip()
                    elif comment.startswith('Base levels:'):
                        header_meta['Base levels'] = comment.split(':', 1)[1].strip()
                    elif 'subset combinations' in comment.lower():
                        header_meta['Subset mode'] = comment
                    elif comment.startswith('Generated on'):
                        header_meta['Generated on'] = comment.split('on', 1)[1].strip()
                    elif comment.startswith('python3 runner.py') or comment.startswith('python runner.py'):
                        header_meta['Command'] = comment
            
            if header_meta:
                logger.info("")
                logger.info("Header parameters:")
                for key, value in header_meta.items():
                    logger.info(f"  {key}: {value}")
            
            log_summary("Load Pre-initialized Header", "SUCCESS", {
                'Source': str(header_source),
                'Destination': str(header_dest),
                **header_meta,
            }, [str(header_dest)])
            logger.info("="*60)
        
        if args.model_coeffs == "microbenchmarks":
            # ============================================================
            # MICROBENCHMARKS MODE: Full benchmark suite (default behavior)
            # ============================================================
        
            # Display microbenchmark execution section
            logger.info("")
            logger.info("="*60)
            logger.info("MICROBENCHMARK EXECUTION")
            logger.info("="*60)
            total_benchmarks = len(alloc_policies) * len(threads_list) * len(k_list)
            logger.info(f"Total benchmark configurations: {total_benchmarks}")
            logger.info(f"  - Policies: {len(alloc_policies)}")
            logger.info(f"  - Thread configs: {len(threads_list)}")
            logger.info(f"  - Consecutive bytes: {len(k_list)}")
            logger.info("="*60)
            
            # Collect data from all allocation policies first
            all_combined_df = pd.DataFrame()
            benchmark_count = 0
            
            for policy_idx, alloc_policy in enumerate(alloc_policies):
                logger.info(f"")
                logger.info(f"Policy {policy_idx + 1}/{len(alloc_policies)}: {alloc_policy.upper()}")
                logger.info(f"  Thread configs: {len(threads_list)} {threads_list}")
                logger.info(f"  Consecutive bytes: {len(k_list)} {k_list}")
                
                logger.debug(f"\n{'='*60}")
                logger.debug(f"Processing allocation policy: {alloc_policy}")
                logger.debug(f"{'='*60}")
                
                combined_df = pd.DataFrame()
                policy_files = []
                
                # Track what we're doing
                running_benchmarks = []
                loading_benchmarks = []
                
                for num_threads in threads_list:
                    for consecutive_bytes in k_list:
                        benchmark_count += 1
                        filename = f"{results_dir}/{system_string_local}/{alloc_policy}/k_{consecutive_bytes}/t_{num_threads}/models/HW_model.csv"
                        
                        if not os.path.exists(filename):
                            running_benchmarks.append(f"{num_threads}t,k={consecutive_bytes}B")
                            logger.debug(f"Running benchmark {benchmark_count}/{total_benchmarks}: {alloc_policy}, {num_threads}t, k={consecutive_bytes}")
                            config = BenchmarkConfig(
                                system_name=system_string_local,
                                consecutive_bytes=consecutive_bytes,
                                num_threads=num_threads,
                                create_sync_level=args.create_sync_level,
                                alloc_policy=alloc_policy,
                                bw_kernel_indices=bw_kernel_indices,
                                kernel_aggregator=args.kernel_aggregator
                            )
                            output = create_hw_tree(config)
                            os.makedirs(os.path.dirname(filename), exist_ok=True)
                            output.to_csv(filename, index=False)
                            track_file(filename)
                            
                            # Generate plots
                            plot_files = plot_PaHiS_submodel(output, system_string_local, alloc_policy, consecutive_bytes, num_threads, config.system_cores,
                                                                r_scalar=config.r_scalar if config.r_scalar else 0.0,
                                                                r_vec=config.r_vec if config.r_vec else 0.0)
                            if plot_files:
                                for plot_file in plot_files:
                                    track_file(plot_file)
                                    policy_files.append(plot_file)
                        else:
                            loading_benchmarks.append(f"{num_threads}t,k={consecutive_bytes}B")
                            logger.debug(f"Loading existing benchmark data from {filename}")
                            output = pd.read_csv(filename)
                        
                        # Keep only rows matching the selected bandwidth kernel indices
                        if 'idx' in output.columns:
                            # Rows with idx=0 (chase) are not in the HW_model.csv;
                            # GLOBAL_SYNC rows have idx=0 and must be kept.
                            bw_rows_mask = output['idx'].isin(bw_kernel_indices)
                            non_bw_rows_mask = ~output['idx'].isin([v for v in KERNEL_NAME_TO_IDX.values() if v >= 2])
                            output = output[bw_rows_mask | non_bw_rows_mask]
                        
                        # Add alloc_policy column to distinguish between policies
                        output['alloc_policy'] = alloc_policy
                        combined_df = pd.concat([combined_df, output], ignore_index=True)
                
                # Show summary of what was done
                if running_benchmarks:
                    logger.info(f"  Running benchmarks: {len(running_benchmarks)} configs")
                    logger.debug(f"  Running: {', '.join(running_benchmarks)}")
                if loading_benchmarks:
                    logger.info(f"  Loading existing data: {len(loading_benchmarks)} configs")
                    logger.debug(f"  Loading: {', '.join(loading_benchmarks)}")
                
                # Log policy completion
                policy_details = {
                    'Threads': f"{len(threads_list)} configs",
                    'Consecutive Bytes': f"{len(k_list)} configs", 
                    'Data Points': len(combined_df),
                    'Files Generated': len(policy_files)
                }
                log_summary(f"Policy {alloc_policy}", "SUCCESS", policy_details, policy_files)
                
                logger.debug(f"\nCombined data for {alloc_policy} policy:")
                logger.debug(combined_df.to_markdown())
                
                # Add to overall combined dataframe
                all_combined_df = pd.concat([all_combined_df, combined_df], ignore_index=True)
        
            globals()['global_sync_boundaries'] = {}
    
            # Handle empirical_poly or empirical_monotony for GLOBAL_SYNC if requested
            if args.create_sync_level == "empirical_poly" or args.create_sync_level == "empirical_monotony":
                logger.debug(f"\n{'='*60}")
                mode_desc = "empirical_monotony (averaging + monotonicity enforcement)" if args.create_sync_level == "empirical_monotony" else "empirical_poly (polynomial fitting)"
                logger.debug(f"Processing {args.create_sync_level} for GLOBAL_SYNC: {mode_desc}")
                logger.debug(f"{'='*60}")
                
                # Extract all GLOBAL_SYNC data from ALL policies
                all_sync_data = all_combined_df[all_combined_df['level_name'] == 'GLOBAL_SYNC']
                
                if all_sync_data.empty:
                    logger.debug(f"  Warning: No GLOBAL_SYNC data found")
                else:
                    # Debug: show what policies are in the data
                    if 'alloc_policy' in all_sync_data.columns:
                        unique_policies_in_data = all_sync_data['alloc_policy'].unique()
                        logger.debug(f"  Policies found in all_combined_df: {list(unique_policies_in_data)}")
                        logger.debug(f"  Total data points: {len(all_sync_data)}")
                        for policy in unique_policies_in_data:
                            policy_data = all_sync_data[all_sync_data['alloc_policy'] == policy]
                            logger.debug(f"    {policy}: {len(policy_data)} data points")
                    else:
                        logger.debug(f"  Warning: 'alloc_policy' column not found in all_sync_data")
                        logger.debug(f"  Total data points: {len(all_sync_data)}")
                
                # Step 1: Identify boundary points and average across ALL policies
                if all_sync_data.empty:
                    logger.debug(f"  ERROR: Cannot proceed - no GLOBAL_SYNC data available")
                    boundary_averages = {}
                else:
                    unique_threads = sorted(all_sync_data['threads'].unique())
                    min_threads = min(unique_threads)
                    max_threads = max(unique_threads)
                    
                    logger.debug(f"  Step 1: Boundary point averaging (across ALL policies)")
                    logger.debug(f"    Boundary threads: min={min_threads}, max={max_threads}")
                    logger.debug(f"    All policies in dataset: {list(all_sync_data['alloc_policy'].unique()) if 'alloc_policy' in all_sync_data.columns else 'N/A'}")
                    
                    # Average boundary points across all policies
                    boundary_averages = {}
                    for boundary_threads in [min_threads, max_threads]:
                        boundary_data = all_sync_data[all_sync_data['threads'] == boundary_threads]
                        boundary_ls_values = boundary_data['Li'].values
                        boundary_avg = np.mean(boundary_ls_values)
                        
                        # Count unique policies and measurements per policy
                        if 'alloc_policy' in boundary_data.columns:
                            boundary_policies = boundary_data['alloc_policy'].unique()
                            policy_counts = {}
                            for policy in boundary_policies:
                                policy_mask = boundary_data['alloc_policy'] == policy
                                policy_counts[policy] = np.sum(policy_mask)
                        else:
                            boundary_policies = ['N/A']
                            policy_counts = {}
                        
                        boundary_counts = len(boundary_ls_values)
                        boundary_averages[boundary_threads] = boundary_avg
                        logger.debug(f"      threads={boundary_threads}: avg_ls={boundary_avg:.6e} (from {boundary_counts} total measurements)")
                        logger.debug(f"        Measurements by policy: {policy_counts}")
                        logger.debug(f"        Policies involved: {list(boundary_policies)} ({len(boundary_policies)} unique)")
                        logger.debug(f"        All measurements: {[f'{ls:.6e}' for ls in boundary_ls_values]}")
    
                    globals()['global_sync_boundaries'] = {
                        'min_threads': min_threads,
                        'max_threads': max_threads,
                        'ls': boundary_averages.copy()
                    }
                    logger.debug(f"  Stored global GLOBAL_SYNC boundary averages for runtime use: {globals()['global_sync_boundaries']}")
                
                # Store adjusted monotonic data for empirical_monotony mode
                adjusted_monotonic_data = {}  # {policy: {threads: [list], ls: [list]}}
                
                # Step 2: For each policy, average non-boundary points and use boundary averages
                logger.debug(f"  Step 2: Policy-specific averaging (non-boundary points only)")
                for alloc_policy in alloc_policies:
                    logger.debug(f"Processing policy: {alloc_policy}")
                    
                    # Filter data for this policy and GLOBAL_SYNC level
                    policy_sync_data = all_combined_df[
                        (all_combined_df.get('alloc_policy', alloc_policy) == alloc_policy) & 
                        (all_combined_df['level_name'] == 'GLOBAL_SYNC')
                    ]
                    
                    if policy_sync_data.empty:
                        logger.debug(f"  Warning: No GLOBAL_SYNC data found for policy {alloc_policy}")
                        continue
                    
                    logger.debug(f"    Raw data points: {len(policy_sync_data)}")
                    logger.debug(f"    Raw data breakdown by thread count:")
                    for threads in sorted(policy_sync_data['threads'].unique()):
                        thread_data = policy_sync_data[policy_sync_data['threads'] == threads]
                        ls_values = thread_data['Li'].values
                        consecutive_bytes = thread_data['consecutive_bytes'].values if 'consecutive_bytes' in thread_data.columns else ['N/A'] * len(ls_values)
                        logger.debug(f"      threads={threads}: {len(ls_values)} measurements")
                        for i, (ls_val, k_val) in enumerate(zip(ls_values, consecutive_bytes)):
                            logger.debug(f"        Point {i+1}: ls={ls_val:.6e}, consecutive_bytes={k_val}")
                        logger.debug(f"        Summary: ls values = {[f'{ls:.6e}' for ls in ls_values]}")
                        if len(ls_values) > 1:
                            logger.debug(f"        Range: min={ls_values.min():.6e}, max={ls_values.max():.6e}, spread={ls_values.max()-ls_values.min():.6e}")
                    
                    # For each thread count ti, average all ls(ti) values (within this policy)
                    # But use boundary averages for min/max threads
                    averaged_data = policy_sync_data.groupby('threads')['Li'].agg(['mean', 'count', 'std']).reset_index()
                    averaged_data.columns = ['threads', 'ls_avg', 'count', 'std']
                    
                    # Replace boundary thread averages with cross-policy averages
                    for boundary_threads in [min_threads, max_threads]:
                        if boundary_threads in averaged_data['threads'].values:
                            idx = averaged_data[averaged_data['threads'] == boundary_threads].index[0]
                            old_avg = averaged_data.loc[idx, 'ls_avg']
                            averaged_data.loc[idx, 'ls_avg'] = boundary_averages[boundary_threads]
                            logger.debug(f"      Replaced boundary threads={boundary_threads}: {old_avg:.6e} -> {boundary_averages[boundary_threads]:.6e} (cross-policy average)")
                    
                    logger.debug(f"    After averaging (non-boundary: policy-only, boundary: cross-policy):")
                    logger.debug(f"    Averaged data points: {len(averaged_data)}")
                    for _, row in averaged_data.iterrows():
                        std_str = f"{row['std']:.6e}" if not pd.isna(row['std']) else 'N/A'
                        threads_val = int(row['threads'])
                        ls_avg_val = row['ls_avg']
                        count_val = int(row['count'])
                        is_boundary = threads_val in [min_threads, max_threads]
                        boundary_tag = " (boundary - cross-policy avg)" if is_boundary else " (policy-only avg)"
                        logger.debug(f"      threads={threads_val}: ls_avg={ls_avg_val:.6e}{boundary_tag}")
                        logger.debug(f"        From {count_val} measurement(s), std={std_str}")
                        # Show what was averaged
                        original_points = policy_sync_data[policy_sync_data['threads'] == threads_val]['Li'].values
                        if len(original_points) > 1:
                            diff_from_avg = original_points - ls_avg_val
                            logger.debug(f"        Original points: {[f'{p:.6e}' for p in original_points]}")
                            logger.debug(f"        Differences from average: {[f'{d:+.6e}' for d in diff_from_avg]}")
                    
                    # Now use the merged data for fitting
                    # Handle endpoints: use boundary averages from Step 1
                    # Handle middle points: use averaged data from this policy
                    
                    combined_threads = []
                    combined_li = []
                    
                    logger.debug(f"  Step 3: Building combined dataset (boundary from cross-policy, middle from policy)")
                    
                    # Add boundary data using cross-policy averages from Step 1
                    for boundary_threads in [min_threads, max_threads]:
                        if boundary_threads in averaged_data['threads'].values:
                            combined_threads.append(boundary_threads)
                            combined_li.append(boundary_averages[boundary_threads])
                            logger.debug(f"    Boundary threads {boundary_threads}: using cross-policy average Li = {boundary_averages[boundary_threads]:.6e}")
                    
                    # Add middle data points from this policy only (excluding boundaries)
                    middle_data = averaged_data[
                        (averaged_data['threads'] != min_threads) & 
                        (averaged_data['threads'] != max_threads)
                    ]
                    
                    if not middle_data.empty:
                        combined_threads.extend(middle_data['threads'].values)
                        combined_li.extend(middle_data['ls_avg'].values)
                        logger.debug(f"    Middle data points: {len(middle_data)} from {alloc_policy} policy (policy-only averages)")
                        for _, row in middle_data.iterrows():
                            logger.debug(f"      threads={int(row['threads'])}: ls={row['ls_avg']:.6e}")
                    
                    # Convert to numpy arrays for monotonicity checking
                    threads_array_pre = np.array(combined_threads)
                    li_sync_array_pre = np.array(combined_li)
                    
                    # Sort by thread count for monotonicity checking
                    sort_idx = np.argsort(threads_array_pre)
                    threads_array_pre = threads_array_pre[sort_idx]
                    li_sync_array_pre = li_sync_array_pre[sort_idx]
                    
                    monotony_factor = 1.2
                    logger.debug(f"  Step 4: Enforcing weighted monotonicity (t_i > t_j -> ls(t_i) >= {monotony_factor:.2f} * ls(t_j))")
                    logger.debug(f"    Pre-monotonicity dataset: {len(threads_array_pre)} points")
                    logger.debug(f"      Threads: {threads_array_pre.tolist()}")
                    logger.debug(f"      LS values: {[f'{ls:.6e}' for ls in li_sync_array_pre]}")
                    
                    # Convert back to DataFrame format for monotonicity processing
                    averaged_data = pd.DataFrame({
                        'threads': threads_array_pre,
                        'ls_avg': li_sync_array_pre,
                        'count': [1] * len(threads_array_pre)  # Count doesn't matter for monotonicity
                    })
                    
                    # Sort by thread count (ascending)
                    averaged_data = averaged_data.sort_values('threads').reset_index(drop=True)
                    
                    # Check for violations and merge groups
                    merged_data = []
                    current_group = []
                    
                    logger.debug(f"    Checking monotonicity violations:")
                    logger.debug(f"    Processing {len(averaged_data)} averaged points in order...")
                    violation_found = False
                    
                    for idx, row in averaged_data.iterrows():
                        threads = int(row['threads'])
                        ls_avg = row['ls_avg']
                        
                        if not current_group:
                            # Start new group
                            current_group = [{'threads': threads, 'ls_avg': ls_avg, 'count': int(row['count'])}]
                            logger.debug(f"      [{idx+1}/{len(averaged_data)}] Starting new group: threads={threads}, ls={ls_avg:.6e}")
                        else:
                            # Calculate the averaged value of the current group (if merged, use average)
                            group_threads_list = [g['threads'] for g in current_group]
                            group_ls_list = [g['ls_avg'] for g in current_group]
                            if len(current_group) == 1:
                                # Single point in group, use its value
                                last_group_ls = current_group[0]['ls_avg']
                                last_group_threads = current_group[0]['threads']
                            else:
                                # Multiple points merged, use averaged value
                                last_group_threads = int(np.mean(group_threads_list))
                                last_group_ls = np.mean(group_ls_list)
                                logger.debug(f"        Current merged group contains {len(current_group)} points: threads={group_threads_list}, ls={group_ls_list}")
                                logger.debug(f"        Using averaged value for comparison: threads={last_group_threads}, ls={last_group_ls:.6e}")
                            
                            logger.debug(f"      [{idx+1}/{len(averaged_data)}] Checking: threads={threads} (ls={ls_avg:.6e}) vs last group threads={last_group_threads} (ls={last_group_ls:.6e})")
                            
                            threshold = last_group_ls * monotony_factor
                            if threads > last_group_threads and ls_avg < threshold:
                                # Violation: threads increased but weighted ls requirement failed
                                violation_found = True
                                ls_diff = ls_avg - last_group_ls
                                logger.debug(f"        ❌ VIOLATION: threads increased ({last_group_threads} -> {threads}) but ls fell below {monotony_factor:.2f}× requirement ({last_group_ls:.6e} -> {ls_avg:.6e}, diff={ls_diff:+.6e})")
                                logger.debug(f"        Merging threads={threads} (ls={ls_avg:.6e}) into current group")
                                # Add to current group (will be averaged later)
                                current_group.append({'threads': threads, 'ls_avg': ls_avg, 'count': int(row['count'])})
                                logger.debug(f"        Current group now contains {len(current_group)} points")
                            else:
                                # No violation with last group, finalize current group and start new one
                                group_adjusted = finalize_weighted_group(
                                    current_group,
                                    merged_data,
                                    monotony_factor,
                                    context_label="        ",
                                    boundary_averages=boundary_averages,
                                    min_threads=min_threads,
                                    max_threads=max_threads,
                                )
                                if group_adjusted:
                                    violation_found = True
                                
                                # Start new group
                                current_group = [{'threads': threads, 'ls_avg': ls_avg, 'count': int(row['count'])}]
                                logger.debug(f"        Starting new group: threads={threads}, ls={ls_avg:.6e}")
                    
                    # Finalize last group
                    if current_group:
                        group_adjusted = finalize_weighted_group(
                            current_group,
                            merged_data,
                            monotony_factor,
                            context_label="      [Final] ",
                            boundary_averages=boundary_averages,
                            min_threads=min_threads,
                            max_threads=max_threads,
                        )
                        if group_adjusted:
                            violation_found = True
                    
                    # Convert merged data to DataFrame and verify monotonicity
                    merged_df = pd.DataFrame(merged_data)
                    merged_df = merged_df.sort_values('threads').reset_index(drop=True)
                    
                    logger.debug(f"    After monotonicity enforcement:")
                    logger.debug(f"      Reduced from {len(averaged_data)} to {len(merged_df)} data points")
                    if violation_found:
                        logger.debug(f"      ✓ Monotonicity violations resolved by merging groups")
                    else:
                        logger.debug(f"      ✓ No violations found, data already monotonic")
                    
                    logger.debug(f"    Final merged data ({len(merged_df)} points):")
                    for idx, (_, row) in enumerate(merged_df.iterrows()):
                        threads_val = int(row['threads'])
                        ls_val = row['ls_avg']
                        count_val = int(row['count'])
                        logger.debug(f"      Point {idx+1}: threads={threads_val}, ls={ls_val:.6e} (from {count_val} measurement(s))")
                    thread_seq = [int(row['threads']) for _, row in merged_df.iterrows()]
                    ls_seq = []
                    for _, row in merged_df.iterrows():
                        ls_seq.append(f"{row['ls_avg']:.6e}")
                    logger.debug(f"    Thread sequence: {thread_seq}")
                    logger.debug(f"    LS sequence: {ls_seq}")
                    
                    # Verify monotonicity of final data
                    logger.debug(f"    Verifying final monotonicity:")
                    monotonic_ok = True
                    for i in range(len(merged_df) - 1):
                        t1 = int(merged_df.iloc[i]['threads'])
                        ls1 = merged_df.iloc[i]['ls_avg']
                        t2 = int(merged_df.iloc[i+1]['threads'])
                        ls2 = merged_df.iloc[i+1]['ls_avg']
                        required = ls1 * monotony_factor
                        if t2 > t1 and ls2 < required:
                            logger.debug(f"      ERROR: Still violates at threads={t2} (ls={ls2:.6e}) < {monotony_factor:.2f}× threads={t1} (ls={ls1:.6e})")
                            monotonic_ok = False
                        else:
                            logger.debug(f"      ✓ threads={t1} (ls={ls1:.6e}) -> threads={t2} (ls={ls2:.6e}) meets {monotony_factor:.2f}× requirement")
                    
                    if not monotonic_ok:
                        logger.debug(f"      ERROR: Monotonicity verification failed")
                        exit(1)
                    else:
                        logger.debug(f"      ✓ All monotonicity constraints satisfied")
                    
                    # After monotonicity enforcement, use the merged data
                    # The merged_df already has boundary points (from cross-policy avg) and middle points (from policy)
                    
                    combined_threads = []
                    combined_li = []
                    
                    logger.debug(f"  Step 5: Building final dataset from monotonicity-enforced data")
                    
                    # Add all points from merged_df (boundaries already have cross-policy averages)
                    for _, row in merged_df.iterrows():
                        threads_val = int(row['threads'])
                        ls_val = row['ls_avg']
                        combined_threads.append(threads_val)
                        combined_li.append(ls_val)
                    
                    # Sort by thread count
                    sort_idx = np.argsort(combined_threads)
                    combined_threads = [combined_threads[i] for i in sort_idx]
                    combined_li = [combined_li[i] for i in sort_idx]
                    
                    logger.debug(f"    Final dataset: {len(combined_threads)} points")
                    for threads, ls in zip(combined_threads, combined_li):
                        is_boundary = threads in [min_threads, max_threads]
                        tag = " (boundary - cross-policy)" if is_boundary else " (policy-only)"
                        logger.debug(f"      threads={threads}: ls={ls:.6e}{tag}")
                    
                    # Convert to numpy arrays
                    threads_array = np.array(combined_threads).reshape(-1, 1)
                    li_sync_array = np.array(combined_li)
                    
                    logger.debug(f"  Final combined dataset for regression: {len(threads_array)} data points for {alloc_policy}")
                    sorted_pairs_for_display = sorted(zip(combined_threads, combined_li), key=lambda x: x[0])
                    logger.debug(f"    Threads: {[t for t, _ in sorted_pairs_for_display]}")
                    logger.debug(f"    Li values: {[f'{li:.6e}' for _, li in sorted_pairs_for_display]}")
                    
                    # Verify monotonicity of final combined dataset
                    sorted_pairs = sorted(zip(combined_threads, combined_li), key=lambda x: x[0])
                    logger.debug(f"    Monotonicity check of final dataset:")
                    for i in range(len(sorted_pairs) - 1):
                        t1, ls1 = sorted_pairs[i]
                        t2, ls2 = sorted_pairs[i+1]
                        required = ls1 * monotony_factor
                        if t2 > t1 and ls2 < required:
                            logger.debug(f"      ERROR: Final dataset violates weighted monotonicity: threads={t2} (ls={ls2:.6e}) < {monotony_factor:.2f}× threads={t1} (ls={ls1:.6e})")
                            exit(1)
                        else:
                            logger.debug(f"      ✓ threads={t1} (ls={ls1:.6e}) -> threads={t2} (ls={ls2:.6e}) meets {monotony_factor:.2f}× requirement")
                    logger.debug(f"    ✓ Final dataset satisfies weighted monotonicity")
                    
                    # For empirical_monotony, store the adjusted data and skip polynomial fitting
                    if args.create_sync_level == "empirical_monotony":
                        # Store the adjusted monotonic data points
                        adjusted_monotonic_data[alloc_policy] = {
                            'threads': combined_threads.copy(),
                            'ls': combined_li.copy()
                        }
                        logger.debug(f"  Stored adjusted monotonic data for {alloc_policy}: {len(combined_threads)} points")
                        logger.debug(f"    Skipping polynomial fitting (using adjusted values directly)")
                        
                        # Store the adjusted data in a global variable for use in parse_hw_parameters
                        globals()[f'adjusted_monotonic_{alloc_policy}'] = {
                            'threads': combined_threads.copy(),
                            'ls': combined_li.copy()
                        }
                        continue  # Skip to next policy
                    
                    # Use Non-Negative Least Squares (NNLS) for polynomial fitting with positive coefficients
                    logger.debug(f"  Step 4: Polynomial regression using NNLS")
                    logger.debug(f"  Using Non-Negative Least Squares (NNLS) for {alloc_policy}...")
                    
                    from scipy.optimize import nnls
                    
                    # Create polynomial features matrix: [x^2, x, 1] for each thread count
                    threads_flat = threads_array.flatten()
                    A = np.column_stack([
                        threads_flat**2,  # x^2 terms
                        threads_flat,      # x terms  
                        np.ones(len(threads_flat))  # constant terms
                    ])
                    
                    logger.debug(f"    Regression matrix A shape: {A.shape}")
                    logger.debug(f"    Input data for regression:")
                    for i, (t, ls) in enumerate(zip(threads_flat, li_sync_array)):
                        logger.debug(f"      Point {i+1}: threads={t}, ls={ls:.6e}, [x²={t**2:.1f}, x={t:.1f}, 1]")
                    
                    # Use NNLS to find non-negative coefficients
                    logger.debug(f"    Running NNLS solver...")
                    coeffs, residual = nnls(A, li_sync_array)
                    a, b, c = coeffs[0], coeffs[1], coeffs[2]
                    
                    logger.debug(f"    NNLS completed:")
                    logger.debug(f"      Coefficients: a={a:.6e}, b={b:.6e}, c={c:.6e}")
                    logger.debug(f"      Residual (sum of squared errors): {residual:.6e}")
                    logger.debug(f"      Model: ls(threads) = {a:.6e} * threads² + {b:.6e} * threads + {c:.6e}")
    
                    alignment_applied = False
                    if boundary_averages:
                        a, b, c, alignment_applied = enforce_polynomial_boundary_alignment(
                            a, b, c, min_threads, max_threads, boundary_averages, alloc_policy
                        )
                        if alignment_applied:
                            logger.debug(f"      ✓ Enforced cross-policy boundary equality for {alloc_policy}")
                
                    # Check if y(threads=1) > 0
                    y_at_one = a * (1**2) + b * 1 + c
                    if y_at_one <= 0:
                        logger.debug(f"  WARNING: y(threads=1) = {y_at_one:.6e} <= 0, adjusting intercept")
                        # Ensure y(1) > 0 by adjusting the intercept
                        c = max(c, 1e-10 - a - b)  # Ensure a + b + c > 0
                        logger.debug(f"  Adjusted coefficients: a={a:.6e}, b={b:.6e}, c={c:.6e}")
                        y_at_one = a * (1**2) + b * 1 + c
                        logger.debug(f"  New y(threads=1) = {y_at_one:.6e}")
                    
                    # Calculate fit quality
                    predicted_values = a * (threads_array.flatten()**2) + b * threads_array.flatten() + c
                    residuals = li_sync_array - predicted_values
                    r2_score = 1 - np.sum(residuals**2) / np.sum((li_sync_array - np.mean(li_sync_array))**2)
                    rmse = np.sqrt(np.mean(residuals**2))
                    mae = np.mean(np.abs(residuals))
                    
                    logger.debug(f"NNLS polynomial fit for {alloc_policy} policy (endpoint-constrained):")
                    logger.debug(f"  Li_sync = {a:.6e} * threads² + {b:.6e} * threads + {c:.6e}")
                    logger.debug(f"  R² = {r2_score:.4f}")
                    logger.debug(f"  RMSE = {rmse:.6e}")
                    logger.debug(f"  MAE = {mae:.6e}")
                    logger.debug(f"  Number of data points = {len(threads_array)}")
                    logger.debug(f"  Thread range: {threads_array.min()} - {threads_array.max()}")
                    logger.debug(f"  Li range: {li_sync_array.min():.6e} - {li_sync_array.max():.6e}")
                    
                    # Debug: show original vs predicted values
                    logger.debug(f"  Data comparison:")
                    for i, (thread_val, orig_val, pred_val) in enumerate(zip(threads_array.flatten(), li_sync_array, predicted_values)):
                        logger.debug(f"    threads={thread_val:2d}: orig={orig_val:.6e}, pred={pred_val:.6e}, diff={abs(orig_val-pred_val):.6e}")
                
                    # Validate that all predicted values are positive
                    logger.debug(f"  Validating positive predictions:")
                    negative_count = 0
                    for thread_val in threads_array.flatten():
                        pred_val = a * (thread_val**2) + b * thread_val + c
                        if pred_val <= 0:
                            negative_count += 1
                            logger.debug(f"    ERROR: threads={thread_val} -> Li_sync={pred_val:.6e} (non-positive)")
                    
                    if negative_count == 0:
                        logger.debug(f"    ✓ All predicted values are positive")
                    else:
                        logger.debug(f"    ERROR: {negative_count} predicted values are non-positive")
                        logger.debug(f"    This violates the constraint that latency must be positive.")
                        logger.debug(f"    Please check your polynomial coefficients or data quality.")
                        exit(1)
                    
                    # Verify endpoint equality
                    y_at_min = a * (min_threads**2) + b * min_threads + c
                    y_at_max = a * (max_threads**2) + b * max_threads + c
                    logger.debug(f"  Endpoint verification:")
                    logger.debug(f"    y({min_threads}) = {y_at_min:.6e}")
                    logger.debug(f"    y({max_threads}) = {y_at_max:.6e}")
                    
                    # Log polynomial fitting summary
                    poly_details = {
                        'Data Points': len(threads_array),
                        'R² Score': f"{r2_score:.4f}",
                        'RMSE': f"{rmse:.6e}",
                        'Coefficients': f"a={a:.6e}, b={b:.6e}, c={c:.6e}",
                        'Endpoint Constrained': 'Yes'
                    }
                    log_summary(f"Polynomial Fit {alloc_policy}", "SUCCESS", poly_details)
                
                    # Validate simple constraints
                    logger.debug(f"  Validating constraints:")
                    
                    # Check non-negative coefficients
                    coeffs_ok = (a >= 0) and (b >= 0) and (c >= 0)
                    logger.debug(f"    Non-negative coefficients: a={a:.6e}, b={b:.6e}, c={c:.6e} - {'✓' if coeffs_ok else '✗'}")
                    
                    # Check y(threads=1) > 0
                    y_at_one = a * (1**2) + b * 1 + c
                    y_at_one_ok = y_at_one > 0
                    logger.debug(f"    y(threads=1) = {y_at_one:.6e} - {'✓' if y_at_one_ok else '✗'}")
                    
                    if not (coeffs_ok and y_at_one_ok):
                        logger.debug(f"    ERROR: Constraints violated")
                        logger.debug(f"    Please check your polynomial coefficients or data quality.")
                        exit(1)
                    else:
                        logger.debug(f"    ✓ All constraints satisfied")
                    
                    # Store the regression coefficients and characteristics for later use
                    globals()[f'lr_coeffs_{alloc_policy}'] = {
                        'a': a, 'b': b, 'c': c,
                        'model_type': 'nnls_endpoint_constrained', 'model_name': 'NNLS with Endpoint Constraints',
                        'r2': r2_score, 'rmse': rmse, 'mae': mae, 'n_points': len(threads_array)
                    }
        
            # Display hardware analysis section
            logger.info("")
            logger.info("="*60)
            logger.info("HARDWARE ANALYSIS")
            logger.info("="*60)
            logger.info(f"Processing {len(alloc_policies)} allocation policies")
            logger.info(f"Generating hardware parameters and models...")
            logger.info("="*60)
            
            # Enforce parameter monotonicity on raw benchmark data before subset generation
            logger.info("")
            logger.info("="*60)
            logger.info("PARAMETER MONOTONICITY ENFORCEMENT")
            logger.info("="*60)
            logger.info("Enforcing monotonicity constraints on raw benchmark data...")
            logger.info("="*60)
            
            all_combined_df = enforce_parameter_monotonicity(all_combined_df, kernel_aggregator=args.kernel_aggregator, keep_levels=keep_levels)
            
            # Generate hardware parameters for all policies combined
            logger.debug(f"\n{'='*60}")
            logger.debug(f"Generating hardware parameters for all policies: {alloc_policies}")
            logger.debug(f"{'='*60}")
            
            # Read r_scalar and r_vec from topology file
            from hw_characterization.scripts.parser import get_topology_file_path as _get_topo_path
            from hw_characterization.scripts.parser import parse_pahis_log_to_df as _parse_topo
            _topo_path = _get_topo_path(system_string_local)
            _, _, _, _, _r_scalar, _r_vec = _parse_topo(_topo_path)
            
            hw_params = parse_hw_parameters(
                system_string_local,
                all_combined_df,
                keep_levels=keep_levels,
                aggregate_method=args.kernel_aggregator,
                alloc_policy=",".join(alloc_policies),  # Pass all policies as comma-separated string
                all_alloc_policies=alloc_policies,
                generate_subsets=args.generate_subsets,
                create_sync_level=args.create_sync_level,
                all_combined_df=all_combined_df,
                r_scalar=_r_scalar if _r_scalar is not None else 0.0,
                r_vec=_r_vec if _r_vec is not None else 0.0,
            )
            
            # Track generated header file
            system_id = system_string_local.replace('-', '_').replace(' ', '_').lower()
            header_file = f"{CPP_MODEL_DIR}/hw_params_{system_id}.hpp"
            track_file(header_file)
            
            hw_details = {
                'Models Generated': len(hw_params),
                'Policies': len(alloc_policies),
                'Header File': header_file
            }
            log_summary("Hardware Parameters", "SUCCESS", hw_details, [header_file])
            
            # Add detailed subset summary to log
            logger.debug(f"Hardware parameter generation completed:")
            logger.debug(f"  Total models: {len(hw_params)}")
            logger.debug(f"  Policies: {alloc_policies}")
            logger.debug(f"  Thread configurations: {len(threads_list)}")
            logger.debug(f"  Total expected models: {len(alloc_policies) * len(threads_list) * 16}")  # 16 is the number of subsets for 5 levels
    
        # Generate hardware parameter summary plots from the header file
        logger.info("")
        logger.info("="*60)
        logger.info("HARDWARE PARAMETER PLOTS")
        logger.info("="*60)
        try:
            from analysis.plotters.hw_params_header_plotter import (
                parse_hpp_header, resolve_configs, plot_configuration
            )
            system_id_plot = system_string_local.replace('-', '_').replace(' ', '_').lower()
            header_path_plot = f"{CPP_MODEL_DIR}/hw_params_{system_id_plot}.hpp"

            if os.path.exists(header_path_plot):
                hw_config, hw_sys_name = parse_hpp_header(header_path_plot)

                # Select thread/policy combos:
                #   1-close, max-close, max/2-close+spread, max/4-close+spread
                avail_threads = sorted(hw_config.threads_options_num)
                max_t = avail_threads[-1]
                half_t = min(avail_threads, key=lambda t: abs(t - max_t / 2))
                quarter_t = min(avail_threads, key=lambda t: abs(t - max_t / 4))

                max_bh = max(hw_config.hw_model_level_bithash)
                plots_dir = f"models/{system_string_local}"
                plots_path = Path(plots_dir)
                plots_path.mkdir(parents=True, exist_ok=True)

                # Remove previous plots in that directory (keep scope tight)
                removed = 0
                for pat in ("hw_params_*.png", "hw_params_*.pdf"):
                    for f in plots_path.glob(pat):
                        try:
                            f.unlink()
                            removed += 1
                        except Exception:
                            # Non-fatal: continue generating new plots
                            pass
                if removed:
                    logger.info(f"Removed {removed} previous plot(s) from: {plots_dir}/")

                # Build explicit (threads, policy) pairs
                plot_specs = [
                    (1, "close"),
                    (max_t, "close"),
                    (half_t, "close"), (half_t, "spread"),
                    (quarter_t, "close"), (quarter_t, "spread"),
                ]
                # Deduplicate while preserving order
                seen = set()
                unique_specs = []
                for spec in plot_specs:
                    if spec not in seen:
                        seen.add(spec)
                        unique_specs.append(spec)

                resolved = []
                for t, p in unique_specs:
                    resolved.extend(resolve_configs(
                        hw_config, threads=[t], policies=[p], bithashes=[max_bh],
                    ))
                logger.info(f"Generating {len(resolved)} hw-parameter plots ...")
                logger.info(f"  Configs: {unique_specs}")
                logger.info(f"  Bithash: {max_bh}")
                logger.info(f"  Output: {plots_dir}/")
                for rc in resolved:
                    plot_path = plot_configuration(rc, hw_sys_name, plots_dir,
                                                   config=hw_config)
                    track_file(plot_path)
                    logger.info(f"  Saved: {plot_path}")
                log_summary("HW Parameter Plots", "SUCCESS",
                            {'Plots': len(resolved), 'Directory': plots_dir})
            else:
                logger.warning(f"Header file not found at {header_path_plot}, skipping plots")
                log_summary("HW Parameter Plots", "WARNING", {'Reason': 'Header not found'})
        except Exception as e:
            logger.warning(f"HW parameter plot generation failed: {e}")
            log_summary("HW Parameter Plots", "WARNING", {'Error': str(e)})
        logger.info("="*60)

        # Display model compilation section
        logger.info("")
        logger.info("="*60)
        logger.info("MODEL COMPILATION")
        logger.info("="*60)
        logger.info("Compiling C++ models with generated hardware parameters...")
        logger.info("="*60)
        
        # Compile models with generated hardware parameters
        logger.debug("\n" + "="*60)
        logger.debug("Compiling models with generated hardware parameters...")
        logger.debug("="*60)
        
        compile_success = compile_models_with_hw_params(system_string_local)
        
        if compile_success:
            # Track compiled executables
            models_build_dir = Path(BUILD_DIR) / "model_build" / "bin"
            if models_build_dir.exists():
                for exe in models_build_dir.glob("*"):
                    if exe.is_file() and exe.stat().st_mode & 0o111:
                        track_file(str(exe))
            
            compile_details = {
                'Status': 'Success',
                'Executables': len(list(models_build_dir.glob("*"))) if models_build_dir.exists() else 0
            }
            log_summary("Model Compilation", "SUCCESS", compile_details)
            
            # Display unit testing section
            logger.info("")
            logger.info("="*60)
            logger.info("UNIT TESTING")
            logger.info("="*60)
            logger.info("Running unit tests for model compilation and execution...")
            logger.info("="*60)
            
            # Run unit tests for model compilation and execution
            logger.debug("\n" + "="*60)
            logger.debug("Running unit tests for model compilation and execution...")
            logger.debug("="*60)
            
            test_success = run_model_unit_tests(system_string_local)
            test_details = {'Status': 'Success' if test_success else 'Failed'}
            log_summary("Unit Tests", "SUCCESS" if test_success else "WARNING", test_details)
        else:
            log_summary("Model Compilation", "ERROR", {'Status': 'Failed'})
            logger.info("⚠️  Model compilation failed, skipping unit tests")
        
        # Generate final summary
        logger.info("")
        logger.info("="*60)
        logger.info("EXECUTION SUMMARY")
        logger.info("="*60)
        
        # Show file generation summary
        files_summary = get_generated_files_summary()
        if files_summary:
            logger.info("Generated Files:")
            for file_type, files in files_summary.items():
                logger.info(f"  {file_type}: {len(files)} files")
                for file_path in files[:3]:  # Show first 3 files
                    logger.info(f"    - {file_path}")
                if len(files) > 3:
                    logger.info(f"    ... and {len(files) - 3} more")
        
        # Show step summary
        logger.info("")
        logger.info("Steps Completed:")
        for step in step_summaries:
            status_symbol = "✓" if step['status'] == "SUCCESS" else "✗" if step['status'] == "ERROR" else "⚠"
            logger.info(f"  {status_symbol} {step['step']}")
        
        logger.info("")
        logger.info(f"Detailed log: {log_file}")
        logger.info("="*60)
        
    except Exception as e:
        logger.error(f"Fatal error: {e}")
        log_summary("Execution", "ERROR", {'Error': str(e)})
        raise

if __name__ == "__main__":
    _main()