.gitlab-ci.yml

@@ -67,7 +67,7 @@ Static Type:
     - poetry run mypy . --exclude 'setup\.py$'
     # We can remove the flag once this is resolved https://github.com/pypa/setuptools/issues/2345
 
-# TODO add more tests before activate automatic release
+# TODO add more tests before activating auto-release
 # Bump_release:
 #   stage: release
 #   script:

.gitlab/issue_templates/bug.md

 ## Summary
 
-(Summarize the bug encountered concisely)
+{Summarize the bug encountered concisely}
+
+I confirm that I have (if relevant):
+- [ ] Read the troubleshooting guide: https://gitlab.com/aliby/aliby/-/wikis/Troubleshooting-(basic)
+- [ ] Updated aliby and aliby-baby.
+- [ ] Tried the unit test.
+- [ ] Tried a scaled-down version of my experiment (distributed=0, filter=0, tps=10)
+- [ ] Tried re-postprocessing.
 
 ## Steps to reproduce
 
-(How one can reproduce the issue - this is very important)
+{How one can reproduce the issue - this is very important}
+
+- aliby version: 0.1.{...}, or if development/unreleased version, commit SHA: {...}
+- platform(s):
+    - [ ] Jura
+    - [ ] Other Linux, please specify distribution and version: {...}
+    - [ ] MacOS, please specify version: {...}
+    - [ ] Windows, please specify version: {...}
+- experiment ID: {...}
+   - Any special things you need to know about this experiment: {...}
 
 ## What is the current bug behavior?
 
@@ -19,6 +35,12 @@
 (Paste any relevant logs - please use code blocks (```) to format console output, logs, and code, as
 it's very hard to read otherwise.)
 
+```
+{PASTE YOUR ERROR MESSAGE HERE!!}
+
+
+```
+
 ## Possible fixes
 
 (If you can, link to the line of code that might be responsible for the problem)

README.md

@@ -11,15 +11,18 @@ End-to-end processing of cell microscopy time-lapses. ALIBY automates segmentati
 ## Quickstart Documentation
 Installation of [VS Studio](https://visualstudio.microsoft.com/downloads/#microsoft-visual-c-redistributable-for-visual-studio-2022) Native MacOS support for is under work, but you can use containers (e.g., Docker, Podman) in the meantime.
 
-For analysing local data
+To analyse local data
  ```bash
-pip install aliby # aliby[network] if you want to access an OMERO server
+pip install aliby 
  ```
-
+ Add any of the optional flags `omero` and `utils` (e.g., `pip install aliby[omero, utils]`). `omero` provides tools to connect with an OMERO server and `utils` provides visualisation, user interface and additional deep learning tools.
+  
 See our [installation instructions]( https://aliby.readthedocs.io/en/latest/INSTALL.html ) for more details.
 
 ### CLI
 
+If installed via poetry, you have access to a Command Line Interface (CLI)
+
  ```bash
 aliby-run --expt_id EXPT_PATH --distributed 4 --tps None
  ```
@@ -31,6 +34,8 @@ And to run Omero servers, the basic arguments are shown:
 
 The output is a folder with the original logfiles and a set of hdf5 files, one with the results of each multidimensional inside.
 
+For more information, including available options, see the page on [running the analysis pipeline](https://aliby.readthedocs.io/en/latest/PIPELINE.html)
+
 ## Using specific components
 
 ### Access raw data
@@ -76,12 +81,18 @@ It fetches the metadata from the Image object, and uses the TilerParameters valu
 ```python
 fpath = "h5/location"
 
-trap_id = 9
-trange = list(range(0, 30))
+tile_id = 9
+trange = range(0, 10)
 ncols = 8
 
 riv = remoteImageViewer(fpath)
-trap_tps = riv.get_trap_timepoints(trap_id, trange, ncols)
+trap_tps = [riv.tiler.get_tiles_timepoint(tile_id, t) for t in trange] 
+
+# You can also access labelled traps
+m_ts = riv.get_labelled_trap(tile_id=0, tps=[0])
+
+# And plot them directly
+riv.plot_labelled_trap(trap_id=0, channels=[0, 1, 2, 3], trange=range(10))
 ```
 
 Depending on the network speed can take several seconds at the moment.
@@ -91,8 +102,8 @@ For a speed-up: take fewer z-positions if you can.
 Alternatively, if you want to get all the traps at a given timepoint:
 
 ```python
-timepoint = 0
-seg_expt.get_tiles_timepoints(timepoint, tile_size=96, channels=None,
+timepoint = (4,6)
+tiler.get_tiles_timepoint(timepoint, channels=None,
                                 z=[0,1,2,3,4])
 ```
 

docs/source/INSTALL.md

@@ -62,7 +62,7 @@ For Windows, the simplest way to install it is using conda (or mamba). You can i
     $ \PATH\TO\POETRY\LOCATION\poetry install
 
   - MacOS
-  Under work (See issue https://github.com/ome/omero-py/issues/317)
+  For local access and processing, follow the same instructions as Linux. Remote access to OMERO servers depends on some issues in one of our depedencies being solved (See issue https://github.com/ome/omero-py/issues/317)
 
 ### Git version
 
@@ -71,9 +71,23 @@ Install [ poetry ](https://python-poetry.org/docs/#installation) for dependency
 In case you want to have local version:
 
     $ git clone git@gitlab.com/aliby/aliby.git
-    $ cd aliby && poetry install --all-extras
+    $ cd aliby
+    
+ and then either
 
-This will automatically install the [ BABY ](https://gitlab.com/aliby/baby) segmentation software. Support for additional segmentation and tracking algorithms is under development.
+    $$ poetry install --all-extras
+
+for everything, including tools to access OMERO servers, or
+
+    $$ poetry install
+
+for a version with only local access, or
+
+    $$ poetry install --with dev
+
+to install with compatible versions of the development tools we use, such as black.
+
+These commands will automatically install the [ BABY ](https://gitlab.com/aliby/baby) segmentation software. Support for additional segmentation and tracking algorithms is under development.
 
 ## Omero Server
 
@@ -111,3 +125,45 @@ docker-compose stop
 Segmentation has been tested on: Mac OSX Mojave, Ubuntu 20.04 and Arch Linux.
 Data processing has been tested on all the above and Windows 11.
 
+### Detailed Windows installation
+#### Create environment
+Open anaconda powershell as administrator
+```shell  script
+conda create -n devaliby2 -c conda-forge python=3.8 omero-py
+conda activate devaliby2
+```
+
+#### Install poetry
+    You may have to specify the python executable to get this to work :
+```shell script
+(Invoke-WebRequest -Uri https://install.python-poetry.org -UseBasicParsing).Content | C:\Users\USERNAME\Anaconda3\envs\devaliby2\python.exe -
+
+```    Also specify full path when running poetry (there must be a way to sort this)
+
+- Clone the repository (Assuming you have ssh properly set up)
+```shell script
+git clone git@gitlab.com:aliby/aliby.git
+cd aliby
+poetry install --all-extras
+```
+
+You may need to run the full poetry path twice - first time gave an error message, worked second time
+
+```shell script
+C:\Users\v1iclar2\AppData\Roaming\Python\Scripts\poetry install --all-extras
+```
+
+confirm installation of aliby - python...import aliby - get no error message
+
+#### Access the virtual environment from the IDE (e.g., PyCharm)
+New project
+In location - navigate to the aliby folder (eg c::/Users/Public/Repos/aliby
+
+- Select the correct python interpreter
+click the interpreter name at the bottom right
+click add local interpreter
+on the left click conda environment
+click the 3 dots to the right of the interpreter path and navigate to the python executable from the environment created above (eg C:\Users\v1iclar2\Anaconda3\envs\devaliby2\python.exe)
+
+#### Potential Windows issues
+- Sometimes the library pywin32 gives trouble, just install it using pip or conda 

docs/source/PIPELINE.md

+# Running the analysis pipeline
+
+You can run the analysis pipeline either via the command line interface (CLI) or using a script that incorporates the `aliby.pipeline.Pipeline` object.
+
+## CLI
+
+On a CLI, you can use the `aliby-run` command.  This command takes options as follows:
+- `--host`: Address of image-hosting server.
+- `--username`: Username to access image-hosting server.
+- `--password`: Password to access image-hosting server.
+- `--expt_id`: Number ID of experiment stored on host server.
+- `--distributed`: Number of distributed cores to use for segmentation and signal processing.  If 0, there is no parallelisation.
+- `--tps`: Optional.  Number of time points from the beginning of the experiment to use.  If not specified, the pipeline processes all time points.
+- `--directory`: Optional.  Parent directory to save the data files (HDF5) generated, `./data` by default; the files will be stored in a child directory whose name is the name of the experiment.
+- `--filter`: Optional.  List of positions to use for analysis.  Alternatively, a regex (regular expression) or list of regexes to search for positions.  **Note: for the CLI, currently it is not able to take a list of strings as input.**
+- `--overwrite`: Optional.  Whether to overwrite an existing data directory.  True by default.
+- `--override_meta`: Optional.  Whether to overwrite an existing data directory.  True by default.
+
+Example usage:
+ ```bash
+aliby-run --expt_id EXPT_PATH --distributed 4 --tps None
+ ```
+
+And to run Omero servers, the basic arguments are shown:
+ ```bash
+ aliby-run --expt_id XXX --host SERVER.ADDRESS --user USER --password PASSWORD 
+ ```
+
+
+## Script
+
+Use the `aliby.pipeline.Pipeline` object and supply a dictionary, following the example below.  The meaning of the parameters are the same as described in the CLI section above.
+
+```python
+#!/usr/bin/env python3
+
+from aliby.pipeline import Pipeline, PipelineParameters
+
+# Specify experiment IDs
+ids = [101, 102]
+
+for i in ids:
+    print(i)
+    try:
+        params = PipelineParameters.default(
+            # Create dictionary to define pipeline parameters.
+            general={
+                "expt_id": i,
+                "distributed": 6,
+                "host": "INSERT ADDRESS HERE",
+                "username": "INSERT USERNAME HERE",
+                "password": "INSERT PASSWORD HERE",
+                # Ensure data will be overwriten
+                "override_meta": True,
+                "overwrite": True,
+            }
+        )
+
+        # Fine-grained control beyond general parameters:
+        # change specific leaf in the extraction tree.
+        # This example tells the pipeline to additionally compute the
+        # nuc_est_conv quantity, which is a measure of the degree of
+        # localisation of a signal in a cell.
+        params = params.to_dict()
+        leaf_to_change = params["extraction"]["tree"]["GFP"]["np_max"]
+        leaf_to_change.add("nuc_est_conv")
+
+        # Regenerate PipelineParameters
+        p = Pipeline(PipelineParameters.from_dict(params))
+
+        # Run pipeline
+        p.run()
+        
+    # Error handling
+    except Exception as e:
+        print(e)
+```
+
+This example code can be the contents of a `run.py` file, and you can run it via
+
+```bash
+python run.py
+```
+
+in the appropriate virtual environment.
+
+Alternatively, the example code can be the contents of a cell in a jupyter notebook.

docs/source/index.rst

@@ -4,15 +4,14 @@
    contain the root `toctree` directive.
 
 .. toctree::
-   :hidden:
 
    Home page <self>
    Installation <INSTALL.md>
+   Pipeline options <PIPELINE.md>
+   Contributing <CONTRIBUTING.md>
    ..
       Examples <examples.rst>
-   Reference <api.rst>
-   ..
-      Contributing <CONTRIBUTING.md>
+      Reference <api.rst>
    ..
 
 .. include:: ../../README.md

docs/source/specifications/io_dependencies.org

+#+title: Input/Output Stage Dependencies
+
+Overview of what fields are required for each consecutive step to run, and
+
+- Registration
+- Tiler
+  - Requires:
+    - None
+  # - Optionally:
+  - Produces:
+    - /trap_info
+- Tiler
+  - Requires:
+    - None
+  - Produces:
+    - /trap_info

docs/source/specifications/roadmap.org

+#+title: ALIBY roadmap
+
+Overview of potential improvements, goals, issues and other thoughts worth keeping in the repository. In general, it is things that the original developer would have liked to implement had there been enough time.
+
+* General goals
+- Simplify code base
+- Reduce dependency on BABY
+- Abstract components beyond cell outlines (i.e, vacuole, or other ROIs)
+- Enable providing metadata defaults (remove dependency of metadata)
+- (Relevant to BABY): Migrate aliby-baby to Pytorch from Keras. Immediately after upgrade h5py to the latest version (we are stuck in 2.10.0 due to Keras).
+
+* Long-term tasks (Soft Eng)
+- Support external segmentation/tracking/lineage/processing tools
+  - Split segmentation, tracking and lineage into independent Steps
+- Implement the pipeline as an acyclic graph
+- Isolate lineage and tracking into a section of aliby or an independent package
+- Abstract cells into "ROIs" or "Outlines"
+- Abstract lineage into "Outline relationships" (this may help study cell-to-cell interactions in the future)
+- Add support to next generation microscopy formats.
+- Make live cell processing great again! (low priority)
+
+* Potential features
+- Flat field correction (requires research on what is the best way to do it)
+- Support for monotiles (e.g., agarose pads)
+- Support the user providing location of tiles (could be a GUI in which the user selects a region)
+- Support multiple neural networks (e.g., vacuole/nucleus in adition to cell segmentation)
+- Use CellPose as a backup for accuracy-first pipelines
+
+* Potential CLI(+matplotlib) interfaces
+The fastest way to get a gui-like interface is by using matplotlib as a panel to update and read keyboard inputs to interact with the data. All of this can be done within matplotlib in a few hundreds of line of code.
+
+- Annotate intracellular contents
+- Interface to adjust the parameters for calibration
+- Basic selection of region of interest in a per-position basis
+
+* Sections in need of refactoring
+** Extraction
+Extraction could easily increase its processing speed. Most of the code was not originally written using casting and vectorised operations.
+- Reducing the use of python loops to the minimum
+- Replacing nested functions with functional mappings (extraction be faster and clearer with a functional programming approach)
+- Replacing the tree with a set of tuples and delegating processing order to dask.
+  Dask can produce its own internal tree and optimise the order of rendering the tree unnecessary
+
+** Postprocessing.
+- Clarify the limits of picking and merging classes: These are temporal procedures; in the future segmentation should become more accurate, making picking Picker redundant; better tracking/lineage assignemnt will make merging redundant.
+- Formalise how lineage and reshaper processes are handled
+- Non-destructive postprocessing.
+  The way postprocessing is done is destructive at the moment. If we aim to perform more complex data analysis automatically an implementation of complementary and tractable sub-pipelines is essential. (low priority, perhaps within scripts)
+- Functionalise parameter-process schema. This schema provides a decent structure, but it requires a lot of boilerplate code. To transition the best option is probably a function that converts Process classes into a function, and another that extracts default values from a Parameters class. This could in theory replace most Process-Parameters pairs. Lineage functions will pose a problem and a common interface to get lineage or outline-to-outline relationships demands to be engineered.
+
+** Compiler/Reporter
+- Remove compiler step, and focus on designing an adequate report, then build it straight after postprocessing ends.
+
+** Writers/Readers
+- Consider storing signals that are similar (e.g., signals arising from each channel) in a single multidimensional array to save storage space. (mid priority)
+- Refactor (Extraction/Postprocessing) Writer to use the DynamicWriter Abstract Base Class.
+
+** Pipeline
+Pipeline is in dire need of refactoring, as it coordinates too many things. The best approach would be to modify the structure to delegate more responsibilities to Steps (such as validation) and Writers (such as writing metadata).
+
+* Testing
+- I/O interfaces
+- Visualisation helpers and other functions
+- Running one pipeline from another
+- Groupers
+
+* Documentation
+- Tutorials and how-to for the usual tasks
+- How to deal with different types of data
+- How to aggregate data from multiple experiments
+- Contribution guidelines (after developing some)
+
+* Tools/alternatives that may be worth considering for the future
+- trio/asyncio/anyio for concurrent processing of individual threads
+- Pandas -> Polars: Reconsider after pandas 2.0; they will become interoperable
+- awkward arrays: Better way to represent data series with different sizes
+- h5py -> zarr: OME-ZARR format is out now, it is possible that the field will move in that direction. This would also make us being stuck in h5py 2.10.0 less egregious.
+- Use CellACDC's work on producing a common interface to access a multitude of segmentation algorithms.
+
+* Secrets in the code
+- As aliby is adapted to future Python versions, keep up with the "FUTURE" statements that enunciate how code can be improved in new python version
+- Track FIXMEs and, if we cannot solve them immediately, open an associated issue
+
+* Minor inconveniences to fix
+- Update CellTracker models by training with current scikit-learn (currently it warns that the models were trained in an older version of sklearn )

examples/parsers/swainlab_logfile_header_example.log

+2022-10-10 15:31:27,350 - INFO 
+Swain Lab microscope experiment log file
+GIT commit: e5d5e33 fix: changes to a few issues with focus control on Batman.
+Microscope name: Batman
+Date: 022-10-10 15:31:27
+Log file path: D:\AcquisitionDataBatman\Swain Lab\Ivan\RAW DATA\2022\Oct\10-Oct-2022\pH_med_to_low00\pH_med_to_low.log
+Micromanager config file: C:\Users\Public\Microscope control\Micromanager config files\Batman_python_15_4_22.cfg
+Omero project: Default project
+Omero tags: 
+Experiment details: Effect on growth and cytoplasmic pH of switch from normal pH (4.25) media to higher pH (5.69). Switching is run using the Oxygen software
+-----Acquisition settings-----
+
+2022-10-10 15:31:27,350 - INFO Image Configs:
+Image config,Channel,Description,Exposure (ms), Number of Z sections,Z spacing (um),Sectioning method
+brightfield1,Brightfield,Default bright field config,30,5,0.6,PIFOC
+pHluorin405_0_4,pHluorin405,Phluorin excitation from 405 LED 0.4v and 10ms exposure,5,1,0.6,PIFOC
+pHluorin488_0_4,GFPFast,Phluorin excitation from 488 LED 0.4v,10,1,0.6,PIFOC
+cy5,cy5,Default cy5,30,1,0.6,PIFOC
+
+Device properties:
+Image config,device,property,value
+pHluorin405_0_4,DTOL-DAC-1,Volts,0.4
+pHluorin488_0_4,DTOL-DAC-2,Volts,0.4
+cy5,DTOL-DAC-3,Volts,4
+
+2022-10-10 15:31:27,353 - INFO 
+group: YST_247 field: position
+Name, X, Y, Z, Autofocus offset
+YST_247_001,-8968,-3319,2731.125040696934,123.25
+YST_247_002,-8953,-3091,2731.3000406995416,123.25
+YST_247_003,-8954,-2849,2731.600040704012,122.8
+YST_247_004,-8941,-2611,2730.7750406917185,122.8
+YST_247_005,-8697,-2541,2731.4500407017767,118.6
+group: YST_247 field: time
+start: 0
+interval: 300
+frames: 180
+
+group: YST_247 field: config
+brightfield1: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin405_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin488_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+cy5: 0xfffffffffffffffffffffffffffffffffffffffffffff
+
+2022-10-10 15:31:27,356 - INFO 
+group: YST_1510 field: position
+Name,X,Y,Z,Autofocus offset
+YST_1510_001,-6450,-230,2343.300034917891,112.55
+YST_1510_002,-6450,-436,2343.350034918636,112.55
+YST_1510_003,-6450,-639,2344.000034928322,116.8
+YST_1510_004,-6450,-831,2344.250034932047,116.8
+YST_1510_005,-6848,-536,2343.3250349182636,110
+group: YST_1510 field: time
+start: 0
+interval: 300
+frames: 180
+
+group: YST_1510 field: config
+brightfield1: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin405_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin488_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+cy5: 0xfffffffffffffffffffffffffffffffffffffffffffff
+2022-10-10 15:31:27,359 - INFO 
+group: YST_1511 field: position
+Name, X, Y, Z, Autofocus offset
+YST_1511_001,-10618,-1675,2716.900040484965,118.7
+YST_1511_002,-10618,-1914,2717.2250404898077,122.45
+YST_1511_003,-10367,-1695,2718.2500405050814,120.95
+YST_1511_004,-10367,-1937,2718.8250405136496,120.95
+YST_1511_005,-10092,-1757,2719.975040530786,119.45
+
+group: YST_1511 field: time
+start: 0
+interval: 300
+frames: 180
+
+group: YST_1511 field: config
+brightfield1: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin405_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin488_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+cy5: 0xfffffffffffffffffffffffffffffffffffffffffffff
+
+2022-10-10 15:31:27,362 - INFO 
+group: YST_1512 field: position
+Name,X,Y,Z,Autofocus offset
+YST_1512_001,-8173,-2510,2339.0750348549336,115.65
+YST_1512_002,-8173,-2718,2338.0250348392874,110.8
+YST_1512_003,-8173,-2963,2336.625034818426,110.8
+YST_1512_004,-8457,-2963,2336.350034814328,110.9
+YST_1512_005,-8481,-2706,2337.575034832582,113.3
+group: YST_1512 field: time
+start: 0
+interval: 300
+frames: 180
+
+group: YST_1512 field: config
+brightfield1: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin405_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin488_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+cy5: 0xfffffffffffffffffffffffffffffffffffffffffffff
+
+2022-10-10 15:31:27,365 - INFO 
+group: YST_1513 field: position
+Name,X,Y,Z,Autofocus offset
+YST_1513_001,-6978,-2596,2339.8750348668545,113.3
+YST_1513_002,-6978,-2380,2340.500034876168,113.3
+YST_1513_003,-6971,-2163,2340.8750348817557,113.3
+YST_1513_004,-6971,-1892,2341.2500348873436,113.3
+YST_1513_005,-6692,-1892,2341.550034891814,113.3
+group: YST_1513 field: time
+start: 0
+interval: 300
+frames: 180
+
+group: YST_1513 field: config
+brightfield1: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin405_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+pHluorin488_0_4: 0xfffffffffffffffffffffffffffffffffffffffffffff
+cy5: 0xfffffffffffffffffffffffffffffffffffffffffffff
+
+2022-10-10 15:31:27,365 - INFO 
+2022-10-10 15:31:27,365 - INFO 
+-----Experiment started-----

pyproject.toml

 [tool.poetry]
 name = "aliby"
-version = "0.1.55"
+version = "0.1.64"
 description = "Process and analyse live-cell imaging data"
 authors = ["Alan Munoz <alan.munoz@ed.ac.uk>"]
 packages = [
@@ -14,6 +14,8 @@ readme = "README.md"
 
 [tool.poetry.scripts]
 aliby-run = "aliby.bin.run:run"
+aliby-annotate = "aliby.bin.annotate:annotate"
+aliby-visualise = "aliby.bin.visualise:napari_overlay"
 
 [build-system]
 requires = ["setuptools", "poetry-core>=1.0.0"]
@@ -25,6 +27,7 @@ PyYAML = "^6.0"
 flatten-dict = "^0.4.2"
 gaussianprocessderivatives = "^0.1.5"
 numpy = ">=1.21.6"
+Bottleneck = "^1.3.5"
 opencv-python = "^4.7.0.72"
 pathos = "^0.2.8" # Lambda-friendly multithreading
 p-tqdm = "^1.3.3"
@@ -40,17 +43,16 @@ requests-toolbelt = "^0.9.1"
 scikit-image = ">=0.18.1"
 tqdm = "^4.62.3" # progress bars
 xmltodict = "^0.13.0" # read ome-tiff metadata
+zarr = "^2.14.0"
+GitPython = "^3.1.27"
+h5py = "2.10" # File I/O
 
 
 # Networking
 omero-py = { version = ">=5.6.2", optional = true } # contact omero server
-zeroc-ice = { version="3.6.5", optional = true } # networking interface, slow to build
-GitPython = "^3.1.27"
-Bottleneck = "^1.3.5"
-zarr = "^2.14.0"
 
-[tool.poetry.extras]
-omero = [ "omero-py" ]
+# Baby segmentation
+aliby-baby = {version = "^0.1.17", optional=true}
 
 # Postprocessing
 [tool.poetry.group.pp.dependencies]
@@ -58,10 +60,9 @@ leidenalg = "^0.8.8"
 more-itertools = "^8.12.0"
 pycatch22 = "^0.4.2"
 
+[tool.poetry.group.pp]
+optional = true
 
-[tool.poetry.group.baby.dependencies]
-aliby-baby = "^0.1.15"
-h5py = "2.10" # File I/O
 
 [tool.poetry.group.dev]
 optional = true
@@ -81,13 +82,7 @@ coverage = "^7.0.4"
 jupytext = "^1.14.4"
 grid-strategy = "^0.0.1"
 readchar = "^4.0.3"
-
-
-[tool.poetry.group.napari]
-optional = true
-
-[tool.poetry.group.napari.dependencies]
-napari = ">=0.4.16"
+ipdb = "^0.13.11"
 
 [tool.poetry.group.docs]
 optional = true
@@ -104,12 +99,21 @@ optional = true
 [tool.poetry.group.test.dependencies]
 pytest = "^6.2.5"
 
-[tool.poetry.group.torch]
+[tool.poetry.group.utils]
 optional = true
 
-[tool.poetry.group.torch.dependencies]
-torch = "^1.13.1"
-pytorch-lightning = "^1.9.3"
+# Dependency groups can only be used by a poetry installation, not pip
+[tool.poetry.group.utils.dependencies]
+napari = {version = ">=0.4.16", optional=true}
+Torch = {version = "^1.13.1", optional=true}
+pytorch-lightning = {version = "^1.9.3", optional=true}
+torchvision = {version = "^0.14.1", optional=true}
+trio = {version = "^0.22.0", optional=true}
+grid-strategy = {version = "^0.0.1", optional=true}
+
+[tool.poetry.extras]
+omero = ["omero-py"]
+baby = ["aliby-baby"]
 
 [tool.black]
 line-length = 79

src/agora/abc.py

@@ -3,7 +3,7 @@ import typing as t
 from abc import ABC, abstractmethod
 from collections.abc import Iterable
 from copy import copy
-from pathlib import Path, PosixPath
+from pathlib import Path
 from time import perf_counter
 from typing import Union
 
@@ -60,14 +60,14 @@ class ParametersABC(ABC):
         else:
             return iterable
 
-    def to_yaml(self, path: Union[PosixPath, str] = None):
+    def to_yaml(self, path: Union[Path, str] = None):
         """
         Returns a yaml stream of the attributes of the class instance.
         If path is provided, the yaml stream is saved there.
 
         Parameters
         ----------
-        path : Union[PosixPath, str]
+        path : Union[Path, str]
             Output path.
         """
         if path:
@@ -80,7 +80,7 @@ class ParametersABC(ABC):
         return cls(**d)
 
     @classmethod
-    def from_yaml(cls, source: Union[PosixPath, str]):
+    def from_yaml(cls, source: Union[Path, str]):
         """
         Returns instance from a yaml filename or stdin
         """
@@ -202,7 +202,7 @@ class ProcessABC(ABC):
     def run(self):
         pass
 
-    def _log(self, message: str, level: str = "warn"):
+    def _log(self, message: str, level: str = "warning"):
         # Log messages in the corresponding level
         logger = logging.getLogger("aliby")
         getattr(logger, level)(f"{self.__class__.__name__}: {message}")
@@ -211,7 +211,7 @@ class ProcessABC(ABC):
 def check_type_recursive(val1, val2):
     same_types = True
     if not isinstance(val1, type(val2)) and not all(
-        type(x) in (PosixPath, str) for x in (val1, val2)  # Ignore str->path
+        type(x) in (Path, str) for x in (val1, val2)  # Ignore str->path
     ):
         return False
     if not isinstance(val1, t.Iterable) and not isinstance(val2, t.Iterable):
@@ -249,5 +249,5 @@ class StepABC(ProcessABC):
         return self._run_tp(tp, **kwargs)
 
     def run(self):
-        # Replace run withn run_tp
+        # Replace run with run_tp
         raise Warning("Steps use run_tp instead of run")

src/agora/io/bridge.py

@@ -162,5 +162,8 @@ def image_creds_from_h5(fpath: str):
     attrs = attrs_from_h5(fpath)
     return (
         attrs["image_id"],
-        yaml.safe_load(attrs["parameters"])["general"]["server_info"],
+        {
+            k: yaml.safe_load(attrs["parameters"])["general"][k]
+            for k in ("username", "password", "host")
+        },
     )

src/agora/io/cells.py

 import logging
 import typing as t
 from itertools import groupby
-from pathlib import Path, PosixPath
+from pathlib import Path
 from functools import lru_cache, cached_property
 
 import h5py
@@ -26,13 +26,13 @@ class Cells:
     """
 
     def __init__(self, filename, path="cell_info"):
-        self.filename: t.Optional[t.Union[str, PosixPath]] = filename
+        self.filename: t.Optional[t.Union[str, Path]] = filename
         self.cinfo_path: t.Optional[str] = path
         self._edgemasks: t.Optional[str] = None
         self._tile_size: t.Optional[int] = None
 
     @classmethod
-    def from_source(cls, source: t.Union[PosixPath, str]):
+    def from_source(cls, source: t.Union[Path, str]):
         return cls(Path(source))
 
     def _log(self, message: str, level: str = "warn"):
@@ -278,7 +278,7 @@ class Cells:
     def ntimepoints(self) -> int:
         return self["timepoint"].max() + 1
 
-    @property
+    @cached_property
     def _cells_vs_tps(self):
         # Binary matrix showing the presence of all cells in all time points
         ncells_per_tile = [len(x) for x in self.labels]
@@ -286,13 +286,29 @@ class Cells:
             (sum(ncells_per_tile), self.ntimepoints), dtype=bool
         )
 
-        cumsum = np.roll(np.cumsum(ncells_per_tile), shift=1)
-        cumsum[0] = 0
         cells_vs_tps[
-            cumsum[self["trap"]] + self["cell_label"] - 1, self["timepoint"]
+            self._cell_cumsum[self["trap"]] + self["cell_label"] - 1,
+            self["timepoint"],
         ] = True
         return cells_vs_tps
 
+    @cached_property
+    def _cell_cumsum(self):
+        # Cumulative sum indicating the number of cells per tile
+        ncells_per_tile = [len(x) for x in self.labels]
+        cumsum = np.roll(np.cumsum(ncells_per_tile), shift=1)
+        cumsum[0] = 0
+        return cumsum
+
+    def _flat_index_to_tuple_location(self, idx: int) -> t.Tuple[int, int]:
+        # Convert a cell index to a tuple
+        # Note that it assumes tiles and cell labels are flattened, but
+        # it is agnostic to tps
+
+        tile_id = int(np.where(idx + 1 > self._cell_cumsum)[0][-1])
+        cell_label = idx - self._cell_cumsum[tile_id] + 1
+        return tile_id, cell_label
+
     @property
     def _tiles_vs_cells_vs_tps(self):
         ncells_mat = np.zeros(
@@ -304,39 +320,28 @@ class Cells:
         ] = True
         return ncells_mat
 
-    def cell_tp_where(self, min_consecutive_tps: int = 15):
+    def cell_tp_where(
+        self,
+        min_consecutive_tps: int = 15,
+        interval: None or t.Tuple[int, int] = None,
+    ):
         window = sliding_window_view(
             self._cells_vs_tps, min_consecutive_tps, axis=1
         )
 
         tp_min = window.sum(axis=-1) == min_consecutive_tps
-        return tp_min
-
-    def matrix_trap_tp_where(
-        self, min_ncells: int = 2, min_consecutive_tps: int = 5
-    ):
-        """
-        Return a matrix of shape (ntraps x ntps - min_consecutive_tps to
-        indicate traps and time-points where min_ncells are available for at least min_consecutive_tps
 
-        Parameters
-        ---------
-            min_ncells: int Minimum number of cells
-            min_consecutive_tps: int
-                Minimum number of time-points a
+        # Apply an interval filter to focucs on a slice
+        if interval is not None:
+            interval = tuple(np.array(interval))
+        else:
+            interval = (0, window.shape[1])
 
-        Returns
-        ---------
-            (ntraps x ( ntps-min_consecutive_tps )) 2D boolean numpy array where rows are trap ids and columns are timepoint windows.
-            If the value in a cell is true its corresponding trap and timepoint contains more than min_ncells for at least min_consecutive time-points.
-        """
+        low_boundary, high_boundary = interval
 
-        window = sliding_window_view(
-            self._tiles_vs_cells_vs_tps, min_consecutive_tps, axis=2
-        )
-        tp_min = window.sum(axis=-1) == min_consecutive_tps
-        ncells_tp_min = tp_min.sum(axis=1) >= min_ncells
-        return ncells_tp_min
+        tp_min[:, :low_boundary] = False
+        tp_min[:, high_boundary:] = False
+        return tp_min
 
     @lru_cache(20)
     def mothers_in_trap(self, trap_id: int):
@@ -572,10 +577,9 @@ class Cells:
     def _sample_tiles_tps(
         self,
         size=1,
-        # min_ncells: int = 2,
-        # max_ncells: int = 2,
-        min_consecutive_ntps: int = 10,
+        min_consecutive_ntps: int = 15,
         seed: int = 0,
+        interval=None,
     ) -> t.Tuple[np.ndarray, np.ndarray]:
         """
         Sample tiles that have a minimum number of cells and are occupied for at least a minimum number of consecutive timepoints.
@@ -590,10 +594,12 @@ class Cells:
             The minimum number of consecutive timepoints a cell must be present in a trap.
         seed: int, optional (default=0)
             Random seed value for reproducibility.
+        interval: None or Tuple(int,int), optional (default=None)
+            Random seed value for reproducibility.
 
         Returns
         -------
-        Tuple[np.ndarray, np.ndarray]
+        Tuple[np.ndarray, np.ndarray,np.ndarray]
             A tuple of 1D numpy arrays containing the indices of the sampled tiles and the corresponding timepoints.
         """
         # cell_availability_matrix = self.matrix_trap_tp_where(
@@ -603,30 +609,43 @@ class Cells:
         # # Find all valid tiles with min_ncells for at least min_tps
         # tile_ids, tps = np.where(cell_availability_matrix)
         cell_availability_matrix = self.cell_tp_where(
-            min_consecutive_tps=min_consecutive_ntps
+            min_consecutive_tps=min_consecutive_ntps,
+            interval=interval,
         )
 
         # Find all valid tiles with min_ncells for at least min_tps
-        index_id, _ = np.where(cell_availability_matrix)
+        index_id, tps = np.where(cell_availability_matrix)
+
+        if interval is None:  # Limit search
+            interval = (0, cell_availability_matrix.shape[1])
 
         np.random.seed(seed)
-        choices = np.random.choice(index_id, size=size)
-        return (
-            self["trap"][choices],
-            self["cell_label"][choices],
-            self["timepoint"][choices],
-        )
+        choices = np.random.randint(len(index_id), size=size)
+
+        linear_indices = np.zeros_like(self["cell_label"], dtype=bool)
+        for cell_index_flat, tp in zip(index_id[choices], tps[choices]):
+            tile_id, cell_label = self._flat_index_to_tuple_location(
+                cell_index_flat
+            )
+            linear_indices[
+                (
+                    (self["cell_label"] == cell_label)
+                    & (self["trap"] == tile_id)
+                    & (self["timepoint"] == tp)
+                )
+            ] = True
+
+        return linear_indices
 
     def _sample_masks(
         self,
-        size=1,
-        min_ncells: int = 2,
-        min_consecutive_ntps: int = 5,
+        size: int = 1,
+        min_consecutive_ntps: int = 15,
         interval: t.Union[None, t.Tuple[int, int]] = None,
-        seed=0,
+        seed: int = 0,
     ) -> t.Tuple[t.Tuple[t.List[int], t.List[int], t.List[int]], np.ndarray]:
         """
-        Sample a number of cells from different tiles each.
+        Sample a number of cells from within an interval.
 
         Parameters
         ----------
@@ -649,27 +668,53 @@ class Cells:
             The second tuple contains:
             - `masks`: A list of 2D numpy arrays representing the binary masks of the sampled cells at each timepoint.
         """
-        tile_ids, _, tps = self._sample_tiles_tps(
+        sampled_bitmask = self._sample_tiles_tps(
             size=size,
             min_consecutive_ntps=min_consecutive_ntps,
             seed=seed,
+            interval=interval,
         )
+
         #  Sort sampled tiles to use automatic cache when possible
-        order = np.argsort(tps)
-        tile_ids = tile_ids[order].tolist()
-        tps = tps[order].tolist()
+        tile_ids = self["trap"][sampled_bitmask]
+        cell_labels = self["cell_label"][sampled_bitmask]
+        tps = self["timepoint"][sampled_bitmask]
 
-        cell_ids = []
         masks = []
-        for tp, tile_id in zip(tps, tile_ids):
-            tile_masks = self.at_time(tp)[tile_id]
+        for tile_id, cell_label, tp in zip(tile_ids, cell_labels, tps):
+            local_idx = self.labels_at_time(tp)[tile_id].index(cell_label)
+            tile_mask = self.at_time(tp)[tile_id][local_idx]
+            masks.append(tile_mask)
+
+        return (tile_ids, cell_labels, tps), np.stack(masks)
 
-            np.random.seed(seed)
-            cell_id = np.random.randint(len(tile_masks))
+    def matrix_trap_tp_where(
+        self, min_ncells: int = 2, min_consecutive_tps: int = 5
+    ):
+        """
+        NOTE CURRENLTY UNUSED WITHIN ALIBY THE MOMENT. MAY BE USEFUL IN THE FUTURE.
+
+        Return a matrix of shape (ntraps x ntps - min_consecutive_tps) to
+        indicate traps and time-points where min_ncells are available for at least min_consecutive_tps
 
-            cell_ids.append(cell_id)
-            masks.append(tile_masks[cell_id])
-        return (tile_ids, cell_ids, tps), np.stack(masks)
+        Parameters
+        ---------
+            min_ncells: int Minimum number of cells
+            min_consecutive_tps: int
+                Minimum number of time-points a
+
+        Returns
+        ---------
+            (ntraps x ( ntps-min_consecutive_tps )) 2D boolean numpy array where rows are trap ids and columns are timepoint windows.
+            If the value in a cell is true its corresponding trap and timepoint contains more than min_ncells for at least min_consecutive time-points.
+        """
+
+        window = sliding_window_view(
+            self._tiles_vs_cells_vs_tps, min_consecutive_tps, axis=2
+        )
+        tp_min = window.sum(axis=-1) == min_consecutive_tps
+        ncells_tp_min = tp_min.sum(axis=1) >= min_ncells
+        return ncells_tp_min
 
 
 def stack_masks_in_tile(

src/agora/io/metadata.py

 """
-Anthology of interfaces for different parsers and lack of them.
+Anthology of interfaces fordispatch_metadata_parse different parsers and lack of them.
 
 ALIBY decides on using different metadata parsers based on two elements:
 
@@ -12,10 +12,11 @@ If there are no metadata files, ALIBY requires indicating indices for tiler, seg
 
 """
 import glob
+import logging
 import os
 import typing as t
 from datetime import datetime
-from pathlib import Path, PosixPath
+from pathlib import Path
 
 import pandas as pd
 from pytz import timezone
@@ -97,7 +98,9 @@ def find_file(root_dir, regex):
         )
         file = [sorted(file)[0]]
     if len(file) == 0:
-        print("Warning:Metadata: No valid swainlab .log found.")
+        logging.getLogger("aliby").log(
+            logging.WARNING, "Metadata: No valid swainlab .log found."
+        )
     else:
         return file[0]
     return None
@@ -173,7 +176,7 @@ def get_meta_from_legacy(parsed_metadata: dict):
     return result
 
 
-def parse_swainlab_metadata(filedir: t.Union[str, PosixPath]):
+def parse_swainlab_metadata(filedir: t.Union[str, Path]):
     """
     Dispatcher function that determines which parser to use based on the file ending.
 
@@ -202,7 +205,7 @@ def parse_swainlab_metadata(filedir: t.Union[str, PosixPath]):
     return minimal_meta
 
 
-def dispatch_metadata_parser(filepath: t.Union[str, PosixPath]):
+def dispatch_metadata_parser(filepath: t.Union[str, Path]):
     """
     Function to dispatch different metadata parsers that convert logfiles into a
     basic metadata dictionary. Currently only contains the swainlab log parsers.
@@ -219,7 +222,7 @@ def dispatch_metadata_parser(filepath: t.Union[str, PosixPath]):
     return parsed_meta
 
 
-def dir_to_meta(path: PosixPath, suffix="tiff"):
+def dir_to_meta(path: Path, suffix="tiff"):
     filenames = list(path.glob(f"*.{suffix}"))
 
     try:

src/agora/io/signal.py

@@ -2,7 +2,7 @@ import logging
 import typing as t
 from copy import copy
 from functools import cached_property, lru_cache
-from pathlib import PosixPath
+from pathlib import Path
 
 import bottleneck as bn
 import h5py
@@ -11,7 +11,7 @@ import pandas as pd
 
 from agora.io.bridge import BridgeH5
 from agora.io.decorators import _first_arg_str_to_df
-from agora.utils.association import validate_association
+from agora.utils.indexing import validate_association
 from agora.utils.kymograph import add_index_levels
 from agora.utils.merge import apply_merges
 
@@ -23,7 +23,7 @@ class Signal(BridgeH5):
     Signal assumes that the metadata and data are accessible to perform time-adjustments and apply previously recorded post-processes.
     """
 
-    def __init__(self, file: t.Union[str, PosixPath]):
+    def __init__(self, file: t.Union[str, Path]):
         """Define index_names for dataframes, candidate fluorescence channels, and composite statistics."""
         super().__init__(file, flag=None)
         self.index_names = (
@@ -47,20 +47,25 @@ class Signal(BridgeH5):
     def __getitem__(self, dsets: t.Union[str, t.Collection]):
         """Get and potentially pre-process data from h5 file and return as a dataframe."""
         if isinstance(dsets, str):  # no pre-processing
-            df = self.get_raw(dsets)
-            return self.add_name(df, dsets)
+            return self.get(dsets)
         elif isinstance(dsets, list):  # pre-processing
             is_bgd = [dset.endswith("imBackground") for dset in dsets]
             # Check we are not comparing tile-indexed and cell-indexed data
             assert sum(is_bgd) == 0 or sum(is_bgd) == len(
                 dsets
             ), "Tile data and cell data can't be mixed"
-            return [
-                self.add_name(self.apply_prepost(dset), dset) for dset in dsets
-            ]
+            return [self.get(dset) for dset in dsets]
         else:
             raise Exception(f"Invalid type {type(dsets)} to get datasets")
 
+    def get(self, dsets: t.Union[str, t.Collection], **kwargs):
+        """Get and potentially pre-process data from h5 file and return as a dataframe."""
+        if isinstance(dsets, str):  # no pre-processing
+            df = self.get_raw(dsets, **kwargs)
+            prepost_applied = self.apply_prepost(dsets, **kwargs)
+
+            return self.add_name(prepost_applied, dsets)
+
     @staticmethod
     def add_name(df, name):
         """Add column of identical strings to a dataframe."""
@@ -129,18 +134,24 @@ class Signal(BridgeH5):
         Returns an array with three columns: the tile id, the mother label, and the daughter label.
         """
         if lineage_location is None:
-            lineage_location = "postprocessing/lineage"
-            if merged:
-                lineage_location += "_merged"
+            lineage_location = "modifiers/lineage_merged"
         with h5py.File(self.filename, "r") as f:
+            # if lineage_location not in f:
+            #     lineage_location = lineage_location.split("_")[0]
+            if lineage_location not in f:
+                lineage_location = "postprocessing/lineage"
             tile_mo_da = f[lineage_location]
-            lineage = np.array(
-                (
-                    tile_mo_da["trap"],
-                    tile_mo_da["mother_label"],
-                    tile_mo_da["daughter_label"],
-                )
-            ).T
+
+            if isinstance(tile_mo_da, h5py.Dataset):
+                lineage = tile_mo_da[()]
+            else:
+                lineage = np.array(
+                    (
+                        tile_mo_da["trap"],
+                        tile_mo_da["mother_label"],
+                        tile_mo_da["daughter_label"],
+                    )
+                ).T
         return lineage
 
     @_first_arg_str_to_df
@@ -171,7 +182,7 @@ class Signal(BridgeH5):
 
         """
         if isinstance(merges, bool):
-            merges: np.ndarray = self.get_merges() if merges else np.array([])
+            merges: np.ndarray = self.load_merges() if merges else np.array([])
         if merges.any():
             merged = apply_merges(data, merges)
         else:
@@ -203,9 +214,8 @@ class Signal(BridgeH5):
                     merged = pd.DataFrame([], index=index)
         return merged
 
-    # Alan: do we need two similar properties - see below?
-    @property
-    def datasets(self):
+    @cached_property
+    def p_available(self):
         """Print data sets available in h5 file."""
         if not hasattr(self, "_available"):
             self._available = []
@@ -214,11 +224,6 @@ class Signal(BridgeH5):
         for sig in self._available:
             print(sig)
 
-    @cached_property
-    def p_available(self):
-        """Print data sets available in h5 file."""
-        self.datasets
-
     @cached_property
     def available(self):
         """Get data sets available in h5 file."""
@@ -298,7 +303,7 @@ class Signal(BridgeH5):
             self._log(f"Could not fetch dataset {dataset}: {e}", "error")
             raise e
 
-    def get_merges(self):
+    def load_merges(self):
         """Get merge events going up to the first level."""
         with h5py.File(self.filename, "r") as f:
             merges = f.get("modifiers/merges", np.array([]))
@@ -315,7 +320,9 @@ class Signal(BridgeH5):
         with h5py.File(self.filename, "r") as f:
             picks = set()
             if path in f:
-                picks = set(zip(*[f[path + name] for name in names]))
+                picks = set(
+                    zip(*[f[path + name] for name in names if name in f[path]])
+                )
             return picks
 
     def dataset_to_df(self, f: h5py.File, path: str) -> pd.DataFrame:

src/agora/io/writer.py

@@ -172,7 +172,6 @@ class DynamicWriter:
                             # append or create new dataset
                             self._append(value, key, hgroup)
                     except Exception as e:
-                        print(key, value)
                         self._log(
                             f"{key}:{value} could not be written: {e}", "error"
                         )
@@ -230,7 +229,6 @@ class LinearBabyWriter(DynamicWriter):
     Assumes the edgemasks are of form ((None, tile_size, tile_size), bool).
     """
 
-    # TODO make this YAML: Alan: why?
     compression = "gzip"
     _default_tile_size = 117
     datatypes = {
@@ -319,11 +317,7 @@ class StateWriter(DynamicWriter):
     @staticmethod
     def format_values_tpback(states: list, val_name: str):
         """Unpacks a dict of state data into tp_back, trap, value."""
-        # initialise as empty lists
-        # Alan: is this initialisation necessary?
-        tp_back, trap, value = [
-            [[] for _ in states[0][val_name]] for _ in range(3)
-        ]
+
         # store results as a list of tuples
         lbl_tuples = [
             (tp_back, trap, cell_label)
@@ -334,6 +328,11 @@ class StateWriter(DynamicWriter):
         # unpack list of tuples to define variables
         if len(lbl_tuples):
             tp_back, trap, value = zip(*lbl_tuples)
+        else:
+            # set as empty lists
+            tp_back, trap, value = [
+                [[] for _ in states[0][val_name]] for _ in range(3)
+            ]
         return tp_back, trap, value
 
     @staticmethod
@@ -409,9 +408,9 @@ class StateWriter(DynamicWriter):
 
 #################### Extraction version ###############################
 class Writer(BridgeH5):
-    """Class to transform data into compatible structures."""
-
-    # Alan: when is this used?
+    """
+    Class to transform data into compatible structures.
+    Used by Extractor and Postprocessor within the pipeline."""
 
     def __init__(self, filename, flag=None, compression="gzip"):
         """
@@ -473,7 +472,7 @@ class Writer(BridgeH5):
             self.write_pd(f, path, data, compression=self.compression)
         # data is a multi-index dataframe
         elif isinstance(data, pd.MultiIndex):
-            # Alan: should we still not compress here?
+            # TODO: benchmark I/O speed when using compression
             self.write_index(f, path, data)  # , compression=self.compression)
         # data is a dictionary of dataframes
         elif isinstance(data, Dict) and np.all(
@@ -550,7 +549,7 @@ class Writer(BridgeH5):
                 compression=kwargs.get("compression", None),
             )
             dset = f[values_path]
-            dset[()] = df.values
+            dset[()] = df.values.astype("float16")
 
             # create dateset and write indices
             if not len(df):  # Only write more if not empty
@@ -567,21 +566,18 @@ class Writer(BridgeH5):
                 )
                 dset = f[indices_path]
                 dset[()] = df.index.get_level_values(level=name).tolist()
-            # create dataset and write columns
-            if (
-                df.columns.dtype == int
-                or df.columns.dtype == np.dtype("uint")
-                or df.columns.name == "timepoint"
-            ):
+
+                # create dataset and write time points as columns
                 tp_path = path + "/timepoint"
-                f.create_dataset(
-                    name=tp_path,
-                    shape=(df.shape[1],),
-                    maxshape=(max_tps,),
-                    dtype="uint16",
-                )
-                tps = list(range(df.shape[1]))
-                f[tp_path][tps] = tps
+                if tp_path not in f:
+                    f.create_dataset(
+                        name=tp_path,
+                        shape=(df.shape[1],),
+                        maxshape=(max_tps,),
+                        dtype="uint16",
+                    )
+                    tps = list(range(df.shape[1]))
+                    f[tp_path][tps] = tps
             else:
                 f[path].attrs["columns"] = df.columns.tolist()
         else:

src/agora/utils/association.py

 #!/usr/bin/env jupyter
-"""
-Utilities based on association are used to efficiently acquire indices of tracklets with some kind of relationship.
-This can be:
-    - Cells that are to be merged
-    - Cells that have a linear relationship
-"""
-
-import numpy as np
-import typing as t
-
-
-def validate_association(
-    association: np.ndarray,
-    indices: np.ndarray,
-    match_column: t.Optional[int] = None,
-) -> t.Tuple[np.ndarray, np.ndarray]:
-
-    """Select rows from the first array that are present in both.
-        We use casting for fast multiindexing, generalising for lineage dynamics
-
-
-        Parameters
-        ----------
-        association : np.ndarray
-            2-D array where columns are (trap, mother, daughter) or 3-D array where
-            dimensions are (X, (trap,mother), (trap,daughter))
-        indices : np.ndarray
-            2-D array where each column is a different level.
-        match_column: int
-            int indicating a specific column is required to match (i.e.
-            0-1 for target-source when trying to merge tracklets or mother-bud for lineage)
-            must be present in indices. If it is false one match suffices for the resultant indices
-            vector to be True.
-
-        Returns
-        -------
-        np.ndarray
-            1-D boolean array indicating valid merge events.
-        np.ndarray
-            1-D boolean array indicating indices with an association relationship.
-
-        Examples
-        --------
-
-        >>> import numpy as np
-        >>> from agora.utils.association import validate_association
-        >>> merges = np.array(range(12)).reshape(3,2,2)
-        >>> indices = np.array(range(6)).reshape(3,2)
-
-        >>> print(merges, indices)
-        >>> print(merges); print(indices)
-        [[[ 0  1]
-          [ 2  3]]
-
-         [[ 4  5]
-          [ 6  7]]
-
-         [[ 8  9]
-          [10 11]]]
-
-        [[0 1]
-         [2 3]
-         [4 5]]
-
-        >>> valid_associations, valid_indices  = validate_association(merges, indices)
-        >>> print(valid_associations, valid_indices)
-    [ True False False] [ True  True False]
-
-    """
-    if association.ndim < 3:
-        # Reshape into 3-D array for broadcasting if neded
-        association = np.stack(
-            (association[:, [0, 1]], association[:, [0, 2]]), axis=1
-        )
-
-    # Compare existing association with available indices
-    # Swap trap and label axes for the association array to correctly cast
-    valid_ndassociation = association[..., None] == indices.T[None, ...]
-
-    # Broadcasting is confusing (but efficient):
-    # First we check the dimension across trap and cell id, to ensure both match
-    valid_cell_ids = valid_ndassociation.all(axis=2)
-
-    if match_column is None:
-        # Then we check the merge tuples to check which cases have both target and source
-        valid_association = valid_cell_ids.any(axis=2).all(axis=1)
-
-        # Finally we check the dimension that crosses all indices, to ensure the pair
-        # is present in a valid merge event.
-        valid_indices = (
-            valid_ndassociation[valid_association].all(axis=2).any(axis=(0, 1))
-        )
-    else:  # We fetch specific indices if we aim for the ones with one present
-        valid_indices = valid_cell_ids[:, match_column].any(axis=0)
-        # Valid association then becomes a boolean array, true means that there is a
-        # match (match_column) between that cell and the index
-        valid_association = (
-            valid_cell_ids[:, match_column] & valid_indices
-        ).any(axis=1)
-
-    return valid_association, valid_indices

src/agora/utils/cast.py

+#!/usr/bin/env jupyter
+
+"""
+Convert some types to others
+"""
+
+
+def _str_to_int(x: str or None):
+    """
+    Cast string as int if possible. If Nonetype return None.
+    """
+    if x is not None:
+        try:
+            return int(x)
+        except:
+            return x

src/agora/utils/indexing.py

+#!/usr/bin/env jupyter
+"""
+Utilities based on association are used to efficiently acquire indices of tracklets with some kind of relationship.
+This can be:
+    - Cells that are to be merged
+    - Cells that have a linear relationship
+"""
+
+import numpy as np
+import typing as t
+
+
+def validate_association(
+    association: np.ndarray,
+    indices: np.ndarray,
+    match_column: t.Optional[int] = None,
+) -> t.Tuple[np.ndarray, np.ndarray]:
+
+    """Select rows from the first array that are present in both.
+        We use casting for fast multiindexing, generalising for lineage dynamics
+
+
+        Parameters
+        ----------
+        association : np.ndarray
+            2-D array where columns are (trap, mother, daughter) or 3-D array where
+            dimensions are (X,trap,2), containing tuples ((trap,mother), (trap,daughter))
+            across the 3rd dimension.
+        indices : np.ndarray
+            2-D array where each column is a different level. This should not include mother_label.
+        match_column: int
+            int indicating a specific column is required to match (i.e.
+            0-1 for target-source when trying to merge tracklets or mother-bud for lineage)
+            must be present in indices. If it is false one match suffices for the resultant indices
+            vector to be True.
+
+        Returns
+        -------
+        np.ndarray
+            1-D boolean array indicating valid merge events.
+        np.ndarray
+            1-D boolean array indicating indices with an association relationship.
+
+        Examples
+        --------
+
+        >>> import numpy as np
+        >>> from agora.utils.indexing import validate_association
+        >>> merges = np.array(range(12)).reshape(3,2,2)
+        >>> indices = np.array(range(6)).reshape(3,2)
+
+        >>> print(merges, indices)
+        >>> print(merges); print(indices)
+        [[[ 0  1]
+          [ 2  3]]
+
+         [[ 4  5]
+          [ 6  7]]
+
+         [[ 8  9]
+          [10 11]]]
+
+        [[0 1]
+         [2 3]
+         [4 5]]
+
+        >>> valid_associations, valid_indices  = validate_association(merges, indices)
+        >>> print(valid_associations, valid_indices)
+    [ True False False] [ True  True False]
+
+    """
+    if association.ndim == 2:
+        # Reshape into 3-D array for broadcasting if neded
+        # association = np.stack(
+        #     (association[:, [0, 1]], association[:, [0, 2]]), axis=1
+        # )
+        association = _assoc_indices_to_3d(association)
+
+    # Compare existing association with available indices
+    # Swap trap and label axes for the association array to correctly cast
+    valid_ndassociation = association[..., None] == indices.T[None, ...]
+
+    # Broadcasting is confusing (but efficient):
+    # First we check the dimension across trap and cell id, to ensure both match
+    valid_cell_ids = valid_ndassociation.all(axis=2)
+
+    if match_column is None:
+        # Then we check the merge tuples to check which cases have both target and source
+        valid_association = valid_cell_ids.any(axis=2).all(axis=1)
+
+        # Finally we check the dimension that crosses all indices, to ensure the pair
+        # is present in a valid merge event.
+        valid_indices = (
+            valid_ndassociation[valid_association].all(axis=2).any(axis=(0, 1))
+        )
+    else:  # We fetch specific indices if we aim for the ones with one present
+        valid_indices = valid_cell_ids[:, match_column].any(axis=0)
+        # Valid association then becomes a boolean array, true means that there is a
+        # match (match_column) between that cell and the index
+        valid_association = (
+            valid_cell_ids[:, match_column] & valid_indices
+        ).any(axis=1)
+
+    return valid_association, valid_indices
+
+
+def _assoc_indices_to_3d(ndarray: np.ndarray):
+    """
+    Convert the last column to a new row while repeating all previous indices.
+
+    This is useful when converting a signal multiindex before comparing association.
+
+    Assumes the input array has shape (N,3)
+    """
+    result = ndarray
+    if len(ndarray) and ndarray.ndim > 1:
+        if ndarray.shape[1] == 3:  # Faster indexing for single positions
+            result = np.transpose(
+                np.hstack((ndarray[:, [0]], ndarray)).reshape(-1, 2, 2),
+                axes=[0, 2, 1],
+            )
+        else:  # 20% slower but more general indexing
+            columns = np.arange(ndarray.shape[1])
+
+            result = np.stack(
+                (
+                    ndarray[:, np.delete(columns, -1)],
+                    ndarray[:, np.delete(columns, -2)],
+                ),
+                axis=1,
+            )
+    return result
+
+
+def _3d_index_to_2d(array: np.ndarray):
+    """
+    Opposite to _assoc_indices_to_3d.
+    """
+    result = array
+    if len(array):
+        result = np.concatenate(
+            (array[:, 0, :], array[:, 1, 1, np.newaxis]), axis=1
+        )
+    return result
+
+
+def compare_indices(x: np.ndarray, y: np.ndarray) -> np.ndarray:
+    """
+    Fetch two 2-D indices and return a binary 2-D matrix
+    where a True value links two cells where all cells are the same
+    """
+    return (x[..., None] == y.T[None, ...]).all(axis=1)