Chapter 5: Functions & Modules - Building Reusable Scientific Code

import math

def calculate_mean(values):
    """
    Calculate the arithmetic mean of a list of values.

    This is our first function - notice the structure!
    """
    # Explicit validation: raise an error for invalid input
    if len(values) == 0:
        raise ValueError("Cannot calculate mean of empty list")

    total = math.fsum(values)  # more numerically stable than sum for mixed magnitudes
    count = len(values)
    mean = total / count
    return mean

# Using the function
measurements = [23.5, 24.1, 23.8, 24.3, 23.9]
avg = calculate_mean(measurements)
print(f"Mean temperature: {avg:.2f} degC")

# The validation would raise an error with empty list:
# empty_mean = calculate_mean([])
# ValueError: Cannot calculate mean of empty list

def report_statistics(data, name="Dataset"):
    """Report basic statistics without returning values."""
    if len(data) == 0:
        raise ValueError("report_statistics() requires at least one value")
    mean = sum(data) / len(data)
    minimum = min(data)
    maximum = max(data)

    print(f"Statistics for {name}:")
    print(f"  Mean: {mean:.3f}")
    print(f"  Range: [{minimum:.3f}, {maximum:.3f}]")
    # No return statement - returns None implicitly

# Report some calculations
temperatures = [20.1, 21.5, 19.8, 22.3, 20.9]
report_statistics(temperatures, "Temperature ( degC)")

import math

def analyze_data(values, ddof=0):
    """
    Calculate multiple statistics at once.

    Returns:
        mean, std_dev, min_val, max_val
    """
    n = len(values)
    if n == 0:
        raise ValueError("Cannot analyze empty list")

    mean = math.fsum(values) / n

    # Calculate standard deviation
    # ddof=0 gives population variance (divide by n)
    # ddof=1 gives sample variance (divide by n-1)
    if n - ddof <= 0:
        raise ValueError(f"Need n > ddof for variance (n={n}, ddof={ddof})")
    variance = math.fsum((x - mean)**2 for x in values) / (n - ddof)
    std_dev = variance ** 0.5

    return mean, std_dev, min(values), max(values)

# Analyze experimental data
data = [9.8, 10.1, 9.9, 10.2, 9.7, 10.0]
mean, std, min_val, max_val = analyze_data(data)

print(f"Mean: {mean:.3f}")
print(f"Std Dev: {std:.3f}")
print(f"Range: [{min_val}, {max_val}]")

def kinetic_energy_cgs(mass_g, velocity_cms):
    """
    Calculate kinetic energy in ergs.

    Parameters:
        mass_g: mass in grams
        velocity_cms: velocity in cm/s

    Returns:
        energy in ergs (g*cm^2/s^2)
    """
    import math

    if not (math.isfinite(mass_g) and math.isfinite(velocity_cms)):
        raise ValueError("mass_g and velocity_cms must be finite")
    if mass_g < 0:
        raise ValueError(f"mass_g must be nonnegative, got {mass_g}")
    energy_ergs = 0.5 * mass_g * velocity_cms**2
    return energy_ergs

# Example: electron moving at 1% speed of light
electron_mass = 9.109e-28  # grams
c_light = 2.998e10  # speed of light in cm/s
electron_velocity = 0.01 * c_light  # 1% of c
ke = kinetic_energy_cgs(electron_mass, electron_velocity)

print(f"Electron kinetic energy: {ke:.2e} ergs")
print(f"In eV: {ke/1.602e-12:.2f} eV")  # 1 eV = 1.602e-12 ergs

"""
DESIGN: Function to validate numerical calculation

PURPOSE: Ensure numerical results are physically reasonable
INPUT: value, expected_range
OUTPUT: boolean (True if valid)
PROCESS:
    - Check if value is finite
    - Check if value is within range
"""

def validate_result(value, min_val=None, max_val=None):
    """
    Validate a numerical result is reasonable.

    Simple validation - we'll expand this pattern later!
    """
    import math

    if min_val is not None and max_val is not None and min_val > max_val:
        raise ValueError("min_val must not exceed max_val")

    # Check if finite
    if not math.isfinite(value):
        return False

    # Check bounds if provided
    if min_val is not None and value < min_val:
        return False
    if max_val is not None and value > max_val:
        return False

    return True

# Test validation
energy = kinetic_energy_cgs(1.0, 100.0)
print(f"Energy = {energy} ergs")
print(f"Valid (positive)? {validate_result(energy, min_val=0)}")
print(f"Valid (in range)? {validate_result(energy, 0, 10000)}")

def escape_velocity(mass_g, radius_cm):
    """
    Calculate escape velocity from a spherical body.

    Parameters:
        mass_g: mass of the body in grams
        radius_cm: radius of the body in cm

    Returns:
        escape velocity in cm/s
    """
    if mass_g < 0:
        raise ValueError("mass_g must be nonnegative")
    if radius_cm <= 0:
        raise ValueError("radius_cm must be positive")
    G = 6.674e-8  # gravitational constant in cm^3/(g*s^2)
    v_escape = (2 * G * mass_g / radius_cm) ** 0.5
    return v_escape

# Earth parameters
earth_mass = 5.972e27    # grams
earth_radius = 6.371e8   # cm

# Correct calculation
v_correct = escape_velocity(earth_mass, earth_radius)

# WRONG: Arguments reversed
v_wrong = escape_velocity(earth_radius, earth_mass)

print(f"Correct escape velocity: {v_correct:.2e} cm/s = {v_correct/1e5:.1f} km/s")
print(f"Wrong (args reversed): {v_wrong:.2e} cm/s")
print("Python accepts the reversed call because both inputs are numbers;")
print("the result is physically meaningless because their units were swapped.")

# Keywords make intent clear
v_safe = escape_velocity(radius_cm=earth_radius, mass_g=earth_mass)
print(f"\nWith keywords: {v_safe:.2e} cm/s")

def simulate_decay(initial_atoms, half_life_years=5730, time_years=0):
    """
    Simulate radioactive decay.

    Parameters:
        initial_atoms: starting number of atoms
        half_life_years: half-life in years (default: C-14 = 5730)
        time_years: elapsed time in years (default: 0)
    """
    import math
    if initial_atoms < 0:
        raise ValueError("initial_atoms must be nonnegative")
    if half_life_years <= 0:
        raise ValueError("half_life_years must be positive")
    if time_years < 0:
        raise ValueError("time_years must be nonnegative in this forward-time model")
    remaining = initial_atoms * (0.5 ** (time_years / half_life_years))
    return remaining

# Different usage patterns
n0 = 1000000  # One million atoms

# Use all defaults (returns initial count)
print(f"At t=0: {simulate_decay(n0):.0f} atoms")

# Override just time (uses default half-life for C-14)
print(f"After 5730 years: {simulate_decay(n0, time_years=5730):.0f} atoms")

# Override all parameters (U-238 example)
print(f"U-238 after 1 billion years: {simulate_decay(n0, 4.5e9, 1e9):.0f} atoms")

# THE TRAP - DON'T DO THIS!
def add_measurement_buggy(value, data=[]):
    """BUGGY VERSION - default list is shared!"""
    data.append(value)
    return data

# Watch the disaster unfold
day1 = add_measurement_buggy(23.5)
print(f"Day 1: {day1}")

day2 = add_measurement_buggy(24.1)  # Surprise!
print(f"Day 2: {day2}")  # Contains BOTH days!

print(f"Same list? {day1 is day2}")  # True - it's the same object!

def add_measurement_fixed(value, data=None):
    """CORRECT VERSION - new list each time if needed."""
    if data is None:
        data = []  # Create new list when not provided
    data.append(value)
    return data

# Now it works correctly
day1 = add_measurement_fixed(23.5)
day2 = add_measurement_fixed(24.1)
print(f"Day 1: {day1}")
print(f"Day 2: {day2}")  # Separate lists!

# Can still provide existing list
combined = []
add_measurement_fixed(23.5, combined)
add_measurement_fixed(24.1, combined)
print(f"Combined: {combined}")

def calculate_statistics(*values, method='all'):
    """
    Calculate statistics for any number of values.

    *args collects positional arguments into a tuple.
    """
    if not values:
        return None

    n = len(values)
    mean = sum(values) / n

    if method == 'mean':
        return mean

    # Calculate variance and std dev
    variance = sum((x - mean)**2 for x in values) / n
    std_dev = variance ** 0.5

    if method == 'all':
        return {
            'n': n,
            'mean': mean,
            'std': std_dev,
            'min': min(values),
            'max': max(values)
        }

# Works with any number of arguments!
print(calculate_statistics(1, 2, 3, 4, 5))
print(f"Mean only: {calculate_statistics(2.5, 3.7, 1.2, method='mean'):.2f}")

def run_experiment(name, **parameters):
    """
    Run experiment with flexible parameters.

    **kwargs collects keyword arguments into a dictionary.
    """
    print(f"=== Experiment: {name} ===")

    # Default parameters in CGS units
    defaults = {
        'temperature_k': 293.15,       # Kelvin
        'pressure_dynes_cm2': 1.01325e6,  # 1 atm in dynes/cm^2
        'trials': 10
    }

    # Update with user parameters
    config = defaults.copy()
    config.update(parameters)

    print("Configuration:")
    for key, value in config.items():
        print(f"  {key}: {value}")

    return f"Completed {config['trials']} trials"

# Simple call
result = run_experiment("Test A")

# Complex call with many options
result = run_experiment(
    "Test B",
    temperature_k=350,
    pressure_dynes_cm2=2e6,
    trials=50,
    catalyst="Platinum"
)

# WRONG - SyntaxError!
# def bad_function(x=10, y):  # Can't have required after optional
#     pass

# CORRECT - Required parameters first, then optional
def good_function(x, y=10):
    """Required 'x' must come before optional 'y'."""
    return x + y

# Usage
print(good_function(5))      # x=5, y=10 (default)
print(good_function(5, 20))  # x=5, y=20 (override)

# Example 1: Mixing positional-only and keyword-only (Python 3.8+)
def divide(numerator, denominator, /, *, precision=2):
    """
    numerator, denominator: must be positional
    precision: must be keyword
    """
    result = numerator / denominator
    return round(result, precision)

print(divide(10, 3))                    # Uses default precision
print(divide(10, 3, precision=4))       # Specifies precision
# divide(numerator=10, denominator=3)   # ERROR! Must be positional

# Example 2: Flexible argument collection
def process(*data_points, normalize=True, **options):
    """
    data_points: any number of positional args
    normalize: keyword-only parameter
    options: any extra keyword arguments
    """
    print(f"Processing {len(data_points)} points")
    print(f"Normalize: {normalize}")
    print(f"Options: {options}")

process(1, 2, 3, 4, normalize=False, threshold=0.5, verbose=True)

# Without keyword-only: unclear what True and False mean
def process_data_unclear(data, normalize, remove_outliers):
    """What do the boolean parameters mean at the call site?"""
    pass

# Call is confusing:
# process_data_unclear(measurements, True, False)  # What's True? What's False?

# With keyword-only: self-documenting
def process_data_clear(data, *, normalize=True, remove_outliers=False):
    """The * forces callers to name the boolean parameters."""
    if normalize:
        # Normalize the data...
        pass
    if remove_outliers:
        # Remove outliers...
        pass

# Call is clear:
# process_data_clear(measurements, normalize=True, remove_outliers=False)

# Pattern 1: Simple function with defaults
def calculate_interest(principal, rate=0.05, years=1):
    """Most common pattern - some required, some optional."""
    return principal * (1 + rate) ** years

# Pattern 2: Collecting variable arguments
def average(*numbers):
    """Accept any number of arguments."""
    if not numbers:
        return 0
    return sum(numbers) / len(numbers)

print(f"Average of 2, 4, 6, 8: {average(2, 4, 6, 8)}")

# Pattern 3: Force clarity with keyword-only
def create_particle(mass, velocity, *, charge=0, spin=0.5):
    """Mass and velocity are required, charge and spin must be named."""
    return {
        'mass': mass,
        'velocity': velocity,
        'charge': charge,
        'spin': spin
    }

# Clear what each parameter means:
electron = create_particle(9.109e-28, 1e6, charge=-1.602e-19, spin=0.5)

# Pattern 4: Flexible keyword arguments
def plot(x, y, **options):
    """Accept x, y, and any number of plot options."""
    print(f"Plotting {len(x)} points")
    for key, value in options.items():
        print(f"  {key}: {value}")

plot([1, 2, 3], [2, 4, 6], color='red', linewidth=2, marker='o')

# The / marker (Python 3.8+): everything before is positional-only
def divide_pos_only(dividend, divisor, /):
    """Both arguments MUST be passed positionally."""
    return dividend / divisor

print(divide_pos_only(10, 2))           # Works
# print(divide_pos_only(dividend=10, divisor=2))  # ERROR! Can't use keywords

# The * marker: everything after is keyword-only
def configure(name, *, debug=False, verbose=False):
    """name can be positional or keyword, others must be keyword."""
    print(f"Configuring {name}")
    if debug:
        print("  Debug mode enabled")
    if verbose:
        print("  Verbose output enabled")

configure("MyApp", debug=True)  # Must name debug
# configure("MyApp", True)  # ERROR! Can't pass debug positionally

# Demonstrating LEGB with scientific context
speed_of_light = 2.998e10  # Global scope (cm/s in CGS)

def calculate_energy():
    # Enclosing scope
    rest_mass = 9.109e-28  # electron mass in grams

    def relativistic_factor(velocity_cms):
        # Local scope
        c = speed_of_light  # Accesses global (cm/s)
        beta = velocity_cms / c

        # Local calculation
        import math
        if beta >= 1:
            return float('inf')  # Cannot exceed speed of light
        gamma = 1 / math.sqrt(1 - beta**2)

        print(f"  Inside relativistic: beta = {beta:.3f}, gamma = {gamma:.3f}")
        return gamma

    # Use inner function
    v = 0.5 * speed_of_light  # Half light speed (cm/s)
    gamma = relativistic_factor(v)

    # Calculate relativistic energy (E = gammamc^2)
    energy_ergs = rest_mass * speed_of_light**2 * gamma
    print(f"Inside calculate: E = {energy_ergs:.2e} ergs")
    return energy_ergs

result = calculate_energy()
print(f"Global scope: c = {speed_of_light:.2e} cm/s")
print(f"Final energy: {result:.2e} ergs = {result/1.602e-12:.2f} eV")

counter = 0  # Global

def increment_wrong():
    """This will crash - Python sees assignment and assumes local!"""
    # counter += 1  # UnboundLocalError!
    print("Uncomment the line above to see the error")

def increment_with_global():
    """Explicitly use global variable."""
    global counter
    counter += 1
    return counter

def increment_better(current_count):
    """Best approach - avoid global state entirely."""
    return current_count + 1

# Test the better approach
count = 0
count = increment_better(count)
count = increment_better(count)
print(f"Count: {count}")

# Global approach (generally avoid)
increment_with_global()
increment_with_global()
print(f"Global counter: {counter}")

def create_integrator(method='rectangle'):
    """
    Factory function that creates specialized integrators.

    The returned function 'remembers' either the ``rectangle`` or ``midpoint``
    method.
    """
    if method not in {'rectangle', 'midpoint'}:
        raise ValueError("method must be 'rectangle' or 'midpoint'")

    def integrate(func, a, b, n=100):
        """Numerically integrate func from a to b using a positive integer n."""
        if not isinstance(n, int) or isinstance(n, bool) or n < 1:
            raise ValueError("n must be a positive integer")
        dx = (b - a) / n
        total = 0

        if method == 'rectangle':
            # Left rectangle rule
            for i in range(n):
                x = a + i * dx
                total += func(x) * dx

        elif method == 'midpoint':
            # Midpoint rule (more accurate)
            for i in range(n):
                x = a + (i + 0.5) * dx
                total += func(x) * dx

        return total

    # Return the customized function
    return integrate

# Create specialized integrators
rect_integrate = create_integrator('rectangle')
mid_integrate = create_integrator('midpoint')

# Test with x^2 from 0 to 1 (exact answer = 1/3)
def f(x):
    return x**2

rect_result = rect_integrate(f, 0, 1, 100)
mid_result = mid_integrate(f, 0, 1, 100)
exact = 1/3

print(f"Rectangle rule: {rect_result:.6f}")
print(f"Midpoint rule:  {mid_result:.6f}")
print(f"Exact:          {exact:.6f}")
print(f"Rectangle error: {abs(rect_result-exact):.6f}")
print(f"Midpoint error:  {abs(mid_result-exact):.6f}")

# Traditional function
def square(x):
    return x**2

# Equivalent lambda
square_lambda = lambda x: x**2

print(f"Traditional: {square(5)}")
print(f"Lambda: {square_lambda(5)}")

# Lambdas shine for sorting and filtering
data_points = [
    {'x': 1.5, 'y': 2.3, 'error': 0.1},
    {'x': 2.1, 'y': 4.7, 'error': 0.05},
    {'x': 0.8, 'y': 1.2, 'error': 0.2},
]

# Sort by x value
by_x = sorted(data_points, key=lambda p: p['x'])
print("\nSorted by x:")
for p in by_x:
    print(f"  x={p['x']}, y={p['y']}")

# Sort by relative error (error/y)
by_rel_error = sorted(data_points, key=lambda p: p['error']/p['y'])
print("\nSorted by relative error:")
for p in by_rel_error:
    rel_err = p['error']/p['y'] * 100
    print(f"  x={p['x']}, relative error={rel_err:.1f}%")

import math
from functools import reduce

# Sample measurements in CGS units
measurements_cms = [981, 979, 983, 9900, 980, 982, 978]  # cm/s^2 (g values)

# FILTER: Remove outliers
mean_estimate = 980  # cm/s^2
tolerance = 10  # cm/s^2

good_data = list(filter(
    lambda x: abs(x - mean_estimate) < tolerance,
    measurements_cms
))

print(f"Original: {measurements_cms}")
print(f"Filtered: {good_data}")
print(f"Removed {len(measurements_cms) - len(good_data)} outliers")

# MAP: Apply calibration factor
calibration = 1.002  # 0.2% correction
calibrated = list(map(lambda x: x * calibration, good_data))
print(f"Calibrated: {[f'{x:.1f}' for x in calibrated]}")

# REDUCE: Calculate mean
# Guard clause: filtering might remove everything
if not calibrated:
    raise ValueError("No data left after filtering; cannot compute mean")

# In real code, prefer math.fsum(...) or sum(...) over reduce(...)
total = reduce(lambda a, b: a + b, calibrated, 0.0)
mean = total / len(calibrated)
print(f"Final mean: {mean:.1f} cm/s^2")

# In array frameworks such as JAX, transformations act most predictably on
# functions whose inputs and outputs make state explicit. Python's built-in
# map() does not itself promise GPU execution or parallelism.

def apply_operator(data, operator_func):
    """
    Apply any operation to data.

    operator_func determines the transformation.
    """
    return [operator_func(x) for x in data]

# Define different operations
def grams_to_kilograms(g):
    """Convert grams to kilograms."""
    return g / 1000.0

def dynes_to_newtons(dynes):
    """Convert dynes to Newtons (1 N = 10^5 dynes)."""
    return dynes / 1e5

def normalize(x, mean=0, std=1):
    """Standardize value."""
    return (x - mean) / std

# Use different operations
masses_g = [100, 200, 300, 400]  # grams
masses_kg = apply_operator(masses_g, grams_to_kilograms)
print(f"Grams: {masses_g}")
print(f"Kilograms: {masses_kg}")

forces_dynes = [1e5, 2e5, 3e5, 4e5]  # dynes
forces_n = apply_operator(forces_dynes, dynes_to_newtons)
print(f"\nDynes: {forces_dynes}")
print(f"Newtons: {forces_n}")

# Create custom normalizer with closure
values = [2, 4, 6, 8, 10]
mean = sum(values) / len(values)
std = (sum((x - mean)**2 for x in values) / len(values)) ** 0.5

normalizer = lambda x: normalize(x, mean, std)
normalized = apply_operator(values, normalizer)
print(f"\nOriginal: {values}")
print(f"Normalized: {[f'{x:.2f}' for x in normalized]}")

def integrate_trapezoid(x_values, y_values):
    """
    Integrate using the trapezoidal rule.

    Parameters
    ----------
    x_values : list or array
        X coordinates of data points (strictly increasing)
    y_values : list or array
        Y coordinates of data points

    Returns
    -------
    float
        Approximate integral of y with respect to x

    Examples
    --------
    >>> x = [0, 1, 2]
    >>> y = [0, 1, 4]  # y = x^2
    >>> integrate_trapezoid(x, y)
    3.0

    With only these two intervals, `3.0` is the trapezoidal *approximation*.
    The exact integral of $x^2$ from 0 to 2 is $8/3 \approx 2.667$.

    Notes
    -----
    Uses the composite trapezoidal rule:
    Integral y dx approx sum(y[i] + y[i+1])/2 * (x[i+1] - x[i])

    More accurate for smooth functions with
    many points.
    """
    if len(x_values) != len(y_values):
        raise ValueError("x and y must have same length")
    if len(x_values) < 2:
        raise ValueError("trapezoid integration requires at least two points")
    if any(right <= left for left, right in zip(x_values, x_values[1:])):
        raise ValueError("x_values must be strictly increasing")

    integral = 0.0
    for i in range(len(x_values) - 1):
        dx = x_values[i + 1] - x_values[i]
        avg_y = (y_values[i] + y_values[i + 1]) / 2
        integral += avg_y * dx

    return integral

# Test the function
x = [0, 0.5, 1.0, 1.5, 2.0]
y = [0, 0.25, 1.0, 2.25, 4.0]  # y = x^2
result = integrate_trapezoid(x, y)
print(f"Integral of x^2 from 0 to 2: {result:.3f}")
print(f"Exact answer: {8/3:.3f}")
print(f"Error: {abs(result - 8/3):.3f}")

def quadratic_formula(a, b, c):
    """
    Solve quadratic equation ax^2 + bx + c = 0.

    Returns:
        tuple: (root1, root2) or (None, None) if no real roots
    """
    import math

    if a == 0:
        raise ValueError("Not a quadratic: a must be nonzero")

    # Calculate discriminant
    discriminant = b**2 - 4*a*c

    # Check if real roots exist
    if discriminant < 0:
        return None, None

    sqrt_disc = math.sqrt(discriminant)
    if discriminant == 0:
        root = -b / (2 * a)
        return root, root

    # Numerically stable quadratic formula:
    # Avoids catastrophic cancellation when |b| is large.
    q = -0.5 * (b + math.copysign(sqrt_disc, b))
    root1 = q / a
    root2 = c / q
    return root1, root2

def test_quadratic():
    """Test quadratic formula solver."""
    import math

    # Test 1: Simple case (x^2 - 5x + 6 = 0)
    # Roots should be 2 and 3
    r1, r2 = quadratic_formula(1, -5, 6)
    roots = sorted([r1, r2])
    assert math.isclose(roots[0], 2.0, rel_tol=0.0, abs_tol=1e-10)
    assert math.isclose(roots[1], 3.0, rel_tol=0.0, abs_tol=1e-10)
    for r in (r1, r2):
        assert math.isclose(1 * r * r + (-5) * r + 6, 0.0, rel_tol=0.0, abs_tol=1e-9)
    print("Test 1: Simple roots - PASSED")

    # Test 2: No real roots (x^2 + 1 = 0)
    r1, r2 = quadratic_formula(1, 0, 1)
    assert r1 is None and r2 is None
    print("Test 2: No real roots - PASSED")

    # Test 3: Single root (x^2 - 2x + 1 = 0)
    # Root should be 1 (twice)
    r1, r2 = quadratic_formula(1, -2, 1)
    assert math.isclose(r1, 1.0, rel_tol=0.0, abs_tol=1e-10)
    assert math.isclose(r2, 1.0, rel_tol=0.0, abs_tol=1e-10)
    print("Test 3: Repeated root - PASSED")

    print("All tests passed!")

# Run tests
test_quadratic()

def kinetic_energy_typed(mass_g: float, velocity_cms: float) -> float:
    """
    Calculate kinetic energy with type hints.

    The ': float' annotations indicate expected types.
    The '-> float' indicates the return type.
    These are documentation only - Python doesn't enforce them!
    """
    return 0.5 * mass_g * velocity_cms**2

# Type hints help IDEs provide better autocomplete and catch potential errors
energy: float = kinetic_energy_typed(9.109e-28, 3e9)
print(f"Energy: {energy:.2e} ergs")

# Python does not check annotations before entering the function. Integers are
# accepted where floats are annotated in the usual numeric type relationship:
energy2 = kinetic_energy_typed(1, 100)  # Still works!
print(f"Integer inputs are accepted: {energy2} ergs")

# For more complex types, import from typing module
from typing import List, Tuple, Optional

def analyze_data_typed(values: List[float]) -> Tuple[float, float, Optional[float]]:
    """
    Analyze data with complex type hints.

    List[float] means a list of floats
    Tuple[float, float, Optional[float]] means it returns
    a tuple with two floats and maybe a third float (or None)
    """
    if not values:
        return 0.0, 0.0, None

    n = len(values)
    mean = sum(values) / n
    variance = sum((x - mean)**2 for x in values) / n
    std_dev = variance ** 0.5

    sorted_vals = sorted(values)
    median = sorted_vals[n // 2] if n % 2 == 1 else (sorted_vals[n//2-1] + sorted_vals[n//2]) / 2

    return mean, std_dev, median

# The IDE now knows exactly what types are returned
mean_val, std_val, median_val = analyze_data_typed([1.0, 2.0, 3.0, 4.0, 5.0])
print(f"Mean: {mean_val:.1f}, Std: {std_val:.2f}, Median: {median_val:.1f}")

import timeit

def simple_add(a, b):
    """Minimal work - overhead visible."""
    return a + b

def complex_calc(x):
    """Heavy computation - overhead negligible."""
    result = 0
    for i in range(100):
        result += (x + i) ** 0.5
    return result / 100

# Measure overhead with timeit (more reliable than manual timing)
simple_time = timeit.timeit("simple_add(1.0, 2.0)", globals=globals(), number=200_000)
inline_time = timeit.timeit("1.0 + 2.0", number=200_000)
complex_time = timeit.timeit("complex_calc(1.0)", globals=globals(), number=5_000)

print("Function Call Analysis:")
print(f"Simple function: {simple_time*1000:.1f} ms")
print(f"Inline addition: {inline_time*1000:.1f} ms")
print(f"Overhead factor: {simple_time/inline_time:.1f}x")
print(f"\nComplex function: {complex_time*1000:.1f} ms")
print("\nThis is one interpreter and workload; profile the calculation that matters.")

from functools import lru_cache
import time
import math

# Fibonacci without memoization (exponentially slow!)
def fib_slow(n):
    if n < 2:
        return n
    return fib_slow(n-1) + fib_slow(n-2)

# Fibonacci with memoization (linear time!)
@lru_cache(maxsize=128)
def fib_fast(n):
    if n < 2:
        return n
    return fib_fast(n-1) + fib_fast(n-2)

# Add astrophysics-specific memoization example
@lru_cache(maxsize=128)
def _planck_function_cached(wavelength_nm, temp_k):
    """
    Cached Planck function for blackbody radiation.
    Expensive calculation cached for repeated calls.
    Common in spectral fitting where same temps recur.
    """
    h = 6.626e-27  # Planck constant in erg*s
    c = 2.998e10   # Speed of light in cm/s
    k = 1.381e-16  # Boltzmann constant in erg/K

    wavelength_cm = wavelength_nm * 1e-7

    # Planck function calculation
    exp_term = (h * c) / (wavelength_cm * k * temp_k)
    if exp_term > 700:  # Prevent overflow
        return 0

    numerator = 2 * h * c**2 / wavelength_cm**5
    denominator = math.expm1(exp_term)  # stable for small exp_term
    return numerator / denominator

# Cache the exact inputs supplied to this function. For a model-defined grid,
# use grid indices as the cache keys and document the grid resolution separately.
def planck_function(wavelength_nm, temp_k):
    return _planck_function_cached(wavelength_nm, temp_k)

# Compare performance
n = 28

start = time.perf_counter()
result_slow = fib_slow(n)
time_slow = time.perf_counter() - start

start = time.perf_counter()
result_fast = fib_fast(n)
time_fast = time.perf_counter() - start

print(f"Fibonacci({n}) = {result_fast}")
print(f"Without memoization: {time_slow*1000:.1f} ms")
print(f"With memoization: {time_fast*1000:.3f} ms")
print(f"Local ratio for this first call: {time_slow/time_fast:.0f}x")

# Check cache statistics
print(f"\nCache info: {fib_fast.cache_info()}")

# Test Planck function caching
print(f"\nPlanck function at 500nm, 5778K: {planck_function(500, 5778):.2e} erg/(s cm^2 cm sr)")
print(f"Cache info: {_planck_function_cached.cache_info()}")