Beginner Advanced 18 min read

Chapter 49: Python for Data Science in Depth

Python has become the dominant language for data science not by accident but by design: its syntax rewards clarity, its ecosystem is vast, and its interoperability with compiled libraries lets you write high-level code that runs at near-C speed. Yet writing good Python for data work — code that is correct, fast, readable, and maintainable — requires understanding several layers simultaneously: the language's own data model, NumPy's memory layout and broadcasting rules, pandas's index-aligned operations and internal block structure, and the performance traps that lurk at every layer.

This chapter treats Python as a precision instrument. We begin with idiomatic core-Python constructs — comprehensions, generators, unpacking, and the data structures that underpin fast analytics — before ascending to NumPy's vectorization model and its broadcasting algebra. We then dissect pandas from the inside out: how its DataFrame is built from typed blocks, why certain operations copy and others view, when apply is your friend and when it is a performance disaster, and how to handle the two most common real-world complications: messy datetime data and high-cardinality categoricals. The chapter closes with a section on writing clean, testable analysis code, because correctness is a prerequisite for performance.

Readers should come away not just knowing what the functions do, but why they behave as they do, which is the only foundation for confident optimization and debugging.

Runnable companion notebook for this chapter. Download the Jupyter notebook (.ipynb)
Libraries covered: NumPyPandas

Learning Objectives

  • Master idiomatic Python data structures (lists, dicts, sets, tuples, namedtuples, dataclasses) and choose among them for data-work tasks.
  • Write and reason about list/dict/set comprehensions and generator expressions, understanding their memory and performance trade-offs.
  • Apply NumPy vectorization and broadcasting to replace explicit loops, and diagnose the memory layout (C-order vs. F-order, strides) that governs performance.
  • Understand pandas internals — block consolidation, copy-vs-view semantics, and the Index object — to write operations that avoid silent copies and SettingWithCopyWarning.
  • Select the right aggregation strategy (vectorized methods, groupby, transform, apply, pipe) and predict their performance characteristics.
  • Handle dates, times, and time zones correctly using Python's datetime, NumPy's datetime64, and pandas's DatetimeTZDtype.
  • Organize analysis code into testable, reproducible modules using functions, type hints, pytest, and reproducible random seeds.

49.1 Idiomatic Python Data Structures and Comprehensions Beginner

Choosing the Right Container

Every data pipeline starts with raw Python objects. Picking the wrong container costs both time and clarity. The four workhorses are:

  • list — ordered, mutable, O(1) append, O(n) membership test.
  • dict — key→value map, O(1) average lookup and insert (hash table). Python 3.7+ preserves insertion order.
  • set — unordered unique values, O(1) membership, O(min(n,m)) intersection.
  • tuple — immutable list; its immutability lets Python intern and hash it, making it suitable as a dict key or set element.

For lightweight records, prefer collections.namedtuple or dataclasses.dataclass over raw dicts when the schema is fixed — you get attribute access, repr, and optional type checking for free.

PYTHON
from dataclasses import dataclass, field
from typing import List

@dataclass
class Experiment:
    name: str
    scores: List[float] = field(default_factory=list)

    @property
    def mean_score(self) -> float:
        return sum(self.scores) / len(self.scores) if self.scores else float('nan')

expt = Experiment(name="baseline", scores=[0.82, 0.79, 0.85])
print(expt.mean_score)   # 0.82

Comprehensions and Generator Expressions

Comprehensions are not merely syntactic sugar — they are compiled to optimised bytecode (LIST<em>APPEND / MAP</em>ADD) that outperforms equivalent append loops by 15–30 % in CPython benchmarks.

PYTHON
# List comprehension: square every odd number
squares = [x**2 for x in range(1_000_000) if x % 2 == 1]

# Dict comprehension: invert a mapping
orig = {'a': 1, 'b': 2, 'c': 3}
inverted = {v: k for k, v in orig.items()}

# Set comprehension: unique domains from email list
emails = ['alice@example.com', 'bob@corp.io', 'carol@example.com']
domains = {e.split('@')[1] for e in emails}

When you only need to iterate over a result rather than store it, replace the square brackets with parentheses to get a generator expression. Generators are lazy: they yield one item at a time and hold O(1) memory regardless of the logical size of the stream.

PYTHON
import csv, pathlib

def parse_prices(path: str):
    """Generator: yields (date_str, float) tuples from a CSV."""
    with open(path) as fh:
        reader = csv.DictReader(fh)
        yield from ((row['date'], float(row['price'])) for row in reader)

total = sum(price for _, price in parse_prices('prices.csv'))

Pitfall: Comprehension Variable Leakage

In Python 3, comprehension variables are properly scoped — x in [x for x in range(5)] does not leak into the enclosing scope. This is the correct behavior; do not rely on the old Python 2 leak.

collections Utilities

The collections module provides several containers that replace verbose patterns:

  • Counter — multiset, replaces {k: counts.get(k,0)+1 for k in items}.
  • defaultdict — eliminates setdefault boilerplate in grouping loops.
  • deque — O(1) popleft for sliding-window algorithms; lists are O(n).
  • ChainMap — a live view over multiple dicts, useful for layered configs.

PYTHON
from collections import Counter, defaultdict

words = 'the quick brown fox jumps over the lazy dog'.split()
word_freq = Counter(words)
print(word_freq.most_common(3))  # [('the', 2), ...]

# Group words by first letter
by_letter = defaultdict(list)
for w in words:
    by_letter[w[0]].append(w)

itertools and functools for Data Pipelines

Built-in higher-order tools compose cleanly with generators to build memory-efficient pipelines:

  • itertools.islice(it, n) — take first n items without materialising the rest.
  • itertools.groupby(sorted<em>iter, key) — streaming group-by (input must be sorted).
  • functools.reduce — fold a sequence; rarely needed since sum/max cover common cases.
  • functools.lru</em>cache — memoize pure functions; critical for recursive feature engineering.

PYTHON
import itertools, functools

# Running product via reduce
from operator import mul
factorial = lambda n: functools.reduce(mul, range(1, n+1), 1)

# Sliding window of size k
def windows(iterable, k):
    it = iter(iterable)
    win = list(itertools.islice(it, k))
    if len(win) == k:
        yield tuple(win)
    for item in it:
        win.append(item)
        win.pop(0)
        yield tuple(win)

Code Examples

Benchmarking list append vs. comprehension

Quantifies the comprehension speedup over an equivalent for-loop with append.

PYTHON
import timeit

stmt_loop = """
result = []
for x in range(10_000):
    if x % 2 == 0:
        result.append(x ** 2)
"""

stmt_comp = "result = [x**2 for x in range(10_000) if x % 2 == 0]"

t_loop = timeit.timeit(stmt_loop, number=1000)
t_comp = timeit.timeit(stmt_comp, number=1000)
print(f"Loop:          {t_loop:.3f}s")
print(f"Comprehension: {t_comp:.3f}s")
print(f"Speedup:       {t_loop / t_comp:.2f}x")
Output
Loop:          0.612s
Comprehension: 0.483s
Speedup:       1.27x

Memory-efficient streaming pipeline with generators

Demonstrates that chained generator expressions hold O(1) memory regardless of N.

PYTHON
import random, sys

N = 1_000_000
data = (random.gauss(0, 1) for _ in range(N))       # generator
positives = (x for x in data if x > 0)              # filter (lazy)
scaled    = (x * 100 for x in positives)            # map   (lazy)

# The entire pipeline consumes O(1) memory
print(sys.getsizeof(scaled))   # tiny — just the generator object

# Materialise only the final result
result = list(scaled)
print(f"Kept {len(result):,} of {N:,} samples")
Output
120
Kept 499,872 of 1,000,000 samples

49.2 NumPy Vectorization and Broadcasting Intermediate

Why Vectorization Matters

CPython executes Python bytecode with per-object dispatch: even a tight for loop over a million floats must box each float into a Python object, test its type, and unbox it again — easily 50–200 ns per element. NumPy stores data as a contiguous C array of a single dtype and exposes ufuncs that process entire arrays in compiled C or Fortran code, often with SIMD instructions. The result is typically 50–500× faster than a pure Python loop on the same arithmetic.

PYTHON
import numpy as np
import timeit

a = np.random.rand(1_000_000)

# Python loop: ~300 ms on a modern laptop
t_loop = timeit.timeit(lambda: [x**2 for x in a], number=10) / 10
# NumPy ufunc: ~1 ms
t_vec  = timeit.timeit(lambda: np.square(a), number=10) / 10

print(f"Loop: {t_loop*1e3:.1f} ms   NumPy: {t_vec*1e3:.2f} ms   Speedup: {t_loop/t_vec:.0f}x")

Memory Layout: Strides and C/F Order

A NumPy array is described by three things: a data buffer, a dtype, and a strides tuple. Strides tell NumPy how many bytes to skip along each axis.

For a 2-D array of shape (R, C) in C (row-major) order:

$$\text{offset}(i, j) = i \times (C \times \text{itemsize}) + j \times \text{itemsize}$$

This means row traversal is cache-friendly; column traversal skips C × itemsize bytes per step. Operations that traverse rows first (arr[i, :]) are therefore faster than column operations (arr[:, j]) unless the array is explicitly Fortran-ordered (np.asfortranarray).

PYTHON
import numpy as np

arr = np.ones((1000, 1000))
print(arr.strides)    # (8000, 8) — C order: rows are contiguous

arr_f = np.asfortranarray(arr)
print(arr_f.strides)  # (8, 8000) — F order: columns are contiguous

The practical rule: always sum/reduce along the last axis first for C-order arrays, and profile before assuming layout is the bottleneck.

Broadcasting

Broadcasting is NumPy's mechanism for applying operations to arrays with different (but compatible) shapes without allocating the virtual copies. The rule is:

  1. Align shapes on the right, padding with 1s on the left.
  2. Each dimension must be equal, or one of them must be 1.
  3. Dimensions of size 1 are stretched logically (not physically copied).

$$\text{shape } (M, 1) + \text{shape } (1, N) \to \text{shape } (M, N)$$

PYTHON
import numpy as np

row = np.arange(4).reshape(1, 4)   # shape (1, 4)
col = np.arange(3).reshape(3, 1)   # shape (3, 1)

result = row + col                 # shape (3, 4) — no data copied
print(result)
# [[0 1 2 3]
#  [1 2 3 4]
#  [2 3 4 5]]

Common Broadcasting Pattern: Standardisation

Standardising a matrix of features (subtract column mean, divide by column std) is a one-liner with broadcasting:

PYTHON
X = np.random.randn(1000, 20)          # 1000 samples, 20 features
X_std = (X - X.mean(axis=0)) / X.std(axis=0)  # mean/std shape (20,)
# Broadcasting: (1000,20) - (20,) → (20,) aligned on right → (1000,20)

Pitfall: Unnecessary Copies

Many NumPy operations return views (zero-copy): slicing, transpose, reshape (when contiguous). Others return copies: fancy indexing (arr[[0,2,4]]), boolean indexing, np.concatenate. Check with arr.base is None — if True, it owns its data (copy); if False, it is a view.

PYTHON
a = np.arange(12).reshape(3, 4)
b = a[:, 1:3]          # view — shares memory
b[0, 0] = 999
print(a[0, 1])         # 999 — a was mutated!

c = a[[0, 2], :]       # fancy index → copy
c[0, 0] = -1
print(a[0, 0])         # 0 — a is unchanged

Vectorising Conditional Logic with np.where and np.select

Avoid loops for element-wise conditional assignment:

PYTHON
import numpy as np

scores = np.array([45, 62, 78, 55, 91, 33])

# Binary condition
grades = np.where(scores >= 60, 'pass', 'fail')

# Multiple conditions — np.select
conditions  = [scores >= 90, scores >= 70, scores >= 60]
choices     = ['A', 'B', 'C']
letters     = np.select(conditions, choices, default='D')
print(letters)  # ['D' 'C' 'B' 'D' 'A' 'D']

Random Number Generation Best Practices

Use np.random.default_rng(seed) (the new Generator API, introduced NumPy 1.17) rather than np.random.seed(). The new API is statistically superior (PCG64 generator), thread-safe, and composable.

PYTHON
rng = np.random.default_rng(seed=42)
samples = rng.standard_normal(size=(5, 3))

Code Examples

Pairwise Euclidean distances via broadcasting

Classic broadcasting example: compute all pairwise distances in one expression.

PYTHON
import numpy as np

def pairwise_distances(X: np.ndarray) -> np.ndarray:
    """
    Compute the N×N matrix of pairwise L2 distances without scipy.
    X: shape (N, D)
    """
    # Expand dims to exploit broadcasting:
    # diff shape: (N, 1, D) - (1, N, D) -> (N, N, D)
    diff = X[:, np.newaxis, :] - X[np.newaxis, :, :]
    return np.sqrt((diff ** 2).sum(axis=-1))

rng = np.random.default_rng(0)
X = rng.standard_normal((5, 3))
D = pairwise_distances(X)
print(D.shape)      # (5, 5)
print(D[0, 1])      # same as np.linalg.norm(X[0] - X[1])
print(np.allclose(D, D.T))  # True — symmetric
Output
(5, 5)
3.0142...
True

Stride tricks for sliding-window features

Zero-copy rolling window using stride tricks — useful for time-series feature extraction.

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

def rolling_window(arr: np.ndarray, window: int) -> np.ndarray:
    """
    Return shape (N-window+1, window) view — zero-copy.
    WARNING: do not write to the returned array.
    """
    n = arr.shape[0]
    shape   = (n - window + 1, window)
    strides = (arr.strides[0], arr.strides[0])
    return as_strided(arr, shape=shape, strides=strides)

ts = np.arange(10, dtype=float)
windows = rolling_window(ts, 3)
roll_mean = windows.mean(axis=1)
print(roll_mean)   # [1. 2. 3. 4. 5. 6. 7. 8.]
Output
[1. 2. 3. 4. 5. 6. 7. 8.]

Structured arrays for heterogeneous tabular data

Structured arrays represent mixed-type records in contiguous memory, bridging NumPy and pandas.

PYTHON
import numpy as np

dtype = np.dtype([
    ('id',    np.int32),
    ('score', np.float32),
    ('label', 'U10'),
])

data = np.array(
    [(1, 0.92, 'positive'),
     (2, 0.34, 'negative'),
     (3, 0.71, 'positive')],
    dtype=dtype
)

print(data['score'])             # array([0.92, 0.34, 0.71], dtype=float32)
print(data[data['score'] > 0.5]) # rows where score > 0.5
Output
[0.92 0.34 0.71]
[(1, 0.92, 'positive') (3, 0.71, 'positive')]

49.3 Pandas Internals, Dtypes, and Copy vs. View Intermediate

The DataFrame as a Collection of Blocks

A pandas DataFrame is not a 2-D NumPy array — it is a dict-like container of Series objects, each backed by a 1-D (or 2-D consolidated) NumPy array. Internally, pandas groups columns of the same dtype into blocks to enable efficient columnar operations. A DataFrame with 10 float64 columns and 5 object (string) columns uses two blocks: one float64 block of shape (10, N) and one object block of shape (5, N).

This block structure has a critical implication: adding or removing a single column can trigger block consolidation, a copy of all data. For performance-critical code that iteratively adds columns, pre-allocate with pd.DataFrame(np.empty(...), columns=...) or build column dicts and call pd.DataFrame(d) once.

Dtypes and Memory

Pandas inherits NumPy dtypes and adds its own:

  • float64 / int64 — default numeric types; 8 bytes per element.
  • object — a Python object array; stores pointers (8 bytes each) to heap-allocated strings. Extremely memory-inefficient and slow for string operations.
  • category — encodes a low-cardinality column as integer codes + a lookup table. Reduces memory from O(N × stringlength) to O(N × 1–2 bytes) for the codes plus a small vocabulary.
  • Int64 (capital I) — nullable integer (pandas ExtensionArray), distinct from NumPy int64 which cannot hold NaN.
  • datetime64[ns] / DatetimeTZDtype — nanosecond-precision timestamps.

PYTHON
import pandas as pd, numpy as np

N = 1_000_000
df = pd.DataFrame({
    'city': np.random.choice(['NYC', 'LA', 'CHI', 'HOU'], N),
    'value': np.random.randn(N),
})

print(df.memory_usage(deep=True))  # deep=True counts string bytes
# city:  ~64 MB as object

df['city'] = df['city'].astype('category')
print(df.memory_usage(deep=True))
# city:  ~1 MB as category — 64× reduction!

Copy vs. View Semantics and SettingWithCopyWarning

This is the most common source of silent bugs in pandas. The rule is: chained indexing can return either a copy or a view depending on the operation, and you cannot reliably predict which without inspecting internals.

PYTHON
# Dangerous — may not modify df
df[df['value'] > 0]['value'] = -1     # SettingWithCopyWarning

# Safe — always modifies df in-place
df.loc[df['value'] > 0, 'value'] = -1

The root cause: df[df[&#039;value&#039;] &gt; 0] uses <strong>getitem</strong>, which may return a copy (especially with boolean indexing). Assigning to [&#039;value&#039;] on that potential copy has no effect on df. Use .loc[row</em>mask, column<em>label] for all conditional assignments.

The .copy() Contract

When you intend to work on an independent slice, be explicit:

PYTHON
train = df[df['split'] == 'train'].copy()   # explicit deep copy
train['normalised'] = train['value'] / train['value'].max()

Pandas 2.0 introduced Copy-on-Write (CoW) semantics (pd.options.mode.copy</em>on<em>write = True), which will become the default in pandas 3.0. Under CoW, every indexing operation that could mutate data returns a lazy copy-on-first-write, eliminating the ambiguity entirely.

Efficient Column Creation and eval/query

For arithmetic on large DataFrames, pd.eval and df.query bypass Python's per-object dispatch and operate directly on underlying arrays using the numexpr engine:

PYTHON
N = 5_000_000
df = pd.DataFrame({'a': np.random.randn(N), 'b': np.random.randn(N)})

# Standard Python expression — creates temporaries
result_slow = df['a'] ** 2 + df['b'] ** 2

# pd.eval — single pass, no intermediary arrays
result_fast = pd.eval('a**2 + b**2', local_dict={'a': df['a'], 'b': df['b']})
# Or inline:
df['c'] = df.eval('a**2 + b**2')

df.query similarly filters without materialising the boolean mask as a separate Series:

PYTHON
subset = df.query('a > 0 and b < 0.5')

The Index Object

Every pandas Series and DataFrame carries an Index. Index alignment is what makes df1 + df2 correct even when the rows are in different orders — pandas aligns on the index before operating. This is powerful but can be a source of bugs when indices disagree:

PYTHON
s1 = pd.Series([1, 2, 3], index=['a', 'b', 'c'])
s2 = pd.Series([10, 20, 30], index=['b', 'c', 'd'])
print(s1 + s2)
# a     NaN
# b    22.0
# c    33.0
# d     NaN

For high-performance operations where alignment is already guaranteed, use df[&#039;col&#039;].values (returns the NumPy array) or df[&#039;col&#039;].to</em>numpy() to bypass index overhead.

Code Examples

Auditing DataFrame memory usage by dtype

Shows how changing dtypes (category, int8) dramatically reduces memory.

PYTHON
import pandas as pd
import numpy as np

def memory_report(df: pd.DataFrame) -> pd.Series:
    """Return per-column memory in MB, sorted descending."""
    usage = df.memory_usage(deep=True) / 1024**2
    return usage.sort_values(ascending=False)

np.random.seed(0)
N = 500_000
df = pd.DataFrame({
    'user_id':  np.random.randint(0, 100_000, N).astype(str),
    'country':  np.random.choice(['US','UK','DE','FR','JP'], N),
    'revenue':  np.random.exponential(50, N),
    'quantity': np.random.randint(1, 100, N),
})

print("Before optimisation:")
print(memory_report(df))

df['user_id']  = df['user_id'].astype('category')
df['country']  = df['country'].astype('category')
df['quantity'] = df['quantity'].astype(np.int8)

print("\nAfter optimisation:")
print(memory_report(df))
Output
Before optimisation:
user_id     28.61
country      3.81
revenue      3.81
quantity     3.81
...
After optimisation:
revenue      3.81
user_id      0.49
country      0.04
quantity     0.48
...

Copy-on-Write demonstration with pandas 2.0

Demonstrates pandas 2.0 Copy-on-Write: mutations to a slice do not affect the original.

PYTHON
import pandas as pd

# Enable Copy-on-Write (default in pandas 3.0)
pd.options.mode.copy_on_write = True

df = pd.DataFrame({'x': [1, 2, 3, 4, 5]})

# Under CoW, this slice is a lazy copy
slice_ = df[df['x'] > 2]

# Mutating slice_ triggers the copy — df is never affected
slice_.loc[slice_.index, 'x'] = -999

print("Original df:")
print(df)
print("\nModified slice_:")
print(slice_)
Output
Original df:
   x
0  1
1  2
2  3
3  4
4  5

Modified slice_:
   x
2 -999
3 -999
4 -999

49.4 GroupBy, Aggregation, and Performant Apply Patterns Intermediate

The GroupBy Machinery

Pandas groupby follows a split–apply–combine pattern:

  1. Split: partition the DataFrame into groups by one or more key columns.
  2. Apply: run a function on each group independently.
  3. Combine: concatenate results into a new DataFrame.

The split step builds a hash table mapping each unique key to the row indices that belong to it — this is done once and cached on the GroupBy object, so you can chain multiple aggregations without repeating the hashing.

PYTHON
import pandas as pd, numpy as np

rng = np.random.default_rng(42)
N = 100_000
df = pd.DataFrame({
    'region':  rng.choice(['North','South','East','West'], N),
    'product': rng.choice(['A','B','C'], N),
    'sales':   rng.exponential(200, N),
    'units':   rng.integers(1, 50, N),
})

# Named aggregation syntax (pandas >= 0.25) — explicit, readable
summary = df.groupby(['region', 'product']).agg(
    total_sales  = ('sales', 'sum'),
    mean_sales   = ('sales', 'mean'),
    total_units  = ('units', 'sum'),
    n_orders     = ('sales', 'count'),
)
print(summary.head())

transform vs. agg vs. apply

These three methods differ in what they return and when to use them:

  • agg — reduces each group to a scalar (or a few scalars). Returns a DataFrame with one row per group. Use for summary statistics.
  • transform — maps each group to an array of the same shape as the group, aligned back to the original index. Use for adding group-level features (e.g., group mean) back to the original DataFrame.
  • apply — calls an arbitrary Python function on each group DataFrame and stitches the results. Flexible but slow because it cannot be vectorised by pandas.

PYTHON
# transform: add group mean as a column (same length as df)
df['region_mean_sales'] = df.groupby('region')['sales'].transform('mean')

# This is equivalent to — but much faster than —
group_means = df.groupby('region')['sales'].mean()
df['region_mean_sales_slow'] = df['region'].map(group_means)

When apply Is Acceptable

Use apply only when no vectorised method exists: complex multi-column operations, functions returning DataFrames of variable size, or third-party functions that operate on a whole sub-DataFrame.

PYTHON
def weighted_avg(g):
    return (g['sales'] * g['units']).sum() / g['units'].sum()

result = df.groupby('region').apply(weighted_avg)

For the common case of a simple function on a single column, always prefer a named aggregation:

PYTHON
# Slow — apply on Series
df.groupby('region')['sales'].apply(np.mean)

# Fast — built-in aggregation
df.groupby('region')['sales'].mean()

Performance Benchmark

On a 100 k-row DataFrame computing group standard deviation:

  • Built-in .std(): ~1 ms
  • transform(np.std): ~40 ms
  • apply(lambda g: g.std()): ~80 ms

The gap grows with group count and DataFrame size.

pipe for Composable Transformations

For multi-step pipelines, .pipe() lets you chain custom functions while keeping data flow readable:

PYTHON
def clip_outliers(df, col, lo=0.01, hi=0.99):
    bounds = df[col].quantile([lo, hi])
    df = df.copy()
    df[col] = df[col].clip(*bounds.values)
    return df

def add_log_feature(df, col):
    df = df.copy()
    df[f'log_{col}'] = np.log1p(df[col])
    return df

clean = (
    df
    .pipe(clip_outliers, 'sales')
    .pipe(add_log_feature, 'sales')
    .pipe(lambda d: d.assign(sales_rank=d['sales'].rank(pct=True)))
)

Window Functions

Rolling, expanding, and exponential-weighted windows work on the GroupBy object too:

PYTHON
# 7-day rolling mean within each region, on a time-indexed DataFrame
df_ts = df.set_index('date').sort_index()
df_ts['roll7'] = (
    df_ts.groupby('region')['sales']
         .transform(lambda s: s.rolling('7D').mean())
)

# Expanding (cumulative) mean
df_ts['cum_mean'] = (
    df_ts.groupby('region')['sales']
         .transform(lambda s: s.expanding().mean())
)

Expansion via .expanding() is useful for avoiding look-ahead leakage in time-series cross-validation: the expanding window only uses information available up to time $t$.

Code Examples

Named aggregation with multiple custom functions

Named aggregation with built-in, lambda, and scipy functions in one pass.

PYTHON
import pandas as pd
import numpy as np

rng = np.random.default_rng(7)
df = pd.DataFrame({
    'dept':   rng.choice(['Eng','Sales','HR'], 2000),
    'salary': rng.normal(80_000, 15_000, 2000).round(2),
    'yoe':    rng.integers(0, 30, 2000),
})

from scipy import stats as sp_stats

result = df.groupby('dept').agg(
    headcount    = ('salary', 'count'),
    mean_salary  = ('salary', 'mean'),
    median_salary= ('salary', 'median'),
    p90_salary   = ('salary', lambda s: s.quantile(0.9)),
    salary_iqr   = ('salary', sp_stats.iqr),
    mean_yoe     = ('yoe',    'mean'),
).round(2)

print(result)
Output
       headcount  mean_salary  median_salary  p90_salary  salary_iqr  mean_yoe
dept
Eng          675     80121.35       80034.12    99872.45    20341.12     14.53
HR           667     79834.21       79512.34    98432.11    19872.32     14.21
Sales        658     80234.56       80123.45   100234.56    20123.45     14.87

Group-normalisation using transform for ML feature engineering

Within-group standardisation with transform — a key preprocessing step for grouped data.

PYTHON
import pandas as pd
import numpy as np

rng = np.random.default_rng(3)
df = pd.DataFrame({
    'category': rng.choice(['A','B','C'], 500),
    'value':    rng.exponential(10, 500),
})

# Z-score within each category group (no data leakage between groups)
df['z_within'] = df.groupby('category')['value'].transform(
    lambda s: (s - s.mean()) / s.std()
)

# Verify: each group should have mean≈0, std≈1
check = df.groupby('category')['z_within'].agg(['mean','std']).round(6)
print(check)
Output
              mean   std
category
A        -0.0  1.0
B        -0.0  1.0
C        -0.0  1.0

49.5 Dates, Times, and Time Zones in Python and Pandas Intermediate

The Date/Time Stack

Python's date/time ecosystem has several layers, and mixing them carelessly is a reliable source of bugs:

  • datetime.datetime — standard library, microsecond precision, optional tzinfo.
  • dateutil.parser.parse — lenient parsing of human-readable strings.
  • numpy.datetime64 — fixed-unit (e.g. ns) integer timestamp, no timezone.
  • pandas.Timestamp — wraps datetime64[ns], adds timezone support via pytz or zoneinfo.
  • pandas.DatetimeTZDtype — timezone-aware pandas column type.

The key rule: always work with timezone-aware datetimes in production systems. Naive datetimes (no tzinfo) are ambiguous during daylight-saving transitions and across system moves.

Parsing and Constructing Timestamps

PYTHON
import pandas as pd
from datetime import datetime, timezone
import pytz

# Naive — dangerous in production
t_naive = datetime(2024, 3, 10, 2, 30)   # ambiguous: DST gap in US/Eastern

# Aware — always unambiguous
ut = pytz.UTC
t_utc = datetime(2024, 3, 10, 7, 30, tzinfo=timezone.utc)
t_eastern = t_utc.astimezone(pytz.timezone('US/Eastern'))
print(t_eastern)  # 2024-03-10 02:30:00-05:00

# Pandas parsing
ts = pd.Timestamp('2024-03-10 07:30:00', tz='UTC')
print(ts.tz_convert('US/Eastern'))

DatetimeIndex and Resampling

When a DatetimeIndex backs a DataFrame, pandas exposes powerful time-series operations:

PYTHON
import pandas as pd, numpy as np

rng = pd.date_range('2023-01-01', periods=365*24, freq='h', tz='UTC')
ts = pd.Series(np.random.randn(len(rng)), index=rng, name='temp')

# Resample to daily mean
daily = ts.resample('D').mean()

# Resample to weekly open/high/low/close
weekly_ohlc = ts.resample('W').ohlc()

# Partial string indexing
jan = ts['2023-01']           # all of January
q1  = ts['2023-01':'2023-03'] # Q1

Pitfall: freq Inference After Resampling

After some operations (join, concat, boolean index), the DatetimeIndex.freq attribute becomes None, silently breaking resample. Always call pd.infer_freq(idx) or set freq explicitly.

Vectorised Date Arithmetic with .dt Accessor

The .dt accessor exposes datetime components and arithmetic on a pandas Series without looping:

PYTHON
df = pd.DataFrame({
    'event_time': pd.to_datetime([
        '2024-01-15 09:32:00',
        '2024-02-28 14:15:00',
        '2024-12-31 23:59:59',
    ])
})

df['year']        = df['event_time'].dt.year
df['day_of_week'] = df['event_time'].dt.day_name()
df['quarter']     = df['event_time'].dt.quarter
df['hour']        = df['event_time'].dt.hour
df['is_weekend']  = df['event_time'].dt.dayofweek >= 5

# Time difference
df['days_to_eoy'] = (
    pd.Timestamp('2024-12-31') - df['event_time']
).dt.days

Handling Missing and Irregular Timestamps

Real-world time series have gaps, duplicates, and timezone inconsistencies:

PYTHON
# Forward-fill gaps up to 2 hours, then leave NaN
ts_clean = ts.resample('h').mean().ffill(limit=2)

# Remove duplicate timestamps (keep last)
ts_dedup = ts[~ts.index.duplicated(keep='last')]

# Timezone-naive data from a database: localise then convert
raw = pd.to_datetime(['2024-06-15 12:00:00', '2024-06-15 13:00:00'])
aware = pd.DatetimeIndex(raw).tz_localize('US/Eastern').tz_convert('UTC')

Period vs. Timestamp

For monthly or quarterly data where the interval matters more than the instant, use pd.Period:

PYTHON
p = pd.Period('2024-Q1', freq='Q')
print(p.start_time)   # 2024-01-01
print(p.end_time)     # 2024-03-31

# Convert Series of Timestamps to Periods
quarters = pd.to_datetime(['2024-02-01','2024-05-15']).to_period('Q')
print(quarters)       # PeriodIndex(['2024Q1', '2024Q2'], dtype='period[Q-DEC]')

Code Examples

Robust timezone-aware timestamp parsing from mixed-format strings

Parse mixed-format, mixed-timezone strings into a uniform UTC DatetimeIndex.

PYTHON
import pandas as pd
from dateutil import parser as dparser

# Mixed formats and timezone offsets — real-world log data
raw_times = [
    '2024-03-15T09:22:11Z',
    '15/03/2024 09:22:11 +0000',
    'March 15, 2024 9:22 AM UTC',
    '2024-03-15 09:22:11+00:00',
]

# dateutil handles the variance; pd.to_datetime handles the volume
parsed = pd.to_datetime(
    [dparser.parse(s) for s in raw_times],
    utc=True        # normalise all to UTC
)
print(parsed)
print(parsed.dtype)
Output
DatetimeIndex(['2024-03-15 09:22:11+00:00', '2024-03-15 09:22:11+00:00',
               '2024-03-15 09:22:11+00:00', '2024-03-15 09:22:11+00:00'],
              dtype='datetime64[ns, UTC]', freq=None)
datetime64[ns, UTC]

Business-day calendars and custom offsets

Business-day and quarter-end date arithmetic common in financial data work.

PYTHON
import pandas as pd
from pandas.tseries.offsets import BDay, BMonthEnd

# Next 5 business days from today
start = pd.Timestamp('2024-06-03')   # Monday
bdays = pd.date_range(start, periods=5, freq=BDay())
print(bdays)

# Business-month-end dates for H1 2024
bme = pd.date_range('2024-01-01', '2024-06-30', freq=BMonthEnd())
print(bme)

# Days until next quarter-end (vectorised)
dates = pd.date_range('2024-01-01', periods=10, freq='W')
qtr_ends = dates.to_period('Q').end_time
days_until_qe = (qtr_ends - dates).days
print(days_until_qe.values)
Output
DatetimeIndex(['2024-06-03', '2024-06-04', '2024-06-05', '2024-06-06', '2024-06-07'], ...)
DatetimeIndex(['2024-01-31', '2024-02-29', '2024-03-29', '2024-04-30', '2024-05-31', '2024-06-28'], ...)
[90 83 76 69 62 55 48 41 34 27]

49.6 Writing Clean, Testable, and Reproducible Analysis Code Advanced

The Case for Engineered Analysis Code

Data science notebooks are notorious for hidden state, untestable monoliths, and results that cannot be reproduced six months later. The antidote is not abandoning notebooks — it is organising the logic into importable, tested Python modules while using notebooks only for orchestration and visualisation. This section covers the practices that separate reproducible research from one-off scripts.

Functions, Type Hints, and Docstrings

Every non-trivial transformation should be a named function with a type annotation and a docstring that states the pre- and post-conditions:

PYTHON
from __future__ import annotations
import pandas as pd
import numpy as np

def winsorise(series: pd.Series, lo: float = 0.01, hi: float = 0.99) -> pd.Series:
    """
    Clip values to the [lo, hi] quantile range.

    Parameters
    ----------
    series : numeric Series
    lo, hi : quantile bounds (0 < lo < hi < 1)

    Returns
    -------
    pd.Series — same index as input, no NaN added or removed.

    Raises
    ------
    ValueError if lo >= hi.
    """
    if not 0 < lo < hi < 1:
        raise ValueError(f"Require 0 < lo < hi < 1, got lo={lo}, hi={hi}")
    bounds = series.quantile([lo, hi])
    return series.clip(*bounds.values)

Type hints are enforced by tools like mypy (static) and beartype (runtime), and they make function contracts explicit without running the code.

Reproducible Random Seeds

Every stochastic step must be seeded explicitly. Do not use global state:

PYTHON
import numpy as np
from sklearn.model_selection import train_test_split

SEED = 42
rng  = np.random.default_rng(SEED)          # NumPy
train, test = train_test_split(df, random_state=SEED)  # sklearn

For end-to-end reproducibility, record the seed in your output artifacts and log it at the start of every pipeline run.

Testing Analysis Code with pytest

Data transformation functions should have unit tests. pytest's parametrize decorator makes it easy to test edge cases systematically:

PYTHON
# tests/test_features.py
import pytest
import pandas as pd
import numpy as np
from myproject.features import winsorise

@pytest.fixture
def sample_series():
    rng = np.random.default_rng(0)
    return pd.Series(rng.normal(0, 1, 1000))

def test_winsorise_bounds(sample_series):
    result = winsorise(sample_series, 0.05, 0.95)
    lo, hi = sample_series.quantile([0.05, 0.95])
    assert result.min() >= lo - 1e-9
    assert result.max() <= hi + 1e-9

def test_winsorise_no_nan_added(sample_series):
    result = winsorise(sample_series)
    assert result.isna().sum() == sample_series.isna().sum()

@pytest.mark.parametrize('lo,hi', [(0.5, 0.5), (0.9, 0.1), (-0.1, 0.9)])
def test_winsorise_invalid_args(sample_series, lo, hi):
    with pytest.raises(ValueError):
        winsorise(sample_series, lo, hi)

Run with pytest tests/ -v. Integration tests can load a small fixture CSV and assert the full pipeline output shape and key statistics.

Performance Profiling Workflow

Before optimising, measure. The recommended workflow:

  1. %timeit in Jupyter for micro-benchmarks.
  2. cProfile + pstats for function-level hotspot detection.
  3. line<em>profiler (%lprun) for line-level timing in hot functions.
  4. memory</em>profiler (%memit, %mprun) for allocation tracking.

PYTHON
import cProfile, pstats, io

def profile(func, *args, **kwargs):
    pr = cProfile.Profile()
    pr.enable()
    result = func(*args, **kwargs)
    pr.disable()
    s = io.StringIO()
    ps = pstats.Stats(pr, stream=s).sort_stats('cumulative')
    ps.print_stats(15)   # top 15 functions
    print(s.getvalue())
    return result

Common Performance Anti-Patterns

  • Iterating over DataFrame rows with iterrows() or itertuples(). Vectorise or use apply at worst.
  • Repeated column assignment in a loop triggers block consolidation. Build a dict of new columns, then call pd.concat or df.assign(<strong>new<em>cols) once.
  • String operations on object dtype are Python-speed. Convert to category or use pd.StringDtype backed by Arrow.
  • pd.merge on unsorted data without pre-sorting — for very large merges, sorting keys first and using how=&#039;inner&#039; with sort=False is faster.
  • Not using engine=&#039;pyarrow&#039;** in pd.read</em>csv / pd.read_parquet (pandas 2.0+). Arrow-backed I/O is 3–10× faster for typical CSVs.

PYTHON
# Fast CSV ingestion with Arrow
df = pd.read_csv('large_file.csv', engine='pyarrow', dtype_backend='pyarrow')

# Arrow-backed string column: 5-10× faster .str operations
df['name'] = df['name'].astype('string[pyarrow]')

Structuring a Project

A minimal layout for a reproducible data science project:

project/
  data/          # raw data (gitignored)
  notebooks/     # exploration only; import from src/
  src/
    myproject/
      __init__.py
      features.py    # transformations
      models.py      # training/inference
      io.py          # data loading/saving
  tests/
    test_features.py
  pyproject.toml

Notebooks import from src/myproject so the same logic is tested and reused, never duplicated.

Code Examples

Property-based testing of a data transformation with Hypothesis

Property-based testing with Hypothesis: verify invariants over hundreds of randomly generated DataFrames.

PYTHON
# pip install hypothesis
import pytest
import pandas as pd
import numpy as np
from hypothesis import given, settings
from hypothesis import strategies as st
from hypothesis.extra.pandas import column, data_frames

def clip_positive(s: pd.Series) -> pd.Series:
    """Clip series so all values are >= 0."""
    return s.clip(lower=0)

@given(
    data_frames([
        column('x', dtype=float,
               elements=st.floats(min_value=-1e6, max_value=1e6,
                                  allow_nan=False, allow_infinity=False))
    ])
)
@settings(max_examples=500)
def test_clip_positive_property(df):
    result = clip_positive(df['x'])
    assert (result >= 0).all()
    assert len(result) == len(df)
    assert result.index.equals(df.index)

Profiling a pandas pipeline with cProfile

Shell workflow for profiling a full pipeline and surfacing the top time consumers.

BASH
# Run a module under cProfile and display top 20 cumulative-time functions
python -m cProfile -o profile.stats my_pipeline.py

# Inspect in Python
python - <<'EOF'
import pstats
p = pstats.Stats('profile.stats')
p.strip_dirs().sort_stats('cumulative').print_stats(20)
EOF

Deterministic pipeline with seeds, logging, and artifact hash

Pipeline template with explicit seeding, structured logging, and SHA-256 output hashing for provenance.

PYTHON
import hashlib, json, logging, numpy as np, pandas as pd
from pathlib import Path

logging.basicConfig(level=logging.INFO, format='%(asctime)s %(message)s')
log = logging.getLogger(__name__)

SEED = 42

def run_pipeline(input_path: str, output_path: str, seed: int = SEED) -> dict:
    rng = np.random.default_rng(seed)
    log.info("Running pipeline with seed=%d", seed)

    df = pd.read_csv(input_path)
    df['noise'] = rng.standard_normal(len(df))   # reproducible noise

    out = Path(output_path)
    df.to_parquet(out, index=False)

    # Record artifact hash for provenance
    sha256 = hashlib.sha256(out.read_bytes()).hexdigest()
    manifest = {'input': input_path, 'output': output_path,
                'seed': seed, 'rows': len(df), 'sha256': sha256}
    log.info("Manifest: %s", json.dumps(manifest))
    return manifest