multielectrode_grasp/code/data_overview_2.py

# -*- coding: utf-8 -*-
"""
Code for generating the second data figure in the manuscript.

Authors: Julia Sprenger, Lyuba Zehl, Michael Denker


Copyright (c) 2017, Institute of Neuroscience and Medicine (INM-6),
Forschungszentrum Juelich, Germany
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
* Neither the names of the copyright holders nor the names of the contributors
may be used to endorse or promote products derived from this software without
specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
"""

# This loads the Neo and odML libraries shipped with this code. For production
# use, please use the newest releases of odML and Neo.
import load_local_neo_odml_elephant

import os

import matplotlib.pyplot as plt
from matplotlib import gridspec, transforms

import quantities as pq
import numpy as np

from neo import (AnalogSignal, SpikeTrain)
from neo.utils import *
from reachgraspio import reachgraspio


from neo_utils import load_segment

# =============================================================================
# Define data and metadata directories and general settings
# =============================================================================


def get_monkey_datafile(monkey):
    if monkey == "Lilou":
        return "l101210-001"  # ns2 (behavior) and ns5 present
    elif monkey == "Nikos2":
        return "i140703-001"  # ns2 and ns6 present
    else:
        return ""


# Enter your dataset directory here
datasetdir = os.path.join('..', 'datasets')

chosen_els = {'Lilou': range(3, 97, 7), 'Nikos2': range(1, 97, 7)}
chosen_el = {
    'Lilou': chosen_els['Lilou'][0],
    'Nikos2': chosen_els['Nikos2'][0]}
chosen_unit = 1
trial_indexes = range(14)
trial_index = trial_indexes[0]
chosen_events = ['TS-ON', 'WS-ON', 'CUE-ON', 'CUE-OFF', 'GO-ON', 'SR-ON',
                 'RW-ON', 'WS-OFF']  # , 'RW-OFF'

# =============================================================================
# Load data and metadata for a monkey
# =============================================================================
# CHANGE this parameter to load data of the different monkeys
# monkey = 'Nikos2'
monkey = 'Lilou'

datafile = get_monkey_datafile(monkey)

session = reachgraspio.ReachGraspIO(
    filename=os.path.join(datasetdir, datafile),
    odml_directory=datasetdir,
    verbose=False)

bl = session.read_block(lazy=True)
    # channels=chosen_els[monkey],
    # units=[1],  # loading only unit_id 1
    # load_waveforms=False,
    # load_events=True,
    # scaling='voltage')

seg = bl.segments[0]

# get start and stop events of trials
start_events = get_events(
    seg, **{
        'name': 'TrialEvents',
        'trial_event_labels': 'TS-ON',
        'performance_in_trial': session.performance_codes['correct_trial']})
stop_events = get_events(
    seg, **{
        'name': 'TrialEvents',
        'trial_event_labels': 'RW-ON',
        'performance_in_trial': session.performance_codes['correct_trial']})

# there should only be one event object for these conditions
assert len(start_events) == 1
assert len(stop_events) == 1

# insert epochs between 10ms before TS to 50ms after RW corresponding to trails
ep = add_epoch(seg,
               start_events[0],
               stop_events[0],
               pre=-250 * pq.ms,
               post=500 * pq.ms,
               segment_type='complete_trials')
ep.array_annotate(trialtype=start_events[0].array_annotations['belongs_to_trialtype'])

# access single epoch of this data_segment
epochs = get_epochs(seg, **{'segment_type': 'complete_trials'})
assert len(epochs) == 1

# remove spiketrains not belonging to chosen_electrode
seg.spiketrains = seg.filter(targdict={'unit_id': chosen_unit},
                             recursive=True, objects='SpikeTrainProxy')
# remove all non-neural signals
seg.analogsignals = seg.filter(targdict={'neural_signal': True},
                               objects='AnalogSignalProxy')

# use most raw data if multiple versions are present
raw_signal = seg.analogsignals[0]
for sig in seg.analogsignals:
    if sig.sampling_rate > raw_signal.sampling_rate:
        raw_signal = sig
seg.analogsignals = [raw_signal]

# replacing the segment with a new segment containing all data
# to speed up cutting of segments
seg = load_segment(seg, load_wavefroms=True)

# only keep the chosen electrode signal in the AnalogSignal object
mask = np.isin(seg.analogsignals[0].array_annotations['channel_ids'], chosen_els[monkey])

seg.analogsignals[0] = seg.analogsignals[0][:, mask]

# cut segments according to inserted 'complete_trials' epochs and reset trial
# times
cut_segments = cut_segment_by_epoch(seg, epochs[0], reset_time=True)

# explicitly adding trial type annotations to cut segments
for i, cut_seg in enumerate(cut_segments):
    cut_seg.annotate(trialtype=epochs[0].array_annotations['trialtype'][i])

# =============================================================================
# Define figure and subplot axis for first data overview
# =============================================================================
fig = plt.figure(facecolor='w')
fig.set_size_inches(7.0, 9.9)  # (w, h) in inches
# #(7.0, 9.9) corresponds to A4 portrait ratio

gs = gridspec.GridSpec(
    nrows=2,
    ncols=2,
    left=0.1,
    bottom=0.05,
    right=0.9,
    top=0.975,
    wspace=0.1,
    hspace=0.1,
    width_ratios=None,
    height_ratios=[2, 1])

ax1 = plt.subplot(gs[0, 0])  # top left
ax2 = plt.subplot(gs[0, 1], sharex=ax1)  # top right
ax3 = plt.subplot(gs[1, 0], sharex=ax1)  # bottom left
ax4 = plt.subplot(gs[1, 1], sharex=ax1)  # bottom right

fontdict_titles = {'fontsize': 9, 'fontweight': 'bold'}
fontdict_axis = {'fontsize': 10, 'fontweight': 'bold'}

# the x coords of the event labels are data, and the y coord are axes
event_label_transform = transforms.blended_transform_factory(ax1.transData,
                                                             ax1.transAxes)

trialtype_colors = {
    'SGHF': 'MediumBlue', 'SGLF': 'Turquoise',
    'PGHF': 'DarkGreen', 'PGLF': 'YellowGreen',
    'LFSG': 'Orange', 'LFPG': 'Yellow',
    'HFSG': 'DarkRed', 'HFPG': 'OrangeRed',
    'SGSG': 'SteelBlue', 'PGPG': 'LimeGreen',
    None: 'black'}

event_colors = {
    'TS-ON': 'indigo', 'TS-OFF': 'indigo',
    'WS-ON': 'purple', 'WS-OFF': 'purple',
    'CUE-ON': 'crimson', 'CUE-OFF': 'crimson',
    'GO-ON': 'orangered', 'GO-OFF': 'orangered',
    'SR-ON': 'darkorange',
    'RW-ON': 'orange', 'RW-OFF': 'orange'}

electrode_cmap = plt.get_cmap('bone')
electrode_colors = [electrode_cmap(x) for x in
                    np.tile(np.array([0.3, 0.7]), int(len(chosen_els[monkey]) / 2))]

time_unit = 'ms'
lfp_unit = 'uV'

# define scaling factors for analogsignals

anasig_std = np.mean(np.std(cut_segments[trial_index].analogsignals[0].rescale(lfp_unit), axis=0))
anasig_offset = 3 * anasig_std


# =============================================================================
# SUPPLEMENTARY PLOTTING functions
# =============================================================================

def add_scalebar(ax, std):
    # the x coords of the scale bar are axis, and the y coord are data
    scalebar_transform = transforms.blended_transform_factory(ax.transAxes,
                                                              ax.transData)
    # adding scalebar
    yscalebar = max(int(std.rescale(lfp_unit)), 1) * getattr(pq, lfp_unit) * 2
    scalebar_offset = -2 * std
    ax.vlines(x=0.4,
              ymin=(scalebar_offset - yscalebar).magnitude,
              ymax=scalebar_offset.magnitude,
              color='k',
              linewidth=4,
              transform=scalebar_transform)
    ax.text(0.4, (scalebar_offset - 0.5 * yscalebar).magnitude,
            ' %i %s' % (yscalebar.magnitude, lfp_unit),
            ha="left", va="center", rotation=0, color='k',
            size=8, transform=scalebar_transform)


# =============================================================================
# PLOT DATA OF SINGLE TRIAL (left plots)
# =============================================================================

# get data of selected trial
selected_trial = cut_segments[trial_index]

# PLOT DATA FOR EACH CHOSEN ELECTRODE
for el_idx, electrode_id in enumerate(chosen_els[monkey]):

    # PLOT ANALOGSIGNALS in upper plot
    chosen_el_idx = np.where(cut_segments[0].analogsignals[0].array_annotations['channel_ids'] == electrode_id)[0][0]
    anasig = selected_trial.analogsignals[0][:, chosen_el_idx]
    ax1.plot(anasig.times.rescale(time_unit),
             np.asarray(anasig.rescale(lfp_unit))
             + anasig_offset.magnitude * el_idx,
             color=electrode_colors[el_idx])

    # PLOT SPIKETRAINS in lower plot
    spiketrains = selected_trial.filter(
        channel_id=electrode_id, objects=SpikeTrain)
    for spiketrain in spiketrains:
        ax3.plot(spiketrain.times.rescale(time_unit),
                 np.zeros(len(spiketrain.times)) + el_idx, 'k|')

# PLOT EVENTS in both plots
for event_type in chosen_events:
    # get events of each chosen event type
    event_data = get_events(selected_trial,
                                      **{'trial_event_labels': event_type})
    for event in event_data:
        event_color = event_colors[event.array_annotations['trial_event_labels'][0]]
        # adding lines
        for ax in [ax1, ax3]:
            ax.axvline(event.times.rescale(time_unit),
                       color=event_color,
                       zorder=0.5)
        # adding labels
        ax1.text(event.times.rescale(time_unit), 0,
                 event.array_annotations['trial_event_labels'][0],
                 ha="center", va="top", rotation=45, color=event_color,
                 size=8, transform=event_label_transform)

# SUBPLOT ADJUSTMENTS
ax1.set_title('single trial', fontdict=fontdict_titles)
ax1.set_ylabel('electrode id', fontdict=fontdict_axis)
ax1.set_yticks(np.arange(len(chosen_els[monkey])) * anasig_offset)
ax1.set_yticklabels(chosen_els[monkey])

ax1.autoscale(enable=True, axis='y')
plt.setp(ax1.get_xticklabels(), visible=False)  # show no xticklabels
ax3.set_ylabel('electrode id', fontdict=fontdict_axis)
ax3.set_yticks(range(0, len(chosen_els[monkey])))
ax3.set_yticklabels(np.asarray(chosen_els[monkey]))
ax3.set_ylim(-1, len(chosen_els[monkey]))
ax3.set_xlabel('time [%s]' % time_unit, fontdict=fontdict_axis)
# ax3.autoscale(axis='y')

# =============================================================================
# PLOT DATA OF SINGLE ELECTRODE
# =============================================================================


# plot data for each chosen trial
chosen_el_idx = np.where(cut_segments[0].analogsignals[0].array_annotations['channel_ids'] == chosen_el[monkey])[0][0]
for trial_idx, trial_id in enumerate(trial_indexes):
    trial_spikes = cut_segments[trial_id].filter(channel_id=chosen_el[monkey], objects='SpikeTrain')
    trial_type = cut_segments[trial_id].annotations['trialtype']
    trial_color = trialtype_colors[trial_type]

    t_signal = cut_segments[trial_id].analogsignals[0][:, chosen_el_idx]
    # PLOT ANALOGSIGNALS in upper plot
    ax2.plot(t_signal.times.rescale(time_unit),
             np.asarray(t_signal.rescale(lfp_unit))
             + anasig_offset.magnitude * trial_idx,
             color=trial_color, zorder=1)

    for t_data in trial_spikes:
        # PLOT SPIKETRAINS in lower plot
        ax4.plot(t_data.times.rescale(time_unit),
                 np.ones(len(t_data.times)) + trial_idx, 'k|')

    # PLOT EVENTS in both plots
    for event_type in chosen_events:
        # get events of each chosen event type
        event_data = get_events(cut_segments[trial_id], **{'trial_event_labels': event_type})
        for event in event_data:
            color = event_colors[event.array_annotations['trial_event_labels'][0]]
            ax2.vlines(x=event.times.rescale(time_unit),
                       ymin=(trial_idx - 0.5) * anasig_offset,
                       ymax=(trial_idx + 0.5) * anasig_offset,
                       color=color,
                       zorder=2)
            ax4.vlines(x=event.times.rescale(time_unit),
                       ymin=trial_idx + 1 - 0.4,
                       ymax=trial_idx + 1 + 0.4,
                       color=color,
                       zorder=0.5)

# SUBPLOT ADJUSTMENTS
ax2.set_title('single electrode', fontdict=fontdict_titles)
ax2.set_ylabel('trial id', fontdict=fontdict_axis)
ax2.set_yticks(np.asarray(trial_indexes) * anasig_offset)
ax2.set_yticklabels(
    [epochs[0].array_annotations['trial_id'][_] for _ in trial_indexes])
ax2.yaxis.set_label_position("right")
ax2.tick_params(direction='in', length=3, labelleft='off', labelright='on')
ax2.autoscale(enable=True, axis='y')
add_scalebar(ax2, anasig_std)
plt.setp(ax2.get_xticklabels(), visible=False)  # show no xticklabels

ax4.set_ylabel('trial id', fontdict=fontdict_axis)
ax4.set_xlabel('time [%s]' % time_unit, fontdict=fontdict_axis)

start, end = ax4.get_xlim()
ax4.xaxis.set_ticks(np.arange(start, end, 1000))
ax4.xaxis.set_ticks(np.arange(start, end, 500), minor=True)
ax4.set_yticks(range(1, len(trial_indexes) + 1))
ax4.set_yticklabels(np.asarray(
    [epochs[0].array_annotations['trial_id'][_] for _ in trial_indexes]))
ax4.yaxis.set_label_position("right")
ax4.tick_params(direction='in', length=3, labelleft='off', labelright='on')
ax4.autoscale(enable=True, axis='y')

# GENERAL PLOT ADJUSTMENTS
# adjust font sizes of ticks
for ax in [ax4.yaxis, ax4.xaxis, ax3.xaxis, ax3.yaxis]:
    for tick in ax.get_major_ticks():
        tick.label.set_fontsize(10)

# adjust time range on x axis
t_min = np.min([cut_segments[tid].t_start.rescale(time_unit)
                for tid in trial_indexes])
t_max = np.max([cut_segments[tid].t_stop.rescale(time_unit)
                for tid in trial_indexes])
ax1.set_xlim(t_min, t_max)
add_scalebar(ax1, anasig_std)

# =============================================================================
# SAVE FIGURE
# =============================================================================
fname = 'data_overview_2_%s' % monkey
for file_format in ['eps', 'pdf', 'png']:
    fig.savefig(fname + '.%s' % file_format, dpi=400, format=file_format)