Chapter 7 — Computational Experiments

# --- Visualization: Algorithm runtime experiment with variance ---
# A proper experiment measures not just mean runtime but distribution.
import numpy as np
import matplotlib.pyplot as plt
import time
plt.style.use('seaborn-v0_8-whitegrid')

def insertion_sort(arr):
    a = arr.copy()
    for i in range(1, len(a)):
        key, j = a[i], i - 1
        while j >= 0 and a[j] > key:
            a[j+1] = a[j]; j -= 1
        a[j+1] = key
    return a

N_SIZES = [50, 100, 200, 400, 800]
N_REPS = 30
INPUT_TYPES = {'random': lambda n: np.random.permutation(n),
               'sorted': lambda n: np.arange(n),
               'reverse': lambda n: np.arange(n)[::-1]}

results = {t: {n: [] for n in N_SIZES} for t in INPUT_TYPES}

for input_type, gen in INPUT_TYPES.items():
    for n in N_SIZES:
        for _ in range(N_REPS):
            arr = gen(n)
            t0 = time.perf_counter()
            insertion_sort(arr)
            results[input_type][n].append(time.perf_counter() - t0)

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

colors = {'random': 'steelblue', 'sorted': 'green', 'reverse': 'coral'}
for input_type, color in colors.items():
    means = [np.mean(results[input_type][n]) * 1e6 for n in N_SIZES]
    stds  = [np.std(results[input_type][n])  * 1e6 for n in N_SIZES]
    axes[0].plot(N_SIZES, means, 'o-', color=color, label=input_type, linewidth=2)
    axes[0].fill_between(N_SIZES,
                          [m-s for m,s in zip(means,stds)],
                          [m+s for m,s in zip(means,stds)],
                          alpha=0.2, color=color)

axes[0].set_xlabel('Input size n'); axes[0].set_ylabel('Runtime (μs)')
axes[0].set_title('Insertion Sort: Runtime by Input Type (mean ± 1 std)')
axes[0].legend()

# Log-log to check O(n²) scaling
for input_type, color in colors.items():
    means = [np.mean(results[input_type][n]) * 1e6 for n in N_SIZES]
    axes[1].loglog(N_SIZES, means, 'o-', color=color, label=input_type, linewidth=2)

n_ref = np.array(N_SIZES, dtype=float)
axes[1].loglog(n_ref, n_ref**2 / n_ref[-1]**2 * np.mean(results['random'][N_SIZES[-1]]) * 1e6,
               'k--', alpha=0.5, label='O(n²) reference')
axes[1].loglog(n_ref, n_ref / n_ref[-1] * np.mean(results['sorted'][N_SIZES[-1]]) * 1e6,
               'k:', alpha=0.5, label='O(n) reference')
axes[1].set_xlabel('n (log)'); axes[1].set_ylabel('Runtime μs (log)')
axes[1].set_title('Log-Log: Confirm O(n²) average, O(n) best case')
axes[1].legend(fontsize=8)

plt.tight_layout()
plt.show()

# --- Implementation: Experiment framework ---
import numpy as np
import time

def run_experiment(algorithm, input_generator, param_values, n_reps=20):
    """
    Run a controlled timing experiment on an algorithm.
    
    Args:
        algorithm:       callable(input) -> any (we measure its runtime)
        input_generator: callable(param) -> input
        param_values:    list of parameter values to test
        n_reps:          number of replications per parameter value
    Returns:
        dict: param -> {'mean': float, 'std': float, 'snr': float}
    """
    results = {}
    for p in param_values:
        times = []
        for _ in range(n_reps):
            inp = input_generator(p)
            t0 = time.perf_counter()
            algorithm(inp)
            times.append(time.perf_counter() - t0)
        times = np.array(times)
        results[p] = {
            'mean': times.mean(),
            'std': times.std(),
            'snr': times.mean() / (times.std() + 1e-15),
        }
    return results

def estimate_complexity(param_values, mean_times):
    """
    Estimate algorithmic complexity exponent from timing data.
    Fits log(T) = alpha * log(n) + const on log-log scale.
    """
    log_n = np.log(param_values)
    log_t = np.log(mean_times)
    alpha = np.polyfit(log_n, log_t, 1)[0]
    return alpha

# Experiment: merge sort
def merge_sort(arr):
    if len(arr) <= 1: return arr
    mid = len(arr) // 2
    L, R = merge_sort(arr[:mid]), merge_sort(arr[mid:])
    result = []
    i = j = 0
    while i < len(L) and j < len(R):
        if L[i] <= R[j]: result.append(L[i]); i += 1
        else: result.append(R[j]); j += 1
    return result + L[i:] + R[j:]

SIZES = [100, 200, 500, 1000, 2000]
np.random.seed(5)
exp_results = run_experiment(
    algorithm=merge_sort,
    input_generator=lambda n: list(np.random.permutation(n)),
    param_values=SIZES,
    n_reps=15
)

means = [exp_results[n]['mean'] for n in SIZES]
alpha = estimate_complexity(SIZES, means)

print(f"Merge sort complexity estimate: O(n^{alpha:.3f})")
print(f"Expected: O(n log n) ≈ O(n^1.0) to O(n^1.1) for practical sizes")
print()
for n in SIZES:
    r = exp_results[n]
    print(f"  n={n:5d}  mean={r['mean']*1e6:.1f}μs  std={r['std']*1e6:.1f}μs  SNR={r['snr']:.1f}")

# --- Experiment: Confirmation bias in algorithm testing ---
# Hypothesis: Testing only "nice" inputs gives systematically optimistic results.
# Try changing: SEED to see how much results vary across random seeds

import numpy as np
SEED = 42  # <-- try: 1, 2, 3, 100
np.random.seed(SEED)

# A flawed binary search that fails on duplicate values
def buggy_binary_search(arr, target):
    """Binary search — fails silently when duplicates exist."""
    lo, hi = 0, len(arr) - 1
    while lo <= hi:
        mid = (lo + hi) // 2
        if arr[mid] == target: return mid
        elif arr[mid] < target: lo = mid + 1
        else: hi = mid - 1
    return -1

# Test 1: Only unique-value arrays ("nice" testing)
n_nice_fails = 0
for _ in range(1000):
    arr = sorted(np.random.choice(10000, 50, replace=False))
    target = arr[np.random.randint(len(arr))]
    idx = buggy_binary_search(arr, target)
    if arr[idx] != target: n_nice_fails += 1

# Test 2: Arrays with duplicates (adversarial testing)
n_dup_fails = 0
for _ in range(1000):
    arr = sorted(np.random.choice(20, 50, replace=True))  # high chance of duplicates
    target = arr[np.random.randint(len(arr))]
    idx = buggy_binary_search(arr, target)
    if idx == -1 or arr[idx] != target: n_dup_fails += 1

print(f"Failure rate, unique inputs:     {n_nice_fails/1000:.1%}")
print(f"Failure rate, duplicate inputs:  {n_dup_fails/1000:.1%}")
print()
print("Testing only nice inputs misses the bug entirely.")
print("A well-designed experiment includes adversarial cases.")

# --- Mini Project: Sorting Algorithm Tournament ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

algorithms = {
    'insertion_sort': insertion_sort,
    'merge_sort': lambda arr: merge_sort(list(arr)),
    'python_sorted': sorted,
}

distributions = {
    'random':        lambda n: np.random.permutation(n),
    'nearly_sorted': lambda n: np.array(sorted(np.random.permutation(n))[:-int(n*0.1)] +
                                        list(np.random.permutation(int(n*0.1)))),
    'reverse':       lambda n: np.arange(n)[::-1],
}

SIZES = [100, 500, 1000]
N_REPS = 15
np.random.seed(7)

print(f"{'Algorithm':<18} {'Distribution':<16}", end='')
for n in SIZES:
    print(f" n={n:>5} (μs)", end='')
print()
print("-" * 80)

for alg_name, alg in algorithms.items():
    for dist_name, dist_gen in distributions.items():
        print(f"{alg_name:<18} {dist_name:<16}", end='')
        for n in SIZES:
            times = []
            for _ in range(N_REPS):
                arr = dist_gen(n)
                t0 = time.perf_counter()
                alg(arr)
                times.append(time.perf_counter() - t0)
            print(f" {np.mean(times)*1e6:>11.1f}", end='')
        print()
    print()