Chained protocols

brainbox/behavior/training.py

@@ -1,4 +1,5 @@
 from one.api import ONE
+from one.alf.exceptions import ALFObjectNotFound
 import datetime
 import re
 import numpy as np
@@ -131,7 +132,10 @@ def get_sessions(subj, date=None, one=None):
     sess_dates = []
     if len(sessions) < 3:
         for n, _ in enumerate(sessions):
-            trials_ = one.load_object(sessions[n]['url'].split('/')[-1], 'trials')
+            try:
+                trials_ = one.load_object(sessions[n]['url'].split('/')[-1], 'trials')
+            except ALFObjectNotFound:
+                trials_ = None
 
             if trials_:
                 task_protocol.append(re.search('tasks_(.*)Choice',
@@ -142,7 +146,11 @@ def get_sessions(subj, date=None, one=None):
     else:
         n = 0
         while len(trials) < 3:
-            trials_ = one.load_object(sessions[n]['url'].split('/')[-1], 'trials')
+            print(sessions[n]['url'].split('/')[-1])
+            try:
+                trials_ = one.load_object(sessions[n]['url'].split('/')[-1], 'trials')
+            except ALFObjectNotFound:
+                trials_ = None
 
             if trials_:
                 task_protocol.append(re.search('tasks_(.*)Choice',

brainbox/io/one.py

@@ -706,6 +706,28 @@ def load_wheel_reaction_times(eid, one=None):
     return firstMove_times - trials['goCue_times']
 
 
+def load_iti(trials):
+    """
+    The inter-trial interval (ITI) time for each trial, defined as the period of open-loop grey
+    screen commencing at stimulus off and lasting until the quiescent period at the start of the
+    following trial.  Note that the ITI for the first trial is the time between the first trial
+    and the next, therefore the last value is NaN.
+
+    Parameters
+    ----------
+    trials : one.alf.io.AlfBunch
+        An ALF trials object containing the keys {'intervals', 'stimOff_times'}.
+
+    Returns
+    -------
+    np.array
+        An array of inter-trial intervals, the last value being NaN.
+    """
+    if not {'intervals', 'stimOff_times'} <= trials.keys():
+        raise ValueError('trials must contain keys {"intervals", "stimOff_times"}')
+    return np.r_[(np.roll(trials['intervals'][:, 0], -1) - trials['stimOff_times'])[:-1], np.nan]
+
+
 def load_channels_from_insertion(ins, depths=None, one=None, ba=None):
 
     PROV_2_VAL = {

brainbox/io/spikeglx.py

@@ -132,7 +132,7 @@ def __init__(self, pid, one, typ='ap', cache_folder=None, remove_cached=False):
             self.chunks['chunk_bounds'] = np.array(self.chunks['chunk_bounds'])
         super(Streamer, self).__init__(meta_file, ignore_warnings=True)
 
-    def read(self, nsel=slice(0, 10000), csel=slice(None), sync=True):
+    def read(self, nsel=slice(0, 10000), csel=slice(None), sync=True, volts=True):
         """
         overload the read function by downloading the necessary chunks
         """
@@ -142,7 +142,11 @@ def read(self, nsel=slice(0, 10000), csel=slice(None), sync=True):
         _logger.debug(f'Streamer: caching sample {n0}, (t={n0 / self.fs})')
         self.cache_folder.mkdir(exist_ok=True, parents=True)
         sr = self._download_raw_partial(first_chunk=first_chunk, last_chunk=last_chunk)
-        data = sr[nsel.start - n0: nsel.stop - n0, csel]
+        if not volts:
+            data = np.copy(sr._raw[nsel.start - n0:nsel.stop - n0, csel])
+        else:
+            data = sr[nsel.start - n0: nsel.stop - n0, csel]
+
         sr.close()
         if self.remove_cached:
             shutil.rmtree(self.target_dir)

brainbox/plot.py

@@ -4,15 +4,15 @@
 
 Run the following to set-up the workspace to run the docstring examples:
 >>> from brainbox import processing
->>> import alf.io as aio
+>>> import one.alf.io as alfio
 >>> import numpy as np
 >>> import matplotlib.pyplot as plt
 >>> import ibllib.ephys.spikes as e_spks
 # (*Note, if there is no 'alf' directory, make 'alf' directory from 'ks2' output directory):
 >>> e_spks.ks2_to_alf(path_to_ks_out, path_to_alf_out)
 # Load the alf spikes bunch and clusters bunch, and get a units bunch.
->>> spks_b = aio.load_object(path_to_alf_out, 'spikes')
->>> clstrs_b = aio.load_object(path_to_alf_out, 'clusters')
+>>> spks_b = alfio.load_object(path_to_alf_out, 'spikes')
+>>> clstrs_b = alfio.load_object(path_to_alf_out, 'clusters')
 >>> units_b = processing.get_units_bunch(spks_b)  # may take a few mins to compute
 """
 
@@ -28,7 +28,7 @@
 from brainbox.metrics import single_units
 from brainbox.processing import bincount2D
 from brainbox.io.spikeglx import extract_waveforms
-from ibllib.io import spikeglx
+import spikeglx
 
 
 def feat_vars(units_b, units=None, feat_name='amps', dist='norm', test='ks', cmap_name='coolwarm',

brainbox/tests/test_io.py

@@ -5,6 +5,7 @@
 import shutil
 
 import numpy as np
+import numpy.testing
 
 from brainbox.io import one as bbone
 from one.api import ONE
@@ -66,5 +67,24 @@ def tearDown(self) -> None:
         shutil.rmtree(self.tmpdir)
 
 
-if __name__ == "__main__":
+class TestIO_ONE(unittest.TestCase):
+    """Tests for brainbox.io.one functions that don't require fixtures on disk."""
+    def test_load_iti(self):
+        """Test for brainbox.io.one.load_iti function."""
+        trials = bbone.alfio.AlfBunch({})
+        trials.intervals = np.array([
+            [114.52487625, 117.88103707],
+            [118.5169474, 122.89742147],
+            [123.49302927, 126.12216664],
+            [126.68107337, 129.53872083],
+            [130.11952807, 133.90539162]
+        ])
+        trials.stimOff_times = [117.38098379, 122.39736201, 125.62210278, 129.03865947, 133.4053633]
+        expected = np.array([1.13596361, 1.09566726, 1.05897059, 1.0808686, np.nan])
+        np.testing.assert_array_almost_equal(bbone.load_iti(trials), expected)
+        _ = trials.pop('stimOff_times')
+        self.assertRaises(ValueError, bbone.load_iti, trials)
+
+
+if __name__ == '__main__':
     unittest.main(exit=False, verbosity=2)

examples/exploring_data/data_download.ipynb

@@ -365,4 +365,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 0
-}
+}

ibllib/__init__.py

@@ -2,7 +2,7 @@
 import logging
 import warnings
 
-__version__ = '2.20.0'
+__version__ = '2.21.0'
 warnings.filterwarnings('always', category=DeprecationWarning, module='ibllib')
 
 # if this becomes a full-blown library we should let the logging configuration to the discretion of the dev

ibllib/atlas/allen_structure_tree.csv

@@ -367,12 +367,12 @@ id,atlas_id,name,acronym,st_level,ontology_id,hemisphere_id,weight,parent_struct
 289,1026,Temporal association areas layer 5,TEa5,,1,3,8690,541,7,1,364,/997/8/567/688/695/315/541/289/,15B0B3,,,f,365,475299804,734881840,Temporal association areas layer 5
 729,1081,Temporal association areas layer 6a,TEa6a,,1,3,8690,541,7,1,365,/997/8/567/688/695/315/541/729/,15B0B3,,,f,366,1289955072,734881840,Temporal association areas layer 6a
 786,1088,Temporal association areas layer 6b,TEa6b,,1,3,8690,541,7,1,366,/997/8/567/688/695/315/541/786/,15B0B3,,,f,367,3588912826,734881840,Temporal association areas layer 6b
-922,256,Perirhinal area,PERI,,1,3,8690,315,6,1,367,/997/8/567/688/695/315/922/,0,,,f,368,173995527,734881840,Perirhinal area
-540,1057,Perirhinal area layer 1,PERI1,,1,3,8690,922,7,1,368,/997/8/567/688/695/315/922/540/,0,,,f,369,161492388,734881840,Perirhinal area layer 1
-888,959,Perirhinal area layer 2/3,PERI2/3,,1,3,8690,922,7,1,369,/997/8/567/688/695/315/922/888/,0,,,f,370,1642584549,734881840,Perirhinal area layer 2/3
-692,1076,Perirhinal area layer 5,PERI5,,1,3,8690,922,7,1,370,/997/8/567/688/695/315/922/692/,0,,,f,371,248375741,734881840,Perirhinal area layer 5
-335,890,Perirhinal area layer 6a,PERI6a,,1,3,8690,922,7,1,371,/997/8/567/688/695/315/922/335/,0,,,f,372,1984229208,734881840,Perirhinal area layer 6a
-368,894,Perirhinal area layer 6b,PERI6b,,1,3,8690,922,7,1,372,/997/8/567/688/695/315/922/368/,0,,,f,373,4014849762,734881840,Perirhinal area layer 6b
+922,256,Perirhinal area,PERI,,1,3,8690,315,6,1,367,/997/8/567/688/695/315/922/,0E9684,,,f,368,173995527,734881840,Perirhinal area
+540,1057,Perirhinal area layer 1,PERI1,,1,3,8690,922,7,1,368,/997/8/567/688/695/315/922/540/,0E9684,,,f,369,161492388,734881840,Perirhinal area layer 1
+888,959,Perirhinal area layer 2/3,PERI2/3,,1,3,8690,922,7,1,369,/997/8/567/688/695/315/922/888/,0E9684,,,f,370,1642584549,734881840,Perirhinal area layer 2/3
+692,1076,Perirhinal area layer 5,PERI5,,1,3,8690,922,7,1,370,/997/8/567/688/695/315/922/692/,0E9684,,,f,371,248375741,734881840,Perirhinal area layer 5
+335,890,Perirhinal area layer 6a,PERI6a,,1,3,8690,922,7,1,371,/997/8/567/688/695/315/922/335/,0E9684,,,f,372,1984229208,734881840,Perirhinal area layer 6a
+368,894,Perirhinal area layer 6b,PERI6b,,1,3,8690,922,7,1,372,/997/8/567/688/695/315/922/368/,0E9684,,,f,373,4014849762,734881840,Perirhinal area layer 6b
 895,111,Ectorhinal area,ECT,,1,3,8690,315,6,1,373,/997/8/567/688/695/315/895/,0D9F91,,,f,374,3399607880,734881840,Ectorhinal area
 836,1094,Ectorhinal area/Layer 1,ECT1,,1,3,8690,895,7,1,374,/997/8/567/688/695/315/895/836/,0D9F91,,,f,375,2240743692,734881840,Ectorhinal area/Layer 1
 427,1043,Ectorhinal area/Layer 2/3,ECT2/3,,1,3,8690,895,7,1,375,/997/8/567/688/695/315/895/427/,0D9F91,,,f,376,993691642,734881840,Ectorhinal area/Layer 2/3

ibllib/atlas/atlas.py

@@ -1058,6 +1058,22 @@ def _ccf_order(ccf_order, reverse=False):
         else:
             ValueError("ccf_order needs to be either 'mlapdv' or 'apdvml'")
 
+    def compute_regions_volume(self):
+        """
+        Sums the number of voxels in the labels volume for each region.
+        Then compute volumes for all of the levels of hierarchy in cubic mm.
+        :return:
+        """
+        nr = self.regions.id.shape[0]
+        count = np.bincount(self.label.flatten(), minlength=nr)
+        self.regions.compute_hierarchy()
+        self.regions.volume = np.zeros_like(count)
+        for i in np.arange(nr):
+            if count[i] == 0:
+                continue
+            self.regions.volume[np.unique(self.regions.hierarchy[:, i])] += count[i]
+        self.regions.volume = self.regions.volume * (self.res_um / 1e3) ** 3
+
 
 def NeedlesAtlas(*args, **kwargs):
     """

ibllib/atlas/flatmaps.py

@@ -3,18 +3,21 @@
 """
 from functools import lru_cache
 import logging
+import json
 
 import numpy as np
 from scipy.interpolate import interp1d
 import matplotlib.pyplot as plt
 import matplotlib.colors
+from matplotlib import cm
 
 from iblutil.numerical import ismember
 from iblutil.util import Bunch
 from iblutil.io.hashfile import md5
 import one.remote.aws as aws
 
 from ibllib.atlas.atlas import AllenAtlas, BrainRegions
+from ibllib.atlas.plots import plot_polygon, plot_polygon_with_hole, coords_for_poly_hole
 
 
 _logger = logging.getLogger(__name__)
@@ -132,6 +135,115 @@ def swanson(filename="swanson2allen.npz"):
     return s2a
 
 
+def swanson_json(filename="swansonpaths.json"):
+
+    OLD_MD5 = ['f848783954883c606ca390ceda9e37d2']
+
+    json_file = AllenAtlas._get_cache_dir().joinpath(filename)
+    if not json_file.exists() or md5(json_file) in OLD_MD5:
+        json_file.parent.mkdir(exist_ok=True, parents=True)
+        _logger.info(f'downloading swanson paths from {aws.S3_BUCKET_IBL} s3 bucket...')
+        aws.s3_download_file(f'atlas/{json_file.name}', json_file)
+
+    with open(json_file) as f:
+        sw_json = json.load(f)
+
+    return sw_json
+
+
+def plot_swanson_vector(acronyms=None, values=None, ax=None, hemisphere=None, br=None, orientation='landscape',
+                        empty_color='silver', vmin=None, vmax=None, cmap='cividis', **kwargs):
+
+    br = BrainRegions() if br is None else br
+    br.compute_hierarchy()
+
+    if ax is None:
+        fig, ax = plt.subplots()
+        ax.set_axis_off()
+
+    if acronyms is not None:
+        ibr, vals = br.propagate_down(acronyms, values)
+        colormap = cm.get_cmap(cmap)
+        vmin = vmin or np.nanmin(vals)
+        vmax = vmax or np.nanmax(vals)
+        norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
+        rgba_color = colormap(norm(vals), bytes=True)
+
+    sw = swanson()
+    sw_json = swanson_json()
+
+    for i, reg in enumerate(sw_json):
+
+        if acronyms is None:
+            color = br.rgba[br.mappings['Swanson'][reg['thisID']]] / 255
+        else:
+            idx = np.where(ibr == reg['thisID'])[0]
+            color = rgba_color[idx[0]] / 255 if len(idx) > 0 else empty_color
+
+        coords = reg['coordsReg']
+
+        if reg['hole']:
+            vertices, codes = coords_for_poly_hole(coords)
+            if orientation == 'portrait':
+                vertices[:, [0, 1]] = vertices[:, [1, 0]]
+                plot_polygon_with_hole(ax, vertices, codes, color, **kwargs)
+                if hemisphere is not None:
+                    color_inv = color if hemisphere == 'mirror' else empty_color
+                    vertices_inv = np.copy(vertices)
+                    vertices_inv[:, 0] = -1 * vertices_inv[:, 0] + (sw.shape[0] * 2)
+                    plot_polygon_with_hole(ax, vertices_inv, codes, color_inv, **kwargs)
+            else:
+                plot_polygon_with_hole(ax, vertices, codes, color, **kwargs)
+                if hemisphere is not None:
+                    color_inv = color if hemisphere == 'mirror' else empty_color
+                    vertices_inv = np.copy(vertices)
+                    vertices_inv[:, 1] = -1 * vertices_inv[:, 1] + (sw.shape[0] * 2)
+                    plot_polygon_with_hole(ax, vertices_inv, codes, color_inv, **kwargs)
+        else:
+            coords = [coords] if type(coords) == dict else coords
+            for c in coords:
+
+                if orientation == 'portrait':
+                    xy = np.c_[c['y'], c['x']]
+                    plot_polygon(ax, xy, color, **kwargs)
+                    if hemisphere is not None:
+                        color_inv = color if hemisphere == 'mirror' else empty_color
+                        xy_inv = np.copy(xy)
+                        xy_inv[:, 0] = -1 * xy_inv[:, 0] + (sw.shape[0] * 2)
+                        plot_polygon(ax, xy_inv, color_inv, **kwargs)
+                else:
+                    xy = np.c_[c['x'], c['y']]
+                    plot_polygon(ax, xy, color, **kwargs)
+                    if hemisphere is not None:
+                        color_inv = color if hemisphere == 'mirror' else empty_color
+                        xy_inv = np.copy(xy)
+                        xy_inv[:, 1] = -1 * xy_inv[:, 1] + (sw.shape[0] * 2)
+                        plot_polygon(ax, xy_inv, color_inv, **kwargs)
+
+    if orientation == 'portrait':
+        ax.set_ylim(0, sw.shape[1])
+        if hemisphere is None:
+            ax.set_xlim(0, sw.shape[0])
+        else:
+            ax.set_xlim(0, 2 * sw.shape[0])
+    else:
+        ax.set_xlim(0, sw.shape[1])
+        if hemisphere is None:
+            ax.set_ylim(0, sw.shape[0])
+        else:
+            ax.set_ylim(0, 2 * sw.shape[0])
+
+    def format_coord(x, y):
+        ind = sw[int(y), int(x)]
+        ancestors = br.ancestors(br.id[ind])['acronym']
+        return f'sw-{ind}, x={x:1.4f}, y={y:1.4f}, aid={br.id[ind]}-{br.acronym[ind]} \n {ancestors}'
+
+    ax.format_coord = format_coord
+
+    ax.invert_yaxis()
+    ax.set_aspect('equal')
+
+
 def plot_swanson(acronyms=None, values=None, ax=None, hemisphere=None, br=None,
                  orientation='landscape', annotate=False, empty_color='silver', **kwargs):
     """
@@ -196,7 +308,7 @@ def plot_swanson(acronyms=None, values=None, ax=None, hemisphere=None, br=None,
     def format_coord(x, y):
         ind = s2a[int(y), int(x)]
         ancestors = br.ancestors(br.id[ind])['acronym']
-        return f'x={x:1.4f}, y={x:1.4f}, {br.acronym[ind]} \n {ancestors}'
+        return f'sw-{ind}, x={x:1.4f}, y={y:1.4f}, aid={br.id[ind]}-{br.acronym[ind]} \n {ancestors}'
 
     ax.format_coord = format_coord
     return ax