Chapter 4: Data Structures - Organizing Scientific Data

Learning Objectives¶

By the end of this chapter, you will be able to:

• (1) Choose optimal data structures (list, tuple, dict, set) based on algorithmic requirements and O(n) vs O(1) performance constraints.

• (2) Predict operation complexity using Big-O notation for common operations like search, insertion, and deletion.

• (3) Understand Python’s reference model and memory layout explaining why NumPy arrays are 100× faster than lists.

• (4) Implement defensive copying strategies using copy.copy() and copy.deepcopy() to prevent aliasing bugs.

• (5) Profile memory usage with sys.getsizeof() and optimize data structure choices for particle simulations.

• (6) Master dictionary operations for O(1) lookups using hash tables for caching expensive calculations.

• (7) Apply set operations (union, intersection, difference) for particle tracking and domain decomposition.

• (8) Debug three common bugs: mutable defaults, aliasing during grid updates, and iteration modification.

• (9) Design hybrid data structures combining lists, dicts, and sets for real physics simulations.

Prerequisites Check¶

Chapter Overview¶

Imagine you’re simulating the interactions between a million particles - whether they’re dark matter particles in a galaxy, molecules in a gas, or nodes in a network. Each timestep, you need to find which particles are close enough to interact strongly. With the wrong data structure, this neighbor search could take hours per timestep. With the right one - a spatial hash table or tree structure - it takes seconds. That’s the difference between a simulation finishing in a day versus running for months. This chapter teaches you to make these critical choices that determine whether your code scales to research problems or remains stuck with toy models.

This chapter transforms you from someone who stores data to someone who orchestrates it strategically for scientific computing. You’ll discover not just that dictionary lookups are fast, but why they’re fast - through hash functions that turn particle positions into array indices. You’ll understand when they might fail - like when hash collisions cluster your data. And you’ll learn how to verify performance yourself - because in computational science, measurement beats assumption every time.

These concepts directly prepare you for the numerical computing ahead. The memory layout discussions explain why NumPy arrays can be 100× more efficient than Python lists for vector operations. The immutability concepts prepare you for functional programming paradigms used in modern frameworks like JAX, where immutable operations enable automatic differentiation. The performance profiling skills will help you identify bottlenecks whether you’re solving differential equations or analyzing experimental data. By chapter’s end, you’ll think about data organization like a computational scientist, architecting for performance from the start.

4.1 What Is a Data Structure?¶

A data structure is fundamentally about organizing information to match your access patterns. Think about an N-body simulation where particles interact through forces - these particles could represent stars in a star cluster, molecules in a gas, or charges in a plasma. You could store particles in order of creation (like a list), organize them by spatial region for fast neighbor finding (like a dictionary of cells), or track unique particle IDs (like a set). Each choice profoundly affects your simulation’s performance - the difference between O(n²) all-pairs checks and O(n log n) tree-based algorithms.

Quick Preview: Python’s Core Data Structures¶

*Python 3.7+ dicts maintain insertion order

Building Intuition: Measuring Speed Empirically¶

Ready to discover something that will transform how you think about organizing data in your simulations? We’re going to measure the dramatic performance difference between different data structures. This empirical approach - measure first, understand why, then optimize - is exactly how computational scientists approach performance optimization.

# Let's discover why data structure choice matters for simulations!
import time
import random

# Create test data representing particle IDs in a simulation
sizes = [100, 1000, 10000, 100000]

print("Searching for a particle NOT in the collection (worst case):")
print("=" * 60)

for n in sizes:
    # Simulate particle IDs
    particle_list = list(range(n))
    particle_set = set(range(n))
    
    # Search for a particle that escaped our simulation boundary
    escaped_particle = -1
    
    # Time list search - what you might naturally try first
    start = time.perf_counter()
    found = escaped_particle in particle_list
    list_time = time.perf_counter() - start
    
    # Time set search - the optimized approach
    start = time.perf_counter()
    found = escaped_particle in particle_set
    set_time = time.perf_counter() - start
    
    print(f"Size {n:6d}: List: {list_time*1e6:8.2f} μs, "
          f"Set: {set_time*1e6:8.2f} μs")
    if list_time > 0:
        print(f"           Set is {list_time/set_time:,.0f}× faster!")

print("\nFor collision detection between 100,000 particles:")
print("List approach: hours per timestep")
print("Set approach: milliseconds per timestep!")

This performance difference isn’t just academic - it determines whether your galaxy simulation can evolve for billions of years or gets stuck after a few million. The secret behind this magic is the hash table, which you’re about to understand completely.

Understanding Big-O Notation¶

Now that you’ve witnessed this dramatic performance difference empirically, let’s understand the mathematical pattern behind it. Big-O notation describes how an algorithm’s runtime scales with input size - it’s the language computational scientists use to discuss whether an algorithm is feasible for large-scale simulations.

import matplotlib.pyplot as plt
import math

# Visualize how different algorithms scale for physics simulations
n = range(1, 101)
constant = [1 for _ in n]
logarithmic = [math.log(x) for x in n]
linear = list(n)
quadratic = [x**2 for x in n]

plt.figure(figsize=(10, 6))
plt.plot(n, constant, label='O(1) - Hash table lookup', linewidth=2)
plt.plot(n, logarithmic, label='O(log n) - Tree-based methods', linewidth=2)
plt.plot(n, linear, label='O(n) - Direct summation', linewidth=2)
plt.plot(n, quadratic, label='O(n²) - All-pairs (naive)', linewidth=2)

plt.xlim(0, 100)
plt.ylim(0, 500)
plt.xlabel('Number of Particles (thousands)')
plt.ylabel('Time (arbitrary units)')
plt.title('Algorithm Scaling: Why Data Structures Matter')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

print("For a million-particle simulation:")
print(f"  O(1): 1 operation - spatial hash lookup")
print(f"  O(log n): {math.log2(1e6):.0f} operations - tree traversal")
print(f"  O(n): 1,000,000 operations - single particle force sum")
print(f"  O(n²): 1,000,000,000,000 operations - all pairs (impossible!)")

print("\nIn astronomy: This enables galaxy simulations with billions of stars")
print("In chemistry: This allows protein folding with millions of atoms")
print("In climate: This permits global models with millions of grid cells")

This graph reveals why algorithmic complexity matters for computational physics. The naive O(n²) approach that checks all particle pairs becomes impossible beyond a few thousand particles. But with smart data structures like trees (O(n log n)) or spatial hashing (O(n)), we can simulate entire galaxies, proteins, or climate systems!

4.2 Lists: Python’s Workhorse Sequence¶

Lists are Python’s most versatile data structure - perfect for particle arrays, time series data, or any ordered collection. But here’s something crucial for scientific computing: Python lists don’t store your numbers directly! Understanding this explains why NumPy arrays are so much faster for numerical work.

How Lists Really Work in Memory¶

Let’s explore how Python actually stores your simulation data in memory. This knowledge will help you understand when to use Python lists versus when you need NumPy arrays.

# Stage 1: Basic memory inspection (10 lines)
import sys

# Position of one particle in 3D space
position = [1.5e10, 2.3e10, 0.8e10]  # cm (typical scale for solar system)

# The list container itself
list_size = sys.getsizeof(position)
print(f"List container: {list_size} bytes")

# Each float is a full Python object!
element_sizes = [sys.getsizeof(x) for x in position]
print(f"Each coordinate: {element_sizes[0]} bytes")

# Stage 2: Calculate total overhead (10 lines)
# Total memory footprint
total = list_size + sum(element_sizes)
print(f"Total: {total} bytes for 3 floats!")
print(f"That's {total/24:.1f}× more than raw floats would need!")

print("\nWhy so much memory?")
print("- List stores pointers, not values")
print("- Each float is a full object with type info")
print("- Python must track reference counts")
print("\nThis is why NumPy arrays (Chapter 7) store raw values contiguously!")

Here’s what’s actually happening in memory when you store particle positions:

Python List of Positions:          Objects in Memory:
┌─────────────────┐               ┌──────────────────┐
│ size: 3         │               │ float: 1.5e10    │
│ capacity: 4     │          ┌───>│ type info        │
├─────────────────┤          │    │ ref count        │
│ pointer ────────┼──────────┘    └──────────────────┘
│ pointer ────────┼──────────────>┌──────────────────┐
│ pointer ────────┼──────┐        │ float: 2.3e10    │
│ (unused)        │    │          └──────────────────┘
└─────────────────┘    └────────>┌──────────────────┐
                                  │ float: 0.8e10    │
                                  └──────────────────┘

NumPy Array (preview):
┌────┬────┬────┐
│1.5 │2.3 │0.8 │  <- Raw values, no pointers!
└────┴────┴────┘

List Operations: Performance for Simulations¶

Different list operations have vastly different costs - critical knowledge when processing particle data or building spatial data structures.

import time

# Simulate a growing particle system
particles = list(range(100_000))

# Adding particles at the END (common pattern)
start = time.perf_counter()
particles.append(100_000)  # New particle enters domain
particles.pop()            # Remove for fair comparison
end_time = time.perf_counter() - start

# Adding particles at the BEGINNING (avoid this!)
start = time.perf_counter()
particles.insert(0, -1)    # Insert at front
particles.pop(0)           # Remove from front
begin_time = time.perf_counter() - start

print(f"Adding particle at END:   {end_time*1e6:.2f} microseconds")
print(f"Adding particle at START: {begin_time*1e6:.2f} microseconds")
print(f"\nBeginning is {begin_time/end_time:.0f}× slower!")

print("\nWhy? Inserting at the beginning shifts ALL particles in memory!")
print("For particle systems, use collections.deque if you need")
print("fast operations at both ends (e.g., boundary conditions).")

List Growth Strategy: Understanding Dynamic Arrays¶

Watch how Python manages memory as your particle system grows. This pattern appears in many languages and understanding it helps you write efficient simulation codes.

# Stage 1: Basic growth observation (12 lines)
import sys

particles = []
print("Watch Python's list growth strategy:")
print("(Critical for understanding simulation performance)")
print("Length → Capacity (overallocation)")
print("-" * 40)

for i in range(10):
    old_size = sys.getsizeof(particles)
    particles.append(i)
    new_size = sys.getsizeof(particles)
    
    if new_size != old_size:
        print(f"{len(particles):4d} → larger allocation")

# Stage 2: Calculate exact capacities (15 lines)
particles = []
previous_size = 0

for i in range(20):
    old_size = sys.getsizeof(particles)
    particles.append(i)
    new_size = sys.getsizeof(particles)
    
    if new_size != old_size:
        # Calculate capacity from size change
        capacity = (new_size - sys.getsizeof([])) // 8
        actual_length = len(particles)
        overalloc = ((capacity - actual_length) / actual_length * 100 
                     if actual_length > 0 else 0)
        print(f"Length {actual_length:4d} → Capacity {capacity:4d} "
              f"({overalloc:5.1f}% extra)")

This growth pattern is why appending to lists is “amortized O(1)” - usually fast, occasionally slow when reallocation happens. For time-critical simulations, pre-allocate your arrays when particle count is known!

4.3 Tuples: The Power of Immutability¶

What if you need to guarantee that initial conditions or physical constants won’t change during your simulation? Enter tuples - Python’s immutable sequences. This isn’t a limitation - it’s protection against an entire category of bugs that plague scientific codes!

Understanding Immutability in Physics Simulations¶

Let’s see how immutability protects your simulations and enables powerful optimizations:

# Physical constants should NEVER change during simulation
class PhysicsConstants:
    """Demonstrating safe vs unsafe constant storage."""
    
    # UNSAFE: Mutable list (can be accidentally modified)
    unsafe_constants = [
        6.674e-8,   # G (cm³/g/s²)
        2.998e10,   # c (cm/s)
        1.381e-16,  # k_B (erg/K)
    ]
    
    # SAFE: Immutable tuple (modification attempts fail loudly)
    safe_constants = (
        6.674e-8,   # G (cm³/g/s²)
        2.998e10,   # c (cm/s)
        1.381e-16,  # k_B (erg/K)
    )

# Demonstration of the danger
def buggy_calculation(constants):
    """This function has a typo that modifies constants!"""
    if isinstance(constants, list):
        # Oops! Used = instead of == in a complex calculation
        constants[0] = constants[0] * 1e10  # BUG: Modifies G!
        return "Calculated (with hidden corruption)"
    else:
        # With tuple, this would raise an error immediately
        try:
            constants[0] = constants[0] * 1e10
        except TypeError as e:
            return f"Caught bug immediately: {e}"

# Test with both
print("With list:", buggy_calculation(PhysicsConstants.unsafe_constants.copy()))
print("With tuple:", buggy_calculation(PhysicsConstants.safe_constants))

print("\nImmutability turns subtle runtime corruption into immediate, obvious errors!")

Tuples as Dictionary Keys: Caching Expensive Calculations¶

Here’s where immutability becomes a superpower for computational physics - tuples can be dictionary keys, enabling powerful memoization patterns for expensive calculations:

# Stage 1: Basic caching concept (15 lines)
# Cache expensive physics calculations
cache = {}
calculation_count = 0

def gravitational_force(pos1, pos2, m1, m2):
    """Calculate gravitational force with caching."""
    global calculation_count
    
    # Create cache key from positions (must be tuples!)
    cache_key = (pos1, pos2, m1, m2)
    
    if cache_key in cache:
        print(f"  Cache hit! Avoided expensive calculation")
        return cache[cache_key]
    
    calculation_count += 1
    print(f"  Computing force (calculation #{calculation_count})")

# Stage 2: Complete implementation with physics (20 lines)
# Continue from Stage 1
def gravitational_force(pos1, pos2, m1, m2):
    """Calculate gravitational force with automatic caching."""
    global calculation_count
    
    cache_key = (pos1, pos2, m1, m2)
    if cache_key in cache:
        print(f"  Cache hit!")
        return cache[cache_key]
    
    calculation_count += 1
    print(f"  Computing (calculation #{calculation_count})")
    
    # Unpack positions
    x1, y1, z1 = pos1
    x2, y2, z2 = pos2
    
    # Calculate distance
    dx, dy, dz = x2 - x1, y2 - y1, z2 - z1
    r = (dx**2 + dy**2 + dz**2) ** 0.5
    
    # Gravitational force
    G = 6.674e-8  # cm³/g/s²
    F = G * m1 * m2 / r**2
    
    cache[cache_key] = F
    return F

# Stage 3: Demonstrate caching benefits (10 lines)
# Simulate repeated force calculations in N-body code
sun_pos = (0.0, 0.0, 0.0)
earth_pos = (1.496e13, 0.0, 0.0)  # 1 AU in cm
m_sun = 1.989e33  # grams
m_earth = 5.972e27  # grams

print("First calculation:")
F1 = gravitational_force(sun_pos, earth_pos, m_sun, m_earth)

print("\nSecond calculation (same positions):")
F2 = gravitational_force(sun_pos, earth_pos, m_sun, m_earth)

print(f"\nForce: {F1:.2e} dynes")
print(f"Calculations performed: {calculation_count}")
print("In astronomy: Caching speeds up N-body codes 10-100×!")

4.4 The Mutable vs Immutable Distinction¶

Time for one of Python’s most important concepts for scientific computing! Understanding mutability is the key to avoiding mysterious bugs where your simulation state changes unexpectedly. This connects directly to Chapter 1’s defensive programming principles - catching errors early prevents corrupted results.

Python’s Reference Model in Physics Simulations¶

Python doesn’t store values in variables - it stores references to objects. This has profound implications for simulation codes:

# Critical concept for simulation state management!

print("With IMMUTABLE objects (safe for constants):")
G_constant = 6.674e-8
G_backup = G_constant
print(f"Initially: G={G_constant:.3e}, backup={G_backup:.3e}")
print(f"Same object in memory? {G_constant is G_backup}")

G_backup = 6.674e-7  # Creates NEW object (typo won't affect original)
print(f"After change: G={G_constant:.3e}, backup={G_backup:.3e}")
print("Original constant is safe!")

print("\n" + "="*50)
print("With MUTABLE objects (dangerous for state!):")

# Particle system state
particles = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]  # Positions
backup = particles  # Think you're making a backup?

print(f"Initially: same object? {particles is backup}")

# Modify what you think is the backup
backup[0][0] = 999.0  # Trying to test something

print(f"Original particles: {particles}")
print(f"'Backup': {backup}")
print("😱 You just corrupted your simulation state!")

print("\nLesson: Use copy.deepcopy() for simulation state backups!")

The Classic Mutable Default Argument Bug¶

This bug is so dangerous that Python linters specifically check for it. It’s particularly insidious in iterative simulations:

# Stage 1: Demonstrate the bug (15 lines)
# THE BUG: Default mutable created ONCE at function definition!
def accumulate_energies_buggy(energy, history=[]):
    """BUGGY: Trying to track energy history."""
    history.append(energy)
    return history

# Simulate multiple independent runs
print("Run 1:")
run1 = accumulate_energies_buggy(100)
run1 = accumulate_energies_buggy(95)
print(f"  Energy history: {run1}")

print("\nRun 2 (should be independent):")
run2 = accumulate_energies_buggy(200)
print(f"  Energy history: {run2}")
print("  Contains Run 1 data! Runs are coupled! 😱")

# Stage 2: Show the fix (15 lines)
print("THE FIX: Use None sentinel pattern")

def accumulate_energies_fixed(energy, history=None):
    """Safe version - creates new list for each run."""
    if history is None:
        history = []  # Fresh list for each simulation
    history.append(energy)
    return history

print("\nRun 3:")
run3 = accumulate_energies_fixed(100)
run3 = accumulate_energies_fixed(95)
print(f"  Energy history: {run3}")

print("\nRun 4 (properly independent):")
run4 = accumulate_energies_fixed(200)
print(f"  Energy history: {run4}")
print("  Runs are independent! ✅")

Shallow vs Deep Copies in Grid Simulations¶

This distinction is critical for grid-based simulations in computational fluid dynamics, stellar atmospheres, or any field-based calculation:

# Stage 1: Set up the problem (12 lines)
import copy

# Create a 2D grid for temperature distribution
print("Simulating heat diffusion on a grid:")
grid = [[20.0, 20.0, 20.0],
        [20.0, 50.0, 20.0],  # Hot spot in center
        [20.0, 20.0, 20.0]]
print(f"Initial grid: {grid[1]}")  # Middle row

# Shallow copy - DANGEROUS for grids!
print("\n--- SHALLOW COPY (aliases inner arrays) ---")
grid_next_shallow = grid.copy()

# Stage 2: Show the shallow copy problem (10 lines)
# Try to update center point for next timestep
grid_next_shallow[1][1] = 35.0  # Cooling

print(f"Original grid center: {grid[1][1]}°C")
print(f"Next step grid center: {grid_next_shallow[1][1]}°C")
print("😱 Modified both grids! Simulation is corrupted!")

# Reset for deep copy demo
grid = [[20.0, 20.0, 20.0],
        [20.0, 50.0, 20.0],
        [20.0, 20.0, 20.0]]

# Stage 3: Show the deep copy solution (12 lines)
# Deep copy - SAFE for grids!
print("\n--- DEEP COPY (independent arrays) ---")
grid_next_deep = copy.deepcopy(grid)
grid_next_deep[1][1] = 35.0  # Cooling

print(f"Original grid center: {grid[1][1]}°C")
print(f"Next step grid center: {grid_next_deep[1][1]}°C")
print("✅ Grids are independent! Simulation is correct!")

print("\nMemory visualization:")
print("Shallow: grid → [ref1, ref2, ref3] → same inner arrays")
print("Deep:    grid → [ref1, ref2, ref3] → independent arrays")

4.5 Dictionaries: O(1) Lookup Magic for Physics¶

Now we explore one of computer science’s most elegant inventions - the dictionary! In computational physics, dictionaries are perfect for particle properties, lookup tables, and caching expensive calculations. You’re about to understand exactly how they achieve near-instantaneous lookups regardless of size.

Organizing Simulation Data with Dictionaries¶

Let’s see how dictionaries transform particle management in an N-body simulation:

# Stage 1: Basic dictionary structure (15 lines)
# N-body simulation: Managing particle properties efficiently
import json

# Individual particle data (could have thousands)
particles = {
    'particle_0001': {
        'mass': 1.989e30,  # grams (0.001 solar masses)
        'position': [1.5e13, 0.0, 0.0],  # cm
        'velocity': [0.0, 3.0e6, 0.0],   # cm/s
        'type': 'star'
    },
    'particle_0002': {
        'mass': 5.972e27,  # grams (Earth mass)
        'position': [2.5e13, 0.0, 0.0],
        'velocity': [0.0, 2.5e6, 0.0],
        'type': 'planet'
    }
}

# Stage 2: O(1) lookup demonstration (15 lines)
# O(1) lookup by particle ID - instant access!
target = 'particle_0002'
print(f"Accessing {target}:")
print(f"  Mass: {particles[target]['mass']:.2e} g")
print(f"  Type: {particles[target]['type']}")

# Organize by type for efficient group operations
by_type = {}
for pid, data in particles.items():
    ptype = data['type']
    if ptype not in by_type:
        by_type[ptype] = []
    by_type[ptype].append(pid)

print(f"\nParticles by type: {by_type}")
print("In astronomy: Perfect for star/planet/dark matter separation")

Understanding Hash Tables for Physics Applications¶

Let’s demystify how dictionaries achieve O(1) lookup - critical for lookup tables in equation of state calculations:

def demonstrate_hashing_for_physics():
    """Show how hash tables enable fast lookups in physics codes."""
    
    # Common lookup table keys in physics simulations
    physics_keys = [
        (1.0e6, 1.0e-3),   # (Temperature, Density) for EOS (equation of state)
        (2.0e6, 2.0e-3),
        (5.0e6, 1.0e-2),
        "H_ionization",     # Reaction rates
        "He_ionization",
        "opacity_table"     # Opacity types
    ]
    
    print("How hash tables accelerate physics lookups:")
    print("=" * 60)
    print("Key                  → Hash        → Bucket")
    print("-" * 60)
    
    for key in physics_keys:
        hash_value = hash(key)
        # Simplified bucket calculation
        bucket = abs(hash_value) % 100
        if isinstance(key, tuple):
            key_str = f"T={key[0]:.0e}, ρ={key[1]:.0e}"
        else:
            key_str = str(key)
        print(f"{key_str:20s} → {hash_value:11d} → {bucket:3d}")
    
    print("\nThe O(1) lookup process:")
    print("1. Hash the key → integer")
    print("2. Map to bucket index")
    print("3. Direct array access!")
    print("\nThis is why equation of state tables are fast!")

demonstrate_hashing_for_physics()

Dictionary Performance for Particle Lookups¶

Time to see the dramatic performance difference for particle system queries:

# Stage 1: Create test data (10 lines)
import time
import random

# Create a large particle system
n_particles = 1_000_000
print(f"Creating system with {n_particles:,} particles...")

# List approach (what beginners might try)
particle_list = [(f"particle_{i:07d}", random.uniform(1e27, 1e30)) 
                 for i in range(n_particles)]

# Stage 2: Create dictionary version (8 lines)
# Dictionary approach (professional solution)
particle_dict = {f"particle_{i:07d}": random.uniform(1e27, 1e30) 
                 for i in range(n_particles)}

# Search for specific particle
target = "particle_0500000"

print(f"Finding {target} among {n_particles:,} particles...")

# Stage 3: Time the lookups (15 lines)
# Time list search - O(n)
start = time.perf_counter()
for pid, mass in particle_list:
    if pid == target:
        mass_list = mass
        break
list_time = time.perf_counter() - start

# Time dict lookup - O(1)
start = time.perf_counter()
mass_dict = particle_dict[target]
dict_time = time.perf_counter() - start

print(f"\nList search: {list_time*1000:.3f} ms")
print(f"Dict lookup: {dict_time*1000:.6f} ms")
print(f"Dictionary is {list_time/dict_time:,.0f}× faster!")
print("\nFor astronomical catalogs with millions of objects:")
print("This difference determines feasibility!")

4.7 Memory and Performance for Scientific Computing¶

Now that we understand how different data structures behave algorithmically with sets, dictionaries, and lists, let’s examine their actual memory footprint and cache performance implications. Understanding memory layout is crucial for scientific computing because it explains why specialized numerical libraries are so much faster than pure Python and will help you write cache-efficient simulation codes.

Memory Usage in Particle Systems¶

Let’s examine memory costs for different data structures in a particle simulation context:

# Stage 1: Basic Memory Measurement (8 lines)
import sys

def measure_particle_memory(n=10000):
    """Compare memory for particle ID storage."""
    ids_list = [f"p_{i:06d}" for i in range(n)]
    ids_set = set(ids_list)
    
    print(f"Storing {n:,} particle IDs:")
    print(f"  List: {sys.getsizeof(ids_list):,} bytes")
    print(f"  Set:  {sys.getsizeof(ids_set):,} bytes")

measure_particle_memory()

# Stage 2: Understanding the Tradeoff (10 lines)
def understand_physics_tradeoffs():
    """Memory vs speed tradeoffs in physics codes."""
    # Small particle system
    n = 100
    positions = [[i*1.0, 0.0, 0.0] for i in range(n)]
    
    list_bytes = sys.getsizeof(positions)
    # If we used spatial hashing (dict of cells)
    dict_bytes = sys.getsizeof({i: pos for i, pos in enumerate(positions)})
    
    print(f"100 particles as list: {list_bytes:,} bytes")
    print(f"With spatial hash: ~{dict_bytes:,} bytes")
    print(f"Extra {dict_bytes - list_bytes:,} bytes enables O(1) neighbor finding!")

understand_physics_tradeoffs()

# Stage 3: Complete Memory Profile (12 lines)
def profile_simulation_structures(n_particles=1000):
    """Memory comparison for complete particle system."""
    import numpy as np  # Preview of Chapter 7
    
    structures = {
        'List of lists': [[i, i*1.0, 0.0, 0.0] for i in range(n_particles)],
        'List of tuples': [(i, i*1.0, 0.0, 0.0) for i in range(n_particles)],
        'Dictionary': {i: [i*1.0, 0.0, 0.0] for i in range(n_particles)},
        'NumPy (preview)': np.zeros((n_particles, 4))
    }
    
    print(f"Memory for {n_particles:,} particles (ID, x, y, z):")
    print("-" * 50)
    for name, struct in structures.items():
        size = sys.getsizeof(struct)
        per_particle = size / n_particles
        print(f"{name:15s}: {size:8,} bytes ({per_particle:.1f} bytes/particle)")

profile_simulation_structures()
print("\nNotice NumPy's dramatic efficiency - that's Chapter 7!")

Cache Efficiency in Grid Computations¶

Modern CPUs are fast, but memory access is slow. Understanding cache efficiency is critical for grid-based simulations:

# Stage 1: Setup (8 lines)
import time

# Demonstrate cache effects in finite difference calculations
size = 500  # Grid size for heat equation
grid = [[20.0 + 0.1*i*j for j in range(size)] 
        for i in range(size)]

print("Simulating heat diffusion finite differences:")
print("(Common in CFD and atmospheric modeling)")

# Stage 2: Row-major access (10 lines)
# Row-major access (cache-friendly in Python)
start = time.perf_counter()
total = 0.0
for i in range(size):
    for j in range(size):
        total += grid[i][j]
row_time = time.perf_counter() - start

print(f"Processing {size}×{size} grid:")
print(f"Row-major (cache-friendly):   {row_time*1000:.1f} ms")

# Stage 3: Column-major comparison (12 lines)
# Column-major access (cache-hostile in Python)
start = time.perf_counter()
total = 0.0
for j in range(size):
    for i in range(size):
        total += grid[i][j]
col_time = time.perf_counter() - start

print(f"Column-major (cache-hostile): {col_time*1000:.1f} ms")
print(f"Column-major is {col_time/row_time:.1f}× slower!")

print("\nNote: Python's object overhead masks some effects.")
print("In NumPy/C/Fortran with contiguous arrays, this difference")
print("can be 10-100×! Critical for CFD and climate codes.")

4.8 Choosing the Right Structure for Physics¶

After exploring all options, how do you choose? Here’s your decision framework for computational physics:

Performance Quick Reference for Physics¶

* Dict searches by key (particle ID)
† Amortized - occasionally O(n) during resize

Real Example: Combining Data Structures¶

Let’s see how different data structures work together in a simple particle tracking system:

# Simple example: Tracking particles with multiple data structures

# Constants (tuple - can't be changed accidentally!)
CONSTANTS = (6.674e-8, 2.998e10, 1.381e-16)  # G, c, k_B

# Particle data (dictionary for O(1) lookup by ID)
particles = {
    'p001': {'mass': 1e30, 'x': 0.0, 'y': 0.0},
    'p002': {'mass': 2e30, 'x': 10.0, 'y': 0.0},
    'p003': {'mass': 1.5e30, 'x': 5.0, 'y': 5.0},
}

# Active particles (set for fast membership testing)
active = {'p001', 'p002'}  # p003 is inactive

# Time series measurements (list - ordered by time)
measurements = [
    (0.0, 100.5),  # (time, energy)
    (1.0, 99.8),
    (2.0, 99.1),
]

# Example: Find total mass of active particles
total_mass = 0
for pid in active:  # O(1) membership test for each
    if pid in particles:  # O(1) lookup
        total_mass += particles[pid]['mass']

print(f"Total active mass: {total_mass:.1e} g")
print(f"Number of measurements: {len(measurements)}")
print(f"Constants are protected: type(CONSTANTS) = {type(CONSTANTS).__name__}")

# This simple combination shows:
# - Tuple for constants (immutable)
# - Dictionary for particle data (O(1) lookup)
# - Set for active particles (O(1) membership)
# - List for time series (ordered)

Main Takeaways¶

You’ve just mastered concepts that will transform your computational physics from toy problems to research-grade simulations! The journey from understanding simple lists to architecting with dictionaries and sets represents a fundamental shift in how you approach scientific computing. You now see data structures not just as containers, but as carefully chosen tools where the right choice can mean the difference between simulations that finish in hours versus weeks.

The most profound insight from this chapter is that data structure choice often matters more than algorithm optimization. We saw how switching from a list to a set for particle lookups gave us a 100,000× speedup - no amount of code optimization could achieve that! This is the secret that separates research codes from student projects: professionals spend more time architecting their data organization than writing physics equations. The best physics insight in the world is useless if your code can’t handle realistic problem sizes.

The mutable versus immutable distinction that seemed abstract at first is actually critical for scientific computing. Every time you use a tuple for physical constants or configuration parameters, you’re preventing bugs that have literally crashed spacecraft and corrupted published results. When you properly use defensive copying for your simulation state, you’re protecting yourself from the aliasing bugs that have plagued computational physics for decades. Remember: in the Cassini example, a simple data structure choice nearly lost a billion-dollar mission.

The performance principles you’ve learned extend far beyond Python. The cache efficiency concepts explain why codes like LAMMPS and GADGET obsess over memory layout - it’s not premature optimization, it’s the difference between simulating thousands versus millions of atoms. The Big-O notation you’ve mastered is the universal language for discussing whether an algorithm scales to galaxy-sized problems. The hash table concept underlying dictionaries appears in every parallel communication library, every adaptive mesh refinement code, and every spatial indexing scheme you’ll encounter.

Most importantly, you now understand the connection between data structures and numerical methods. That O(n²) all-pairs particle interaction that’s impossible for large systems? It becomes O(n log n) with tree-based structures. The O(n) searching that makes timesteps crawl? It becomes O(1) with spatial hashing. These aren’t just computer science concepts - they’re the difference between simulating a few hundred particles and modeling entire galaxies with billions of stars.

Remember: every data structure makes trade-offs. Lists give you ordering but slow searches. Dictionaries trade memory for lightning-fast lookups. Sets provide uniqueness and mathematical operations but lose ordering. There’s no universally “best” structure - only the best structure for your specific physics problem. As you tackle research simulations, you’ll combine multiple structures: NumPy arrays for number crunching, dictionaries for particle properties, sets for collision detection, and spatial data structures for neighbor finding.

With this foundation, you’re ready to architect simulations that scale to astronomical proportions. The next chapter on functions and modules will show you how to organize this knowledge into reusable, testable components. After that, NumPy will supercharge your numerical operations using the memory layout principles you now understand. You’re no longer just writing code - you’re engineering solutions to the computational challenges of modern physics!

Definitions¶

Aliasing: When two or more variables refer to the same object in memory, causing modifications through one to affect the others.

Amortized O(1): An operation that is usually O(1) but occasionally O(n), where the average over many operations remains O(1).

Big-O Notation: Mathematical notation describing how an algorithm’s runtime scales with input size, focusing on the dominant term.

Cache: Small, fast memory close to the CPU that stores recently accessed data for quick retrieval.

Data Structure: A way of organizing data in computer memory to enable efficient access and modification.

Deep Copy: Creating a completely independent copy of an object and all objects it contains, recursively.

Dictionary: A mutable mapping type storing key-value pairs with O(1) average lookup time using hash tables.

Hash Function: A function mapping data of arbitrary size to fixed-size values, enabling fast lookups.

Hash Table: The underlying implementation for dictionaries and sets, enabling O(1) average-case lookups.

Immutable: Objects whose state cannot be modified after creation (tuples, strings, numbers).

List: Python’s built-in mutable sequence type that stores references to objects in order.

Memoization: Caching technique storing function results to avoid recalculating for the same inputs.

Mutable: Objects whose state can be modified after creation (lists, dictionaries, sets).

Named Tuple: A tuple subclass allowing element access by name as well as index.

O(1) - Constant Time: An operation whose runtime doesn’t depend on input size.

O(n) - Linear Time: An operation whose runtime grows proportionally with input size.

O(n²) - Quadratic Time: An operation whose runtime grows with the square of input size.

O(log n) - Logarithmic Time: An operation whose runtime grows logarithmically with input size.

O(n log n): Common complexity for efficient sorting algorithms and tree-based operations.

Reference: A variable that points to an object in memory rather than containing the value directly.

Set: A mutable collection of unique, unordered elements with O(1) membership testing.

Shallow Copy: Creating a new container with references to the same contained objects.

Spatial Hashing: Organizing particles by spatial region for O(1) neighbor finding.

Tuple: An immutable sequence type that cannot be changed after creation.

Key Takeaways¶

✓ Data structure choice can change performance by factors of 1,000,000× or more for large-scale simulations

✓ Lists are versatile but have O(n) search - use for ordered particle arrays that you’ll vectorize with NumPy

✓ Dictionaries and sets provide O(1) lookup through hash tables - essential for particle lookups and caching

✓ Tuples prevent modification bugs - use for physical constants and configuration parameters

✓ The shallow vs deep copy distinction is critical for grid-based simulations

✓ Python stores references to objects, explaining why NumPy’s contiguous arrays are faster

✓ Memory layout affects cache performance by 2-10× even in Python

✓ Every data structure trades something (memory, speed, flexibility) for something else

✓ Mutable default arguments are dangerous - always use the None sentinel pattern (Chapter 1 callback)

✓ Spatial data structures transform O(n²) physics problems into O(n) or O(n log n)

✓ Sets provide elegant mathematical operations for domain decomposition and particle tracking

✓ Dictionaries enable memoization that can speed up equation of state calculations 100×

Python Module & Method Reference¶

This reference section catalogs all Python modules, functions, and methods introduced in this chapter. Keep this as a quick lookup guide for your physics simulations.

Standard Library Modules¶

time module - High-resolution timing for performance measurement

import time

• time.perf_counter() - Returns float seconds with highest available resolution

  • Use for timing code segments: start = time.perf_counter()

  • Always use for benchmarking (not time.time() which can go backwards!)

sys module - System-specific parameters (expanded from Chapter 1)

import sys

• sys.getsizeof(object) - Returns memory size of object in bytes

  • Includes object overhead but not contained objects

  • Essential for memory profiling in simulations

copy module - Create object copies with control over depth

import copy

• copy.copy(x) - Creates shallow copy (new container, same contents)

• copy.deepcopy(x) - Creates deep copy (recursively copies all contents)

  • Critical for grid simulations to avoid aliasing bugs

  • Use for simulation state backups

random module - Generate random numbers for Monte Carlo

import random

• random.random() - Random float in [0.0, 1.0)

• random.uniform(a, b) - Random float in [a, b]

• random.choice(seq) - Random element from sequence

• random.sample(population, k) - k unique random elements

json module - Save/load structured data

import json

• json.dump(obj, file, indent=2) - Write object to file with formatting

• json.load(file) - Read object from file

• json.dumps(obj) - Convert object to JSON string

• json.loads(string) - Parse JSON string to object

math module - Mathematical functions (review from Chapter 2)

import math

• math.sqrt(x) - Square root (use for distances)

• math.log(x) - Natural logarithm

• math.log2(x) - Base-2 logarithm (for Big-O analysis)

Collections Module¶

collections module - Specialized container datatypes

from collections import OrderedDict, Counter, defaultdict, deque, namedtuple

OrderedDict - Dictionary that remembers insertion order

• .move_to_end(key, last=True) - Move key to end (or beginning if last=False)

• .popitem(last=True) - Remove and return (key, value) pair from end/beginning

• Use for LRU caches in physics calculations

Counter - Dictionary subclass for counting hashable objects

• Counter(iterable) - Create from iterable

• .most_common(n) - Returns n most common (element, count) pairs

• .update(iterable) - Add counts from iterable

• Perfect for histogram analysis of particle properties

defaultdict - Dictionary with default value factory

• defaultdict(list) - Missing keys default to empty list

• defaultdict(int) - Missing keys default to 0

• Eliminates KeyError in particle grouping

deque - Double-ended queue with O(1) operations at both ends

• deque(maxlen=n) - Fixed-size queue (old elements dropped)

• .append(x) / .appendleft(x) - Add to right/left end

• .pop() / .popleft() - Remove from right/left end

• Essential for boundary conditions in simulations

namedtuple - Tuple subclass with named fields

Particle = namedtuple('Particle', ['id', 'mass', 'x', 'y', 'z'])
p = Particle(1, 1.0e30, 0, 0, 0)
print(p.mass)  # Access by name

Functools Module¶

functools module - Higher-order functions and operations

from functools import lru_cache

@lru_cache(maxsize=128) - Decorator for automatic memoization

@lru_cache(maxsize=1000)
def expensive_calculation(T, rho):
    return complex_physics_calculation(T, rho)

• Caches function results automatically

• maxsize=None for unlimited cache

• .cache_info() - Returns cache statistics

• .cache_clear() - Clear the cache

Built-in Functions for Data Structures¶

Type Checking

• type(obj) - Returns object’s type

• isinstance(obj, type) - Check if obj is instance of type

• id(obj) - Returns unique identifier (memory address)

• obj1 is obj2 - Check if same object (identity)

Container Operations

• len(container) - Number of elements

• x in container - Membership test (O(1) for sets/dicts, O(n) for lists)

• sorted(iterable) - Returns new sorted list

• reversed(sequence) - Returns reverse iterator

• enumerate(iterable, start=0) - Returns (index, value) pairs

• zip(iter1, iter2, ...) - Parallel iteration

Copying

• list(iterable) - Create new list from iterable

• tuple(iterable) - Create new tuple from iterable

• set(iterable) - Create new set from iterable

• dict(pairs) - Create dictionary from (key, value) pairs

Data Structure Methods Summary¶

List Methods (Expanded)

• .append(x) - Add to end (amortized O(1))

• .extend(iterable) - Add all elements from iterable

• .insert(i, x) - Insert at position i (O(n))

• .remove(x) - Remove first x (O(n))

• .pop(i=-1) - Remove and return element at i

• .clear() - Remove all elements

• .index(x) - Find position of x (O(n))

• .count(x) - Count occurrences of x

• .sort() - Sort in-place

• .reverse() - Reverse in-place

• .copy() - Create shallow copy

Dictionary Methods (Expanded)

• .get(key, default=None) - Safe access with default

• .setdefault(key, default) - Get or set with default

• .pop(key, default) - Remove and return value

• .popitem() - Remove and return arbitrary (key, value)

• .update(other) - Update with other dict/pairs

• .keys() - View of keys

• .values() - View of values

• .items() - View of (key, value) pairs

• .clear() - Remove all items

• .copy() - Shallow copy

Set Methods (Expanded)

• .add(x) - Add element

• .remove(x) - Remove element (raises KeyError if missing)

• .discard(x) - Remove element (no error if missing)

• .pop() - Remove and return arbitrary element

• .clear() - Remove all elements

• .union(other) or | - Elements in either set

• .intersection(other) or & - Elements in both sets

• .difference(other) or - - Elements in first but not second

• .symmetric_difference(other) or ^ - Elements in either but not both

• .update(other) or |= - Add elements from other

• .intersection_update(other) or &= - Keep only elements in both

• .difference_update(other) or -= - Remove elements in other

• .symmetric_difference_update(other) or ^= - Keep elements in either but not both

• .issubset(other) - Check if all elements in other

• .issuperset(other) - Check if contains all elements of other

• .isdisjoint(other) - Check if no elements in common

• .copy() - Create shallow copy

Quick Reference Tables¶

Data Structure Operations¶

Performance Quick Reference¶

* Dict searches by key, not value
† Amortized - occasionally O(n) during resize

Common Methods for Physics¶

When to Use Each Structure¶

Next Chapter Preview¶

With data structures mastered, Chapter 5 will explore functions and modules - how to organize code for reusability, testing, and collaboration. You’ll learn how Python’s function model, with first-class functions and closure support, enables powerful patterns like decorators and functional programming techniques. These concepts prepare you for building modular physics engines, creating reusable analysis pipelines, and understanding the functional programming paradigm used in modern frameworks like JAX where functions transform into automatically differentiable computational graphs. Get ready to transform your scripts into professional-quality scientific libraries!