Chapter 147 — NumPy Vector Operations

# --- Visualization: axis conventions for reductions ---

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches

plt.style.use('seaborn-v0_8-whitegrid')

A = np.array([[1, 2, 3, 4],
              [5, 6, 7, 8],
              [9,10,11,12]], dtype=float)

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

def draw_matrix_with_highlight(ax, M, highlight_axis, title, result):
    rows, cols = M.shape
    ax.imshow(M, cmap='Blues', vmin=0, vmax=M.max()+2, aspect='auto')
    for i in range(rows):
        for j in range(cols):
            ax.text(j, i, f'{M[i,j]:.0f}', ha='center', va='center', fontsize=12)

    if highlight_axis == 0:  # collapse rows -> highlight columns
        for j in range(cols):
            rect = patches.Rectangle((j-0.5, -0.5), 1, rows, lw=2,
                                      edgecolor='red', facecolor='none')
            ax.add_patch(rect)
        for j, val in enumerate(result):
            ax.text(j, rows+0.3, f'{val:.0f}', ha='center', va='center',
                    fontsize=11, color='red', fontweight='bold')
        ax.text(-0.8, rows+0.3, 'sum→', fontsize=10, color='red')
    else:  # collapse cols -> highlight rows
        for i in range(rows):
            rect = patches.Rectangle((-0.5, i-0.5), cols, 1, lw=2,
                                      edgecolor='green', facecolor='none')
            ax.add_patch(rect)
        for i, val in enumerate(result):
            ax.text(cols+0.3, i, f'{val:.0f}', ha='center', va='center',
                    fontsize=11, color='green', fontweight='bold')
        ax.text(cols+0.3, -0.9, 'sum→', fontsize=10, color='green')

    ax.set_xlim(-1, cols + 0.8)
    ax.set_ylim(rows + 0.8, -1.2)
    ax.set_xticks([]); ax.set_yticks([])
    ax.set_title(title, fontsize=10)

axes[0].imshow(A, cmap='Blues', aspect='auto')
for i in range(3):
    for j in range(4):
        axes[0].text(j, i, f'{A[i,j]:.0f}', ha='center', va='center', fontsize=12)
axes[0].set_xticks([]); axes[0].set_yticks([])
axes[0].set_title('Original matrix\nshape (3,4)', fontsize=10)

draw_matrix_with_highlight(axes[1], A, 0, 'sum(axis=0)\nresult shape (4,)\n[collapse rows]', A.sum(axis=0))
draw_matrix_with_highlight(axes[2], A, 1, 'sum(axis=1)\nresult shape (3,)\n[collapse columns]', A.sum(axis=1))

plt.suptitle('Axis conventions: axis=0 collapses rows, axis=1 collapses columns', fontsize=12)
plt.tight_layout()
plt.show()

# --- Section 5a: Array creation ---

import numpy as np

print("=== Array creation ===")
print(np.zeros((3, 4)))                     # zeros matrix
print(np.ones(5))                            # ones vector
print(np.eye(3))                             # identity matrix
print(np.arange(0, 10, 2))                  # [0,2,4,6,8]
print(np.linspace(0, 1, 5))                 # [0, .25, .5, .75, 1]
print(np.random.default_rng(0).standard_normal((2, 3)))  # random

# Useful constructors
print(np.diag([1, 2, 3]))                   # diagonal matrix from vector
print(np.diag(np.array([[1,2],[3,4]])))      # extract diagonal from matrix
print(np.full((2, 3), 7.0))                 # constant matrix

# --- Section 5b: Indexing and slicing ---

import numpy as np

A = np.arange(12).reshape(3, 4).astype(float)
print("A =", A)
print()

print("A[0]       =", A[0])       # first row (shape (4,))
print("A[0:1]     =", A[0:1])     # first row (shape (1,4)) — keeps dims
print("A[:,1]     =", A[:,1])     # second column (shape (3,))
print("A[1:, 2:]  =", A[1:, 2:]) # submatrix (shape (2,2))
print("A[[0,2], :]=", A[[0,2],:]  ) # rows 0 and 2 (fancy indexing)
print("A[A > 5]   =", A[A > 5])  # boolean mask (shape (N,) flat)

# View vs copy
b = A[0]          # VIEW: modifying b modifies A
c = A[0].copy()   # COPY: independent
b[0] = 999
print(f"\nAfter b[0]=999: A[0,0]={A[0,0]} (view modified A)")
c[0] = 999
print(f"After c[0]=999: A[0,0]={A[0,0]} (copy did NOT modify A)")

# --- Section 5c: Shape manipulation ---

import numpy as np

v = np.arange(12)
print("Original:", v.shape)             # (12,)
print("Reshape (3,4):", v.reshape(3,4).shape)
print("Reshape (2,2,3):", v.reshape(2,2,3).shape)
print("Reshape (4,-1):", v.reshape(4,-1).shape)  # -1 inferred

print()
a = np.array([1., 2., 3.])  # shape (3,)
print("a shape:", a.shape)
print("a[:,None] shape:", a[:,None].shape)  # (3,1) — column vector
print("a[None,:] shape:", a[None,:].shape)  # (1,3) — row vector
print("a.reshape(-1,1) shape:", a.reshape(-1,1).shape)  # (3,1)

# Stack and concatenate
b = np.array([4., 5., 6.])
print("\nnp.stack([a,b], axis=0):", np.stack([a,b], axis=0).shape)  # (2,3)
print("np.stack([a,b], axis=1):", np.stack([a,b], axis=1).shape)  # (3,2)
print("np.concatenate([a,b]):", np.concatenate([a,b]).shape)      # (6,)

# --- Section 5d: Reductions and linear algebra ---

import numpy as np

rng = np.random.default_rng(42)
A = rng.standard_normal((4, 5))

print("=== Reductions ===")
print("sum()     :", A.sum())                    # scalar
print("sum(ax=0) :", A.sum(axis=0).shape)        # (5,) — per column
print("sum(ax=1) :", A.sum(axis=1).shape)        # (4,) — per row
print("mean(ax=0):", A.mean(axis=0).round(3))
print("max()     :", A.max())
print("argmax(ax=1):", A.argmax(axis=1))          # index of max in each row
print("any(A>2)  :", (A>2).any())
print("all(A>-3) :", (A>-3).all())

print("\n=== Linear algebra ===")
# These are the workhorses of all Part VI operations
v = rng.standard_normal(5)
print("A @ v shape       :", (A @ v).shape)       # (4,)
print("norm(v)           :", np.linalg.norm(v).round(4))
print("norm(A, axis=1)   :", np.linalg.norm(A, axis=1).round(3))
print("det(A[:4,:4])     :", np.linalg.det(A[:4,:4]).round(4))
print("rank(A)           :", np.linalg.matrix_rank(A))

U, s, Vt = np.linalg.svd(A, full_matrices=False)
print(f"SVD shapes: U{U.shape}, s{s.shape}, Vt{Vt.shape}")
print("Reconstruction error:", np.linalg.norm(A - U @ np.diag(s) @ Vt).round(10))

# --- Experiment 1: View vs copy pitfall ---
# Hypothesis: Slicing returns a view; fancy indexing returns a copy.

import numpy as np

A = np.arange(9, dtype=float).reshape(3,3)
print("Original A:")
print(A)

# Slice: VIEW
B = A[0:2, 0:2]
B[0,0] = 999
print("\nAfter B=A[0:2,0:2], B[0,0]=999 (slice = view):")
print(A)  # A is modified!

# Reset
A = np.arange(9, dtype=float).reshape(3,3)

# Fancy indexing: COPY
C = A[[0,1], :][:, [0,1]]
C[0,0] = 999
print("\nAfter C=A[[0,1]][:,[0,1]], C[0,0]=999 (fancy = copy):")
print(A)  # A is NOT modified

print()
print("Rule: use .copy() whenever you want an independent array.")

# --- Experiment 2: dtype overflow and precision ---
# Hypothesis: Integer arrays silently overflow; float64 is the safe default.

import numpy as np

# Integer overflow
a_int8 = np.array([127], dtype=np.int8)
print(f"int8 max = 127. Adding 1: {a_int8 + 1}  (overflow, wraps to -128)")

# Float precision
a = np.array([0.1, 0.2], dtype=np.float32)
b = np.array([0.1, 0.2], dtype=np.float64)
print(f"float32: 0.1 + 0.2 = {a.sum():.10f}")
print(f"float64: 0.1 + 0.2 = {b.sum():.10f}")
print(f"Exact:               0.30000000000000")

# Division dtype change
x = np.array([1, 2, 3])  # int64
print(f"\nInteger array: {x.dtype}")
print(f"After /2.0: {(x/2.0).dtype}  values={x/2.0}")
print(f"After //2:  {(x//2).dtype}   values={x//2}  (floor div, stays int)")

# --- Experiment 3: keepdims for broadcasting compatibility ---
# Hypothesis: Without keepdims, reductions drop a dimension and break broadcasting.

import numpy as np

X = np.random.default_rng(0).standard_normal((4, 5))

# Without keepdims
means = X.mean(axis=1)             # shape (4,) — 1D
try:
    result = X - means             # ERROR: (4,5) - (4,) — ambiguous broadcast
    print("Without keepdims: subtracted (unexpected broadcast)")
    print("Result shape:", result.shape)
except:
    print("Without keepdims may broadcast unexpectedly")

# With keepdims
means_kd = X.mean(axis=1, keepdims=True)  # shape (4,1)
result_correct = X - means_kd             # (4,5) - (4,1) -> (4,5) ✓
print(f"With keepdims: means shape={means_kd.shape}, result shape={result_correct.shape}")
print(f"Row means after subtraction: {result_correct.mean(axis=1).round(10)}")
print("(all ~0 — rows are zero-centered)")

# --- Mini Project: NumPy Array Operations Cheat Sheet Validator ---
#
# Build a self-checking reference notebook: implement each operation
# and verify it against a known result or mathematical property.
# This serves as both a learning tool and a regression test suite.

import numpy as np

rng = np.random.default_rng(0)
PASS = "✓"
FAIL = "✗"

def check(name, condition):
    status = PASS if condition else FAIL
    print(f"  [{status}] {name}")

A = rng.standard_normal((4, 5))
v = rng.standard_normal(5)
w = rng.standard_normal(5)

print("=== Linear algebra identities ===")
check("(AB)ᵀ = BᵀAᵀ",
      np.allclose((A[:3] @ A[:3].T).T, A[:3] @ A[:3].T))
check("||Av|| ≤ σ_max ||v|| (spectral norm)",
      np.linalg.norm(A @ v) <= np.linalg.svd(A, compute_uv=False)[0] * np.linalg.norm(v) + 1e-10)
check("SVD reconstruction A = UΣVᵀ",
      np.allclose(A, np.linalg.svd(A, full_matrices=False)[0] @
                  np.diag(np.linalg.svd(A, full_matrices=False)[1]) @
                  np.linalg.svd(A, full_matrices=False)[2]))
check("v · w = w · v (commutativity)",
      np.isclose(v @ w, w @ v))
check("||v||² = v · v",
      np.isclose(np.linalg.norm(v)**2, v @ v))

print("\n=== Broadcasting ===")
means = A.mean(axis=1, keepdims=True)
A_centered = A - means
check("Row centering: row means ≈ 0",
      np.allclose(A_centered.mean(axis=1), 0, atol=1e-12))

norms = np.linalg.norm(A, axis=1, keepdims=True)
A_normalized = A / norms
check("Row normalization: row norms ≈ 1",
      np.allclose(np.linalg.norm(A_normalized, axis=1), 1, atol=1e-12))

print("\n=== Reductions ===")
check("sum(axis=0) + sum(axis=1) consistency",
      np.isclose(A.sum(axis=0).sum(), A.sum(axis=1).sum()))
check("argmax gives correct index",
      A[2, A[2].argmax()] == A[2].max())

print("\n=== Shape manipulation ===")
B = A.reshape(2, 10)
check("reshape preserves total elements", B.size == A.size)
check("reshape is a view (shares memory)", np.shares_memory(A, B))
check("transpose swaps shape", A.T.shape == (5, 4))

print("\nAll checks complete.")