5. Custom analysis with the Python API

Notebook 2 covered the aggregates: how long each region took in total, on average, per rank. This notebook goes one level down, to the individual calls, which is where you go when the question you have is not one the built-in summaries answer.

Three tools do most of the work:

Method	Gives you
results.events()	one entry per recorded call
results.to_events_dataframe()	the same, as a pandas DataFrame
results.call_stack()	the calls with their nesting reconstructed

As everywhere in this API, times are in seconds.

import tempfile
import time
from pathlib import Path

from scope_profiler import ProfileManager, read_h5

WORKDIR = Path(tempfile.mkdtemp(prefix="scope-profiler-tutorial-"))
DATA_FILE = WORKDIR / "profiling_data.h5"

ProfileManager.setup(file_path=str(DATA_FILE))


@ProfileManager.profile("assemble")
def assemble():
    time.sleep(0.002)


with ProfileManager.profile_region("setup"):
    time.sleep(0.01)

for step in range(6):
    with ProfileManager.profile_region("timestep"):
        assemble()
        with ProfileManager.profile_region("solve"):
            # The solve gets slower as the run progresses.
            time.sleep(0.002 + 0.001 * step)
    if step % 3 == 2:
        with ProfileManager.profile_region("checkpoint"):
            time.sleep(0.004)

ProfileManager.finalize(verbose=False)

Reading the results back in the same script

ProfileManager.read_results() opens the file the current configuration just wrote, so you do not have to repeat the path. Under MPI only rank 0 writes the merged file, so guard the call with a rank check there.

results = ProfileManager.read_results()
print(results)
print("same file:", results.file_path == DATA_FILE)

<ProfilingResults 'profiling_data.h5': 5 region(s), 1 rank(s)>
same file: True

One entry per call

events() returns the long-form (“tidy”) view of the run: every recorded call, in region order, then rank, then call order.

events = results.events()

print(f"{len(events)} calls recorded\n")
for event in events[:4]:
    print(event)

21 calls recorded

{'name': 'setup', 'rank': 0, 'call_index': 0, 'start': 0.001959806000002118, 'end': 0.012039456000024984, 'duration': 0.010079650000022866, 'call_id': 0, 'parent_id': -1}
{'name': 'timestep', 'rank': 0, 'call_index': 0, 'start': 0.012221035999999685, 'end': 0.016385643000006667, 'duration': 0.004164607000006981, 'call_id': 1, 'parent_id': -1}
{'name': 'timestep', 'rank': 0, 'call_index': 1, 'start': 0.016389048000007733, 'end': 0.021536933000021463, 'duration': 0.00514788500001373, 'call_id': 4, 'parent_id': -1}
{'name': 'timestep', 'rank': 0, 'call_index': 2, 'start': 0.021540118000018538, 'end': 0.02768315900001994, 'duration': 0.0061430410000014035, 'call_id': 7, 'parent_id': -1}

Each entry carries the region name, the rank that recorded it, its call_index within that rank, and start / end / duration in seconds.

The timeline starts at zero

By default the timestamps are measured from the start of the run — ProfileManager.setup() records that instant, and ProfilingResults exposes it as run_start_time. The raw values come from perf_counter_ns(), a monotonic clock with an arbitrary origin that is comparable only within one run, and you can still get them with relative=False.

print("relative :", [round(e["start"], 4) for e in results.events(include="timestep")])
print(
    "absolute :",
    [round(e["start"], 4) for e in results.events(include="timestep", relative=False)],
)

print("\nrun start    :", results.run_start_time)
print("first entry  :", results.minimum_start_time)
print("last exit    :", results.maximum_end_time)
print("wall clock   :", results.time_span, "s")

relative : [0.0122, 0.0164, 0.0215, 0.0318, 0.0389, 0.0471]
absolute : [163.4414, 163.4456, 163.4508, 163.461, 163.4681, 163.4763]

run start    : 163.429213445
first entry  : 163.431173251
last exit    : 163.489485019
wall clock   : 0.058311768000010034 s

The time before the first region

setup() runs before any region is entered, so the gap between the two is work the instrumentation never saw — imports, reading input, building a mesh. It is often the first surprise a profile delivers. The run’s origin is the moment setup() is called, so calling it as early as possible is what makes that gap visible.

print(
    f"{results.startup_time * 1e3:.2f} ms elapsed before the first region was entered"
)

1.96 ms elapsed before the first region was entered

Choosing a different origin

time_origin is the zero point in use: the registered start time, or the first region entry for files that have none. Override it per call with origin.

The plot_* functions are the one exception — they frame their x axis on the first region entry, so a long startup does not fill a chart with empty space. The last line below reproduces exactly what a chart’s axis shows.

None of this is required reading for older files: one written before setup() recorded a start time reports run_start_time as None, falls back to the first region entry, and behaves exactly as it always did.

print("origin       :", results.time_origin)
print("is run start :", results.time_origin == results.run_start_time)

print(
    "\nfrom run start   :",
    [round(e["start"], 4) for e in results.events(include="solve")],
)
print(
    "from first entry :",
    [
        round(e["start"], 4)
        for e in results.events(include="solve", origin=results.minimum_start_time)
    ],
)

origin       : 163.429213445
is run start : True

from run start   : [0.0143, 0.0185, 0.0236, 0.0339, 0.041, 0.0491]
from first entry : [0.0124, 0.0165, 0.0217, 0.0319, 0.039, 0.0472]

time_span is the width of the profiled window, which is what you want when asking what fraction of the run went into this region — the sum of all region durations double-counts nested regions and can exceed the wall clock.

for region in results.get_regions():
    share = region.total_duration / results.time_span
    print(f"  {region.name:<11} {share:6.1%} of the run")

  setup        17.3% of the run
  timestep     68.4% of the run
  assemble     21.3% of the run
  solve        47.0% of the run
  checkpoint   14.0% of the run

Selecting what you get

events() takes the same include / exclude regexes as the rest of the API, plus ranks:

results.events(include="solve")             # one region
results.events(include=["solve", "assemble"])
results.events(exclude="checkpoint")
results.events(ranks=0)                     # one rank, or ranks=[0, 2]

solve_events = results.events(include="solve")
print("calls  :", len(solve_events))
print("indices:", [e["call_index"] for e in solve_events])
print("slowest:", max(solve_events, key=lambda e: e["duration"]))

calls  : 6
indices: [0, 1, 2, 3, 4, 5]
slowest: {'name': 'solve', 'rank': 0, 'call_index': 5, 'start': 0.04913167200001567, 'end': 0.05619917399999963, 'duration': 0.0070675019999839606, 'call_id': 19, 'parent_id': 17}

As a DataFrame

to_events_dataframe() returns the same data with one row per call, which is usually the shortest path from a question to an answer. It needs pandas (part of the plot extra).

frame = results.to_events_dataframe()
frame.head()

	name	rank	call_index	start	end	duration
0	setup	0	0	0.001960	0.012039	0.010080
1	timestep	0	0	0.012221	0.016386	0.004165
2	timestep	0	1	0.016389	0.021537	0.005148
3	timestep	0	2	0.021540	0.027683	0.006143
4	timestep	0	3	0.031777	0.038923	0.007147

From here every question is a groupby. A few that come up constantly:

# How consistent is each region, call to call?
frame.groupby("name")["duration"].agg(["count", "mean", "std", "max"]).sort_values(
    "mean", ascending=False
)

	count	mean	std	max
name
setup	1	0.010080	NaN	0.010080
timestep	6	0.006646	0.001861	0.009136
solve	6	0.004568	0.001870	0.007068
checkpoint	2	0.004069	0.000003	0.004071
assemble	6	0.002067	0.000002	0.002071

# The ten slowest individual calls in the whole run.
frame.nlargest(10, "duration")[["name", "rank", "call_index", "start", "duration"]]

	name	rank	call_index	start	duration
0	setup	0	0	0.001960	0.010080
6	timestep	0	5	0.047064	0.009136
5	timestep	0	4	0.038926	0.008135
4	timestep	0	3	0.031777	0.007147
18	solve	0	5	0.049132	0.007068
3	timestep	0	2	0.021540	0.006143
17	solve	0	4	0.040995	0.006065
2	timestep	0	1	0.016389	0.005148
16	solve	0	3	0.033850	0.005072
1	timestep	0	0	0.012221	0.004165

# Time spent per region, per rank - the shape load-imbalance analysis wants.
frame.pivot_table(index="rank", columns="name", values="duration", aggfunc="sum")

name	assemble	checkpoint	setup	solve	timestep
rank
0	0.012401	0.008137	0.01008	0.027411	0.039874

Down to one region

Region and MPIRegion expose the same view for a single region, so you can narrow before you widen. Region also hands back the stored integers directly, for the rare case where the nanosecond values matter more than the convenience of seconds.

solve = results["solve"]  # an MPIRegion: the region across all ranks
print("all ranks   :", len(solve.events()), "calls")
print("rank 0 only :", len(solve.events(ranks=0)), "calls")

rank0 = solve[0]  # a Region: one rank
print("\nseconds     :", rank0.durations[:3])
print("nanoseconds :", rank0.durations_ns[:3])
print("start (ns)  :", rank0.start_times_ns[:3])

all ranks   : 6 calls
rank 0 only : 6 calls

seconds     : [0.00206877 0.00306926 0.00406826]
nanoseconds : [2068768 3069262 4068264]
start (ns)  : [163443529098 163447675037 163452827088]

The call stack

Regions record no call graph of their own — each call is just a (start, end) pair. Nesting is therefore reconstructed from containment: a call that starts while another is still open is treated as its child. That is what the flame graph draws, and call_stack() hands you the same structure as plain dicts.

calls = results.call_stack(rank=0)

for call in calls[:8]:
    indent = "  " * call["depth"]
    print(f"{indent}{call['name']:<12} {call['duration'] * 1e3:6.2f} ms")

setup         10.08 ms
timestep       4.16 ms
  assemble       2.07 ms
  solve          2.07 ms
timestep       5.15 ms
  assemble       2.07 ms
  solve          3.07 ms
timestep       6.14 ms

Every call carries its depth and the index of its parent in the returned list. Indices rather than names, because a region called repeatedly — or recursively — contributes several entries under one name.

call_stack_roots() and call_stack_children() turn that flat list into a tree you can walk:

from scope_profiler import call_stack_children, call_stack_roots

children = call_stack_children(calls)

for root in call_stack_roots(calls):
    names = [calls[child]["name"] for child in children[root]]
    print(f"{calls[root]['name']:<12} contains {names}")

setup        contains []
timestep     contains ['assemble', 'solve']
timestep     contains ['assemble', 'solve']
timestep     contains ['assemble', 'solve']
checkpoint   contains []
timestep     contains ['assemble', 'solve']
timestep     contains ['assemble', 'solve']
timestep     contains ['assemble', 'solve']
checkpoint   contains []

Self time

A useful thing the aggregates cannot give you: exclusive (self) time — how long a region took on its own, with the time spent inside its children removed. It falls straight out of the tree.

from collections import defaultdict

total = defaultdict(float)
self_time = defaultdict(float)

for index, call in enumerate(calls):
    nested = sum(calls[child]["duration"] for child in children[index])
    total[call["name"]] += call["duration"]
    self_time[call["name"]] += call["duration"] - nested

print(f"{'region':<12}{'total [ms]':>12}{'self [ms]':>12}")
for name in sorted(total, key=total.get, reverse=True):
    print(f"{name:<12}{total[name] * 1e3:12.2f}{self_time[name] * 1e3:12.2f}")

region        total [ms]   self [ms]
timestep           39.87        0.06
solve              27.41       27.41
assemble           12.40       12.40
setup              10.08       10.08
checkpoint          8.14        8.14

timestep is nearly all children — it is a container. The regions with large self time are the ones actually doing work, and the ones worth optimising.

Several ranks

Everything above is rank-aware; a serial run just happens to have one rank. To show the shape without launching mpirun, here is a two-rank file written by hand — it also documents the on-disk layout, which is plain HDF5:

import h5py
import numpy as np

MPI_FILE = WORKDIR / "two_ranks.h5"
NS = 1_000_000_000  # timestamps are stored as int64 nanoseconds

with h5py.File(MPI_FILE, "w") as h5file:
    for rank, factor in enumerate([1.0, 1.6]):  # rank 1 is the slow one
        regions = h5file.create_group(f"rank{rank}").create_group("regions")
        starts = np.array([0.0, 0.5, 1.0]) * NS
        durations = np.array([0.2, 0.25, 0.2]) * factor * NS
        group = regions.create_group("solve")
        group.create_dataset("start_times", data=starts.astype(np.int64))
        group.create_dataset("end_times", data=(starts + durations).astype(np.int64))

mpi_reader = read_h5(MPI_FILE)
print(mpi_reader)

mpi_reader.to_events_dataframe()

<ProfilingResults 'two_ranks.h5': 1 region(s), 2 rank(s)>

	name	rank	call_index	start	end	duration
0	solve	0	0	0.0	0.20	0.20
1	solve	0	1	0.5	0.75	0.25
2	solve	0	2	1.0	1.20	0.20
3	solve	1	0	0.0	0.32	0.32
4	solve	1	1	0.5	0.90	0.40
5	solve	1	2	1.0	1.32	0.32

# Load imbalance: how much longer did the slowest rank take than the fastest?
per_rank = mpi_reader["solve"].total_durations()
print("total per rank:", per_rank)
print(f"imbalance     : {max(per_rank.values()) / min(per_rank.values()):.2f}x")

total per rank: {0: 0.65, 1: 1.04}
imbalance     : 1.60x

Next

6. Building your own plots takes the same three tools and draws charts with them.