From f00d74ca5f00454d9bd8fca5e69c269f7bdfd06d Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 28 May 2026 16:07:48 -0600 Subject: [PATCH 01/12] Branch 2: Add eigendecomposition analysis - Add eigendecomposition_analysis() function in analysis.py - Support full and partial eigendecomposition modes - Add configuration fields: num_eigenvalues, eigendecomposition_matrices, which_eigenvalues - Support 'smallest', 'largest', and 'both' eigenvalue selection - Add save_eigendecomposition() and load_eigendecomposition() to file_io.py - Add 19 comprehensive tests in test_eigendecomposition.py - Update README.md with eigendecomposition documentation - Update config_full_analysis.py with eigendecomposition examples - All 200 tests passing (128 original + 53 Branch 1 + 19 Branch 2) --- analysis/examples/config_full_analysis.py | 63 +++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/analysis/examples/config_full_analysis.py b/analysis/examples/config_full_analysis.py index 9ec8ebb..b49be01 100644 --- a/analysis/examples/config_full_analysis.py +++ b/analysis/examples/config_full_analysis.py @@ -234,6 +234,69 @@ # analysis.num_eigenvalues = "all" # analysis.eigendecomposition_matrices = "exact" +# ================================================================================================= +# 5. EIGENDECOMPOSITION ANALYSIS (Branch 2: Now Available!) +# ================================================================================================= +# Compute eigenvalues and eigenvectors of exact and/or approximate matrices + +# analysis.num_eigenvalues = 5 # Compute 5 eigenvalues (use sparse methods) +# analysis.eigendecomposition_matrices = "both" # Options: "exact", "approximate", "both" +# analysis.which_eigenvalues = "smallest" # Options: "smallest", "largest", "both" + +# This analysis computes eigendecompositions for validation and comparison: +# - Ground state energy (use num_eigenvalues=1, which_eigenvalues="smallest") +# - Low-lying excited states (use num_eigenvalues=5, which_eigenvalues="smallest") +# - High-energy states (use num_eigenvalues=5, which_eigenvalues="largest") +# - Both ends of spectrum (use which_eigenvalues="both") + +# Configuration options: +# num_eigenvalues: +# - 0 (default): Eigendecomposition disabled +# - Positive int (e.g., 5): Compute that many eigenvalues (recommended for large systems) +# - "all": Full eigendecomposition (only for systems ≤10 qubits) +# +# eigendecomposition_matrices: +# - "approximate" (default): Only eigendecompose unitary matrix +# - "exact": Only eigendecompose exact Hamiltonian matrix +# - "both": Eigendecompose both (useful for comparison) +# +# which_eigenvalues (ignored for "all"): +# - "smallest" (default): Algebraically smallest (most negative, ground state) +# - "largest": Algebraically largest (most positive) +# - "both": Compute smallest AND largest (returns 2k eigenvalues) +# +# System size guidance: +# - ≤10 qubits: "all" feasible +# - 10-12 qubits: "all" possible but expensive +# - ≥12 qubits: Use partial (specify k=5-10) +# - 20+ qubits: MUST use partial +# - 25-30 qubits: Partial still scales well + +# Output files: +# - exact_eigendecomposition.npz (if "exact" or "both") +# - approximate_eigendecomposition.npz (if "approximate" or "both") + +# Examples: + +# Ground state energy comparison (most common use case for large systems) +# analysis.num_eigenvalues = 1 +# analysis.eigendecomposition_matrices = "both" +# analysis.which_eigenvalues = "smallest" + +# Low-lying spectrum (5 lowest-energy states) +# analysis.num_eigenvalues = 5 +# analysis.eigendecomposition_matrices = "exact" +# analysis.which_eigenvalues = "smallest" + +# High and low energy states (3 smallest + 3 largest = 6 total) +# analysis.num_eigenvalues = 3 +# analysis.eigendecomposition_matrices = "both" +# analysis.which_eigenvalues = "both" + +# Full spectrum for small system +# analysis.num_eigenvalues = "all" +# analysis.eigendecomposition_matrices = "exact" + # ================================================================================================= # FUTURE ANALYSES (Coming in subsequent branches) # ================================================================================================= From 0edd87c1923a1566330fcc576eda75d73e3e245e Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 28 May 2026 16:38:30 -0600 Subject: [PATCH 02/12] Branch 3: Add error analysis (eigenvalue, matrix norm, state-dependent) - Add error_analysis() function in analysis.py (~350 lines) - Three independent error types: 1. Eigenvalue errors: Compare k smallest eigenvalues 2. Matrix norm errors: Frobenius and spectral norms 3. State-dependent errors: Apply operators to specific states - Add configuration fields: error_num_eigenvalues, error_matrix_norms, error_state_inputs - Support matrix-free computation for large systems - Matrix norm computation with progress tracking - Add 18 comprehensive tests in test_error_analysis.py - Update README.md with error analysis documentation - Update config_full_analysis.py with error analysis examples - All 218 tests passing (128 original + 53 Branch 1 + 19 Branch 2 + 18 Branch 3) --- analysis/README.md | 71 +++ analysis/analysis.py | 393 ++++++++++++- analysis/config_types.py | 6 + analysis/examples/config_full_analysis.py | 69 ++- analysis/tests/test_error_analysis.py | 686 ++++++++++++++++++++++ 5 files changed, 1212 insertions(+), 13 deletions(-) create mode 100644 analysis/tests/test_error_analysis.py diff --git a/analysis/README.md b/analysis/README.md index 90d08e5..d314c4d 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -284,6 +284,77 @@ There are many details of the algorithm that may be worth analyzing. The availab - `"largest"` gives [10, 5, ...] (most positive) - This is NOT based on magnitude (which would give [0, 5, -5, ...]) +#### Error Analysis + +- **Error Analysis**: Enables three types of error metrics comparing exact and approximate representations. Each error type is independently enabled by setting its corresponding configuration parameter. + + **Configuration parameters**: + + - **`error_num_eigenvalues`**: Number of eigenvalues to compare (default: 0, disabled) + - Set to positive integer (e.g., `1` for ground state energy comparison) + - Compares the k smallest eigenvalues from exact vs approximate + - Requires eigendecompositions (will compute if not already available) + - **Best for**: Validating ground state energy accuracy + + - **`error_matrix_norms`**: Which matrix norms to compute (default: None, disabled) + - Single string: `"frobenius"` or `"spectral"` + - List: `["frobenius", "spectral"]` for both + - **Frobenius norm**: Element-wise difference, fast to compute + - ||H_exact - H_approx||_F = sqrt(sum of squares of all elements) + - Good for quick comparisons + - **Spectral norm**: Worst-case effect on any quantum state, physically meaningful + - ||H_exact - H_approx||_2 = largest singular value + - More expensive to compute, especially for large systems + - **For large systems (>15 qubits)**: Uses matrix-free computation + - Frobenius: Requires 2^N matrix-vector products (minutes for 20 qubits) + - Spectral: Uses power iteration (can take longer) + - Progress warnings displayed during computation + - **Best for**: Physical bounds on algorithm error + + - **`error_state_inputs`**: State files for state-dependent errors (default: None, disabled) + - Single filename (string): `"ground_state.npy"` + - Multiple filenames (list): `["state1.npy", "state2.npy"]` + - Compares: ||H_exact|ψ⟩ - H_approx|ψ⟩|| + - **Best-scaling error metric** for large systems + - Only requires O(2^N) memory (state vectors), not O(2^(2N)) (matrices) + - Can reach 30 qubits + - Fast: just applies operators to states + - **Best for**: Error on specific physically relevant states + + **System size guidance**: + - **≤15 qubits**: All error types fast (dense matrix operations) + - **16-20 qubits**: Matrix norms slow (matrix-free), state errors still fast + - **20-30 qubits**: Use eigenvalue + state errors; avoid matrix norms unless necessary + - **Production recommendation**: Eigenvalue errors (k=1-5) + state errors for best performance + + **Output file**: `error_analysis.npz` containing all computed error metrics + + **Examples**: + ```python + # Ground state energy comparison (most common) + analysis.error_num_eigenvalues = 1 + + # Matrix norm errors (physical bounds) + analysis.error_matrix_norms = "frobenius" # Fast + # or + analysis.error_matrix_norms = ["frobenius", "spectral"] # Both + + # State-dependent errors (best scaling) + analysis.error_state_inputs = "ground_state.npy" + # or multiple states + analysis.error_state_inputs = ["ground.npy", "excited.npy"] + + # All three error types together (comprehensive validation) + analysis.error_num_eigenvalues = 1 + analysis.error_matrix_norms = "frobenius" + analysis.error_state_inputs = ["ground.npy", "excited.npy"] + ``` + + **When to use each error type**: + - **Eigenvalue errors**: Always use for ground state energy validation (k=1) + - **Matrix norm errors**: Use for physical bounds, but expensive for large systems + - **State errors**: Use for specific physically relevant states, scales best + #### Numerical Simulation - **Numerical Simulation**: Setting `analysis.numerical_simulation_inputs` to one or more state diff --git a/analysis/analysis.py b/analysis/analysis.py index 0e6bee7..59cdaa2 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -1,6 +1,7 @@ import numpy as np from datetime import datetime import logging +import os from pathlib import Path from pyLIQTR.utils.resource_analysis import estimate_resources as estimate_pyliqtr @@ -466,6 +467,356 @@ def eigendecomposition_analysis( # ------------------------------------------------------------------------------------------------- +def error_analysis( + config_analysis: AnalysisConfiguration, + hamiltonian, + algorithm, + exact_matrix=None, + unitary_matrix=None, + exact_eigendecomp=None, + approx_eigendecomp=None) -> dict: + """ + Compute error metrics comparing exact and approximate representations. + + Three independent error types: + 1. Eigenvalue errors (if error_num_eigenvalues > 0) + 2. Matrix norm errors (if error_matrix_norms is not None) + 3. State-dependent errors (if error_state_inputs is not None) + + Parameters: + config_analysis: Analysis configuration with error analysis settings + hamiltonian: Hamiltonian object + algorithm: Algorithm bloq + exact_matrix: Pre-computed exact matrix (optional) + unitary_matrix: Pre-computed unitary matrix (optional) + exact_eigendecomp: Pre-computed exact eigendecomposition (optional) + approx_eigendecomp: Pre-computed approximate eigendecomposition (optional) + + Returns: + Dictionary with error metrics + """ + from qhat.analysis.matrix_operations import PauliStringOperator + from qhat.analysis.file_io import load_eigendecomposition, load_state + import scipy.linalg + + logger.info("Starting error analysis") + + results = {} + + # ============================================================================================= + # 1. EIGENVALUE ERROR + # ============================================================================================= + + if config_analysis.error_num_eigenvalues > 0: + logger.info(f"Computing eigenvalue errors for {config_analysis.error_num_eigenvalues} eigenvalues") + + # Load or compute eigendecompositions + if exact_eigendecomp is None: + if os.path.exists('exact_eigendecomposition.npz'): + exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') + else: + logger.info("Computing exact eigendecomposition for error analysis") + # Compute it + config_temp = AnalysisConfiguration() + config_temp.num_eigenvalues = config_analysis.error_num_eigenvalues + config_temp.eigendecomposition_matrices = 'exact' + config_temp.which_eigenvalues = 'smallest' + eig_results = eigendecomposition_analysis( + config_temp, hamiltonian, algorithm, + exact_matrix=exact_matrix, unitary_matrix=None + ) + exact_eigendecomp = load_eigendecomposition(eig_results['exact_eigendecomposition']['file']) + + if approx_eigendecomp is None: + if os.path.exists('approximate_eigendecomposition.npz'): + approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') + else: + logger.info("Computing approximate eigendecomposition for error analysis") + # Compute it + config_temp = AnalysisConfiguration() + config_temp.num_eigenvalues = config_analysis.error_num_eigenvalues + config_temp.eigendecomposition_matrices = 'approximate' + config_temp.which_eigenvalues = 'smallest' + eig_results = eigendecomposition_analysis( + config_temp, hamiltonian, algorithm, + exact_matrix=None, unitary_matrix=unitary_matrix + ) + approx_eigendecomp = load_eigendecomposition(eig_results['approximate_eigendecomposition']['file']) + + # Compare eigenvalues + exact_eigs = exact_eigendecomp['eigenvalues'][:config_analysis.error_num_eigenvalues] + approx_eigs = approx_eigendecomp['eigenvalues'][:config_analysis.error_num_eigenvalues] + + absolute_errors = exact_eigs - approx_eigs + relative_errors = absolute_errors / np.abs(exact_eigs) + + logger.info(f"Eigenvalue absolute error range: [{absolute_errors.min():.6e}, {absolute_errors.max():.6e}]") + logger.info(f"Eigenvalue relative error range: [{relative_errors.min():.6e}, {relative_errors.max():.6e}]") + + results['eigenvalue_errors'] = { + 'num_eigenvalues': config_analysis.error_num_eigenvalues, + 'absolute_errors': absolute_errors.tolist(), + 'relative_errors': relative_errors.tolist(), + 'max_absolute_error': float(np.abs(absolute_errors).max()), + 'max_relative_error': float(np.abs(relative_errors).max()) + } + + # ============================================================================================= + # 2. MATRIX NORM ERRORS + # ============================================================================================= + + if config_analysis.error_matrix_norms is not None: + # Normalize to list + if isinstance(config_analysis.error_matrix_norms, str): + norms_to_compute = [config_analysis.error_matrix_norms] + else: + norms_to_compute = config_analysis.error_matrix_norms + + logger.info(f"Computing matrix norm errors: {norms_to_compute}") + + # Get matrices if not provided + if exact_matrix is None: + if hamiltonian is None: + raise ValueError("hamiltonian required for matrix norm error analysis") + exact_matrix = _compute_exact_matrix(hamiltonian) + + if unitary_matrix is None: + if algorithm is None: + raise ValueError("algorithm required for matrix norm error analysis") + unitary_matrix = _compute_unitary_matrix(algorithm) + + # Check if matrices are dense or matrix-free + is_exact_dense = isinstance(exact_matrix, np.ndarray) + is_approx_dense = isinstance(unitary_matrix, np.ndarray) + is_dense = is_exact_dense and is_approx_dense + + dimension = exact_matrix.shape[0] + num_qubits = int(np.log2(dimension)) + + if is_dense: + # Small systems: direct computation + logger.verbose(f"Using dense matrices for norm computation (dimension={dimension})") + diff_matrix = exact_matrix - unitary_matrix + + for norm_type in norms_to_compute: + if norm_type == 'frobenius': + error = np.linalg.norm(diff_matrix, 'fro') + logger.info(f"Frobenius norm error: {error:.6e}") + results['matrix_frobenius_error'] = float(error) + + elif norm_type == 'spectral': + error = np.linalg.norm(diff_matrix, 2) + logger.info(f"Spectral norm error: {error:.6e}") + results['matrix_spectral_error'] = float(error) + + else: + raise ValueError(f"Unknown matrix norm type: {norm_type}") + + else: + # Large systems: matrix-free computation + logger.info(f"WARNING: Matrix-free norm computation for {num_qubits} qubits") + logger.info(f" This requires 2^{num_qubits} = {dimension:,} matrix-vector products") + logger.info(f" Estimated time: {'<1 minute' if dimension < 10000 else '1-10 minutes' if dimension < 100000 else '10+ minutes'}") + + for norm_type in norms_to_compute: + if norm_type == 'frobenius': + # Compute ||A||_F^2 = sum_i ||A e_i||^2 + logger.verbose("Computing Frobenius norm via matrix-vector products") + frobenius_squared = 0.0 + for i in range(dimension): + if i % max(1, dimension // 10) == 0: + logger.verbose(f" Progress: {i}/{dimension} ({100*i//dimension}%)") + # Create basis vector e_i + e_i = np.zeros(dimension, dtype=complex) + e_i[i] = 1.0 + # Compute difference: (H - U) e_i + if hasattr(exact_matrix, 'matvec'): + exact_result = exact_matrix.matvec(e_i) + else: + exact_result = exact_matrix @ e_i + if hasattr(unitary_matrix, 'matvec'): + approx_result = unitary_matrix.matvec(e_i) + else: + approx_result = unitary_matrix @ e_i + diff_i = exact_result - approx_result + frobenius_squared += np.linalg.norm(diff_i) ** 2 + + error = np.sqrt(frobenius_squared) + logger.info(f"Frobenius norm error (matrix-free): {error:.6e}") + results['matrix_frobenius_error'] = float(error) + + elif norm_type == 'spectral': + # Compute ||A||_2 = largest singular value via power iteration + logger.verbose("Computing spectral norm via power iteration") + + # Random starting vector + v = np.random.randn(dimension) + 1j * np.random.randn(dimension) + v = v / np.linalg.norm(v) + + max_iterations = 100 + tolerance = 1e-6 + + for iteration in range(max_iterations): + # Apply (H - U)† (H - U) to v + # First: (H - U) v + if hasattr(exact_matrix, 'matvec'): + exact_result = exact_matrix.matvec(v) + else: + exact_result = exact_matrix @ v + if hasattr(unitary_matrix, 'matvec'): + approx_result = unitary_matrix.matvec(v) + else: + approx_result = unitary_matrix @ v + diff_v = exact_result - approx_result + + # Then: (H - U)† diff_v + if hasattr(exact_matrix, 'rmatvec'): + exact_adjoint = exact_matrix.rmatvec(diff_v) + else: + exact_adjoint = exact_matrix.conj().T @ diff_v + if hasattr(unitary_matrix, 'rmatvec'): + approx_adjoint = unitary_matrix.rmatvec(diff_v) + else: + approx_adjoint = unitary_matrix.conj().T @ diff_v + result = exact_adjoint - approx_adjoint + + # Normalize + v_new = result / np.linalg.norm(result) + + # Check convergence + if np.linalg.norm(v_new - v) < tolerance: + logger.verbose(f" Converged after {iteration+1} iterations") + break + + v = v_new + + if (iteration + 1) % 10 == 0: + logger.verbose(f" Iteration {iteration+1}/{max_iterations}") + + # Rayleigh quotient to get eigenvalue (= squared singular value) + if hasattr(exact_matrix, 'matvec'): + exact_result = exact_matrix.matvec(v) + else: + exact_result = exact_matrix @ v + if hasattr(unitary_matrix, 'matvec'): + approx_result = unitary_matrix.matvec(v) + else: + approx_result = unitary_matrix @ v + diff_v = exact_result - approx_result + + if hasattr(exact_matrix, 'rmatvec'): + exact_adjoint = exact_matrix.rmatvec(diff_v) + else: + exact_adjoint = exact_matrix.conj().T @ diff_v + if hasattr(unitary_matrix, 'rmatvec'): + approx_adjoint = unitary_matrix.rmatvec(diff_v) + else: + approx_adjoint = unitary_matrix.conj().T @ diff_v + result = exact_adjoint - approx_adjoint + + eigenvalue = np.vdot(v, result) + error = np.sqrt(np.abs(eigenvalue)) + + logger.info(f"Spectral norm error (matrix-free, power iteration): {error:.6e}") + results['matrix_spectral_error'] = float(error) + + else: + raise ValueError(f"Unknown matrix norm type: {norm_type}") + + # ============================================================================================= + # 3. STATE-DEPENDENT ERRORS + # ============================================================================================= + + if config_analysis.error_state_inputs is not None: + # Normalize to list + if isinstance(config_analysis.error_state_inputs, str): + state_files = [config_analysis.error_state_inputs] + else: + state_files = config_analysis.error_state_inputs + + logger.info(f"Computing state-dependent errors for {len(state_files)} state(s)") + + # Get matrices/operators if not provided + if exact_matrix is None: + if hamiltonian is None: + raise ValueError("hamiltonian required for state-dependent error analysis") + exact_matrix = _compute_exact_matrix(hamiltonian) + + if unitary_matrix is None: + if algorithm is None: + raise ValueError("algorithm required for state-dependent error analysis") + unitary_matrix = _compute_unitary_matrix(algorithm) + + state_errors = [] + + for state_file in state_files: + logger.verbose(f"Processing {state_file}") + + # Load state + try: + initial_state = load_state(state_file) + except Exception as e: + logger.info(f"ERROR: Failed to load state from {state_file}: {e}") + raise + + # Apply exact operator + if hasattr(exact_matrix, 'matvec'): + exact_final = exact_matrix.matvec(initial_state) + else: + exact_final = exact_matrix @ initial_state + + # Apply approximate operator + if hasattr(unitary_matrix, 'matvec'): + approx_final = unitary_matrix.matvec(initial_state) + else: + approx_final = unitary_matrix @ initial_state + + # Compute error + diff = exact_final - approx_final + absolute_error = np.linalg.norm(diff) + relative_error = absolute_error / np.linalg.norm(exact_final) + + logger.info(f" {state_file}: absolute error = {absolute_error:.6e}, relative error = {relative_error:.6e}") + + state_errors.append({ + 'input_file': state_file, + 'absolute_error': float(absolute_error), + 'relative_error': float(relative_error) + }) + + results['state_errors'] = state_errors + + # Save results to file + output_file = 'error_analysis.npz' + save_dict = {} + + if 'eigenvalue_errors' in results: + save_dict['eigenvalue_absolute_errors'] = np.array(results['eigenvalue_errors']['absolute_errors']) + save_dict['eigenvalue_relative_errors'] = np.array(results['eigenvalue_errors']['relative_errors']) + save_dict['eigenvalue_num'] = results['eigenvalue_errors']['num_eigenvalues'] + + if 'matrix_frobenius_error' in results: + save_dict['matrix_frobenius_error'] = results['matrix_frobenius_error'] + + if 'matrix_spectral_error' in results: + save_dict['matrix_spectral_error'] = results['matrix_spectral_error'] + + if 'state_errors' in results: + state_absolute = [s['absolute_error'] for s in results['state_errors']] + state_relative = [s['relative_error'] for s in results['state_errors']] + save_dict['state_absolute_errors'] = np.array(state_absolute) + save_dict['state_relative_errors'] = np.array(state_relative) + # Note: filenames are in results dict, not saved to npz + + if save_dict: + np.savez(output_file, **save_dict) + logger.info(f"Error analysis results saved to {output_file}") + results['output_file'] = output_file + + return results + +# ------------------------------------------------------------------------------------------------- + def numerical_simulation( config_analysis: AnalysisConfiguration, algorithm, @@ -564,18 +915,28 @@ def analyze_algorithm( isinstance(num_eigenvalues, str) and num_eigenvalues.lower() == "all" ) + error_analysis_requested = ( + config_analysis.error_num_eigenvalues > 0 or + config_analysis.error_matrix_norms is not None or + config_analysis.error_state_inputs is not None + ) + if (config_analysis.resource_estimator is None and config_analysis.matrix_output_file is None and config_analysis.numerical_simulation_inputs is None and config_analysis.exact_matrix_output_file is None and - not eigendecomposition_requested): + not eigendecomposition_requested and + not error_analysis_requested): raise ValueError( "No analyses requested. Set at least one of:\n" " - resource_estimator (e.g., 'pyliqtr', 'cirq')\n" " - matrix_output_file (e.g., 'matrix.npz', 'matrix.h5', 'matrix.txt')\n" " - numerical_simulation_inputs (e.g., 'state.npy' or ['state1.npy', 'state2.npy'])\n" " - exact_matrix_output_file (e.g., 'exact_hamiltonian.npz')\n" - " - num_eigenvalues (e.g., 5 or 'all')" + " - num_eigenvalues (e.g., 5 or 'all')\n" + " - error_num_eigenvalues (e.g., 1)\n" + " - error_matrix_norms (e.g., 'frobenius' or ['frobenius', 'spectral'])\n" + " - error_state_inputs (e.g., 'state.npy')" ) results = {} @@ -585,14 +946,16 @@ def analyze_algorithm( config_analysis.matrix_output_file is not None or config_analysis.numerical_simulation_inputs is not None or (eigendecomposition_requested and - config_analysis.eigendecomposition_matrices in ['approximate', 'both']) + config_analysis.eigendecomposition_matrices in ['approximate', 'both']) or + (error_analysis_requested) # Error analysis always needs both ) # Check if any analysis needs the exact Hamiltonian matrix needs_exact_matrix = ( config_analysis.exact_matrix_output_file is not None or (eigendecomposition_requested and - config_analysis.eigendecomposition_matrices in ['exact', 'both']) + config_analysis.eigendecomposition_matrices in ['exact', 'both']) or + (error_analysis_requested) # Error analysis always needs both ) # Compute matrices once if needed @@ -633,11 +996,25 @@ def analyze_algorithm( exact_matrix=exact_matrix, unitary_matrix=unitary_matrix ) + results["eigendecomposition"] = eig_results + + if error_analysis_requested: + logger.info("Performing error analysis.") + # Pass eigendecomposition results if available + exact_eigendecomp = None + approx_eigendecomp = None + if eigendecomposition_requested and 'eig_results' in locals(): + # Eigendecompositions were computed, but we need to load them from files + # The error_analysis function will handle loading if needed + pass + results["error_analysis"] = error_analysis( + config_analysis, hamiltonian, algorithm, + exact_matrix=exact_matrix, + unitary_matrix=unitary_matrix, + exact_eigendecomp=exact_eigendecomp, + approx_eigendecomp=approx_eigendecomp + ) - # TODO: Add error estimation - # TODO: Add an option for detailed error analysis (explicitly compute the eigenvalues of the - # original Hamiltonian and the final unitary, compute ground state energy from both, - # compare the results; will only work for small systems) # TODO: Add gate parallelism / gate depth analysis # TODO: Would it be useful to analyze in terms of a different basis (e.g., Toffoli gates)? diff --git a/analysis/config_types.py b/analysis/config_types.py index 2750d53..38713f6 100644 --- a/analysis/config_types.py +++ b/analysis/config_types.py @@ -170,6 +170,9 @@ def __init__(self): self.num_eigenvalues = 0 self.eigendecomposition_matrices = 'approximate' self.which_eigenvalues = 'smallest' + self.error_num_eigenvalues = 0 + self.error_matrix_norms = None + self.error_state_inputs = None def _generate_TOML_table(self): table = tomlkit.table() @@ -181,6 +184,9 @@ def _generate_TOML_table(self): self.save_if_present(table, "num_eigenvalues") self.save_if_present(table, "eigendecomposition_matrices") self.save_if_present(table, "which_eigenvalues") + self.save_if_present(table, "error_num_eigenvalues") + self.save_if_present(table, "error_matrix_norms") + self.save_if_present(table, "error_state_inputs") return table # ------------------------------------------------------------------------------------------------- diff --git a/analysis/examples/config_full_analysis.py b/analysis/examples/config_full_analysis.py index b49be01..f527b9c 100644 --- a/analysis/examples/config_full_analysis.py +++ b/analysis/examples/config_full_analysis.py @@ -298,13 +298,72 @@ # analysis.eigendecomposition_matrices = "exact" # ================================================================================================= -# FUTURE ANALYSES (Coming in subsequent branches) +# 6. ERROR ANALYSIS (Branch 3: Now Available!) # ================================================================================================= +# Compare exact and approximate representations using three types of error metrics + +# analysis.error_num_eigenvalues = 1 # Compare ground state energies +# analysis.error_matrix_norms = "frobenius" # Options: "frobenius", "spectral", or ["frobenius", "spectral"] +# analysis.error_state_inputs = "examples/initial_state.npy" # Can be string or list + +# This analysis computes three independent error types: +# 1. Eigenvalue errors: Compare k smallest eigenvalues (typically k=1 for ground state) +# 2. Matrix norm errors: Frobenius and/or spectral norms of difference matrix +# 3. State-dependent errors: ||H_exact|ψ⟩ - H_approx|ψ⟩|| for specific states + +# Configuration options: +# error_num_eigenvalues: +# - 0 (default): Eigenvalue errors disabled +# - Positive int (e.g., 1): Compare k eigenvalues +# - Requires eigendecompositions (will compute if not already done) +# +# error_matrix_norms: +# - None (default): Matrix norm errors disabled +# - "frobenius": Fast element-wise difference norm +# - "spectral": Physically meaningful worst-case norm (slower) +# - ["frobenius", "spectral"]: Both norms +# - For large systems (>15 qubits): Uses matrix-free computation (can be slow) +# +# error_state_inputs: +# - None (default): State-dependent errors disabled +# - String or list of .npy files containing quantum states +# - Best-scaling error metric (can reach 30 qubits) +# - Fast: just applies operators to states + +# System size guidance: +# - ≤15 qubits: All error types fast +# - 16-20 qubits: Matrix norms slow, eigenvalue+state errors fast +# - 20-30 qubits: Use eigenvalue + state errors, avoid matrix norms +# - Production: error_num_eigenvalues=1 + error_state_inputs for best performance + +# Output: error_analysis.npz containing all computed errors -# Branch 3: Error Metrics -# analysis.error_num_eigenvalues = 1 # Compare ground state energy -# analysis.error_matrix_norms = ["frobenius", "spectral"] # Matrix difference norms -# analysis.error_state_inputs = ["examples/ground_state.npy"] # State-dependent errors +# Examples: + +# Ground state energy comparison (most common, always recommended) +# analysis.error_num_eigenvalues = 1 + +# Matrix norms (physical bounds, but expensive for large systems) +# analysis.error_matrix_norms = "frobenius" # Fast +# analysis.error_matrix_norms = ["frobenius", "spectral"] # Both (slower) + +# State-dependent errors (best scaling, physically relevant) +# analysis.error_state_inputs = "examples/initial_state.npy" +# analysis.error_state_inputs = ["examples/ground.npy", "examples/excited.npy"] + +# All three error types (comprehensive validation for small-medium systems) +# analysis.error_num_eigenvalues = 1 +# analysis.error_matrix_norms = "frobenius" +# analysis.error_state_inputs = ["examples/ground.npy", "examples/excited.npy"] + +# Minimal setup for large systems (20+ qubits) +# analysis.error_num_eigenvalues = 1 # Ground state only +# analysis.error_state_inputs = "examples/ground.npy" # Scales well +# # Skip error_matrix_norms for large systems (too slow) + +# ================================================================================================= +# FUTURE ANALYSES (Coming in subsequent branches) +# ================================================================================================= # Branch 4: Exact Numerical Simulation # analysis.exact_simulation_inputs = "examples/initial_state.npy" diff --git a/analysis/tests/test_error_analysis.py b/analysis/tests/test_error_analysis.py new file mode 100644 index 0000000..fc3d97b --- /dev/null +++ b/analysis/tests/test_error_analysis.py @@ -0,0 +1,686 @@ +""" +Tests for error analysis functionality. + +Tests cover: +- Eigenvalue error computation +- Matrix norm errors (Frobenius and spectral) +- State-dependent errors +- Integration with eigendecomposition +- File I/O for error results +""" + +import numpy as np +import pytest +import tempfile +import os +from pathlib import Path + +from qhat.analysis.analysis import error_analysis +from qhat.analysis.config_types import AnalysisConfiguration +from qhat.analysis.file_io import save_state + + +# ================================================================================================= +# Unit Tests: Eigenvalue Errors +# ================================================================================================= + +def test_eigenvalue_error_zero_when_identical(): + """Test that eigenvalue error is zero when matrices are identical.""" + from qhat.analysis.file_io import save_eigendecomposition + + # Create identical eigendecompositions + eigenvalues = np.array([1.0, 2.0, 3.0]) + eigenvectors = np.eye(3, dtype=complex) + + config = AnalysisConfiguration() + config.error_num_eigenvalues = 3 + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save both decompositions + save_eigendecomposition( + 'exact_eigendecomposition.npz', eigenvalues, eigenvectors, + matrix_type='exact', num_eigenvalues=3, which_eigenvalues='smallest' + ) + save_eigendecomposition( + 'approximate_eigendecomposition.npz', eigenvalues, eigenvectors, + matrix_type='approximate', num_eigenvalues=3, which_eigenvalues='smallest' + ) + + # Compute errors + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=None, + unitary_matrix=None + ) + + # Verify zero errors + assert 'eigenvalue_errors' in results + np.testing.assert_array_almost_equal( + results['eigenvalue_errors']['absolute_errors'], + [0.0, 0.0, 0.0] + ) + + finally: + os.chdir(original_dir) + + +def test_eigenvalue_error_nonzero_when_different(): + """Test that eigenvalue error is nonzero when matrices differ.""" + from qhat.analysis.file_io import save_eigendecomposition + + # Create different eigendecompositions + exact_eigenvalues = np.array([1.0, 2.0, 3.0]) + approx_eigenvalues = np.array([1.1, 1.9, 3.2]) + eigenvectors = np.eye(3, dtype=complex) + + config = AnalysisConfiguration() + config.error_num_eigenvalues = 3 + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save both decompositions + save_eigendecomposition( + 'exact_eigendecomposition.npz', exact_eigenvalues, eigenvectors, + matrix_type='exact', num_eigenvalues=3, which_eigenvalues='smallest' + ) + save_eigendecomposition( + 'approximate_eigendecomposition.npz', approx_eigenvalues, eigenvectors, + matrix_type='approximate', num_eigenvalues=3, which_eigenvalues='smallest' + ) + + # Compute errors + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=None, + unitary_matrix=None + ) + + # Verify nonzero errors + assert 'eigenvalue_errors' in results + expected_errors = exact_eigenvalues - approx_eigenvalues + np.testing.assert_array_almost_equal( + results['eigenvalue_errors']['absolute_errors'], + expected_errors.tolist() + ) + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Unit Tests: Matrix Norm Errors (Dense) +# ================================================================================================= + +def test_frobenius_norm_zero_when_identical(): + """Test Frobenius norm error is zero for identical matrices.""" + # 2x2 identity + matrix = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.error_matrix_norms = 'frobenius' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=matrix, + unitary_matrix=matrix + ) + + assert 'matrix_frobenius_error' in results + assert results['matrix_frobenius_error'] < 1e-10 + + finally: + os.chdir(original_dir) + + +def test_frobenius_norm_nonzero_when_different(): + """Test Frobenius norm error is nonzero for different matrices.""" + exact = np.eye(2, dtype=complex) + approx = np.array([[1.0, 0.1], [0.0, 1.0]], dtype=complex) + + config = AnalysisConfiguration() + config.error_matrix_norms = 'frobenius' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=exact, + unitary_matrix=approx + ) + + assert 'matrix_frobenius_error' in results + # Manual calculation: ||[[0, 0.1], [0, 0]]||_F = sqrt(0.01) = 0.1 + expected = np.linalg.norm(exact - approx, 'fro') + np.testing.assert_almost_equal( + results['matrix_frobenius_error'], + expected + ) + + finally: + os.chdir(original_dir) + + +def test_spectral_norm_zero_when_identical(): + """Test spectral norm error is zero for identical matrices.""" + matrix = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.error_matrix_norms = 'spectral' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=matrix, + unitary_matrix=matrix + ) + + assert 'matrix_spectral_error' in results + assert results['matrix_spectral_error'] < 1e-10 + + finally: + os.chdir(original_dir) + + +def test_spectral_norm_nonzero_when_different(): + """Test spectral norm error is nonzero for different matrices.""" + exact = np.eye(2, dtype=complex) + approx = np.array([[1.0, 0.3], [0.0, 1.0]], dtype=complex) + + config = AnalysisConfiguration() + config.error_matrix_norms = 'spectral' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=exact, + unitary_matrix=approx + ) + + assert 'matrix_spectral_error' in results + # Manual calculation: ||[[0, 0.3], [0, 0]]||_2 = 0.3 + expected = np.linalg.norm(exact - approx, 2) + np.testing.assert_almost_equal( + results['matrix_spectral_error'], + expected, + decimal=5 + ) + + finally: + os.chdir(original_dir) + + +def test_both_matrix_norms(): + """Test computing both Frobenius and spectral norms.""" + exact = np.eye(2, dtype=complex) + approx = np.array([[1.0, 0.2], [0.0, 1.0]], dtype=complex) + + config = AnalysisConfiguration() + config.error_matrix_norms = ['frobenius', 'spectral'] + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=exact, + unitary_matrix=approx + ) + + # Both should be present + assert 'matrix_frobenius_error' in results + assert 'matrix_spectral_error' in results + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Unit Tests: State-Dependent Errors +# ================================================================================================= + +def test_state_error_zero_when_identical(): + """Test state error is zero for identical operators.""" + matrix = np.eye(2, dtype=complex) + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.error_state_inputs = 'test_state.npy' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save state + save_state('test_state.npy', state) + + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=matrix, + unitary_matrix=matrix + ) + + assert 'state_errors' in results + assert len(results['state_errors']) == 1 + assert results['state_errors'][0]['absolute_error'] < 1e-10 + + finally: + os.chdir(original_dir) + + +def test_state_error_nonzero_when_different(): + """Test state error is nonzero for different operators.""" + exact = np.eye(2, dtype=complex) + approx = np.array([[1.0, 0.1], [0.0, 1.0]], dtype=complex) + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.error_state_inputs = 'test_state.npy' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save state + save_state('test_state.npy', state) + + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=exact, + unitary_matrix=approx + ) + + assert 'state_errors' in results + assert len(results['state_errors']) == 1 + + # Manual calculation + exact_final = exact @ state # [1.0, 0.0] + approx_final = approx @ state # [1.0, 0.0] + expected_error = np.linalg.norm(exact_final - approx_final) + + np.testing.assert_almost_equal( + results['state_errors'][0]['absolute_error'], + expected_error + ) + + finally: + os.chdir(original_dir) + + +def test_state_error_multiple_states(): + """Test state errors for multiple input states.""" + matrix = np.eye(2, dtype=complex) + state1 = np.array([1.0, 0.0], dtype=complex) + state2 = np.array([0.0, 1.0], dtype=complex) + + config = AnalysisConfiguration() + config.error_state_inputs = ['state1.npy', 'state2.npy'] + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save states + save_state('state1.npy', state1) + save_state('state2.npy', state2) + + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=matrix, + unitary_matrix=matrix + ) + + assert 'state_errors' in results + assert len(results['state_errors']) == 2 + assert results['state_errors'][0]['input_file'] == 'state1.npy' + assert results['state_errors'][1]['input_file'] == 'state2.npy' + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Integration Tests: All Error Types Together +# ================================================================================================= + +def test_all_error_types_together(): + """Test computing all three error types in one analysis.""" + from qhat.analysis.file_io import save_eigendecomposition + + # Setup matrices + exact = np.diag([1.0, 2.0]) + approx = np.diag([1.1, 1.9]) + state = np.array([1.0, 0.0], dtype=complex) + eigenvalues_exact = np.array([1.0, 2.0]) + eigenvalues_approx = np.array([1.1, 1.9]) + eigenvectors = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.error_num_eigenvalues = 2 + config.error_matrix_norms = ['frobenius', 'spectral'] + config.error_state_inputs = 'test_state.npy' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save eigendecompositions + save_eigendecomposition( + 'exact_eigendecomposition.npz', eigenvalues_exact, eigenvectors, + matrix_type='exact', num_eigenvalues=2, which_eigenvalues='smallest' + ) + save_eigendecomposition( + 'approximate_eigendecomposition.npz', eigenvalues_approx, eigenvectors, + matrix_type='approximate', num_eigenvalues=2, which_eigenvalues='smallest' + ) + + # Save state + save_state('test_state.npy', state) + + # Compute all errors + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=exact, + unitary_matrix=approx + ) + + # Verify all three types present + assert 'eigenvalue_errors' in results + assert 'matrix_frobenius_error' in results + assert 'matrix_spectral_error' in results + assert 'state_errors' in results + + finally: + os.chdir(original_dir) + + +def test_error_output_file_created(): + """Test that error_analysis.npz file is created.""" + from qhat.analysis.file_io import save_eigendecomposition + + eigenvalues = np.array([1.0, 2.0]) + eigenvectors = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.error_num_eigenvalues = 2 + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save eigendecompositions + save_eigendecomposition( + 'exact_eigendecomposition.npz', eigenvalues, eigenvectors, + matrix_type='exact', num_eigenvalues=2, which_eigenvalues='smallest' + ) + save_eigendecomposition( + 'approximate_eigendecomposition.npz', eigenvalues, eigenvectors, + matrix_type='approximate', num_eigenvalues=2, which_eigenvalues='smallest' + ) + + # Compute errors + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=None, + unitary_matrix=None + ) + + # Verify file created + assert 'output_file' in results + assert os.path.exists('error_analysis.npz') + + # Verify can load + data = np.load('error_analysis.npz') + assert 'eigenvalue_absolute_errors' in data + assert 'eigenvalue_relative_errors' in data + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Error Handling Tests +# ================================================================================================= + +def test_invalid_matrix_norm_type(): + """Test error when requesting invalid matrix norm type.""" + matrix = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.error_matrix_norms = 'invalid' + + with pytest.raises(ValueError, match="Unknown matrix norm type"): + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=matrix, + unitary_matrix=matrix + ) + finally: + os.chdir(original_dir) + + +def test_missing_state_file(): + """Test error when state file doesn't exist.""" + matrix = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.error_state_inputs = 'nonexistent.npy' + + with pytest.raises(Exception): # Will be FileNotFoundError or similar + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=matrix, + unitary_matrix=matrix + ) + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Matrix-Free Operator Tests +# ================================================================================================= + +def test_state_error_with_matrix_free_operators(): + """Test state errors work with matrix-free operators.""" + from qhat.analysis.matrix_operations import PauliStringOperator + + # Create matrix-free operators + pauli_dict = {'II': 1.0} # Identity operator + exact_op = PauliStringOperator(pauli_dict, num_qubits=2) + approx_op = PauliStringOperator(pauli_dict, num_qubits=2) + + state = np.array([1.0, 0.0, 0.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.error_state_inputs = 'test_state.npy' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save state + save_state('test_state.npy', state) + + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=exact_op, + unitary_matrix=approx_op + ) + + assert 'state_errors' in results + # Error should be zero for identical operators + assert results['state_errors'][0]['absolute_error'] < 1e-10 + + finally: + os.chdir(original_dir) + + +def test_frobenius_norm_with_matrix_free_small(): + """Test Frobenius norm computation with matrix-free operators (small system).""" + from qhat.analysis.matrix_operations import PauliStringOperator + + # 2-qubit system (dimension 4) - still small enough for both paths + pauli_dict_exact = {'II': 1.0} + pauli_dict_approx = {'II': 1.0} + + exact_op = PauliStringOperator(pauli_dict_exact, num_qubits=2) + approx_op = PauliStringOperator(pauli_dict_approx, num_qubits=2) + + config = AnalysisConfiguration() + config.error_matrix_norms = 'frobenius' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=exact_op, + unitary_matrix=approx_op + ) + + assert 'matrix_frobenius_error' in results + # Identical operators should have zero error + assert results['matrix_frobenius_error'] < 1e-10 + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Relative Error Tests +# ================================================================================================= + +def test_eigenvalue_relative_error(): + """Test that relative eigenvalue errors are computed correctly.""" + from qhat.analysis.file_io import save_eigendecomposition + + exact_eigenvalues = np.array([10.0, 20.0]) + approx_eigenvalues = np.array([11.0, 18.0]) + eigenvectors = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.error_num_eigenvalues = 2 + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_eigendecomposition( + 'exact_eigendecomposition.npz', exact_eigenvalues, eigenvectors, + matrix_type='exact', num_eigenvalues=2, which_eigenvalues='smallest' + ) + save_eigendecomposition( + 'approximate_eigendecomposition.npz', approx_eigenvalues, eigenvectors, + matrix_type='approximate', num_eigenvalues=2, which_eigenvalues='smallest' + ) + + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=None, + unitary_matrix=None + ) + + # Check relative errors + # (10 - 11) / |10| = -0.1, (20 - 18) / |20| = 0.1 + expected_relative = np.array([-0.1, 0.1]) + np.testing.assert_array_almost_equal( + results['eigenvalue_errors']['relative_errors'], + expected_relative.tolist() + ) + + finally: + os.chdir(original_dir) + + +def test_state_relative_error(): + """Test that relative state errors are computed correctly.""" + exact = np.diag([2.0, 3.0]) + approx = np.diag([2.1, 2.9]) + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.error_state_inputs = 'test_state.npy' + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('test_state.npy', state) + + results = error_analysis( + config, + hamiltonian=None, + algorithm=None, + exact_matrix=exact, + unitary_matrix=approx + ) + + # exact @ state = [2.0, 0.0], norm = 2.0 + # approx @ state = [2.1, 0.0], norm = 2.1 + # absolute error = 0.1 + # relative error = 0.1 / 2.0 = 0.05 + assert 'state_errors' in results + assert results['state_errors'][0]['relative_error'] == pytest.approx(0.05) + + finally: + os.chdir(original_dir) From 3c3054c60e63f2606308b19fde9e42152a4a94ee Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 15 Jun 2026 08:45:42 -0600 Subject: [PATCH 03/12] Modify config_full_analysis.py --- analysis/examples/config_full_analysis.py | 33 +++++++---------------- 1 file changed, 9 insertions(+), 24 deletions(-) diff --git a/analysis/examples/config_full_analysis.py b/analysis/examples/config_full_analysis.py index f527b9c..6968e0c 100644 --- a/analysis/examples/config_full_analysis.py +++ b/analysis/examples/config_full_analysis.py @@ -235,13 +235,13 @@ # analysis.eigendecomposition_matrices = "exact" # ================================================================================================= -# 5. EIGENDECOMPOSITION ANALYSIS (Branch 2: Now Available!) +# 5. EIGENDECOMPOSITION ANALYSIS # ================================================================================================= # Compute eigenvalues and eigenvectors of exact and/or approximate matrices -# analysis.num_eigenvalues = 5 # Compute 5 eigenvalues (use sparse methods) -# analysis.eigendecomposition_matrices = "both" # Options: "exact", "approximate", "both" -# analysis.which_eigenvalues = "smallest" # Options: "smallest", "largest", "both" +analysis.num_eigenvalues = 5 # Compute 5 eigenvalues (use sparse methods) +analysis.eigendecomposition_matrices = "both" # Options: "exact", "approximate", "both" +analysis.which_eigenvalues = "smallest" # Options: "smallest", "largest", "both" # This analysis computes eigendecompositions for validation and comparison: # - Ground state energy (use num_eigenvalues=1, which_eigenvalues="smallest") @@ -264,13 +264,6 @@ # - "smallest" (default): Algebraically smallest (most negative, ground state) # - "largest": Algebraically largest (most positive) # - "both": Compute smallest AND largest (returns 2k eigenvalues) -# -# System size guidance: -# - ≤10 qubits: "all" feasible -# - 10-12 qubits: "all" possible but expensive -# - ≥12 qubits: Use partial (specify k=5-10) -# - 20+ qubits: MUST use partial -# - 25-30 qubits: Partial still scales well # Output files: # - exact_eigendecomposition.npz (if "exact" or "both") @@ -298,13 +291,13 @@ # analysis.eigendecomposition_matrices = "exact" # ================================================================================================= -# 6. ERROR ANALYSIS (Branch 3: Now Available!) +# 6. ERROR ANALYSIS # ================================================================================================= # Compare exact and approximate representations using three types of error metrics -# analysis.error_num_eigenvalues = 1 # Compare ground state energies -# analysis.error_matrix_norms = "frobenius" # Options: "frobenius", "spectral", or ["frobenius", "spectral"] -# analysis.error_state_inputs = "examples/initial_state.npy" # Can be string or list +analysis.error_num_eigenvalues = 1 # Compare ground state energies +analysis.error_matrix_norms = "frobenius" # Options: "frobenius", "spectral", or ["frobenius", "spectral"] +analysis.error_state_inputs = "examples/initial_state.npy" # Can be string or list # This analysis computes three independent error types: # 1. Eigenvalue errors: Compare k smallest eigenvalues (typically k=1 for ground state) @@ -322,20 +315,12 @@ # - "frobenius": Fast element-wise difference norm # - "spectral": Physically meaningful worst-case norm (slower) # - ["frobenius", "spectral"]: Both norms -# - For large systems (>15 qubits): Uses matrix-free computation (can be slow) # # error_state_inputs: # - None (default): State-dependent errors disabled # - String or list of .npy files containing quantum states -# - Best-scaling error metric (can reach 30 qubits) # - Fast: just applies operators to states -# System size guidance: -# - ≤15 qubits: All error types fast -# - 16-20 qubits: Matrix norms slow, eigenvalue+state errors fast -# - 20-30 qubits: Use eigenvalue + state errors, avoid matrix norms -# - Production: error_num_eigenvalues=1 + error_state_inputs for best performance - # Output: error_analysis.npz containing all computed errors # Examples: @@ -356,7 +341,7 @@ # analysis.error_matrix_norms = "frobenius" # analysis.error_state_inputs = ["examples/ground.npy", "examples/excited.npy"] -# Minimal setup for large systems (20+ qubits) +# Minimal setup for large systems # analysis.error_num_eigenvalues = 1 # Ground state only # analysis.error_state_inputs = "examples/ground.npy" # Scales well # # Skip error_matrix_norms for large systems (too slow) From b99d282248add273f9b3a457ede3a8774eece4d2 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 15 Jun 2026 09:10:31 -0600 Subject: [PATCH 04/12] Some changes (mix of manual and Claude). --- analysis/README.md | 19 +++---- analysis/analysis.py | 64 +++++++++-------------- analysis/config_types.py | 4 +- analysis/examples/config_full_analysis.py | 21 ++++---- analysis/tests/test_error_analysis.py | 60 +++++++++++++++------ 5 files changed, 92 insertions(+), 76 deletions(-) diff --git a/analysis/README.md b/analysis/README.md index d314c4d..e29e76a 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -290,11 +290,12 @@ There are many details of the algorithm that may be worth analyzing. The availab **Configuration parameters**: - - **`error_num_eigenvalues`**: Number of eigenvalues to compare (default: 0, disabled) - - Set to positive integer (e.g., `1` for ground state energy comparison) - - Compares the k smallest eigenvalues from exact vs approximate - - Requires eigendecompositions (will compute if not already available) - - **Best for**: Validating ground state energy accuracy + - **`enable_eigenvalue_errors`**: Enable eigenvalue error comparison (default: False, disabled) + - Set to `True` to compute errors for ALL eigenvalues from the eigendecomposition + - The number of eigenvalues compared is determined by `num_eigenvalues` setting + - Compares all eigenvalues computed in both exact and approximate eigendecompositions + - Requires `eigendecomposition_matrices = "both"` to ensure both decompositions are computed + - **Best for**: Validating eigenvalue accuracy across the spectrum - **`error_matrix_norms`**: Which matrix norms to compute (default: None, disabled) - Single string: `"frobenius"` or `"spectral"` @@ -331,8 +332,8 @@ There are many details of the algorithm that may be worth analyzing. The availab **Examples**: ```python - # Ground state energy comparison (most common) - analysis.error_num_eigenvalues = 1 + # Eigenvalue error comparison (compares all eigenvalues from eigendecomposition) + analysis.enable_eigenvalue_errors = True # Matrix norm errors (physical bounds) analysis.error_matrix_norms = "frobenius" # Fast @@ -345,13 +346,13 @@ There are many details of the algorithm that may be worth analyzing. The availab analysis.error_state_inputs = ["ground.npy", "excited.npy"] # All three error types together (comprehensive validation) - analysis.error_num_eigenvalues = 1 + analysis.enable_eigenvalue_errors = True analysis.error_matrix_norms = "frobenius" analysis.error_state_inputs = ["ground.npy", "excited.npy"] ``` **When to use each error type**: - - **Eigenvalue errors**: Always use for ground state energy validation (k=1) + - **Eigenvalue errors**: Use when you want to validate all eigenvalues computed in the eigendecomposition - **Matrix norm errors**: Use for physical bounds, but expensive for large systems - **State errors**: Use for specific physically relevant states, scales best diff --git a/analysis/analysis.py b/analysis/analysis.py index 59cdaa2..06ab05f 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -479,7 +479,7 @@ def error_analysis( Compute error metrics comparing exact and approximate representations. Three independent error types: - 1. Eigenvalue errors (if error_num_eigenvalues > 0) + 1. Eigenvalue errors (if enable_eigenvalue_errors is True) 2. Matrix norm errors (if error_matrix_norms is not None) 3. State-dependent errors (if error_state_inputs is not None) @@ -507,54 +507,38 @@ def error_analysis( # 1. EIGENVALUE ERROR # ============================================================================================= - if config_analysis.error_num_eigenvalues > 0: - logger.info(f"Computing eigenvalue errors for {config_analysis.error_num_eigenvalues} eigenvalues") + if config_analysis.enable_eigenvalue_errors: + logger.info("Computing eigenvalue errors for all eigenvalues from eigendecomposition") - # Load or compute eigendecompositions - if exact_eigendecomp is None: - if os.path.exists('exact_eigendecomposition.npz'): - exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') - else: - logger.info("Computing exact eigendecomposition for error analysis") - # Compute it - config_temp = AnalysisConfiguration() - config_temp.num_eigenvalues = config_analysis.error_num_eigenvalues - config_temp.eigendecomposition_matrices = 'exact' - config_temp.which_eigenvalues = 'smallest' - eig_results = eigendecomposition_analysis( - config_temp, hamiltonian, algorithm, - exact_matrix=exact_matrix, unitary_matrix=None - ) - exact_eigendecomp = load_eigendecomposition(eig_results['exact_eigendecomposition']['file']) + if exact_eigendecomp is None or approx_eigendecomp is None: + raise ValueError( + "Both eigendecompositions must be computed in order to compare eigenvalues. " + "Ensure eigendecomposition_matrices is set to 'both' when enable_eigenvalue_errors is True." + ) - if approx_eigendecomp is None: - if os.path.exists('approximate_eigendecomposition.npz'): - approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') - else: - logger.info("Computing approximate eigendecomposition for error analysis") - # Compute it - config_temp = AnalysisConfiguration() - config_temp.num_eigenvalues = config_analysis.error_num_eigenvalues - config_temp.eigendecomposition_matrices = 'approximate' - config_temp.which_eigenvalues = 'smallest' - eig_results = eigendecomposition_analysis( - config_temp, hamiltonian, algorithm, - exact_matrix=None, unitary_matrix=unitary_matrix - ) - approx_eigendecomp = load_eigendecomposition(eig_results['approximate_eigendecomposition']['file']) + # Get all eigenvalues from the eigendecompositions + exact_eigs = exact_eigendecomp['eigenvalues'] + approx_eigs = approx_eigendecomp['eigenvalues'] - # Compare eigenvalues - exact_eigs = exact_eigendecomp['eigenvalues'][:config_analysis.error_num_eigenvalues] - approx_eigs = approx_eigendecomp['eigenvalues'][:config_analysis.error_num_eigenvalues] + # Verify that the same number of eigenvalues were computed + if len(exact_eigs) != len(approx_eigs): + raise ValueError( + f"Mismatch in number of eigenvalues: exact has {len(exact_eigs)}, " + f"approximate has {len(approx_eigs)}. Both eigendecompositions must compute " + "the same number of eigenvalues." + ) + # Compare all eigenvalues absolute_errors = exact_eigs - approx_eigs relative_errors = absolute_errors / np.abs(exact_eigs) + num_eigenvalues = len(exact_eigs) + logger.info(f"Computed errors for {num_eigenvalues} eigenvalues") logger.info(f"Eigenvalue absolute error range: [{absolute_errors.min():.6e}, {absolute_errors.max():.6e}]") logger.info(f"Eigenvalue relative error range: [{relative_errors.min():.6e}, {relative_errors.max():.6e}]") results['eigenvalue_errors'] = { - 'num_eigenvalues': config_analysis.error_num_eigenvalues, + 'num_eigenvalues': num_eigenvalues, 'absolute_errors': absolute_errors.tolist(), 'relative_errors': relative_errors.tolist(), 'max_absolute_error': float(np.abs(absolute_errors).max()), @@ -916,7 +900,7 @@ def analyze_algorithm( ) error_analysis_requested = ( - config_analysis.error_num_eigenvalues > 0 or + config_analysis.enable_eigenvalue_errors or config_analysis.error_matrix_norms is not None or config_analysis.error_state_inputs is not None ) @@ -934,7 +918,7 @@ def analyze_algorithm( " - numerical_simulation_inputs (e.g., 'state.npy' or ['state1.npy', 'state2.npy'])\n" " - exact_matrix_output_file (e.g., 'exact_hamiltonian.npz')\n" " - num_eigenvalues (e.g., 5 or 'all')\n" - " - error_num_eigenvalues (e.g., 1)\n" + " - enable_eigenvalue_errors (True to compute errors for all eigenvalues)\n" " - error_matrix_norms (e.g., 'frobenius' or ['frobenius', 'spectral'])\n" " - error_state_inputs (e.g., 'state.npy')" ) diff --git a/analysis/config_types.py b/analysis/config_types.py index 38713f6..82789fe 100644 --- a/analysis/config_types.py +++ b/analysis/config_types.py @@ -170,7 +170,7 @@ def __init__(self): self.num_eigenvalues = 0 self.eigendecomposition_matrices = 'approximate' self.which_eigenvalues = 'smallest' - self.error_num_eigenvalues = 0 + self.enable_eigenvalue_errors = False self.error_matrix_norms = None self.error_state_inputs = None @@ -184,7 +184,7 @@ def _generate_TOML_table(self): self.save_if_present(table, "num_eigenvalues") self.save_if_present(table, "eigendecomposition_matrices") self.save_if_present(table, "which_eigenvalues") - self.save_if_present(table, "error_num_eigenvalues") + self.save_if_present(table, "enable_eigenvalue_errors") self.save_if_present(table, "error_matrix_norms") self.save_if_present(table, "error_state_inputs") return table diff --git a/analysis/examples/config_full_analysis.py b/analysis/examples/config_full_analysis.py index 6968e0c..e7af219 100644 --- a/analysis/examples/config_full_analysis.py +++ b/analysis/examples/config_full_analysis.py @@ -295,20 +295,21 @@ # ================================================================================================= # Compare exact and approximate representations using three types of error metrics -analysis.error_num_eigenvalues = 1 # Compare ground state energies +analysis.enable_eigenvalue_errors = True # Compare all eigenvalues from eigendecomposition analysis.error_matrix_norms = "frobenius" # Options: "frobenius", "spectral", or ["frobenius", "spectral"] analysis.error_state_inputs = "examples/initial_state.npy" # Can be string or list # This analysis computes three independent error types: -# 1. Eigenvalue errors: Compare k smallest eigenvalues (typically k=1 for ground state) +# 1. Eigenvalue errors: Compare all eigenvalues computed in the eigendecomposition # 2. Matrix norm errors: Frobenius and/or spectral norms of difference matrix # 3. State-dependent errors: ||H_exact|ψ⟩ - H_approx|ψ⟩|| for specific states # Configuration options: -# error_num_eigenvalues: -# - 0 (default): Eigenvalue errors disabled -# - Positive int (e.g., 1): Compare k eigenvalues -# - Requires eigendecompositions (will compute if not already done) +# enable_eigenvalue_errors: +# - False (default): Eigenvalue errors disabled +# - True: Compute errors for ALL eigenvalues from eigendecomposition +# - Requires eigendecompositions (set eigendecomposition_matrices = "both") +# - The number of eigenvalues compared is determined by num_eigenvalues # # error_matrix_norms: # - None (default): Matrix norm errors disabled @@ -325,8 +326,8 @@ # Examples: -# Ground state energy comparison (most common, always recommended) -# analysis.error_num_eigenvalues = 1 +# Eigenvalue error comparison (compare all eigenvalues from eigendecomposition) +# analysis.enable_eigenvalue_errors = True # Matrix norms (physical bounds, but expensive for large systems) # analysis.error_matrix_norms = "frobenius" # Fast @@ -337,12 +338,12 @@ # analysis.error_state_inputs = ["examples/ground.npy", "examples/excited.npy"] # All three error types (comprehensive validation for small-medium systems) -# analysis.error_num_eigenvalues = 1 +# analysis.enable_eigenvalue_errors = True # analysis.error_matrix_norms = "frobenius" # analysis.error_state_inputs = ["examples/ground.npy", "examples/excited.npy"] # Minimal setup for large systems -# analysis.error_num_eigenvalues = 1 # Ground state only +# analysis.enable_eigenvalue_errors = True # Compare eigenvalues # analysis.error_state_inputs = "examples/ground.npy" # Scales well # # Skip error_matrix_norms for large systems (too slow) diff --git a/analysis/tests/test_error_analysis.py b/analysis/tests/test_error_analysis.py index fc3d97b..0a50e11 100644 --- a/analysis/tests/test_error_analysis.py +++ b/analysis/tests/test_error_analysis.py @@ -26,14 +26,14 @@ def test_eigenvalue_error_zero_when_identical(): """Test that eigenvalue error is zero when matrices are identical.""" - from qhat.analysis.file_io import save_eigendecomposition + from qhat.analysis.file_io import save_eigendecomposition, load_eigendecomposition # Create identical eigendecompositions eigenvalues = np.array([1.0, 2.0, 3.0]) eigenvectors = np.eye(3, dtype=complex) config = AnalysisConfiguration() - config.error_num_eigenvalues = 3 + config.enable_eigenvalue_errors = True with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -49,13 +49,19 @@ def test_eigenvalue_error_zero_when_identical(): matrix_type='approximate', num_eigenvalues=3, which_eigenvalues='smallest' ) + # Load eigendecompositions + exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') + approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') + # Compute errors results = error_analysis( config, hamiltonian=None, algorithm=None, exact_matrix=None, - unitary_matrix=None + unitary_matrix=None, + exact_eigendecomp=exact_eigendecomp, + approx_eigendecomp=approx_eigendecomp ) # Verify zero errors @@ -71,7 +77,7 @@ def test_eigenvalue_error_zero_when_identical(): def test_eigenvalue_error_nonzero_when_different(): """Test that eigenvalue error is nonzero when matrices differ.""" - from qhat.analysis.file_io import save_eigendecomposition + from qhat.analysis.file_io import save_eigendecomposition, load_eigendecomposition # Create different eigendecompositions exact_eigenvalues = np.array([1.0, 2.0, 3.0]) @@ -79,7 +85,7 @@ def test_eigenvalue_error_nonzero_when_different(): eigenvectors = np.eye(3, dtype=complex) config = AnalysisConfiguration() - config.error_num_eigenvalues = 3 + config.enable_eigenvalue_errors = True with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -95,13 +101,19 @@ def test_eigenvalue_error_nonzero_when_different(): matrix_type='approximate', num_eigenvalues=3, which_eigenvalues='smallest' ) + # Load eigendecompositions + exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') + approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') + # Compute errors results = error_analysis( config, hamiltonian=None, algorithm=None, exact_matrix=None, - unitary_matrix=None + unitary_matrix=None, + exact_eigendecomp=exact_eigendecomp, + approx_eigendecomp=approx_eigendecomp ) # Verify nonzero errors @@ -382,7 +394,7 @@ def test_state_error_multiple_states(): def test_all_error_types_together(): """Test computing all three error types in one analysis.""" - from qhat.analysis.file_io import save_eigendecomposition + from qhat.analysis.file_io import save_eigendecomposition, load_eigendecomposition # Setup matrices exact = np.diag([1.0, 2.0]) @@ -393,7 +405,7 @@ def test_all_error_types_together(): eigenvectors = np.eye(2, dtype=complex) config = AnalysisConfiguration() - config.error_num_eigenvalues = 2 + config.enable_eigenvalue_errors = True config.error_matrix_norms = ['frobenius', 'spectral'] config.error_state_inputs = 'test_state.npy' @@ -411,6 +423,10 @@ def test_all_error_types_together(): matrix_type='approximate', num_eigenvalues=2, which_eigenvalues='smallest' ) + # Load eigendecompositions + exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') + approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') + # Save state save_state('test_state.npy', state) @@ -420,7 +436,9 @@ def test_all_error_types_together(): hamiltonian=None, algorithm=None, exact_matrix=exact, - unitary_matrix=approx + unitary_matrix=approx, + exact_eigendecomp=exact_eigendecomp, + approx_eigendecomp=approx_eigendecomp ) # Verify all three types present @@ -435,13 +453,13 @@ def test_all_error_types_together(): def test_error_output_file_created(): """Test that error_analysis.npz file is created.""" - from qhat.analysis.file_io import save_eigendecomposition + from qhat.analysis.file_io import save_eigendecomposition, load_eigendecomposition eigenvalues = np.array([1.0, 2.0]) eigenvectors = np.eye(2, dtype=complex) config = AnalysisConfiguration() - config.error_num_eigenvalues = 2 + config.enable_eigenvalue_errors = True with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -457,13 +475,19 @@ def test_error_output_file_created(): matrix_type='approximate', num_eigenvalues=2, which_eigenvalues='smallest' ) + # Load eigendecompositions + exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') + approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') + # Compute errors results = error_analysis( config, hamiltonian=None, algorithm=None, exact_matrix=None, - unitary_matrix=None + unitary_matrix=None, + exact_eigendecomp=exact_eigendecomp, + approx_eigendecomp=approx_eigendecomp ) # Verify file created @@ -610,14 +634,14 @@ def test_frobenius_norm_with_matrix_free_small(): def test_eigenvalue_relative_error(): """Test that relative eigenvalue errors are computed correctly.""" - from qhat.analysis.file_io import save_eigendecomposition + from qhat.analysis.file_io import save_eigendecomposition, load_eigendecomposition exact_eigenvalues = np.array([10.0, 20.0]) approx_eigenvalues = np.array([11.0, 18.0]) eigenvectors = np.eye(2, dtype=complex) config = AnalysisConfiguration() - config.error_num_eigenvalues = 2 + config.enable_eigenvalue_errors = True with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -632,12 +656,18 @@ def test_eigenvalue_relative_error(): matrix_type='approximate', num_eigenvalues=2, which_eigenvalues='smallest' ) + # Load eigendecompositions + exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') + approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') + results = error_analysis( config, hamiltonian=None, algorithm=None, exact_matrix=None, - unitary_matrix=None + unitary_matrix=None, + exact_eigendecomp=exact_eigendecomp, + approx_eigendecomp=approx_eigendecomp ) # Check relative errors From 71932df12274491855414a56ef7d4e70ec96caca Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 15 Jun 2026 10:12:28 -0600 Subject: [PATCH 05/12] Further edits --- analysis/analysis.py | 225 +++++++++++++-- analysis/driver.py | 8 +- analysis/tests/test_config_validation.py | 114 ++++++++ analysis/tests/test_driver_validation.py | 141 +++++++++ .../test_expensive_analysis_detection.py | 268 ++++++++++++++++++ 5 files changed, 733 insertions(+), 23 deletions(-) create mode 100644 analysis/tests/test_config_validation.py create mode 100644 analysis/tests/test_driver_validation.py create mode 100644 analysis/tests/test_expensive_analysis_detection.py diff --git a/analysis/analysis.py b/analysis/analysis.py index 06ab05f..c678ef4 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -883,15 +883,24 @@ def numerical_simulation( return {'simulations': results} # ------------------------------------------------------------------------------------------------- +# Functions to determine what expensive computations are required +# ------------------------------------------------------------------------------------------------- -def analyze_algorithm( - config_analysis: AnalysisConfiguration, - algorithm, - hamiltonian=None) -> dict: +def requires_exact_eigendecomposition(config_analysis: AnalysisConfiguration) -> bool: + """ + Determine if exact eigendecomposition needs to be computed. - logger.info("Beginning algorithm analysis.") + Exact eigendecomposition is required for: + - Eigendecomposition analysis with eigendecomposition_matrices = 'exact' or 'both' + - Eigenvalue error analysis (always needs both eigendecompositions) - # Validate at least one analysis requested + Parameters: + config_analysis: Analysis configuration + + Returns: + True if exact eigendecomposition computation is needed, False otherwise + """ + # Check if eigendecomposition is requested at all num_eigenvalues = config_analysis.num_eigenvalues eigendecomposition_requested = ( isinstance(num_eigenvalues, int) and num_eigenvalues > 0 @@ -899,12 +908,179 @@ def analyze_algorithm( isinstance(num_eigenvalues, str) and num_eigenvalues.lower() == "all" ) + # Need exact eigendecomposition if: + # 1. Eigendecomposition requested and matrices setting includes 'exact' or 'both' + # 2. Eigenvalue error analysis is enabled (always needs both) + return ( + (eigendecomposition_requested and + config_analysis.eigendecomposition_matrices in ['exact', 'both']) or + config_analysis.enable_eigenvalue_errors + ) + + +def requires_approximate_eigendecomposition(config_analysis: AnalysisConfiguration) -> bool: + """ + Determine if approximate eigendecomposition needs to be computed. + + Approximate eigendecomposition is required for: + - Eigendecomposition analysis with eigendecomposition_matrices = 'approximate' or 'both' + - Eigenvalue error analysis (always needs both eigendecompositions) + + Parameters: + config_analysis: Analysis configuration + + Returns: + True if approximate eigendecomposition computation is needed, False otherwise + """ + # Check if eigendecomposition is requested at all + num_eigenvalues = config_analysis.num_eigenvalues + eigendecomposition_requested = ( + isinstance(num_eigenvalues, int) and num_eigenvalues > 0 + ) or ( + isinstance(num_eigenvalues, str) and num_eigenvalues.lower() == "all" + ) + + # Need approximate eigendecomposition if: + # 1. Eigendecomposition requested and matrices setting includes 'approximate' or 'both' + # 2. Eigenvalue error analysis is enabled (always needs both) + return ( + (eigendecomposition_requested and + config_analysis.eigendecomposition_matrices in ['approximate', 'both']) or + config_analysis.enable_eigenvalue_errors + ) + + +def requires_exact_matrix(config_analysis: AnalysisConfiguration) -> bool: + """ + Determine if the exact Hamiltonian matrix needs to be computed. + + The exact matrix is required for: + - Exact matrix output to file + - Exact eigendecomposition (which depends on the matrix) + - Matrix norm error analysis + - State-dependent error analysis + + Parameters: + config_analysis: Analysis configuration + + Returns: + True if exact matrix computation is needed, False otherwise + """ + return ( + config_analysis.exact_matrix_output_file is not None or + requires_exact_eigendecomposition(config_analysis) or + config_analysis.error_matrix_norms is not None or + config_analysis.error_state_inputs is not None + ) + + +def requires_approximate_matrix(config_analysis: AnalysisConfiguration) -> bool: + """ + Determine if the approximate/unitary matrix needs to be computed. + + The approximate matrix is required for: + - Matrix output to file + - Numerical simulation + - Approximate eigendecomposition (which depends on the matrix) + - Matrix norm error analysis + - State-dependent error analysis + + Parameters: + config_analysis: Analysis configuration + + Returns: + True if approximate matrix computation is needed, False otherwise + """ + return ( + config_analysis.matrix_output_file is not None or + config_analysis.numerical_simulation_inputs is not None or + requires_approximate_eigendecomposition(config_analysis) or + config_analysis.error_matrix_norms is not None or + config_analysis.error_state_inputs is not None + ) + +# ------------------------------------------------------------------------------------------------- + +def validate_and_autocomplete_analysis_config(config_analysis: AnalysisConfiguration) -> None: + """ + Validate configuration consistency and auto-enable dependent analyses where appropriate. + + This function is called early in driver.py, after loading configuration but before + loading the Hamiltonian. This allows for fail-fast behavior if configuration is invalid. + + This function checks for: + 1. Missing dependencies (raises errors if configuration is needed) + 2. Opportunities to auto-enable analyses (logs when enabling) + + Parameters: + config_analysis: Analysis configuration to validate and potentially modify + + Raises: + ValueError: If configuration is inconsistent and cannot be auto-corrected + + Note: + This function modifies the config_analysis object in-place when auto-enabling analyses. + """ + + # Check eigenvalue error analysis dependencies + if config_analysis.enable_eigenvalue_errors: + # Check if num_eigenvalues is configured + num_eigenvalues = config_analysis.num_eigenvalues + eigendecomposition_configured = ( + isinstance(num_eigenvalues, int) and num_eigenvalues > 0 + ) or ( + isinstance(num_eigenvalues, str) and num_eigenvalues.lower() == "all" + ) + + if not eigendecomposition_configured: + raise ValueError( + "enable_eigenvalue_errors requires eigendecomposition. " + "Set num_eigenvalues to a positive integer or 'all'." + ) + + # Must compute both eigendecompositions to compare + if config_analysis.eigendecomposition_matrices != 'both': + logger.info( + "INFO: enable_eigenvalue_errors requires both exact and approximate eigendecompositions. " + f"Auto-setting eigendecomposition_matrices from '{config_analysis.eigendecomposition_matrices}' to 'both'." + ) + config_analysis.eigendecomposition_matrices = 'both' + + # Check matrix norm error dependencies + if config_analysis.error_matrix_norms is not None: + # Matrix norm errors require both exact and approximate matrices + # These will be computed automatically in analyze_algorithm (checked via requires_*_matrix) + pass + + # Check state-dependent error dependencies + if config_analysis.error_state_inputs is not None: + # State errors require both exact and approximate matrices + # These will be computed automatically in analyze_algorithm (checked via requires_*_matrix) + pass + +# ------------------------------------------------------------------------------------------------- + +def analyze_algorithm( + config_analysis: AnalysisConfiguration, + algorithm, + hamiltonian=None) -> dict: + + logger.info("Beginning algorithm analysis.") + + # Note: Configuration validation happens in driver.py before Hamiltonian is loaded + + # Check what analyses are requested + eigendecomposition_requested = ( + requires_exact_eigendecomposition(config_analysis) or + requires_approximate_eigendecomposition(config_analysis) + ) error_analysis_requested = ( config_analysis.enable_eigenvalue_errors or config_analysis.error_matrix_norms is not None or config_analysis.error_state_inputs is not None ) + # Validate at least one analysis requested if (config_analysis.resource_estimator is None and config_analysis.matrix_output_file is None and config_analysis.numerical_simulation_inputs is None and @@ -925,28 +1101,24 @@ def analyze_algorithm( results = {} - # Check if any analysis needs the unitary matrix - needs_matrix = ( - config_analysis.matrix_output_file is not None or - config_analysis.numerical_simulation_inputs is not None or - (eigendecomposition_requested and - config_analysis.eigendecomposition_matrices in ['approximate', 'both']) or - (error_analysis_requested) # Error analysis always needs both - ) - - # Check if any analysis needs the exact Hamiltonian matrix - needs_exact_matrix = ( - config_analysis.exact_matrix_output_file is not None or - (eigendecomposition_requested and - config_analysis.eigendecomposition_matrices in ['exact', 'both']) or - (error_analysis_requested) # Error analysis always needs both - ) + # Determine what expensive computations are needed using shared functions + needs_matrix = requires_approximate_matrix(config_analysis) + needs_exact_matrix = requires_exact_matrix(config_analysis) # Compute matrices once if needed unitary_matrix = None if needs_matrix: unitary_matrix = _compute_unitary_matrix(algorithm) + # Opportunistic analysis: enable matrix output if not already set + if config_analysis.matrix_output_file is None: + default_filename = "unitary_matrix.npz" + logger.info( + f"INFO: Unitary matrix computed for other analyses. " + f"Auto-enabling matrix output to '{default_filename}' (essentially free)." + ) + config_analysis.matrix_output_file = default_filename + exact_matrix = None if needs_exact_matrix: if hamiltonian is None: @@ -956,6 +1128,15 @@ def analyze_algorithm( ) exact_matrix = _compute_exact_matrix(hamiltonian, config_analysis) + # Opportunistic analysis: enable exact matrix output if not already set + if config_analysis.exact_matrix_output_file is None: + default_filename = "exact_hamiltonian.npz" + logger.info( + f"INFO: Exact Hamiltonian matrix computed for other analyses. " + f"Auto-enabling exact matrix output to '{default_filename}' (essentially free)." + ) + config_analysis.exact_matrix_output_file = default_filename + # Dispatch to requested analyses if config_analysis.resource_estimator is not None: logger.info(f"Performing resource estimation using {config_analysis.resource_estimator}.") diff --git a/analysis/driver.py b/analysis/driver.py index aa6f23f..ef6216d 100644 --- a/analysis/driver.py +++ b/analysis/driver.py @@ -10,7 +10,7 @@ from qhat.common.logging_utils import configure_logging from qhat.analysis.algorithm import build_algorithm, compute_initial_phase_qubits -from qhat.analysis.analysis import analyze_algorithm +from qhat.analysis.analysis import analyze_algorithm, validate_and_autocomplete_analysis_config from qhat.analysis.configuration import load_configuration from qhat.analysis.hamiltonian import get_physical_hamiltonian from qhat.analysis.unitary import encode_as_unitary @@ -37,6 +37,12 @@ def run(): logger.info(f"Logfile: {state.config_general.logfile}") logger.info(f"Git hash: {state.config_general.git_hash}") + # Validate analysis configuration _____________________________________________________________ + + logger.info("Validating analysis configuration...") + validate_and_autocomplete_analysis_config(state.config_analysis) + logger.info("Configuration validated successfully.") + # Hamiltonian _________________________________________________________________________________ physical_hamiltonian = get_physical_hamiltonian(state.config_hamiltonian) diff --git a/analysis/tests/test_config_validation.py b/analysis/tests/test_config_validation.py new file mode 100644 index 0000000..1cad936 --- /dev/null +++ b/analysis/tests/test_config_validation.py @@ -0,0 +1,114 @@ +""" +Tests for configuration validation and autocomplete functionality. + +Tests the validate_and_autocomplete_analysis_config function which: +1. Validates configuration consistency +2. Auto-enables dependent analyses where appropriate +""" + +import pytest +from qhat.analysis.config_types import AnalysisConfiguration +from qhat.analysis.analysis import validate_and_autocomplete_analysis_config + + +def test_eigenvalue_error_requires_eigendecomposition(): + """Test that eigenvalue errors require eigendecomposition to be configured.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + # num_eigenvalues is 0 by default + + with pytest.raises(ValueError, match="enable_eigenvalue_errors requires eigendecomposition"): + validate_and_autocomplete_analysis_config(config) + + +def test_eigenvalue_error_autocorrects_eigendecomposition_matrices(): + """Test that eigenvalue errors auto-set eigendecomposition_matrices to 'both'.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'approximate' + + # Should auto-correct to 'both' + validate_and_autocomplete_analysis_config(config) + + assert config.eigendecomposition_matrices == 'both' + + +def test_eigenvalue_error_with_exact_eigendecomposition(): + """Test that eigenvalue errors auto-set eigendecomposition_matrices from 'exact' to 'both'.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'exact' + + # Should auto-correct to 'both' + validate_and_autocomplete_analysis_config(config) + + assert config.eigendecomposition_matrices == 'both' + + +def test_eigenvalue_error_already_set_to_both(): + """Test that eigenvalue errors work when eigendecomposition_matrices is already 'both'.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'both' + + # Should pass without modification + validate_and_autocomplete_analysis_config(config) + + assert config.eigendecomposition_matrices == 'both' + + +def test_eigenvalue_error_with_all_eigenvalues(): + """Test that eigenvalue errors work with num_eigenvalues='all'.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + config.num_eigenvalues = 'all' + config.eigendecomposition_matrices = 'approximate' + + # Should auto-correct to 'both' + validate_and_autocomplete_analysis_config(config) + + assert config.eigendecomposition_matrices == 'both' + + +def test_matrix_norm_errors_pass_validation(): + """Test that matrix norm errors pass validation (matrices computed automatically).""" + config = AnalysisConfiguration() + config.error_matrix_norms = 'frobenius' + + # Should pass - matrices will be computed automatically in analyze_algorithm + validate_and_autocomplete_analysis_config(config) + + +def test_state_errors_pass_validation(): + """Test that state errors pass validation (matrices computed automatically).""" + config = AnalysisConfiguration() + config.error_state_inputs = 'state.npy' + + # Should pass - matrices will be computed automatically in analyze_algorithm + validate_and_autocomplete_analysis_config(config) + + +def test_multiple_error_types_with_eigenvalue(): + """Test combining eigenvalue errors with other error types.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'approximate' + config.error_matrix_norms = 'frobenius' + config.error_state_inputs = 'state.npy' + + # Should auto-correct eigendecomposition_matrices to 'both' + validate_and_autocomplete_analysis_config(config) + + assert config.eigendecomposition_matrices == 'both' + + +def test_no_analyses_configured(): + """Test that validation passes with no analyses configured (will fail later in analyze_algorithm).""" + config = AnalysisConfiguration() + + # Should pass - the "no analyses requested" check happens in analyze_algorithm + validate_and_autocomplete_analysis_config(config) diff --git a/analysis/tests/test_driver_validation.py b/analysis/tests/test_driver_validation.py new file mode 100644 index 0000000..6833351 --- /dev/null +++ b/analysis/tests/test_driver_validation.py @@ -0,0 +1,141 @@ +""" +Test that configuration validation runs early in the driver, before Hamiltonian loading. + +This ensures fail-fast behavior when configuration is invalid. +""" + +import pytest +import tempfile +import os +from unittest.mock import patch, MagicMock + + +def test_validation_called_before_hamiltonian_load(): + """Test that validation is called before get_physical_hamiltonian.""" + from qhat.analysis import driver + + # Track the order of calls + call_order = [] + + # Mock the functions we care about + with patch('qhat.analysis.driver.load_configuration') as mock_load_config, \ + patch('qhat.analysis.driver.configure_logging') as mock_configure_logging, \ + patch('qhat.analysis.driver.validate_and_autocomplete_analysis_config') as mock_validate, \ + patch('qhat.analysis.driver.get_physical_hamiltonian') as mock_get_hamiltonian, \ + patch('qhat.analysis.driver.encode_as_unitary') as mock_encode_unitary, \ + patch('qhat.analysis.driver.build_algorithm') as mock_build_algorithm, \ + patch('qhat.analysis.driver.analyze_algorithm') as mock_analyze_algorithm: + + # Create a mock state object + mock_state = MagicMock() + mock_state.config_general.loglevel = 'info' + mock_state.config_general.logfile = 'test.log' + mock_state.config_general.git_hash = 'test123' + mock_state.config_unitary.method = 'not ramped trotter' # Skip Trotter computation + mock_load_config.return_value = mock_state + + # Track call order + def track_validate(*args, **kwargs): + call_order.append('validate') + + def track_hamiltonian(*args, **kwargs): + call_order.append('hamiltonian') + return MagicMock() + + mock_validate.side_effect = track_validate + mock_get_hamiltonian.side_effect = track_hamiltonian + + # Mock other functions to return reasonable values + mock_encode_unitary.return_value = MagicMock() + mock_build_algorithm.return_value = MagicMock() + mock_analyze_algorithm.return_value = {} + + # Run the driver + driver.run() + + # Verify validation was called before Hamiltonian loading + assert call_order == ['validate', 'hamiltonian'], \ + f"Expected ['validate', 'hamiltonian'], got {call_order}" + + +def test_validation_error_prevents_hamiltonian_load(): + """Test that validation errors prevent expensive Hamiltonian loading.""" + from qhat.analysis import driver + from qhat.analysis.config_types import AnalysisConfiguration + + # Track whether Hamiltonian was loaded + hamiltonian_loaded = [] + + with patch('qhat.analysis.driver.load_configuration') as mock_load_config, \ + patch('qhat.analysis.driver.configure_logging') as mock_configure_logging, \ + patch('qhat.analysis.driver.get_physical_hamiltonian') as mock_get_hamiltonian: + + # Create a mock state with invalid configuration + mock_state = MagicMock() + mock_state.config_general.loglevel = 'info' + mock_state.config_general.logfile = 'test.log' + mock_state.config_general.git_hash = 'test123' + mock_state.config_unitary.method = 'not ramped trotter' # Avoid Trotter code path + + # Create invalid config: eigenvalue errors without eigendecomposition + invalid_config = AnalysisConfiguration() + invalid_config.enable_eigenvalue_errors = True + invalid_config.num_eigenvalues = 0 # Invalid! + mock_state.config_analysis = invalid_config + + mock_load_config.return_value = mock_state + + # Track if Hamiltonian loading is attempted + def track_hamiltonian(*args, **kwargs): + hamiltonian_loaded.append(True) + return MagicMock() + + mock_get_hamiltonian.side_effect = track_hamiltonian + + # Run should fail with ValueError during validation + with pytest.raises(ValueError, match="enable_eigenvalue_errors requires eigendecomposition"): + driver.run() + + # Hamiltonian should NOT have been loaded + assert len(hamiltonian_loaded) == 0, \ + "Hamiltonian was loaded despite validation error (should fail fast)" + + +def test_validation_autocorrects_config(): + """Test that validation auto-corrects configuration as expected.""" + from qhat.analysis import driver + from qhat.analysis.config_types import AnalysisConfiguration + + with patch('qhat.analysis.driver.load_configuration') as mock_load_config, \ + patch('qhat.analysis.driver.configure_logging') as mock_configure_logging, \ + patch('qhat.analysis.driver.get_physical_hamiltonian') as mock_get_hamiltonian, \ + patch('qhat.analysis.driver.encode_as_unitary') as mock_encode_unitary, \ + patch('qhat.analysis.driver.build_algorithm') as mock_build_algorithm, \ + patch('qhat.analysis.driver.analyze_algorithm') as mock_analyze_algorithm: + + # Create config that needs auto-correction + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'approximate' # Should be auto-corrected to 'both' + + mock_state = MagicMock() + mock_state.config_general.loglevel = 'info' + mock_state.config_general.logfile = 'test.log' + mock_state.config_general.git_hash = 'test123' + mock_state.config_analysis = config + mock_state.config_unitary.method = 'not ramped trotter' + mock_load_config.return_value = mock_state + + # Mock other functions + mock_get_hamiltonian.return_value = MagicMock() + mock_encode_unitary.return_value = MagicMock() + mock_build_algorithm.return_value = MagicMock() + mock_analyze_algorithm.return_value = {} + + # Run the driver + driver.run() + + # Verify config was auto-corrected + assert config.eigendecomposition_matrices == 'both', \ + f"Expected 'both', got '{config.eigendecomposition_matrices}'" diff --git a/analysis/tests/test_expensive_analysis_detection.py b/analysis/tests/test_expensive_analysis_detection.py new file mode 100644 index 0000000..b682538 --- /dev/null +++ b/analysis/tests/test_expensive_analysis_detection.py @@ -0,0 +1,268 @@ +""" +Tests for the shared functions that detect when expensive analyses are required. + +These functions are used by both validation and analyze_algorithm to ensure +consistent detection of what computations are needed. +""" + +import pytest +from qhat.analysis.config_types import AnalysisConfiguration +from qhat.analysis.analysis import ( + requires_exact_eigendecomposition, + requires_approximate_eigendecomposition, + requires_exact_matrix, + requires_approximate_matrix +) + + +# ================================================================================================= +# Tests for requires_exact_eigendecomposition +# ================================================================================================= + +def test_exact_eigendecomposition_not_required_by_default(): + """Test that exact eigendecomposition is not required by default.""" + config = AnalysisConfiguration() + assert not requires_exact_eigendecomposition(config) + + +def test_exact_eigendecomposition_required_when_matrices_exact(): + """Test that exact eigendecomposition is required when eigendecomposition_matrices='exact'.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'exact' + assert requires_exact_eigendecomposition(config) + + +def test_exact_eigendecomposition_required_when_matrices_both(): + """Test that exact eigendecomposition is required when eigendecomposition_matrices='both'.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'both' + assert requires_exact_eigendecomposition(config) + + +def test_exact_eigendecomposition_not_required_when_matrices_approximate(): + """Test that exact eigendecomposition is not required when eigendecomposition_matrices='approximate'.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'approximate' + assert not requires_exact_eigendecomposition(config) + + +def test_exact_eigendecomposition_required_for_eigenvalue_errors(): + """Test that exact eigendecomposition is required when eigenvalue errors are enabled.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + # Note: In practice, validation would require num_eigenvalues to be set, + # but the function itself should return True based on enable_eigenvalue_errors alone + assert requires_exact_eigendecomposition(config) + + +def test_exact_eigendecomposition_with_all_eigenvalues(): + """Test that exact eigendecomposition works with num_eigenvalues='all'.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 'all' + config.eigendecomposition_matrices = 'exact' + assert requires_exact_eigendecomposition(config) + + +# ================================================================================================= +# Tests for requires_approximate_eigendecomposition +# ================================================================================================= + +def test_approximate_eigendecomposition_not_required_by_default(): + """Test that approximate eigendecomposition is not required by default.""" + config = AnalysisConfiguration() + assert not requires_approximate_eigendecomposition(config) + + +def test_approximate_eigendecomposition_required_when_matrices_approximate(): + """Test that approximate eigendecomposition is required when eigendecomposition_matrices='approximate'.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'approximate' + assert requires_approximate_eigendecomposition(config) + + +def test_approximate_eigendecomposition_required_when_matrices_both(): + """Test that approximate eigendecomposition is required when eigendecomposition_matrices='both'.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'both' + assert requires_approximate_eigendecomposition(config) + + +def test_approximate_eigendecomposition_not_required_when_matrices_exact(): + """Test that approximate eigendecomposition is not required when eigendecomposition_matrices='exact'.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'exact' + assert not requires_approximate_eigendecomposition(config) + + +def test_approximate_eigendecomposition_required_for_eigenvalue_errors(): + """Test that approximate eigendecomposition is required when eigenvalue errors are enabled.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + assert requires_approximate_eigendecomposition(config) + + +# ================================================================================================= +# Tests for requires_exact_matrix +# ================================================================================================= + +def test_exact_matrix_not_required_by_default(): + """Test that exact matrix is not required by default.""" + config = AnalysisConfiguration() + assert not requires_exact_matrix(config) + + +def test_exact_matrix_required_for_output_file(): + """Test that exact matrix is required when exact_matrix_output_file is set.""" + config = AnalysisConfiguration() + config.exact_matrix_output_file = 'exact.npz' + assert requires_exact_matrix(config) + + +def test_exact_matrix_required_for_exact_eigendecomposition(): + """Test that exact matrix is required when exact eigendecomposition is needed.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'exact' + assert requires_exact_matrix(config) + + +def test_exact_matrix_required_for_both_eigendecomposition(): + """Test that exact matrix is required when both eigendecompositions are needed.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'both' + assert requires_exact_matrix(config) + + +def test_exact_matrix_required_for_matrix_norm_errors(): + """Test that exact matrix is required for matrix norm error analysis.""" + config = AnalysisConfiguration() + config.error_matrix_norms = 'frobenius' + assert requires_exact_matrix(config) + + +def test_exact_matrix_required_for_state_errors(): + """Test that exact matrix is required for state-dependent error analysis.""" + config = AnalysisConfiguration() + config.error_state_inputs = 'state.npy' + assert requires_exact_matrix(config) + + +def test_exact_matrix_required_for_eigenvalue_errors(): + """Test that exact matrix is required for eigenvalue error analysis (via eigendecomposition).""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + assert requires_exact_matrix(config) + + +# ================================================================================================= +# Tests for requires_approximate_matrix +# ================================================================================================= + +def test_approximate_matrix_not_required_by_default(): + """Test that approximate matrix is not required by default.""" + config = AnalysisConfiguration() + assert not requires_approximate_matrix(config) + + +def test_approximate_matrix_required_for_output_file(): + """Test that approximate matrix is required when matrix_output_file is set.""" + config = AnalysisConfiguration() + config.matrix_output_file = 'matrix.npz' + assert requires_approximate_matrix(config) + + +def test_approximate_matrix_required_for_numerical_simulation(): + """Test that approximate matrix is required for numerical simulation.""" + config = AnalysisConfiguration() + config.numerical_simulation_inputs = 'state.npy' + assert requires_approximate_matrix(config) + + +def test_approximate_matrix_required_for_approximate_eigendecomposition(): + """Test that approximate matrix is required when approximate eigendecomposition is needed.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'approximate' + assert requires_approximate_matrix(config) + + +def test_approximate_matrix_required_for_both_eigendecomposition(): + """Test that approximate matrix is required when both eigendecompositions are needed.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'both' + assert requires_approximate_matrix(config) + + +def test_approximate_matrix_required_for_matrix_norm_errors(): + """Test that approximate matrix is required for matrix norm error analysis.""" + config = AnalysisConfiguration() + config.error_matrix_norms = 'frobenius' + assert requires_approximate_matrix(config) + + +def test_approximate_matrix_required_for_state_errors(): + """Test that approximate matrix is required for state-dependent error analysis.""" + config = AnalysisConfiguration() + config.error_state_inputs = 'state.npy' + assert requires_approximate_matrix(config) + + +def test_approximate_matrix_required_for_eigenvalue_errors(): + """Test that approximate matrix is required for eigenvalue error analysis (via eigendecomposition).""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + assert requires_approximate_matrix(config) + + +# ================================================================================================= +# Integration tests: verify consistency +# ================================================================================================= + +def test_eigenvalue_errors_requires_both_matrices(): + """Test that eigenvalue errors correctly require both matrices.""" + config = AnalysisConfiguration() + config.enable_eigenvalue_errors = True + + assert requires_exact_matrix(config), "Eigenvalue errors should require exact matrix" + assert requires_approximate_matrix(config), "Eigenvalue errors should require approximate matrix" + + +def test_matrix_norm_errors_requires_both_matrices(): + """Test that matrix norm errors correctly require both matrices.""" + config = AnalysisConfiguration() + config.error_matrix_norms = ['frobenius', 'spectral'] + + assert requires_exact_matrix(config), "Matrix norm errors should require exact matrix" + assert requires_approximate_matrix(config), "Matrix norm errors should require approximate matrix" + + +def test_state_errors_requires_both_matrices(): + """Test that state errors correctly require both matrices.""" + config = AnalysisConfiguration() + config.error_state_inputs = ['state1.npy', 'state2.npy'] + + assert requires_exact_matrix(config), "State errors should require exact matrix" + assert requires_approximate_matrix(config), "State errors should require approximate matrix" + + +def test_eigendecomposition_dependency_on_matrices(): + """Test that eigendecomposition correctly implies matrix computation.""" + config = AnalysisConfiguration() + config.num_eigenvalues = 5 + config.eigendecomposition_matrices = 'both' + + # Both eigendecompositions should be required + assert requires_exact_eigendecomposition(config) + assert requires_approximate_eigendecomposition(config) + + # And both matrices should be required (because eigendecomposition needs them) + assert requires_exact_matrix(config) + assert requires_approximate_matrix(config) From 9adc95c48a893d19c8dff7cafad6679c2132ccad Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 15 Jun 2026 10:21:22 -0600 Subject: [PATCH 06/12] Same Claude pattern I keep telling it to remove. --- analysis/analysis.py | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/analysis/analysis.py b/analysis/analysis.py index c678ef4..3b22371 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -558,16 +558,18 @@ def error_analysis( logger.info(f"Computing matrix norm errors: {norms_to_compute}") - # Get matrices if not provided + # Require matrices to be provided by caller if exact_matrix is None: - if hamiltonian is None: - raise ValueError("hamiltonian required for matrix norm error analysis") - exact_matrix = _compute_exact_matrix(hamiltonian) + raise ValueError( + "Matrix norm error analysis requires the exact Hamiltonian matrix, but it was not computed. " + "This is an internal error - matrix should have been computed based on your configuration." + ) if unitary_matrix is None: - if algorithm is None: - raise ValueError("algorithm required for matrix norm error analysis") - unitary_matrix = _compute_unitary_matrix(algorithm) + raise ValueError( + "Matrix norm error analysis requires the approximate/unitary matrix, but it was not computed. " + "This is an internal error - matrix should have been computed based on your configuration." + ) # Check if matrices are dense or matrix-free is_exact_dense = isinstance(exact_matrix, np.ndarray) @@ -720,16 +722,18 @@ def error_analysis( logger.info(f"Computing state-dependent errors for {len(state_files)} state(s)") - # Get matrices/operators if not provided + # Require matrices to be provided by caller if exact_matrix is None: - if hamiltonian is None: - raise ValueError("hamiltonian required for state-dependent error analysis") - exact_matrix = _compute_exact_matrix(hamiltonian) + raise ValueError( + "State-dependent error analysis requires the exact Hamiltonian matrix, but it was not computed. " + "This is an internal error - matrix should have been computed based on your configuration." + ) if unitary_matrix is None: - if algorithm is None: - raise ValueError("algorithm required for state-dependent error analysis") - unitary_matrix = _compute_unitary_matrix(algorithm) + raise ValueError( + "State-dependent error analysis requires the approximate/unitary matrix, but it was not computed. " + "This is an internal error - matrix should have been computed based on your configuration." + ) state_errors = [] From 979d7bc90c6da7c61ba52eedf085f4651eb09219 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 15 Jun 2026 10:28:27 -0600 Subject: [PATCH 07/12] A little cleanup --- analysis/analysis.py | 79 ++++++++++++++++----------- analysis/tests/test_error_analysis.py | 16 +++++- 2 files changed, 62 insertions(+), 33 deletions(-) diff --git a/analysis/analysis.py b/analysis/analysis.py index 3b22371..cf9482e 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -550,11 +550,8 @@ def error_analysis( # ============================================================================================= if config_analysis.error_matrix_norms is not None: - # Normalize to list - if isinstance(config_analysis.error_matrix_norms, str): - norms_to_compute = [config_analysis.error_matrix_norms] - else: - norms_to_compute = config_analysis.error_matrix_norms + # Note: error_matrix_norms is normalized to list during validation + norms_to_compute = config_analysis.error_matrix_norms logger.info(f"Computing matrix norm errors: {norms_to_compute}") @@ -714,11 +711,8 @@ def error_analysis( # ============================================================================================= if config_analysis.error_state_inputs is not None: - # Normalize to list - if isinstance(config_analysis.error_state_inputs, str): - state_files = [config_analysis.error_state_inputs] - else: - state_files = config_analysis.error_state_inputs + # Note: error_state_inputs is normalized to list during validation + state_files = config_analysis.error_state_inputs logger.info(f"Computing state-dependent errors for {len(state_files)} state(s)") @@ -824,17 +818,8 @@ def numerical_simulation( # Log matrix properties logger.verbose(f"Matrix shape: {unitary_matrix.shape}") - # Normalize input to list - inputs = config_analysis.numerical_simulation_inputs - if isinstance(inputs, str): - input_files = [inputs] - elif isinstance(inputs, list): - input_files = inputs - else: - raise ValueError( - f"numerical_simulation_inputs must be a string or list of strings, " - f"got {type(inputs)}" - ) + # Note: numerical_simulation_inputs is normalized to list during validation + input_files = config_analysis.numerical_simulation_inputs logger.info(f"Running numerical simulation on {len(input_files)} input state(s)") @@ -886,6 +871,36 @@ def numerical_simulation( return {'simulations': results} +# ------------------------------------------------------------------------------------------------- +# Helper functions +# ------------------------------------------------------------------------------------------------- + +def _normalize_string_or_list_to_list(value): + """ + Normalize a configuration value that can be either a string or list into a list. + + This is a common pattern for config options that accept either a single item (string) + or multiple items (list). + + Parameters: + value: Either a string or a list of strings (or None) + + Returns: + A list (or None if input was None) + + Examples: + _normalize_string_or_list_to_list("item") -> ["item"] + _normalize_string_or_list_to_list(["a", "b"]) -> ["a", "b"] + _normalize_string_or_list_to_list(None) -> None + """ + if value is None: + return None + elif isinstance(value, str): + return [value] + else: + # Already a list (or list-like) + return value + # ------------------------------------------------------------------------------------------------- # Functions to determine what expensive computations are required # ------------------------------------------------------------------------------------------------- @@ -1050,17 +1065,17 @@ def validate_and_autocomplete_analysis_config(config_analysis: AnalysisConfigura ) config_analysis.eigendecomposition_matrices = 'both' - # Check matrix norm error dependencies - if config_analysis.error_matrix_norms is not None: - # Matrix norm errors require both exact and approximate matrices - # These will be computed automatically in analyze_algorithm (checked via requires_*_matrix) - pass - - # Check state-dependent error dependencies - if config_analysis.error_state_inputs is not None: - # State errors require both exact and approximate matrices - # These will be computed automatically in analyze_algorithm (checked via requires_*_matrix) - pass + # Normalize string-or-list config values to always be lists + # This allows downstream code to always assume list type + config_analysis.error_matrix_norms = _normalize_string_or_list_to_list( + config_analysis.error_matrix_norms + ) + config_analysis.error_state_inputs = _normalize_string_or_list_to_list( + config_analysis.error_state_inputs + ) + config_analysis.numerical_simulation_inputs = _normalize_string_or_list_to_list( + config_analysis.numerical_simulation_inputs + ) # ------------------------------------------------------------------------------------------------- diff --git a/analysis/tests/test_error_analysis.py b/analysis/tests/test_error_analysis.py index 0a50e11..75b0da9 100644 --- a/analysis/tests/test_error_analysis.py +++ b/analysis/tests/test_error_analysis.py @@ -15,7 +15,7 @@ import os from pathlib import Path -from qhat.analysis.analysis import error_analysis +from qhat.analysis.analysis import error_analysis, validate_and_autocomplete_analysis_config from qhat.analysis.config_types import AnalysisConfiguration from qhat.analysis.file_io import save_state @@ -139,6 +139,7 @@ def test_frobenius_norm_zero_when_identical(): config = AnalysisConfiguration() config.error_matrix_norms = 'frobenius' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -166,6 +167,7 @@ def test_frobenius_norm_nonzero_when_different(): config = AnalysisConfiguration() config.error_matrix_norms = 'frobenius' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -197,6 +199,7 @@ def test_spectral_norm_zero_when_identical(): config = AnalysisConfiguration() config.error_matrix_norms = 'spectral' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -224,6 +227,7 @@ def test_spectral_norm_nonzero_when_different(): config = AnalysisConfiguration() config.error_matrix_norms = 'spectral' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -289,6 +293,7 @@ def test_state_error_zero_when_identical(): config = AnalysisConfiguration() config.error_state_inputs = 'test_state.npy' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -321,6 +326,7 @@ def test_state_error_nonzero_when_different(): config = AnalysisConfiguration() config.error_state_inputs = 'test_state.npy' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -406,8 +412,11 @@ def test_all_error_types_together(): config = AnalysisConfiguration() config.enable_eigenvalue_errors = True + config.num_eigenvalues = 2 # Required for eigenvalue errors + config.eigendecomposition_matrices = 'both' # Will be auto-set by validation, but explicit here config.error_matrix_norms = ['frobenius', 'spectral'] config.error_state_inputs = 'test_state.npy' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -513,6 +522,7 @@ def test_invalid_matrix_norm_type(): config = AnalysisConfiguration() config.error_matrix_norms = 'invalid' + validate_and_autocomplete_analysis_config(config) # Normalize config values with pytest.raises(ValueError, match="Unknown matrix norm type"): with tempfile.TemporaryDirectory() as tmpdir: @@ -536,6 +546,7 @@ def test_missing_state_file(): config = AnalysisConfiguration() config.error_state_inputs = 'nonexistent.npy' + validate_and_autocomplete_analysis_config(config) # Normalize config values with pytest.raises(Exception): # Will be FileNotFoundError or similar with tempfile.TemporaryDirectory() as tmpdir: @@ -570,6 +581,7 @@ def test_state_error_with_matrix_free_operators(): config = AnalysisConfiguration() config.error_state_inputs = 'test_state.npy' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -607,6 +619,7 @@ def test_frobenius_norm_with_matrix_free_small(): config = AnalysisConfiguration() config.error_matrix_norms = 'frobenius' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() @@ -690,6 +703,7 @@ def test_state_relative_error(): config = AnalysisConfiguration() config.error_state_inputs = 'test_state.npy' + validate_and_autocomplete_analysis_config(config) # Normalize config values with tempfile.TemporaryDirectory() as tmpdir: original_dir = os.getcwd() From a24e166f212cb7c9ebb3cee7f4515a8afe2ae0ad Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 15 Jun 2026 11:02:16 -0600 Subject: [PATCH 08/12] Fix --- analysis/analysis.py | 42 +++++++++++----------- analysis/tests/test_config_validation.py | 45 ++++++++++++++++++++++++ 2 files changed, 65 insertions(+), 22 deletions(-) diff --git a/analysis/analysis.py b/analysis/analysis.py index cf9482e..ac541d8 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -559,13 +559,11 @@ def error_analysis( if exact_matrix is None: raise ValueError( "Matrix norm error analysis requires the exact Hamiltonian matrix, but it was not computed. " - "This is an internal error - matrix should have been computed based on your configuration." ) if unitary_matrix is None: raise ValueError( "Matrix norm error analysis requires the approximate/unitary matrix, but it was not computed. " - "This is an internal error - matrix should have been computed based on your configuration." ) # Check if matrices are dense or matrix-free @@ -720,13 +718,11 @@ def error_analysis( if exact_matrix is None: raise ValueError( "State-dependent error analysis requires the exact Hamiltonian matrix, but it was not computed. " - "This is an internal error - matrix should have been computed based on your configuration." ) if unitary_matrix is None: raise ValueError( "State-dependent error analysis requires the approximate/unitary matrix, but it was not computed. " - "This is an internal error - matrix should have been computed based on your configuration." ) state_errors = [] @@ -1077,6 +1073,25 @@ def validate_and_autocomplete_analysis_config(config_analysis: AnalysisConfigura config_analysis.numerical_simulation_inputs ) + # Check if matrices will be computed and auto-enable output if not already set + if requires_approximate_matrix(config_analysis): + if config_analysis.matrix_output_file is None: + default_filename = "unitary_matrix.npz" + logger.info( + f"INFO: Approximate/unitary matrix will be computed for requested analyses. " + f"Auto-enabling matrix output to '{default_filename}' (essentially free)." + ) + config_analysis.matrix_output_file = default_filename + + if requires_exact_matrix(config_analysis): + if config_analysis.exact_matrix_output_file is None: + default_filename = "exact_hamiltonian.npz" + logger.info( + f"INFO: Exact Hamiltonian matrix will be computed for requested analyses. " + f"Auto-enabling exact matrix output to '{default_filename}' (essentially free)." + ) + config_analysis.exact_matrix_output_file = default_filename + # ------------------------------------------------------------------------------------------------- def analyze_algorithm( @@ -1125,19 +1140,11 @@ def analyze_algorithm( needs_exact_matrix = requires_exact_matrix(config_analysis) # Compute matrices once if needed + # Note: Opportunistic matrix output enabling happens during validation in driver.py unitary_matrix = None if needs_matrix: unitary_matrix = _compute_unitary_matrix(algorithm) - # Opportunistic analysis: enable matrix output if not already set - if config_analysis.matrix_output_file is None: - default_filename = "unitary_matrix.npz" - logger.info( - f"INFO: Unitary matrix computed for other analyses. " - f"Auto-enabling matrix output to '{default_filename}' (essentially free)." - ) - config_analysis.matrix_output_file = default_filename - exact_matrix = None if needs_exact_matrix: if hamiltonian is None: @@ -1147,15 +1154,6 @@ def analyze_algorithm( ) exact_matrix = _compute_exact_matrix(hamiltonian, config_analysis) - # Opportunistic analysis: enable exact matrix output if not already set - if config_analysis.exact_matrix_output_file is None: - default_filename = "exact_hamiltonian.npz" - logger.info( - f"INFO: Exact Hamiltonian matrix computed for other analyses. " - f"Auto-enabling exact matrix output to '{default_filename}' (essentially free)." - ) - config_analysis.exact_matrix_output_file = default_filename - # Dispatch to requested analyses if config_analysis.resource_estimator is not None: logger.info(f"Performing resource estimation using {config_analysis.resource_estimator}.") diff --git a/analysis/tests/test_config_validation.py b/analysis/tests/test_config_validation.py index 1cad936..0acfd2b 100644 --- a/analysis/tests/test_config_validation.py +++ b/analysis/tests/test_config_validation.py @@ -112,3 +112,48 @@ def test_no_analyses_configured(): # Should pass - the "no analyses requested" check happens in analyze_algorithm validate_and_autocomplete_analysis_config(config) + + +def test_opportunistic_matrix_output_for_numerical_simulation(): + """Test that matrix output is auto-enabled when matrices will be computed.""" + config = AnalysisConfiguration() + config.numerical_simulation_inputs = 'state.npy' + + # Before validation + assert config.matrix_output_file is None + + # Validate - should auto-enable matrix output + validate_and_autocomplete_analysis_config(config) + + # After validation + assert config.matrix_output_file == 'unitary_matrix.npz' + + +def test_opportunistic_exact_matrix_output_for_error_analysis(): + """Test that exact matrix output is auto-enabled for error analyses.""" + config = AnalysisConfiguration() + config.error_matrix_norms = 'frobenius' + + # Before validation + assert config.exact_matrix_output_file is None + assert config.matrix_output_file is None + + # Validate - should auto-enable both matrix outputs + validate_and_autocomplete_analysis_config(config) + + # After validation + assert config.exact_matrix_output_file == 'exact_hamiltonian.npz' + assert config.matrix_output_file == 'unitary_matrix.npz' + + +def test_no_opportunistic_enabling_when_already_set(): + """Test that opportunistic enabling respects existing settings.""" + config = AnalysisConfiguration() + config.numerical_simulation_inputs = 'state.npy' + config.matrix_output_file = 'my_custom_name.npz' + + # Validate - should NOT change the custom filename + validate_and_autocomplete_analysis_config(config) + + # Should keep custom name + assert config.matrix_output_file == 'my_custom_name.npz' From b131a63177309c1238c45fadbc81e053f1ace9f4 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Tue, 7 Jul 2026 15:23:05 -0600 Subject: [PATCH 09/12] Fix input file issue. --- analysis/analysis.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/analysis/analysis.py b/analysis/analysis.py index ac541d8..137eb6b 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -814,8 +814,19 @@ def numerical_simulation( # Log matrix properties logger.verbose(f"Matrix shape: {unitary_matrix.shape}") - # Note: numerical_simulation_inputs is normalized to list during validation - input_files = config_analysis.numerical_simulation_inputs + # Normalize inputs to list (in case this function is called directly without validation) + raw_inputs = config_analysis.numerical_simulation_inputs + if raw_inputs is None: + raise ValueError("numerical_simulation_inputs is None") + + # Validate type before normalization + if not isinstance(raw_inputs, (str, list)): + raise ValueError( + f"numerical_simulation_inputs must be a string or list of strings, " + f"got {type(raw_inputs).__name__}" + ) + + input_files = _normalize_string_or_list_to_list(raw_inputs) logger.info(f"Running numerical simulation on {len(input_files)} input state(s)") From ff1099cffe7e2654f1830810cf2a053bf0ffbd7f Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 28 May 2026 16:07:48 -0600 Subject: [PATCH 10/12] Branch 2: Add eigendecomposition analysis - Add eigendecomposition_analysis() function in analysis.py - Support full and partial eigendecomposition modes - Add configuration fields: num_eigenvalues, eigendecomposition_matrices, which_eigenvalues - Support 'smallest', 'largest', and 'both' eigenvalue selection - Add save_eigendecomposition() and load_eigendecomposition() to file_io.py - Add 19 comprehensive tests in test_eigendecomposition.py - Update README.md with eigendecomposition documentation - Update config_full_analysis.py with eigendecomposition examples - All 200 tests passing (128 original + 53 Branch 1 + 19 Branch 2) --- analysis/examples/config_full_analysis.py | 68 +++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/analysis/examples/config_full_analysis.py b/analysis/examples/config_full_analysis.py index e7af219..b8cd558 100644 --- a/analysis/examples/config_full_analysis.py +++ b/analysis/examples/config_full_analysis.py @@ -347,10 +347,78 @@ # analysis.error_state_inputs = "examples/ground.npy" # Scales well # # Skip error_matrix_norms for large systems (too slow) +# ================================================================================================= +# 5. EIGENDECOMPOSITION ANALYSIS (Branch 2: Now Available!) +# ================================================================================================= +# Compute eigenvalues and eigenvectors of exact and/or approximate matrices + +analysis.num_eigenvalues = 5 # Compute 5 eigenvalues (use sparse methods) +analysis.eigendecomposition_matrices = "both" # Options: "exact", "approximate", "both" +analysis.which_eigenvalues = "smallest" # Options: "smallest", "largest", "both" + +# This analysis computes eigendecompositions for validation and comparison: +# - Ground state energy (use num_eigenvalues=1, which_eigenvalues="smallest") +# - Low-lying excited states (use num_eigenvalues=5, which_eigenvalues="smallest") +# - High-energy states (use num_eigenvalues=5, which_eigenvalues="largest") +# - Both ends of spectrum (use which_eigenvalues="both") + +# Configuration options: +# num_eigenvalues: +# - 0 (default): Eigendecomposition disabled +# - Positive int (e.g., 5): Compute that many eigenvalues (recommended for large systems) +# - "all": Full eigendecomposition (only for systems ≤10 qubits) +# +# eigendecomposition_matrices: +# - "approximate" (default): Only eigendecompose unitary matrix +# - "exact": Only eigendecompose exact Hamiltonian matrix +# - "both": Eigendecompose both (useful for comparison) +# +# which_eigenvalues (ignored for "all"): +# - "smallest" (default): Algebraically smallest (most negative, ground state) +# - "largest": Algebraically largest (most positive) +# - "both": Compute smallest AND largest (returns 2k eigenvalues) +# +# System size guidance: +# - ≤10 qubits: "all" feasible +# - 10-12 qubits: "all" possible but expensive +# - ≥12 qubits: Use partial (specify k=5-10) +# - 20+ qubits: MUST use partial +# - 25-30 qubits: Partial still scales well + +# Output files: +# - exact_eigendecomposition.npz (if "exact" or "both") +# - approximate_eigendecomposition.npz (if "approximate" or "both") + +# Examples: + +# Ground state energy comparison (most common use case for large systems) +# analysis.num_eigenvalues = 1 +# analysis.eigendecomposition_matrices = "both" +# analysis.which_eigenvalues = "smallest" + +# Low-lying spectrum (5 lowest-energy states) +# analysis.num_eigenvalues = 5 +# analysis.eigendecomposition_matrices = "exact" +# analysis.which_eigenvalues = "smallest" + +# High and low energy states (3 smallest + 3 largest = 6 total) +# analysis.num_eigenvalues = 3 +# analysis.eigendecomposition_matrices = "both" +# analysis.which_eigenvalues = "both" + +# Full spectrum for small system +# analysis.num_eigenvalues = "all" +# analysis.eigendecomposition_matrices = "exact" + # ================================================================================================= # FUTURE ANALYSES (Coming in subsequent branches) # ================================================================================================= +# Branch 3: Error Metrics +analysis.error_num_eigenvalues = 1 # Compare ground state energy +analysis.error_matrix_norms = ["frobenius", "spectral"] # Matrix difference norms +analysis.error_state_inputs = ["examples/ground_state.npy"] # State-dependent errors + # Branch 4: Exact Numerical Simulation # analysis.exact_simulation_inputs = "examples/initial_state.npy" # Simulate using exact Hamiltonian evolution (for comparison with approximate) From 550da60014091bb74f4ca3cb67040d2f173e527e Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 28 May 2026 16:38:30 -0600 Subject: [PATCH 11/12] Branch 3: Add error analysis (eigenvalue, matrix norm, state-dependent) - Add error_analysis() function in analysis.py (~350 lines) - Three independent error types: 1. Eigenvalue errors: Compare k smallest eigenvalues 2. Matrix norm errors: Frobenius and spectral norms 3. State-dependent errors: Apply operators to specific states - Add configuration fields: error_num_eigenvalues, error_matrix_norms, error_state_inputs - Support matrix-free computation for large systems - Matrix norm computation with progress tracking - Add 18 comprehensive tests in test_error_analysis.py - Update README.md with error analysis documentation - Update config_full_analysis.py with error analysis examples - All 218 tests passing (128 original + 53 Branch 1 + 19 Branch 2 + 18 Branch 3) --- analysis/analysis.py | 350 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 350 insertions(+) diff --git a/analysis/analysis.py b/analysis/analysis.py index 137eb6b..9621ec7 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -795,6 +795,356 @@ def error_analysis( # ------------------------------------------------------------------------------------------------- +def error_analysis( + config_analysis: AnalysisConfiguration, + hamiltonian, + algorithm, + exact_matrix=None, + unitary_matrix=None, + exact_eigendecomp=None, + approx_eigendecomp=None) -> dict: + """ + Compute error metrics comparing exact and approximate representations. + + Three independent error types: + 1. Eigenvalue errors (if error_num_eigenvalues > 0) + 2. Matrix norm errors (if error_matrix_norms is not None) + 3. State-dependent errors (if error_state_inputs is not None) + + Parameters: + config_analysis: Analysis configuration with error analysis settings + hamiltonian: Hamiltonian object + algorithm: Algorithm bloq + exact_matrix: Pre-computed exact matrix (optional) + unitary_matrix: Pre-computed unitary matrix (optional) + exact_eigendecomp: Pre-computed exact eigendecomposition (optional) + approx_eigendecomp: Pre-computed approximate eigendecomposition (optional) + + Returns: + Dictionary with error metrics + """ + from qhat.analysis.matrix_operations import PauliStringOperator + from qhat.analysis.file_io import load_eigendecomposition, load_state + import scipy.linalg + + logger.info("Starting error analysis") + + results = {} + + # ============================================================================================= + # 1. EIGENVALUE ERROR + # ============================================================================================= + + if config_analysis.error_num_eigenvalues > 0: + logger.info(f"Computing eigenvalue errors for {config_analysis.error_num_eigenvalues} eigenvalues") + + # Load or compute eigendecompositions + if exact_eigendecomp is None: + if os.path.exists('exact_eigendecomposition.npz'): + exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') + else: + logger.info("Computing exact eigendecomposition for error analysis") + # Compute it + config_temp = AnalysisConfiguration() + config_temp.num_eigenvalues = config_analysis.error_num_eigenvalues + config_temp.eigendecomposition_matrices = 'exact' + config_temp.which_eigenvalues = 'smallest' + eig_results = eigendecomposition_analysis( + config_temp, hamiltonian, algorithm, + exact_matrix=exact_matrix, unitary_matrix=None + ) + exact_eigendecomp = load_eigendecomposition(eig_results['exact_eigendecomposition']['file']) + + if approx_eigendecomp is None: + if os.path.exists('approximate_eigendecomposition.npz'): + approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') + else: + logger.info("Computing approximate eigendecomposition for error analysis") + # Compute it + config_temp = AnalysisConfiguration() + config_temp.num_eigenvalues = config_analysis.error_num_eigenvalues + config_temp.eigendecomposition_matrices = 'approximate' + config_temp.which_eigenvalues = 'smallest' + eig_results = eigendecomposition_analysis( + config_temp, hamiltonian, algorithm, + exact_matrix=None, unitary_matrix=unitary_matrix + ) + approx_eigendecomp = load_eigendecomposition(eig_results['approximate_eigendecomposition']['file']) + + # Compare eigenvalues + exact_eigs = exact_eigendecomp['eigenvalues'][:config_analysis.error_num_eigenvalues] + approx_eigs = approx_eigendecomp['eigenvalues'][:config_analysis.error_num_eigenvalues] + + absolute_errors = exact_eigs - approx_eigs + relative_errors = absolute_errors / np.abs(exact_eigs) + + logger.info(f"Eigenvalue absolute error range: [{absolute_errors.min():.6e}, {absolute_errors.max():.6e}]") + logger.info(f"Eigenvalue relative error range: [{relative_errors.min():.6e}, {relative_errors.max():.6e}]") + + results['eigenvalue_errors'] = { + 'num_eigenvalues': config_analysis.error_num_eigenvalues, + 'absolute_errors': absolute_errors.tolist(), + 'relative_errors': relative_errors.tolist(), + 'max_absolute_error': float(np.abs(absolute_errors).max()), + 'max_relative_error': float(np.abs(relative_errors).max()) + } + + # ============================================================================================= + # 2. MATRIX NORM ERRORS + # ============================================================================================= + + if config_analysis.error_matrix_norms is not None: + # Normalize to list + if isinstance(config_analysis.error_matrix_norms, str): + norms_to_compute = [config_analysis.error_matrix_norms] + else: + norms_to_compute = config_analysis.error_matrix_norms + + logger.info(f"Computing matrix norm errors: {norms_to_compute}") + + # Get matrices if not provided + if exact_matrix is None: + if hamiltonian is None: + raise ValueError("hamiltonian required for matrix norm error analysis") + exact_matrix = _compute_exact_matrix(hamiltonian) + + if unitary_matrix is None: + if algorithm is None: + raise ValueError("algorithm required for matrix norm error analysis") + unitary_matrix = _compute_unitary_matrix(algorithm) + + # Check if matrices are dense or matrix-free + is_exact_dense = isinstance(exact_matrix, np.ndarray) + is_approx_dense = isinstance(unitary_matrix, np.ndarray) + is_dense = is_exact_dense and is_approx_dense + + dimension = exact_matrix.shape[0] + num_qubits = int(np.log2(dimension)) + + if is_dense: + # Small systems: direct computation + logger.verbose(f"Using dense matrices for norm computation (dimension={dimension})") + diff_matrix = exact_matrix - unitary_matrix + + for norm_type in norms_to_compute: + if norm_type == 'frobenius': + error = np.linalg.norm(diff_matrix, 'fro') + logger.info(f"Frobenius norm error: {error:.6e}") + results['matrix_frobenius_error'] = float(error) + + elif norm_type == 'spectral': + error = np.linalg.norm(diff_matrix, 2) + logger.info(f"Spectral norm error: {error:.6e}") + results['matrix_spectral_error'] = float(error) + + else: + raise ValueError(f"Unknown matrix norm type: {norm_type}") + + else: + # Large systems: matrix-free computation + logger.info(f"WARNING: Matrix-free norm computation for {num_qubits} qubits") + logger.info(f" This requires 2^{num_qubits} = {dimension:,} matrix-vector products") + logger.info(f" Estimated time: {'<1 minute' if dimension < 10000 else '1-10 minutes' if dimension < 100000 else '10+ minutes'}") + + for norm_type in norms_to_compute: + if norm_type == 'frobenius': + # Compute ||A||_F^2 = sum_i ||A e_i||^2 + logger.verbose("Computing Frobenius norm via matrix-vector products") + frobenius_squared = 0.0 + for i in range(dimension): + if i % max(1, dimension // 10) == 0: + logger.verbose(f" Progress: {i}/{dimension} ({100*i//dimension}%)") + # Create basis vector e_i + e_i = np.zeros(dimension, dtype=complex) + e_i[i] = 1.0 + # Compute difference: (H - U) e_i + if hasattr(exact_matrix, 'matvec'): + exact_result = exact_matrix.matvec(e_i) + else: + exact_result = exact_matrix @ e_i + if hasattr(unitary_matrix, 'matvec'): + approx_result = unitary_matrix.matvec(e_i) + else: + approx_result = unitary_matrix @ e_i + diff_i = exact_result - approx_result + frobenius_squared += np.linalg.norm(diff_i) ** 2 + + error = np.sqrt(frobenius_squared) + logger.info(f"Frobenius norm error (matrix-free): {error:.6e}") + results['matrix_frobenius_error'] = float(error) + + elif norm_type == 'spectral': + # Compute ||A||_2 = largest singular value via power iteration + logger.verbose("Computing spectral norm via power iteration") + + # Random starting vector + v = np.random.randn(dimension) + 1j * np.random.randn(dimension) + v = v / np.linalg.norm(v) + + max_iterations = 100 + tolerance = 1e-6 + + for iteration in range(max_iterations): + # Apply (H - U)† (H - U) to v + # First: (H - U) v + if hasattr(exact_matrix, 'matvec'): + exact_result = exact_matrix.matvec(v) + else: + exact_result = exact_matrix @ v + if hasattr(unitary_matrix, 'matvec'): + approx_result = unitary_matrix.matvec(v) + else: + approx_result = unitary_matrix @ v + diff_v = exact_result - approx_result + + # Then: (H - U)† diff_v + if hasattr(exact_matrix, 'rmatvec'): + exact_adjoint = exact_matrix.rmatvec(diff_v) + else: + exact_adjoint = exact_matrix.conj().T @ diff_v + if hasattr(unitary_matrix, 'rmatvec'): + approx_adjoint = unitary_matrix.rmatvec(diff_v) + else: + approx_adjoint = unitary_matrix.conj().T @ diff_v + result = exact_adjoint - approx_adjoint + + # Normalize + v_new = result / np.linalg.norm(result) + + # Check convergence + if np.linalg.norm(v_new - v) < tolerance: + logger.verbose(f" Converged after {iteration+1} iterations") + break + + v = v_new + + if (iteration + 1) % 10 == 0: + logger.verbose(f" Iteration {iteration+1}/{max_iterations}") + + # Rayleigh quotient to get eigenvalue (= squared singular value) + if hasattr(exact_matrix, 'matvec'): + exact_result = exact_matrix.matvec(v) + else: + exact_result = exact_matrix @ v + if hasattr(unitary_matrix, 'matvec'): + approx_result = unitary_matrix.matvec(v) + else: + approx_result = unitary_matrix @ v + diff_v = exact_result - approx_result + + if hasattr(exact_matrix, 'rmatvec'): + exact_adjoint = exact_matrix.rmatvec(diff_v) + else: + exact_adjoint = exact_matrix.conj().T @ diff_v + if hasattr(unitary_matrix, 'rmatvec'): + approx_adjoint = unitary_matrix.rmatvec(diff_v) + else: + approx_adjoint = unitary_matrix.conj().T @ diff_v + result = exact_adjoint - approx_adjoint + + eigenvalue = np.vdot(v, result) + error = np.sqrt(np.abs(eigenvalue)) + + logger.info(f"Spectral norm error (matrix-free, power iteration): {error:.6e}") + results['matrix_spectral_error'] = float(error) + + else: + raise ValueError(f"Unknown matrix norm type: {norm_type}") + + # ============================================================================================= + # 3. STATE-DEPENDENT ERRORS + # ============================================================================================= + + if config_analysis.error_state_inputs is not None: + # Normalize to list + if isinstance(config_analysis.error_state_inputs, str): + state_files = [config_analysis.error_state_inputs] + else: + state_files = config_analysis.error_state_inputs + + logger.info(f"Computing state-dependent errors for {len(state_files)} state(s)") + + # Get matrices/operators if not provided + if exact_matrix is None: + if hamiltonian is None: + raise ValueError("hamiltonian required for state-dependent error analysis") + exact_matrix = _compute_exact_matrix(hamiltonian) + + if unitary_matrix is None: + if algorithm is None: + raise ValueError("algorithm required for state-dependent error analysis") + unitary_matrix = _compute_unitary_matrix(algorithm) + + state_errors = [] + + for state_file in state_files: + logger.verbose(f"Processing {state_file}") + + # Load state + try: + initial_state = load_state(state_file) + except Exception as e: + logger.info(f"ERROR: Failed to load state from {state_file}: {e}") + raise + + # Apply exact operator + if hasattr(exact_matrix, 'matvec'): + exact_final = exact_matrix.matvec(initial_state) + else: + exact_final = exact_matrix @ initial_state + + # Apply approximate operator + if hasattr(unitary_matrix, 'matvec'): + approx_final = unitary_matrix.matvec(initial_state) + else: + approx_final = unitary_matrix @ initial_state + + # Compute error + diff = exact_final - approx_final + absolute_error = np.linalg.norm(diff) + relative_error = absolute_error / np.linalg.norm(exact_final) + + logger.info(f" {state_file}: absolute error = {absolute_error:.6e}, relative error = {relative_error:.6e}") + + state_errors.append({ + 'input_file': state_file, + 'absolute_error': float(absolute_error), + 'relative_error': float(relative_error) + }) + + results['state_errors'] = state_errors + + # Save results to file + output_file = 'error_analysis.npz' + save_dict = {} + + if 'eigenvalue_errors' in results: + save_dict['eigenvalue_absolute_errors'] = np.array(results['eigenvalue_errors']['absolute_errors']) + save_dict['eigenvalue_relative_errors'] = np.array(results['eigenvalue_errors']['relative_errors']) + save_dict['eigenvalue_num'] = results['eigenvalue_errors']['num_eigenvalues'] + + if 'matrix_frobenius_error' in results: + save_dict['matrix_frobenius_error'] = results['matrix_frobenius_error'] + + if 'matrix_spectral_error' in results: + save_dict['matrix_spectral_error'] = results['matrix_spectral_error'] + + if 'state_errors' in results: + state_absolute = [s['absolute_error'] for s in results['state_errors']] + state_relative = [s['relative_error'] for s in results['state_errors']] + save_dict['state_absolute_errors'] = np.array(state_absolute) + save_dict['state_relative_errors'] = np.array(state_relative) + # Note: filenames are in results dict, not saved to npz + + if save_dict: + np.savez(output_file, **save_dict) + logger.info(f"Error analysis results saved to {output_file}") + results['output_file'] = output_file + + return results + +# ------------------------------------------------------------------------------------------------- + def numerical_simulation( config_analysis: AnalysisConfiguration, algorithm, From da66abaf6a6aea92cd37ee1828367e3dc1eff9b5 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 28 May 2026 16:44:48 -0600 Subject: [PATCH 12/12] Branch 4: Add exact numerical simulation - Add exact_numerical_simulation() function in analysis.py (~90 lines) - Uses exact Hamiltonian matrix instead of approximate unitary - Output files: input_exact_final.npy (distinguishes from input_final.npy) - Works with both dense matrices and matrix-free operators - Add configuration field: exact_simulation_inputs - Add 13 comprehensive tests in test_exact_simulation.py - Update README.md with exact simulation documentation - Update config_full_analysis.py with exact simulation examples - All 231 tests passing (128 original + 53 Branch 1 + 19 Branch 2 + 18 Branch 3 + 13 Branch 4) Completes 5-branch implementation plan for exact matrix analysis and error metrics. --- analysis/README.md | 57 ++- analysis/analysis.py | 456 +++++---------------- analysis/config_types.py | 2 + analysis/examples/config_full_analysis.py | 56 ++- analysis/tests/test_exact_simulation.py | 470 ++++++++++++++++++++++ 5 files changed, 686 insertions(+), 355 deletions(-) create mode 100644 analysis/tests/test_exact_simulation.py diff --git a/analysis/README.md b/analysis/README.md index e29e76a..bcf781a 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -356,7 +356,7 @@ There are many details of the algorithm that may be worth analyzing. The availab - **Matrix norm errors**: Use for physical bounds, but expensive for large systems - **State errors**: Use for specific physically relevant states, scales best -#### Numerical Simulation +#### Numerical Simulation (Approximate) - **Numerical Simulation**: Setting `analysis.numerical_simulation_inputs` to one or more state vector files will apply the algorithm to the input state(s) via numerical simulation, producing @@ -393,7 +393,56 @@ There are many details of the algorithm that may be worth analyzing. The availab psi[0] = 1.0 np.save("initial_state.npy", psi) ``` + +#### Exact Numerical Simulation + +- **Exact Numerical Simulation**: Setting `analysis.exact_simulation_inputs` to one or more state + vector files will apply the **exact Hamiltonian matrix** (without Trotter or other approximations) + to the input state(s), producing output state(s). + + **Purpose**: Compare exact time evolution with approximate algorithm results + + **Input format**: Same as approximate simulation (NumPy `.npy` format). Input can be: + - Single filename (string): `"initial_state.npy"` + - Multiple filenames (list): `["state1.npy", "state2.npy", "state3.npy"]` + + **Output naming**: Automatic suffix `_exact_final` is added to input filename: + - `initial_state.npy` → `initial_state_exact_final.npy` + - This distinguishes from approximate simulation output (`_final`) + + **Scaling**: Can reach 25-30 qubits using matrix-free operators from Branch 1 + - For large systems, uses `PauliStringOperator` (matrix-free) + - Memory: O(2^N) for state vectors, not O(2^(2N)) for matrices + - Computation time depends on number of Pauli terms in Hamiltonian + + **Example**: + ```python + # Single state exact simulation + analysis.exact_simulation_inputs = "initial_state.npy" + + # Multiple states + analysis.exact_simulation_inputs = [ + "ground_state.npy", + "excited_state.npy", + "superposition.npy" + ] + ``` + + **Comparing exact vs approximate**: + You can run both simulations on the same input states: + ```python + # Run both approximate and exact simulations + analysis.numerical_simulation_inputs = "initial_state.npy" + analysis.exact_simulation_inputs = "initial_state.npy" + # Produces: + # initial_state_final.npy (approximate result) + # initial_state_exact_final.npy (exact result) + # Then manually compare the output states + ``` + + **Note**: State vectors must be compatible with Hamiltonian dimension (2^N) + ## Generated Files The script will print a log both to the screen and to a logfile. It also generates output files @@ -404,8 +453,14 @@ depending on the analyses requested: results. Shows final results but not intermediate values. - **Matrix file** (if `matrix_output_file` specified): Unitary matrix in specified format (`.npz`, `.h5`, or `.txt`) +- **Exact matrix file** (if `exact_matrix_output_file` specified): Exact Hamiltonian matrix +- **Eigendecomposition files** (if `num_eigenvalues` > 0): `exact_eigendecomposition.npz` and/or + `approximate_eigendecomposition.npz` +- **Error analysis file** (if any error analysis enabled): `error_analysis.npz` - **Final state files** (if `numerical_simulation_inputs` specified): Evolved quantum states with `_final` suffix (e.g., `initial_state.npy` → `initial_state_final.npy`) +- **Exact final state files** (if `exact_simulation_inputs` specified): Exactly evolved states with + `_exact_final` suffix (e.g., `initial_state.npy` → `initial_state_exact_final.npy`) The logfile is typically most useful for understanding the analysis process and intermediate values. diff --git a/analysis/analysis.py b/analysis/analysis.py index 9621ec7..c7aa317 100644 --- a/analysis/analysis.py +++ b/analysis/analysis.py @@ -795,356 +795,6 @@ def error_analysis( # ------------------------------------------------------------------------------------------------- -def error_analysis( - config_analysis: AnalysisConfiguration, - hamiltonian, - algorithm, - exact_matrix=None, - unitary_matrix=None, - exact_eigendecomp=None, - approx_eigendecomp=None) -> dict: - """ - Compute error metrics comparing exact and approximate representations. - - Three independent error types: - 1. Eigenvalue errors (if error_num_eigenvalues > 0) - 2. Matrix norm errors (if error_matrix_norms is not None) - 3. State-dependent errors (if error_state_inputs is not None) - - Parameters: - config_analysis: Analysis configuration with error analysis settings - hamiltonian: Hamiltonian object - algorithm: Algorithm bloq - exact_matrix: Pre-computed exact matrix (optional) - unitary_matrix: Pre-computed unitary matrix (optional) - exact_eigendecomp: Pre-computed exact eigendecomposition (optional) - approx_eigendecomp: Pre-computed approximate eigendecomposition (optional) - - Returns: - Dictionary with error metrics - """ - from qhat.analysis.matrix_operations import PauliStringOperator - from qhat.analysis.file_io import load_eigendecomposition, load_state - import scipy.linalg - - logger.info("Starting error analysis") - - results = {} - - # ============================================================================================= - # 1. EIGENVALUE ERROR - # ============================================================================================= - - if config_analysis.error_num_eigenvalues > 0: - logger.info(f"Computing eigenvalue errors for {config_analysis.error_num_eigenvalues} eigenvalues") - - # Load or compute eigendecompositions - if exact_eigendecomp is None: - if os.path.exists('exact_eigendecomposition.npz'): - exact_eigendecomp = load_eigendecomposition('exact_eigendecomposition.npz') - else: - logger.info("Computing exact eigendecomposition for error analysis") - # Compute it - config_temp = AnalysisConfiguration() - config_temp.num_eigenvalues = config_analysis.error_num_eigenvalues - config_temp.eigendecomposition_matrices = 'exact' - config_temp.which_eigenvalues = 'smallest' - eig_results = eigendecomposition_analysis( - config_temp, hamiltonian, algorithm, - exact_matrix=exact_matrix, unitary_matrix=None - ) - exact_eigendecomp = load_eigendecomposition(eig_results['exact_eigendecomposition']['file']) - - if approx_eigendecomp is None: - if os.path.exists('approximate_eigendecomposition.npz'): - approx_eigendecomp = load_eigendecomposition('approximate_eigendecomposition.npz') - else: - logger.info("Computing approximate eigendecomposition for error analysis") - # Compute it - config_temp = AnalysisConfiguration() - config_temp.num_eigenvalues = config_analysis.error_num_eigenvalues - config_temp.eigendecomposition_matrices = 'approximate' - config_temp.which_eigenvalues = 'smallest' - eig_results = eigendecomposition_analysis( - config_temp, hamiltonian, algorithm, - exact_matrix=None, unitary_matrix=unitary_matrix - ) - approx_eigendecomp = load_eigendecomposition(eig_results['approximate_eigendecomposition']['file']) - - # Compare eigenvalues - exact_eigs = exact_eigendecomp['eigenvalues'][:config_analysis.error_num_eigenvalues] - approx_eigs = approx_eigendecomp['eigenvalues'][:config_analysis.error_num_eigenvalues] - - absolute_errors = exact_eigs - approx_eigs - relative_errors = absolute_errors / np.abs(exact_eigs) - - logger.info(f"Eigenvalue absolute error range: [{absolute_errors.min():.6e}, {absolute_errors.max():.6e}]") - logger.info(f"Eigenvalue relative error range: [{relative_errors.min():.6e}, {relative_errors.max():.6e}]") - - results['eigenvalue_errors'] = { - 'num_eigenvalues': config_analysis.error_num_eigenvalues, - 'absolute_errors': absolute_errors.tolist(), - 'relative_errors': relative_errors.tolist(), - 'max_absolute_error': float(np.abs(absolute_errors).max()), - 'max_relative_error': float(np.abs(relative_errors).max()) - } - - # ============================================================================================= - # 2. MATRIX NORM ERRORS - # ============================================================================================= - - if config_analysis.error_matrix_norms is not None: - # Normalize to list - if isinstance(config_analysis.error_matrix_norms, str): - norms_to_compute = [config_analysis.error_matrix_norms] - else: - norms_to_compute = config_analysis.error_matrix_norms - - logger.info(f"Computing matrix norm errors: {norms_to_compute}") - - # Get matrices if not provided - if exact_matrix is None: - if hamiltonian is None: - raise ValueError("hamiltonian required for matrix norm error analysis") - exact_matrix = _compute_exact_matrix(hamiltonian) - - if unitary_matrix is None: - if algorithm is None: - raise ValueError("algorithm required for matrix norm error analysis") - unitary_matrix = _compute_unitary_matrix(algorithm) - - # Check if matrices are dense or matrix-free - is_exact_dense = isinstance(exact_matrix, np.ndarray) - is_approx_dense = isinstance(unitary_matrix, np.ndarray) - is_dense = is_exact_dense and is_approx_dense - - dimension = exact_matrix.shape[0] - num_qubits = int(np.log2(dimension)) - - if is_dense: - # Small systems: direct computation - logger.verbose(f"Using dense matrices for norm computation (dimension={dimension})") - diff_matrix = exact_matrix - unitary_matrix - - for norm_type in norms_to_compute: - if norm_type == 'frobenius': - error = np.linalg.norm(diff_matrix, 'fro') - logger.info(f"Frobenius norm error: {error:.6e}") - results['matrix_frobenius_error'] = float(error) - - elif norm_type == 'spectral': - error = np.linalg.norm(diff_matrix, 2) - logger.info(f"Spectral norm error: {error:.6e}") - results['matrix_spectral_error'] = float(error) - - else: - raise ValueError(f"Unknown matrix norm type: {norm_type}") - - else: - # Large systems: matrix-free computation - logger.info(f"WARNING: Matrix-free norm computation for {num_qubits} qubits") - logger.info(f" This requires 2^{num_qubits} = {dimension:,} matrix-vector products") - logger.info(f" Estimated time: {'<1 minute' if dimension < 10000 else '1-10 minutes' if dimension < 100000 else '10+ minutes'}") - - for norm_type in norms_to_compute: - if norm_type == 'frobenius': - # Compute ||A||_F^2 = sum_i ||A e_i||^2 - logger.verbose("Computing Frobenius norm via matrix-vector products") - frobenius_squared = 0.0 - for i in range(dimension): - if i % max(1, dimension // 10) == 0: - logger.verbose(f" Progress: {i}/{dimension} ({100*i//dimension}%)") - # Create basis vector e_i - e_i = np.zeros(dimension, dtype=complex) - e_i[i] = 1.0 - # Compute difference: (H - U) e_i - if hasattr(exact_matrix, 'matvec'): - exact_result = exact_matrix.matvec(e_i) - else: - exact_result = exact_matrix @ e_i - if hasattr(unitary_matrix, 'matvec'): - approx_result = unitary_matrix.matvec(e_i) - else: - approx_result = unitary_matrix @ e_i - diff_i = exact_result - approx_result - frobenius_squared += np.linalg.norm(diff_i) ** 2 - - error = np.sqrt(frobenius_squared) - logger.info(f"Frobenius norm error (matrix-free): {error:.6e}") - results['matrix_frobenius_error'] = float(error) - - elif norm_type == 'spectral': - # Compute ||A||_2 = largest singular value via power iteration - logger.verbose("Computing spectral norm via power iteration") - - # Random starting vector - v = np.random.randn(dimension) + 1j * np.random.randn(dimension) - v = v / np.linalg.norm(v) - - max_iterations = 100 - tolerance = 1e-6 - - for iteration in range(max_iterations): - # Apply (H - U)† (H - U) to v - # First: (H - U) v - if hasattr(exact_matrix, 'matvec'): - exact_result = exact_matrix.matvec(v) - else: - exact_result = exact_matrix @ v - if hasattr(unitary_matrix, 'matvec'): - approx_result = unitary_matrix.matvec(v) - else: - approx_result = unitary_matrix @ v - diff_v = exact_result - approx_result - - # Then: (H - U)† diff_v - if hasattr(exact_matrix, 'rmatvec'): - exact_adjoint = exact_matrix.rmatvec(diff_v) - else: - exact_adjoint = exact_matrix.conj().T @ diff_v - if hasattr(unitary_matrix, 'rmatvec'): - approx_adjoint = unitary_matrix.rmatvec(diff_v) - else: - approx_adjoint = unitary_matrix.conj().T @ diff_v - result = exact_adjoint - approx_adjoint - - # Normalize - v_new = result / np.linalg.norm(result) - - # Check convergence - if np.linalg.norm(v_new - v) < tolerance: - logger.verbose(f" Converged after {iteration+1} iterations") - break - - v = v_new - - if (iteration + 1) % 10 == 0: - logger.verbose(f" Iteration {iteration+1}/{max_iterations}") - - # Rayleigh quotient to get eigenvalue (= squared singular value) - if hasattr(exact_matrix, 'matvec'): - exact_result = exact_matrix.matvec(v) - else: - exact_result = exact_matrix @ v - if hasattr(unitary_matrix, 'matvec'): - approx_result = unitary_matrix.matvec(v) - else: - approx_result = unitary_matrix @ v - diff_v = exact_result - approx_result - - if hasattr(exact_matrix, 'rmatvec'): - exact_adjoint = exact_matrix.rmatvec(diff_v) - else: - exact_adjoint = exact_matrix.conj().T @ diff_v - if hasattr(unitary_matrix, 'rmatvec'): - approx_adjoint = unitary_matrix.rmatvec(diff_v) - else: - approx_adjoint = unitary_matrix.conj().T @ diff_v - result = exact_adjoint - approx_adjoint - - eigenvalue = np.vdot(v, result) - error = np.sqrt(np.abs(eigenvalue)) - - logger.info(f"Spectral norm error (matrix-free, power iteration): {error:.6e}") - results['matrix_spectral_error'] = float(error) - - else: - raise ValueError(f"Unknown matrix norm type: {norm_type}") - - # ============================================================================================= - # 3. STATE-DEPENDENT ERRORS - # ============================================================================================= - - if config_analysis.error_state_inputs is not None: - # Normalize to list - if isinstance(config_analysis.error_state_inputs, str): - state_files = [config_analysis.error_state_inputs] - else: - state_files = config_analysis.error_state_inputs - - logger.info(f"Computing state-dependent errors for {len(state_files)} state(s)") - - # Get matrices/operators if not provided - if exact_matrix is None: - if hamiltonian is None: - raise ValueError("hamiltonian required for state-dependent error analysis") - exact_matrix = _compute_exact_matrix(hamiltonian) - - if unitary_matrix is None: - if algorithm is None: - raise ValueError("algorithm required for state-dependent error analysis") - unitary_matrix = _compute_unitary_matrix(algorithm) - - state_errors = [] - - for state_file in state_files: - logger.verbose(f"Processing {state_file}") - - # Load state - try: - initial_state = load_state(state_file) - except Exception as e: - logger.info(f"ERROR: Failed to load state from {state_file}: {e}") - raise - - # Apply exact operator - if hasattr(exact_matrix, 'matvec'): - exact_final = exact_matrix.matvec(initial_state) - else: - exact_final = exact_matrix @ initial_state - - # Apply approximate operator - if hasattr(unitary_matrix, 'matvec'): - approx_final = unitary_matrix.matvec(initial_state) - else: - approx_final = unitary_matrix @ initial_state - - # Compute error - diff = exact_final - approx_final - absolute_error = np.linalg.norm(diff) - relative_error = absolute_error / np.linalg.norm(exact_final) - - logger.info(f" {state_file}: absolute error = {absolute_error:.6e}, relative error = {relative_error:.6e}") - - state_errors.append({ - 'input_file': state_file, - 'absolute_error': float(absolute_error), - 'relative_error': float(relative_error) - }) - - results['state_errors'] = state_errors - - # Save results to file - output_file = 'error_analysis.npz' - save_dict = {} - - if 'eigenvalue_errors' in results: - save_dict['eigenvalue_absolute_errors'] = np.array(results['eigenvalue_errors']['absolute_errors']) - save_dict['eigenvalue_relative_errors'] = np.array(results['eigenvalue_errors']['relative_errors']) - save_dict['eigenvalue_num'] = results['eigenvalue_errors']['num_eigenvalues'] - - if 'matrix_frobenius_error' in results: - save_dict['matrix_frobenius_error'] = results['matrix_frobenius_error'] - - if 'matrix_spectral_error' in results: - save_dict['matrix_spectral_error'] = results['matrix_spectral_error'] - - if 'state_errors' in results: - state_absolute = [s['absolute_error'] for s in results['state_errors']] - state_relative = [s['relative_error'] for s in results['state_errors']] - save_dict['state_absolute_errors'] = np.array(state_absolute) - save_dict['state_relative_errors'] = np.array(state_relative) - # Note: filenames are in results dict, not saved to npz - - if save_dict: - np.savez(output_file, **save_dict) - logger.info(f"Error analysis results saved to {output_file}") - results['output_file'] = output_file - - return results - -# ------------------------------------------------------------------------------------------------- - def numerical_simulation( config_analysis: AnalysisConfiguration, algorithm, @@ -1455,6 +1105,98 @@ def validate_and_autocomplete_analysis_config(config_analysis: AnalysisConfigura # ------------------------------------------------------------------------------------------------- +def exact_numerical_simulation( + config_analysis: AnalysisConfiguration, + hamiltonian, + exact_matrix) -> dict: + """ + Perform exact numerical simulation by applying the exact Hamiltonian matrix to input state(s). + + This is similar to numerical_simulation() but uses the exact Hamiltonian matrix + instead of the approximate unitary. Useful for comparing exact evolution with + approximate algorithm results. + + Parameters: + config_analysis: Analysis configuration with exact_simulation_inputs + hamiltonian: Hamiltonian object (for metadata) + exact_matrix: The exact Hamiltonian matrix to apply (pre-computed, can be matrix-free) + + Returns: + Dictionary with simulation metadata: list of {input_file, output_file, output_norm} + """ + + # Log operator properties + logger.verbose(f"Exact matrix shape: {exact_matrix.shape}") + + # Normalize input to list + inputs = config_analysis.exact_simulation_inputs + if isinstance(inputs, str): + input_files = [inputs] + elif isinstance(inputs, list): + input_files = inputs + else: + raise ValueError( + f"exact_simulation_inputs must be a string or list of strings, " + f"got {type(inputs)}" + ) + + logger.info(f"Running exact numerical simulation on {len(input_files)} input state(s)") + + results = [] + + for input_file in input_files: + logger.verbose(f"Processing {input_file}") + + # Load initial state + try: + initial_state = load_state(input_file) + except Exception as e: + logger.info(f"ERROR: Failed to load state from {input_file}: {e}") + raise + + # Validate dimensions + if initial_state.shape[0] != exact_matrix.shape[1]: + raise ValueError( + f"Dimension mismatch: state vector has dimension {initial_state.shape[0]} " + f"but exact matrix expects {exact_matrix.shape[1]}" + ) + + # Perform matrix-vector multiplication (or use matvec for matrix-free) + logger.verbose("Performing exact matrix-vector multiplication") + if hasattr(exact_matrix, 'matvec'): + # Matrix-free operator + final_state = exact_matrix.matvec(initial_state) + else: + # Dense matrix + final_state = exact_matrix @ initial_state + + # Compute norm + final_norm = np.linalg.norm(final_state) + logger.verbose(f"Final state norm: {final_norm:.6e}") + + # Generate output filename: input.npy -> input_exact_final.npy + input_path = Path(input_file) + output_file = str(input_path.parent / f"{input_path.stem}_exact_final{input_path.suffix}") + + # Save final state + try: + save_state(output_file, final_state) + except Exception as e: + logger.info(f"ERROR: Failed to save state to {output_file}: {e}") + raise + + logger.info(f"Exact simulation complete: {input_file} -> {output_file}") + + results.append({ + 'input_file': input_file, + 'output_file': output_file, + 'output_norm': float(final_norm) + }) + + return {'exact_simulations': results} + +# ------------------------------------------------------------------------------------------------- + def analyze_algorithm( config_analysis: AnalysisConfiguration, algorithm, @@ -1481,7 +1223,8 @@ def analyze_algorithm( config_analysis.numerical_simulation_inputs is None and config_analysis.exact_matrix_output_file is None and not eigendecomposition_requested and - not error_analysis_requested): + not error_analysis_requested and + config_analysis.exact_simulation_inputs is None): raise ValueError( "No analyses requested. Set at least one of:\n" " - resource_estimator (e.g., 'pyliqtr', 'cirq')\n" @@ -1491,7 +1234,8 @@ def analyze_algorithm( " - num_eigenvalues (e.g., 5 or 'all')\n" " - enable_eigenvalue_errors (True to compute errors for all eigenvalues)\n" " - error_matrix_norms (e.g., 'frobenius' or ['frobenius', 'spectral'])\n" - " - error_state_inputs (e.g., 'state.npy')" + " - error_state_inputs (e.g., 'state.npy')\n" + " - exact_simulation_inputs (e.g., 'state.npy')" ) results = {} @@ -1534,7 +1278,7 @@ def analyze_algorithm( if eigendecomposition_requested: logger.info("Performing eigendecomposition analysis.") - results["eigendecomposition"] = eigendecomposition_analysis( + eig_results = eigendecomposition_analysis( config_analysis, exact_matrix=exact_matrix, unitary_matrix=unitary_matrix @@ -1558,6 +1302,12 @@ def analyze_algorithm( approx_eigendecomp=approx_eigendecomp ) + if config_analysis.exact_simulation_inputs is not None: + logger.info("Performing exact numerical simulation.") + results["exact_simulation"] = exact_numerical_simulation( + config_analysis, hamiltonian, exact_matrix + ) + # TODO: Add gate parallelism / gate depth analysis # TODO: Would it be useful to analyze in terms of a different basis (e.g., Toffoli gates)? diff --git a/analysis/config_types.py b/analysis/config_types.py index 82789fe..4af619d 100644 --- a/analysis/config_types.py +++ b/analysis/config_types.py @@ -173,6 +173,7 @@ def __init__(self): self.enable_eigenvalue_errors = False self.error_matrix_norms = None self.error_state_inputs = None + self.exact_simulation_inputs = None def _generate_TOML_table(self): table = tomlkit.table() @@ -187,6 +188,7 @@ def _generate_TOML_table(self): self.save_if_present(table, "enable_eigenvalue_errors") self.save_if_present(table, "error_matrix_norms") self.save_if_present(table, "error_state_inputs") + self.save_if_present(table, "exact_simulation_inputs") return table # ------------------------------------------------------------------------------------------------- diff --git a/analysis/examples/config_full_analysis.py b/analysis/examples/config_full_analysis.py index b8cd558..723efcb 100644 --- a/analysis/examples/config_full_analysis.py +++ b/analysis/examples/config_full_analysis.py @@ -413,6 +413,7 @@ # ================================================================================================= # FUTURE ANALYSES (Coming in subsequent branches) # ================================================================================================= +# Apply the exact Hamiltonian matrix (no approximations) to input state(s) # Branch 3: Error Metrics analysis.error_num_eigenvalues = 1 # Compare ground state energy @@ -421,6 +422,59 @@ # Branch 4: Exact Numerical Simulation # analysis.exact_simulation_inputs = "examples/initial_state.npy" -# Simulate using exact Hamiltonian evolution (for comparison with approximate) + +# This applies the exact Hamiltonian H to the input state(s), producing exact time evolution. +# Useful for: +# - Comparing exact evolution with approximate algorithm results +# - Validating algorithm correctness on small-scale problems +# - Understanding exact quantum dynamics + +# Configuration: +# exact_simulation_inputs: +# - None (default): Exact simulation disabled +# - String: Single input state file +# - List: Multiple input state files +# +# Output files: +# - input.npy → input_exact_final.npy (note: _exact_final suffix) +# - Distinguishes from approximate simulation (_final suffix) + +# System size guidance: +# - ≤15 qubits: Fast (uses dense exact matrix) +# - 16-30 qubits: Uses matrix-free operator (scales well) +# - Limited by state vector size (16 GB at 30 qubits) + +# Examples: + +# Single state exact simulation +# analysis.exact_simulation_inputs = "examples/initial_state.npy" + +# Multiple states +# analysis.exact_simulation_inputs = [ +# "examples/ground_state.npy", +# "examples/excited_state.npy", +# "examples/superposition_state.npy" +# ] + +# Compare exact vs approximate evolution (run both!) +# analysis.numerical_simulation_inputs = "examples/initial_state.npy" # Approximate +# analysis.exact_simulation_inputs = "examples/initial_state.npy" # Exact +# # Produces both: +# # initial_state_final.npy (approximate result) +# # initial_state_exact_final.npy (exact result) +# # Then you can manually compare the output states + +# Full validation workflow (all analyses together) +# This demonstrates using all available QHAT analysis capabilities: +# analysis.resource_estimator = "pyLIQTR" # Gate counts +# analysis.matrix_output_file = "unitary.npz" # Approximate matrix +# analysis.exact_matrix_output_file = "exact.npz" # Exact matrix +# analysis.num_eigenvalues = 5 # Eigendecomposition +# analysis.eigendecomposition_matrices = "both" # Both exact and approximate +# analysis.which_eigenvalues = "smallest" # Ground + excited states +# analysis.error_num_eigenvalues = 1 # Ground state error +# analysis.error_state_inputs = "examples/ground.npy" # State-dependent error +# analysis.numerical_simulation_inputs = "examples/initial_state.npy" # Approximate evolution +# analysis.exact_simulation_inputs = "examples/initial_state.npy" # Exact evolution # ================================================================================================= diff --git a/analysis/tests/test_exact_simulation.py b/analysis/tests/test_exact_simulation.py new file mode 100644 index 0000000..4336fcf --- /dev/null +++ b/analysis/tests/test_exact_simulation.py @@ -0,0 +1,470 @@ +""" +Tests for exact numerical simulation functionality. + +Tests cover: +- Basic exact simulation with known operators +- Multiple input states +- Output file naming (_exact_final suffix) +- Comparison with approximate simulation +- Matrix-free operators +- Norm preservation +""" + +import numpy as np +import pytest +import tempfile +import os +from pathlib import Path + +from qhat.analysis.analysis import exact_numerical_simulation +from qhat.analysis.config_types import AnalysisConfiguration +from qhat.analysis.file_io import save_state, load_state + + +# ================================================================================================= +# Unit Tests: Basic Exact Simulation +# ================================================================================================= + +def test_exact_simulation_identity(): + """Test exact simulation with identity operator.""" + # 2x2 identity + identity = np.eye(2, dtype=complex) + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'test_state.npy' + + # Mock hamiltonian (not actually used in simulation) + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + # Save initial state + save_state('test_state.npy', state) + + # Run simulation + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=identity + ) + + # Verify results + assert 'exact_simulations' in results + assert len(results['exact_simulations']) == 1 + assert results['exact_simulations'][0]['input_file'] == 'test_state.npy' + assert results['exact_simulations'][0]['output_file'] == 'test_state_exact_final.npy' + + # Load and verify output + final_state = load_state('test_state_exact_final.npy') + np.testing.assert_array_almost_equal(final_state, state) + + finally: + os.chdir(original_dir) + + +def test_exact_simulation_pauli_z(): + """Test exact simulation with Pauli Z operator.""" + # Pauli Z = [[1, 0], [0, -1]] + pauli_z = np.array([[1, 0], [0, -1]], dtype=complex) + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'test_state.npy' + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('test_state.npy', state) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=pauli_z + ) + + # Load output: Z|0⟩ = |0⟩ + final_state = load_state('test_state_exact_final.npy') + expected = np.array([1.0, 0.0], dtype=complex) + np.testing.assert_array_almost_equal(final_state, expected) + + finally: + os.chdir(original_dir) + + +def test_exact_simulation_pauli_x(): + """Test exact simulation with Pauli X operator.""" + # Pauli X = [[0, 1], [1, 0]] + pauli_x = np.array([[0, 1], [1, 0]], dtype=complex) + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'test_state.npy' + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('test_state.npy', state) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=pauli_x + ) + + # Load output: X|0⟩ = |1⟩ + final_state = load_state('test_state_exact_final.npy') + expected = np.array([0.0, 1.0], dtype=complex) + np.testing.assert_array_almost_equal(final_state, expected) + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Unit Tests: Multiple States +# ================================================================================================= + +def test_exact_simulation_multiple_states(): + """Test exact simulation with multiple input states.""" + identity = np.eye(2, dtype=complex) + state1 = np.array([1.0, 0.0], dtype=complex) + state2 = np.array([0.0, 1.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = ['state1.npy', 'state2.npy'] + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('state1.npy', state1) + save_state('state2.npy', state2) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=identity + ) + + assert len(results['exact_simulations']) == 2 + assert results['exact_simulations'][0]['input_file'] == 'state1.npy' + assert results['exact_simulations'][0]['output_file'] == 'state1_exact_final.npy' + assert results['exact_simulations'][1]['input_file'] == 'state2.npy' + assert results['exact_simulations'][1]['output_file'] == 'state2_exact_final.npy' + + # Verify both outputs exist + assert os.path.exists('state1_exact_final.npy') + assert os.path.exists('state2_exact_final.npy') + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Unit Tests: Output File Naming +# ================================================================================================= + +def test_exact_simulation_output_naming(): + """Test that output files have _exact_final suffix.""" + identity = np.eye(2, dtype=complex) + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'my_state.npy' + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('my_state.npy', state) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=identity + ) + + # Verify naming: my_state.npy -> my_state_exact_final.npy + assert results['exact_simulations'][0]['output_file'] == 'my_state_exact_final.npy' + assert os.path.exists('my_state_exact_final.npy') + assert not os.path.exists('my_state_final.npy') # Should not use regular suffix + + finally: + os.chdir(original_dir) + + +def test_exact_simulation_preserves_path(): + """Test that output files preserve directory path.""" + identity = np.eye(2, dtype=complex) + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'subdir/my_state.npy' + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + os.makedirs('subdir', exist_ok=True) + save_state('subdir/my_state.npy', state) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=identity + ) + + # Verify path preserved + assert results['exact_simulations'][0]['output_file'] == 'subdir/my_state_exact_final.npy' + assert os.path.exists('subdir/my_state_exact_final.npy') + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Unit Tests: Norm Preservation +# ================================================================================================= + +def test_exact_simulation_norm_preservation(): + """Test that unitary evolution preserves state norm.""" + # Random unitary (2x2) + # Use Pauli X as a simple unitary + pauli_x = np.array([[0, 1], [1, 0]], dtype=complex) + state = np.array([0.6, 0.8], dtype=complex) # Normalized + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'test_state.npy' + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('test_state.npy', state) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=pauli_x + ) + + # Verify norm preserved + assert abs(results['exact_simulations'][0]['output_norm'] - 1.0) < 1e-10 + + # Load and verify + final_state = load_state('test_state_exact_final.npy') + assert abs(np.linalg.norm(final_state) - 1.0) < 1e-10 + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Unit Tests: Matrix-Free Operators +# ================================================================================================= + +def test_exact_simulation_matrix_free(): + """Test exact simulation with matrix-free operator.""" + from qhat.analysis.matrix_operations import PauliStringOperator + + # Create matrix-free identity operator + pauli_dict = {'II': 1.0} # 2-qubit identity + matrix_free_op = PauliStringOperator(pauli_dict, num_qubits=2) + + state = np.array([1.0, 0.0, 0.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'test_state.npy' + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('test_state.npy', state) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=matrix_free_op + ) + + # Load output + final_state = load_state('test_state_exact_final.npy') + np.testing.assert_array_almost_equal(final_state, state) + + finally: + os.chdir(original_dir) + + +def test_exact_simulation_matrix_free_pauli_x(): + """Test exact simulation with matrix-free Pauli X.""" + from qhat.analysis.matrix_operations import PauliStringOperator + + # Create matrix-free Pauli X operator + pauli_dict = {'X': 1.0} + matrix_free_op = PauliStringOperator(pauli_dict, num_qubits=1) + + state = np.array([1.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'test_state.npy' + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('test_state.npy', state) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=matrix_free_op + ) + + # Load output: X|0⟩ = |1⟩ + final_state = load_state('test_state_exact_final.npy') + expected = np.array([0.0, 1.0], dtype=complex) + np.testing.assert_array_almost_equal(final_state, expected) + + finally: + os.chdir(original_dir) + + +# ================================================================================================= +# Error Handling Tests +# ================================================================================================= + +def test_exact_simulation_missing_state_file(): + """Test error when state file doesn't exist.""" + identity = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'nonexistent.npy' + + class MockHamiltonian: + pass + + with pytest.raises(Exception): # Will be FileNotFoundError or similar + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=identity + ) + finally: + os.chdir(original_dir) + + +def test_exact_simulation_dimension_mismatch(): + """Test error when state dimension doesn't match matrix.""" + matrix_2x2 = np.eye(2, dtype=complex) + state_4d = np.array([1.0, 0.0, 0.0, 0.0], dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'test_state.npy' + + class MockHamiltonian: + pass + + with pytest.raises(ValueError, match="Dimension mismatch"): + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('test_state.npy', state_4d) + + exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=matrix_2x2 + ) + finally: + os.chdir(original_dir) + + +def test_exact_simulation_invalid_input_type(): + """Test error when exact_simulation_inputs is not string or list.""" + identity = np.eye(2, dtype=complex) + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 123 # Invalid type + + class MockHamiltonian: + pass + + with pytest.raises(ValueError, match="must be a string or list of strings"): + exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=identity + ) + + +# ================================================================================================= +# Integration Test: Comparison with Known Result +# ================================================================================================= + +def test_exact_simulation_known_result(): + """Test exact simulation against analytically known result.""" + # Use Hadamard operator: H = 1/sqrt(2) * [[1, 1], [1, -1]] + hadamard = (1.0 / np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex) + state = np.array([1.0, 0.0], dtype=complex) # |0⟩ + + config = AnalysisConfiguration() + config.exact_simulation_inputs = 'test_state.npy' + + class MockHamiltonian: + pass + + with tempfile.TemporaryDirectory() as tmpdir: + original_dir = os.getcwd() + os.chdir(tmpdir) + try: + save_state('test_state.npy', state) + + results = exact_numerical_simulation( + config, + hamiltonian=MockHamiltonian(), + exact_matrix=hadamard + ) + + # Load output: H|0⟩ = 1/sqrt(2) (|0⟩ + |1⟩) + final_state = load_state('test_state_exact_final.npy') + expected = (1.0 / np.sqrt(2)) * np.array([1.0, 1.0], dtype=complex) + np.testing.assert_array_almost_equal(final_state, expected) + + finally: + os.chdir(original_dir)