NumPy Loop vs Vectorisation: 45-Minute Training Bottleneck

# io.thecodeforge.numpy.loop_cost
# Demonstrates the interpreter overhead cost of Python loops vs vectorised operations
# Run with: python loop_cost.py or in a Jupyter notebook with %timeit

import numpy as np
import time

arr = np.random.randn(1_000_000)  # 1M elements, float64

# ---- Python loop: pays interpreter overhead on every element ----
start = time.perf_counter()
total_loop = 0.0
for x in arr:
    total_loop += x
elapsed_loop = time.perf_counter() - start
print(f'Python loop:     {elapsed_loop:.3f}s  (result={total_loop:.4f})')

# ---- Vectorised: hands a C pointer to a compiled SIMD loop ----
start = time.perf_counter()
total_vec = arr.sum()
elapsed_vec = time.perf_counter() - start
print(f'NumPy vectorised: {elapsed_vec:.4f}s  (result={total_vec:.4f})')

print(f'Speedup: {elapsed_loop / elapsed_vec:.0f}x')
# Typical output on a modern CPU:
# Python loop:      0.287s
# NumPy vectorised: 0.0019s
# Speedup: ~150x

# ---- Per-element overhead breakdown ----
# Python loop: ~300ns per element
#   Boxing scalar to Python float:    ~80ns
#   Bytecode execution (BINARY_ADD):  ~100ns
#   Reference count increment:         ~40ns
#   Temporary object deallocation:     ~50ns
#   Loop bookkeeping (JUMP_ABSOLUTE):  ~30ns
#
# Vectorised: ~0.2ns per element
#   Raw memory access from C pointer:  ~0.1ns
#   SIMD add (AVX2 processes 4 doubles per instruction): ~0.1ns
#   Zero Python overhead per element

# ---- The overhead is fixed, not proportional ----
# For arr + 1 (trivial), interpreter overhead is 99% of runtime
# For np.sin(x) (complex), overhead shrinks because C computation dominates

# Simple operation: loop vs vectorised gap is largest
start = time.perf_counter()
result_loop = np.array([x + 1.0 for x in arr])
elapsed_add_loop = time.perf_counter() - start

start = time.perf_counter()
result_vec = arr + 1.0
elapsed_add_vec = time.perf_counter() - start

print(f'\nElement-wise add:')
print(f'  Loop:        {elapsed_add_loop:.3f}s')
print(f'  Vectorised:  {elapsed_add_vec:.4f}s')
print(f'  Speedup:     {elapsed_add_loop / elapsed_add_vec:.0f}x')
# Speedup for addition is typically 100-150x — overhead dominates

# Complex operation: overhead shrinks in proportion
start = time.perf_counter()
result_sin_loop = np.array([float(np.sin(x)) for x in arr])
elapsed_sin_loop = time.perf_counter() - start

start = time.perf_counter()
result_sin_vec = np.sin(arr)
elapsed_sin_vec = time.perf_counter() - start

print(f'\nElement-wise np.sin:')
print(f'  Loop:        {elapsed_sin_loop:.3f}s')
print(f'  Vectorised:  {elapsed_sin_vec:.4f}s')
print(f'  Speedup:     {elapsed_sin_loop / elapsed_sin_vec:.0f}x')
# Speedup for sin is typically 10-20x — C computation now significant

# io.thecodeforge.numpy.patterns
# Direct replacements for the most common Python loop patterns over NumPy arrays
# Each pattern shown as: loop version (slow) -> vectorised version (fast)

import numpy as np

arr = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])

# ---- Pattern 1: Adjacent differences ----
# Loop (slow):
diffs_loop = [arr[i+1] - arr[i] for i in range(len(arr) - 1)]
# Vectorised:
diffs = np.diff(arr)
print(f'np.diff:    {diffs}')  # [1. 1. 1. 1. 1. 1. 1.]

# ---- Pattern 2: Cumulative sum ----
# Loop (slow):
cumulative_loop = [sum(arr[:i+1]) for i in range(len(arr))]
# Vectorised:
cumulative = np.cumsum(arr)
print(f'np.cumsum:  {cumulative}')  # [ 1.  3.  6. 10. 15. 21. 28. 36.]

# ---- Pattern 3: Value clamping ----
# Loop (slow):
clipped_loop = [min(max(x, 2.0), 6.0) for x in arr]
# Vectorised:
clipped = np.clip(arr, 2.0, 6.0)
print(f'np.clip:    {clipped}')  # [2. 2. 3. 4. 5. 6. 6. 6.]

# ---- Pattern 4: Conditional assignment — two branches ----
# Loop (slow):
result_loop = np.empty_like(arr)
for i in range(len(arr)):
    if arr[i] > 4.0:
        result_loop[i] = arr[i] * 10
    else:
        result_loop[i] = arr[i] * 0.5

# Vectorised with np.where:
result_where = np.where(arr > 4.0, arr * 10, arr * 0.5)
print(f'np.where:   {result_where}')  # [ 0.5  1.   1.5  2.   50.  60.  70.  80. ]

# Vectorised with boolean masking (preferred when branches have very different costs):
# Boolean masking only evaluates the branch for the subset that matches
result_mask = np.empty_like(arr)
mask = arr > 4.0
result_mask[mask] = arr[mask] * 10
result_mask[~mask] = arr[~mask] * 0.5
print(f'bool mask:  {result_mask}')   # same result

# Note: np.where(cond, x, y) evaluates BOTH x and y for all elements,
# then selects. Boolean masking evaluates each branch only for the matching
# subset — more efficient when one branch is expensive and rarely true.

# ---- Pattern 5: Multiple branches — np.select ----
# Extends np.where to N conditions
conditions = [arr < 2.0, (arr >= 2.0) & (arr < 5.0), arr >= 5.0]
choices = [arr * 0.0, arr * 1.0, arr * 2.0]  # zero, keep, double
result_select = np.select(conditions, choices, default=arr)
print(f'np.select:  {result_select}')  # [0. 0. 3. 4. 10. 12. 14. 16.]

# ---- Pattern 6: Rolling window statistics (the key production pattern) ----
# Loop (slow) — creates 1 temporary array per iteration:
# for i in range(window, len(arr)):
#     rolling_mean[i] = arr[i-window:i].mean()

def rolling_stats(arr, window):
    """Rolling mean and std using cumulative sums. O(n), no Python loop."""
    # Pad with zero at position 0 so indexing arithmetic works cleanly
    cs = np.cumsum(np.insert(arr, 0, 0.0))    # prefix sum
    cs2 = np.cumsum(np.insert(arr**2, 0, 0.0))  # prefix sum of squares

    rolling_sum = cs[window:] - cs[:-window]
    rolling_sum2 = cs2[window:] - cs2[:-window]

    mean = rolling_sum / window
    # Computational formula: Var = E[X^2] - E[X]^2
    # Clip to 0 to avoid negative values from floating-point precision loss
    var = np.maximum(rolling_sum2 / window - mean**2, 0.0)
    std = np.sqrt(var)

    return mean, std

mean, std = rolling_stats(arr, window=3)
print(f'Rolling mean (w=3): {mean}')  # [2. 3. 4. 5. 6. 7.]
print(f'Rolling std  (w=3): {std}')   # [~0.816 ~0.816 ...]

# ---- Pattern 7: Normalisation across rows of a 2D array ----
# Loop (slow):
matrix = np.random.randn(1000, 50)  # 1000 rows, 50 features
norm_loop = np.empty_like(matrix)
for i in range(len(matrix)):
    norm_loop[i] = (matrix[i] - matrix[i].mean()) / matrix[i].std()

# Vectorised with broadcasting — keepdims preserves shape for broadcasting:
means = matrix.mean(axis=1, keepdims=True)  # shape (1000, 1)
stds = matrix.std(axis=1, keepdims=True)    # shape (1000, 1)
norm_vec = (matrix - means) / stds          # broadcasts across 50 columns
print(f'Row normalisation matches: {np.allclose(norm_loop, norm_vec)}')

# io.thecodeforge.numpy.vectorize_vs_numba
# Demonstrates:
# 1. np.vectorize — what it actually does (convenience, not speed)
# 2. Numba @njit — genuine compiled performance for sequential logic
# 3. Side-by-side benchmark showing the difference

import numpy as np
import time

rng = np.random.default_rng(42)
arr = rng.standard_normal(10_000_000).astype(np.float64)  # 10M elements

# ---- np.vectorize: the convenience wrapper that is NOT faster ----
# Under the hood: Python loop + per-element type coercion
# Use this when you want array broadcasting on a scalar function
# Do NOT use this when you want speed

def clip_scalar(x, lo, hi):
    """Works on scalars. np.vectorize will call this once per element."""
    return lo if x < lo else (hi if x > hi else x)

vclip = np.vectorize(clip_scalar)

# For a small demo array:
small = np.array([-0.5, 0.3, 0.7, 1.5, -0.1])
print(f'np.vectorize result: {vclip(small, 0.0, 1.0)}')
# [0.  0.3 0.7 1.  0. ]

# Benchmark: np.vectorize vs raw loop vs numpy built-in
start = time.perf_counter()
_ = vclip(arr, 0.0, 1.0)
time_vectorize = time.perf_counter() - start

start = time.perf_counter()
result_py = np.empty_like(arr)
for i, x in enumerate(arr):
    result_py[i] = 0.0 if x < 0.0 else (1.0 if x > 1.0 else x)
time_loop = time.perf_counter() - start

start = time.perf_counter()
_ = np.clip(arr, 0.0, 1.0)  # built-in ufunc
time_numpy = time.perf_counter() - start

print(f'np.vectorize: {time_vectorize:.2f}s')
print(f'Python loop:  {time_loop:.2f}s')
print(f'np.clip:      {time_numpy:.4f}s')
# Typical output:
# np.vectorize: 4.87s  <- same order as the loop, sometimes worse
# Python loop:  4.61s  <- the baseline we're trying to beat
# np.clip:      0.031s <- 150x faster — this is what vectorisation means

# ---- Numba @njit: compiled sequential logic ----
# Use this when the loop has genuine sequential dependencies
# (result[i] depends on result[i-1]) and cannot be vectorised

try:
    from numba import njit

    @njit
    def running_max(arr):
        """
        Running maximum: result[i] = max(arr[0], arr[1], ..., arr[i]).
        Genuinely sequential — result[i] depends on result[i-1].
        Cannot be vectorised with standard NumPy.
        """
        result = np.empty_like(arr)
        result[0] = arr[0]
        for i in range(1, len(arr)):
            result[i] = result[i-1] if result[i-1] > arr[i] else arr[i]
        return result

    @njit
    def exponential_smooth(arr, alpha):
        """
        Exponential weighted moving average (EWMA).
        result[i] = alpha * arr[i] + (1 - alpha) * result[i-1]
        Classic sequential recurrence — each output depends on the prior output.
        """
        result = np.empty_like(arr)
        result[0] = arr[0]
        for i in range(1, len(arr)):
            result[i] = alpha * arr[i] + (1.0 - alpha) * result[i-1]
        return result

    # Warmup call — triggers LLVM compilation (~0.5-2s one-time cost)
    _ = running_max(arr[:100])
    _ = exponential_smooth(arr[:100], 0.1)

    # Benchmark compiled Numba vs equivalent Python loop
    start = time.perf_counter()
    result_nb = running_max(arr)
    time_numba = time.perf_counter() - start

    # Python equivalent for comparison
    start = time.perf_counter()
    result_py2 = np.empty_like(arr)
    result_py2[0] = arr[0]
    for i in range(1, len(arr)):
        result_py2[i] = max(result_py2[i-1], arr[i])
    time_python2 = time.perf_counter() - start

    print(f'\nRunning max (10M elements):')
    print(f'  Python loop: {time_python2:.2f}s')
    print(f'  Numba @njit: {time_numba:.4f}s')
    print(f'  Speedup:     {time_python2 / time_numba:.0f}x')
    print(f'  Results match: {np.allclose(result_nb, result_py2)}')
    # Typical output:
    # Python loop: 3.84s
    # Numba @njit: 0.022s
    # Speedup:     ~175x

except ImportError:
    print('Numba not installed — pip install numba')
    print('For environments without Numba, Cython is an alternative')
    print('For pandas integration, consider bottleneck or scipy alternatives')

# ---- Summary: what to use when ----
# np.vectorize: when you want to apply a scalar Python function element-wise
#               for API convenience — accept Python loop speed
# np.clip / np.where / ufuncs: when the operation has a NumPy equivalent
#               — 100-300x faster than any loop
# Numba @njit: when the loop has genuine sequential dependencies
#               — 50-200x faster than Python loop, near-C speed
# Cython: when team already uses it and prefers explicit type annotations
#               — similar performance to Numba, more setup overhead

# io.thecodeforge.numpy.broadcasting_memory
# Broadcasting efficiency, memory layout effects, and in-place operations
# These factors affect performance after you've already eliminated the Python loop

import numpy as np
import time

# ---- Broadcasting: no copies for shape-compatible operations ----
rng = np.random.default_rng(0)
matrix = rng.standard_normal((10_000, 500))  # 10K rows, 500 columns

# Row normalisation via broadcasting — no temporary (1000, 500) copy of means
means = matrix.mean(axis=1, keepdims=True)  # shape: (10000, 1)
stds = matrix.std(axis=1, keepdims=True)    # shape: (10000, 1)
# The subtraction broadcasts (10000, 1) across 500 columns without materialising
# a full (10000, 500) copy of the means
start = time.perf_counter()
norm = (matrix - means) / stds
time_broadcast = time.perf_counter() - start
print(f'Broadcasting normalisation: {time_broadcast:.4f}s')

# Equivalent with explicit expansion — forces a copy allocation
means_expanded = np.tile(means, (1, 500))   # shape: (10000, 500) — full copy
start = time.perf_counter()
norm_tiled = (matrix - means_expanded) / stds
time_tiled = time.perf_counter() - start
print(f'Tiled (explicit copy):      {time_tiled:.4f}s')
print(f'Broadcasting speedup:       {time_tiled / time_broadcast:.1f}x')
# Broadcasting is typically 1.5-3x faster because it avoids the tile allocation

# ---- Memory layout: C-order vs Fortran-order ----
# C-order: rows are contiguous — efficient for row-wise operations
# F-order: columns are contiguous — efficient for column-wise operations

c_arr = np.ascontiguousarray(rng.standard_normal((5000, 5000)))  # C-order
f_arr = np.asfortranarray(c_arr.copy())                          # F-order

# Row sum: accessing contiguous elements in C-order
start = time.perf_counter()
for _ in range(5):
    _ = c_arr.sum(axis=1)  # sum across columns per row
time_c_row = (time.perf_counter() - start) / 5

start = time.perf_counter()
for _ in range(5):
    _ = f_arr.sum(axis=1)  # same operation, Fortran layout — cache-unfriendly
time_f_row = (time.perf_counter() - start) / 5

print(f'\nRow sum, C-order:        {time_c_row:.4f}s')
print(f'Row sum, Fortran-order:  {time_f_row:.4f}s')
print(f'Layout speedup (row op): {time_f_row / time_c_row:.1f}x')
# C-order is typically 2-5x faster for row operations due to cache locality

# Verify memory layout:
print(f'\nc_arr is C-contiguous: {c_arr.flags["C_CONTIGUOUS"]}')
print(f'f_arr is F-contiguous: {f_arr.flags["F_CONTIGUOUS"]}')

# ---- In-place operations: avoid intermediate allocations ----
large = rng.standard_normal(10_000_000).astype(np.float64)

# Out-of-place: allocates a 80MB temporary array
start = time.perf_counter()
result_oop = large + 1.0  # allocates new array
time_oop = time.perf_counter() - start

# In-place: no allocation — modifies large directly
large_copy = large.copy()  # preserve original for fair comparison
start = time.perf_counter()
large_copy += 1.0  # modifies in place, no allocation
time_inplace = time.perf_counter() - start

print(f'\nOut-of-place (large + 1.0): {time_oop:.4f}s')
print(f'In-place (large += 1.0):    {time_inplace:.4f}s')
print(f'In-place speedup:           {time_oop / time_inplace:.1f}x')
# In-place is typically 1.5-2x faster for large arrays
# because it avoids memory allocation and halves the memory bandwidth required

# ---- Practical rule: check array flags before assuming layout ----
print(f'\nCheck any array layout with arr.flags:')
print(large.flags)

// io.thecodeforge
import numpy as np
import time

# Contiguous array
arr = np.random.rand(10_000_000)

# Non-contiguous view via slicing every other element
sliced = arr[::2]  # strides increase, not contiguous

# Time sum on contiguous
t0 = time.time()
_ = arr.sum()
t1 = time.time()
print(f"Contiguous sum: {t1-t0:.4f}s")

# Time sum on non-contiguous
t0 = time.time()
_ = sliced.sum()
t1 = time.time()
print(f"Sliced view sum: {t1-t0:.4f}s")

# Fix: copy to force contiguity
fixed = sliced.copy()
t0 = time.time()
_ = fixed.sum()
t1 = time.time()
print(f"Copied contiguous: {t1-t0:.4f}s")

// io.thecodeforge
import numpy as np

# Two arrays for sine computation
angles = np.arange(0, np.pi, 0.1, dtype=np.float32)

# Float32 sine – watch precision loss
sin32 = np.sin(angles)
print(f"float32 output sample: {sin32[0]:.10f}")

# Force float64 output
sin64 = np.sin(angles, dtype=np.float64)
print(f"float64 output sample: {sin64[0]:.10f}")

# Reuse buffer to avoid allocation
buffer = np.empty_like(angles, dtype=np.float64)
np.sin(angles, out=buffer)
print(f"Buffer reuse matches: {np.allclose(sin64, buffer)}")