diff --git a/brainbox/CONTRIBUTING.md b/brainbox/CONTRIBUTING.md
index 2c293c455..03560beea 100644
--- a/brainbox/CONTRIBUTING.md
+++ b/brainbox/CONTRIBUTING.md
@@ -24,7 +24,9 @@ We suggest using [Anaconda](https://www.anaconda.com/distribution/), which is de
 Once you have installed Anaconda, the next step is to create an environment for working with brainbox. This requires you to have the `environment.yml` file which lives in the top directory of this repository. We will just clone the whole repository now though, since you will need it later, using the following command on *nix systems:
 
 ```bash
-git clone https://github.com/int-brain-lab/brainbox
+git clone https://github.com/int-brain-lab/ibllib/
+cd ./ibllib
+git checkout brainbox
 ```
 
 Note: please navigate to the folder where you want to run this command beforehand, e.g. `/home/username/Documents` if you want the `brainbox` repository to live in your Documents folder
@@ -47,7 +49,7 @@ if you are developing from the terminal, in order to activate the environment yo
 
 # Git, GitFlow, and you (15 minutes)
 
-**TL;DR: We use [Git](https://rogerdudler.github.io/git-guide/) with a [GitFlow](https://www.atlassian.com/git/tutorials/comparing-workflows/gitflow-workflow) workflow to develop Brainbox. Please create new feature branches of `develop` for writing code and then make a pull request to have it added to `develop`.**
+**TL;DR: We use [Git](https://rogerdudler.github.io/git-guide/) with a [GitFlow](https://www.atlassian.com/git/tutorials/comparing-workflows/gitflow-workflow) workflow to develop Brainbox. Please create new feature branches of `brainbox` for writing code and then make a pull request to have it added to `brainbox`.**
 
 For those unfamiliar with it, Git is a system for *version control*, which allows you to make changes to whatever you put into it (Git isn't limited to just code!) that are:
 
@@ -64,17 +66,11 @@ This way you can keep a version of your code that you know works (called `master
 
 For an explanation of the basics of Git, [this guide by Roger Dudler](http://git.huit.harvard.edu/guide/) is a necessary five-minute read on the basics of manipulating a repository.
 
-Brainbox uses [GitFlow](https://www.atlassian.com/git/tutorials/comparing-workflows/gitflow-workflow) as a model for how to organize our repository. This means that there are two main branches that always exist, `master` and `develop`, the latter of which is the basis for all development of the toolbox. If you want to incorporate a new feature into the repository, e.g. a raster plot, you can run the following git command in your repository:
-
-```bash
-git flow feature start rasterplot
-```
-
-and Git will automatically create a new branch for you to work on and make it active. Once you've created a few commits and feel confident that your code is working well, you can create a pull request on GitHub so we can incorporate your code into the `develop` branch for future release.
+Brainbox uses [GitFlow](https://www.atlassian.com/git/tutorials/comparing-workflows/gitflow-workflow) as a model for how to organize our repository. This means that there are two main branches that always exist, `master` and `develop`, the latter of which we have renamed `brainbox`, as brainbox is a part of `ibllib`. `brainbox` is the basis for all development of the toolbox.
 
 
 # Writing code for Brainbox
 
-We require all code in Brainbox to conform to [PEP8](https://www.python.org/dev/peps/pep-0008/) guidelines, with [Numpy-style](https://numpydoc.readthedocs.io/en/latest/format.html) docstrings. We require all contributors to use `flake8` as a linter to check their code before a pull request.
+We require all code in Brainbox to conform to [PEP8](https://www.python.org/dev/peps/pep-0008/) guidelines, with [Numpy-style](https://numpydoc.readthedocs.io/en/latest/format.html) docstrings. We require all contributors to use `flake8` as a linter to check their code before a pull request. Please check the `.flake8` file in the top level of `ibllib` for an exact specification for how to set up your particular instance of flake8.
 
 [MORE GUIDELINES HERE PLEASE]
\ No newline at end of file
diff --git a/brainbox/core/core.py b/brainbox/core/core.py
index a199e7bb4..04b127d51 100644
--- a/brainbox/core/core.py
+++ b/brainbox/core/core.py
@@ -24,7 +24,8 @@ def __init__(self, times, values, columns=None, *args, **kwargs):
         a time stamp associated. TS objects have obligatory 'times' and 'values' entries which
         must be passed at construction, the length of both of which must match. TimeSeries takes an
         optional 'columns' argument, which defaults to None, that is a set of labels for the
-        columns in 'values'.
+        columns in 'values'. These are also exposed via the dot syntax as pointers to the specific
+        columns which they reference.
 
         :param times: an ordered object containing a list of timestamps for the time series data
         :param values: an ordered object containing the associated measurements for each time stamp
diff --git a/brainbox/environment.yml b/brainbox/environment.yml
new file mode 100644
index 000000000..26b8aba8d
--- /dev/null
+++ b/brainbox/environment.yml
@@ -0,0 +1,190 @@
+name: brainboxdev
+channels:
+  - defaults
+dependencies:
+  - _libgcc_mutex=0.1=main
+  - _tflow_select=2.3.0=mkl
+  - absl-py=0.7.1=py37_0
+  - alabaster=0.7.12=py37_0
+  - asn1crypto=0.24.0=py37_0
+  - astor=0.7.1=py37_0
+  - astroid=2.2.5=py37_0
+  - attrs=19.1.0=py37_1
+  - babel=2.7.0=py_0
+  - backcall=0.1.0=py37_0
+  - blas=1.0=mkl
+  - bleach=3.1.0=py37_0
+  - bzip2=1.0.8=h7b6447c_0
+  - c-ares=1.15.0=h7b6447c_1
+  - ca-certificates=2019.5.15=0
+  - certifi=2019.6.16=py37_0
+  - cffi=1.12.2=py37h2e261b9_1
+  - chardet=3.0.4=py37_1
+  - cloudpickle=1.2.1=py_0
+  - conda-package-handling=1.3.11=py37_0
+  - cryptography=2.6.1=py37h1ba5d50_0
+  - cycler=0.10.0=py37_0
+  - dbus=1.13.6=h746ee38_0
+  - decorator=4.4.0=py37_1
+  - defusedxml=0.6.0=py_0
+  - docutils=0.14=py37_0
+  - entrypoints=0.3=py37_0
+  - expat=2.2.6=he6710b0_0
+  - fontconfig=2.13.0=h9420a91_0
+  - freetype=2.9.1=h8a8886c_1
+  - gast=0.2.2=py37_0
+  - glib=2.56.2=hd408876_0
+  - gmp=6.1.2=h6c8ec71_1
+  - google-pasta=0.1.7=py_0
+  - grpcio=1.16.1=py37hf8bcb03_1
+  - gst-plugins-base=1.14.0=hbbd80ab_1
+  - gstreamer=1.14.0=hb453b48_1
+  - h5py=2.9.0=py37h7918eee_0
+  - hdf5=1.10.4=hb1b8bf9_0
+  - icu=58.2=h9c2bf20_1
+  - idna=2.8=py37_0
+  - imagesize=1.1.0=py37_0
+  - intel-openmp=2019.4=243
+  - ipykernel=5.1.1=py37h39e3cac_0
+  - ipython=7.6.1=py37h39e3cac_0
+  - ipython_genutils=0.2.0=py37_0
+  - ipywidgets=7.5.0=py_0
+  - isort=4.3.21=py37_0
+  - jedi=0.13.3=py37_0
+  - jeepney=0.4=py37_0
+  - jinja2=2.10.1=py37_0
+  - joblib=0.13.2=py37_0
+  - jpeg=9b=h024ee3a_2
+  - jsonschema=3.0.1=py37_0
+  - jupyter=1.0.0=py37_7
+  - jupyter_client=5.3.1=py_0
+  - jupyter_console=6.0.0=py37_0
+  - jupyter_core=4.5.0=py_0
+  - keras-applications=1.0.8=py_0
+  - keras-preprocessing=1.1.0=py_1
+  - keyring=18.0.0=py37_0
+  - kiwisolver=1.1.0=py37he6710b0_0
+  - lazy-object-proxy=1.4.1=py37h7b6447c_0
+  - libarchive=3.3.3=h5d8350f_5
+  - libedit=3.1.20181209=hc058e9b_0
+  - libffi=3.2.1=hd88cf55_4
+  - libgcc-ng=8.2.0=hdf63c60_1
+  - libgfortran-ng=7.3.0=hdf63c60_0
+  - libpng=1.6.37=hbc83047_0
+  - libprotobuf=3.8.0=hd408876_0
+  - libsodium=1.0.16=h1bed415_0
+  - libstdcxx-ng=8.2.0=hdf63c60_1
+  - libuuid=1.0.3=h1bed415_2
+  - libxcb=1.13=h1bed415_1
+  - libxml2=2.9.9=hea5a465_1
+  - lz4-c=1.8.1.2=h14c3975_0
+  - lzo=2.10=h49e0be7_2
+  - markdown=3.1.1=py37_0
+  - markupsafe=1.1.1=py37h7b6447c_0
+  - matplotlib=3.1.0=py37h5429711_0
+  - mccabe=0.6.1=py37_1
+  - mistune=0.8.4=py37h7b6447c_0
+  - mkl=2019.4=243
+  - mkl-service=2.0.2=py37h7b6447c_0
+  - mkl_fft=1.0.12=py37ha843d7b_0
+  - mkl_random=1.0.2=py37hd81dba3_0
+  - nbconvert=5.5.0=py_0
+  - nbformat=4.4.0=py37_0
+  - ncurses=6.1=he6710b0_1
+  - notebook=5.7.8=py37_0
+  - numpy=1.16.4=py37h7e9f1db_0
+  - numpy-base=1.16.4=py37hde5b4d6_0
+  - numpydoc=0.9.1=py_0
+  - openssl=1.1.1c=h7b6447c_1
+  - packaging=19.0=py37_0
+  - pandas=0.24.2=py37he6710b0_0
+  - pandoc=2.2.3.2=0
+  - pandocfilters=1.4.2=py37_1
+  - parso=0.5.0=py_0
+  - patsy=0.5.1=py37_0
+  - pcre=8.43=he6710b0_0
+  - pexpect=4.7.0=py37_0
+  - pickleshare=0.7.5=py37_0
+  - pip=19.0.3=py37_0
+  - prometheus_client=0.7.1=py_0
+  - prompt_toolkit=2.0.9=py37_0
+  - protobuf=3.8.0=py37he6710b0_0
+  - psutil=5.6.3=py37h7b6447c_0
+  - ptyprocess=0.6.0=py37_0
+  - pycodestyle=2.5.0=py37_0
+  - pycosat=0.6.3=py37h14c3975_0
+  - pycparser=2.19=py37_0
+  - pyflakes=2.1.1=py37_0
+  - pygments=2.4.2=py_0
+  - pylint=2.3.1=py37_0
+  - pyopenssl=19.0.0=py37_0
+  - pyparsing=2.4.0=py_0
+  - pyqt=5.9.2=py37h05f1152_2
+  - pyrsistent=0.14.11=py37h7b6447c_0
+  - pysocks=1.6.8=py37_0
+  - python=3.7.3=h0371630_0
+  - python-dateutil=2.8.0=py37_0
+  - python-libarchive-c=2.8=py37_10
+  - pytz=2019.1=py_0
+  - pyzmq=18.0.0=py37he6710b0_0
+  - qt=5.9.7=h5867ecd_1
+  - qtawesome=0.5.7=py37_1
+  - qtconsole=4.5.1=py_0
+  - qtpy=1.8.0=py_0
+  - readline=7.0=h7b6447c_5
+  - requests=2.22.0=py37_0
+  - rope=0.14.0=py_0
+  - ruamel_yaml=0.15.46=py37h14c3975_0
+  - scikit-learn=0.21.2=py37hd81dba3_0
+  - seaborn=0.9.0=py37_0
+  - secretstorage=3.1.1=py37_0
+  - send2trash=1.5.0=py37_0
+  - setuptools=41.0.0=py37_0
+  - sip=4.19.8=py37hf484d3e_0
+  - six=1.12.0=py37_0
+  - snowballstemmer=1.9.0=py_0
+  - sphinx=2.1.2=py_0
+  - sphinxcontrib-applehelp=1.0.1=py_0
+  - sphinxcontrib-devhelp=1.0.1=py_0
+  - sphinxcontrib-htmlhelp=1.0.2=py_0
+  - sphinxcontrib-jsmath=1.0.1=py_0
+  - sphinxcontrib-qthelp=1.0.2=py_0
+  - sphinxcontrib-serializinghtml=1.1.3=py_0
+  - spyder=3.3.6=py37_0
+  - spyder-kernels=0.5.1=py37_0
+  - sqlite=3.27.2=h7b6447c_0
+  - statsmodels=0.10.0=py37hdd07704_0
+  - tensorboard=1.14.0=py37hf484d3e_0
+  - tensorflow=1.14.0=mkl_py37h45c423b_0
+  - tensorflow-base=1.14.0=mkl_py37h7ce6ba3_0
+  - tensorflow-estimator=1.14.0=py_0
+  - termcolor=1.1.0=py37_1
+  - terminado=0.8.2=py37_0
+  - testpath=0.4.2=py37_0
+  - tk=8.6.8=hbc83047_0
+  - tornado=6.0.3=py37h7b6447c_0
+  - tqdm=4.32.1=py_0
+  - traitlets=4.3.2=py37_0
+  - urllib3=1.24.1=py37_0
+  - wcwidth=0.1.7=py37_0
+  - webencodings=0.5.1=py37_1
+  - werkzeug=0.15.4=py_0
+  - wheel=0.33.1=py37_0
+  - widgetsnbextension=3.5.0=py37_0
+  - wrapt=1.11.2=py37h7b6447c_0
+  - wurlitzer=1.0.2=py37_0
+  - xz=5.2.4=h14c3975_4
+  - yaml=0.1.7=had09818_2
+  - zeromq=4.3.1=he6710b0_3
+  - zlib=1.2.11=h7b6447c_3
+  - zstd=1.3.7=h0b5b093_0
+  - pip:
+    - autopep8==1.4.4
+    - colorlog==4.0.2
+    - dataclasses==0.6
+    - flake8==3.7.8
+    - globus-sdk==1.8.0
+    - ibllib==1.0.6
+    - pyjwt==1.7.1
+    - scipy==1.3.0
+
diff --git a/brainbox/plot/__init__.py b/brainbox/plot/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/brainbox/population/__init__.py b/brainbox/population/__init__.py
index e69de29bb..241d62821 100644
--- a/brainbox/population/__init__.py
+++ b/brainbox/population/__init__.py
@@ -0,0 +1 @@
+from .population import *
\ No newline at end of file
diff --git a/brainbox/population/population.py b/brainbox/population/population.py
new file mode 100644
index 000000000..8cb919f6c
--- /dev/null
+++ b/brainbox/population/population.py
@@ -0,0 +1,141 @@
+'''
+Population functions.
+
+Code from https://github.com/cortex-lab/phylib/blob/master/phylib/stats/ccg.py by C. Rossant.
+'''
+
+import numpy as np
+
+
+def _index_of(arr, lookup):
+    """Replace scalars in an array by their indices in a lookup table.
+
+    Implicitely assume that:
+
+    * All elements of arr and lookup are non-negative integers.
+    * All elements or arr belong to lookup.
+
+    This is not checked for performance reasons.
+
+    """
+    # Equivalent of np.digitize(arr, lookup) - 1, but much faster.
+    # TODO: assertions to disable in production for performance reasons.
+    # TODO: np.searchsorted(lookup, arr) is faster on small arrays with large
+    # values
+    lookup = np.asarray(lookup, dtype=np.int32)
+    m = (lookup.max() if len(lookup) else 0) + 1
+    tmp = np.zeros(m + 1, dtype=np.int)
+    # Ensure that -1 values are kept.
+    tmp[-1] = -1
+    if len(lookup):
+        tmp[lookup] = np.arange(len(lookup))
+    return tmp[arr]
+
+
+def _increment(arr, indices):
+    """Increment some indices in a 1D vector of non-negative integers.
+    Repeated indices are taken into account."""
+    bbins = np.bincount(indices)
+    arr[:len(bbins)] += bbins
+    return arr
+
+
+def _diff_shifted(arr, steps=1):
+    return arr[steps:] - arr[:len(arr) - steps]
+
+
+def _create_correlograms_array(n_clusters, winsize_bins):
+    return np.zeros((n_clusters, n_clusters, winsize_bins // 2 + 1), dtype=np.int32)
+
+
+def _symmetrize_correlograms(correlograms):
+    """Return the symmetrized version of the CCG arrays."""
+
+    n_clusters, _, n_bins = correlograms.shape
+    assert n_clusters == _
+
+    # We symmetrize c[i, j, 0].
+    # This is necessary because the algorithm in correlograms()
+    # is sensitive to the order of identical spikes.
+    correlograms[..., 0] = np.maximum(
+        correlograms[..., 0], correlograms[..., 0].T)
+
+    sym = correlograms[..., 1:][..., ::-1]
+    sym = np.transpose(sym, (1, 0, 2))
+
+    return np.dstack((sym, correlograms))
+
+
+def xcorr(spike_times, spike_clusters, bin_size=None, window_size=None):
+    """Compute all pairwise cross-correlograms among the clusters appearing in `spike_clusters`.
+
+    Parameters
+    ----------
+
+    :param spike_times: Spike times in seconds.
+    :type spike_times: array-like
+    :param spike_clusters: Spike-cluster mapping.
+    :type spike_clusters: array-like
+    :param bin_size: Size of the bin, in seconds.
+    :type bin_size: float
+    :param window_size: Size of the window, in seconds.
+    :type window_size: float
+
+    Returns an `(n_clusters, n_clusters, winsize_samples)` array with all pairwise
+    cross-correlograms.
+
+    """
+    assert np.all(np.diff(spike_times) >= 0), ("The spike times must be increasing.")
+    assert spike_times.ndim == 1
+    assert spike_times.shape == spike_clusters.shape
+
+    # Find `binsize`.
+    bin_size = np.clip(bin_size, 1e-5, 1e5)  # in seconds
+
+    # Find `winsize_bins`.
+    window_size = np.clip(window_size, 1e-5, 1e5)  # in seconds
+    winsize_bins = 2 * int(.5 * window_size / bin_size) + 1
+
+    # Take the cluster order into account.
+    clusters = np.unique(spike_clusters)
+    n_clusters = len(clusters)
+
+    # Like spike_clusters, but with 0..n_clusters-1 indices.
+    spike_clusters_i = _index_of(spike_clusters, clusters)
+
+    # Shift between the two copies of the spike trains.
+    shift = 1
+
+    # At a given shift, the mask precises which spikes have matching spikes
+    # within the correlogram time window.
+    mask = np.ones_like(spike_times, dtype=np.bool)
+
+    correlograms = _create_correlograms_array(n_clusters, winsize_bins)
+
+    # The loop continues as long as there is at least one spike with
+    # a matching spike.
+    while mask[:-shift].any():
+        # Interval between spike i and spike i+shift.
+        spike_diff = _diff_shifted(spike_times, shift)
+
+        # Binarize the delays between spike i and spike i+shift.
+        spike_diff_b = np.round(spike_diff / bin_size).astype(np.int64)
+
+        # Spikes with no matching spikes are masked.
+        mask[:-shift][spike_diff_b > (winsize_bins / 2)] = False
+
+        # Cache the masked spike delays.
+        m = mask[:-shift].copy()
+        d = spike_diff_b[m]
+
+        # Find the indices in the raveled correlograms array that need
+        # to be incremented, taking into account the spike clusters.
+        indices = np.ravel_multi_index(
+            (spike_clusters_i[:-shift][m], spike_clusters_i[+shift:][m], d), correlograms.shape)
+
+        # Increment the matching spikes in the correlograms array.
+        _increment(correlograms.ravel(), indices)
+
+        shift += 1
+
+    return _symmetrize_correlograms(correlograms)
diff --git a/brainbox/processing/processing.py b/brainbox/processing/processing.py
index 81e16e0aa..4e093a41c 100644
--- a/brainbox/processing/processing.py
+++ b/brainbox/processing/processing.py
@@ -9,7 +9,7 @@
 from brainbox import core
 
 
-def sync(dt, timeseries=None, times=None, values=None, offsets=None, interp='zero',
+def sync(dt, times=None, values=None, timeseries=None, offsets=None, interp='zero',
          fillval=np.nan):
     """
     Function for resampling a single or multiple time series to a single, evenly-spaced, delta t
diff --git a/brainbox/singlecell/__init__.py b/brainbox/singlecell/__init__.py
index e69de29bb..4d5935db6 100644
--- a/brainbox/singlecell/__init__.py
+++ b/brainbox/singlecell/__init__.py
@@ -0,0 +1 @@
+from .singlecell import *
diff --git a/brainbox/singlecell/singlecell.py b/brainbox/singlecell/singlecell.py
new file mode 100644
index 000000000..f66880369
--- /dev/null
+++ b/brainbox/singlecell/singlecell.py
@@ -0,0 +1,26 @@
+'''
+Single-cell functions.
+'''
+
+import numpy as np
+from brainbox.population import xcorr
+
+
+def acorr(spike_times, bin_size=None, window_size=None):
+    """Compute the auto-correlogram of a neuron.
+
+    Parameters
+    ----------
+
+    :param spike_times: Spike times in seconds.
+    :type spike_times: array-like
+    :param bin_size: Size of the bin, in seconds.
+    :type bin_size: float
+    :param window_size: Size of the window, in seconds.
+    :type window_size: float
+
+    Returns an `(winsize_samples,)` array with the auto-correlogram.
+
+    """
+    xc = xcorr(spike_times, np.zeros_like(spike_times), bin_size=bin_size, window_size=window_size)
+    return xc[0, 0, :]
diff --git a/env_ibl.yml b/env_ibl.yml
new file mode 100644
index 000000000..8025667f0
--- /dev/null
+++ b/env_ibl.yml
@@ -0,0 +1,16 @@
+name: ibl
+channels:
+  - defaults
+dependencies:
+  - numpy
+  - scipy
+  - seaborn
+  - flake8
+  - matplotlib
+  - pandas
+  - requests
+  - pip:
+    - colorlog>=4.0.2
+    - dataclasses>=0.6
+    - globus-sdk>=1.8.0
+
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 000000000..28b3de2ab
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,21 @@
+name: ibldev
+channels:
+  - defaults
+dependencies:
+  - numpy
+  - scipy
+  - seaborn
+  - flake8
+  - matplotlib
+  - pandas
+  - requests
+  - jupyter
+  - tensorflow
+  - keras
+  - scikit-learn
+  - ipython
+  - plotly
+  - pip:
+    - colorlog>=4.0.2
+    - dataclasses>=0.6
+    - globus-sdk>=1.8.0
diff --git a/examples/brainbox/plot_peths.py b/examples/brainbox/plot_peths.py
new file mode 100644
index 000000000..5208375b2
--- /dev/null
+++ b/examples/brainbox/plot_peths.py
@@ -0,0 +1,208 @@
+import numpy as np
+
+from pathlib import Path
+from oneibl.one import ONE
+import alf.io as ioalf
+
+
+def filter_trials(trials, choice, stim_side, stim_contrast):
+    """
+    Filter trials by choice and presented stimulus constrast/side
+
+    :param trials: dict with keys `choice`, `contrastLeft`, `contrastRight`
+    :type trials: dict
+    :param choice: subject's choice (1=left, -1=right)
+    :type choice: int
+    :param stim_side: where stimulus was presented (1=left, -1=right)
+    :type stim_side: int
+    :param stim_contrast: contrast of grating stimulus
+    :type stim_contrast: float
+    :return: subset of trials filtered by choice and stimulus contrast/side
+    :rtype: np.ndarray
+    """
+    contrast = 'contrastLeft' if stim_side == 1 else 'contrastRight'
+    trial_ids = np.where(
+        (trials['choice'] == choice) & (trials[contrast] == stim_contrast))[0]
+    return trial_ids
+
+
+def calculate_peths(
+        spike_times, spike_clusters, cluster_ids, align_times, pre_time=0.2,
+        post_time=0.5, bin_size=0.025, smoothing=0.025):
+    """
+    Calcluate peri-event time histograms; return means and standard deviations
+    for each time point across specified clusters
+
+    :param spike_times: spike times (in seconds)
+    :type spike_times: array-like
+    :param spike_clusters: cluster ids corresponding to each event in `spikes`
+    :type spike_clusters: array-like
+    :param cluster_ids: subset of cluster ids for calculating peths
+    :type cluster_ids: array-like
+    :param align_times: times (in seconds) to align peths to
+    :type align_times: array-like
+    :param pre_time: time (in seconds) to precede align times in peth
+    :type pre_time: float
+    :param post_time: time (in seconds) to follow align times in peth
+    :type post_time: float
+    :param bin_size: width of time windows (in seconds) to bin spikes
+    :type bin_size: float
+    :param smoothing: standard deviation (in seconds) of Gaussian kernel for
+        smoothing peths; use `smoothing=0` to skip smoothing
+    :type smoothing: float
+    :return: (psth_means, psth_stds)
+    :rtype: tuple with two elements, each of shape `(n_trials, n_clusters, n_bins)`
+    """
+
+    from scipy.signal import gaussian
+    from scipy.signal import convolve
+
+    # initialize containers
+    n_bins_pre = int(np.ceil(pre_time / bin_size))
+    n_bins_post = int(np.ceil(post_time / bin_size))
+    n_bins = n_bins_pre + n_bins_post
+    binned_spikes = np.zeros(shape=(len(align_times), len(cluster_ids), n_bins))
+
+    # build gaussian kernel if requested
+    if smoothing > 0:
+        w = n_bins - 1 if n_bins % 2 == 0 else n_bins
+        window = gaussian(w, std=smoothing / bin_size)
+
+    ids = np.unique(cluster_ids)
+
+    # bin spikes
+    for i, t_0 in enumerate(align_times):
+
+        # define bin edges
+        ts_pre = t_0 - np.arange(n_bins_pre, 0, -1) * bin_size
+        ts_post = t_0 + np.arange(n_bins_post + 1) * bin_size
+        ts = np.concatenate([ts_pre, ts_post])
+
+        # filter spikes
+        idxs = (spike_times > ts[0]) & \
+               (spike_times <= ts[-1]) & np.isin(spike_clusters, cluster_ids)
+        i_spikes = spike_times[idxs]
+        i_clusters = spike_clusters[idxs]
+
+        # bin spikes similar to bincount2D: x = spike times, y = spike clusters
+        xscale = ts
+        xind = (np.floor((i_spikes - ts[0]) / bin_size)).astype(np.int64)
+        yscale, yind = np.unique(i_clusters, return_inverse=True)
+        nx, ny = [xscale.size, yscale.size]
+        ind2d = np.ravel_multi_index(np.c_[yind, xind].transpose(), dims=(ny, nx))
+        r = np.bincount(ind2d, minlength=nx * ny, weights=None).reshape(ny, nx)
+
+        # smooth
+        if smoothing > 0:
+            for j in range(r.shape[0]):
+                r[j, :] = convolve(r[j, :], window, mode='same', method='auto')
+
+        # store (ts represent bin edges, so there are one fewer bins)
+        bs_idxs = np.isin(ids, yscale)
+        binned_spikes[i, bs_idxs, :] = r[:, :-1]
+
+    # average
+    peth_means = np.mean(binned_spikes, axis=0)
+    peth_stds = np.std(binned_spikes, axis=0)
+
+    return peth_means, peth_stds
+
+
+def plot_peths(
+        means, stds=None, marker_idx=0, bin_size=1.0, linewidth=1,
+        onset_label='event', **kwargs):
+    """
+    Plot peths with optional errorbars
+
+    :param means: peth means
+    :type means: np.ndarray of shape `(n_clusters, n_bins)`
+    :param stds: peth standard deviations
+    :type stds: np.ndarray of shape `(n_clusters, n_bins)`
+    :param marker_idx: index into bin dimension for drawing dashed line
+    :type marker_idx: int
+    :param bin_size: bin size of peth (in seconds) for labeling time axis
+    :type bin_size: float
+    :param linewidth: linewidth of peth means
+    :type linewidth: float
+    :param onset_label: name of alignment type, i.e. 'go cue', 'stimulus', etc.
+        the x-axis is labeled as 'Time from `onset_label` onset (s)'
+    :type onset_label: str
+    :param kwargs: additional arguments into matplotlib.axis.Axes.plot function
+    :return: figure handle
+    :rtype: matplotlib figure handle
+    """
+
+    import matplotlib.pyplot as plt
+    from ibllib.dsp import rms
+
+    scale = 1.0 / rms(means.flatten())
+    n_clusts, n_bins = means.shape
+    ts = (np.arange(n_bins) - marker_idx) * bin_size
+
+    fig = plt.figure(figsize=(3, 0.5 * n_clusts))
+    ax = plt.gca()
+    ax.spines['top'].set_visible(False)
+    ax.spines['left'].set_visible(False)
+    ax.spines['right'].set_visible(False)
+
+    # plot individual traces
+    for i in range(n_clusts):
+        # extract mean and offset
+        ys = i + means[i, :] * scale
+        # plot mean
+        ax.plot(ts, ys, linewidth=linewidth, **kwargs)
+        # plot standard error
+        if stds is not None:
+            ax.fill_between(
+                ts, ys - stds[i, :] * scale, ys + stds[i, :] * scale,
+                linewidth=linewidth, alpha=0.5)
+    # plot vertical line at marker_idx (time = 0 seconds)
+    ax.axvline(x=0, ymin=0.02, ymax=0.98, linestyle='--', color='k')
+    # add labels
+    ax.set_xlabel('Time from %s onset (s)' % onset_label, fontsize=12)
+    ax.set_ylim([-0.5, n_clusts - 0.5])
+    ax.set_yticks([])
+    ax.set_yticklabels([])
+
+    plt.show()
+
+    return fig
+
+
+if __name__ == '__main__':
+
+    BIN_SIZE = 0.025  # seconds
+    SMOOTH_SIZE = 0.025  # seconds; standard deviation of gaussian kernel
+    PRE_TIME = 0.2  # seconds to plot before event onset
+    POST_TIME = 0.5  # seconds to plot after event onset
+
+    # get the data from flatiron
+    one = ONE()
+    eid = one.search(subject='ZM_1735', date='2019-08-01', number=1)
+    D = one.load(eid[0], clobber=False, download_only=True)
+    session_path = Path(D.local_path[0]).parent
+
+    # load objects
+    spikes = ioalf.load_object(session_path, 'spikes')
+    trials = ioalf.load_object(session_path, '_ibl_trials')
+
+    # filter trials by choice and contrast/side
+    trial_ids = filter_trials(trials, choice=1, stim_contrast=1.0, stim_side=1)
+
+    # align to go cue
+    align_times = trials['goCue_times'][trial_ids]
+
+    # define subset of clusters to plot
+    cluster_ids = np.unique(spikes['clusters'])[np.arange(10)]
+
+    # calculate peths of specified trials/clusters, aligned to desired event
+    peth_means, peth_stds = calculate_peths(
+        spikes['times'], spikes['clusters'], cluster_ids, align_times,
+        pre_time=PRE_TIME, post_time=POST_TIME, bin_size=BIN_SIZE,
+        smoothing=SMOOTH_SIZE)
+
+    # plot peths along with standard error
+    fig = plot_peths(
+        peth_means, peth_stds / np.sqrt(len(trial_ids)),
+        marker_idx=int(PRE_TIME / BIN_SIZE), bin_size=BIN_SIZE,
+        onset_label='go cue')
diff --git a/examples/brainbox/raster_cluster_ordered.py b/examples/brainbox/raster_cluster_ordered.py
new file mode 100644
index 000000000..d5d845659
--- /dev/null
+++ b/examples/brainbox/raster_cluster_ordered.py
@@ -0,0 +1,50 @@
+from pathlib import Path
+import matplotlib.pyplot as plt
+import numpy as np
+import rastermap
+from oneibl.one import ONE
+import alf.io as ioalf
+import ibllib.plots as iblplt
+from brainbox.processing import bincount2D
+
+T_BIN = 0.01
+
+# get the data from flatiron and the current folder
+one = ONE()
+eid = one.search(lab='wittenlab', date='2019-08-04')
+D = one.load(eid[0])
+session_path = Path(D.local_path[0]).parent
+
+# load objects
+spikes = ioalf.load_object(session_path, 'spikes')
+clusters = ioalf.load_object(session_path, 'clusters')
+channels = ioalf.load_object(session_path, 'channels')
+trials = ioalf.load_object(session_path, '_ibl_trials')
+
+# compute raster map as a function of cluster number
+R, times, clusters = bincount2D(spikes['times'], spikes['clusters'], T_BIN)
+
+
+# Using rastermap defaults to order activity matrix
+# by similarity of activity (requires R to contain floats)
+model = rastermap.mapping.Rastermap().fit(R.astype(float))
+isort = np.argsort(model.embedding[:, 0])
+R = R[isort, :]
+
+# Alternatively, order activity by cortical depth of neurons
+# d=dict(zip(spikes['clusters'],spikes['depths']))
+# y=sorted([[i,d[i]] for i in d])
+# isort=argsort([x[1] for x in y])
+# R=R[isort,:]
+
+# plot raster map
+plt.imshow(R, aspect='auto', cmap='binary', vmax=T_BIN / 0.001 / 4,
+           extent=np.r_[times[[0, -1]], clusters[[0, -1]]], origin='lower')
+# plot trial start and reward time
+reward = trials['goCue_times']
+
+iblplt.vertical_lines(reward, ymin=0, ymax=clusters[-1], color='m', linewidth=0.5,
+                      label='valve openings')
+plt.xlabel('Time (s)')
+plt.ylabel('Cluster #')
+plt.legend()
diff --git a/examples/brainbox/raster_clusters.py b/examples/brainbox/raster_clusters.py
index 006192fb1..996068800 100644
--- a/examples/brainbox/raster_clusters.py
+++ b/examples/brainbox/raster_clusters.py
@@ -6,7 +6,7 @@
 import alf.io as ioalf
 import ibllib.plots as iblplt
 
-from brainbox.misc import bincount2D
+from brainbox.processing import bincount2D
 
 T_BIN = 0.01
 
diff --git a/ibllib_env.yml b/ibllib_env.yml
index a4e2763e0..26d13b60f 100644
--- a/ibllib_env.yml
+++ b/ibllib_env.yml
@@ -4,68 +4,21 @@ channels:
   - conda-forge
   - conda-canary
 dependencies:
-  - asn1crypto=0.24.0=py37_1003
-  - blas=1.0=mkl
-  - ca-certificates=2019.6.16=hecc5488_0
-  - certifi=2019.6.16=py37_1
-  - cffi=1.12.3=py37hb32ad35_0
-  - chardet=3.0.4=py37_1003
-  - colorama=0.4.1=py_0
-  - colorlog=4.0.2=py37_1000
-  - cryptography=2.7=py37hb32ad35_0
-  - cycler=0.10.0=py_1
-  - entrypoints=0.3=py37_1000
-  - flake8=3.7.8=py37_1
-  - freetype=2.10.0=h563cfd7_1
-  - globus-sdk=1.8.0=py_0
-  - icc_rt=2019.0.0=h0cc432a_1
-  - icu=64.2=he025d50_0
-  - idna=2.8=py37_1000
-  - intel-openmp=2019.4=245
-  - jpeg=9c=hfa6e2cd_1001
-  - kiwisolver=1.1.0=py37he980bc4_0
-  - libblas=3.8.0=12_mkl
-  - libcblas=3.8.0=12_mkl
-  - liblapack=3.8.0=12_mkl
-  - libpng=1.6.37=h7602738_0
-  - matplotlib=3.1.1=py37_1
-  - matplotlib-base=3.1.1=py37h2852a4a_1
-  - mccabe=0.6.1=py_1
-  - mkl=2019.4=245
-  - mkl-service=2.1.0=py37hfa6e2cd_0
-  - numpy=1.17.0=py37hc71023c_0
-  - openssl=1.1.1c=hfa6e2cd_0
-  - pandas=0.25.0=py37he350917_0
-  - patsy=0.5.1=py_0
-  - pip:
-    - colorlog>=4.0.2
-    - dataclasses>=0.6
-    - globus-sdk>=1.8.0 
-  - pycodestyle=2.5.0=py_0
-  - pycparser=2.19=py37_1
-  - pyflakes=2.1.1=py_0
-  - pyjwt=1.7.1=py_0
-  - pyopenssl=19.0.0=py37_0
-  - pyparsing=2.4.2=py_0
-  - pyqt=5.9.2=py37h6538335_2
-  - pysocks=1.7.0=py37_0
-  - python=3.7.3=h510b542_1
-  - python-dateutil=2.8.0=py_0
-  - pytz=2019.2=py_0
-  - qt=5.9.7=h506e8af_3
-  - requests=2.22.0=py37_1
-  - scipy=1.3.1=py37h29ff71c_0
-  - seaborn=0.9.0=py_1
-  - setuptools=41.0.1=py37_0
-  - sip=4.19.8=py37h6538335_1000
-  - six=1.12.0=py37_1000
-  - sqlite=3.29.0=hfa6e2cd_0
-  - statsmodels=0.10.1=py37hfa6e2cd_0
-  - tornado=6.0.3=py37hfa6e2cd_0
-  - urllib3=1.25.3=py37_0
-  - vc=14.1=h0510ff6_4
-  - vs2015_runtime=14.15.26706=h3a45250_4
-  - wheel=0.33.6=py37_0
-  - win_inet_pton=1.1.0=py37_0
-  - wincertstore=0.2=py37_1002
-  - zlib=1.2.11=h2fa13f4_1005
+ - numpy >=1.16.4
+ - scipy >=1.3.0
+ - seaborn >=0.9.0
+ - flake8 >=3.7.8
+ - matplotlib >=3.0.3
+ - pandas >=0.24.2
+ - requests >=2.22.0
+ - jupyter
+ - tensorflow
+ - keras
+ - scikit-learn
+ - ipython
+ - plotly
+ - pip >=19.0
+ - pip:
+   - colorlog >=4.0.2
+   - dataclasses >=0.6
+   - globus-sdk >=1.8.0
diff --git a/tests/brainbox/test_population.py b/tests/brainbox/test_population.py
new file mode 100644
index 000000000..7bd3a070b
--- /dev/null
+++ b/tests/brainbox/test_population.py
@@ -0,0 +1,62 @@
+from brainbox.population import xcorr
+import unittest
+import numpy as np
+
+
+def _random_data(max_cluster):
+    nspikes = 10000
+    spike_times = np.cumsum(np.random.exponential(scale=.025, size=nspikes))
+    spike_clusters = np.random.randint(0, max_cluster, nspikes)
+    return spike_times, spike_clusters
+
+
+class TestPopulation(unittest.TestCase):
+    def test_xcorr_0(self):
+        # 0: 0, 10
+        # 1: 10, 20
+        spike_times = np.array([0, 10, 10, 20])
+        spike_clusters = np.array([0, 1, 0, 1])
+        bin_size = 1
+        winsize_bins = 2 * 3 + 1
+
+        c_expected = np.zeros((2, 2, 7), dtype=np.int32)
+
+        c_expected[1, 0, 3] = 1
+        c_expected[0, 1, 3] = 1
+
+        c = xcorr(spike_times, spike_clusters, bin_size=bin_size, window_size=winsize_bins)
+
+        self.assertTrue(np.allclose(c, c_expected))
+
+    def test_xcorr_1(self):
+        # 0: 2, 10, 12, 30
+        # 1: 3, 24
+        # 2: 20
+        spike_times = np.array([2, 3, 10, 12, 20, 24, 30, 40], dtype=np.uint64)
+        spike_clusters = np.array([0, 1, 0, 0, 2, 1, 0, 2])
+        bin_size = 1
+        winsize_bins = 2 * 3 + 1
+
+        c_expected = np.zeros((3, 3, 7), dtype=np.int32)
+        c_expected[0, 1, 4] = 1
+        c_expected[1, 0, 2] = 1
+        c_expected[0, 0, 1] = 1
+        c_expected[0, 0, 5] = 1
+
+        c = xcorr(spike_times, spike_clusters, bin_size=bin_size, window_size=winsize_bins)
+
+        self.assertTrue(np.allclose(c, c_expected))
+
+    def test_xcorr_2(self):
+        max_cluster = 10
+        spike_times, spike_clusters = _random_data(max_cluster)
+        bin_size, winsize_bins = .001, .05
+
+        c = xcorr(spike_times, spike_clusters, bin_size=bin_size, window_size=winsize_bins)
+
+        self.assertEqual(c.shape, (max_cluster, max_cluster, 51))
+
+
+if __name__ == "__main__":
+    np.random.seed(0)
+    unittest.main(exit=False)
diff --git a/tests/brainbox/test_processing.py b/tests/brainbox/test_processing.py
index 589271159..3c553cf87 100644
--- a/tests/brainbox/test_processing.py
+++ b/tests/brainbox/test_processing.py
@@ -13,7 +13,7 @@ def test_sync(self):
         samples = times**3
         # Use cubic interpolation to resample to uniform interval
         cubes = core.TimeSeries(times=times, values=samples, columns=('cubic',))
-        resamp = processing.sync(0.1, cubes, interp='cubic', fillval='extrapolate')
+        resamp = processing.sync(0.1, timeseries=cubes, interp='cubic', fillval='extrapolate')
         # Check that the sync function is returning a new time series object
         self.assertTrue(isinstance(resamp, core.TimeSeries))
         # Test that all returned sample times are uniformly spaced
@@ -29,7 +29,8 @@ def test_sync(self):
         samples2 = times2**2
         squares = core.TimeSeries(times=times2, values=samples2, columns=('square',))
         # Use cubic interpolation again, this time on both timeseries
-        resamp2 = processing.sync(0.1, [squares, cubes], interp='cubic', fillval='extrapolate')
+        resamp2 = processing.sync(0.1, timeseries=[squares, cubes], interp='cubic',
+                                  fillval='extrapolate')
         # Check that the new TS has both squares and cubes as keys and attribs
         self.assertTrue(hasattr(resamp2, 'cubic'))
         self.assertTrue(hasattr(resamp2, 'square'))
diff --git a/tests/brainbox/test_singlecell.py b/tests/brainbox/test_singlecell.py
new file mode 100644
index 000000000..29ff90353
--- /dev/null
+++ b/tests/brainbox/test_singlecell.py
@@ -0,0 +1,22 @@
+from brainbox.singlecell import acorr
+import unittest
+import numpy as np
+
+
+class TestPopulation(unittest.TestCase):
+    def test_acorr_0(self):
+        spike_times = np.array([0, 10, 10, 20])
+        bin_size = 1
+        winsize_bins = 2 * 3 + 1
+
+        c_expected = np.zeros(7, dtype=np.int32)
+        c_expected[3] = 1
+
+        c = acorr(spike_times, bin_size=bin_size, window_size=winsize_bins)
+
+        self.assertTrue(np.allclose(c, c_expected))
+
+
+if __name__ == "__main__":
+    np.random.seed(0)
+    unittest.main(exit=False)