.flake8

+[flake8]
+ignore = E203, E266, E501, W503, F403, F401
+max-line-length = 79
+select = B,C,E,F,W,T4,B9
+exclude =
+    # No need to traverse our git directory
+    .git,
+    # There's no value in checking cache directories
+    __pycache__,
+    # Ignore virtual environment contents
+    .venv
+    # The conf file is mostly autogenerated, ignore it
+    docs/source/conf.py,
+    # The old directory contains Flake8 2.0
+    old,
+    # This contains our built documentation
+    build,
+    # This contains builds of flake8 that we don't want to check
+    dist,
+    # Any data produced inside the folder during development
+    data/,
+    # Temporarily ignore tests
+    tests/
+max-complexity = 18

.gitignore

@@ -110,7 +110,6 @@ venv.bak/
 omero_py/omeroweb/
 omero_py/pipeline/
 **.ipynb
-data/
 notebooks/
 *.pdf
 *.h5

.gitlab-ci.yml

-image: python:3.7
+image: python:3.8
 
 cache:
     key: "project-${CI_JOB_NAME}"
     paths:
-      - .cache/pip
       - .venv
     key:
       files:
         - poetry.lock
 
+variables:
+  TRIGGER_PYPI_NAME: ""
+
 stages:
-  - test
-  - check
-  - release
+  - tests
+  - checks
+  # - release
 
 before_script:
   - test -e $HOME/.poetry/bin/ || curl -sSL https://install.python-poetry.org | python3 -
@@ -20,52 +22,68 @@ before_script:
   - poetry --version
   - poetry config virtualenvs.in-project true
   - pip install --upgrade pip
-  - git remote rm origin && git remote add origin https://${ACCESS_TOKEN_NAME}:${ACCESS_TOKEN}@${CI_SERVER_HOST}/${CI_PROJECT_PATH}.git
-  - git config pull.rebase false
-  - git pull origin HEAD:master
-  - if [ ${var+TRIGGER_PYPI_NAME} ]; then echo "Pipeline triggered by ${TRIGGER_PYPI_NAME}"; poetry update ${TRIGGER_PYPI_NAME}; fi
-  - rm -f poetry.lock
-  - poetry install -vv
-  # - poetry export -f requirements.txt --output requirements.txt
-  # - sed '/^zeroc-ice/,/^[z.*]/{/z.*/!d}' requirements.txt | grep -v zeroc-ice > nozice.txt
-  # - pip install -r nozice.txt --no-deps && pip install pytest black mypy
+  # - git remote rm origin && git remote add origin https://${ACCESS_TOKEN_NAME}:${ACCESS_TOKEN}@${CI_SERVER_HOST}/${CI_PROJECT_PATH}.git
+  # - git config pull.rebase false
+  # - git pull origin HEAD:master
+  # - rm -rf ~/.cache/pypoetry
+  # - if [ ${var+TRIGGER_PYPI_NAME} ]; then echo "Pipeline triggered by ${TRIGGER_PYPI_NAME}"; poetry add ${TRIGGER_PYPI_NAME}@latest; fi
+  # - export WITHOUT="docs,network";
+  - export ARGS="--with test,dev";
+  - if [[ "$CI_STAGE_NAME" == "tests" ]]; then echo "Installing system dependencies for ${CI_STAGE_NAME}"; apt update && apt install -y ffmpeg libsm6 libxext6; fi
+  - if [[ "$CI_JOB_NAME" == "Static Type" ]]; then echo "Activating development group"; export ARGS="${ARGS},dev"; fi
+  - if [[ "$CI_JOB_NAME" == "Network Tools Tests" ]]; then echo "Setting flag to compile zeroc-ice"; export ARGS="${ARGS} --all-extras"; fi
+  - poetry install -vv $ARGS
 
-Unit test:
-  stage: test
-  # allow_failure: true
+Local Tests:
+  stage: tests
   script:
-    - apt update && apt install ffmpeg libsm6 libxext6  -y
-    - poetry run pytest ./tests/
+    # - poetry install -vv
+    - poetry run coverage run -m --branch pytest ./tests --ignore ./tests/aliby/network --ignore ./tests/aliby/pipeline
+    - poetry run coverage report -m
+    - poetry run coverage xml
+  coverage: '/(?i)total.*? (100(?:\.0+)?\%|[1-9]?\d(?:\.\d+)?\%)$/'
+  artifacts:
+    reports:
+      coverage_report:
+        coverage_format: cobertura
+        path: coverage.xml
 
-Python Code Lint:
-  stage: check
+Network Tools Tests:
+  stage: tests
   script:
-    - poetry run black .
+    - poetry run pytest ./tests/aliby/network
+    - DIRNAME="test_datasets"
+    - curl https://zenodo.org/record/7513194/files/test_datasets.tar.gz\?download\=1 -o "test_datasets.tar.gz"
+    - mkdir -p $DIRNAME
+    - tar xvf test_datasets.tar.gz -C $DIRNAME
+    - poetry run pytest -s tests/aliby/pipeline --file $DIRNAME/560_2022_11_30_pypipeline_unit_test_reconstituted_00
 
 Static Type:
-  stage: check
+  stage: checks
   allow_failure: true
   script:
+    - poetry run black .
+    - poetry run isort .
     - poetry run mypy . --exclude 'setup\.py$'
     # We can remove the flag once this is resolved https://github.com/pypa/setuptools/issues/2345
 
-Bump_release:
-  stage: release
-  script:
-    - poetry version ${BUMP_RULE}
-    - git config --global user.email ${GITLAB_USER_EMAIL}
-    - git config --global user.name ${GITLAB_USER_NAME}
-    - git add pyproject.toml poetry.lock
-    - git commit -m "Bump version"
-    # - git pull  && git push -o ci.skip origin HEAD:master && poetry publish --build --username ${PYPI_USER} --password ${PYPI_PASSWORD}
-    # - git pull  && git push -o ci.skip origin HEAD:master
-    - echo "TRIGGER_PYPI_NAME=$(cat pyproject.toml | grep '^name =' | head -n 1 | cut -f3 -d' ' | tr -d \")" >> build.env
-  artifacts:
-    reports:
-      dotenv: build.env
-  only:
-    - master
-  needs:
-    job: Unit test
+# TODO add more tests before activating auto-release
+# Bump_release:
+#   stage: release
+#   script:
+#     - git config --global user.email ${GITLAB_USER_EMAIL}
+#     - git config --global user.name ${GITLAB_USER_NAME}
+#     # - git pull origin HEAD:MASTER && poetry version ${BUMP_RULE} && git add poetry.lock add pyproject.toml poetry.lock && git commit -m "Bump version" && git push -o ci.skip origin HEAD:master && poetry publish --build --username ${PYPI_USER} --password ${PYPI_PASSWORD}
+#     - echo "TRIGGER_PYPI_NAME=$(cat pyproject.toml | grep '^name =' | head -n 1 | cut -f3 -d' ' | tr -d \")" >> build.env
+#     - echo "Exporting TRIGGER_PYPI_NAME as ${TRIGGER_PYPI_NAME}"
+#   only:
+#     - master
+#   except:
+#     changes:
+#       - tests/
+#       - .*
+
+#   needs:
+#     job: Unit test
 
 ## Custom stages ##

.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,7 +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)
-

.pre-commit-config.yaml

+exclude: '^examples/'
+repos:
+  - repo: https://github.com/pre-commit/pre-commit-hooks
+    rev: v4.3.0
+    hooks:
+      # - id: check-yaml
+      - id: end-of-file-fixer
+      - id: trailing-whitespace
+        args: [--markdown-linebreak-ext=md]
+      - id: check-added-large-files
+      - id: check-docstring-first
+      - id: name-tests-test
+        args: [--pytest-test-first]
+  - repo: https://github.com/psf/black
+    rev: 22.6.0
+    hooks:
+      - id: black
+  - repo: https://github.com/PyCQA/isort
+    rev: 5.10.1
+    hooks:
+      - id: isort
+  - repo: https://github.com/PyCQA/flake8
+    rev: 4.0.1
+    hooks:
+      - id: flake8
+        exclude: ^tests/ # Temporarily exclude tests

.readthedocs.yaml

+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-20.04
+  tools:
+    python: "3.8"
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+   configuration: docs/source/conf.py
+
+# If using Sphinx, optionally build your docs in additional formats such as PDF
+# formats:
+#    - pdf
+
+# Optionally declare the Python requirements required to build your docs
+python:
+   install:
+   - requirements: docs/requirements.txt

CONTRIBUTING.md

+# Contributing
+
+We focus our work on python 3.8 due to the current neural network being developed on tensorflow 1. In the near future we will migrate the network to pytorch to support more recent versions of all packages.
+
+## Issues
+All issues are managed within the gitlab [ repository ](https://gitlab.com/aliby/aliby/-/issues), if you don't have an account on the University of Edinburgh's gitlab instance and would like to submit issues please get in touch with [Prof. Peter Swain](mailto:peter.swain@ed.ac.uk ).
+
+## Data aggregation
+
+ALIBY has been tested by a few research groups, but we welcome new data sources for the models and pipeline to be as general as possible. Please get in touch with [ us ](mailto:peter.swain@ed.ac.uk ) if you are interested in testing it on your data.

INSTALL.md

-# Installation 
-
-Tested on: Mac OSX Mojave
-
-## Requirements
-We strongly recommend installing within a python environment as there are many dependencies that you may not want polluting your regular python environment. 
-Make sure you are using python 3.
-
-An environment can be created with using the conda package manager:
-
-    $ conda create --name <env>
-    $ conda activate <env>
-
-Which you can deactivate with: 
-
-    $ conda deactivate
-
-Or using virtualenv: 
-
-    $ python -m virtualenv /path/to/venv/
-    $ source /path/to/venv/bin/activate
-    
-This will download all of your packages under `/path/to/venv` and then activate it. 
-Deactivate using 
-
-    $ deactivate
-
-You will also need to make sure you have a recent version of pip.
-In your local environment, run: 
-
-    $ pip install --upgrade pip
-
-## Pipeline installation
-
-Once you have created your local environment, run:
-
-    $ cd pipeline-core
-    $ pip install -e ./
-
-You will be asked to put in your GitLab credentials for a couple of the packages installed as dependencies.
-The `-e` option will install in 'editable' mode: all of the changes that are made to the code in the repository will immediately be reflected in the installation. 
-This is very useful to keep up with any changes and branching that we make. 
-That way, if you want to keep up with the most recent updates, just run: 
-
-    $ git pull
-
-If you would rather be more in charge of which updates you install and which you don't, remove `-e` from your installation command. 
-In this case (or if you run into a dependency error) in order to update your installed version you will have to run: 
-
-    $ cd pipeline-core
-    $ git pull
-    $ pip install ./
-

README.md

-# Pipeline core
+# ALIBY (Analyser of Live-cell Imaging for Budding Yeast)
 
-The core classes and methods for the python microfluidics, microscopy, and 
-analysis pipeline.
-
-### Installation
-See [INSTALL.md](./INSTALL.md) for installation instructions.
+[![docs](https://readthedocs.org/projects/aliby/badge/?version=master)](https://aliby.readthedocs.io/en/latest)
+[![PyPI version](https://badge.fury.io/py/aliby.svg)](https://badge.fury.io/py/aliby)
+[![pipeline](https://gitlab.com/aliby/aliby/badges/master/pipeline.svg?key_text=master)](https://gitlab.com/aliby/aliby/-/pipelines)
+[![dev pipeline](https://gitlab.com/aliby/aliby/badges/dev/pipeline.svg?key_text=dev)](https://gitlab.com/aliby/aliby/-/commits/dev)
+[![coverage](https://gitlab.com/aliby/aliby/badges/dev/coverage.svg)](https://gitlab.com/aliby/aliby/-/commits/dev)
 
+End-to-end processing of cell microscopy time-lapses. ALIBY automates segmentation, tracking, lineage predictions, post-processing and report production. It leverages the existing Python ecosystem and open-source scientific software available to produce seamless and standardised pipelines.
 
 ## Quickstart Documentation
-### Setting up a server
-For testing and development, the easiest way to set up an OMERO server is by
-using Docker images. 
-[The software carpentry](https://software-carpentry.org/) and the [Open
- Microscopy Environment](https://www.openmicroscopy.org), have provided
-[instructions](https://ome.github.io/training-docker/) to do this.
-
-The `docker-compose.yml` file can be used to create an OMERO server with an
-accompanying PostgreSQL database, and an OMERO web server.
-It is described in detail 
-[here](https://ome.github.io/training-docker/12-dockercompose/).
-
-Our version of the `docker-compose.yml` has been adapted from the above to
-use version 5.6 of OMERO.
-
-To start these containers (in background):
-```shell script
-cd pipeline-core
-docker-compose up -d
-```
-Omit the `-d` to run in foreground.
+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.
 
-To stop them, in the same directory, run:
-```shell script
-docker-compose stop
-```
+To analyse local data
+ ```bash
+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
+ ```
+
+And to run Omero servers, the basic arguments are shown:
+ ```bash
+ aliby-run --expt_id XXX --host SERVER.ADDRESS --user USER --password PASSWORD 
+ ```
+
+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)
 
-### Raw data access
+## Using specific components
 
+### Access raw data
+
+ALIBY's tooling can also be used as an interface to OMERO servers, for example, to fetch a brightfield channel.
  ```python
 from aliby.io.omero import Dataset, Image
 
@@ -56,129 +60,53 @@ with Image(list(image_ids.values())[0], **server_info) as image:
     imgs = dimg[tps, image.metadata["channels"].index("Brightfield"), 2, ...].compute()
     # tps timepoints, Brightfield channel, z=2, all x,y
 ```
- 
+
 ### Tiling the raw data
 
-A `Tiler` object performs trap registration. It is built in different ways, the easiest one is using an image and a the default parameters set.
+A `Tiler` object performs trap registration. It may be built in different ways but the simplest one is using an image and a the default parameters set.
 
 ```python
 from aliby.tile.tiler import Tiler, TilerParameters
 with Image(list(image_ids.values())[0], **server_info) as image:
     tiler = Tiler.from_image(image, TilerParameters.default())
+    tiler.run_tp(0)
 ```
 
 The initialisation should take a few seconds, as it needs to align the images
-in time. 
+in time.
 
 It fetches the metadata from the Image object, and uses the TilerParameters values (all Processes in aliby depend on an associated Parameters class, which is in essence a dictionary turned into a class.)
 
-#### Get a timelapse for a given trap
-TODO: Update this
+#### Get a timelapse for a given tile (remote connection)
 ```python
-channels = [0] #Get only the first channel, this is also the default
-z = [0, 1, 2, 3, 4] #Get all z-positions
-trap_id = 0
-tile_size = 117
-
-# Get a timelapse of the trap
-# The default trap size is 96 by 96
-# The trap is in the center of the image, except for edge cases
-# The output has shape (C, T, X, Y, Z), so in this example: (1, T, 96, 96, 5)
-timelapse = seg_expt.get_trap_timelapse(trap_id, tile_size=tile_size, 
-                                        channels=channels, z=z)
-```
+fpath = "h5/location"
 
-This can take several seconds at the moment.
-For a speed-up: take fewer z-positions if you can.
+tile_id = 9
+trange = range(0, 10)
+ncols = 8
 
-If you're not sure what indices to use:
-```python
-seg_expt.channels # Get a list of channels
-channel = 'Brightfield'
-ch_id = seg_expt.get_channel_index(channel)
+riv = remoteImageViewer(fpath)
+trap_tps = [riv.tiler.get_tiles_timepoint(tile_id, t) for t in trange] 
 
-n_traps = seg_expt.n_traps # Get the number of traps 
-```
-
-#### Get the traps for a given time point
-Alternatively, if you want to get all the traps at a given timepoint:
+# You can also access labelled traps
+m_ts = riv.get_labelled_trap(tile_id=0, tps=[0])
 
-```python
-timepoint = 0
-seg_expt.get_traps_timepoints(timepoint, tile_size=96, channels=None, 
-                                z=[0,1,2,3,4])
-```
-
-## Reading MATLAB files
-*Disclaimer: this is very much still in development so it may not always
- work for you case. If you run into any problems please let me know, or even
- better start an Issue on the project describing your problem.*
-
-At the moment the best/only way to read matlab files is through a `matObject`:
-
-```python
-from aliby.io.matlab import matObject
-cTimelapse = matObject('/path/to/cTimelapse.mat')
-```
-
-You can see an overview of what's in the object:
-```python
-cTimelapse.describe()
+# And plot them directly
+riv.plot_labelled_trap(trap_id=0, channels=[0, 1, 2, 3], trange=range(10))
 ```
 
-The `matObject` has some dictionary-like features although it is *not* a
-dictionary (yet). You can access different parts of the object using keys
-, though, and can use the `keys()` function to do so. This will usually
-work at the first few levels, but if it doesn't you may have run into an
-object that's actually a list or a numpy array.
-
-```python
-cTimelapse.keys()
-```
+Depending on the network speed can take several seconds at the moment.
+For a speed-up: take fewer z-positions if you can.
 
-This should return an iterable of the upper level keys. For example, a
-timelapse object will usually have a `timelapseTrapsOmero` key which you
-can look deeper into in the same manner. Once you've found what you want
-you can usually access it as you would a nested dictionary, for instance:  
+#### Get the tiles for a given time point
+Alternatively, if you want to get all the traps at a given timepoint:
 
 ```python
-cTimelapse['timelapseTrapsOmero']['cTimepoint']['trapLocations']
-```
-
-For more information about using MATLAB files in python objects, please see
- [this page](docs/matlab.md).
-
-## Development guidelines
-In order to separate the python2, python3, and "currently working" versions 
-(\#socialdistancing) of the pipeline, please use the branches:
-* python2.7: for any development on the 2 version
-* python3.6-dev: for any added features for the python3 version
-* master: very sparingly and only for changes that need to be made in both
- versions as I will be merging changes from master into the development
- branches frequently
-    * Ideally for adding features into any branch, espeically master, create
-     a new branch first, then create a pull request (from within Gitlab) before 
-     merging it back so we can check each others' code. This is just to make
-     sure that we can always use the code that is in the master branch without
-     any issues.
- 
-Branching cheat-sheet:
-```git
-git branch my_branch # Create a new branch called branch_name from master
-git branch my_branch another_branch #Branch from another_branch, not master
-git checkout -b my_branch # Create my_branch and switch to it
-
-# Merge changes from master into your branch
-git pull #get any remote changes in master
-git checkout my_branch
-git merge master
-
-# Merge changes from your branch into another branch
-git checkout another_branch
-git merge my_branch #check the doc for --no-ff option, you might want to use it
+timepoint = (4,6)
+tiler.get_tiles_timepoint(timepoint, channels=None,
+                                z=[0,1,2,3,4])
 ```
 
-## TODO
-### Tests
-* test full pipeline with OMERO experiment (no download.)
 
+### Contributing
+See [CONTRIBUTING](https://aliby.readthedocs.io/en/latest/INSTALL.html) on how to help out or get involved.

aliby/experiment.py

-"""Core classes for the pipeline"""
-import atexit
-import itertools
-import os
-import abc
-import glob
-import json
-import warnings
-from getpass import getpass
-from pathlib import Path
-import re
-import logging
-from typing import Union
-
-import h5py
-from tqdm import tqdm
-import pandas as pd
-
-import omero
-from omero.gateway import BlitzGateway
-from logfile_parser import Parser
-
-from agora.io.utils import accumulate
-from aliby.timelapse import TimelapseOMERO, TimelapseLocal
-
-from agora.io.writer import Writer
-
-logger = logging.getLogger(__name__)
-
-########################### Dask objects ###################################
-##################### ENVIRONMENT INITIALISATION ################
-import omero
-from omero.gateway import BlitzGateway, PixelsWrapper
-from omero.model import enums as omero_enums
-import numpy as np
-
-# Set up the pixels so that we can reuse them across sessions (?)
-PIXEL_TYPES = {
-    omero_enums.PixelsTypeint8: np.int8,
-    omero_enums.PixelsTypeuint8: np.uint8,
-    omero_enums.PixelsTypeint16: np.int16,
-    omero_enums.PixelsTypeuint16: np.uint16,
-    omero_enums.PixelsTypeint32: np.int32,
-    omero_enums.PixelsTypeuint32: np.uint32,
-    omero_enums.PixelsTypefloat: np.float32,
-    omero_enums.PixelsTypedouble: np.float64,
-}
-
-
-class NonCachedPixelsWrapper(PixelsWrapper):
-    """Extend gateway.PixelWrapper to override _prepareRawPixelsStore."""
-
-    def _prepareRawPixelsStore(self):
-        """
-        Creates RawPixelsStore and sets the id etc
-        This overrides the superclass behaviour to make sure that
-        we don't re-use RawPixelStore in multiple processes since
-        the Store may be closed in 1 process while still needed elsewhere.
-        This is needed when napari requests may planes simultaneously,
-        e.g. when switching to 3D view.
-        """
-        ps = self._conn.c.sf.createRawPixelsStore()
-        ps.setPixelsId(self._obj.id.val, True, self._conn.SERVICE_OPTS)
-        return ps
-
-
-omero.gateway.PixelsWrapper = NonCachedPixelsWrapper
-# Update the BlitzGateway to use our NonCachedPixelsWrapper
-omero.gateway.refreshWrappers()
-
-
-######################  DATA ACCESS ###################
-import dask.array as da
-from dask import delayed
-
-
-def get_data_lazy(image) -> da.Array:
-    """Get 5D dask array, with delayed reading from OMERO image."""
-    nt, nc, nz, ny, nx = [getattr(image, f"getSize{x}")() for x in "TCZYX"]
-    pixels = image.getPrimaryPixels()
-    dtype = PIXEL_TYPES.get(pixels.getPixelsType().value, None)
-    get_plane = delayed(lambda idx: pixels.getPlane(*idx))
-
-    def get_lazy_plane(zct):
-        return da.from_delayed(get_plane(zct), shape=(ny, nx), dtype=dtype)
-
-    # 5D stack: TCZXY
-    t_stacks = []
-    for t in range(nt):
-        c_stacks = []
-        for c in range(nc):
-            z_stack = []
-            for z in range(nz):
-                z_stack.append(get_lazy_plane((z, c, t)))
-            c_stacks.append(da.stack(z_stack))
-        t_stacks.append(da.stack(c_stacks))
-    return da.stack(t_stacks)
-
-
-# Metadata writer
-from aliby.io.metadata_parser import parse_logfiles
-
-
-class MetaData:
-    """Small metadata Process that loads log."""
-
-    def __init__(self, log_dir, store):
-        self.log_dir = log_dir
-        self.store = store
-
-    def load_logs(self):
-        parsed_flattened = parse_logfiles(self.log_dir)
-        return parsed_flattened
-
-    def run(self):
-        metadata_writer = Writer(self.store)
-        metadata_dict = self.load_logs()
-        metadata_writer.write(path="/", meta=metadata_dict, overwrite=False)
-
-
-########################### Old Objects ####################################
-
-
-class Experiment(abc.ABC):
-    """
-    Abstract base class for experiments.
-    Gives all the functions that need to be implemented in both the local
-    version and the Omero version of the Experiment class.
-
-    As this is an abstract class, experiments can not be directly instantiated
-    through the usual `__init__` function, but must be instantiated from a
-    source.
-    >>> expt = Experiment.from_source(root_directory)
-    Data from the current timelapse can be obtained from the experiment using
-    colon and comma separated slicing.
-    The order of data is C, T, X, Y, Z
-    C, T and Z can have any slice
-    X and Y will only consider the beginning and end as we want the images
-    to be continuous
-    >>> bf_1 = expt[0, 0, :, :, :] # First channel, first timepoint, all x,y,z
-    """
-
-    __metaclass__ = abc.ABCMeta
-
-    # metadata_parser = AcqMetadataParser()
-
-    def __init__(self):
-        self.exptID = ""
-        self._current_position = None
-        self.position_to_process = 0
-
-    def __getitem__(self, item):
-        return self.current_position[item]
-
-    @property
-    def shape(self):
-        return self.current_position.shape
-
-    @staticmethod
-    def from_source(*args, **kwargs):
-        """
-        Factory method to construct an instance of an Experiment subclass (
-        either ExperimentOMERO or ExperimentLocal).
-
-        :param source: Where the data is stored (OMERO server or directory
-        name)
-        :param kwargs: If OMERO server, `user` and `password` keyword
-        arguments are required. If the data is stored locally keyword
-        arguments are ignored.
-        """
-        if len(args) > 1:
-            logger.debug("ExperimentOMERO: {}".format(args, kwargs))
-            return ExperimentOMERO(*args, **kwargs)
-        else:
-            logger.debug("ExperimentLocal: {}".format(args, kwargs))
-            return ExperimentLocal(*args, **kwargs)
-
-    @property
-    @abc.abstractmethod
-    def positions(self):
-        """Returns list of available position names"""
-        return
-
-    @abc.abstractmethod
-    def get_position(self, position):
-        return
-
-    @property
-    def current_position(self):
-        return self._current_position
-
-    @property
-    def channels(self):
-        return self._current_position.channels
-
-    @current_position.setter
-    def current_position(self, position):
-        self._current_position = self.get_position(position)
-
-    def get_hypercube(self, x, y, z_positions, channels, timepoints):
-        return self.current_position.get_hypercube(
-            x, y, z_positions, channels, timepoints
-        )
-
-
-# Todo: cache images like in ExperimentLocal
-class ExperimentOMERO(Experiment):
-    """
-    Experiment class to organise different timelapses.
-    Connected to a Dataset object which handles database I/O.
-    """
-
-    def __init__(self, omero_id, host, port=4064, **kwargs):
-        super(ExperimentOMERO, self).__init__()
-        self.exptID = omero_id
-        # Get annotations
-        self.use_annotations = kwargs.get("use_annotations", True)
-        self._files = None
-        self._tags = None
-
-        # Create a connection
-        self.connection = BlitzGateway(
-            kwargs.get("username") or input("Username: "),
-            kwargs.get("password") or getpass("Password: "),
-            host=host,
-            port=port,
-        )
-        connected = self.connection.connect()
-        try:
-            assert connected is True, "Could not connect to server."
-        except AssertionError as e:
-            self.connection.close()
-            raise (e)
-        try:  # Run everything that could cause the initialisation to fail
-            self.dataset = self.connection.getObject("Dataset", self.exptID)
-            self.name = self.dataset.getName()
-            # Create positions objects
-            self._positions = {
-                img.getName(): img.getId()
-                for img in sorted(
-                    self.dataset.listChildren(), key=lambda x: x.getName()
-                )
-            }
-            # Set up local cache
-            self.root_dir = Path(kwargs.get("save_dir", "./")) / self.name
-            if not self.root_dir.exists():
-                self.root_dir.mkdir(parents=True)
-            self.compression = kwargs.get("compression", None)
-            self.image_cache = h5py.File(self.root_dir / "images.h5", "a")
-
-            # Set up the current position as the first in the list
-            self._current_position = self.get_position(self.positions[0])
-            self.running_tp = 0
-        except Exception as e:
-            # Close the connection!
-            print("Error in initialisation, closing connection.")
-            self.connection.close()
-            print(self.connection.isConnected())
-            raise e
-        atexit.register(self.close)  # Close everything if program ends
-
-    def close(self):
-        print("Clean-up on exit.")
-        self.image_cache.close()
-        self.connection.close()
-
-    @property
-    def files(self):
-        if self._files is None:
-            self._files = {
-                x.getFileName(): x
-                for x in self.dataset.listAnnotations()
-                if isinstance(x, omero.gateway.FileAnnotationWrapper)
-            }
-        return self._files
-
-    @property
-    def tags(self):
-        if self._tags is None:
-            self._tags = {
-                x.getName(): x
-                for x in self.dataset.listAnnotations()
-                if isinstance(x, omero.gateway.TagAnnotationWrapper)
-            }
-        return self._tags
-
-    @property
-    def positions(self):
-        return list(self._positions.keys())
-
-    def _get_position_annotation(self, position):
-        # Get file annotations filtered by position name and ordered by
-        # creation date
-        r = re.compile(position)
-        wrappers = sorted(
-            [self.files[key] for key in filter(r.match, self.files)],
-            key=lambda x: x.creationEventDate(),
-            reverse=True,
-        )
-        # Choose newest file
-        if len(wrappers) < 1:
-            return None
-        else:
-            # Choose the newest annotation and cache it
-            annotation = wrappers[0]
-            filepath = self.root_dir / annotation.getFileName().replace("/", "_")
-            if not filepath.exists():
-                with open(str(filepath), "wb") as fd:
-                    for chunk in annotation.getFileInChunks():
-                        fd.write(chunk)
-            return filepath
-
-    def get_position(self, position):
-        """Get a Timelapse object for a given position by name"""
-        # assert position in self.positions, "Position not available."
-        img = self.connection.getObject("Image", self._positions[position])
-        if self.use_annotations:
-            annotation = self._get_position_annotation(position)
-        else:
-            annotation = None
-        return TimelapseOMERO(img, annotation, self.image_cache)
-
-    def cache_locally(
-        self,
-        root_dir="./",
-        positions=None,
-        channels=None,
-        timepoints=None,
-        z_positions=None,
-    ):
-        """
-        Save the experiment locally.
-
-        :param root_dir: The directory in which the experiment will be
-        saved. The experiment will be a subdirectory of "root_directory"
-        and will be named by its id.
-        """
-        logger.warning("Saving experiment {}; may take some time.".format(self.name))
-
-        if positions is None:
-            positions = self.positions
-        if channels is None:
-            channels = self.current_position.channels
-        if timepoints is None:
-            timepoints = range(self.current_position.size_t)
-        if z_positions is None:
-            z_positions = range(self.current_position.size_z)
-
-        save_dir = Path(root_dir) / self.name
-        if not save_dir.exists():
-            save_dir.mkdir()
-        # Save the images
-        for pos_name in tqdm(positions):
-            pos = self.get_position(pos_name)
-            pos_dir = save_dir / pos_name
-            if not pos_dir.exists():
-                pos_dir.mkdir()
-            self.cache_set(pos, range(pos.size_t))
-
-        self.cache_logs(save_dir)
-        # Save the file annotations
-        cache_config = dict(
-            positions=positions,
-            channels=channels,
-            timepoints=timepoints,
-            z_positions=z_positions,
-        )
-        with open(str(save_dir / "cache.config"), "w") as fd:
-            json.dump(cache_config, fd)
-        logger.info("Downloaded experiment {}".format(self.exptID))
-
-    def cache_logs(self, **kwargs):
-        # Save the file annotations
-        tags = dict()  # and the tag annotations
-        for annotation in self.dataset.listAnnotations():
-            if isinstance(annotation, omero.gateway.FileAnnotationWrapper):
-                filepath = self.root_dir / annotation.getFileName().replace("/", "_")
-                if str(filepath).endswith("txt") and not filepath.exists():
-                    # Save only the text files
-                    with open(str(filepath), "wb") as fd:
-                        for chunk in annotation.getFileInChunks():
-                            fd.write(chunk)
-            if isinstance(annotation, omero.gateway.TagAnnotationWrapper):
-                key = annotation.getDescription()
-                if key == "":
-                    key = "misc. tags"
-                if key in tags:
-                    if not isinstance(tags[key], list):
-                        tags[key] = [tags[key]]
-                    tags[key].append(annotation.getValue())
-                else:
-                    tags[key] = annotation.getValue()
-        with open(str(self.root_dir / "omero_tags.json"), "w") as fd:
-            json.dump(tags, fd)
-        return
-
-    def run(self, keys: Union[list, int], store, **kwargs):
-        if self.running_tp == 0:
-            self.cache_logs(**kwargs)
-            self.running_tp = 1  # Todo rename based on annotations
-        run_tps = dict()
-        for pos, tps in accumulate(keys):
-            position = self.get_position(pos)
-            run_tps[pos] = position.run(tps, store, save_dir=self.root_dir)
-        # Update the keys to match what was actually run
-        keys = [(pos, tp) for pos in run_tps for tp in run_tps[pos]]
-        return keys
-
-
-class ExperimentLocal(Experiment):
-    def __init__(self, root_dir, finished=True):
-        super(ExperimentLocal, self).__init__()
-        self.root_dir = Path(root_dir)
-        self.exptID = self.root_dir.name
-        self._pos_mapper = dict()
-        # Fixme: Made the assumption that the Acq file gets saved before the
-        #  experiment is run and that the information in that file is
-        #  trustworthy.
-        acq_file = self._find_acq_file()
-        acq_parser = Parser("multiDGUI_acq_format")
-        with open(acq_file, "r") as fd:
-            metadata = acq_parser.parse(fd)
-        self.metadata = metadata
-        self.metadata["finished"] = finished
-        self.files = [f for f in self.root_dir.iterdir() if f.is_file()]
-        self.image_cache = h5py.File(self.root_dir / "images.h5", "a")
-        if self.finished:
-            cache = self._find_cache()
-            # log = self._find_log() # Todo: add log metadata
-            if cache is not None:
-                with open(cache, "r") as fd:
-                    cache_config = json.load(fd)
-                self.metadata.update(**cache_config)
-        self._current_position = self.get_position(self.positions[0])
-
-    def _find_file(self, regex):
-        file = glob.glob(os.path.join(str(self.root_dir), regex))
-        if len(file) != 1:
-            return None
-        else:
-            return file[0]
-
-    def _find_acq_file(self):
-        file = self._find_file("*[Aa]cq.txt")
-        if file is None:
-            raise ValueError(
-                "Cannot load this experiment. There are either "
-                "too many or too few acq files."
-            )
-        return file
-
-    def _find_cache(self):
-        return self._find_file("cache.config")
-
-    @property
-    def finished(self):
-        return self.metadata["finished"]
-
-    @property
-    def running(self):
-        return not self.metadata["finished"]
-
-    @property
-    def positions(self):
-        return self.metadata["positions"]["posname"]
-
-    def _get_position_annotation(self, position):
-        r = re.compile(position)
-        files = list(filter(lambda x: r.match(x.stem), self.files))
-        if len(files) == 0:
-            return None
-        files = sorted(files, key=lambda x: x.lstat().st_ctime, reverse=True)
-        # Get the newest and return as string
-        return files[0]
-
-    def get_position(self, position):
-        if position not in self._pos_mapper:
-            annotation = self._get_position_annotation(position)
-            self._pos_mapper[position] = TimelapseLocal(
-                position,
-                self.root_dir,
-                finished=self.finished,
-                annotation=annotation,
-                cache=self.image_cache,
-            )
-        return self._pos_mapper[position]
-
-    def run(self, keys, store, **kwargs):
-        """
-
-        :param keys: List of (position, time point) tuples to process.
-        :return:
-        """
-        run_tps = dict()
-        for pos, tps in accumulate(keys):
-            run_tps[pos] = self.get_position(pos).run(tps, store)
-        # Update the keys to match what was actually run
-        keys = [(pos, tp) for pos in run_tps for tp in run_tps[pos]]
-        return keys

aliby/extract.py

-"""
-A module to extract data from a processed experiment.
-"""
-import h5py
-import numpy as np
-from tqdm import tqdm
-
-from core.io.matlab import matObject
-from growth_rate.estimate_gr import estimate_gr
-
-
-class Extracted:
-    # TODO write the filtering functions.
-    def __init__(self):
-        self.volume = None
-        self._keep = None
-
-    def filter(self, filename=None, **kwargs):
-        """
-        1. Filter out small non-growing tracks. This means:
-            a. the cell size never reaches beyond a certain size-threshold
-            volume_thresh or
-            b. the cell's volume doesn't increase by at least a minimum
-            amount over its lifetime
-        2. Join daughter tracks that are contiguous and within a volume
-           threshold of each other
-        3. Discard tracks that are shorter than a threshold number of
-           timepoints
-
-        This function is used to fix tracking/bud errors in post-processing.
-        The parameters define the thresholds used to determine which cells are
-        discarded.
-        FIXME Ideally we get to a point where this is no longer needed.
-        :return:
-        """
-        #self.join_tracks()
-        filter_out = self.filter_size(**kwargs)
-        filter_out += self.filter_lifespan(**kwargs)
-        # TODO save data or just filtering parameters?
-        #self.to_hdf(filename)
-        self.keep = ~filter_out
-
-    def filter_size(self, volume_thresh=7, growth_thresh=10, **kwargs):
-        """Filter out small and non-growing cells.
-        :param volume_thresh: Size threshold for small cells
-        :param growth_thresh: Size difference threshold for non-growing cells
-        """
-        filter_out = np.where(np.max(self.volume, axis=1) < volume_thresh,
-                              True, False)
-        growth = [v[v > 0] for v in self.volume]
-        growth = np.array([v[-1] - v[0] if len(v) > 0 else 0 for v in growth])
-        filter_out += np.where(growth < growth_thresh, True, False)
-        return filter_out
-
-    def filter_lifespan(self, min_time=5, **kwargs):
-        """Remove daughter cells that have a small life span.
-
-        :param min_time: The minimum life span, under which cells are removed.
-        """
-        # TODO What if there are nan values?
-        filter_out = np.where(np.count_nonzero(self.volume, axis=1) <
-                              min_time, True, False)
-        return filter_out
-
-    def join_tracks(self, threshold=7):
-        """ Join contiguous tracks that are within a certain volume
-        threshold of each other.
-
-        :param threshold: Maximum volume difference to join contiguous tracks.
-        :return:
-        """
-        # For all pairs of cells
-        #
-        pass
-
-
-class ExtractedHDF(Extracted):
-    # TODO pull all the data out of the HFile and filter!
-    def __init__(self, file):
-        # We consider the data to be read-only
-        self.hfile = h5py.File(file, 'r')
-
-
-class ExtractedMat(Extracted):
-    """ Pulls the extracted data out of the MATLAB cTimelapse file.
-
-    This is mostly a convenience function in order to run the
-    gaussian-processes growth-rate estimation
-    """
-    def __init__(self, file, debug=False):
-        ct = matObject(file)
-        self.debug = debug
-        # Pre-computed data
-        # TODO what if there is no timelapseTrapsOmero?
-        self.metadata = ct['timelapseTrapsOmero']['metadata']
-        self.extracted_data = ct['timelapseTrapsOmero']['extractedData']
-        self.channels = ct['timelapseTrapsOmero']['extractionParameters'][
-            'functionParameters']['channels'].tolist()
-        self.time_settings = ct['timelapseTrapsOmero']['metadata']['acq'][
-            'times']
-        # Get filtering information
-        n_cells = self.extracted_data['cellNum'][0].shape
-        self.keep = np.full(n_cells, True)
-        # Not yet computed data
-        self._growth_rate = None
-        self._daughter_index = None
-
-
-    def get_channel_index(self, channel):
-        """Get index of channel based on name. This only considers
-        fluorescence channels."""
-        return self.channels.index(channel)
-
-    @property
-    def trap_num(self):
-        return self.extracted_data['trapNum'][0][self.keep]
-
-    @property
-    def cell_num(self):
-        return self.extracted_data['cellNum'][0][self.keep]
-
-    def identity(self, cell_idx):
-        """Get the (position), trap, and cell label given a cell's global
-        index."""
-        # Todo include position when using full strain
-        trap = self.trap_num[cell_idx]
-        cell = self.cell_num[cell_idx]
-        return trap, cell
-
-    def global_index(self, trap_id, cell_label):
-        """Get the global index of a cell given it's trap/cellNum
-        combination."""
-        candidates = np.where(np.logical_and(
-                            (self.trap_num == trap_id), # +1?
-                            (self.cell_num == cell_label)
-                        ))[0]
-        # TODO raise error if number of candidates != 1
-        if len(candidates) == 1:
-            return candidates[0]
-        elif len(candidates) == 0:
-            return -1
-        else:
-            raise(IndexError("No such cell/trap combination"))
-
-    @property
-    def daughter_label(self):
-        """Returns the cell label of the daughters of each cell over the
-        timelapse.
-
-        0 corresponds to no daughter. This *not* the index of the daughter
-        cell within the data. To get this, use daughter_index.
-        """
-        return self.extracted_data['daughterLabel'][0][self.keep]
-
-    def _single_daughter_idx(self, mother_idx, daughter_labels):
-        trap_id, _ = self.identity(mother_idx)
-        daughter_index = [self.global_index(trap_id, cell_label) for
-                          cell_label
-                          in daughter_labels]
-        return daughter_index
-
-    @property
-    def daughter_index(self):
-        """Returns the global index of the daughters of each cell.
-
-        This is different from the daughter label because it corresponds to
-        the index of the daughter when counting all of the cells. This can
-        be used to index within the data arrays.
-        """
-        if self._daughter_index is None:
-            daughter_index = [self._single_daughter_idx(i, daughter_labels)
-                          for i, daughter_labels in enumerate(
-                                  self.daughter_label)]
-            self._daughter_index = np.array(daughter_index)
-        return self._daughter_index
-
-    @property
-    def births(self):
-        return np.array(self.extracted_data['births'][0].todense())[self.keep]
-
-    @property
-    def volume(self):
-        """Get the volume of all of the cells"""
-        return np.array(self.extracted_data['volume'][0].todense())[self.keep]
-
-    def _gr_estimation(self):
-        dt = self.time_settings['interval'] / 360  # s to h conversion
-        results = []
-        for v in tqdm(self.volume):
-            results.append(estimate_gr(v, dt))
-        merged = {k: np.stack([x[k] for x in results]) for k in results[0]}
-        self._gr_results = merged
-        return
-
-    @property
-    def growth_rate(self):
-        """Get the growth rate for all cells.
-
-        Note that this uses the gaussian processes method of estimating
-        growth rate by default. If there is no growth rate in the given file
-        (usually the case for MATLAB), it needs to run estimation first.
-        This can take a while.
-        """
-        # TODO cache the results of growth rate estimation.
-        if self._gr_results is None:
-            dt = self.time_settings['interval'] / 360  # s to h conversion
-            self._growth_rate = [estimate_gr(v, dt) for v in self.volume]
-        return self._gr_results['growth_rate']
-
-    def _fluo_attribute(self, channel, attribute):
-        channel_id = self.get_channel_index(channel)
-        res = np.array(self.extracted_data[attribute][channel_id].todense())
-        return res[self.keep]
-
-    def protein_localisation(self, channel, method='nucEstConv'):
-        """Returns protein localisation data for a given channel.
-
-        Uses the 'nucEstConv' by default. Alternatives are 'smallPeakConv',
-        'max5px', 'max2p5pc'
-        """
-        return self._fluo_attribute(channel, method)
-
-    def background_fluo(self, channel):
-        return self._fluo_attribute(channel, 'imBackground')
-
-    def mean(self, channel):
-        return self._fluo_attribute(channel, 'mean')
-
-    def median(self, channel):
-        return self._fluo_attribute(channel, 'median')
-
-    def filter(self, filename=None):
-        """Filters and saves results to and HDF5 file.
-
-        This is necessary because we cannot write to the MATLAB file,
-        so the results of the filter cannot be saved in the object.
-        """
-        super().filter(filename=filename)
-        self._growth_rate = None  # reset growth rate so it is recomputed
-
-    def to_hdf(self, filename):
-        """Store the current results, including any filtering done, to a file.
-
-        TODO Should we save filtered results or just re-do?
-        :param filename:
-        :return:
-        """
-        store = h5py.File(filename, 'w')
-        try:
-            # Store (some of the) metadata
-            for meta in ['experiment', 'username', 'microscope',
-                              'comments', 'project', 'date', 'posname',
-                              'exptid']:
-                store.attrs[meta] = self.metadata[meta]
-            # TODO store timing information?
-            store.attrs['time_interval'] = self.time_settings['interval']
-            store.attrs['timepoints'] = self.time_settings['ntimepoints']
-            store.attrs['total_duration'] = self.time_settings['totalduration']
-            # Store volume, births, daughterLabel, trapNum, cellNum
-            for key in ['volume', 'births', 'daughter_label', 'trap_num',
-                        'cell_num']:
-                store[key] = getattr(self, key)
-            # Store growth rate results
-            if self._gr_results:
-                grp = store.create_group('gaussian_process')
-                for key, val in self._gr_results.items():
-                    grp[key] = val
-            for channel in self.channels:
-                # Create a group for each channel
-                grp = store.create_group(channel)
-                # Store protein_localisation, background fluorescence, mean, median
-                # for each channel
-                grp['protein_localisation'] = self.protein_localisation(channel)
-                grp['background_fluo'] = self.background_fluo(channel)
-                grp['mean'] = self.mean(channel)
-                grp['median'] = self.median(channel)
-        finally:
-            store.close()
-

aliby/grouper.py

-#!/usr/bin/env python3
-
-from abc import ABC, abstractmethod, abstractproperty
-from pathlib import Path
-from pathos.multiprocessing import Pool
-
-import h5py
-import numpy as np
-import pandas as pd
-
-from agora.io.signal import Signal
-
-
-class Grouper(ABC):
-    """
-    Base grouper class
-    """
-
-    files = []
-
-    def __init__(self, dir):
-        self.files = list(Path(dir).glob("*.h5"))
-        self.load_signals()
-
-    def load_signals(self):
-        self.signals = {f.name[:-3]: Signal(f) for f in self.files}
-
-    @property
-    def fsignal(self):
-        return list(self.signals.values())[0]
-
-    @property
-    def siglist(self):
-        return self.fsignal.datasets
-
-    @abstractproperty
-    def group_names():
-        pass
-
-    def concat_signal(self, path, reduce_cols=None, axis=0, pool=8):
-        group_names = self.group_names
-        sitems = self.signals.items()
-        if pool:
-            with Pool(pool) as p:
-                signals = p.map(
-                    lambda x: concat_signal_ind(path, group_names, x[0], x[1]),
-                    sitems,
-                )
-        else:
-            signals = [
-                concat_signal_ind(path, group_names, name, signal)
-                for name, signal in sitems
-            ]
-
-        signals = [s for s in signals if s is not None]
-        sorted = pd.concat(signals, axis=axis).sort_index()
-        if reduce_cols:
-            sorted = sorted.apply(np.nanmean, axis=1)
-            spath = path.split("/")
-            sorted.name = "_".join([spath[1], spath[-1]])
-
-        return sorted
-
-    @property
-    def ntraps(self):
-        for pos, s in self.signals.items():
-            with h5py.File(s.filename, "r") as f:
-                print(pos, f["/trap_info/trap_locations"].shape[0])
-
-    def traplocs(self):
-        d = {}
-        for pos, s in self.signals.items():
-            with h5py.File(s.filename, "r") as f:
-                d[pos] = f["/trap_info/trap_locations"][()]
-        return d
-
-
-class MetaGrouper(Grouper):
-    """Group positions using metadata's 'group' number"""
-
-    pass
-
-
-class NameGrouper(Grouper):
-    """
-    Group a set of positions using a subsection of the name
-    """
-
-    def __init__(self, dir, by=None):
-        super().__init__(dir=dir)
-
-        if by is None:
-            by = (0, -4)
-        self.by = by
-
-    @property
-    def group_names(self):
-        if not hasattr(self, "_group_names"):
-            self._group_names = {}
-            for name in self.signals.keys():
-                self._group_names[name] = name[self.by[0] : self.by[1]]
-
-        return self._group_names
-
-    def aggregate_multisignals(self, paths=None, **kwargs):
-
-        aggregated = pd.concat(
-            [
-                self.concat_signal(path, reduce_cols=np.nanmean, **kwargs)
-                for path in paths
-            ],
-            axis=1,
-        )
-        # ph = pd.Series(
-        #     [
-        #         self.ph_from_group(x[list(aggregated.index.names).index("group")])
-        #         for x in aggregated.index
-        #     ],
-        #     index=aggregated.index,
-        #     name="media_pH",
-        # )
-        # self.aggregated = pd.concat((aggregated, ph), axis=1)
-
-        return aggregated
-
-
-class phGrouper(NameGrouper):
-    """
-    Grouper for pH calibration experiments where all surveyed media pH values
-    are within a single experiment.
-    """
-
-    def __init__(self, dir, by=(3, 7)):
-        super().__init__(dir=dir, by=by)
-
-    def get_ph(self):
-        self.ph = {gn: self.ph_from_group(gn) for gn in self.group_names}
-
-    @staticmethod
-    def ph_from_group(group_name):
-        if group_name.startswith("ph_"):
-            group_name = group_name[3:]
-
-        return float(group_name.replace("_", "."))
-
-    def aggregate_multisignals(self, paths):
-
-        aggregated = pd.concat(
-            [self.concat_signal(path, reduce_cols=np.nanmean) for path in paths], axis=1
-        )
-        ph = pd.Series(
-            [
-                self.ph_from_group(x[list(aggregated.index.names).index("group")])
-                for x in aggregated.index
-            ],
-            index=aggregated.index,
-            name="media_pH",
-        )
-        aggregated = pd.concat((aggregated, ph), axis=1)
-
-        return aggregated
-
-
-def concat_signal_ind(path, group_names, group, signal):
-    print("Looking at ", group)
-    # try:
-    combined = signal[path]
-    combined["position"] = group
-    combined["group"] = group_names[group]
-    combined.set_index(["group", "position"], inplace=True, append=True)
-    combined.index = combined.index.swaplevel(-2, 0).swaplevel(-1, 1)
-
-    return combined
-    # except:
-    #     return None

aliby/io/matlab.py

-"""Read and convert MATLAB files from Swain Lab platform.
-
-TODO: Information that I need from lab members esp J and A
-    * Lots of examples to try
-    * Any ideas on what these Map objects are?
-
-TODO: Update Swain Lab wiki
-
-All credit to Matt Bauman for
-the reverse engineering at https://nbviewer.jupyter.org/gist/mbauman/9121961
-"""
-
-import re
-import struct
-import sys
-from collections import Iterable
-from io import BytesIO
-
-import h5py
-import numpy as np
-import pandas as pd
-import scipy
-from numpy.compat import asstr
-
-# TODO only use this if scipy>=1.6 or so
-from scipy.io import matlab
-from scipy.io.matlab.mio5 import MatFile5Reader
-from scipy.io.matlab.mio5_params import mat_struct
-
-from aliby.io.utils import read_int, read_string, read_delim
-
-
-def read_minimat_vars(rdr):
-    rdr.initialize_read()
-    mdict = {"__globals__": []}
-    i = 0
-    while not rdr.end_of_stream():
-        hdr, next_position = rdr.read_var_header()
-        name = asstr(hdr.name)
-        if name == "":
-            name = "var_%d" % i
-            i += 1
-        res = rdr.read_var_array(hdr, process=False)
-        rdr.mat_stream.seek(next_position)
-        mdict[name] = res
-        if hdr.is_global:
-            mdict["__globals__"].append(name)
-    return mdict
-
-
-def read_workspace_vars(fname):
-    fp = open(fname, "rb")
-    rdr = MatFile5Reader(fp, struct_as_record=True, squeeze_me=True)
-    vars = rdr.get_variables()
-    fws = vars["__function_workspace__"]
-    ws_bs = BytesIO(fws.tostring())
-    ws_bs.seek(2)
-    rdr.mat_stream = ws_bs
-    # Guess byte order.
-    mi = rdr.mat_stream.read(2)
-    rdr.byte_order = mi == b"IM" and "<" or ">"
-    rdr.mat_stream.read(4)  # presumably byte padding
-    mdict = read_minimat_vars(rdr)
-    fp.close()
-    return mdict
-
-
-class matObject:
-    """A python read-out of MATLAB objects
-    The objects pulled out of the
-    """
-
-    def __init__(self, filepath):
-        self.filepath = filepath  # For record
-        self.classname = None
-        self.object_name = None
-        self.buffer = None
-        self.version = None
-        self.names = None
-        self.segments = None
-        self.heap = None
-        self.attrs = dict()
-        self._init_buffer()
-        self._init_heap()
-        self._read_header()
-        self.parse_file()
-
-    def __getitem__(self, item):
-        return self.attrs[item]
-
-    def keys(self):
-        """Returns the names of the available properties"""
-        return self.attrs.keys()
-
-    def get(self, item, default=None):
-        return self.attrs.get(item, default)
-
-    def _init_buffer(self):
-        fp = open(self.filepath, "rb")
-        rdr = MatFile5Reader(fp, struct_as_record=True, squeeze_me=True)
-        vars = rdr.get_variables()
-        self.classname = vars["None"]["s2"][0].decode("utf-8")
-        self.object_name = vars["None"]["s0"][0].decode("utf-8")
-        fws = vars["__function_workspace__"]
-        self.buffer = BytesIO(fws.tostring())
-        fp.close()
-
-    def _init_heap(self):
-        super_data = read_workspace_vars(self.filepath)
-        elem = super_data["var_0"][0, 0]
-        if isinstance(elem, mat_struct):
-            self.heap = elem.MCOS[0]["arr"]
-        else:
-            self.heap = elem["MCOS"][0]["arr"]
-
-    def _read_header(self):
-        self.buffer.seek(248)  # the start of the header
-        version = read_int(self.buffer)
-        n_str = read_int(self.buffer)
-
-        offsets = read_int(self.buffer, n=6)
-
-        # check that the next two are zeros
-        reserved = read_int(self.buffer, n=2)
-        assert all(
-            [x == 0 for x in reserved]
-        ), "Non-zero reserved header fields: {}".format(reserved)
-        # check that we are at the right place
-        assert self.buffer.tell() == 288, "String elemnts begin at 288"
-        hdrs = []
-        for i in range(n_str):
-            hdrs.append(read_string(self.buffer))
-        self.names = hdrs
-        self.version = version
-        # The offsets are actually STARTING FROM 248 as well
-        self.segments = [x + 248 for x in offsets]  # list(offsets)
-        return
-
-    def parse_file(self):
-        # Get class attributes from segment 1
-        self.buffer.seek(self.segments[0])
-        classes = self._parse_class_attributes(self.segments[1])
-        # Get first set of properties from segment 2
-        self.buffer.seek(self.segments[1])
-        props1 = self._parse_properties(self.segments[2])
-        # Get the property description from segment 3
-        self.buffer.seek(self.segments[2])
-        object_info = self._parse_prop_description(classes, self.segments[3])
-        # Get more properties from segment 4
-        self.buffer.seek(self.segments[3])
-        props2 = self._parse_properties(self.segments[4])
-        # Check that the last segment is empty
-        self.buffer.seek(self.segments[4])
-        seg5_length = (self.segments[5] - self.segments[4]) // 8
-        read_delim(self.buffer, seg5_length)
-        props = (props1, props2)
-        self._to_attrs(object_info, props)
-
-    def _to_attrs(self, object_info, props):
-        """Re-organise the various classes and subclasses into a nested
-        dictionary.
-        :return:
-        """
-        for pkg_clss, indices, idx in object_info:
-            pkg, clss = pkg_clss
-            idx = max(indices)
-            which = indices.index(idx)
-            obj = flatten_obj(props[which][idx])
-            subdict = self.attrs
-            if pkg != "":
-                subdict = self.attrs.setdefault(pkg, {})
-            if clss in subdict:
-                if isinstance(subdict[clss], list):
-                    subdict[clss].append(obj)
-                else:
-                    subdict[clss] = [subdict[clss]]
-                    subdict[clss].append(obj)
-            else:
-                subdict[clss] = obj
-
-    def describe(self):
-        describe(self.attrs)
-
-    def _parse_class_attributes(self, section_end):
-        """Read the Class attributes = the first segment"""
-        read_delim(self.buffer, 4)
-        classes = []
-        while self.buffer.tell() < section_end:
-            package_index = read_int(self.buffer) - 1
-            package = self.names[package_index] if package_index > 0 else ""
-            name_idx = read_int(self.buffer) - 1
-            name = self.names[name_idx] if name_idx > 0 else ""
-            classes.append((package, name))
-            read_delim(self.buffer, 2)
-        return classes
-
-    def _parse_prop_description(self, classes, section_end):
-        """Parse the description of each property = the third segment"""
-        read_delim(self.buffer, 6)
-        object_info = []
-        while self.buffer.tell() < section_end:
-            class_idx = read_int(self.buffer) - 1
-            class_type = classes[class_idx]
-            read_delim(self.buffer, 2)
-            indices = [x - 1 for x in read_int(self.buffer, 2)]
-            obj_id = read_int(self.buffer)
-            object_info.append((class_type, indices, obj_id))
-        return object_info
-
-    def _parse_properties(self, section_end):
-        """
-        Parse the actual values of the attributes == segments 2 and 4
-        """
-        read_delim(self.buffer, 2)
-        props = []
-        while self.buffer.tell() < section_end:
-            n_props = read_int(self.buffer)
-            d = parse_prop(n_props, self.buffer, self.names, self.heap)
-            if not d:  # Empty dictionary
-                break
-            props.append(d)
-            # Move to next 8-byte aligned offset
-            self.buffer.seek(self.buffer.tell() + self.buffer.tell() % 8)
-        return props
-
-    def to_hdf(self, filename):
-        f = h5py.File(filename, mode="w")
-        save_to_hdf(f, "/", self.attrs)
-
-
-def describe(d, indent=0, width=4, out=None):
-    for key, value in d.items():
-        print(f'{"": <{width * indent}}' + str(key), file=out)
-        if isinstance(value, dict):
-            describe(value, indent + 1, out=out)
-        elif isinstance(value, np.ndarray):
-            print(
-                f'{"": <{width * (indent + 1)}} {value.shape} array '
-                f"of type {value.dtype}",
-                file=out,
-            )
-        elif isinstance(value, scipy.sparse.csc.csc_matrix):
-            print(
-                f'{"": <{width * (indent + 1)}} {value.shape} '
-                f"sparse matrix of type {value.dtype}",
-                file=out,
-            )
-        elif isinstance(value, Iterable) and not isinstance(value, str):
-            print(
-                f'{"": <{width * (indent + 1)}} {type(value)} of len ' f"{len(value)}",
-                file=out,
-            )
-        else:
-            print(f'{"": <{width * (indent + 1)}} {value}', file=out)
-
-
-def parse_prop(n_props, buff, names, heap):
-    d = dict()
-    for i in range(n_props):
-        name_idx, flag, heap_idx = read_int(buff, 3)
-        if flag not in [0, 1, 2] and name_idx == 0:
-            n_props = flag
-            buff.seek(buff.tell() - 1)  # go back on one byte
-            d = parse_prop(n_props, buff, names, heap)
-        else:
-            item_name = names[name_idx - 1]
-            if flag == 0:
-                d[item_name] = names[heap_idx]
-            elif flag == 1:
-                d[item_name] = heap[heap_idx + 2]  # Todo: what is the heap?
-            elif flag == 2:
-                assert 0 <= heap_idx <= 1, (
-                    "Boolean flag has a value other " "than 0 or 1 "
-                )
-                d[item_name] = bool(heap_idx)
-            else:
-                raise ValueError(
-                    "unknown flag {} for property {} with heap "
-                    "index {}".format(flag, item_name, heap_idx)
-                )
-    return d
-
-
-def is_object(x):
-    """Checking object dtype for structured numpy arrays"""
-    if x.dtype.names is not None and len(x.dtype.names) > 1:  # Complex obj
-        return all(x.dtype[ix] == np.object for ix in range(len(x.dtype)))
-    else:  # simple object
-        return x.dtype == np.object
-
-
-def flatten_obj(arr):
-    # TODO turn structured arrays into nested dicts of lists rather that
-    #  lists of dicts
-    if isinstance(arr, np.ndarray):
-        if arr.dtype.names:
-            arrdict = dict()
-            for fieldname in arr.dtype.names:
-                arrdict[fieldname] = flatten_obj(arr[fieldname])
-            arr = arrdict
-        elif arr.dtype == np.object and arr.ndim == 0:
-            arr = flatten_obj(arr[()])
-        elif arr.dtype == np.object and arr.ndim > 0:
-            try:
-                arr = np.stack(arr)
-                if arr.dtype.names:
-                    d = {k: flatten_obj(arr[k]) for k in arr.dtype.names}
-                    arr = d
-            except:
-                arr = [flatten_obj(x) for x in arr.tolist()]
-    elif isinstance(arr, dict):
-        arr = {k: flatten_obj(v) for k, v in arr.items()}
-    elif isinstance(arr, list):
-        try:
-            arr = flatten_obj(np.stack(arr))
-        except:
-            arr = [flatten_obj(x) for x in arr]
-    return arr
-
-
-def save_to_hdf(h5file, path, dic):
-    """
-    Saving a MATLAB object to HDF5
-    """
-    if isinstance(dic, list):
-        dic = {str(i): v for i, v in enumerate(dic)}
-    for key, item in dic.items():
-        if isinstance(item, (int, float, str)):
-            h5file[path].attrs.create(key, item)
-        elif isinstance(item, list):
-            if len(item) == 0 and path + key not in h5file:  # empty list empty group
-                h5file.create_group(path + key)
-            if all(isinstance(x, (int, float, str)) for x in item):
-                if path not in h5file:
-                    h5file.create_group(path)
-                h5file[path].attrs.create(key, item)
-            else:
-                if path + key not in h5file:
-                    h5file.create_group(path + key)
-                save_to_hdf(
-                    h5file, path + key + "/", {str(i): x for i, x in enumerate(item)}
-                )
-        elif isinstance(item, scipy.sparse.csc.csc_matrix):
-            try:
-                h5file.create_dataset(
-                    path + key, data=item.todense(), compression="gzip"
-                )
-            except Exception as e:
-                print(path + key)
-                raise e
-        elif isinstance(item, (np.ndarray, np.int64, np.float64)):
-            if item.dtype == np.dtype("<U1"):  # Strings to 'S' type for HDF5
-                item = item.astype("S")
-            try:
-                h5file.create_dataset(path + key, data=item, compression="gzip")
-            except Exception as e:
-                print(path + key)
-                raise e
-        elif isinstance(item, dict):
-            if path + key not in h5file:
-                h5file.create_group(path + key)
-            save_to_hdf(h5file, path + key + "/", item)
-        elif item is None:
-            continue
-        else:
-            raise ValueError(f"Cannot save {type(item)} type at key {path + key}")
-
-
-## NOT YET FULLY IMPLEMENTED!
-
-
-class _Info:
-    def __init__(self, info):
-        self.info = info
-        self._identity = None
-
-    def __getitem__(self, item):
-        val = self.info[item]
-        if val.shape[0] == 1:
-            val = val[0]
-        if 0 in val[1].shape:
-            val = val[0]
-        if isinstance(val, scipy.sparse.csc.csc_matrix):
-            return np.asarray(val.todense())
-        if val.dtype == np.dtype("O"):
-            # 3d "sparse matrix"
-            if all(isinstance(x, scipy.sparse.csc.csc_matrix) for x in val):
-                val = np.array([x.todense() for x in val])
-            # TODO: The actual object data
-        equality = val[0] == val[1]
-        if isinstance(equality, scipy.sparse.csc.csc_matrix):
-            equality = equality.todense()
-        if equality.all():
-            val = val[0]
-        return np.squeeze(val)
-
-    @property
-    def categories(self):
-        return self.info.dtype.names
-
-
-class TrapInfo(_Info):
-    def __init__(self, info):
-        """
-        The information on all of the traps in a given position.
-
-        :param info: The TrapInfo structure, can be found in the heap of
-        the CTimelapse at index 7
-        """
-        super().__init__(info)
-
-
-class CellInfo(_Info):
-    def __init__(self, info):
-        """
-        The extracted information of all cells in a given position.
-        :param info: The CellInfo structure, can be found in the heap
-        of the CTimelapse at index 15.
-        """
-        super().__init__(info)
-
-    @property
-    def identity(self):
-        if self._identity is None:
-            self._identity = pd.DataFrame(
-                zip(self["trapNum"], self["cellNum"]), columns=["trapNum", "cellNum"]
-            )
-        return self._identity
-
-    def index(self, trapNum, cellNum):
-        query = "trapNum=={} and cellNum=={}".format(trapNum, cellNum)
-        try:
-            result = self.identity.query(query).index[0]
-        except Exception as e:
-            print(query)
-            raise e
-        return result
-
-    @property
-    def nucEstConv1(self):
-        return np.asarray(self.info["nuc_est_conv"][0][0].todense())
-
-    @property
-    def nucEstConv2(self):
-        return np.asarray(self.info["nuc_est_conv"][0][1].todense())
-
-    @property
-    def mothers(self):
-        return np.where((self["births"] != 0).any(axis=1))[0]
-
-    def daughters(self, mother_index):
-        """
-        Get daughters of cell with index `mother_index`.
-
-        :param mother_index: the index of the mother within the data. This is
-        different from the mother's cell/trap identity.
-        """
-        daughter_ids = np.unique(self["daughterLabel"][mother_index]).tolist()
-        daughter_ids.remove(0)
-        mother_trap = self.identity["trapNum"].loc[mother_index]
-        daughters = [self.index(mother_trap, cellNum) for cellNum in daughter_ids]
-        return daughters
-
-
-def _todict(matobj):
-    """
-    A recursive function which constructs from matobjects nested dictionaries
-    """
-    if not hasattr(matobj, "_fieldnames"):
-        return matobj
-    d = {}
-    for strg in matobj._fieldnames:
-        elem = matobj.__dict__[strg]
-        if isinstance(elem, matlab.mio5_params.mat_struct):
-            d[strg] = _todict(elem)
-        elif isinstance(elem, np.ndarray):
-            d[strg] = _toarray(elem)
-        else:
-            d[strg] = elem
-    return d
-
-
-def _toarray(ndarray):
-    """
-    A recursive function which constructs ndarray from cellarrays
-    (which are loaded as numpy ndarrays), recursing into the elements
-    if they contain matobjects.
-    """
-    if ndarray.dtype != "float64":
-        elem_list = []
-        for sub_elem in ndarray:
-            if isinstance(sub_elem, matlab.mio5_params.mat_struct):
-                elem_list.append(_todict(sub_elem))
-            elif isinstance(sub_elem, np.ndarray):
-                elem_list.append(_toarray(sub_elem))
-            else:
-                elem_list.append(sub_elem)
-        return np.array(elem_list)
-    else:
-        return ndarray
-
-
-from pathlib import Path
-
-
-class Strain:
-    """The cell info for all the positions of a strain."""
-
-    def __init__(self, origin, strain):
-        self.origin = Path(origin)
-        self.files = [x for x in origin.iterdir() if strain in str(x)]
-        self.cts = [matObject(x) for x in self.files]
-        self.cinfos = [CellInfo(x.heap[15]) for x in self.cts]
-        self._identity = None
-
-    def __getitem__(self, item):
-        try:
-            return np.concatenate([c[item] for c in self.cinfos])
-        except ValueError:  # If first axis is the channel
-            return np.concatenate([c[item] for c in self.cinfos], axis=1)
-
-    @property
-    def categories(self):
-        return set.union(*[set(c.categories) for c in self.cinfos])
-
-    @property
-    def identity(self):
-        if self._identity is None:
-            identities = []
-            for pos_id, cinfo in enumerate(self.cinfos):
-                identity = cinfo.identity
-                identity["position"] = pos_id
-                identities.append(identity)
-            self._identity = pd.concat(identities, ignore_index=True)
-        return self._identity
-
-    def index(self, posNum, trapNum, cellNum):
-        query = "position=={} and trapNum=={} and cellNum=={}".format(
-            posNum, trapNum, cellNum
-        )
-        try:
-            result = self.identity.query(query).index[0]
-        except Exception as e:
-            raise e
-        return result
-
-    @property
-    def mothers(self):
-        # At least two births are needed to be considered a mother cell
-        return np.where(np.count_nonzero(self["births"], axis=1) > 3)[0]
-
-    def daughters(self, mother_index):
-        """
-        Get daughters of cell with index `mother_index`.
-
-        :param mother_index: the index of the mother within the data. This is
-        different from the mother's pos/trap/cell identity.
-        """
-        daughter_ids = np.unique(self["daughterLabel"][mother_index]).tolist()
-        if 0 in daughter_ids:
-            daughter_ids.remove(0)
-        mother_pos_trap = self.identity[["position", "trapNum"]].loc[mother_index]
-        daughters = []
-        for cellNum in daughter_ids:
-            try:
-                daughters.append(self.index(*mother_pos_trap, cellNum))
-            except IndexError:
-                continue
-        return daughters

aliby/io/omero.py

-import h5py
-import omero
-from omero.gateway import BlitzGateway
-from aliby.experiment import get_data_lazy
-from aliby.cells import CellsHDF
-
-
-class Argo:
-    # TODO use the one in extraction?
-    def __init__(
-        self, host="islay.bio.ed.ac.uk", username="upload", password="***REMOVED***"
-    ):
-        self.conn = None
-        self.host = host
-        self.username = username
-        self.password = password
-
-    def get_meta(self):
-        pass
-
-    def __enter__(self):
-        self.conn = BlitzGateway(
-            host=self.host, username=self.username, passwd=self.password
-        )
-        self.conn.connect()
-        return self
-
-    def __exit__(self, *exc):
-        self.conn.close()
-        return False
-
-
-class Dataset(Argo):
-    def __init__(self, expt_id, **server_info):
-        super().__init__(**server_info)
-        self.expt_id = expt_id
-        self._files = None
-
-    @property
-    def dataset(self):
-        return self.conn.getObject("Dataset", self.expt_id)
-
-    @property
-    def name(self):
-        return self.dataset.getName()
-
-    @property
-    def date(self):
-        return self.dataset.getDate()
-
-    @property
-    def unique_name(self):
-        return "_".join((self.date.strftime("%Y_%m_%d").replace("/", "_"), self.name))
-
-    def get_images(self):
-        return {im.getName(): im.getId() for im in self.dataset.listChildren()}
-
-    @property
-    def files(self):
-        if self._files is None:
-            self._files = {
-                x.getFileName(): x
-                for x in self.dataset.listAnnotations()
-                if isinstance(x, omero.gateway.FileAnnotationWrapper)
-            }
-        if not len(self._files):
-            raise Exception("Exception:Metadata: Experiment has no annotation files.")
-        return self._files
-
-    @property
-    def tags(self):
-        if self._tags is None:
-            self._tags = {
-                x.getName(): x
-                for x in self.dataset.listAnnotations()
-                if isinstance(x, omero.gateway.TagAnnotationWrapper)
-            }
-        return self._tags
-
-    def cache_logs(self, root_dir):
-        for name, annotation in self.files.items():
-            filepath = root_dir / annotation.getFileName().replace("/", "_")
-            if str(filepath).endswith("txt") and not filepath.exists():
-                # Save only the text files
-                with open(str(filepath), "wb") as fd:
-                    for chunk in annotation.getFileInChunks():
-                        fd.write(chunk)
-        return True
-
-
-class Image(Argo):
-    def __init__(self, image_id, **server_info):
-        super().__init__(**server_info)
-        self.image_id = image_id
-        self._image_wrap = None
-
-    @property
-    def image_wrap(self):
-        # TODO check that it is alive/ connected
-        if self._image_wrap is None:
-            self._image_wrap = self.conn.getObject("Image", self.image_id)
-        return self._image_wrap
-
-    @property
-    def name(self):
-        return self.image_wrap.getName()
-
-    @property
-    def data(self):
-        return get_data_lazy(self.image_wrap)
-
-    @property
-    def metadata(self):
-        meta = dict()
-        meta["size_x"] = self.image_wrap.getSizeX()
-        meta["size_y"] = self.image_wrap.getSizeY()
-        meta["size_z"] = self.image_wrap.getSizeZ()
-        meta["size_c"] = self.image_wrap.getSizeC()
-        meta["size_t"] = self.image_wrap.getSizeT()
-        meta["channels"] = self.image_wrap.getChannelLabels()
-        meta["name"] = self.image_wrap.getName()
-        return meta
-
-
-class Cells(CellsHDF):
-    def __init__(self, filename):
-        file = h5py.File(filename, "r")
-        super().__init__(file)
-
-    def __enter__(self):
-        return self
-
-    def __exit__(self, *exc):
-        self.close
-        return False

aliby/io/utils.py

-import re
-import struct
-
-
-def clean_ascii(text):
-    return re.sub(r'[^\x20-\x7F]', '.', text)
-
-
-def xxd(x, start=0, stop=None):
-    if stop is None:
-        stop = len(x)
-    for i in range(start, stop, 8):
-        # Row number
-        print("%04d" % i, end="   ")
-        # Hexadecimal bytes
-        for r in range(i, i + 8):
-            print("%02x" % x[r], end="")
-            if (r + 1) % 4 == 0:
-                print("  ", end="")
-        # ASCII
-        print("   ", clean_ascii(x[i:i + 8].decode('utf-8', errors='ignore')),
-              "   ", end="")
-        # Int32
-        print('{:>10} {:>10}'.format(*struct.unpack('II', x[i: i + 8])),
-              end="   ")
-        print("")  # Newline
-    return
-
-
-# Buffer reading functions
-def read_int(buffer, n=1):
-    res = struct.unpack('I' * n, buffer.read(4 * n))
-    if n == 1:
-        res = res[0]
-    return res
-
-
-def read_string(buffer):
-    return ''.join([x.decode() for x in iter(lambda: buffer.read(1), b'\x00')])
-
-
-def read_delim(buffer, n):
-    delim = read_int(buffer, n)
-    assert all([x == 0 for x in delim]), "Unknown nonzero value in delimiter"

aliby/multiexperiment.py

-from pathos.multiprocessing import Pool
-
-from aliby.pipeline import PipelineParameters, Pipeline
-
-
-class MultiExp:
-    """
-    Manages cases when you need to segment several different experiments with a single
-    position (e.g. pH calibration).
-    """
-
-    def __init__(self, expt_ids, npools=8, *args, **kwargs):
-
-        self.expt_ids = expt_ids
-
-    def run(self):
-        run_expt = lambda expt: Pipeline(
-            PipelineParameters.default(general={"expt_id": expt, "distributed": 0})
-        ).run()
-        with Pool(npools) as p:
-            results = p.map(lambda x: self.create_pipeline(x), self.exp_ids)
-
-    @classmethod
-    def default(self):
-        return cls(expt_ids=list(range(20448, 20467 + 1)))

aliby/pipeline.py

-"""
-Pipeline and chaining elements.
-"""
-import logging
-import os
-from abc import ABC, abstractmethod
-from typing import List
-from pathlib import Path
-import traceback
-
-import itertools
-import yaml
-from tqdm import tqdm
-from time import perf_counter
-
-import numpy as np
-import pandas as pd
-from pathos.multiprocessing import Pool
-
-from aliby.experiment import MetaData
-from aliby.haystack import initialise_tf
-from aliby.baby_client import BabyRunner, BabyParameters
-from aliby.tile.tiler import Tiler, TilerParameters
-from aliby.io.omero import Dataset, Image
-from agora.abc import ParametersABC, ProcessABC
-from agora.io.writer import TilerWriter, BabyWriter
-from agora.io.signal import Signal
-from extraction.core.extractor import Extractor, ExtractorParameters
-from extraction.core.functions.defaults import exparams_from_meta
-from postprocessor.core.processor import PostProcessor, PostProcessorParameters
-
-
-class PipelineParameters(ParametersABC):
-    def __init__(self, general, tiler, baby, extraction, postprocessing):
-        self.general = general
-        self.tiler = tiler
-        self.baby = baby
-        self.extraction = extraction
-        self.postprocessing = postprocessing
-
-    @classmethod
-    def default(
-        cls,
-        general={},
-        tiler={},
-        baby={},
-        extraction={},
-        postprocessing={},
-    ):
-        """
-        Load unit test experiment
-        :expt_id: Experiment id
-        :directory: Output directory
-
-        Provides default parameters for the entire pipeline. This downloads the logfiles and sets the default
-        timepoints and extraction parameters from there.
-        """
-        expt_id = general.get("expt_id", 19993)
-        directory = Path(general.get("directory", "../data"))
-        with Dataset(int(expt_id), **general.get("server_info")) as conn:
-            directory = directory / conn.unique_name
-            if not directory.exists():
-                directory.mkdir(parents=True)
-                # Download logs to use for metadata
-            conn.cache_logs(directory)
-        meta = MetaData(directory, None).load_logs()
-        tps = meta["time_settings/ntimepoints"][0]
-        defaults = {
-            "general": dict(
-                id=expt_id,
-                distributed=0,
-                tps=tps,
-                directory=directory,
-                strain="",
-                earlystop=dict(
-                    min_tp=180,
-                    thresh_pos_clogged=0.3,
-                    thresh_trap_clogged=7,
-                    ntps_to_eval=5,
-                ),
-            )
-        }
-        defaults["tiler"] = TilerParameters.default().to_dict()
-        defaults["baby"] = BabyParameters.default().to_dict()
-        defaults["extraction"] = exparams_from_meta(meta)
-        defaults["postprocessing"] = PostProcessorParameters.default().to_dict()
-        for k in defaults.keys():
-            exec("defaults[k].update(" + k + ")")
-        return cls(**{k: v for k, v in defaults.items()})
-
-    def load_logs(self):
-        parsed_flattened = parse_logfiles(self.log_dir)
-        return parsed_flattened
-
-
-class Pipeline(ProcessABC):
-    """
-    A chained set of Pipeline elements connected through pipes.
-    """
-
-    # Tiling, Segmentation,Extraction and Postprocessing should use their own default parameters
-
-    # Early stop for clogging
-    earlystop = {
-        "min_tp": 180,
-        "thresh_pos_clogged": 0.3,
-        "thresh_trap_clogged": 7,
-        "ntps_to_eval": 5,
-    }
-
-    def __init__(self, parameters: PipelineParameters):
-        super().__init__(parameters)
-        self.store = self.parameters.general["directory"]
-
-    @classmethod
-    def from_yaml(cls, fpath):
-        # This is just a convenience function, think before implementing
-        # for other processes
-        return cls(parameters=PipelineParameters.from_yaml(fpath))
-
-    def run(self):
-        # Config holds the general information, use in main
-        # Steps holds the description of tasks with their parameters
-        # Steps: all holds general tasks
-        # steps: strain_name holds task for a given strain
-        config = self.parameters.to_dict()
-        expt_id = config["general"]["id"]
-        distributed = config["general"]["distributed"]
-        strain_filter = config["general"]["strain"]
-        root_dir = config["general"]["directory"]
-        root_dir = Path(root_dir)
-
-        print("Searching OMERO")
-        # Do all initialis
-        with Dataset(int(expt_id), **self.general["server_info"]) as conn:
-            image_ids = conn.get_images()
-            directory = root_dir / conn.unique_name
-            if not directory.exists():
-                directory.mkdir(parents=True)
-                # Download logs to use for metadata
-            conn.cache_logs(directory)
-
-        # Modify to the configuration
-        self.parameters.general["directory"] = directory
-        config["general"]["directory"] = directory
-
-        # Filter TODO integrate filter onto class and add regex
-        if isinstance(strain_filter, str):
-            image_ids = {
-                k: v for k, v in image_ids.items() if k.startswith(strain_filter)
-            }
-        elif isinstance(strain_filter, int):
-            image_ids = {
-                k: v for i, (k, v) in enumerate(image_ids.items()) if i == strain_filter
-            }
-
-        if distributed != 0:  # Gives the number of simultaneous processes
-            with Pool(distributed) as p:
-                results = p.map(lambda x: self.create_pipeline(x), image_ids.items())
-            return results
-        else:  # Sequential
-            results = []
-            for k, v in image_ids.items():
-                r = self.create_pipeline((k, v))
-                results.append(r)
-
-    def create_pipeline(self, image_id):
-        config = self.parameters.to_dict()
-        name, image_id = image_id
-        general_config = config["general"]
-        session = None
-        earlystop = general_config["earlystop"]
-        try:
-            directory = general_config["directory"]
-            with Image(image_id, **self.general["server_info"]) as image:
-                filename = f"{directory}/{image.name}.h5"
-                try:
-                    os.remove(filename)
-                except:
-                    pass
-
-                # Run metadata first
-                process_from = 0
-                # if True:  # not Path(filename).exists():
-                meta = MetaData(directory, filename)
-                meta.run()
-                tiler = Tiler.from_image(
-                    image, TilerParameters.from_dict(config["tiler"])
-                )
-                # else: TODO add support to continue local experiments?
-                #     tiler = Tiler.from_hdf5(image.data, filename)
-                #     s = Signal(filename)
-                #     process_from = s["/general/None/extraction/volume"].columns[-1]
-                #     if process_from > 2:
-                #         process_from = process_from - 3
-                #         tiler.n_processed = process_from
-
-                writer = TilerWriter(filename)
-                session = initialise_tf(2)
-                runner = BabyRunner.from_tiler(
-                    BabyParameters.from_dict(config["baby"]), tiler
-                )
-                bwriter = BabyWriter(filename)
-
-                # Limit extraction parameters during run using the available channels in tiler
-                av_channels = set((*tiler.channels, "general"))
-                config["extraction"]["tree"] = {
-                    k: v
-                    for k, v in config["extraction"]["tree"].items()
-                    if k in av_channels
-                }
-                config["extraction"]["sub_bg"] = av_channels.intersection(
-                    config["extraction"]["sub_bg"]
-                )
-
-                av_channels_wsub = av_channels.union(
-                    [c + "_bgsub" for c in config["extraction"]["sub_bg"]]
-                )
-                for op in config["extraction"]["multichannel_ops"]:
-                    config["extraction"]["multichannel_ops"][op] = [
-                        x
-                        for x in config["extraction"]["multichannel_ops"]
-                        if len(x[0]) == len(av_channels_wsub.intersection(x[0]))
-                    ]
-                config["extraction"]["multichannel_ops"] = {
-                    k: v
-                    for k, v in config["extraction"]["multichannel_ops"].items()
-                    if len(v)
-                }
-
-                exparams = ExtractorParameters.from_dict(config["extraction"])
-                ext = Extractor.from_tiler(exparams, store=filename, tiler=tiler)
-
-                # RUN
-                tps = general_config["tps"]
-                frac_clogged_traps = 0
-                for i in tqdm(
-                    range(process_from, tps), desc=image.name, initial=process_from
-                ):
-                    if (
-                        frac_clogged_traps < earlystop["thresh_pos_clogged"]
-                        or i < earlystop["min_tp"]
-                    ):
-                        t = perf_counter()
-                        trap_info = tiler.run_tp(i)
-                        logging.debug(f"Timing:Trap:{perf_counter() - t}s")
-                        t = perf_counter()
-                        writer.write(trap_info, overwrite=[])
-                        logging.debug(f"Timing:Writing-trap:{perf_counter() - t}s")
-                        t = perf_counter()
-                        seg = runner.run_tp(i)
-                        logging.debug(f"Timing:Segmentation:{perf_counter() - t}s")
-                        # logging.debug(
-                        #     f"Segmentation failed:Segmentation:{perf_counter() - t}s"
-                        # )
-                        t = perf_counter()
-                        bwriter.write(seg, overwrite=["mother_assign"])
-                        logging.debug(f"Timing:Writing-baby:{perf_counter() - t}s")
-                        t = perf_counter()
-
-                        tmp = ext.run(tps=[i])
-                        logging.debug(f"Timing:Extraction:{perf_counter() - t}s")
-                    else:  # Stop if more than X% traps are clogged
-                        logging.debug(
-                            f"EarlyStop:{earlystop['thresh_pos_clogged']*100}% traps clogged at time point {i}"
-                        )
-                        print(
-                            f"Stopping analysis at time {i} with {frac_clogged_traps} clogged traps"
-                        )
-                        break
-
-                    if (
-                        i > earlystop["min_tp"]
-                    ):  # Calculate the fraction of clogged traps
-                        frac_clogged_traps = self.check_earlystop(filename, earlystop)
-                        logging.debug(f"Quality:Clogged_traps:{frac_clogged_traps}")
-                        print("Frac clogged traps: ", frac_clogged_traps)
-
-                # Run post processing
-                post_proc_params = PostProcessorParameters.from_dict(
-                    self.parameters.postprocessing
-                ).to_dict()
-                PostProcessor(filename, post_proc_params).run()
-                return True
-        except Exception as e:  # bug in the trap getting
-            print(f"Caught exception in worker thread (x = {name}):")
-            # This prints the type, value, and stack trace of the
-            # current exception being handled.
-            traceback.print_exc()
-            print()
-            raise e
-        finally:
-            if session:
-                session.close()
-
-    def check_earlystop(self, filename, es_parameters):
-        s = Signal(filename)
-        df = s["/extraction/general/None/area"]
-        frac_clogged_traps = (
-            df[df.columns[-1 - es_parameters["ntps_to_eval"] : -1]]
-            .dropna(how="all")
-            .notna()
-            .groupby("trap")
-            .apply(sum)
-            .apply(np.mean, axis=1)
-            > es_parameters["thresh_trap_clogged"]
-        ).mean()
-        return frac_clogged_traps

aliby/post_processing.py

-"""
-Post-processing utilities
-
-Notes: I don't have statistics on ranges of radii for each of the knots in
-the radial spline representation, but we regularly extract the average of
-these radii for each cell. So, depending on camera/lens, we get:
-    * 60x evolve: mean radii of 2-14 pixels (and measured areas of 30-750
-    pixels^2)
-    * 60x prime95b: mean radii of 3-24 pixels (and measured areas of 60-2000
-	pixels^2)
-
-And I presume that for a 100x lens we would get an ~5/3 increase over those
-values.
-
-In terms of the current volume estimation method, it's currently only
-implemented in the AnalysisToolbox repository, but it's super simple:
-
-mVol = 4/3*pi*sqrt(mArea/pi).^3
-
-where mArea is simply the sum of pixels for that cell.
-"""
-import matplotlib.pyplot as plt
-import numpy as np
-from mpl_toolkits.mplot3d.art3d import Poly3DCollection
-from scipy import ndimage
-from skimage.morphology import erosion, ball
-from skimage import measure, draw
-
-
-def my_ball(radius):
-    """Generates a ball-shaped structuring element.
-
-    This is the 3D equivalent of a disk.
-    A pixel is within the neighborhood if the Euclidean distance between
-    it and the origin is no greater than radius.
-
-    Parameters
-    ----------
-    radius : int
-        The radius of the ball-shaped structuring element.
-
-    Other Parameters
-    ----------------
-    dtype : data-type
-        The data type of the structuring element.
-
-    Returns
-    -------
-    selem : ndarray
-        The structuring element where elements of the neighborhood
-        are 1 and 0 otherwise.
-    """
-    n = 2 * radius + 1
-    Z, Y, X = np.mgrid[-radius:radius:n * 1j,
-              -radius:radius:n * 1j,
-              -radius:radius:n * 1j]
-    X **= 2
-    Y **= 2
-    Z **= 2
-    X += Y
-    X += Z
-    # s = X ** 2 + Y ** 2 + Z ** 2
-    return X <= radius * radius
-
-def circle_outline(r):
-    return ellipse_perimeter(r, r)
-
-def ellipse_perimeter(x, y):
-    im_shape = int(2*max(x, y) + 1)
-    img = np.zeros((im_shape, im_shape), dtype=np.uint8)
-    rr, cc = draw.ellipse_perimeter(int(im_shape//2), int(im_shape//2),
-                                    int(x), int(y))
-    img[rr, cc] = 1
-    return np.pad(img, 1)
-
-def capped_cylinder(x, y):
-    max_size = (y + 2*x + 2)
-    pixels = np.zeros((max_size, max_size))
-
-    rect_start = ((max_size-x)//2, x + 1)
-    rr, cc = draw.rectangle_perimeter(rect_start, extent=(x, y),
-                                     shape=(max_size, max_size))
-    pixels[rr, cc] = 1
-    circle_centres = [(max_size//2 - 1, x),
-                      (max_size//2 - 1, max_size - x - 1 )]
-    for r, c in circle_centres:
-        rr, cc = draw.circle_perimeter(r, c, (x + 1)//2,
-                                       shape=(max_size, max_size))
-        pixels[rr, cc] = 1
-    pixels = ndimage.morphology.binary_fill_holes(pixels)
-    pixels ^= erosion(pixels)
-    return pixels
-
-def volume_of_sphere(radius):
-    return 4 / 3 * np.pi * radius**3
-
-def plot_voxels(voxels):
-    verts, faces, normals, values = measure.marching_cubes_lewiner(
-        voxels, 0)
-    fig = plt.figure(figsize=(10, 10))
-    ax = fig.add_subplot(111, projection='3d')
-    mesh = Poly3DCollection(verts[faces])
-    mesh.set_edgecolor('k')
-    ax.add_collection3d(mesh)
-    ax.set_xlim(0, voxels.shape[0])
-    ax.set_ylim(0, voxels.shape[1])
-    ax.set_zlim(0, voxels.shape[2])
-    plt.tight_layout()
-    plt.show()
-
-# Volume estimation
-def union_of_spheres(outline, shape='my_ball', debug=False):
-    filled = ndimage.binary_fill_holes(outline)
-    nearest_neighbor = ndimage.morphology.distance_transform_edt(
-        outline == 0) * filled
-    voxels = np.zeros((filled.shape[0], filled.shape[1], max(filled.shape)))
-    c_z = voxels.shape[2] // 2
-    for x,y in zip(*np.where(filled)):
-        radius = nearest_neighbor[(x,y)]
-        if radius > 0:
-            if shape == 'ball':
-                b = ball(radius)
-            elif shape == 'my_ball':
-                b = my_ball(radius)
-            else:
-                raise ValueError(f"{shape} is not an accepted value for "
-                                 f"shape.")
-            centre_b = ndimage.measurements.center_of_mass(b)
-
-            I,J,K = np.ogrid[:b.shape[0], :b.shape[1], :b.shape[2]]
-            voxels[I + int(x - centre_b[0]), J + int(y - centre_b[1]),
-                   K + int(c_z - centre_b[2])] += b
-    if debug:
-        plot_voxels(voxels)
-    return voxels.astype(bool).sum()
-
-def improved_uos(outline, shape='my_ball', debug=False):
-    filled = ndimage.binary_fill_holes(outline)
-    nearest_neighbor = ndimage.morphology.distance_transform_edt(
-        outline == 0) * filled
-    voxels = np.zeros((filled.shape[0], filled.shape[1], max(filled.shape)))
-    c_z = voxels.shape[2] // 2
-
-    while np.any(nearest_neighbor != 0):
-        radius = np.max(nearest_neighbor)
-        x, y = np.argwhere(nearest_neighbor == radius)[0]
-        if shape == 'ball':
-            b = ball(np.ceil(radius))
-        elif shape == 'my_ball':
-            b = my_ball(np.ceil(radius))
-        else:
-            raise ValueError(f"{shape} is not an accepted value for shape")
-        centre_b = ndimage.measurements.center_of_mass(b)
-
-        I, J, K = np.ogrid[:b.shape[0], :b.shape[1], :b.shape[2]]
-        voxels[I + int(x - centre_b[0]), J + int(y - centre_b[1]),
-               K + int(c_z - centre_b[2])] += b
-
-        # Use the central disk of the ball from voxels to get the circle
-        # = 0 if nn[x,y] < r else nn[x,y]
-        rr, cc = draw.circle(x, y, np.ceil(radius), nearest_neighbor.shape)
-        nearest_neighbor[rr, cc] = 0
-    if debug:
-        plot_voxels(voxels)
-    return voxels.astype(bool).sum()
-
-def conical(outline, debug=False):
-    nearest_neighbor = ndimage.morphology.distance_transform_edt(
-        outline == 0) * ndimage.binary_fill_holes(outline)
-    if debug:
-        hf = plt.figure()
-        ha = hf.add_subplot(111, projection='3d')
-
-        X, Y = np.meshgrid(np.arange(nearest_neighbor.shape[0]),
-                           np.arange(nearest_neighbor.shape[1]))
-        ha.plot_surface(X, Y, nearest_neighbor)
-        plt.show()
-    return 4 * nearest_neighbor.sum()
-
-def volume(outline, method='spheres'):
-    if method=='conical':
-        return conical(outline)
-    elif method=='spheres':
-        return union_of_spheres(outline)
-    else:
-        raise ValueError(f"Method {method} not implemented.")
-
-def circularity(outline):
-    pass
\ No newline at end of file

aliby/results.py

-"""Pipeline results classes and utilities"""
-
-
-class SegmentationResults:
-    """
-    Object storing the data from the Segmentation pipeline.
-    Everything is stored as an `AttributeDict`, which is a `defaultdict` where
-    you can get elements as attributes.
-
-    In addition, it implements:
-     - IO functionality (read from file, write to file)
-    """
-    def __init__(self, raw_expt):
-        pass
-
-
-
-
-class CellResults:
-    """
-    Results on a set of cells TODO: what set of cells, how?
-
-    Contains:
-    * cellInf describing which cells are taken into account
-    * annotations on the cell
-    * segmentation maps of the cell TODO: how to define and save this?
-    * trapLocations TODO: why is this not part of cellInf?
-    """
-
-    def __init__(self, cellInf=None, annotations=None, segmentation=None,
-                 trapLocations=None):
-        self._cellInf = cellInf
-        self._annotations = annotations
-        self._segmentation = segmentation
-        self._trapLocations = trapLocations