.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,49 +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 pull
-  - if [ ${var+TRIGGER_PYPI_NAME} ]; then echo "Pipeline triggered by ${TRIGGER_PYPI_NAME}"; poetry update ${TRIGGER_PYPI_NAME}; fi
-  - poetry install -vv
+  # - 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:
-    - 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 version and 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
-  except:
-    changes:
-      - tests/*
-  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.segment import Tiler, TilerParameters
+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 core.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/cells.py

-import logging
-from pathlib import Path, PosixPath
-from time import perf_counter
-from typing import Union
-from itertools import groupby
-from collections.abc import Iterable
-
-from utils_find_1st import find_1st, cmp_equal
-import h5py
-import numpy as np
-from scipy import ndimage
-from scipy.sparse.base import isdense
-
-from aliby.io.matlab import matObject
-from aliby.utils import timed
-from aliby.io.writer import load_complex
-
-
-def cell_factory(store, type="matlab"):
-    if isinstance(store, matObject):
-        return CellsMat(store)
-    if type == "matlab":
-        mat_object = matObject(store)
-        return CellsMat(mat_object)
-    elif type == "hdf5":
-        return CellsHDF(store)
-    else:
-        raise TypeError(
-            "Could not get cells for type {}:" "valid types are matlab and hdf5"
-        )
-
-
-class Cells:
-    """An object that gathers information about all the cells in a given
-    trap.
-    This is the abstract object, used for type testing
-    """
-
-    def __init__(self):
-        pass
-
-    @staticmethod
-    def from_source(source: Union[PosixPath, str], kind: str = None):
-        if isinstance(source, str):
-            source = Path(source)
-        if kind is None:  # Infer kind from filename
-            kind = "matlab" if source.suffix == ".mat" else "hdf5"
-        return cell_factory(source, kind)
-
-    @staticmethod
-    def _asdense(array):
-        if not isdense(array):
-            array = array.todense()
-        return array
-
-    @staticmethod
-    def _astype(array, kind):
-        # Convert sparse arrays if needed and if kind is 'mask' it fills the outline
-        array = Cells._asdense(array)
-        if kind == "mask":
-            array = ndimage.binary_fill_holes(array).astype(int)
-        return array
-
-    @classmethod
-    def hdf(cls, fpath):
-        return CellsHDF(fpath)
-
-    @classmethod
-    def mat(cls, path):
-        return CellsMat(matObject(store))
-
-
-class CellsHDF(Cells):
-    def __init__(self, filename, path="cell_info"):
-        self.filename = filename
-        self.cinfo_path = path
-        self._edgem_indices = None
-        self._edgemasks = None
-        self._tile_size = None
-
-    def __getitem__(self, item):
-        if item == "edgemasks":
-            return self.edgemasks
-        _item = "_" + item
-        if not hasattr(self, _item):
-            setattr(self, _item, self._fetch(item))
-        return getattr(self, _item)
-
-    def _get_idx(self, cell_id, trap_id):
-        return (self["cell_label"] == cell_id) & (self["trap"] == trap_id)
-
-    def _fetch(self, path):
-        with h5py.File(self.filename, mode="r") as f:
-            return f[self.cinfo_path][path][()]
-
-    @property
-    def ntraps(self):
-        with h5py.File(self.filename, mode="r") as f:
-            return len(f["/trap_info/trap_locations"][()])
-
-    @property
-    def traps(self):
-        return list(set(self["trap"]))
-
-    @property
-    def tile_size(self):  # TODO read from metadata
-        if self._tile_size is None:
-            with h5py.File(self.filename, mode="r") as f:
-                self._tile_size == f["trap_info/tile_size"][0]
-        return self._tile_size
-
-    @property
-    def edgem_indices(self):
-        if self._edgem_indices is None:
-            edgem_path = "edgemasks/indices"
-            self._edgem_indices = load_complex(self._fetch(edgem_path))
-        return self._edgem_indices
-
-    @property
-    def edgemasks(self):
-        if self._edgemasks is None:
-            edgem_path = "edgemasks/values"
-            self._edgemasks = self._fetch(edgem_path)
-
-        return self._edgemasks
-
-    def _edgem_where(self, cell_id, trap_id):
-        ix = trap_id + 1j * cell_id
-        return find_1st(self.edgem_indices == ix, True, cmp_equal)
-
-    @property
-    def labels(self):
-        """
-        Return all cell labels in object
-        We use mother_assign to list traps because it is the only propriety that appears even
-        when no cells are found"""
-        return [self.labels_in_trap(trap) for trap in self.traps]
-
-    def where(self, cell_id, trap_id):
-        """
-        Returns
-        Parameters
-        ----------
-            cell_id: int
-                Cell index
-            trap_id: int
-                Trap index
-
-        Returns
-        ----------
-            indices int array
-            boolean mask array
-            edge_ix int array
-        """
-        indices = self._get_idx(cell_id, trap_id)
-        edgem_ix = self._edgem_where(cell_id, trap_id)
-        return (
-            self["timepoint"][indices],
-            indices,
-            edgem_ix,
-        )  # FIXME edgem_ix makes output different to matlab's Cell
-
-    def outline(self, cell_id, trap_id):
-        times, indices, cell_ix = self.where(cell_id, trap_id)
-        return times, self["edgemasks"][cell_ix, times]
-
-    def mask(self, cell_id, trap_id):
-        times, outlines = self.outline(cell_id, trap_id)
-        return times, np.array(
-            [ndimage.morphology.binary_fill_holes(o) for o in outlines]
-        )
-
-    def at_time(self, timepoint, kind="mask"):
-        ix = self["timepoint"] == timepoint
-        cell_ix = self["cell_label"][ix]
-        traps = self["trap"][ix]
-        indices = traps + 1j * cell_ix
-        choose = np.in1d(self.edgem_indices, indices)
-        edgemasks = self["edgemasks"][choose, timepoint]
-        masks = [
-            self._astype(edgemask, kind) for edgemask in edgemasks if edgemask.any()
-        ]
-        return self.group_by_traps(traps, masks)
-
-    def group_by_traps(self, traps, data):
-        # returns a dict with traps as keys and labels as value
-        iterator = groupby(zip(traps, data), lambda x: x[0])
-        d = {key: [x[1] for x in group] for key, group in iterator}
-        d = {i: d.get(i, []) for i in self.traps}
-        return d
-
-    def labels_in_trap(self, trap_id):
-        # Return set of cell ids in a trap.
-        return set((self["cell_label"][self["trap"] == trap_id]))
-
-    def labels_at_time(self, timepoint):
-        labels = self["cell_label"][self["timepoint"] == timepoint]
-        traps = self["trap"][self["timepoint"] == timepoint]
-        return self.group_by_traps(traps, labels)
-
-
-class CellsMat(Cells):
-    def __init__(self, mat_object):
-        super(CellsMat, self).__init__()
-        # TODO add __contains__ to the matObject
-        timelapse_traps = mat_object.get(
-            "timelapseTrapsOmero", mat_object.get("timelapseTraps", None)
-        )
-        if timelapse_traps is None:
-            raise NotImplementedError(
-                "Could not find a timelapseTraps or "
-                "timelapseTrapsOmero object. Cells "
-                "from cellResults not implemented"
-            )
-        else:
-            self.trap_info = timelapse_traps["cTimepoint"]["trapInfo"]
-
-            if isinstance(self.trap_info, list):
-                self.trap_info = {
-                    k: list([res.get(k, []) for res in self.trap_info])
-                    for k in self.trap_info[0].keys()
-                }
-
-    def where(self, cell_id, trap_id):
-        times, indices = zip(
-            *[
-                (tp, np.where(cell_id == x)[0][0])
-                for tp, x in enumerate(self.trap_info["cellLabel"][:, trap_id].tolist())
-                if np.any(cell_id == x)
-            ]
-        )
-        return times, indices
-
-    def outline(self, cell_id, trap_id):
-        times, indices = self.where(cell_id, trap_id)
-        info = self.trap_info["cell"][times, trap_id]
-
-        def get_segmented(cell, index):
-            if cell["segmented"].ndim == 0:
-                return cell["segmented"][()].todense()
-            else:
-                return cell["segmented"][index].todense()
-
-        segmentation_outline = [
-            get_segmented(cell, idx) for idx, cell in zip(indices, info)
-        ]
-        return times, np.array(segmentation_outline)
-
-    def mask(self, cell_id, trap_id):
-        times, outlines = self.outline(cell_id, trap_id)
-        return times, np.array(
-            [ndimage.morphology.binary_fill_holes(o) for o in outlines]
-        )
-
-    def at_time(self, timepoint, kind="outline"):
-
-        """Returns the segmentations for all the cells at a given timepoint.
-
-        FIXME: this is extremely hacky and accounts for differently saved
-            results in the matlab object. Deprecate ASAP.
-        """
-        # Case 1: only one cell per trap: trap_info['cell'][timepoint] is a
-        # structured array
-        if isinstance(self.trap_info["cell"][timepoint], dict):
-            segmentations = [
-                self._astype(x, "outline")
-                for x in self.trap_info["cell"][timepoint]["segmented"]
-            ]
-        # Case 2: Multiple cells per trap: it becomes a list of arrays or
-        # dictionaries,  one for each trap
-        # Case 2.1 : it's a dictionary
-        elif isinstance(self.trap_info["cell"][timepoint][0], dict):
-            segmentations = []
-            for x in self.trap_info["cell"][timepoint]:
-                seg = x["segmented"]
-                if not isinstance(seg, np.ndarray):
-                    seg = [seg]
-                segmentations.append([self._astype(y, "outline") for y in seg])
-        # Case 2.2 : it's an array
-        else:
-            segmentations = [
-                [self._astype(y, type) for y in x["segmented"]] if x.ndim != 0 else []
-                for x in self.trap_info["cell"][timepoint]
-            ]
-            # Return dict for compatibility with hdf5 output
-        return {i: v for i, v in enumerate(segmentations)}
-
-    def labels_at_time(self, tp):
-        labels = self.trap_info["cellLabel"]
-        labels = [_aslist(x) for x in labels[tp]]
-        labels = {i: [lbl for lbl in lblset] for i, lblset in enumerate(labels)}
-        return labels
-
-    @property
-    def ntraps(self):
-        return len(self.trap_info["cellLabel"][0])
-
-    @property
-    def tile_size(self):
-        pass
-
-
-class ExtractionRunner:
-    """An object to run extraction of fluorescence, and general data out of
-    segmented data.
-
-    Configure with what extraction we want to run.
-    Cell selection criteria.
-    Filtering criteria.
-    """
-
-    def __init__(self, tiler, cells):
-        pass
-
-    def run(self, keys, store, **kwargs):
-        pass
-
-
-def _aslist(x):
-    if isinstance(x, Iterable):
-        if hasattr(x, "tolist"):
-            x = x.tolist()
-    else:
-        x = [x]
-    return x

aliby/core.py

-"""Barebones implementation of the structure/organisation of experiments."""
-
-
-class Experiment:
-    def __init__(self):
-        self.strains = dict()
-        self._metadata = None
-
-    def add_strains(self, name, strain):
-        self.strains[name] = strain
-
-
-class Strain:
-    def __init__(self):
-        self.positions = dict()
-
-    def add_position(self, name, position):
-        self.positions[name] = position
-
-
-class Position:
-    def __init__(self):
-        self.traps = []
-
-    def add_trap(self, trap):
-        self.traps.append(trap)
-
-
-class Trap:  # TODO Name this Tile?
-    def __init__(self):
-        self.cells = []
-
-    def add_cell(self, cell):
-        self.cells.append(cell)

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 aliby.timelapse import TimelapseOMERO, TimelapseLocal
-from aliby.utils import accumulate
-
-from aliby.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 aliby.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/metadata_parser.py

-"""
-Parse microscopy log files according to specified JSON grammars.
-Produces dictionary to include in HDF5
-"""
-import glob
-import os
-import numpy as np
-import pandas as pd
-from datetime import datetime, timezone
-from pytz import timezone
-
-from logfile_parser import Parser
-
-# Paradigm: able to do something with all datatypes present in log files,
-# then pare down on what specific information is really useful later.
-
-# Needed because HDF5 attributes do not support dictionaries
-def flatten_dict(nested_dict, separator="/"):
-    """
-    Flattens nested dictionary
-    """
-    df = pd.json_normalize(nested_dict, sep=separator)
-    return df.to_dict(orient="records")[0]
-
-
-# Needed because HDF5 attributes do not support datetime objects
-# Takes care of time zones & daylight saving
-def datetime_to_timestamp(time, locale="Europe/London"):
-    """
-    Convert datetime object to UNIX timestamp
-    """
-    return timezone(locale).localize(time).timestamp()
-
-
-def find_file(root_dir, regex):
-    file = glob.glob(os.path.join(str(root_dir), regex))
-    if len(file) != 1:
-        return None
-    else:
-        return file[0]
-
-
-# TODO: re-write this as a class if appropriate
-# WARNING: grammars depend on the directory structure of a locally installed
-# logfile_parser repo
-def parse_logfiles(
-    root_dir,
-    acq_grammar="multiDGUI_acq_format.json",
-    log_grammar="multiDGUI_log_format.json",
-):
-    """
-    Parse acq and log files depending on the grammar specified, then merge into
-    single dict.
-    """
-    # Both acq and log files contain useful information.
-    # ACQ_FILE = 'flavin_htb2_glucose_long_ramp_DelftAcq.txt'
-    # LOG_FILE = 'flavin_htb2_glucose_long_ramp_Delftlog.txt'
-    log_parser = Parser(log_grammar)
-    try:
-        log_file = find_file(root_dir, "*log.txt")
-    except:
-        raise ValueError("Experiment log file not found.")
-    with open(log_file, "r") as f:
-        log_parsed = log_parser.parse(f)
-
-    acq_parser = Parser(acq_grammar)
-    try:
-        acq_file = find_file(root_dir, "*[Aa]cq.txt")
-    except:
-        raise ValueError("Experiment acq file not found.")
-    with open(acq_file, "r") as f:
-        acq_parsed = acq_parser.parse(f)
-
-    parsed = {**acq_parsed, **log_parsed}
-
-    for key, value in parsed.items():
-        if isinstance(value, datetime):
-            parsed[key] = datetime_to_timestamp(value)
-
-    parsed_flattened = flatten_dict(parsed)
-    for k, v in parsed_flattened.items():
-        if isinstance(v, list):
-            parsed_flattened[k] = [0 if el is None else el for el in v]
-
-    return parsed_flattened

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/signal.py

-import numpy as np
-from copy import copy
-from itertools import accumulate
-
-from numpy import ndarray
-
-# from more_itertools import first_true
-
-import h5py
-import pandas as pd
-from utils_find_1st import find_1st, cmp_larger
-
-from aliby.io.base import BridgeH5
-
-
-class Signal(BridgeH5):
-    """
-    Class that fetches data from the hdf5 storage for post-processing
-    """
-
-    def __init__(self, file):
-        super().__init__(file, flag=None)
-
-        self.names = ["experiment", "position", "trap"]
-
-    @staticmethod
-    def add_name(df, name):
-        df.name = name
-        return df
-
-    def mothers(self, signal, cutoff=0.8):
-        df = self[signal]
-        get_mothers = lambda df: df.loc[df.notna().sum(axis=1) > df.shape[1] * cutoff]
-        if isinstance(df, pd.DataFrame):
-            return get_mothers(df)
-        elif isinstance(df, list):
-            return [get_mothers(d) for d in df]
-
-    def __getitem__(self, dsets):
-
-        if isinstance(dsets, str) and (
-            dsets.startswith("postprocessing")
-            or dsets.startswith("/postprocessing")
-            or dsets.endswith("imBackground")
-        ):
-            df = self.get_raw(dsets)
-
-        elif isinstance(dsets, str):
-            df = self.apply_prepost(dsets)
-
-        elif isinstance(dsets, list):
-            is_bgd = [dset.endswith("imBackground") for dset in dsets]
-            assert sum(is_bgd) == 0 or sum(is_bgd) == len(
-                dsets
-            ), "Trap data and cell data can't be mixed"
-            with h5py.File(self.filename, "r") as f:
-                return [self.add_name(self.apply_prepost(dset), dset) for dset in dsets]
-
-        return self.add_name(df, dsets)
-
-    def apply_prepost(self, dataset: str):
-        merges = self.get_merges()
-        with h5py.File(self.filename, "r") as f:
-            df = self.dset_to_df(f, dataset)
-
-            merged = df
-            if merges.any():
-                # Split in two dfs, one with rows relevant for merging and one without them
-                mergable_ids = pd.MultiIndex.from_arrays(
-                    np.unique(merges.reshape(-1, 2), axis=0).T,
-                    names=df.index.names,
-                )
-                merged = self.apply_merge(df.loc[mergable_ids], merges)
-
-                nonmergable_ids = df.index.difference(mergable_ids)
-
-                merged = pd.concat(
-                    (merged, df.loc[nonmergable_ids]), names=df.index.names
-                )
-
-            search = lambda a, b: np.where(
-                np.in1d(
-                    np.ravel_multi_index(a.T, a.max(0) + 1),
-                    np.ravel_multi_index(b.T, a.max(0) + 1),
-                )
-            )
-            if "modifiers/picks" in f:
-                picks = self.get_picks(names=merged.index.names)
-                missing_cells = [i for i in picks if tuple(i) not in set(merged.index)]
-
-                if picks:
-                    # return merged.loc[
-                    #     set(picks).intersection([tuple(x) for x in merged.index])
-                    # ]
-                    return merged.loc[picks]
-                else:
-                    if isinstance(merged.index, pd.MultiIndex):
-                        empty_lvls = [[] for i in merged.index.names]
-                        index = pd.MultiIndex(
-                            levels=empty_lvls,
-                            codes=empty_lvls,
-                            names=merged.index.names,
-                        )
-                    else:
-                        index = pd.Index([], name=merged.index.name)
-                    merged = pd.DataFrame([], index=index)
-            return merged
-
-    @property
-    def datasets(self):
-        with h5py.File(self.filename, "r") as f:
-            dsets = f.visititems(self._if_ext_or_post)
-        return dsets
-
-    def get_merged(self, dataset):
-        return self.apply_prepost(dataset, skip_pick=True)
-
-    @property
-    def merges(self):
-        with h5py.File(self.filename, "r") as f:
-            dsets = f.visititems(self._if_merges)
-        return dsets
-
-    @property
-    def n_merges(self):
-        print("{} merge events".format(len(self.merges)))
-
-    @property
-    def merges(self):
-        with h5py.File(self.filename, "r") as f:
-            dsets = f.visititems(self._if_merges)
-        return dsets
-
-    @property
-    def picks(self):
-        with h5py.File(self.filename, "r") as f:
-            dsets = f.visititems(self._if_picks)
-        return dsets
-
-    def apply_merge(self, df, changes):
-        if len(changes):
-
-            for target, source in changes:
-                df.loc[tuple(target)] = self.join_tracks_pair(
-                    df.loc[tuple(target)], df.loc[tuple(source)]
-                )
-                df.drop(tuple(source), inplace=True)
-
-        return df
-
-    def get_raw(self, dataset):
-        if isinstance(dataset, str):
-            with h5py.File(self.filename, "r") as f:
-                return self.dset_to_df(f, dataset)
-        elif isinstance(dataset, list):
-            return [self.get_raw(dset) for dset in dataset]
-
-    def get_merges(self):
-        # fetch merge events going up to the first level
-        with h5py.File(self.filename, "r") as f:
-            merges = f.get("modifiers/merges", np.array([]))
-            if not isinstance(merges, np.ndarray):
-                merges = merges[()]
-
-        return merges
-
-    # def get_picks(self, levels):
-    def get_picks(self, names, path="modifiers/picks/"):
-        with h5py.File(self.filename, "r") as f:
-            if path in f:
-                return list(zip(*[f[path + name] for name in names]))
-                # return f["modifiers/picks"]
-            else:
-                return None
-
-    def dset_to_df(self, f, dataset):
-        dset = f[dataset]
-        names = copy(self.names)
-        if not dataset.endswith("imBackground"):
-            names.append("cell_label")
-        lbls = {lbl: dset[lbl][()] for lbl in names if lbl in dset.keys()}
-        index = pd.MultiIndex.from_arrays(
-            list(lbls.values()), names=names[-len(lbls) :]
-        )
-
-        columns = (
-            dset["timepoint"][()] if "timepoint" in dset else dset.attrs["columns"]
-        )
-
-        df = pd.DataFrame(dset[("values")][()], index=index, columns=columns)
-
-        return df
-
-    @staticmethod
-    def dataset_to_df(f: h5py.File, path: str, mode: str = "h5py"):
-
-        if mode is "h5py":
-            all_indices = ["experiment", "position", "trap", "cell_label"]
-            indices = {k: f[path][k][()] for k in all_indices if k in f[path].keys()}
-            return pd.DataFrame(
-                f[path + "/values"][()],
-                index=pd.MultiIndex.from_arrays(
-                    list(indices.values()), names=indices.keys()
-                ),
-                columns=f[path + "/timepoint"][()],
-            )
-
-    @staticmethod
-    def _if_ext_or_post(name, *args):
-        flag = False
-        if name.startswith("extraction") and len(name.split("/")) == 4:
-            flag = True
-        elif name.startswith("postprocessing") and len(name.split("/")) == 3:
-            flag = True
-
-        if flag:
-            print(name)
-
-    @staticmethod
-    def _if_merges(name: str, obj):
-        if isinstance(obj, h5py.Dataset) and name.startswith("modifiers/merges"):
-            return obj[()]
-
-    @staticmethod
-    def _if_picks(name: str, obj):
-        if isinstance(obj, h5py.Group) and name.endswith("picks"):
-            return obj[()]
-
-    @staticmethod
-    def join_tracks_pair(target, source):
-        tgt_copy = copy(target)
-        end = find_1st(target.values[::-1], 0, cmp_larger)
-        tgt_copy.iloc[-end:] = source.iloc[-end:].values
-        return tgt_copy

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"