5. Custom analysis with the Python API

Download notebookView on GitHubRun when the site was built.

Install for this tutorial: pip install "scope-profiler[pproc]"; pandas is used for event analysis.

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 time

from scope_profiler import ProfileManager

ProfileManager.setup(file_path="profiling_data.h5")


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


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

for step in range(6):
    with ProfileManager.region("timestep"):
        assemble()
        with ProfileManager.region("solve"):
            # The solve gets slower as the run progresses.
            time.sleep(0.002 + 0.001 * step)
    if step % 3 == 2:
        with ProfileManager.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.name == "profiling_data.h5")

<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.002269132000009222, 'end': 0.013287058000003071, 'duration': 0.01101792599999385, 'call_id': 0, 'parent_id': -1}
{'name': 'timestep', 'rank': 0, 'call_index': 0, 'start': 0.013516113000008545, 'end': 0.017687955000013744, 'duration': 0.004171842000005199, 'call_id': 1, 'parent_id': -1}
{'name': 'timestep', 'rank': 0, 'call_index': 1, 'start': 0.017690228000020625, 'end': 0.02281599799999867, 'duration': 0.005125769999978047, 'call_id': 4, 'parent_id': -1}
{'name': 'timestep', 'rank': 0, 'call_index': 2, 'start': 0.022816990000023907, 'end': 0.028939613000005693, 'duration': 0.006122622999981786, '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. At finalize() the profiler also assigns every call a call_id and the parent_id of the call it ran inside (-1 for a top level call), so a saved profile carries its call graph explicitly.

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.

rel = results.events(include="^timestep$")
raw = results.events(include="^timestep$", relative=False)

print("relative start of the first call [s]:", round(rel[0]["start"], 4))
print(
    "raw start minus run start         [s]:",
    round(raw[0]["start"] - results.run_start_time, 4),
)
print("raw start is a perf_counter value, so it is only comparable within one run")

print("\nwall clock   :", round(results.time_span, 3), "s")

relative start of the first call [s]: 0.0135
raw start minus run start         [s]: 0.0135
raw start is a perf_counter value, so it is only comparable within one run

wall clock   : 0.059 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"
)

2.27 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 is the 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 is the run start : True

from run start   : [0.0156, 0.0198, 0.0249, 0.0351, 0.0422, 0.0503]
from first entry : [0.0134, 0.0175, 0.0226, 0.0328, 0.0399, 0.0481]

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        18.6% of the run
  timestep     67.2% of the run
  assemble     20.9% of the run
  solve        46.2% of the run
  checkpoint   13.7% 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.05033543000001828, 'end': 0.05739506700001584, 'duration': 0.007059636999997565, '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 pproc extra).

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

	name	rank	call_index	start	end	duration
0	setup	0	0	0.002269	0.013287	0.011018
1	timestep	0	0	0.013516	0.017688	0.004172
2	timestep	0	1	0.017690	0.022816	0.005126
3	timestep	0	2	0.022817	0.028940	0.006123
4	timestep	0	3	0.033029	0.040150	0.007121

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.011018	NaN	0.011018
timestep	6	0.006631	0.001857	0.009121
solve	6	0.004561	0.001870	0.007060
checkpoint	2	0.004065	0.000006	0.004069
assemble	6	0.002061	0.000002	0.002064

# 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.002269	0.011018
6	timestep	0	5	0.048274	0.009121
5	timestep	0	4	0.040151	0.008123
4	timestep	0	3	0.033029	0.007121
18	solve	0	5	0.050335	0.007060
3	timestep	0	2	0.022817	0.006123
17	solve	0	4	0.042212	0.006061
2	timestep	0	1	0.017690	0.005126
16	solve	0	3	0.035089	0.005060
1	timestep	0	0	0.013516	0.004172

# 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.012363	0.00813	0.011018	0.027365	0.039785

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] - rank0.start_times_ns[0])

all ranks   : 6 calls
rank 0 only : 6 calls

seconds     : [0.0020629  0.00306078 0.00406066]
nanoseconds : [2062901 3060776 4060664]
start (ns)  : [      0 4130499 9254316]

The call stack

Every call is stored as a (start, end) pair, so call_stack() reconstructs the nesting 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 — each with its depth, its inclusive and exclusive duration, and the call_path of region names leading to it.

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         11.02 ms
timestep       4.17 ms
  assemble       2.06 ms
  solve          2.06 ms
timestep       5.13 ms
  assemble       2.06 ms
  solve          3.06 ms
timestep       6.12 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 per call: exclusive (self) time — how long a region took on its own, with the time spent inside its children removed. Every call carries it as exclusive_duration, next to inclusive_duration; it is what the tree above gives you by subtracting the children.

from collections import defaultdict

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

for call in calls:
    total[call["name"]] += call["inclusive_duration"]
    self_time[call["name"]] += call["exclusive_duration"]

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.78        0.06
solve              27.37       27.37
assemble           12.36       12.36
setup              11.02       11.02
checkpoint          8.13        8.13

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 numpy as np

from scope_profiler import MPIRegion, ProfilingResults, Region

NS = 1_000_000_000  # timestamps are stored as int64 nanoseconds
starts = (np.array([0.0, 0.5, 1.0]) * NS).astype(np.int64)
ranks = {}
for rank, factor in enumerate([1.0, 1.6]):  # rank 1 is the slow one
    durations = (np.array([0.2, 0.25, 0.2]) * factor * NS).astype(np.int64)
    ranks[rank] = Region(starts, starts + durations)

mpi_reader = ProfilingResults({"solve": MPIRegion("solve", ranks)}, num_ranks=2)
print(mpi_reader)

mpi_reader.to_events_dataframe()

<ProfilingResults '': 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.