Skip to content

[pyserial] Implement part of the reconstruction in Python - #437

Open
fwyzard wants to merge 21 commits into
masterfrom
pyserial
Open

[pyserial] Implement part of the reconstruction in Python#437
fwyzard wants to merge 21 commits into
masterfrom
pyserial

Conversation

@fwyzard

@fwyzard fwyzard commented Sep 12, 2026

Copy link
Copy Markdown
Contributor

Build on top of #436 (and currently includes it) and adds the pyserial application.

Most of the modules come in two versions: the original C++ implementation, and a Python implementation.

The DataFormats are basically unchanged, except for some minor quality of life improvements.

The modules ported to Python are the beam spot, the raw to cluster, the rechit producer, and the vertex finder.
For the moment the track finder (doublets, ntuplets, fishbone, broken line fit) is only available in C++.

Under these conditions the python version is about a factor 2× slower than the C++ version - but the proof of principle works: there are no memory copies, the C++ and Python modules can inter-operate without issues, and the results are identical:
throughput-04-weak-scaling
ratio-04-weak-scaling

fwyzard and others added 21 commits September 12, 2026 01:34
.gitignore covered *~ but not the swap files an editor leaves beside a file
while it is open, which is how one of them ended up committed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
A copy of src/fwtest under src/pytest, with the two lines that name the
backend changed.  Everything after this commit changes the fork rather than
fwtest, which stays available as the control: the product registry change in
particular is not backward compatible, and every backend in this repository
already carries its own Framework/, bin/ and DataFormats/.

Two headers come over from serial's Framework instead, because fwtest is the
minimal backend and does not carry them: ReusableObjectHolder.h and
RunningAverage.h, which two of serial's reconstruction modules use.  Nothing
in pytest uses them either.  They are here because a later backend will, and
the two Framework/ trees have to stay byte-identical -- a property worth
having from the start rather than restoring later by editing this commit.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
Both were inherited from fwtest and both were silent.

Wrapper<T> brace-initialises its payload, so event.emplace(token, n, 0.0f)
on a std::vector<float> built a *two-element* vector -- {n, 0.0f} read as an
initializer_list -- rather than n elements, and moving one vector into
another did not compile at all.  A tag-selected parenthesised constructor
is added for the cases that need it; the brace-initialising one stays,
because an aggregate product needs it.

Event was copyable in declaration only: it owns its products through
unique_ptr, so the implicit copy constructor existed but could never
compile.  Nothing copies an Event, so nothing noticed.  Deleting it says
what was always true.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
fwtest keys products on std::type_index alone, so produces<T>() throws if
that type is already registered and one event can hold at most one product
of a type.  Comparing two implementations of an algorithm needs two products
of the same type -- one from each -- so the type alone cannot identify a
product.

The alternative is a distinct C++ type per producer, which is what the rest
of pixeltrack-standalone does.  That makes the type system carry which
implementation produced a value, so a consumer cannot be written once and
pointed at either, and a third implementation means a third type.  Keying on
(type, label) instead is what the demo and CMSSW both do, and it is what
lets a Python module stand in for a C++ one later.

consumes<T>() is now ambiguous when several modules produce T.  Rather than
pick one it throws and names the count; a module that wants a particular one
calls consumes<T>(label).  The unlabelled form stays because most products
are unique, and keeping it means the ported plugins compile unchanged.

The label is the module's own name.  Until the next commits it is the plugin
type, which already separates two different plugins publishing one type.

An empty label means "the only product of that type", so consumes<T>("") and
consumes<T>() are the same lookup.  That is what lets a module take the label
of its input from its configuration and pass it on unconditionally: naming an
input becomes something a job does only when it has more than one candidate.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
The schedule and the modules' parameters come from an INI file, parsed with
boost property_tree: [options] modules = a, b, c gives the path, and one
[label] section per module names its plugin with @type and carries the rest
as that module's parameters.  test.ini reproduces the workflow fwtest runs
from its hard-coded path.

A label is now a module's own name rather than its plugin type, so one
plugin can be scheduled twice under different labels and its two products
tell apart by the (type, label) key.

Every module is constructed as Module(ModuleConfig const&, ProductRegistry&).
One signature rather than two -- with the configuration optional -- costs an
ignored argument in the modules that take no parameters, and buys a compiler
error, rather than a silent fallback, when a module that means to read its
section declares the wrong constructor.

The command line still wins: run-scan.py and the test_* targets drive
--numberOfThreads, --numberOfStreams and --maxEvents, so [options] supplies
defaults and an explicit switch overrides them.

Two mistakes that used to be silent are now errors: an unknown key in
[options], and a scheduled label with no section of its own.

LIBNAMES selected every directory entry that was not a plugin, which would
have turned test.ini into an empty library and linked against it.  A library
is now a directory that holds C++ sources.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
A module may now be a Python class.  PythonProducer takes `script` and
`factory` from its section and hands the whole section to that factory as a
dict, along with the very ProductRegistry it is being constructed against:

    def create(config: dict[str, str], registry) -> object   # with .produce(event)

So a Python module declares exactly as a C++ one does, against the same
object, and the tokens it gets back are the only way it can reach a product.
One instance is constructed per stream, each holding its own Python object,
which is what lets the streams run their Python modules concurrently.

The Event offers three ways to publish one, named for what they do:

    event.allocate(token)          construct it in the Event, nothing set,
                                   and hand back a reference to fill
    event.put(token, value)        copy in a product that already exists
    event.emplace(token, *values)  build it from Python values and move it in

allocate() is the zero-copy path: the reference is the product's own storage.
put() and emplace() return nothing, because what is in the Event is not the
object the caller is holding.  What emplace() can build is decided per
product from the dictionary -- a scalar, a std::vector of scalars, or an
aggregate whose members are all scalars -- and anything else says so and
names allocate().  A product that can be neither copied nor built says that
too, at run time: the SoA wrappers hold a unique_ptr and are move-only.

Free-threading is not optional: with the GIL enabled the Python modules of
different streams serialise, and a throughput compared against an all-C++ run
would measure that rather than anything about Python.  The build refuses to
generate python-config.mk from a GIL-enabled interpreter, and the job prints
the GIL state next to the stream count -- the failure is silent by nature, so
it is made loud in two places.  A free-threaded system python3.14t is used
when there is one, otherwise uv installs CPython under external/; either way
what gets embedded is a venv under external/venv.

No binding is written by hand.  rootcling reflects the product classes named
in REFLECT_PRODUCTS, and tools/generate_bindings.C walks that reflection
through TClass and writes the nanobind registrations for every product and
everything reachable from one -- public data members, fixed-size arrays as
zero-copy numpy views, public methods -- together with the Event's dispatch
and the name -> type_index table.  The product list is written once, in the
Makefile.  Nothing under Framework/ names a product type, which is what will
let that directory be identical in every backend.

Three shapes the reflection does not describe on its own, and what the
generator does about them:

- A column that is a pointer and a length.  A format keeping its storage as a
  unique_ptr<T[]> with the length in another member hands it out as a bare
  pointer, and a bare pointer cannot say how far it runs.  A public method
  returning a std::span -- or an edm::Span, since CMSSW has one -- is a pointer
  and a length together, and becomes a numpy view over the column.  Getting
  std::span past rootcling takes one environment variable: its own option
  parser reads -std as -s, and EXTRA_CLING_ARGS is how a flag reaches cling.
- A class reachable only through a method.  Following public data members alone
  stops at the first wrapper that keeps its storage private, and every method
  returning an inner class is then skipped as unbindable, so the generator
  follows what public methods hand back as well.
- Two classes wanting the same Python name.  The second registration would
  quietly shadow the first, so a colliding name keeps its scope.

ROOT also normalises the standard library's namespace away -- a member of type
std::vector<unsigned int> is reported as "vector<unsigned int>", which is not
C++ that compiles.  The qualification is put back where the name arrives from
ROOT rather than at each of the places that emit one, so that the two spellings
cannot be collected as two classes and bound twice either.

Methods are bound through lambdas returning decltype(auto).  Plain auto deduces
by value, which binds every method returning a reference to a copy: it reads
correctly, and only a write through the reference shows the difference, which
is the kind of bug that surfaces a long way from its cause.

A product that cannot be default-constructed but has a public constructor
wiring it to another product and to a column is built by allocate(token, ...),
from the constructor the reflection reports.

The EventSetup is reached the same way, by type and with no token, through a
dispatch generated from the same product list.  The reference it hands out is
not const -- nanobind has no per-instance constness -- and that is worse here
than for the Event, since an EventSetup product is shared by every stream; it
is noted where it is bound.  edm::EventSetup turned out to be copyable in
declaration only, exactly as the Event was and for the same reason: it owns its
products through unique_ptr, so the implicit copy constructor existed but could
never compile, and nothing noticed until nanobind tried to instantiate it.

A scalar product needs no dictionary: nanobind's own casters carry it across
by value, exactly as they already carry a class's scalar members, and ROOT
has no TClass for a fundamental type.  Those are kept out of the LinkDef,
since rootcling has nothing to generate for them and says so once per rule.

Three things the build gets right rather than silently: the dictionary depends
on every header the product header reaches, taken from the compiler, so an
edit to a class it includes regenerates the bindings; it is linked against the
product libraries, because ROOT loads it on its own and needs the typeinfo and
the vtable of every class it describes -- which is also why the ROOT steps
carry the external library directories themselves, rather than relying on the
environment a job is run in; and a generator that fails stops the build with
its own message instead of leaving the previous output in place to be compiled
against.

ROOT is a build-time dependency only: rootcling and one macro, with nothing
linking against it, so a system installation is used rather than downloading
one.  ROOT_HIST=0 keeps the interpreter from writing to $HOME.

Three shared objects, because there are three different things to keep apart.
The generated bindings name every product type, so they are compiled into
edm_core, a Python extension in the library directory that links against the
product libraries and sits above them; the interpreter imports it by name, and
PythonRuntime puts that directory on sys.path first.  Nothing under Framework/
names a product -- checked rather than asserted: libFramework has no undefined
symbol mentioning a format.  nanobind's amalgamation is compiled exactly once,
into libnanobind.so -- two copies would each keep their own type registry and
casts would stop working across the boundary -- which is what NB_SHARED and
NB_BUILD are for: the first turns its core API from hidden symbols into
exported ones, the second marks the one object that defines them.  And
libFramework uses nanobind without owning a module of its own: each shared
object keeps its own pointer to the backend state, filled when *that* object
initialises a module, so Framework's would stay null and the first import from
PythonModuleRunner would segfault.  register_module() is the way in -- it joins
the domain edm_core created rather than making a second one.

Everything that touches nanobind is behind PythonModuleRunner's pimpl, so no
plugin translation unit includes it.  Nothing requires that -- with the core API
exported, a plugin could use nanobind directly -- but it keeps a future
Python-backed module in any plugin costing nothing extra.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
Five modules and square.ini: RandomGenerator publishes a vector of random
floats, SquareCxx squares it, a PythonProducer running python/square.py
squares the same vector with numpy, Compare checks the two element by
element and Printer reports.

Both squarers publish a std::vector<float> and are told apart only by the
label of the module that produced them, which is what the (type, label) key
was added for; fwtest could not express this workflow at all.

Compare runs with tolerance = 0 and strict = true, so the job fails if the
two implementations disagree by a single bit.  They do not: x*x is one IEEE
multiplication either way, and the Python side is a view of the same C++
storage rather than a copy of it, so there is nothing in between to round.

Comparison is a user-defined class rather than a tuple of scalars on
purpose: it is what the rootcling/TClass step reflects into bindings, so the
generator is exercised by this backend and not only by a real workflow.

emplace.ini exercises the other way a Python module publishes: emplace_demo.py
builds a std::vector<float> from a list and a pytest::Comparison from four
numbers, the C++ Printer reads the aggregate knowing nothing about where it
came from, and emplace_check.py reads the vector back.  square.py uses
allocate() instead, which is the zero-copy path -- it fills the product's own
storage through a numpy view -- so between them the two configurations cover
both.

square.py also reads the EventSetup, by type and with no token: IntESProducer
puts a bare int there, and a fundamental type has no TClass, so it crosses
through nanobind's own caster exactly as a scalar member of a product does.
pytest has had an EventSetup since it forked fwtest and no Python module had
ever touched it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
No binding in this backend is written by hand, which makes the generator the
thing that has to be right: a mistake in tools/generate_bindings.C is a mistake
in every product at once.  What there was to check it with was a
std::vector<float> and a four-field aggregate -- a scalar member and a
container, a small corner of what it emits.

Two products carry the shapes a real format has, with none of the content.
pytest::Calibration is conditions, read both one channel at a time -- a method
taking an argument, a class member, a fixed-size array, a static constant --
and as columns.  pytest::Samples is a structure of arrays with std::span
columns, a view handed out by reference and a method with a side effect, built
by a constructor that wires it to the conditions and to a column.

bindings.ini runs CalibrateCxx and python/calibrate.py on the same input and
the same conditions, compares the two products column by column with no
tolerance, and reads both back from a second Python module.  The arithmetic is
one subtraction and one multiplication, never a multiply-add that
-march=native could contract into an FMA, so numpy and C++ agree bit for bit.

The three interface paths nothing had used are covered here too -- put(),
emplace() of a scalar, and reading a scalar product back -- and so are the two
that have to fail: allocate() of a scalar, put() of something that cannot be
copied.

One of these is a regression test with a story behind it.  A method bound as
[](X& o) { return o.method(); } deduces its return type by value, so a method
returning a reference is bound to a copy.  Reading through the copy gives the
right answer -- the columns are pointers and still land in the right memory --
so only a write to the view itself shows the difference.  samples.view()
followed by view.filled = n is that write, and SamplesCompare fails on the
first event if it goes into a temporary.

The configurations are the tests.  test_cpu ran ./pytest --maxEvents 2 with no
configuration at all, which has not been able to do anything since the schedule
started coming from a file; it now runs all four.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
Framework/, bin/, DataFormats/, bindings/, nanobind/, tools/ and the plugin
that runs a Python module, under src/pyserial.  No test modules and no configurations: this is
the framework, and what it is for arrives in the next commit.

The fork is of pytest rather than of serial on purpose.  What pyserial needs
first is the framework -- labelled products, the configuration language, the
Python modules and the generated bindings -- and only then serial's
reconstruction on top of it.  Taking it by copying the directory keeps the
two Framework/ trees byte-identical, which is the property that stops a fix
in one from having to be re-applied by hand to the other.

The result builds and links; it cannot run anything, because a job needs a
configuration file and there is none.  ./pyserial says so and stops.

std::vector<float> is in the product list because a backend cannot have an
empty one: rootcling rejects a LinkDef with no selection rule.  The
reconstruction's own products replace it in the next commit.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
All of serial's reconstruction -- CUDACore, CUDADataFormats, CondFormats,
Geometry, the remaining DataFormats, and the seven producer and validation
plugins -- with reco.ini driving the same path serial hard-codes, and
reco-validate.ini the same job with the count validation scheduled rather than
asked for on the command line.

The plugins are serial's, changed in one respect: each takes the configuration
section its constructor is now handed.  Three of them read a label out of it --
SiPixelRecHitCUDA takes `beamSpot` and `input`, CAHitNtupletCUDA and
PixelVertexProducerCUDA take `input` -- and pass it to consumes<T>() as it
stands, empty or not.  An empty label is "the only product of that type", so
neither of the two configurations here names anything; a job running two
implementations of a module side by side is what makes a label necessary, and
that is what these keys are for.

The backend Makefile becomes serial's -- it carries the CUDA test rules and
the preprocessed-source target that fwtest's does not -- with the Python and
reflection sections kept, and main.cc becomes serial's, which has
--histogram and no --transfer.  Everything under Framework/ and the rest of
bin/ is untouched and stays byte-identical to pytest's.

The two products reflected are PixelTrackHeterogeneous and
ZVertexHeterogeneous, in place of the placeholder the fork carried.  Their
bindings are generated like any others, though the SoA columns behind them
are private at this commit, so the generator skips them and warns; the next
commit is what needs them.

test_cpu runs the configurations, since a bare ./pyserial has no schedule to
run: reco-validate.ini today, and the comparisons as they arrive.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
vertex_finder.py is a port of gpuVertexFinder::Producer::make(): loadTracks,
clusterTracksByDensity, fitVertices, splitVertices, fitVertices again,
sortByPt2, with the same parameters.  reco-python.ini runs it in place of
PixelVertexProducerCUDA, reco-python-validate.ini validates that path, and
reco-compare.ini runs both while VertexCompare checks one against the other,
which is only expressible because products are keyed by label.  Those three
join reco.ini and reco-validate.ini to make five configurations; porting the
next module changes what they say rather than adding a sixth.

Over 200 events the two agree exactly: same vertex count, worst |dz| = 0,
and CountValidator passes on the Python path.  The two loops that are
inherently sequential are written as loops here too, because their results
depend on the order in which they scan; everything else is vectorised.

Nothing is copied at the boundary and no binding here is written by hand:
every attribute the module touches is generated from the ROOT dictionary of
the product it belongs to, and the vertex SoA is allocated inside the Event
and filled through views of its own storage.

That is what eigenSoA's storage being public buys.  It was private behind a
data() returning a bare pointer with no extent, which is not something a
zero-copy view can be built from, so the generator skipped it.  This is the
one place a data format is changed to suit the generator, and it should be
read as a property of this format: eigenSoA is ancient, and the SoA layouts
that replaced it in CMSSW expose each column as a std::span, which carries
its own extent and would need no change at all.

doc/PythonBackends.md records the decisions taken across all of this --
what the alternatives were, which was chosen, and why -- including the
measurements, which only exist now that there is something to measure.  It
lands whole rather than in pieces because much of it only becomes true at
the end; the reasoning for each step is also in that step's commit message.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
BeamSpotToPOD is one line of C++ -- it copies eleven floats from the
EventSetup into the Event -- which makes it the module to port when the
question is what the boundary costs rather than what an algorithm costs.

Measured with a job holding only the beam spot, 20000 events, one stream,
three repeats:

  BeamSpotToPOD, C++                     0.0477 ms/event
  PythonProducer, empty produce()        0.0477 ms/event   +0.000
  PythonProducer, beam_spot.py           0.0503 ms/event   +0.0027

The crossing is free: attaching a thread state, building the Event and
EventSetup arguments and calling into Python cost nothing above the scatter,
under half a microsecond.  The whole difference is the eleven reads and
eleven writes, about 120 ns each through getattr by name, 35 ns by direct
attribute access.  Three thousand times less than the vertex finder, so it
changes nothing there -- but it is why bulk data has to cross as a view: a
module touching 10^4 hits one field at a time would spend milliseconds in
nanobind's attribute protocol and measure that instead.

This is the first module that needs the EventSetup, which Python reaches by
type and with no token, exactly as a C++ module writes es.get<T>().  square.py
already reads the scalar case; this is the class case, and the beam spot
crossing it is what the numbers above measure.

BeamSpotCompare checks the two beam spots field by field.  They are copies of
the same EventSetup product, so the comparison is exact: 200 events, no
difference in any field, and the full chain with the Python beam spot passes
CountValidator.

The three Python configurations gain a section rather than a file:
reco-python.ini and reco-python-validate.ini run the beam spot in Python along
with the vertices, and reco-compare.ini gains beamSpotPy next to beamSpot,
which is what makes recHits name the beam spot it wants.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
pixel_rechits.py is a port of gpuPixelRecHits.h: the extent and charge
distribution of each cluster over its digis, then the CPE position and its
error, then the local-to-global transform corrected for the beam spot.  This
is the module that moves real data -- 48000-67000 digis and 12000-21000
clusters in, 12000-21000 rec hits out, three orders of magnitude more objects
than the vertex finder sees, and every column a numpy view over the C++
storage with nothing copied in either direction.

It is bit-exact.  Over 200 events and 3045548 hits no hit differs in any of
the thirteen columns and the worst absolute difference is 0, and the full
chain built on the Python hits passes CountValidator for tracks and vertices.
The port had room to drift: the position correction branches on cluster size
and edge pixels, and phi is a degree-7 polynomial approximation of atan2
rather than atan2, so the port had to be the same polynomial.

Measured against reco.ini, per event:

              1 stream    16 streams
  C++          25.66 ms     25.40 ms
  beam spot    +0.00        -0.00
  rec hits     +4.36        +5.38
  vertices     +9.02        +9.38

The rec hits cost less than the vertex finder while handling a thousand times
more data, which says what the cost is made of: not the payload, but how much
work stays at Python level.  The vertex finder has two loops that cannot be
vectorised; the rec hits are a segmented reduction and a page of array
arithmetic, thirty-odd numpy calls whatever the number of hits.  Payload
shows up in the scaling instead -- 23% more per event from 1 to 16 streams
against the vertex finder's 4% -- where sixteen streams each touching ~10 MB
of columns compete for memory bandwidth.

The detector geometry is gathered once per stream rather than per event: it
is conditions, and re-reading it would be 37000 scalar attribute reads, 1.3
ms per event by D20's measurement, to fetch constants.

Two things in the framework this needed.  TrackingRecHit2D gains a
constructor taking the CPE parameters by reference and the module offsets as
a column, because a generated binding cannot supply a uint32_t const*; which is what
allocate(token, ...) forwards its arguments to, as a product that must be sized
before it can be filled needs.  The per-layer offsets and the phi index
move into a buildIndex() the C++ producer now calls too: they are bookkeeping
over hits that already exist, and the phi index is a bucketed structure that
nobody should rewrite in numpy.

The formats had to say how long their columns are.  Each keeps its storage as
a unique_ptr<T[]> with the length in a separate member and hands it out as a
bare pointer, which cannot be bound: nine columns of digis and clusters and
thirteen of rec hits were skipped for want of an extent.  Each now also offers
the column as a std::span, which the generator turns into a numpy view; the
private storage and the pointer accessors are untouched, and C++ callers see
no difference.

Binding every public method also instantiates every public method, and one had
a latent bug: HistoContainer::countDirect and friends compare a signed bin
against an unsigned nbins(), which is a -Werror=sign-compare.  Nothing had ever
instantiated them for the signed phi binner, so the warning had nowhere to
appear.

The three Python configurations gain a section again: reco-python.ini and
reco-python-validate.ini build the rec hits in Python, and reco-compare.ini
gains recHitsPy and RecHitCompare, so all three ported modules are now compared
against their C++ counterparts in one job.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
620 runs: threads = streams = 1..31, pinned with taskset to that many distinct
physical cores, five repeats per point with the first discarded, in two sets --
strong scaling (10000 events whatever the thread count) and weak scaling (1000
events per thread).

The Python chain runs at 0.66 of the C++ throughput and that ratio does not
move: 0.606 to 0.699 across the whole scan, with no trend, in both sets.  The
boundary does not get worse under concurrency, which is what free-threading had
to deliver.  Parallel efficiency at 31 threads is 0.91 for C++ and 0.92 for
Python.

Adds the four plots, all 620 measurements as a CSV, and D22 describing them.

The files are named for the scan rather than for the protocol -- scan-01.csv,
throughput-01-weak-scaling.png -- because "weak scaling" says how a job was run
and not what it ran, and there will be more than one weak-scaling scan here.
scan-01.md is the description that goes with them: which modules were C++ and
which were Python, how the job was pinned and repeated, and what each file
holds.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
pixel_clusters.py is a port of SiPixelRawToClusterCUDA and everything behind
it: the FED walk in the producer, RawToDigi, gpuCalibPixel's calibDigis,
gpuClustering's countModules, findClus and clusterChargeCut, and the capped
prefix sum that gives the rec hits their module offsets.

This is the first module of the chain, and the only one whose input is not
already a column: it reads the FED buffers as bytes.

The kernels are translated rather than rethought.  Where the C++ loops over
pixels this computes a column -- which is what numpy is for, and what every
other module here already does -- but the quantities, their order and the
branches are still the C++'s.  findClus in particular is the kernel's
algorithm: the pixels themselves as nodes, the neighbours of a pixel the ones
that come after it in the column histogram within one row, and the minimum
propagated over every pair every round until a round changes nothing, with the
odd rounds following each label to its root.  The histogram becomes a sort, the
neighbour list two searchsorted ranges per pixel, the iteration a reduction
over all the pairs at once and the numbering a cumulative sum; nothing about
which pixels end up in which cluster is decided differently, which is what
makes the products identical rather than merely close.

It agrees exactly.  Over 200 events: 1464992 digis in 45243 modules, with no
difference in any of the seven digi columns or the four per-module ones, and
reco-python-validate.ini then reconstructs the whole chain from the Python
clusters with CountValidator passing -- which is the check that matters, since
the digi, cluster, track and vertex counts now all come out of a Python chain
that starts at the raw data.

Making the inputs reachable took the same three things D21 needed for the
formats: a column as a pointer and a length rather than a bare pointer, a
reference rather than a pointer to a class, and a scalar accessor where the
declaration is something reflection cannot describe -- the gain calibration
keeps its ranges as an array of nested std::pairs.

One build-level detail came with them: allocate() learned to call a constructor
that takes only a size, which is what the digis and the clusters have -- they
are sized for the largest event the detector can produce and then filled.
Nothing else about these formats changes.  Their constructors and destructors
stay in the .cc files they were always in, which is what edm_core linking
against the product libraries is for; only the accessors the Python modules
read through -- the spans, the reference to the cabling map, the gain ranges --
are added to the headers.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
count_validator.py is a port of CountValidator: it reads the digi, cluster,
track and vertex counts the source carries next to the raw data and checks what
the chain produced against them, exactly for the first three and within the
same tolerances for the last two.

It is the first module whose state does not belong to a stream.  A C++ module
says that with a static counter and an endJob(); a Python module now says it
the same way, since a script is imported once and its module-level names are
therefore shared by every stream's instance.  The runner calls endJob() if the
object has one -- most modules do not -- and the framework calls it on the
first stream's module only, as it does for a C++ one.

The locking is the one thing that is genuinely different.  A C++ module gets
std::atomic counters; on a free-threaded interpreter `+=` on a module-level
integer is not atomic either, so the updates go under a threading.Lock, held
for the counter update and nothing else.

`validation = true` puts it in the path: the source starts producing the counts
and the label `countValidator` is scheduled.  The label is what matters, not
the plugin behind it, so reco-python-validate.ini names a PythonProducer there
and the same flag runs the Python validator.

Over the same events it reports the same verdict and the same average
differences as the C++ one, to the digit.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
histo_validator.py is a port of HistoValidator: thirty-three histograms of the
digis, clusters, hits, tracks and vertices of every event, written to a file at
the end of the job.

The file is a format, so this either reproduces it exactly or it is no use.
That means reproducing SimpleAtomicHisto's binning, which has two quirks worth
keeping -- nbins + 2 bins with underflow and overflow at the ends, and a value
that scales exactly to nbins folded back into the last real bin rather than
into the overflow -- and ostream's default formatting of a float, which "%g"
matches.

Both are a line of numpy over a whole column at a time, and np.bincount is what
turns a column of bin indices into counts; the per-value fill loop of the C++
module does not appear at all.

The two files are byte for byte identical.

The counters are the job's rather than a stream's, as in the count validator
and for the same reason: every stream fills the same histograms.

Neither validator is on the measured path -- reco.ini and reco-python.ini run
no validation -- which is what keeps a throughput difference between those two
the reconstruction and nothing else.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
Weak scaling at 1, 11, 21 and 31 cores, the protocol of D22's second set: the
Python chain runs at 0.33 to 0.36 of the C++ throughput, still flat, down from
0.66 with three modules in Python.

All of the difference is the clusterizer: three Python modules cost about 13 ms
per event between them, and the clusterizer alone about 35 ms more, most of it
the labelling -- the one algorithm in this chain that is a graph rather than an
array.

The measurements are scan-02.csv, with scan-02.md alongside them saying what
was run: the same five-module schedule as scan-01 with the clusterizer moved
across too, so the only C++ producer left in the measured arm is the track
finder.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
The module of the previous commit translates findClus; this one replaces it,
and the two are kept side by side so that the difference between them is the
algorithm and nothing else.

pixel_clusters_optimised.py groups the pixels into clusters with two ideas the
kernel does not have:

  - a run of touching pixels within one column is one cluster by construction,
    so the labelling can take runs as its nodes.  There are only about a
    quarter fewer of them -- 1.38 pixels to a run, a cluster being longer along
    the column than across it -- but what they remove is exactly the chains the
    minimum would otherwise walk one edge at a time, and the rounds drop from
    27 an event to 12;
  - those rounds are wildly uneven.  Three of them settle all but a few dozen
    runs of forty thousand and a stubborn cluster drags the rest out to
    nineteen, each still reducing over every edge.  A run's label can change
    only if a neighbour's did, so each round hands the next one only the runs
    next to those that moved.

    Two details decide whether that last one is a gain at all: the front is
    collected with a mask rather than np.unique, whose sort costs more than the
    round it saves, and a round whose front is still wider than a quarter of
    the runs reduces over everything anyway, since laying out a subset is
    itself work.  The version that did neither was 25% slower than the one it
    replaced.

Everything else is one copy of the code, in pixel_clusters_common.py: the FED
walk, the decoding, the calibration, the module boundaries, the charge cut, the
numbering and the products.  Neither module is written in terms of the other --
what they share is a third file, not a base class -- so a fix to any of it
reaches both, and the measurement is of the labelling alone:

  the labelling      35.1 ms   against   8.7 ms
  the whole module   44.8 ms            17.7 ms

Both are exact against SiPixelRawToClusterCUDA over 200 events.

reco-optimised.ini, reco-optimised-validate.ini and reco-optimised-compare.ini
are the three Python configurations with this module in place of the other; the
chain runs at 0.49 to 0.54 of the C++ throughput with it, against 0.35 to 0.38
with the translation.

The weak-scaling numbers for this module are scan-03.csv, with scan-03.md
alongside them.  It is scan-02 with one file swapped, and it ran earlier the
same morning, so each carries its own C++ control and the two are compared as
ratios to those rather than to each other.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
D27 claimed the runs and the active set are worth much less in C++ than in
numpy.  This measures it.

The C++ clusterizer is 4.44 ms an event, of which findClus is 2.363 -- the same
shape as the Python module, where the labelling is 8.7 ms of 17.7.  Inside
findClus: 0.69 ms of histogram, 0.61 of neighbour list, 0.78 of `min` iteration
and 0.33 of numbering, over 1812 modules and 60351 pixels.

The number that settles it is the round count: 5580 rounds an event over 1812
modules, 3.1 per module.  The kernel works one module at a time, so it
converges in three passes over thirty-odd pixels and has no tail for an active
set to skip -- while the same labelling over one graph of the whole event takes
27 rounds translated and 12 with the runs, because a batched graph's round
count is set by the worst cluster anywhere in it.  Batching is not optional in
numpy, so the reshaping buys back a cost the translation itself introduced.

What is left for the runs in C++ is 28% fewer nodes, bounded by 0.4 ms of the
2.36 and paid for by sorting each column by row, which the kernel does not
currently need.  Not worth doing.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
D22 scanned the whole thread range but with the clusterizer still in C++; D26
and D27 had the clusterizer in Python but only four points and three repeats.
Neither measures what the backend now is.  This is that scan: reco.ini against
reco-optimised.ini -- every ported module in Python, only the track finder left
in C++ -- weak scaling, 1000 events per thread after 1000 of warm up, 1 to 31
threads pinned to that many physical cores of one socket, five repeats per
point with the first discarded.  310 runs, no failures.

The chain runs at 0.501 of the C++ throughput, and the ratio does not move:
0.472 to 0.539 over the thirty-one points, with a fitted slope of -0.0001 per
thread against a scatter of 0.014.  That is D22's result at a lower value, and
it is the property that matters -- a fourth Python module costs a constant
factor, not a growing one.

Two things four points could not show.  The one-thread point is the odd one
out, and it belongs to the C++ arm: per-thread throughput *rises* from one
thread to two, 16% for C++ against 6% for Python, because at one stream nothing
overlaps.  That single artefact is the whole spread of the ratio, and it is
also what D22 records as a turbo bump -- which it cannot be, since turbo does
the opposite.  D22 is corrected here.

And measured from any other baseline the two arms degrade alike: 21.5% against
20.1% from two threads, 16.5% against 16.4% from four.  The tempting
explanation for the apparent difference is memory, so that is measured rather
than asserted -- peak RSS, and cache counters differenced between two event
counts so that startup cancels.  At 31 threads the C++ job holds 4.25 GB
against Python's 2.32 GB, so the Python arm uses less memory, not more;
last-level cache misses per event grow 9.7% in C++ and 9.6% in Python; IPC
falls 14.0% and 15.9%.  What is true is only that the Python chain touches six
times the cache references and three times the misses per event: more memory
traffic, not more sensitivity to sharing the machine.

scan-04.csv, the two plots and scan-04.md, which carries the module table, the
protocol and those counters.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WP5x6r6bsRbe8kEQEUJk1S
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant