Typed experiment data, with the spectrum first

A measurement is rarely one array. A spectrum is an intensity and the wavelength axis it was measured on, plus the exposure time, the grating and the file it came from. Until this month a PLESTY module passed those around as loose arrays and dictionaries. plesty-lib now has one contract for such data, called ExperimentData, and ships its first real type, SpectrumData, together with an analyzer that fits every emission line in a stack of spectra.

What changed

ExperimentData (in plesty.lib.data) is an abstract base class: a Python class that other classes build on and that says which methods they must provide. Every typed piece of data built on it carries:

• a kind: a unique lower-case name such as "spectrum". It is the key under which the class is registered, the name written into every saved file, and a type name a schema can use;
• metadata: the acquisition context (exposure, grating, source file, the processing steps applied so far);
• provenance: who produced the data, meaning the device, or the analyzer and the settings it ran with.

With that contract, the same object travels through the whole PLESTY pipeline. Device modules talk to hardware, experiments orchestrate a measurement, and analyzers turn raw data into results:

The work landed in three steps on plesty-lib's development line during September: the contract with SpectrumData, then the SpectrumFit analyzer, then quality diagnostics for the fit. It is available in the 0.5.0.dev snapshot releases.

Using it

1. Load a spectrum, from a file or a whole run

SpectrumData holds a 1-D wavelength axis (default unit nm) and an intensity that is either one spectrum or a stack of spectra, one per step of a scan. Each spectrum in a stack is called a frame. Optional coordinates record the value that was scanned for each frame, such as a waveplate angle.

from plesty.lib.data import ExperimentData, SpectrumData

# One file: LightField .spe 3.0, LightField CSV exports and plain tables are read
spectrum = SpectrumData.load("qd1.spe")

# A whole experiment run: one frame per step that saved a spectrum,
# with the scanned values of that step kept as coordinates
spectrum = SpectrumData.from_run(run_dir, coordinates=("hwp_deg", "power_w"))
spectrum.coordinates["hwp_deg"]      # one angle per frame, unit deg

A coordinate takes its unit from the end of its name: _deg, _w, _nm, _s and so on.

2. Prepare it: every step is recorded

Processing methods return a new spectrum and append the step to metadata["processing"], so a saved result always says how it was prepared.

prepared = spectrum.crop(770, 780).subtract_background("median")
prepared.find_peaks(prominence=100)       # locate lines, NumPy only
prepared.energy()                         # the axis in eV
prepared.save("qd1.npz")                   # lossless; .json, .csv and .txt also work

same = ExperimentData.load("qd1.npz")     # a file plesty wrote names its own type: comes back as a SpectrumData

3. Fit every line in every frame with SpectrumFit

SpectrumFit is an analyzer that ships with plesty-lib. It finds the emission lines, fits each one on a background in every frame, and returns positions, widths, heights and areas with their uncertainties. It needs SciPy, which comes with the fit extra:

uv add "plesty-lib[fit]"

from plesty.lib.analyzer import SpectrumFit

outputs = SpectrumFit(prominence=250, window=0.05, center_tolerance=0.1)(spectrum=spectrum)
fit = outputs["fit"]                      # SpectrumFitData
fit["center"], fit["center_err"]         # arrays of shape (frames, lines)
fit.peak(0)                               # one line across all frames
fit.coordinates["hwp_deg"]               # the scanned angle of each row
outputs["curve"]                         # the fitted model as a SpectrumData

How the fit works, in short:

1. Lines are found in the average of the frames (peaks="auto"), or taken from positions you give. average="median" keeps a cosmic ray that hits a single frame from being mistaken for a line.
2. Groups: lines closer than overlap times their width are fitted together.
3. Tracking: with track=True (the default) each frame starts from the previous frame's result, so a line that drifts across a scan is followed.
4. Models: gaussian, lorentzian or voigt, on a linear or constant background.

Building on it

Add your own data type

A new kind of measurement (an image stack, a time trace, a power curve) becomes a first-class type by subclassing ExperimentData. You declare a kind and implement three methods. Saving, loading, the registry and the schema type come with the base class.

from plesty.lib.data import ExperimentData, PlestyArray

class PowerCurve(ExperimentData):
    kind = "power_curve"            # unique, lower-case; registered on class creation
    schema_version = 1              # bump when the arrays or their meaning change
    # __init__(power, current, metadata, provenance) omitted: it calls
    # super().__init__(metadata, provenance) and stores the two arrays

    def to_arrays(self) -> dict[str, PlestyArray]:
        return {"power": self.power, "current": self.current}

    @classmethod
    def from_arrays(cls, arrays, metadata, provenance):
        return cls(arrays["power"], arrays["current"], metadata, provenance)

    @classmethod
    def example(cls):                # small valid instance for docs, mocks and tests
        ...

Take typed data in an analyzer

An analyzer declares its inputs and outputs in a schema. A registered kind is now a valid type there, and the analyzer framework checks the type before your code runs:

class LineShift(Analyzer):
    input_schema = {"spectrum": {"dtype": "spectrum", "description": "Spectrum stack."}}
    ...

A SpectrumData returned by an analyzer keeps the provenance it already had and gains the analyzer's name, version and settings. Settings must therefore be constructor arguments, never module-level constants.

Fit something other than a spectrum stack

The fitting functions under SpectrumFit are public. An analyzer that fits, for example, a line centre across a polarization sweep can call fit_spectrum(x, y, [Peak(...), ...], model, background, sigma=...), which fits several lines with bounds and per-point errors. fit_peak is the one-line form. Both return a SpectrumFitResult with each line's centre, width, height and area, each with its error, plus r_squared and evaluate(x).