Utilities and I/O

def qsave(obj, name: str, format: str = 'pickle') -> None:
    """
    Save QuTiP object to file.
    
    Parameters:
    - obj: QuTiP object to save (Qobj, Result, etc.)
    - name: Filename (without extension)
    - format: File format ('pickle', 'numpy')
    """

def qload(name: str, format: str = 'pickle') -> object:
    """
    Load QuTiP object from file.
    
    Parameters:
    - name: Filename (without extension)
    - format: File format ('pickle', 'numpy')
    
    Returns:
    - object: Loaded QuTiP object
    """

def file_data_read(filename: str) -> np.ndarray:
    """
    Read numerical data from file.
    
    Parameters:
    - filename: Path to data file
    
    Returns:
    - ndarray: Loaded data array
    """

def file_data_store(filename: str, data: np.ndarray, numtype: str = "complex", 
                   numformat: str = "decimal") -> None:
    """
    Store numerical data to file.
    
    Parameters:
    - filename: Output filename
    - data: Data array to store
    - numtype: 'real' or 'complex'
    - numformat: 'decimal' or 'exp'
    """

def about() -> None:
    """Display information about QuTiP installation and dependencies."""

def cite() -> str:
    """
    Return citation information for QuTiP in BibTeX format.
    
    Returns:
    - str: BibTeX citation string
    """

def n_thermal(w: float, w_th: float) -> float:
    """
    Calculate average thermal occupation number.
    
    Parameters:
    - w: Frequency
    - w_th: Thermal frequency scale (kB*T/ħ)
    
    Returns:
    - float: Thermal occupation ⟨n⟩ = 1/(exp(ω/ωₜₕ) - 1)
    """

def clebsch(j1: float, j2: float, j3: float, m1: float, m2: float, m3: float) -> float:
    """
    Calculate Clebsch-Gordan coefficient ⟨j₁m₁j₂m₂|j₃m₃⟩.
    
    Parameters:
    - j1, j2, j3: Angular momentum quantum numbers
    - m1, m2, m3: Magnetic quantum numbers
    
    Returns:
    - float: Clebsch-Gordan coefficient
    """

def convert_unit(value: float, unit_from: str, unit_to: str) -> float:
    """
    Convert between different physical units.
    
    Parameters:
    - value: Numerical value to convert
    - unit_from: Source unit
    - unit_to: Target unit
    
    Returns:
    - float: Converted value
    """

def iterated_fit(func, data, init_params: list, maxiter: int = 1000, 
                 tol: float = 1e-8) -> tuple:
    """
    Iterative fitting routine for complex functions.
    
    Parameters:
    - func: Function to fit
    - data: Data to fit to
    - init_params: Initial parameter guess
    - maxiter: Maximum iterations
    - tol: Convergence tolerance
    
    Returns:
    - tuple: (fitted_params, residuals, convergence_info)
    """

def polynomial_fit(x: np.ndarray, y: np.ndarray, degree: int) -> np.ndarray:
    """
    Fit polynomial to data points.
    
    Parameters:
    - x: X-coordinate data
    - y: Y-coordinate data  
    - degree: Polynomial degree
    
    Returns:
    - ndarray: Polynomial coefficients
    """

def assert_equal(obj1, obj2, tol: float = 1e-12) -> None:
    """
    Assert two QuTiP objects are equal within tolerance.
    
    Parameters:
    - obj1: First object
    - obj2: Second object
    - tol: Numerical tolerance
    
    Raises:
    - AssertionError: If objects are not equal
    """

def assert_hermitian(op: Qobj, tol: float = 1e-12) -> None:
    """
    Assert that operator is Hermitian.
    
    Parameters:
    - op: Operator to test
    - tol: Numerical tolerance
    
    Raises:
    - AssertionError: If operator is not Hermitian
    """

def assert_unitary(U: Qobj, tol: float = 1e-12) -> None:
    """
    Assert that operator is unitary.
    
    Parameters:
    - U: Operator to test
    - tol: Numerical tolerance
    
    Raises:
    - AssertionError: If operator is not unitary
    """

def run_tests() -> None:
    """Run QuTiP test suite."""

class Settings:
    """
    Global settings object for QuTiP configuration.
    
    Attributes:
        atol: Absolute tolerance for numerical comparisons
        rtol: Relative tolerance for numerical comparisons
        auto_tidyup: Automatically clean small matrix elements
        auto_tidyup_atol: Absolute tolerance for auto tidyup
        num_cpus: Number of CPUs for parallel processing
        debug: Enable debug mode
    """
    def __init__(self):
        self.core = {
            'atol': 1e-12,
            'rtol': 1e-12,
            'auto_tidyup': True,
            'auto_tidyup_atol': 1e-12
        }
        self.compile = {
            'num_cpus': 0  # 0 = auto-detect
        }

# Global settings instance
settings = Settings()

def profile_solver(solver_func, *args, **kwargs) -> dict:
    """
    Profile solver performance.
    
    Parameters:
    - solver_func: Solver function to profile
    - *args, **kwargs: Arguments for solver
    
    Returns:
    - dict: Profiling results with timing information
    """

def memory_usage() -> dict:
    """
    Get current memory usage information.
    
    Returns:
    - dict: Memory usage statistics
    """

import qutip as qt
import numpy as np
import tempfile
import os

# File I/O operations
# Create some quantum objects
psi = qt.bell_state('00')
rho = psi * psi.dag()
H = qt.tensor(qt.sigmaz(), qt.sigmaz())

# Save objects to files
with tempfile.TemporaryDirectory() as tmpdir:
    # Save quantum objects
    qt.qsave(psi, os.path.join(tmpdir, 'bell_state'))
    qt.qsave(rho, os.path.join(tmpdir, 'bell_dm'))
    qt.qsave(H, os.path.join(tmpdir, 'hamiltonian'))
    
    # Load them back
    psi_loaded = qt.qload(os.path.join(tmpdir, 'bell_state'))
    rho_loaded = qt.qload(os.path.join(tmpdir, 'bell_dm'))
    H_loaded = qt.qload(os.path.join(tmpdir, 'hamiltonian'))
    
    # Verify they're the same
    print(f"State fidelity: {qt.fidelity(psi, psi_loaded):.6f}")
    print(f"DM fidelity: {qt.fidelity(rho, rho_loaded):.6f}")
    print(f"Hamiltonian difference: {(H - H_loaded).norm():.2e}")

# Package information
print("QuTiP Information:")
qt.about()

print("\nCitation:")
citation = qt.cite()
print(citation[:200] + "...")  # Show first 200 characters

# Physical utilities
# Thermal occupation number
frequency = 1.0  # Oscillator frequency
temperature = 0.1  # Temperature in same units
n_avg = qt.n_thermal(frequency, temperature)
print(f"Average thermal occupation at T={temperature}: {n_avg:.3f}")

# Clebsch-Gordan coefficients for spin-1/2 + spin-1/2 = spin-1
j1, j2, j3 = 0.5, 0.5, 1.0
m1, m2, m3 = 0.5, 0.5, 1.0
cg_coeff = qt.clebsch(j1, j2, j3, m1, m2, m3)
print(f"Clebsch-Gordan ⟨1/2,1/2; 1/2,1/2|1,1⟩: {cg_coeff:.3f}")

# All possible CG coefficients for this coupling
print("All CG coefficients for j1=1/2, j2=1/2 -> j3=1:")
for m3_val in [1, 0, -1]:
    for m1_val in [0.5, -0.5]:
        for m2_val in [0.5, -0.5]:
            if abs(m1_val + m2_val - m3_val) < 1e-10:  # Selection rule
                cg = qt.clebsch(0.5, 0.5, 1.0, m1_val, m2_val, m3_val)
                print(f"  ⟨1/2,{m1_val:3.1f}; 1/2,{m2_val:3.1f}|1,{m3_val:2.0f}⟩ = {cg:6.3f}")

# Testing utilities
# Create test objects
op1 = qt.sigmax()
op2 = qt.sigmax()
op_hermitian = qt.sigmaz()
op_unitary = qt.rx(np.pi/4)

# Test equality
try:
    qt.assert_equal(op1, op2)
    print("Operators are equal ✓")
except AssertionError:
    print("Operators are not equal ✗")

# Test if operator is Hermitian
try:
    qt.assert_hermitian(op_hermitian)
    print("Operator is Hermitian ✓")
except AssertionError:
    print("Operator is not Hermitian ✗")

# Test if operator is unitary
try:
    qt.assert_unitary(op_unitary)
    print("Operator is unitary ✓")
except AssertionError:
    print("Operator is not unitary ✗")

# Settings configuration
print(f"Current absolute tolerance: {qt.settings.core['atol']}")
print(f"Auto tidyup enabled: {qt.settings.core['auto_tidyup']}")

# Temporarily change settings
old_atol = qt.settings.core['atol']
qt.settings.core['atol'] = 1e-10
print(f"New absolute tolerance: {qt.settings.core['atol']}")

# Restore settings
qt.settings.core['atol'] = old_atol

# Data file operations
# Generate some data
times = np.linspace(0, 10, 100)
expectation_values = np.exp(-0.1 * times) * np.cos(times)

# Save to file
with tempfile.NamedTemporaryFile(mode='w', delete=False, suffix='.dat') as f:
    data_file = f.name

qt.file_data_store(data_file, np.column_stack([times, expectation_values]), 
                   numtype="real")

# Read back
loaded_data = qt.file_data_read(data_file)
print(f"Original data shape: {np.column_stack([times, expectation_values]).shape}")
print(f"Loaded data shape: {loaded_data.shape}")

# Clean up
os.unlink(data_file)

# Mathematical utilities example
def exponential_decay(t, A, gamma, omega, phi):
    """Damped oscillation function."""
    return A * np.exp(-gamma * t) * np.cos(omega * t + phi)

# Generate noisy data
t_data = np.linspace(0, 5, 50)
true_params = [1.0, 0.2, 2.0, 0.0]  # A, gamma, omega, phi
y_clean = exponential_decay(t_data, *true_params)
y_noisy = y_clean + 0.05 * np.random.randn(len(t_data))

# Fit using scipy.optimize (QuTiP's iterated_fit wrapper)
try:
    from scipy.optimize import curve_fit
    fitted_params, covariance = curve_fit(exponential_decay, t_data, y_noisy, 
                                         p0=[1.1, 0.25, 1.8, 0.1])
    print(f"True parameters: {true_params}")
    print(f"Fitted parameters: {fitted_params}")
    print(f"Parameter errors: {np.sqrt(np.diag(covariance))}")
except ImportError:
    print("SciPy not available for fitting example")

# Performance monitoring
def example_calculation():
    """Simple calculation for profiling."""
    N = 50
    psi = qt.rand_ket(N)
    H = qt.rand_herm(N)
    return qt.expect(H, psi)

# Time the calculation
import time
start_time = time.time()
result = example_calculation()
end_time = time.time()
print(f"Calculation result: {result:.6f}")
print(f"Execution time: {end_time - start_time:.4f} seconds")

class Settings:
    """
    Global configuration object.
    
    Attributes:
        core: Core numerical settings (tolerances, etc.)
        compile: Compilation and performance settings
    """

# Utility functions return various types:
# - File I/O functions return None or loaded objects
# - Physical functions return float values
# - Testing functions return None or raise AssertionError
# - Information functions return None (print to console) or str