NumPy チートシート

# NumPy: the fundamental package for numerical computing in Python.
# Core data structure: the n-dimensional array (ndarray).

# Key advantages over pure Python lists:
# - Homogeneous, contiguous memory -> fast vectorized ops
# - Broadcasting: operate on arrays of different shapes
# - Rich set of math / linear algebra / FFT / random routines
# - C/Fortran backend -> orders of magnitude faster than Python loops

import numpy as np
print(np.__version__)   # e.g. 1.26.0

# Install NumPy
pip install numpy
# Or with conda
conda install numpy

# The universal convention is to alias as np
import numpy as np

a = np.array([1, 2, 3])
print(a)            # [1 2 3]
print(type(a))      # <class 'numpy.ndarray'>

import numpy as np

a = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.int64)

a.ndim      # 2          -> number of axes (dimensions)
a.shape     # (2, 3)     -> tuple of axis sizes
a.size      # 6          -> total elements = prod(shape)
a.dtype     # int64      -> element type
a.itemsize  # 8          -> bytes per element
a.nbytes    # 48         -> total bytes = size * itemsize
a.T         # shape (3, 2) transposed view
a.flat      # flat iterator over all elements

import numpy as np

# Specifying dtype at creation
a = np.array([1, 2, 3], dtype=np.float32)
b = np.array([1, 2, 3], dtype="float32")   # string form works too

# Integer types: int8 int16 int32 int64 (and unsigned uint*)
# Float types:   float16 float32 float64 (default)
# Complex:       complex64 complex128
# Others:        bool, string_, unicode_, object

# Cast an existing array
c = a.astype(np.int32)

# Default int/float match the platform
print(np.dtype(int).itemsize)        # 4 or 8
print(np.dtype(float).itemsize)      # 8

# Beware overflow:
small = np.array([200], dtype=np.int8)   # wraps around!

import numpy as np
import time

# Python list: element-wise multiply needs a loop or comprehension
py_list = list(range(1_000_000))
start = time.perf_counter()
py_sq = [x * 2 for x in py_list]
print("list:", time.perf_counter() - start)

# NumPy: vectorized, runs in C
arr = np.arange(1_000_000)
start = time.perf_counter()
np_sq = arr * 2
print("numpy:", time.perf_counter() - start)

# NumPy is typically 20-100x faster for numeric workloads.
# Lists are heterogeneous & flexible; arrays are homogeneous & fast.

import numpy as np

# From a list (1-D)
a = np.array([1, 2, 3])

# From a list of lists (2-D)
b = np.array([[1, 2, 3],
              [4, 5, 6]])

# From nested lists (3-D)
c = np.array([[[1, 2], [3, 4]],
              [[5, 6], [7, 8]]])

# Specify dtype up front
d = np.array([1, 2, 3], dtype=np.float32)

# From any iterable
e = np.array(range(5))         # [0 1 2 3 4]
f = np.array((1, 2, 3))        # from a tuple

import numpy as np

# arange: like range but returns array (watch float edge!)
np.arange(0, 10, 2)         # [0 2 4 6 8]
np.arange(0, 1, 0.1)        # may or may not include 1.0 (float!)

# linspace: exactly num points, endpoints included by default
np.linspace(0, 1, 5)        # [0.   0.25 0.5  0.75 1.  ]

# logspace: 10**start .. 10**stop
np.logspace(0, 3, 4)        # [   1.   10.  100. 1000.]

# Meshgrids for evaluating functions on a grid
x = np.linspace(-2, 2, 5)
y = np.linspace(-2, 2, 5)
X, Y = np.meshgrid(x, y)
Z = X**2 + Y**2             # 5x5 grid of values

import numpy as np

np.zeros((2, 3))            # 2x3 of float64 zeros
np.zeros((2, 3), dtype=int) # 2x3 of int zeros
np.ones((3, 3))             # 3x3 of ones
np.full((2, 2), 7)          # 2x2 filled with 7
np.full((2, 2), np.pi)      # fill with any constant

# Identity & eye
np.eye(3)                   # 3x3 identity matrix
np.eye(3, k=1)              # identity shifted up by one diagonal
np.diag([1, 2, 3])          # diagonal matrix from a vector

# Empty: uninitialized memory (fastest, but garbage values)
np.empty((2, 2))            # allocate without filling

import numpy as np

rng = np.random.default_rng(seed=42)   # modern Generator (NumPy >= 1.17)

rng.random((2, 3))           # uniform in [0, 1), shape (2, 3)
rng.standard_normal((2, 3))  # standard normal (mean 0, std 1)
rng.normal(100, 15, size=5)  # normal(mu=100, sigma=15)
rng.uniform(0, 10, size=5)   # uniform in [low, high)
rng.integers(0, 10, size=5)  # random ints in [low, high)

# Reproducibility: always pass a seed
rng2 = np.random.default_rng(42)
print(rng2.random(3))        # same every run

# Legacy API (still common but prefer default_rng)
# np.random.rand(2, 3)      # uniform [0,1)
# np.random.randn(2, 3)     # standard normal

import numpy as np
a = np.array([1, 2, 3])

# View (shares memory!) - shape change only
v = a.view()
v[0] = 99
print(a)   # [99  2  3]   <- a is mutated!

# Deep copy (independent)
b = a.copy()
b[0] = 0
print(a)   # [99  2  3]   <- a unchanged

# np.array defaults to copying
c = np.array(a)
c is a    # False

# Repeat / tile to expand arrays
np.repeat([1, 2, 3], 2)   # [1 1 2 2 3 3]
np.tile([1, 2, 3], 2)     # [1 2 3 1 2 3]

import numpy as np
a = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

a[0]          # 0      (scalar)
a[-1]         # 9      (last element)
a[2:5]        # [2 3 4]    start:stop (stop excluded)
a[:5]         # [0 1 2 3 4]
a[5:]         # [5 6 7 8 9]
a[::2]        # [0 2 4 6 8]   step=2
a[::-1]       # [9 8 7 6 5 4 3 2 1 0]  reversed
a[1:7:2]      # [1 3 5]       start:stop:step

# Slices return VIEWS (not copies) -> editing changes the original
b = a[2:5]
b[0] = 99
print(a[2])   # 99   <- a was modified!

import numpy as np
a = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

a[1, 2]       # 6      (row 1, col 2) - equivalent to a[1][2] but faster
a[0]          # [1 2 3]   (row 0)
a[0, :]       # [1 2 3]   (row 0, all cols)
a[:, 1]       # [2 5 8]   (all rows, col 1)
a[1:3, 0:2]   # [[4 5]
              #  [7 8]]     rows 1-2, cols 0-1
a[-1]         # [7 8 9]    last row
a[:, -1]      # [3 6 9]    last column

# ... (ellipsis) selects any remaining full axes
b = np.zeros((2, 3, 4))
b[..., 0]     # shape (2, 3) - last axis indexed, others kept
b[1, ...]     # shape (3, 4) - first axis indexed

import numpy as np
a = np.array([1, 2, 3, 4, 5, 6])

# Build a boolean mask
mask = a > 3
# [False False False  True  True  True]

a[mask]            # [4 5 6]   (select elements where mask is True)
a[a > 3]           # same, in one expression
a[(a > 2) & (a < 6)]   # [3 4 5]   combine with & (and), | (or), ~ (not)

# Modify via mask
a[a > 3] = 0       # set all matching elements
a[a == 0] = -1     # conditional assignment

# Count matching elements
(a > 3).sum()      # 3   (True == 1)
np.count_nonzero(a > 3)

# np.where: positions where condition is True
np.where(a > 3)    # (array([3, 4, 5]),) indices

import numpy as np
a = np.array([10, 20, 30, 40, 50])

# Index with a list/array of indices
a[[0, 2, 4]]         # [10 30 50]
a[np.array([0, 2, 4])]   # same

# Indices can repeat
a[[0, 0, 1]]         # [10 10 20]

# Negative indices work
a[[-1, -2]]          # [50 40]

# 2-D fancy indexing
b = np.arange(12).reshape(3, 4)
rows = np.array([0, 2])
cols = np.array([1, 3])
b[rows, cols]        # [1 11]   pairs (0,1) and (2,3)

# Mix a slice and an index array
b[1:, [0, 2]]        # rows 1-2, cols 0 and 2

# np.take: same as fancy indexing, more flexible
np.take(a, [0, 2])   # [10 30]

import numpy as np
a = np.array([1, 2, 3, 4, 5, 6])

# Conditional assignment with a mask
a[a % 2 == 0] = -1     # even -> -1: [1 -1 3 -1 5 -1]

# np.where(cond, x, y): vectorized ternary
b = np.where(a > 3, "big", "small")
# ['small' 'small' 'small' 'big' 'big' 'big']

# np.where with one argument -> indices of True
idx = np.where(a == -1)
# (array([1, 3, 5]),)

# np.select: multiple conditions, multiple choices
conds = [a < 3, a > 3, a == 3]
choices = ["low", "high", "mid"]
np.select(conds, choices, default="ok")

# np.clip: clamp values to a range
np.clip(a, 0, 4)       # values < 0 -> 0, > 4 -> 4

import numpy as np
a = np.arange(12)        # [0 1 2 ... 11], shape (12,)

# reshape to (3, 4)
b = a.reshape(3, 4)
# [[ 0  1  2  3]
#  [ 4  5  6  7]
#  [ 8  9 10 11]]

# Use -1 to infer one dimension
a.reshape(2, -1)         # (2, 6)  - -1 means "figure it out"
a.reshape(-1, 3)         # (4, 3)
a.reshape(2, -1, 2)      # (2, 3, 2)  only one -1 allowed

# Total size must match -> ValueError otherwise
# a.reshape(5, 3)        # 5*3 = 15 != 12

# reshape returns a view when possible (no copy)
b[0, 0] = 99
print(a[0])              # 99   <- a was modified

import numpy as np
a = np.array([[1, 2, 3],
              [4, 5, 6]])

# ravel: returns a VIEW (C order by default) - efficient
a.ravel()              # [1 2 3 4 5 6]
a.ravel(order='F')     # [1 4 2 5 3 6]   Fortran (column) order

# flatten: always returns a COPY - safe to mutate
flat = a.flatten()
flat[0] = 99
print(a[0, 0])         # 1   <- a unchanged

# Order conventions:
# 'C' (default): last axis varies fastest (row-major)
# 'F':           first axis varies fastest (column-major)
# 'A':           Fortran order if array is F-contiguous, else C
# 'K':           preserve memory order

import numpy as np
a = np.arange(12).reshape(3, 4)

a.T                 # shape (4, 3)   transpose (2-D)
a.transpose()       # same

# For 3+ dimensions, specify axis order
b = np.arange(24).reshape(2, 3, 4)
b.transpose(1, 0, 2)   # shape (3, 2, 4)  axes reordered
b.swapaxes(0, 1)       # shape (3, 2, 4)  swap two axes
b.T                   # shape (4, 3, 2)  reverse all axes

# Move axes to new positions
np.moveaxis(b, 0, -1)  # shape (3, 4, 2)

# All of these return VIEWS (no data copy)

import numpy as np
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

# Stack along a NEW axis
np.stack([a, b])            # shape (2, 3)  - default axis=0
np.stack([a, b], axis=1)    # shape (3, 2)  column stack

# Concatenate along an EXISTING axis
np.concatenate([a, b])              # [1 2 3 4 5 6]
np.vstack([a, b])                  # shape (2, 3)  vertical
np.hstack([a, b])                  # [1 2 3 4 5 6] horizontal
np.dstack([a, b])                  # shape (1, 3, 2) depth

# 2-D arrays
A = np.ones((2, 2))
B = np.zeros((2, 2))
np.hstack([A, B])   # shape (2, 4)
np.vstack([A, B])   # shape (4, 2)

# np.block: build from nested blocks
np.block([[A, B], [B, A]])   # shape (4, 4)

import numpy as np
a = np.arange(12)

# Split into equal parts (must divide evenly)
np.split(a, 3)        # [array([0,1,2,3]), array([4,5,6,7]), array([8,9,10,11])]

# Split at specific indices
np.split(a, [3, 5])   # [array([0,1,2]), array([3,4]), array([5,6,7,8,9,10,11])]

# 2-D
b = np.arange(16).reshape(4, 4)
np.vsplit(b, 2)       # split along rows -> 2 arrays of shape (2, 4)
np.hsplit(b, 2)       # split along cols -> 2 arrays of shape (4, 2)

# array_split: allows uneven division
np.array_split(np.arange(10), 3)  # sizes [4, 3, 3] - no error

# dsplit: split along depth (axis 2)
c = np.arange(24).reshape(2, 2, 6)
np.dsplit(c, 3)       # 3 arrays of shape (2, 2, 2)

import numpy as np

# Broadcasting: operate on arrays of different shapes without copying.
# Two shapes are compatible if, for each trailing axis, they are
# equal OR one of them is 1.

# Scalar broadcast to any array
a = np.array([1, 2, 3])
a + 10            # [11 12 13]   10 treated as shape ()

# 1-D with 2-D
A = np.ones((3, 4))
row = np.array([10, 20, 30, 40])
A + row           # row (4,) broadcast to (3, 4) -> added to each row

col = np.array([[100], [200], [300]])   # shape (3, 1)
A + col           # col broadcast across columns -> added to each column

# Incompatible shapes raise ValueError
# A + np.array([1, 2, 3])   # (3,4) vs (3,) -> error

import numpy as np
a = np.array([1, 2, 3])          # shape (3,)
b = np.array([10, 20, 30, 40])   # shape (4,)

# a + b would fail (3,) vs (4,)
# Insert a new axis so shapes become (3, 1) and (4,) -> broadcasts to (3, 4)

a_col = a[:, np.newaxis]   # shape (3, 1)  - same as a.reshape(-1, 1)
b_row = b[np.newaxis, :]   # shape (1, 4)

outer = a_col + b_row      # shape (3, 4)  outer sum
# [[11 21 31 41]
#  [12 22 32 42]
#  [13 23 33 43]]

# np.newaxis / None are equivalent
a[:, None]                 # same as a[:, np.newaxis]

# np.expand_dims is the explicit function form
np.expand_dims(a, axis=1)  # shape (3, 1)

import numpy as np
a = np.array([1, 2, 3])
b = np.array([10, 20, 30, 40])

# Outer product (3,) x (4,) -> (3, 4)
np.outer(a, b)
# [[10 20 30 40]
#  [20 40 60 80]
#  [30 60 90 120]]

# Manual outer via broadcasting
a[:, None] * b[None, :]      # same as np.outer(a, b)

# Pairwise difference: (4,) - (3,) -> (4, 3) or (3, 4)?
diff = b[:, None] - a[None, :]   # shape (4, 3): diff[i, j] = b[i] - a[j]

# Pairwise equality (used in set ops, distance)
x = np.array([1, 2, 3, 4])
y = np.array([2, 4, 6])
x[:, None] == y[None, :]     # shape (4, 3) boolean matrix

# Add a constant vector to each row of a matrix
M = np.zeros((5, 3))
v = np.array([1, 2, 3])
M + v                        # broadcasts v to every row

import numpy as np

# Pitfall 1: (N,) is NOT a row or column vector - it has no axis direction
a = np.array([1, 2, 3])
M = np.ones((3, 3))
M + a            # OK: (3,3) + (3,) -> a broadcast as a ROW
M + a[:, None]   # OK: a as a COLUMN

# Pitfall 2: accidentally broadcasting when you wanted element-wise
x = np.array([1, 2, 3])
y = np.array([10, 20, 30])
x + y            # [11 22 33]   element-wise (both (3,))
x[:, None] + y   # (3,3) matrix  - probably not what you wanted!

# Pitfall 3: shape (1,) vs () - both broadcast but mean different things
np.array([5]).shape   # (1,)
np.array(5).shape     # ()  - 0-D, broadcasts to ANYTHING

# Pitfall 4: memory blow-up from over-broadcasting
# Big arrays broadcast together can produce enormous intermediates
A = np.ones((1000, 1000))
B = np.ones((1000, 1000))
(A * B).sum()    # OK - intermediate is 8MB
# A[:, None, :] * B[None, :, :]   # (1000,1000,1000) = 8GB! crash

import numpy as np

# np.broadcast_to: explicit view of an array broadcast to a new shape
a = np.array([1, 2, 3])            # shape (3,)
b = np.broadcast_to(a, (4, 3))     # shape (4, 3) - VIEW, no copy
# [[1 2 3]
#  [1 2 3]
#  [1 2 3]
#  [1 2 3]]
b[0, 0] = 99   # ERROR: broadcast_to returns read-only view

# np.broadcast_arrays: broadcast several arrays together (as views)
x = np.array([1, 2, 3])           # (3,)
y = np.array([[10], [20]])        # (2, 1)
bx, by = np.broadcast_arrays(x, y)
bx.shape   # (2, 3)
by.shape   # (2, 3)

# np.broadcast: low-level iterator object
b = np.broadcast(x, y)
b.shape    # (2, 3)
b.index    # current position in the iteration
b.iters    # list of iterators over each input

import numpy as np
a = np.array([1, 2, 3, 4])
b = np.array([10, 20, 30, 40])

a + b       # [11 22 33 44]
a - b       # [-9 -18 -27 -36]
a * b       # [10 40 90 160]   element-wise (NOT matrix product)
a / b       # [0.1 0.1 0.1 0.1]  always float
a // b      # [0 0 0 0]        floor division
a ** 2      # [1 4 9 16]       power
a % 2       # [1 0 1 0]        modulo
-a          # [-1 -2 -3 -4]    negation

# In-place ops (modify a, save memory)
a += 10
a *= 2

# Scalar ops broadcast
a + 100     # [102 104 106 108]

# Universal functions (ufuncs) - same operations as function form
np.add(a, b)
np.multiply(a, b)
np.power(a, 2)

import numpy as np
a = np.array([-1.5, 0.0, 1.5, 2.5])

# Math ufuncs
np.abs(a)         # [1.5 0.  1.5 2.5]
np.sqrt(np.abs(a))
np.exp(a)         # e^a
np.log(np.abs(a) + 1)
np.log2(a + 5)
np.log10(a + 5)

# Rounding
np.round(a)       # [-2.  0.  2.  2.]
np.floor(a)       # [-2.  0.  1.  2.]
np.ceil(a)        # [-1.  0.  2.  3.]
np.trunc(a)       # [-1.  0.  1.  2.]  truncate toward zero

# Trigonometry (radians)
np.sin(np.pi/2)   # 1.0
np.cos(0)         # 1.0
np.degrees(np.pi) # 180.0
np.radians(180)   # 3.14159...

# Comparison ufuncs
np.greater(a, 0)  # [False False True True]
np.maximum(a, 0)  # clip negatives to 0  (element-wise max)
np.minimum(a, 0)

import numpy as np
a = np.array([[1, 2, 3],
              [4, 5, 6]])

# Global reductions
a.sum()         # 21
a.prod()        # 720
a.min(), a.max()
a.mean()        # 3.5
a.std(), a.var()
a.argmin(), a.argmax()   # indices of min/max (flattened)

# Axis-wise reductions
a.sum(axis=0)   # [5 7 9]   sum down each column
a.sum(axis=1)   # [6 15]    sum across each row
a.min(axis=1)   # [1 4]     min of each row
a.cumsum(axis=0)   # cumulative sum down columns
a.cumprod(axis=1)  # cumulative product across rows

# Keep dimensions
a.sum(axis=0, keepdims=True)   # shape (1, 3) - keeps the axis

# Specify dtype for accumulation (avoid overflow!)
a.sum(dtype=np.float64)

# NaN-aware versions
b = np.array([1.0, np.nan, 3.0])
np.nansum(b)    # 4.0   treats NaN as 0

import numpy as np
a = np.array([1, 2, 3, 4, 5])

# Cumulative sum and product
np.cumsum(a)         # [ 1  3  6 10 15]
np.cumprod(a)        # [  1   2   6  24 120]

# 2-D: specify axis
b = np.array([[1, 2], [3, 4]])
np.cumsum(b, axis=0)   # [[1 2]
                       #  [4 6]]
np.cumsum(b, axis=1)   # [[1 3]
                       #  [3 7]]

# Differences (inverse of cumulative sum)
np.diff(a)            # [1 1 1 1]   a[i+1] - a[i]
np.diff(a, n=2)       # second-order diff [0 0 0]

# Edge: cumsum on integers can overflow silently
big = np.array([2**30, 2**30, 2**30], dtype=np.int32)
np.cumsum(big)        # wraps around - use int64 or float!

# Cumulative max / min
np.maximum.accumulate(a)   # [1 2 3 4 5]
np.minimum.accumulate(a)   # [1 1 1 1 1]

import numpy as np
a = np.array([1, 2, 3, 4, 5])

# Element-wise comparisons return boolean arrays
a > 3            # [False False False  True  True]
a == 3           # [False False  True False False]
a != 3
np.greater(a, 3)

# Logical ops on boolean arrays
b = (a > 2)
c = (a < 5)
np.logical_and(b, c)   # [False False  True  True False]
np.logical_or(b, c)
np.logical_not(b)
np.logical_xor(b, c)

# Array-level truth testing
np.any(a > 4)    # True
np.all(a > 0)    # True
np.all(a > 0, axis=0)  # per-axis
(a == 3).any()   # True

# np.isclose for float comparisons
np.isclose(0.1 + 0.2, 0.3)        # True  (handles float error)
np.allclose([0.1, 0.2], [0.1, 0.2 + 1e-9], atol=1e-8)  # True

import numpy as np
a = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

a.mean()          # 5.5    arithmetic mean
a.median() if hasattr(a, 'median') else np.median(a)  # 5.5
a.var()           # 8.25   variance (population, ddof=0)
a.std()           # 2.87   standard deviation
a.min(), a.max()  # 1, 10
a.ptp()           # 9      peak-to-peak = max - min

# Sample vs population variance
a.var(ddof=1)     # sample variance (divide by N-1, like pandas default)

# Percentiles & quantiles
np.percentile(a, 50)    # 5.5   median
np.percentile(a, [25, 50, 75])  # [3.25 5.5  7.75]  quartiles
np.quantile(a, 0.5)     # 5.5   same as percentile / 100

# Correlation
x = np.array([1, 2, 3, 4])
y = np.array([2, 4, 5, 4])
np.corrcoef(x, y)       # 2x2 correlation matrix

import numpy as np
data = np.random.default_rng(42).standard_normal(1000)

# 1-D histogram
counts, edges = np.histogram(data, bins=10)
# counts: array of element counts per bin
# edges: array of bin edge positions (length = bins + 1)

# Custom bin edges
np.histogram(data, bins=[-3, -1, 0, 1, 3])

# Density (integrates to 1) vs counts
np.histogram(data, bins=20, density=True)

# 2-D histogram (for heatmaps)
x = np.random.default_rng(1).standard_normal(1000)
y = np.random.default_rng(2).standard_normal(1000)
H, xedges, yedges = np.histogram2d(x, y, bins=20)

# Digitize: which bin does each value fall into?
bins = [-2, -1, 0, 1, 2]
indices = np.digitize(data, bins)   # array of bin indices (1-based)

# Bincount: count occurrences of each non-negative int
np.bincount(np.array([0, 1, 1, 3, 3, 3]))   # [1 2 0 3]

import numpy as np
a = np.array([3, 1, 4, 1, 5, 9, 2, 6])

# Sorting
np.sort(a)             # [1 1 2 3 4 5 6 9]   returns a copy
a.sort()               # in-place sort (modifies a)
np.sort(a)[::-1]       # descending
np.argsort(a)          # indices that would sort a (no copy)

# 2-D sorting along axis
b = np.array([[3, 1, 2], [6, 5, 4]])
np.sort(b, axis=0)     # sort down each column
np.sort(b, axis=1)     # sort across each row (default)

# Partial sorting: top-k smallest (no full sort)
np.partition(a, 3)     # 3 smallest first (unordered), rest after
np.argpartition(a, 3)  # indices of partition

# k-th smallest element (selection algorithm)
np.partition(a, 2)[2]  # 3rd smallest (index 2)

# unique values (sorted)
np.unique([1, 2, 1, 3, 2])   # [1 2 3]
np.unique(a, return_counts=True)  # values and counts

import numpy as np
# Sample data: 3 variables, 5 observations
data = np.array([[1, 2, 3, 4, 5],
                 [2, 4, 1, 3, 5],
                 [5, 4, 3, 2, 1]], dtype=float)

# Covariance matrix (3x3)
# Each ROW is a variable, so set rowvar=False to treat COLUMNS as variables
cov = np.cov(data, rowvar=True)   # 3x3 covariance between rows
# Diagonal = variance of each variable
# Off-diagonal = covariance between pairs

# Correlation matrix
corr = np.corrcoef(data)   # 3x3, values in [-1, 1]

# Correlation between two specific variables
np.corrcoef(data[0], data[1])[0, 1]   # scalar

# Weighted covariance: compute manually
weights = np.array([1, 1, 1, 2, 2])
mean = np.average(data, axis=1, weights=weights, keepdims=True)
centered = data - mean
weighted_cov = (centered * weights) @ centered.T / weights.sum()

# Pearson correlation coefficient (manual)
def pearson(x, y):
    xm, ym = x - x.mean(), y - y.mean()
    return np.dot(xm, ym) / (np.sqrt((xm**2).sum()) * np.sqrt((ym**2).sum()))

import numpy as np
a = np.array([1.0, 2.0, np.nan, 4.0, 5.0])

# Standard functions propagate NaN
a.mean()      # nan
a.sum()       # nan

# NaN-aware equivalents
np.nanmean(a)     # 3.0   ignores NaN
np.nansum(a)      # 12.0
np.nanstd(a)
np.nanmedian(a)
np.nanvar(a)
np.nanmin(a), np.nanmax(a)
np.nanpercentile(a, 50)

# Count non-NaN elements
np.count_nonzero(~np.isnan(a))   # 4

# Find NaN positions
np.isnan(a)            # [False False  True False False]
np.where(np.isnan(a))  # (array([2]),)

# Replace NaN with a value (imputation)
np.nan_to_num(a, nan=0.0)              # replace NaN with 0
np.where(np.isnan(a), 0, a)            # same, more explicit
a_filled = np.copy(a)
a_filled[np.isnan(a)] = np.nanmean(a)  # fill with mean

# Check if any NaN
np.isnan(a).any()      # True

import numpy as np

A = np.array([[1, 2],
              [3, 4]])
B = np.array([[5, 6],
              [7, 8]])

# Matrix product (2-D)
A @ B                  # [[19 22]
                       #  [43 50]]    preferred operator (Py 3.5+)
A.dot(B)               # same, older method
np.dot(A, B)           # same for 2-D

# 1-D: dot is inner product
v = np.array([1, 2, 3])
w = np.array([4, 5, 6])
v @ w                  # 32   (1*4 + 2*5 + 3*6)
np.dot(v, w)           # 32

# Mixed 1-D and 2-D
v = np.array([1, 2])
A @ v                  # [5 11]   matrix-vector product

# Batched matmul (3-D)
X = np.random.rand(10, 3, 4)   # 10 matrices of shape (3, 4)
Y = np.random.rand(10, 4, 2)
(X @ Y).shape          # (10, 3, 2)

# Element-wise product (Hadamard) - NOT matrix product
A * B                  # [[5 12]
                       #  [21 32]]

import numpy as np
A = np.array([[4, 3],
              [6, 3]])

# SVD: A = U @ diag(s) @ Vt
U, s, Vt = np.linalg.svd(A)
# U: (2,2) left singular vectors
# s: (2,) singular values (sorted descending)
# Vt: (2,2) right singular vectors (already transposed)

# Reconstruct A from SVD
A_reconstructed = U @ np.diag(s) @ Vt
np.allclose(A, A_reconstructed)   # True

# QR decomposition: A = Q @ R
Q, R = np.linalg.qr(A)
# Q: orthonormal columns
# R: upper triangular

# Cholesky (for symmetric positive-definite matrices)
S = A @ A.T   # symmetric positive-definite
L = np.linalg.cholesky(S)
np.allclose(S, L @ L.T)   # True

# Eigen-decomposition (square matrices)
w, V = np.linalg.eig(A)
# w: eigenvalues
# V: eigenvectors as columns

# Eigenvalues only (faster)
np.linalg.eigvals(A)

import numpy as np

A = np.array([[3, 2],
              [1, 4]])
b = np.array([7, 6])

# Solve Ax = b (preferred - no explicit inverse)
x = np.linalg.solve(A, b)        # [1.6 1.1]
np.allclose(A @ x, b)            # True

# Multiple right-hand sides
B = np.array([[7, 1], [6, 2]])   # 2x2
X = np.linalg.solve(A, B)        # 2x2 solution matrix

# Inverse (avoid if you only need to solve Ax=b)
A_inv = np.linalg.inv(A)
A_inv @ b    # same result as solve, but slower & less stable

# Pseudo-inverse (works for non-square / singular matrices)
A_pinv = np.linalg.pinv(A)       # Moore-Penrose pseudo-inverse

# Determinant
np.linalg.det(A)    # 10.0

# Matrix rank
np.linalg.matrix_rank(A)   # 2

# Condition number (large = ill-conditioned)
np.linalg.cond(A)   # ~1.9

import numpy as np
v = np.array([3, 4])

# Vector norms
np.linalg.norm(v)         # 5.0    L2 (Euclidean) - default
np.linalg.norm(v, ord=1)  # 7.0    L1 (sum of |x|)
np.linalg.norm(v, ord=np.inf)  # 4.0  L-infinity (max |x|)

# Matrix norms
A = np.array([[1, 2], [3, 4]])
np.linalg.norm(A)         # Frobenius norm (default for matrices)
np.linalg.norm(A, ord='fro')
np.linalg.norm(A, ord=2)  # spectral norm (largest singular value)
np.linalg.norm(A, ord=np.inf)  # max row sum

# Normalizing a vector (unit vector)
v / np.linalg.norm(v)

# Distance between two vectors
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
np.linalg.norm(a - b)         # Euclidean distance
np.linalg.norm(a - b, ord=1)  # Manhattan distance

# Outer product, cross product
np.cross([1, 0, 0], [0, 1, 0])   # [0 0 1]   (3-D only)
np.outer([1, 2], [10, 20, 30])   # (2, 3)

import numpy as np

# Linear least squares: solve Ax = b for over-determined A (more rows than cols)
# Minimizes ||Ax - b||^2
A = np.array([[1, 1],
              [1, 2],
              [1, 3],
              [1, 4]])         # design matrix (4x2)
b = np.array([6, 5, 7, 10])    # observations

x, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)
# x: best-fit parameters [intercept, slope] = [3.5, 1.4]
# residuals: sum of squared errors
# rank: effective rank of A
# sv: singular values of A

# Polynomial fitting (least squares under the hood)
x_data = np.array([0, 1, 2, 3, 4])
y_data = np.array([1, 3, 7, 13, 21])
coeffs = np.polyfit(x_data, y_data, deg=2)   # quadratic fit
# coeffs[0]*x^2 + coeffs[1]*x + coeffs[2]

# Evaluate polynomial
poly = np.poly1d(coeffs)
poly(2.5)   # predict y at x=2.5

# Polynomial roots (zeros)
np.roots([1, -3, 2])   # roots of x^2 - 3x + 2 = 0 -> [2, 1]

import numpy as np

# Create a Generator with a seed (modern API, NumPy >= 1.17)
rng = np.random.default_rng(seed=42)

# Uniform [0, 1) floats
rng.random()              # scalar
rng.random(size=5)        # 1-D array of 5
rng.random((3, 4))        # 3x4 array

# Standard normal (mean 0, std 1)
rng.standard_normal(5)
rng.standard_normal((2, 3))

# Integers in [low, high)
rng.integers(0, 10, size=5)       # 5 random ints in [0, 10)
rng.integers(1, 7, size=(2, 3), endpoint=True)  # include 7

# Choice: sample from an array
rng.choice([10, 20, 30, 40], size=5)
rng.choice([10, 20, 30, 40], size=5, replace=True)   # default
rng.choice(range(100), size=10, replace=False)        # without replacement

# Why default_rng? Faster, better statistical properties, no global state.
# Legacy np.random.rand/randn/seed still works but is discouraged.

import numpy as np
rng = np.random.default_rng(42)

# Continuous distributions
rng.uniform(0, 10, size=5)        # uniform in [low, high)
rng.normal(loc=100, scale=15, size=5)   # normal(mu, sigma)
rng.standard_normal(size=5)       # standard normal
rng.exponential(scale=2.0, size=5)   # exponential
rng.gamma(shape=2.0, scale=1.5, size=5)
rng.beta(a=2, b=5, size=5)
rng.lognormal(mean=0.0, sigma=1.0, size=5)
rng.chisquare(df=3, size=5)
rng.standard_t(df=10, size=5)

# Discrete distributions
rng.binomial(n=10, p=0.5, size=5)
rng.poisson(lam=3.0, size=5)
rng.geometric(p=0.1, size=5)
rng.negative_binomial(n=5, p=0.5, size=5)

# Multi-variate normal
mean = [0, 0]
cov = [[1, 0.5], [0.5, 1]]
rng.multivariate_normal(mean, cov, size=1000)   # shape (1000, 2)

# Shuffle & permutation (in-place / copy)
arr = np.arange(10)
rng.shuffle(arr)        # in-place shuffle
perm = rng.permutation(10)   # returns a new array

import numpy as np

# Same seed -> same sequence of random numbers
r1 = np.random.default_rng(42)
r2 = np.random.default_rng(42)
print(r1.random(3))    # [0.7739 0.4389 0.8586]
print(r2.random(3))    # same exact values

# Different seeds -> independent sequences
r3 = np.random.default_rng(123)
print(r3.random(3))    # different values

# Seed from entropy (non-reproducible) - default if seed omitted
r4 = np.random.default_rng()    # seeded from OS entropy

# Spawn independent child generators (for parallelism)
parent = np.default_rng(42)
children = parent.spawn(4)   # 4 independent Generators
# Each child produces an independent, reproducible stream

# SeedSequence for advanced control
ss = np.random.SeedSequence(42)
child_seeds = ss.spawn(4)
streams = [np.random.default_rng(s) for s in child_seeds]

# Reproducibility across versions: pin numpy version too!

import numpy as np
rng = np.random.default_rng(42)

# Shuffle (in-place)
a = np.arange(10)
rng.shuffle(a)         # a is now shuffled

# Permutation (returns a new array, original untouched)
b = np.arange(10)
perm = rng.permutation(b)
# b is unchanged, perm is the shuffled version

# Permutation of indices (common pattern for shuffling a dataset)
idx = rng.permutation(100)
X_shuffled = X[idx]
y_shuffled = y[idx]

# Random choice
rng.choice(100, size=10, replace=False)   # 10 unique indices
rng.choice([1, 2, 3], size=10, p=[0.1, 0.3, 0.6])  # weighted

# Sample rows of a 2-D array
data = np.arange(20).reshape(10, 2)
row_idx = rng.choice(10, size=3, replace=False)
sample = data[row_idx]    # 3 random rows

# Bootstrap sampling (with replacement)
boot_idx = rng.integers(0, len(data), size=len(data))
boot_sample = data[boot_idx]

import numpy as np

# Legacy global-state API (still common in old code)
np.random.seed(42)
a = np.random.rand(3)      # [0.3745 0.9507 0.7320]

np.random.seed(42)
b = np.random.rand(3)      # same as a - global state was reset

# Functions operate on the global RandomState
np.random.rand(2, 3)       # uniform [0, 1)
np.random.randn(2, 3)      # standard normal
np.random.randint(0, 10, 5)
np.random.choice([1, 2, 3], 2)
np.random.shuffle(arr)

# Create an independent RandomState (better than global)
rs = np.random.RandomState(42)
rs.rand(3)                 # independent of np.random's global state

# Migration: replace np.random.X(args) with rng.X(args)
# rng = np.random.default_rng(seed)
# rand()       -> random()
# randn()      -> standard_normal()
# randint(a,b) -> integers(a, b)   # note: endpoint differs
# seed(s)      -> default_rng(s)   # explicit Generator

import numpy as np
a = np.arange(10, 20)   # [10 11 12 13 14 15 16 17 18 19]

# Index with integer arrays
indices = np.array([0, 3, 5, 8])
a[indices]              # [10 13 15 18]

# Negative indices
a[[-1, -2, -3]]         # [19 18 17]

# Index arrays can repeat (creates duplicates)
a[[0, 0, 1, 1]]         # [10 10 11 11]

# Assign through fancy indexing
a[[0, 3, 5]] = -1       # set those positions to -1
a[[1, 2]] = [100, 200]  # different value per position

# 2-D: pair (row, col) indices element-wise
b = np.arange(12).reshape(3, 4)
rows = np.array([0, 1, 2])
cols = np.array([1, 2, 3])
b[rows, cols]           # [1 6 11]  points (0,1), (1,2), (2,3)

# For a cross-product (rectangle) use np.ix_
b[np.ix_([0, 2], [1, 3])]
# [[ 1  3]
#  [ 9 11]]

import numpy as np
a = np.array([5, 10, 15, 20, 25, 30])

# Boolean mask
mask = a > 15
a[mask]                  # [20 25 30]
a[a > 15]                # equivalent

# Multiple conditions: use & | ~ (NOT and/or)
a[(a > 10) & (a < 30)]   # [15 20 25]
a[(a < 10) | (a > 25)]   # [5 30]
a[~(a > 15)]             # [5 10 15]  negation

# Count True
mask.sum()               # 3   (True == 1)

# Boolean indexing returns a COPY (not a view)
b = a[a > 15]
b[0] = 999
print(a[3])              # 20  - a unchanged

# 2-D boolean mask (must match array shape)
M = np.arange(12).reshape(3, 4)
mask2d = M > 5
M[mask2d]                # 1-D array of matching values (flattened)

# Row selection with a 1-D boolean mask
row_mask = np.array([True, False, True])
M[row_mask]              # rows 0 and 2

import numpy as np
a = np.array([1, 2, 3, 4, 5])

# np.where(cond) -> indices where cond is True
np.where(a > 3)             # (array([3, 4]),)   tuple of index arrays

# np.where(cond, x, y) -> vectorized ternary
np.where(a > 3, "big", "small")
# ['small' 'small' 'small' 'big' 'big']

# np.where with arrays for x and y
np.where(a % 2 == 0, a * 10, a)   # double the evens, keep odds

# Multi-dim where
M = np.array([[1, 2], [3, 4]])
np.where(M > 2)             # (array([1, 1]), array([0, 1]))  (row, col) pairs

# np.choose: pick from multiple arrays based on integer index
choices = [[10, 20, 30], [100, 200, 300], [1000, 2000, 3000]]
np.choose([0, 2, 1], choices)   # [10 2000 300]

# np.select: multiple conditions
conds = [a < 3, a > 4, (a >= 3) & (a <= 4)]
choices = ["low", "high", "mid"]
np.select(conds, choices)   # ['low' 'mid' 'mid' 'mid' 'high']

import numpy as np
a = np.arange(20).reshape(4, 5)

# Select rows 0 and 2, columns 1 and 3 (cross product)
rows = [0, 2]
cols = [1, 3]
a[np.ix_(rows, cols)]
# [[ 1  3]
#  [11 13]]

# Without np.ix_, fancy indexing pairs the arrays element-wise:
a[np.array(rows), np.array(cols)]
# ERROR (or wrong result): pairs (0,1) and (2,3) - points, not rectangle

# Boolean masks also work with np.ix_
row_mask = np.array([True, False, True, False])
col_mask = np.array([False, True, False, True, False])
a[np.ix_(row_mask, col_mask)]   # rows where True, cols where True

# 3-D indexing with ix_
b = np.arange(24).reshape(2, 3, 4)
b[np.ix_([0], [0, 2], [1, 3])].shape   # (1, 2, 2)

# np.ix_ converts 1-D index arrays into broadcastable open grids
# so the result is a (Cartesian) rectangle, not points

import numpy as np
a = np.arange(12).reshape(3, 4)

# Most slicing returns VIEWS (no data copy)
v = a[1:, 1:]   # view
v[0, 0] = 99
print(a[1, 1])  # 99   - a was modified

# Fancy and boolean indexing always return COPIES
b = a[[0, 1, 2]]
b[0] = -1
print(a[0, 0])  # 0   - a unchanged

# Check if two arrays share memory
np.shares_memory(a, v)   # True
np.shares_memory(a, b)   # False

# Force a copy with .copy()
independent = a[1:, 1:].copy()

# Memory layout flags
a.flags['C_CONTIGUOUS']    # True  - row-major (C order)
a.flags['F_CONTIGUOUS']    # False - column-major (Fortran)

# Convert to a guaranteed-contiguous array (sometimes speeds up ops)
a_contig = np.ascontiguousarray(a)
a_fortran = np.asfortranarray(a)

# may_share_memory: weaker, faster check
np.may_share_memory(a, v)  # True

import numpy as np
a = np.arange(12).reshape(3, 4)
b = np.array([1, 2, 3])

# Save a single array (binary .npy format)
np.save("a.npy", a)
loaded = np.load("a.npy")
np.array_equal(a, loaded)   # True

# Save multiple arrays (zipped .npz)
np.savez("data.npz", a=a, b=b, c=np.zeros(3))
data = np.load("data.npz")
print(data["a"])   # the saved array a
print(data["b"])   # array b
print(list(data.keys()))   # ['a', 'b', 'c']

# Compressed savez
np.savez_compressed("data.npz", a=a, b=b)   # smaller file, slower

# Loading returns a lazy NpzFile - arrays load on access
# Close it explicitly if you have many files open
with np.load("data.npz") as data:
    a = data["a"]
    b = data["b"]
# file handle auto-closed after the with block

import numpy as np

# Save as plain text (1-D or 2-D only)
a = np.arange(12).reshape(3, 4)
np.savetxt("a.csv", a, delimiter=",", fmt="%d",
           header="col1,col2,col3,col4", comments="")

# Load text
loaded = np.loadtxt("a.csv", delimiter=",", skiprows=1, dtype=int)

# genfromtxt: handles missing values and irregular data
data = np.genfromtxt("data.csv", delimiter=",",
                     names=True,          # first row = field names
                     dtype=None,          # infer types per column
                     missing_values="",
                     filling_values=0)
# data is a structured array; access columns by name
print(data["age"])

# More flexible: pandas.read_csv, then convert to NumPy
import pandas as pd
df = pd.read_csv("data.csv")
arr = df.to_numpy()      # 2-D ndarray
ages = df["age"].to_numpy()   # 1-D

import numpy as np

# Create a large array on disk (without loading it all into RAM)
shape = (10000, 10000)
dtype = np.float64
mm = np.memmap("big.dat", dtype=dtype, mode="w+", shape=shape)

# Write to it as if it were a normal array
mm[:] = np.arange(shape[0])[:, None] + np.arange(shape[1])[None, :]
mm.flush()   # flush changes to disk

# Re-open for reading without loading the whole file
mm_read = np.memmap("big.dat", dtype=dtype, mode="r", shape=shape)
print(mm_read[0, 0])    # access a small slice -> only that page loads
print(mm_read[-1, -1])

# Mode: 'r' (read-only), 'r+' (read/write), 'w+' (create/overwrite)
# Slicing a memmap returns a normal ndarray (copy) or a memmap view
# depending on the operation.

# Useful for processing huge datasets that don't fit in RAM
# (e.g. iterating over chunks of a multi-GB array)

import numpy as np
import pandas as pd

# NumPy <-> pandas
df = pd.DataFrame({"a": [1, 2, 3], "b": [4.0, 5.0, 6.0]})
arr = df.to_numpy()      # 2-D float array
col = df["a"].to_numpy() # 1-D int array
df2 = pd.DataFrame(arr, columns=["x", "y"])

# NumPy <-> PyTorch
import torch
t = torch.from_numpy(arr)     # shares memory (no copy)
back = t.numpy()              # only works for CPU tensors

# NumPy <-> TensorFlow
import tensorflow as tf
tensor = tf.constant(arr)
back = tensor.numpy()

# NumPy <-> JAX
import jax.numpy as jnp
j_arr = jnp.asarray(arr)      # JAX array (immutable)

# NumPy <-> bytes (for serialization / network)
buf = arr.tobytes()
restored = np.frombuffer(buf, dtype=arr.dtype).reshape(arr.shape)

# Convert list of arrays to a single array (must be same shape)
list_of_arrays = [np.array([1, 2]), np.array([3, 4])]
np.array(list_of_arrays)   # shape (2, 2)

import numpy as np

# Define a structured dtype
dt = np.dtype([("name", "U20"), ("age", "i4"), ("score", "f4")])
data = np.array([("Alice", 30, 9.5),
                 ("Bob", 25, 8.0),
                 ("Carol", 35, 9.8)], dtype=dt)

# Save / load (binary preserves structure)
np.save("people.npy", data)
loaded = np.load("people.npy")
print(loaded["name"])   # ['Alice' 'Bob' 'Carol']
print(loaded["age"].mean())   # 30.0

# CSV with mixed types: better to use pandas
import pandas as pd
df = pd.DataFrame(data)
df.to_csv("people.csv", index=False)

# Convert a structured array to a regular 2-D array
arr_2d = data.view([("age", "i4"), ("score", "f4")]).reshape(-1)  # view
# or just stack the numeric columns
numeric = np.column_stack([data["age"], data["score"]])

# Sort by a field
sorted_data = np.sort(data, order="score")   # ascending by score
sorted_data[-1]["name"]   # 'Carol'

import numpy as np

# From ISO strings
dates = np.array(["2024-01-15", "2024-02-20", "2024-03-25"],
                 dtype="datetime64[D]")
# dtype [D] = day precision

# Different precisions
np.datetime64("2024-01-15")              # day precision
np.datetime64("2024-01-15T10:30")        # minute precision
np.datetime64("2024-01-15T10:30:45.123") # ms precision

# Array of datetimes
np.array(["2024-01", "2024-02", "2024-03"], dtype="datetime64[M]")

# Arange with datetime
np.arange("2024-01", "2024-06", dtype="datetime64[M]")
# ['2024-01' '2024-02' '2024-03' '2024-04' '2024-05']

# Today (UTC)
np.datetime64("today", "D")
np.datetime64("now", "s")

import numpy as np

start = np.datetime64("2024-01-01")
end = np.datetime64("2024-12-31")

# Difference is a timedelta64
diff = end - start
print(diff, diff.dtype)   # 365 days, timedelta64[D]

# Add / subtract timedeltas
start + np.timedelta64(7, "D")       # 2024-01-08
start + np.timedelta64(1, "M")       # 2024-02-01 (month added)
start + np.timedelta64(2, "h")       # 2024-01-01T02:00

# Business day arithmetic
np.busday_offset("2024-01-15", 1)     # next business day
np.is_busday("2024-01-15")            # Monday -> True
np.is_busday("2024-01-13")            # Saturday -> False

# Count business days in a range
np.busday_count("2024-01-01", "2024-02-01")   # 23

# Array of business days
np.arange("2024-01-01", "2024-01-10", dtype="datetime64[D]")
# [Mon Tue Wed Thu Fri Sat Sun Mon Tue]

import numpy as np

# Generate a range of dates
daily = np.arange("2024-01-01", "2024-02-01", dtype="datetime64[D]")
# 31 daily dates

monthly = np.arange("2024-01", "2025-01", dtype="datetime64[M]")
# 12 monthly dates

hourly = np.arange("2024-01-01T00", "2024-01-02T00", dtype="datetime64[h]")
# 24 hourly dates

# Resample daily to monthly: take last value of each month
# NumPy itself doesn't have a resample API - use pandas
import pandas as pd
ts = pd.Series(np.random.rand(31), index=daily)
monthly_mean = ts.resample("ME").mean()    # one value per month

# Compute the weekday of each date
weekdays = (daily.astype("datetime64[D]").view("int64") - 4) % 7
# 0 = Monday, ..., 6 = Sunday

# Or use pandas for richer datetime features
pd.DatetimeIndex(daily).day_name()   # ['Monday' 'Tuesday' ...]

import numpy as np
import datetime as dt

# NumPy datetime64 -> Python datetime
np_dt = np.datetime64("2024-01-15T10:30:00")
py_dt = np_dt.astype(dt.datetime)   # Python datetime object
# datetime.datetime(2024, 1, 15, 10, 30)

# Python datetime -> NumPy datetime64
np_again = np.datetime64(py_dt)

# Convert a NumPy datetime array to Python datetimes
arr = np.arange("2024-01-01", "2024-01-05", dtype="datetime64[D]")
py_list = arr.astype(dt.datetime).tolist()
# [datetime.date(2024, 1, 1), ...]

# Pandas Timestamp (most flexible)
import pandas as pd
ts = pd.Timestamp("2024-01-15 10:30:00", tz="US/Eastern")
ts.to_numpy()   # converts to UTC datetime64 (tz info dropped!)

# Unix timestamp interconversion
ts_unix = (np_dt - np.datetime64("1970-01-01", "s")) / np.timedelta64(1, "s")
# 1705314600.0   (seconds since Unix epoch)
back = np.datetime64(int(ts_unix), "s").astype("datetime64[s]") + \
       np.timedelta64(int((ts_unix % 1) * 1e9), "ns")

import numpy as np

# Compute age in years from birthdates
birthdays = np.array(["1990-05-15", "1985-11-20", "2000-03-10"],
                     dtype="datetime64[D]")
today = np.datetime64("today", "D")
ages_days = (today - birthdays).astype("timedelta64[D]")
ages_years = (ages_days / 365.25).astype(int)
# array([33, 38, 24]) approximately

# Find dates within a window
dates = np.arange("2024-01-01", "2024-06-01", dtype="datetime64[D]")
mask = (dates >= "2024-02-01") & (dates < "2024-03-01")
february = dates[mask]   # all February dates

# Compute day of year
jan1 = dates.astype("datetime64[Y]").astype("datetime64[D]")
# Hmm, this loses month/day; instead:
day_of_year = (dates - dates.astype("datetime64[M]") \
               .astype("datetime64[D]") + 1).astype(int)

# Group by month (manual; pandas is easier)
months = dates.astype("datetime64[M]")
unique_months = np.unique(months)
for m in unique_months:
    mask = months == m
    print(m, dates[mask].size, "days")

import numpy as np
names = np.array(["Alice", "Bob", "Carol", "dave"])

# np.char module: vectorized string operations
np.char.upper(names)      # ['ALICE' 'BOB' 'CAROL' 'DAVE']
np.char.lower(names)      # ['alice' 'bob' 'carol' 'dave']
np.char.title(names)      # ['Alice' 'Bob' 'Carol' 'Dave']
np.char.capitalize(names)
np.char.swapcase(names)

# Strip whitespace
s = np.array(["  hello  ", "world  "])
np.char.strip(s)          # ['hello' 'world']
np.char.lstrip(s)
np.char.rstrip(s)

# Splitting and joining
np.char.split(np.array(["a,b,c", "d,e,f"]), sep=",")
# [list(['a', 'b', 'c']) list(['d', 'e', 'f'])]

np.char.join("-", np.array(["abc", "def"]))
# ['a-b-c' 'd-e-f']

# Replace
np.char.replace(np.array(["foo bar", "baz foo"]), "foo", "XXX")
# ['XXX bar' 'baz XXX']

import numpy as np
words = np.array(["apple", "banana", "cherry", "date"])

# Element-wise comparison
np.char.equal(words, "apple")        # [ True False False False]
np.char.not_equal(words, "apple")

# Lexicographic comparison
np.char.greater(words, "c")          # [False False  True  True]

# Substring search
np.char.find(words, "a")             # index of 'a', -1 if not found
np.char.count(np.array(["banana", "apple"]), "a")
# [3 1]

# Startswith / endswith
np.char.startswith(words, "a")       # [ True False False False]
np.char.endswith(words, "e")         # [ True False False  True]

# Contains (regex-capable in newer versions)
np.char.find(words, "er") >= 0       # [False False  True False]

# String length
np.char.str_len(words)               # [5 6 6 4]

import numpy as np
first = np.array(["Alice", "Bob", "Carol"])
last = np.array(["Smith", "Jones", "Brown"])

# Element-wise concatenation
np.char.add(first, last)              # ['AliceSmith' 'BobJones' 'CarolBrown']
np.char.add(first, " " + last)        # nope, this broadcasts weirdly

# Proper way to add a separator
np.char.add(np.char.add(first, " "), last)
# ['Alice Smith' 'Bob Jones' 'Carol Brown']

# Multiply (repeat)
np.char.multiply(np.array(["ab", "cd"]), 3)
# ['ababab' 'cdcdcd']

# Center / ljust / rjust with fill character
np.char.center(np.array(["hi", "yo"]), width=6, fillchar="-")
# ['--hi--' '--yo--']
np.char.rjust(np.array(["1", "12", "123"]), width=5, fillchar="0")
# ['00001' '00012' '00123']

# Format strings
np.char.mod("%05d", np.array([1, 23, 456]))
# ['00001' '00023' '00456']

# Modern Python f-strings don't vectorize - use np.char.mod or list comp
[f"{n:05d}" for n in [1, 23, 456]]   # Python alternative

import numpy as np

# Fixed-length unicode strings (default)
s = np.array(["hello", "world"])
s.dtype    # dtype('<U5')   <- max length 5 unicode chars

# Fixed-width string dtype (saves memory for ASCII)
s_ascii = np.array(["hello", "world"], dtype="S5")
s_ascii.dtype   # dtype('|S5')  <- bytes, length 5

# If a string exceeds the fixed width, it's truncated
np.array(["abcdef"], dtype="U3")   # ['abc']
np.array(["abcdef"], dtype="S3")   # [b'abc']

# Object dtype for variable-length strings (no truncation, slower)
s_obj = np.array(["a" * 100, "b"], dtype=object)
s_obj.dtype   # dtype('O')

# Convert between unicode and bytes
np.char.encode(np.array(["hello"]), "utf-8")   # array([b'hello'])
np.char.decode(np.array([b"hello"]), "utf-8")  # array(['hello'])

# When loading text data, dtype=None infers per column
# but for strings use U/S for fixed-width or object for variable

import numpy as np

# Build a categorical label array
labels = np.array(["cat", "dog", "cat", "bird", "dog"])

# Encode strings to integers
unique, encoded = np.unique(labels, return_inverse=True)
print(unique)     # ['bird' 'cat' 'dog']
print(encoded)    # [1 2 1 0 2]   integer codes

# Decode integers back to strings
labels_back = unique[encoded]
# ['cat' 'dog' 'cat' 'bird' 'dog']

# Filter rows of a 2-D array by a string column
data = np.array([["Alice", "30"], ["Bob", "25"], ["Carol", "35"]])
mask = np.char.startswith(data[:, 0], "A")
data[mask]   # [['Alice' '30']]

# Pad numbers as zero-padded strings (for sorting filenames)
nums = np.array([1, 10, 100, 2])
padded = np.char.zfill(nums.astype(str), 3)
# ['001' '010' '100' '002']   now sorts lexicographically

# Compute a sorted argsort using string keys
keys = np.array(["banana", "apple", "cherry"])
order = np.argsort(keys)
keys[order]   # ['apple' 'banana' 'cherry']

# Replace multiple patterns (must chain - no vectorized replace-all)
arr = np.array(["a-b", "c_d", "e.f"])
for ch in ["-", "_", "."]:
    arr = np.char.replace(arr, ch, " ")
print(arr)   # ['a b' 'c d' 'e f']

import numpy as np
import time

# BAD: Python loop (slow)
def slow_sum(arr):
    total = 0
    for x in arr:
        total += x
    return total

# GOOD: NumPy reduction (vectorized)
def fast_sum(arr):
    return arr.sum()

big = np.arange(10_000_000)

start = time.perf_counter()
slow_sum(big)
print("loop:", time.perf_counter() - start)   # ~0.6s

start = time.perf_counter()
fast_sum(big)
print("numpy:", time.perf_counter() - start)  # ~0.005s   ~100x faster

# Rule of thumb: if you're writing a for loop over array elements,
# there's probably a vectorized NumPy operation that's 50-100x faster.
# Look for: reductions (sum, mean, max), ufuncs, broadcasting, einsum.

import numpy as np

# einsum: Einstein summation - the Swiss Army knife of tensor ops
A = np.random.rand(3, 4)
B = np.random.rand(4, 5)

# Matrix multiplication: C[i,k] = sum_j A[i,j] * B[j,k]
C = np.einsum("ij,jk->ik", A, B)   # same as A @ B

# Transpose
np.einsum("ij->ji", A)             # A.T

# Sum all elements
np.einsum("ij->", A)               # A.sum()

# Sum along an axis
np.einsum("ij->i", A)              # A.sum(axis=1)
np.einsum("ij->j", A)              # A.sum(axis=0)

# Row-wise dot product with a vector
v = np.random.rand(4)
np.einsum("ij,j->i", A, v)         # A @ v

# Outer product
a = np.random.rand(3)
b = np.random.rand(4)
np.einsum("i,j->ij", a, b)         # np.outer(a, b)

# Batched matrix multiply
X = np.random.rand(10, 3, 4)
Y = np.random.rand(10, 4, 5)
np.einsum("bij,bjk->bik", X, Y)    # batched matmul

# optimize=True lets einsum choose a better evaluation order
np.einsum("ij,jk,kl->il", A, B, A.T, optimize=True)

import numpy as np
import time

# Many ufuncs accept out= to write into an existing buffer (no allocation)
a = np.arange(1_000_000)
out = np.empty(1_000_000)

# Slow: creates a new array
start = time.perf_counter()
b = np.sin(a)
print("alloc:", time.perf_counter() - start)

# Fast: reuses the existing buffer
start = time.perf_counter()
np.sin(a, out=out)
print("reuse:", time.perf_counter() - start)   # faster, no allocation

# In-place operations also avoid copies
a += 10            # in-place add (no temp array)
a *= 2             # in-place multiply
np.sqrt(a, out=a)  # in-place sqrt

# Views vs copies
b = a.reshape(-1, 1000)        # view, no copy
c = a[a > 0.5]                  # copy (fancy indexing)
d = a.transpose()               # view (but maybe non-contiguous)

# ascontiguousarray: copy only if needed
e = np.ascontiguousarray(a)     # no copy if already C-contiguous

import numpy as np
import time

# C order (row-major, default): last axis varies fastest
# F order (column-major): first axis varies fastest

# Iterating along the LAST axis is cache-friendly (C order)
a_c = np.zeros((1000, 1000), order="C")
a_f = np.zeros((1000, 1000), order="F")

# Sum along rows (axis 1): fast for C-order (contiguous in memory)
start = time.perf_counter()
a_c.sum(axis=1)
print("C order, axis=1:", time.perf_counter() - start)

start = time.perf_counter()
a_f.sum(axis=1)
print("F order, axis=1:", time.perf_counter() - start)   # slower

# Rule: for C-order arrays, the last axis is contiguous.
# Iterating a[..., i] is fast; iterating a[i, ...] is slow.

# Force a particular layout (copies only if necessary)
a_c2 = np.ascontiguousarray(a_f)   # now C-contiguous
a_f2 = np.asfortranarray(a_c)      # now F-contiguous

# Check layout
a_c.flags["C_CONTIGUOUS"]    # True
a_c.flags["F_CONTIGUOUS"]    # False

# Strides: bytes to step along each axis
a_c.strides    # (8000, 8)   - 8 bytes per float, 1000 per row
a_f.strides    # (8, 8000)   - first axis is contiguous

import numpy as np
import time

# Quick micro-benchmark
def bench(fn, *args, repeat=5):
    times = []
    for _ in range(repeat):
        start = time.perf_counter()
        fn(*args)
        times.append(time.perf_counter() - start)
    return min(times)

a = np.random.rand(1000, 1000)
print("sum axis=0:", bench(lambda x: x.sum(axis=0), a))
print("sum axis=1:", bench(lambda x: x.sum(axis=1), a))

# Memory profiling: check array size
print(f"size: {a.size} elements")
print(f"memory: {a.nbytes / 1e6:.1f} MB")

# Time individual lines with %timeit (in IPython/Jupyter)
# %timeit a.sum()
# %timeit np.einsum("ij->", a)
# %prun my_function()    # line-by-line profiling

# Find bottlenecks with cProfile
import cProfile
cProfile.run("a @ a", sort="cumtime")

# Check if BLAS is using multiple threads
np.show_config()    # shows BLAS/LAPACK info

# Control thread count (sometimes single-thread is faster for small arrays)
import os
os.environ["OMP_NUM_THREADS"] = "4"
os.environ["MKL_NUM_THREADS"] = "4"

import numpy as np

# Generate a signal: two sine waves + noise
t = np.linspace(0, 1, 500, endpoint=False)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
signal += 0.2 * np.random.default_rng(42).standard_normal(500)

# Compute the FFT
fft_vals = np.fft.fft(signal)
# Complex array: magnitude = amplitude, phase = angle

# Frequencies corresponding to each FFT bin
freqs = np.fft.fftfreq(len(signal), d=t[1] - t[0])

# Magnitude spectrum (absolute value)
magnitude = np.abs(fft_vals)

# Power spectral density
power = magnitude ** 2

# Plot the spectrum (only positive frequencies for real signals)
import matplotlib.pyplot as plt
plt.plot(freqs[:len(freqs)//2], magnitude[:len(freqs)//2])
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")

import numpy as np

# Original signal
t = np.linspace(0, 1, 500, endpoint=False)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)

# Forward FFT
fft_vals = np.fft.fft(signal)

# Inverse FFT reconstructs the original signal
reconstructed = np.fft.ifft(fft_vals)
np.allclose(signal, reconstructed.real)   # True

# Modify the spectrum (e.g. zero out high frequencies) then invert
fft_filtered = fft_vals.copy()
freqs = np.fft.fftfreq(len(signal), d=t[1] - t[0])
fft_filtered[np.abs(freqs) > 10] = 0   # low-pass filter

filtered_signal = np.fft.ifft(fft_filtered).real
# Now only the 5 Hz component remains; the 20 Hz is removed

# rfft / irfft: optimized for real input (half the storage)
rfft_vals = np.fft.rfft(signal)   # only non-negative frequencies
reconstructed2 = np.fft.irfft(rfft_vals, n=len(signal))

import numpy as np

# 2-D FFT for image processing
image = np.random.rand(64, 64)   # treat as a grayscale image

# Forward 2-D FFT
fft2d = np.fft.fft2(image)

# Shift zero frequency to the center (standard convention)
fft2d_shifted = np.fft.fftshift(fft2d)

# Magnitude spectrum (log-scaled for visualization)
magnitude = np.log(np.abs(fft2d_shifted) + 1)

# Inverse: unshift then ifft2
fft2d_unshifted = np.fft.ifftshift(fft2d_shifted)
reconstructed = np.fft.ifft2(fft2d_unshifted).real
np.allclose(image, reconstructed)   # True

# Low-pass filter: keep only central low frequencies
rows, cols = image.shape
crow, ccol = rows // 2, cols // 2
mask = np.zeros((rows, cols), dtype=bool)
mask[crow-10:crow+10, ccol-10:ccol+10] = True   # central 20x20 block
fft2d_shifted[~mask] = 0
blurred = np.fft.ifft2(np.fft.ifftshift(fft2d_shifted)).real

import numpy as np

# Applying a window reduces spectral leakage at the edges of the FFT
N = 1024
t = np.linspace(0, 1, N, endpoint=False)
signal = np.sin(2 * np.pi * 50 * t)

# Rectangular window (no window) - lots of leakage
rect_fft = np.abs(np.fft.fft(signal))

# Hann window - smooth taper to zero at both ends
hann = np.hanning(N)
signal_hann = signal * hann
hann_fft = np.abs(np.fft.fft(signal_hann))

# Other common windows
hamming = np.hamming(N)
blackman = np.blackman(N)
bartlett = np.bartlett(N)

# Compare spectral leakage: Hann has wider main lobe but lower sidelobes
# Trade-off: frequency resolution vs leakage

# For overlap-add processing, use windows that sum to 1 when overlapped:
# Hann at 50% overlap -> COLA (constant overlap-add) property

# Short-time Fourier transform (STFT) via windowed FFTs
def stft(signal, window_size=256, hop=128):
    n_frames = (len(signal) - window_size) // hop + 1
    window = np.hanning(window_size)
    frames = []
    for i in range(n_frames):
        start = i * hop
        frame = signal[start:start + window_size] * window
        frames.append(np.fft.rfft(frame))
    return np.array(frames).T   # shape (n_freqs, n_frames)

import numpy as np

# Convolution theorem: convolution in time = multiplication in frequency
# For large signals, FFT-based convolution is O(N log N) vs O(N^2) direct

signal = np.random.rand(10000)
kernel = np.random.rand(100)

# Method 1: np.convolve (direct, O(N*K))
direct_conv = np.convolve(signal, kernel, mode="same")

# Method 2: FFT-based convolution (fast for large kernels)
def fft_convolve(signal, kernel):
    n = len(signal) + len(kernel) - 1   # full convolution length
    n_fft = 1 << int(np.ceil(np.log2(n)))   # next power of 2 (faster FFT)
    fft_signal = np.fft.rfft(signal, n_fft)
    fft_kernel = np.fft.rfft(kernel, n_fft)
    fft_conv = fft_signal * fft_kernel
    conv = np.fft.irfft(fft_conv, n_fft)[:n]
    return conv

fft_result = fft_convolve(signal, kernel)
# Trim or pad to match 'same' mode if needed

# scipy.signal.fftconvolve is a production-ready version
from scipy.signal import fftconvolve
scipy_conv = fftconvolve(signal, kernel, mode="same")

# Cross-correlation (flip the kernel first)
xcorr = fftconvolve(signal, kernel[::-1], mode="same")

# Autocorrelation
autocorr = fftconvolve(signal, signal[::-1], mode="same")

import numpy as np
import numpy.ma as ma

# A masked array = data + a boolean mask (True = invalid/missing)
data = np.array([1, 2, -999, 4, 5])
mask = np.array([False, False, True, False, False])
m = ma.masked_array(data, mask=mask)
# masked_array(data=[1, 2, --, 4, 5],
#              mask=[False, False,  True, False, False])

# Common constructors
ma.masked_equal(data, -999)        # mask all elements == -999
ma.masked_greater(data, 3)         # mask elements > 3
ma.masked_where(data > 3, data)    # general: mask where condition
ma.masked_invalid([1.0, np.nan, 3.0])   # mask NaN/inf

# Access components
m.data       # the underlying data (including masked values)
m.mask       # the boolean mask
m.filled(0)  # replace masked values with 0, return regular ndarray

# Check if any value is masked
m.mask.any()    # True

import numpy as np
import numpy.ma as ma

m = ma.masked_array([1, 2, 3, 4, 5],
                    mask=[False, False, True, False, False])

# Reductions ignore masked values
m.sum()      # 12   (skips the masked 3)
m.mean()     # 3.0  (mean of unmasked only)
m.max()      # 5
m.min()      # 1

# Arithmetic propagates the mask
m + 10       # [11 12 -- 14 15]
m * 2        # [2 4 -- 8 10]

# Combine masks from two masked arrays
m2 = ma.masked_array([10, 20, 30, 40, 50],
                     mask=[True, False, False, False, True])
m + m2       # [-- 22 -- 44 --]   union of masks

# Logical operations
m < 4        # masked_array, masked entries stay masked

# Compress: extract only the unmasked values (returns 1-D array)
m.compressed()   # [1 2 4 5]   (regular ndarray, no mask)

# Count valid entries
ma.count(m)      # 4
m.count()        # 4

import numpy as np
import numpy.ma as ma

m = ma.masked_array([1.0, 2.0, 3.0, 4.0],
                    mask=[False, True, False, True])

# filled(): replace masked entries with a fill value
m.filled(0)         # [1. 0. 3. 0.]   - regular ndarray
m.filled(np.nan)    # [ 1. nan  3. nan]
m.filled(m.mean())  # fill with mean of unmasked

# The default fill value depends on dtype
m.fill_value    # 1e+20 for floats, 999999 for ints, etc.

# Set a custom fill value
m2 = ma.masked_array([1, 2, 3], mask=[0, 1, 0], fill_value=-1)

# Convert to a plain ndarray (loses mask info!)
arr = np.asarray(m)         # includes masked values! [1. 2. 3. 4.]
arr_filled = np.asarray(m.filled(np.nan))   # NaN where masked

# Convert a regular array with sentinel values to masked
data = np.array([1, -999, 3, -999, 5])
m3 = ma.masked_equal(data, -999)

# Convert masked back to array with NaN (for plotting)
arr_nan = ma.filled(m, np.nan)

import numpy as np
import numpy.ma as ma

m = ma.masked_array(np.arange(12).reshape(3, 4),
                    mask=[[0, 0, 1, 0],
                          [0, 1, 0, 0],
                          [1, 0, 0, 0]])

# Indexing works like regular arrays
m[0]            # first row (mask preserved)
m[1, 2]         # scalar value (or masked if masked)
m[:, 1]         # column 1

# Boolean mask on top of masked array
m[m.data > 5]   # values > 5 (mask still applied)

# Set values (and unmask)
m[0, 0] = 99
m[0, 0] = ma.masked   # re-mask a specific element

# Reset the mask entirely
m.mask = ma.nomask     # no values masked
m[1, 1] = ma.masked    # mask a single value

# Check element-wise
m[0, 2] is ma.masked   # True if that element is masked

# Logical operations keep the mask
(m > 5).data    # boolean array
(m > 5).mask    # propagated mask

# Iteration skips masked values (kind of)
list(m.compressed())   # all unmasked values, flattened

import numpy as np
import numpy.ma as ma

# Scenario 1: Sentinel values in data (e.g. -999 for missing)
temps = np.array([23.5, -999, 24.0, -999, 22.5])
m = ma.masked_equal(temps, -999)
m.mean()   # 23.33...  (ignores -999 entries)

# Scenario 2: Quality flags / validity
values = np.array([10, 12, 8, 15, 11])
quality = np.array([1, 1, 0, 1, 0])   # 0 = bad
m = ma.masked_array(values, mask=(quality == 0))
m.mean()   # (10 + 12 + 15) / 3 = 12.33

# Scenario 3: Grid data with land/sea mask
sst = np.random.rand(5, 5) * 30   # sea surface temperature
land_mask = np.array([
    [1, 1, 0, 0, 0],
    [1, 0, 0, 0, 1],
    [0, 0, 1, 0, 0],
    [0, 1, 0, 0, 1],
    [1, 0, 0, 1, 1],
], dtype=bool)
sst_masked = ma.masked_array(sst, mask=land_mask)
sst_masked.mean()   # average over ocean only

# When NOT to use: if you just have NaN, regular arrays work fine
# (and np.nanmean etc. handle them). Use masked arrays when:
# - You need integer dtype (NaN forces float)
# - You have a non-NaN sentinel (-999, 0, etc.)
# - The mask has semantic meaning (quality flag, region)

import numpy as np

# Define a structured dtype: list of (name, format) tuples
dt = np.dtype([("name", "U20"),
               ("age", "i4"),
               ("weight", "f8"),
               ("active", "?")])    # ? = bool

# Create an array with this dtype
people = np.array([("Alice", 30, 65.5, True),
                   ("Bob", 25, 80.0, False),
                   ("Carol", 35, 55.0, True)], dtype=dt)

print(people)
# [('Alice', 30, 65.5,  True) ('Bob', 25, 80. , False) ...]

# Each element is a record (tuple-like)
people[0]              # ('Alice', 30, 65.5, True)
people[0]["name"]      # 'Alice'

# Type codes: i1/i2/i4/i8 (ints), u1..u8 (unsigned),
#             f2/f4/f8 (floats), U/n (strings), ? (bool), M8 (datetime)

import numpy as np
dt = np.dtype([("name", "U20"), ("age", "i4"), ("score", "f4")])
data = np.array([("Alice", 30, 9.5),
                 ("Bob", 25, 8.0),
                 ("Carol", 35, 9.8)], dtype=dt)

# Access a single field across all records (returns a view)
data["name"]    # ['Alice' 'Bob' 'Carol']
data["age"]     # [30 25 35]
data["age"].mean()   # 30.0

# Modify a field
data["age"] += 1    # increment everyone's age
data["score"] = 0   # set all scores to 0

# Access multiple fields at once
data[["name", "score"]]   # new structured array with just those fields

# Field assignment with another array
data["score"] = [9.5, 8.0, 9.8]

# Iterate records (each record is a tuple)
for rec in data:
    print(rec["name"], rec["age"])

# Iterate fields
for name in data.dtype.names:
    print(name, data[name]

import numpy as np
dt = np.dtype([("name", "U20"), ("age", "i4"), ("score", "f4")])
data = np.array([("Alice", 30, 9.5),
                 ("Bob", 25, 8.0),
                 ("Carol", 35, 9.8),
                 ("Dave", 28, 7.5)], dtype=dt)

# Sort by a single field
np.sort(data, order="age")         # ascending by age
np.sort(data, order="score")[::-1] # descending by score

# Sort by multiple fields (lexicographic)
np.sort(data, order=["score", "age"])  # primary: score, secondary: age

# argsort with order
idx = np.argsort(data, order="score")
data[idx]   # records sorted by score

# Filter by a field
mask = data["age"] > 28
data[mask]   # Alice, Carol

# Complex conditions
mask = (data["score"] > 8) & (data["age"] < 35)
data[mask]   # Alice

# Find unique values of a field
np.unique(data["age"])   # [25 28 30 35]

# Aggregate per field
data["age"].mean()
data["score"].std()

import numpy as np

# Record arrays: access fields as ATTRIBUTES (data.age vs data["age"])
dt = np.dtype([("name", "U20"), ("age", "i4"), ("score", "f4")])

# Method 1: view as recarray
data = np.array([("Alice", 30, 9.5),
                 ("Bob", 25, 8.0)], dtype=dt).view(np.recarray)

data.name    # ['Alice' 'Bob']   - attribute access!
data.age     # [30 25]
data.age.mean()

# Method 2: create directly with np.rec.array
recs = np.rec.array([("Alice", 30, 9.5),
                     ("Bob", 25, 8.0)],
                    dtype=[("name", "U20"), ("age", "i4"), ("score", "f4")])
recs.name    # attribute access
recs[0].name # 'Alice'

# Trade-off: recarray is slightly slower and uses a bit more memory
# than a plain structured array, but the syntax is cleaner.

# In practice, prefer plain structured arrays for performance,
# recarray only when you really want attribute access (e.g. in REPL).

import numpy as np
import pandas as pd

# Same data as a structured array vs a pandas DataFrame
dt = np.dtype([("name", "U20"), ("age", "i4"), ("score", "f4")])
sa = np.array([("Alice", 30, 9.5),
               ("Bob", 25, 8.0),
               ("Carol", 35, 9.8)], dtype=dt)

df = pd.DataFrame(sa)   # easy conversion
#    name  age  score
# 0  Alice   30    9.5
# 1    Bob   25    8.0
# 2  Carol   35    9.8

# Conversion back
sa2 = df.to_records(index=False)

# When to use structured arrays (NumPy):
# - Memory is tight (no pandas overhead, ~2-10x smaller)
# - Interoperating with C/Fortran code that expects structs
# - Simple, fixed-schema numeric data
# - Embedded / scientific computing without pandas available

# When to use pandas:
# - Heterogeneous data with mixed types
# - Need groupby, join, resample, time series
# - Handling missing data (NaN handling is richer)
# - Reading/writing CSV, Excel, SQL, Parquet
# - Most data analysis workflows

# Pandas is built on top of NumPy - usually the better default.
# Structured arrays are a low-level escape hatch.

import numpy as np
A = np.random.rand(3, 4)
B = np.random.rand(4, 5)
v = np.random.rand(4)

# Matrix multiply
np.einsum("ij,jk->ik", A, B)   # = A @ B

# Matrix-vector multiply
np.einsum("ij,j->i", A, v)     # = A @ v

# Vector inner product (dot)
np.einsum("i,i->", v, v)       # = v @ v = sum(v**2)

# Outer product
np.einsum("i,j->ij", v, v)     # = np.outer(v, v)

# Trace
np.einsum("ii->", A) if A.shape[0] == A.shape[1] else None

# Diagonal
M = np.random.rand(5, 5)
np.einsum("ii->i", M)          # = np.diag(M)

# Transpose
np.einsum("ij->ji", A)         # = A.T

# Sum all elements
np.einsum("ij->", A)           # = A.sum()

# Sum along axis 0
np.einsum("ij->j", A)          # = A.sum(axis=0)

# Element-wise (Hadamard) product
np.einsum("ij,ij->ij", A, A)   # = A * A

# Batched matrix multiply
X = np.random.rand(10, 3, 4)
Y = np.random.rand(10, 4, 5)
np.einsum("bij,bjk->bik", X, Y)   # batched A @ B

import numpy as np

# np.lib.stride_tricks.as_strided: create views with custom strides
# Extremely powerful but DANGEROUS - can read out-of-bounds memory

a = np.arange(10)
print(a.strides)   # (8,)   - 8 bytes per int64

# Create a sliding-window view WITHOUT copying data
from numpy.lib.stride_tricks import sliding_window_view, as_strided

# Modern, safe way (NumPy >= 1.20):
windows = sliding_window_view(a, 3)
# shape (8, 3): [[0,1,2], [1,2,3], ..., [7,8,9]]
# This is a VIEW - shares memory with a

# Compute moving average using sliding windows
data = np.random.rand(1000)
window = 5
ma = sliding_window_view(data, window).mean(axis=1)
# shape (996,) - moving average

# The low-level (DANGEROUS) as_strided
# shape = (8, 3), strides = (8, 8)  -> step 1 element along both axes
windows2 = as_strided(a, shape=(8, 3), strides=(8, 8))
# Same as sliding_window_view(a, 3) but you must compute shape/strides
# yourself - get it wrong and you read garbage or crash.

# ALWAYS prefer sliding_window_view - it computes strides safely.

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view

# 1-D rolling window
data = np.arange(10)
windows = sliding_window_view(data, 4)
# shape (7, 4): [[0,1,2,3], [1,2,3,4], ..., [6,7,8,9]]

# Rolling mean (window=4)
rolling_mean = windows.mean(axis=1)
# [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5]

# Rolling max / min / sum
rolling_max = windows.max(axis=1)
rolling_min = windows.min(axis=1)
rolling_sum = windows.sum(axis=1)

# 2-D sliding window (e.g. for image processing)
img = np.random.rand(10, 10)
patches = sliding_window_view(img, (3, 3))
# shape (8, 8, 3, 3): each (i,j) is a 3x3 patch centered around (i+1, j+1)

# Compute patch-wise mean (blur effect)
blurred = patches.mean(axis=(-2, -1))
# shape (8, 8)

# Convolution-like operation (no padding, stride 1)
kernel = np.array([[1, 0, -1], [1, 0, -1], [1, 0, -1]])
convolved = (patches * kernel).sum(axis=(-2, -1))

# Window with step (downsample): slice the windows
windowed_stride2 = sliding_window_view(data, 4)[::2]
# shape (4, 4)

import numpy as np
from numpy.lib.stride_tricks import broadcast_arrays, broadcast_to

# broadcast_to: view of an array broadcast to a new shape (read-only)
a = np.array([1, 2, 3])
b = np.broadcast_to(a, (4, 3))
# shape (4, 3), all rows = [1, 2, 3], no data copied
# b[0, 0] = 99   # ERROR - read-only

# broadcast_arrays: broadcast several arrays together (as views)
x = np.array([1, 2, 3])           # (3,)
y = np.array([[10], [20]])        # (2, 1)
bx, by = broadcast_arrays(x, y)
bx.shape   # (2, 3)
by.shape   # (2, 3)
# bx, by share memory with x, y (just with new strides)

# Use case: build coordinate grids without copying
xs = np.array([0, 1, 2])
ys = np.array([0, 1, 2, 3])
X, Y = np.meshgrid(xs, ys)        # creates copies
X2, Y2 = broadcast_arrays(xs, ys) # views - no copies!
# X2, Y2 are read-only views

# ogrid / mgrid: open / dense meshgrids
ox, oy = np.ogrid[0:3, 0:4]       # column + row vectors (views)
ox.shape   # (3, 1)
oy.shape   # (1, 4)
# Use for broadcasting computations: X = ox + oy  (shape (3, 4))

mx, my = np.mgrid[0:3, 0:4]       # dense (copies), shape (3, 4) each

import numpy as np
from numpy.lib.stride_tricks import as_strided

# Use case: compute pairwise distances without forming the full matrix
# (still O(N^2) compute, but at least no extra copy of the data)
points = np.random.rand(1000, 3)   # 1000 points in 3-D

# Method 1: explicit broadcasting (creates intermediates)
# diff = points[:, None, :] - points[None, :, :]   # (1000, 1000, 3) = 24MB
# dists = np.linalg.norm(diff, axis=-1)            # (1000, 1000)

# Method 2: use as_strided to view points as (1000, 1, 3) and (1, 1000, 3)
# without copying - same compute, less peak memory
n, d = points.shape
p1 = as_strided(points, shape=(n, 1, d),
                strides=(points.strides[0], 0, points.strides[1]))
p2 = as_strided(points, shape=(1, n, d),
                strides=(0, points.strides[0], points.strides[1]))
# p1 and p2 share memory with points (zero-copy broadcast views)
# diff = p1 - p2  still allocates (n, n, d), but the inputs didn't

# Compute a Gram matrix (X @ X.T) with as_strided
X = np.random.rand(100, 5)
gram = np.einsum("ij,ik->jk", X, X)   # = X.T @ X, shape (5, 5)

# Tile an array virtually (without copying)
a = np.array([1, 2, 3])
tiled_view = as_strided(a, shape=(4, 3),
                        strides=(0, a.strides[0]))
# Each row is [1, 2, 3], no data copied