from #149

Small problem I see is for this sequence:

>>> seq = MDASequence(axis_order="gtc",time_plan={"interval": 1, "loops": 2},channels=[{"config": "DAPI", "exposure": 1}],grid_plan={"rows": 2, "columns": 2})
>>> for event in seq:
...    print(event)
...
index=mappingproxy({'t': 0, 'g': 0, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=0.0 x_pos=-0.5 y_pos=0.5
index=mappingproxy({'t': 1, 'g': 0, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=1.0 x_pos=-0.5 y_pos=0.5
index=mappingproxy({'t': 0, 'g': 1, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=0.0 x_pos=0.5 y_pos=0.5
index=mappingproxy({'t': 1, 'g': 1, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=1.0 x_pos=0.5 y_pos=0.5
index=mappingproxy({'t': 0, 'g': 2, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=0.0 x_pos=0.5 y_pos=-0.5
index=mappingproxy({'t': 1, 'g': 2, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=1.0 x_pos=0.5 y_pos=-0.5
index=mappingproxy({'t': 0, 'g': 3, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=0.0 x_pos=-0.5 y_pos=-0.5
index=mappingproxy({'t': 1, 'g': 3, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=1.0 x_pos=-0.5 y_pos=-0.5

Might no be a very typical setup with axis order "gtc, but note how the min_start_time resets for every g index. This leads to the first time series at position 0 respecting the interval, while all later positions are acquired at max speed.

as pointed out by @henrypinkard in micro-manager/pycro-manager#266 (comment), the concept of autofocus should be somewhere

Should we add something like requested_destination to MDASequence?

I don't (immediately) want to add a whole lot of specific saving stuff, but a single field requesting a destination path might be very useful for possible downstream writers: It would provide a directory, base name, and extension. Writers could then increment that base name as desired (until any further spec added instructions on that sort of thing)

@ianhi ?

Hi @tlambert03

I really like the idea of useq-schema!

Over the last couple of days I've been working on putting together a basic C++ library that uses a combination of:

• MMCore I/O handling.
• OpenCV basic image processing.
• TensorRT inference on NVIDIA GPUs and deep learning accelerators.
• Imgui bloat-free Graphical User interface for C++ with minimal dependencies.
• Polyscope viewer for 3D data like meshes and point clouds extracted by OpenCV and TensorRT.

I primarily work in in situ sequencing method development (https://www.biorxiv.org/content/10.1101/722819v2).

One of the things we will see this year in in situ sequencing is an explosion of different commercial vendors offering "hardware boxes" that are essentially just a box with automated fluidics and an epifluorescent microscope to perform some iterative smFISH.

My idea with embarking on the above C++ library is rather that any microscope and hardware can be turned into a sequencing machine operating in a similar fashion to the Illumina Local Run Manager (keeping similar file conventions, user interface etc.).

I really want to use useq-schema to specify and make acquisition runs reproducible. Before embarking on this and digging deeper into the useq-schema I had one question:

Question: Do you have any current plan on how closed-loop events could be handled within the useq-schema. In situ sequencing is full of closed-loop events. Let me give you an example:

Closed-loop histogram equalization: Traditionally each nucleotide can be represented by a single non-overlapping fluorophore. For example A, T, C, G could be represented by Alexa488, Alexa555, Alexa594, Alexa647. Before each sequencing cycle is started the microscope first goes to a region not imaged and takes some snapshots with a given exposure time and light power. The images from these snapshots are then used to compute histograms and extract statistics that are then used to inform the exposure time and light power etc. of each channel setting such that the histograms of each channel are as closely equalized as possible before starting the true acquisition.

There are tons of similar events like this (checking that fluorophores have cleaved before starting the next cycle, adjusting laser power across z-axis etc.).

Do you have any good idea if events like this should already be represented at the level of useq-schema?

I'm wondering whether the autofocus plan has leaked a bit too much implementation details into the model.

I think in the vast majority of cases:

• there will only be one (or no) autofocus device.
• in most cases, I'm guessing that that device will generally be "auto-detectable" by the engine consuming the event (for example, i suspect we could auto-determine what the autofocus device was for micro-manager in the engine for most if not all supported device adapters)

So, that leads me to wonder whether we should update the autofocus plan to be, at minimimum, nothing more than a boolean that says "do autofocus for these axes", allowing autofocus_device_name: str | None. If None, it would imply "do your best to find it".

This thought comes from looking at pymmcore-plus/pymmcore-widgets#201 and wishing we didn't need the user to enter/select a string device name

It would be nice to be able to easily fill in an xarrray array's dims and coords based on a sequence. Or even to be able to directly create an xarray DataArrray from the sequence which will then be filled in via MDA.

Basic approach could be something like this

from typing import List, Dict
def sequence_to_xarray_dims_coords(sequence: MDASequence) -> [List, Dict]:
    mapping = {"t": "time_plan", "c": "channels", "z": "z_plan", "p": "positions"}
    coords = {}
    dims = []
    for ax in axis_order:
        try:
            if ax == "c":
                coords['C'] = [c.config for c in sequence.channels]
            else:
                coords[ax.upper()] = list(getattr(sequence, mapping[ax]))
            dims.append(ax.upper())
        except AttributeError:
            pass
    dims.extend(['Y', "X"])
    return dims, coords

It would be nice to in general accept numpy arrays instead of specifying lists. For example specifying a z_plan with a numpy array breaks:

import numpy as np
from useq import ZRelativePositions

# works
ZRelativePositions(relative=np.arange(-1.5, 1.5, 0.5).tolist())

# breaks
ZRelativePositions(relative=np.arange(-1.5, 1.5, 0.5))

ref: https://numpy.org/devdocs/reference/typing.html

Hi @tlambert03 any chance of a small version bump here to pick up the new metadata change?

tracking upstream issue at: pydantic/pydantic#7645

it's not helpful

The private _previous attribute on the AFAxesPlan has a bug: if the same AFplan object is used to iterate multiple times, the second time(s) will be different than the first. As a rather ugly workaround, we could have iter_sequence clear the _previous attribute...

@fdrgsp, we could add a new Tile specification to useq that would handle most of the logic needed for creating positions in the sample explorer. It could definitely be considered to be an additional dimension (since one could definitely imagine doing a large 3x3 image at every stage position)
I'm thinking something like:

class MDASequence:
    tile_plan: AnyTilePlan  # or a different name... perhaps grid... not sure

class TilePlan(FrozenModel):
    overlap: float = 0.1

class TileBoundingBox(TilePlan):
    x1: float
    x2: float
    ...

class TileRelative(TilePlan):
    nrows: int
    ncols: int
    relative_to: Literal['center', 'top_left', ...]

AnyTilePlan = Union[TileBoundingBox, TileRelative]

etc... basically, thinking of all of the non-overlapping ways to describe a tile in space, and what the bare minimum set of parameters are. @fdrgsp, would you want to take a stab at that? following the patterns in Zplan and TimePlan?

axis naming is a hard issue. it's a necessary evil. There are obvious conventions out there (e.g. XYZCT) but extensions are always needed (e.g. P position, G grid, row/col etc...). Ultimately useq-schema should only care about mapping an arbitrary string to an iterable of coordinates for that dimension. This is mostly what MDASequence does, albeit with hard-coded assumptions about the dimensions for now, so as to make it easy to accomplish the vast majority of use cases. The underlying code still just treats it as an iterable (mostly arbitrary) axes. but more could be done

• first, we definitely need to get rid of axis_order as a string. it should be Sequence[str], not assuming single-character keys #138
• we should make an additional more flexible class

  class NewSeq:
      axes: dict[str, Iterable[coord]]

  MDASequence is essentially a subclass of that more general pattern, with "known" axes PGZCT

we can't use pydantic's extra='forbid' setting, since it makes it very hard to add new fields without breaking backwards compatibility ... but it would be nice to have warnings to show that a kwarg is unrecognized

Very handy when trying to prepare a datastructure for receiving MDA data as noted here: https://github.com/tlambert03/pymmcore-plus/pull/29#discussion_r679267817

A quick function that gives this is:

    def get_axis_order(seq: MDASequence) -> Tuple[str]:
        """Get the axis order using only axes that are present in events."""
        event = next(seq.iter_events())
        event_axes = list(event.index.keys())
        return tuple(a for a in seq.axis_order if a in event_axes)

as pointed out by @henrypinkard in micro-manager/pycro-manager#266 (comment),

In the Channel object in MDASequence, dimension decimation (acquire every n...) should be generalized to arbitrary axes. one possible (but not ideal) pattern:

class Channel:
    ...
    decimate_axes: Dict[str, int]
    offset_axes: Dict[str, int]

# e.g
{'config': 'DAPI', 'decimate_axes': {'t': 5}}  # acquire DAPI every 5th frame
{'config': 'DAPI', 'decimate_axes': {'z': -1}}  # or something like that ... e.g. don't acquire the full stack
 {'config': 'DAPI', 'offset_axes': {'z': 0.5}}  # focal shift of 0.5 in this channel

There should be a way to specify something like Keep Shutter Open for Z stacks like in micro-manager.

In my use case of combining single cell measurements (with a laser) with traditional fluoresence microscopy I'm finding it awkward to incorporate when to collect a laser measurement into the current framework. What I'd like to do is (per position/time point):

• Collect BF z stack
• Aim my laser to a few (~5) points and collect spectra
• Collect GFP z-stack

I could add a raman channel - but there's no easy way to specify that that channel should not be z-stacked, and critically this will change the mda.shape which napari-micromanager relies on display the images. So instead maybe we could add NonImageMDAEvent that can be processed like an MDAEvent, but does not contribute to the shape of the sequence.

It should be possible to specify a time plan with no delay but some known duration of number of loops (e.g. time_plan=TIntervalDuration(interval=0, duration=3)) ... and it is possible to instantiate it... but when you iterate over it in an MDASequence, you'll get an error:

File ~/dev/self/useq-schema/src/useq/_time.py:72, in TimePlan.deltas(self)
     70 def deltas(self) -> Iterator[datetime.timedelta]:
     71     current = timedelta(0)
---> 72     for _ in range(self.loops):  # type: ignore  # TODO
     73         yield current
     74         current += self.interval

File ~/dev/self/useq-schema/src/useq/_time.py:145, in TIntervalDuration.loops(self)
    143 @property
    144 def loops(self) -> int:
--> 145     return self.duration // self.interval + 1

ZeroDivisionError: integer division or modulo by zero

For this to work, I think the consumer (such as MDARunner.run) might actually need to special case an MDASequence with a TIntervalDuration plan of interval 0... or we should add something to MDAEvent like a kill_time: a total duration after which the experiment should be stopped.

in pycromanager, there is an event key called keep_shutter_open that instructs the hardware to keep the shutter(s?) open across C/Z. should add that concept somewhere to the schema

see discussion here micro-manager/pycro-manager#88

useq-schema 0.4.7

I think the sequence setup in the documentation does not work anymore:

>>> from useq import MDASequence, Position, Channel, TIntervalDuration
>>> seq = MDASequence(
...     time_plan={"interval": 0.1, "loops": 2},
...     stage_positions=[(1, 1, 1)],
...     grid_plan={"rows": 2, "cols": 2},
...     z_plan={"range": 3, "step": 1},
...     channels=[{"config": "DAPI", "exposure": 1}]
... )

Works for me if I put "columns" instead of "cols".

Another small problem I see is for this sequence:

>>> seq = MDASequence(axis_order="gtc",time_plan={"interval": 1, "loops": 2},channels=[{"config": "DAPI", "exposure": 1}],grid_plan={"rows": 2, "columns": 2})
>>> for event in seq:
...    print(event)
...
index=mappingproxy({'t': 0, 'g': 0, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=0.0 x_pos=-0.5 y_pos=0.5
index=mappingproxy({'t': 1, 'g': 0, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=1.0 x_pos=-0.5 y_pos=0.5
index=mappingproxy({'t': 0, 'g': 1, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=0.0 x_pos=0.5 y_pos=0.5
index=mappingproxy({'t': 1, 'g': 1, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=1.0 x_pos=0.5 y_pos=0.5
index=mappingproxy({'t': 0, 'g': 2, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=0.0 x_pos=0.5 y_pos=-0.5
index=mappingproxy({'t': 1, 'g': 2, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=1.0 x_pos=0.5 y_pos=-0.5
index=mappingproxy({'t': 0, 'g': 3, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=0.0 x_pos=-0.5 y_pos=-0.5
index=mappingproxy({'t': 1, 'g': 3, 'c': 0}) channel=Channel(config='DAPI') exposure=1.0 min_start_time=1.0 x_pos=-0.5 y_pos=-0.5

Might no be a very typical setup with axis order "gtc, but note how the min_start_time resets for every g index. This leads to the first time series at position 0 respecting the interval, while all later positions are acquired at max speed.

I'm thinking about how to implement a keep shutter open option for pymmcore-plus and finding that there is not a convenient way to interpret the for event in seq.iter_events() to check for how many events we should keep the shutter open.

I think there needs to be a way to implement generating a chunk of events that match some criteria - e.g. give me the next set of events that have channel=='FITC'.

One way to do this maybe would be to introduce a new object MDAEventGroup that could serve as the intermediate, or perhaps to allow MDASequences to contain MDASequences.

The simpler workaround that I'm currently going with is to expand on the iter_ methods with a new iter_chunked method that goes like so:

def iter_chunked(seq: MDASequence, over: str) -> Iterator[List[MDAEvent]]:
    """
    Generate events in chunks where only the given *over* index changes

    Parameters
    ----------
    seq : MDASequence
    over : str
        The index keys over which to allow grouping. E.g. 'z' to return
        all events in the z stack for a given time point, position and channel.

    yields
    ------
    chunk : list of events
    """

    over = over.lower()

    def extract_index(event):
        return tuple(event.index.get(key) for key in event.index.keys() if key not in over)

    events = seq.iter_events()
    e = next(events)
    cur_chunk = extract_index(e)

    chunk = [e]
    for e in events:
        if extract_index(e) == cur_chunk:
            chunk.append(e)
        else:
            yield chunk
            chunk = [e]
            cur_chunk = extract_index(e)
    if len(chunk) > 0:
        # ensure we always return all events
        yield chunk