1
1

data_overview_1.py 28 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744
  1. # -*- coding: utf-8 -*-
  2. """
  3. Code for generating the first data figure in the manuscript.
  4. Authors: Julia Sprenger, Lyuba Zehl, Michael Denker
  5. Copyright (c) 2017, Institute of Neuroscience and Medicine (INM-6),
  6. Forschungszentrum Juelich, Germany
  7. All rights reserved.
  8. Redistribution and use in source and binary forms, with or without
  9. modification, are permitted provided that the following conditions are met:
  10. * Redistributions of source code must retain the above copyright notice, this
  11. list of conditions and the following disclaimer.
  12. * Redistributions in binary form must reproduce the above copyright notice,
  13. this list of conditions and the following disclaimer in the documentation
  14. and/or other materials provided with the distribution.
  15. * Neither the names of the copyright holders nor the names of the contributors
  16. may be used to endorse or promote products derived from this software without
  17. specific prior written permission.
  18. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  19. ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  20. WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  21. DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
  22. FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  23. DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  24. SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  25. CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  26. OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  27. OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  28. """
  29. # This loads the Neo and odML libraries shipped with this code. For production
  30. # use, please use the newest releases of odML and Neo.
  31. import load_local_neo_odml_elephant
  32. import os
  33. import numpy as np
  34. from scipy import stats
  35. import quantities as pq
  36. import matplotlib.pyplot as plt
  37. from matplotlib import gridspec, ticker
  38. from reachgraspio import reachgraspio
  39. import odml.tools
  40. from neo import utils as neo_utils
  41. from neo_utils import load_segment
  42. import odml_utils
  43. # =============================================================================
  44. # Define data and metadata directories
  45. # =============================================================================
  46. def get_monkey_datafile(monkey):
  47. if monkey == "Lilou":
  48. return "l101210-001" # ns2 (behavior) and ns5 present
  49. elif monkey == "Nikos2":
  50. return "i140703-001" # ns2 and ns6 present
  51. else:
  52. return ""
  53. # Enter your dataset directory here
  54. datasetdir = "../datasets/"
  55. trialtype_colors = {
  56. 'SGHF': 'MediumBlue', 'SGLF': 'Turquoise',
  57. 'PGHF': 'DarkGreen', 'PGLF': 'YellowGreen',
  58. 'LFSG': 'Orange', 'LFPG': 'Yellow',
  59. 'HFSG': 'DarkRed', 'HFPG': 'OrangeRed',
  60. 'SGSG': 'SteelBlue', 'PGPG': 'LimeGreen',
  61. 'NONE': 'k', 'PG': 'k', 'SG': 'k', 'LF': 'k', 'HF': 'k'}
  62. event_colors = {
  63. 'TS-ON': 'Gray', # 'TS-OFF': 'Gray',
  64. 'WS-ON': 'Gray', # 'WS-OFF': 'Gray',
  65. 'CUE-ON': 'Gray',
  66. 'CUE-OFF': 'Gray',
  67. 'GO-ON': 'Gray', # 'GO-OFF': 'Gray',
  68. # 'GO/RW-OFF': 'Gray',
  69. 'SR': 'Gray', # 'SR-REP': 'Gray',
  70. 'RW-ON': 'Gray', # 'RW-OFF': 'Gray',
  71. 'STOP': 'Gray'}
  72. # =============================================================================
  73. # Plot helper functions
  74. # =============================================================================
  75. def force_aspect(ax, aspect=1):
  76. ax.set_aspect(abs(
  77. (ax.get_xlim()[1] - ax.get_xlim()[0]) /
  78. (ax.get_ylim()[1] - ax.get_ylim()[0])) / aspect)
  79. def get_arraygrid(signals, chosen_el):
  80. array_grid = np.ones((10, 10)) * 0.7
  81. rejections = np.logical_or(signals.array_annotations['electrode_reject_HFC'],
  82. signals.array_annotations['electrode_reject_LFC'],
  83. signals.array_annotations['electrode_reject_IFC'])
  84. for sig_idx in range(signals.shape[-1]):
  85. connector_aligned_id = signals.array_annotations['connector_aligned_ids'][sig_idx]
  86. x, y = int((connector_aligned_id -1)// 10), int((connector_aligned_id - 1) % 10)
  87. if signals.array_annotations['channel_ids'][sig_idx] == chosen_el:
  88. array_grid[x, y] = -0.7
  89. elif rejections[sig_idx]:
  90. array_grid[x, y] = -0.35
  91. else:
  92. array_grid[x, y] = 0
  93. return np.ma.array(array_grid, mask=np.isnan(array_grid))
  94. # =============================================================================
  95. # Load data and metadata for a monkey
  96. # =============================================================================
  97. # CHANGE this parameter to load data of the different monkeys
  98. monkey = 'Nikos2'
  99. # monkey = 'Lilou'
  100. chosen_el = {'Lilou': 71, 'Nikos2': 63}
  101. chosen_units = {'Lilou': range(1, 5), 'Nikos2': range(1, 5)}
  102. datafile = get_monkey_datafile(monkey)
  103. session = reachgraspio.ReachGraspIO(
  104. filename=os.path.join(datasetdir, datafile),
  105. odml_directory=datasetdir,
  106. verbose=False)
  107. block = session.read_block(lazy=True)
  108. segment = block.segments[0]
  109. # Displaying loaded data structure as string output
  110. print("\nBlock")
  111. print('Attributes ', block.__dict__.keys())
  112. print('Annotations', block.annotations)
  113. print("\nSegment")
  114. print('Attributes ', segment.__dict__.keys())
  115. print('Annotations', segment.annotations)
  116. print("\nEvents")
  117. for x in segment.events:
  118. print('\tEvent with name', x.name)
  119. print('\t\tAttributes ', x.__dict__.keys())
  120. print('\t\tAnnotation keys', x.annotations.keys())
  121. print('\t\ttimes', x.times[:20])
  122. if x.name == 'TrialEvents':
  123. for anno_key in ['trial_id', 'trial_timestamp_id', 'trial_event_labels',
  124. 'trial_reject_IFC']:
  125. print('\t\t'+anno_key, x.array_annotations[anno_key][:20])
  126. print("\nGroups")
  127. for x in block.groups:
  128. print('\tGroup with name', x.name)
  129. print('\t\tAttributes ', x.__dict__.keys())
  130. print('\t\tAnnotations', x.annotations)
  131. print("\nSpikeTrains")
  132. for x in segment.spiketrains:
  133. print('\tSpiketrain with name', x.name)
  134. print('\t\tAttributes ', x.__dict__.keys())
  135. print('\t\tAnnotations', x.annotations)
  136. print('\t\tchannel_id', x.annotations['channel_id'])
  137. print('\t\tunit_id', x.annotations['unit_id'])
  138. print('\t\tis sua', x.annotations['sua'])
  139. print('\t\tis mua', x.annotations['mua'])
  140. print("\nAnalogSignals")
  141. for x in segment.analogsignals:
  142. print('\tAnalogSignal with name', x.name)
  143. print('\t\tAttributes ', x.__dict__.keys())
  144. print('\t\tAnnotations', x.annotations)
  145. print('\t\tchannel_ids', x.array_annotations['channel_ids'])
  146. # get start and stop events of trials
  147. start_events = neo_utils.get_events(
  148. segment,
  149. **{
  150. 'name': 'TrialEvents',
  151. 'trial_event_labels': 'TS-ON',
  152. 'performance_in_trial': 255})
  153. stop_events = neo_utils.get_events(
  154. segment,
  155. **{
  156. 'name': 'TrialEvents',
  157. 'trial_event_labels': 'STOP',
  158. 'performance_in_trial': 255})
  159. # there should only be one event object for these conditions
  160. assert len(start_events) == 1
  161. assert len(stop_events) == 1
  162. # insert epochs between 10ms before TS to 50ms after RW corresponding to trails
  163. ep = neo_utils.add_epoch(
  164. segment,
  165. start_events[0],
  166. stop_events[0],
  167. pre=-250 * pq.ms,
  168. post=500 * pq.ms,
  169. trial_status='complete_trials')
  170. ep.array_annotate(trial_type=start_events[0].array_annotations['belongs_to_trialtype'],
  171. trial_performance=start_events[0].array_annotations['performance_in_trial'])
  172. # access single epoch of this data_segment
  173. epochs = neo_utils.get_epochs(segment, **{'trial_status': 'complete_trials'})
  174. assert len(epochs) == 1
  175. # remove spiketrains not belonging to chosen_electrode
  176. segment.spiketrains = segment.filter(targdict={'channel_id': chosen_el[monkey]},
  177. recursive=True, objects='SpikeTrainProxy')
  178. segment.spiketrains = [st for st in segment.spiketrains if st.annotations['unit_id'] in range(1, 5)]
  179. # replacing the segment with a new segment containing all data
  180. # to speed up cutting of segments
  181. segment = load_segment(segment, load_wavefroms=True, channel_indexes=[chosen_el[monkey]])
  182. # use most raw neuronal data if multiple versions are present
  183. max_sampling_rate = max([a.sampling_rate for a in segment.analogsignals])
  184. idx = 0
  185. while idx < len(segment.analogsignals):
  186. signal = segment.analogsignals[idx]
  187. if signal.annotations['neural_signal'] and signal.sampling_rate < max_sampling_rate:
  188. segment.analogsignals.pop(idx)
  189. else:
  190. idx += 1
  191. # neural_signals = []
  192. # behav_signals = []
  193. # for sig in segment.analogsignals:
  194. # if sig.annotations['neural_signal']:
  195. # neural_signals.append(sig)
  196. # else:
  197. # behav_signals.append(sig)
  198. #
  199. # chosen_raw = neural_signals[0]
  200. # for sig in neural_signals:
  201. # if sig.sampling_rate > chosen_raw.sampling_rate:
  202. # chosen_raw = sig
  203. #
  204. # segment.analogsignals = behav_signals + [chosen_raw]
  205. # cut segments according to inserted 'complete_trials' epochs and reset trial times
  206. cut_segments = neo_utils.cut_segment_by_epoch(segment, epochs[0], reset_time=True)
  207. # =============================================================================
  208. # Define data for overview plots
  209. # =============================================================================
  210. trial_index = {'Lilou': 0, 'Nikos2': 6}
  211. trial_segment = cut_segments[trial_index[monkey]]
  212. blackrock_elid_list = block.annotations['avail_electrode_ids']
  213. # get 'TrialEvents'
  214. event = trial_segment.events[2]
  215. start = np.where(event.array_annotations['trial_event_labels'] == 'TS-ON')[0][0]
  216. trialx_trty = event.array_annotations['belongs_to_trialtype'][start]
  217. trialx_trtimeid = event.array_annotations['trial_timestamp_id'][start]
  218. trialx_color = trialtype_colors[trialx_trty]
  219. # find trial index for next trial with opposite force type (for ax5b plot)
  220. if 'LF' in trialx_trty:
  221. trialz_trty = trialx_trty.replace('LF', 'HF')
  222. else:
  223. trialz_trty = trialx_trty.replace('HF', 'LF')
  224. for i, tr in enumerate(cut_segments):
  225. eventz = tr.events[2]
  226. nextft = np.where(eventz.array_annotations['trial_event_labels'] == 'TS-ON')[0][0]
  227. if eventz.array_annotations['belongs_to_trialtype'][nextft] == trialz_trty:
  228. trialz_trtimeid = eventz.array_annotations['trial_timestamp_id'][nextft]
  229. trialz_color = trialtype_colors[trialz_trty]
  230. trialz_seg = tr
  231. break
  232. # =============================================================================
  233. # Define figure and subplot axis for first data overview
  234. # =============================================================================
  235. fig = plt.figure()
  236. fig.set_size_inches(6.5, 10.) # (w, h) in inches
  237. gs = gridspec.GridSpec(
  238. nrows=5,
  239. ncols=4,
  240. left=0.05,
  241. bottom=0.07,
  242. right=0.9,
  243. top=0.975,
  244. wspace=0.3,
  245. hspace=0.5,
  246. width_ratios=None,
  247. height_ratios=[1, 3, 3, 6, 3])
  248. ax1 = plt.subplot(gs[0, :]) # top row / odml data
  249. # second row
  250. ax2a = plt.subplot(gs[1, 0]) # electrode overview plot
  251. ax2b = plt.subplot(gs[1, 1]) # waveforms unit 1
  252. ax2c = plt.subplot(gs[1, 2]) # waveforms unit 2
  253. ax2d = plt.subplot(gs[1, 3]) # waveforms unit 3
  254. ax3 = plt.subplot(gs[2, :]) # third row / spiketrains
  255. ax4 = plt.subplot(gs[3, :], sharex=ax3) # fourth row / raw signal
  256. ax5a = plt.subplot(gs[4, 0:3]) # fifth row / behavioral signals
  257. ax5b = plt.subplot(gs[4, 3])
  258. fontdict_titles = {'fontsize': 'small', 'fontweight': 'bold'}
  259. fontdict_axis = {'fontsize': 'x-small'}
  260. wf_time_unit = pq.ms
  261. wf_signal_unit = pq.microvolt
  262. plotting_time_unit = pq.s
  263. raw_signal_unit = wf_signal_unit
  264. behav_signal_unit = pq.V
  265. # =============================================================================
  266. # PLOT TRIAL SEQUENCE OF SUBSESSION
  267. # =============================================================================
  268. # load complete metadata collection
  269. odmldoc = odml.load(datasetdir + datafile + '.odml')
  270. # get total trial number
  271. trno_tot = odml_utils.get_TrialCount(odmldoc)
  272. trno_ctr = odml_utils.get_TrialCount(odmldoc, performance_code=255)
  273. trno_ertr = trno_tot - trno_ctr
  274. # get trial id of chosen trial (and next trial with opposite force)
  275. trtimeids = odml_utils.get_TrialIDs(odmldoc, idtype='TrialTimestampID')
  276. trids = odml_utils.get_TrialIDs(odmldoc)
  277. trialx_trid = trids[trtimeids.index(trialx_trtimeid)]
  278. trialz_trid = trids[trtimeids.index(trialz_trtimeid)]
  279. # get all trial ids for grip error trials
  280. trids_pc191 = odml_utils.get_trialids_pc(odmldoc, 191)
  281. # get all trial ids for correct trials
  282. trids_pc255 = odml_utils.get_trialids_pc(odmldoc, 255)
  283. # get occurring trial types
  284. octrty = odml_utils.get_OccurringTrialTypes(odmldoc, code=False)
  285. # Subplot 1: Trial sequence
  286. boxes, labels = [], []
  287. for tt in octrty:
  288. # Plot trial ids of current trial type into trial sequence bar plot
  289. left = odml_utils.get_trialids_trty(odmldoc, tt)
  290. height = np.ones_like(left)
  291. width = 1.
  292. if tt in ['NONE', 'PG', 'SG', 'LF', 'HF']:
  293. color = 'w'
  294. else:
  295. color = trialtype_colors[tt]
  296. B = ax1.bar(
  297. x=left, height=height, width=width, color=color, linewidth=0.001, align='edge')
  298. # Mark trials of current trial type (left) if a grip error occurred
  299. x = [i for i in list(set(left) & set(trids_pc191))]
  300. y = np.ones_like(x) * 2.0
  301. ax1.scatter(x, y, s=5, color='k', marker='*')
  302. # Mark trials of current trial type (left) if any other error occurred
  303. x = [i for i in list(
  304. set(left) - set(trids_pc255) - set(trids_pc191))]
  305. y = np.ones_like(x) * 2.0
  306. ax1.scatter(x, y, s=5, color='gray', marker='*')
  307. # Collect information for trial type legend
  308. if tt not in ['PG', 'SG', 'LF', 'HF']:
  309. boxes.append(B[0])
  310. if tt == 'NONE':
  311. # use errors for providing total trial number
  312. labels.append('total: # %i' % trno_tot)
  313. # add another box and label for error numbers
  314. boxes.append(B[0])
  315. labels.append('* errors: # %i' % trno_ertr)
  316. else:
  317. # trial type trial numbers
  318. labels.append(tt + ': # %i' % len(left))
  319. # mark chosen trial
  320. x = [trialx_trid]
  321. y = np.ones_like(x) * 2.0
  322. ax1.scatter(x, y, s=5, marker='D', color='Red', edgecolors='Red')
  323. # mark next trial with opposite force
  324. x = [trialz_trid]
  325. y = np.ones_like(x) * 2.0
  326. ax1.scatter(x, y, s=5, marker='D', color='orange', edgecolors='orange')
  327. # Generate trial type legend; bbox: (left, bottom, width, height)
  328. leg = ax1.legend(
  329. boxes, labels, bbox_to_anchor=(0., 1., 0.5, 0.1), loc=3, handlelength=1.1,
  330. ncol=len(labels), borderaxespad=0., handletextpad=0.4,
  331. prop={'size': 'xx-small'})
  332. leg.draw_frame(False)
  333. # adjust x and y axis
  334. xticks = list(range(1, 101, 10)) + [100]
  335. ax1.set_xticks(xticks)
  336. ax1.set_xticklabels([str(int(t)) for t in xticks], size='xx-small')
  337. ax1.set_xlabel('trial ID', size='x-small')
  338. ax1.set_xlim(1.-width/2., 100.+width/2.)
  339. ax1.yaxis.set_visible(False)
  340. ax1.set_ylim(0, 3)
  341. ax1.spines['top'].set_visible(False)
  342. ax1.spines['left'].set_visible(False)
  343. ax1.spines['right'].set_visible(False)
  344. ax1.tick_params(direction='out', top=False, left=False, right=False)
  345. ax1.set_title('sequence of the first 100 trials', fontdict_titles, y=2)
  346. ax1.set_aspect('equal')
  347. # =============================================================================
  348. # PLOT ELECTRODE POSITION of chosen electrode
  349. # =============================================================================
  350. neural_signals = [sig for sig in trial_segment.analogsignals if sig.annotations['neural_signal']]
  351. assert len(neural_signals) == 1
  352. neural_signals = neural_signals[0]
  353. arraygrid = get_arraygrid(neural_signals, chosen_el[monkey])
  354. cmap = plt.cm.RdGy
  355. ax2a.pcolormesh(
  356. arraygrid, vmin=-1, vmax=1, lw=1, cmap=cmap, edgecolors='k',
  357. #shading='faceted'
  358. )
  359. force_aspect(ax2a, aspect=1)
  360. ax2a.tick_params(
  361. bottom=False, top=False, left=False, right=False,
  362. labelbottom=False, labeltop=False, labelleft=False, labelright=False)
  363. ax2a.set_title('electrode pos.', fontdict_titles)
  364. # =============================================================================
  365. # PLOT WAVEFORMS of units of the chosen electrode
  366. # =============================================================================
  367. unit_ax_translator = {1: ax2b, 2: ax2c, 3: ax2d}
  368. unit_type = {1: '', 2: '', 3: ''}
  369. wf_lim = []
  370. # plotting waveform for all spiketrains available
  371. for spiketrain in trial_segment.spiketrains:
  372. unit_id = spiketrain.annotations['unit_id']
  373. # get unit type
  374. if spiketrain.annotations['sua']:
  375. unit_type[unit_id] = 'SUA'
  376. elif spiketrain.annotations['mua']:
  377. unit_type[unit_id] = 'MUA'
  378. elif unit_id in [0, 255]:
  379. continue
  380. else:
  381. raise ValueError(f'Found unit with id {unit_id}, that is not SUA or MUA.')
  382. # get correct ax
  383. ax = unit_ax_translator[unit_id]
  384. # get wf sampling time before threshold crossing
  385. left_sweep = spiketrain.left_sweep
  386. # plot waveforms in subplots according to unit id
  387. for st_id, st in enumerate(spiketrain):
  388. wf = spiketrain.waveforms[st_id]
  389. wf_lim.append((np.min(wf), np.max(wf)))
  390. wf_color = str(
  391. (st / spiketrain.t_stop).rescale('dimensionless').magnitude)
  392. times = range(len(wf[0])) * spiketrain.units - left_sweep
  393. ax.plot(
  394. times.rescale(wf_time_unit), wf[0].rescale(wf_signal_unit),
  395. color=wf_color)
  396. ax.set_xlim(
  397. times.rescale(wf_time_unit)[0], times.rescale(wf_time_unit)[-1])
  398. # adding xlabels and titles
  399. for unit_id, ax in unit_ax_translator.items():
  400. ax.set_title('unit %i (%s)' % (unit_id, unit_type[unit_id]),
  401. fontdict_titles)
  402. ax.tick_params(direction='in', length=3, labelsize='xx-small',
  403. labelleft=False, labelright=False)
  404. ax.set_xlabel(wf_time_unit.dimensionality.latex, fontdict_axis)
  405. xticklocator = ticker.MaxNLocator(nbins=5)
  406. ax.xaxis.set_major_locator(xticklocator)
  407. ax.set_ylim(np.min(wf_lim), np.max(wf_lim))
  408. force_aspect(ax, aspect=1)
  409. # adding ylabel
  410. ax2d.tick_params(labelsize='xx-small', labelright=True)
  411. ax2d.set_ylabel(wf_signal_unit.dimensionality.latex, fontdict_axis)
  412. ax2d.yaxis.set_label_position("right")
  413. # =============================================================================
  414. # PLOT SPIKETRAINS of units of chosen electrode
  415. # =============================================================================
  416. plotted_unit_ids = []
  417. # plotting all available spiketrains
  418. for st in trial_segment.spiketrains:
  419. unit_id = st.annotations['unit_id']
  420. plotted_unit_ids.append(unit_id)
  421. ax3.plot(st.times.rescale(plotting_time_unit),
  422. np.zeros(len(st.times)) + unit_id,
  423. 'k|')
  424. # setting layout of spiketrain plot
  425. ax3.set_ylim(min(plotted_unit_ids) - 0.5, max(plotted_unit_ids) + 0.5)
  426. ax3.set_ylabel(r'unit ID', fontdict_axis)
  427. ax3.yaxis.set_major_locator(ticker.MultipleLocator(base=1))
  428. ax3.yaxis.set_label_position("right")
  429. ax3.tick_params(axis='y', direction='in', length=3, labelsize='xx-small',
  430. labelleft=False, labelright=True)
  431. ax3.invert_yaxis()
  432. ax3.set_title('spiketrains', fontdict_titles)
  433. # =============================================================================
  434. # PLOT "raw" SIGNAL of chosen trial of chosen electrode
  435. # =============================================================================
  436. # get "raw" data from chosen electrode
  437. el_raw_sig = [a for a in trial_segment.analogsignals if a.annotations['neural_signal']]
  438. assert len(el_raw_sig) == 1
  439. el_raw_sig = el_raw_sig[0]
  440. # plotting raw signal trace of chosen electrode
  441. chosen_el_idx = np.where(el_raw_sig.array_annotations['channel_ids'] == chosen_el[monkey])[0][0]
  442. ax4.plot(el_raw_sig.times.rescale(plotting_time_unit),
  443. el_raw_sig[:, chosen_el_idx].squeeze().rescale(raw_signal_unit),
  444. color='k')
  445. # setting layout of raw signal plot
  446. ax4.set_ylabel(raw_signal_unit.units.dimensionality.latex, fontdict_axis)
  447. ax4.yaxis.set_label_position("right")
  448. ax4.tick_params(axis='y', direction='in', length=3, labelsize='xx-small',
  449. labelleft=False, labelright=True)
  450. ax4.set_title('"raw" signal', fontdict_titles)
  451. ax4.set_xlim(trial_segment.t_start.rescale(plotting_time_unit),
  452. trial_segment.t_stop.rescale(plotting_time_unit))
  453. ax4.xaxis.set_major_locator(ticker.MultipleLocator(base=1))
  454. # =============================================================================
  455. # PLOT EVENTS across ax3 and ax4 and add time bar
  456. # =============================================================================
  457. # find trial relevant events
  458. startidx = np.where(event.array_annotations['trial_event_labels'] == 'TS-ON')[0][0]
  459. stopidx = np.where(event.array_annotations['trial_event_labels'][startidx:] == 'STOP')[0][0] + startidx + 1
  460. for ax in [ax3, ax4]:
  461. xticks = []
  462. xticklabels = []
  463. for ev_id, ev in enumerate(event[startidx:stopidx]):
  464. ev_labels = event.array_annotations['trial_event_labels'][startidx:stopidx]
  465. if ev_labels[ev_id] in event_colors.keys():
  466. ev_color = event_colors[ev_labels[ev_id]]
  467. ax.axvline(
  468. ev.rescale(plotting_time_unit), color=ev_color, zorder=0.5)
  469. xticks.append(ev.rescale(plotting_time_unit))
  470. if ev_labels[ev_id] == 'CUE-OFF':
  471. xticklabels.append('-OFF')
  472. elif ev_labels[ev_id] == 'GO-ON':
  473. xticklabels.append('GO')
  474. else:
  475. xticklabels.append(ev_labels[ev_id])
  476. ax.set_xticks(xticks)
  477. ax.set_xticklabels(xticklabels)
  478. ax.tick_params(axis='x', direction='out', length=3, labelsize='xx-small',
  479. labeltop=False, top=False)
  480. timebar_ypos = ax4.get_ylim()[0] + np.diff(ax4.get_ylim())[0] / 10
  481. timebar_labeloffset = np.diff(ax4.get_ylim())[0] * 0.01
  482. timebar_xmin = xticks[-2] + ((xticks[-1] - xticks[-2]) / 2 - 0.25 * pq.s)
  483. timebar_xmax = timebar_xmin + 0.5 * pq.s
  484. ax4.plot([timebar_xmin, timebar_xmax], [timebar_ypos, timebar_ypos], '-',
  485. linewidth=3, color='k')
  486. ax4.text(timebar_xmin + 0.25 * pq.s, timebar_ypos + timebar_labeloffset,
  487. '500 ms', ha='center', va='bottom', size='xx-small', color='k')
  488. # =============================================================================
  489. # PLOT BEHAVIORAL SIGNALS of chosen trial
  490. # =============================================================================
  491. # get behavioral signals
  492. ainp_signals = [nsig for nsig in trial_segment.analogsignals if not nsig.annotations['neural_signal']][0]
  493. force_channel_idx = np.where(ainp_signals.array_annotations['channel_ids'] == 141)[0][0]
  494. ainp_trialz_signals = [a for a in trialz_seg.analogsignals if not a.annotations['neural_signal']]
  495. assert len(ainp_trialz_signals)
  496. ainp_trialz = ainp_trialz_signals[0][:, force_channel_idx]
  497. # find out what signal to use
  498. trialx_sec = odmldoc['Recording']['TaskSettings']['Trial_%03i' % trialx_trid]
  499. # get correct channel id
  500. trialx_chids = [143]
  501. FSRi = trialx_sec['AnalogEvents'].properties['UsedForceSensor'].values[0]
  502. FSRinfosec = odmldoc['Setup']['Apparatus']['TargetObject']['FSRSensor']
  503. if 'SG' in trialx_trty:
  504. sgchids = FSRinfosec.properties['SGChannelIDs'].values
  505. trialx_chids.append(min(sgchids) if FSRi == 1 else max(sgchids))
  506. else:
  507. pgchids = FSRinfosec.properties['PGChannelIDs'].values
  508. trialx_chids.append(min(pgchids) if FSRi == 1 else max(pgchids))
  509. # define time epoch
  510. startidx = np.where(event.array_annotations['trial_event_labels'] == 'SR')[0][0]
  511. stopidx = np.where(event.array_annotations['trial_event_labels'] == 'OBB')[0][0]
  512. sr = event[startidx].rescale(plotting_time_unit)
  513. stop = event[stopidx].rescale(plotting_time_unit) + 0.050 * pq.s
  514. startidx = np.where(event.array_annotations['trial_event_labels'] == 'FSRplat-ON')[0][0]
  515. stopidx = np.where(event.array_annotations['trial_event_labels'] == 'FSRplat-OFF')[0][0]
  516. fplon = event[startidx].rescale(plotting_time_unit)
  517. fploff = event[stopidx].rescale(plotting_time_unit)
  518. # define time epoch trialz
  519. startidx = np.where(eventz.array_annotations['trial_event_labels'] == 'FSRplat-ON')[0][0]
  520. stopidx = np.where(eventz.array_annotations['trial_event_labels'] == 'FSRplat-OFF')[0][0]
  521. fplon_trz = eventz[startidx].rescale(plotting_time_unit)
  522. fploff_trz = eventz[stopidx].rescale(plotting_time_unit)
  523. # plotting grip force and object displacement
  524. ai_legend = []
  525. ai_legend_txt = []
  526. for chidx, chid in enumerate(ainp_signals.array_annotations['channel_ids']):
  527. ainp = ainp_signals[:, chidx]
  528. if ainp.array_annotations['channel_ids'][0] in trialx_chids:
  529. ainp_times = ainp.times.rescale(plotting_time_unit)
  530. mask = (ainp_times > sr) & (ainp_times < stop)
  531. ainp_ampli = stats.zscore(ainp.magnitude[mask])
  532. if ainp.array_annotations['channel_ids'][0] != 143:
  533. color = 'gray'
  534. ai_legend_txt.append('grip force')
  535. else:
  536. color = 'k'
  537. ai_legend_txt.append('object disp.')
  538. ai_legend.append(
  539. ax5a.plot(ainp_times[mask], ainp_ampli, color=color)[0])
  540. # get force load of this trial for next plot
  541. elif ainp.array_annotations['channel_ids'][0] == 141:
  542. ainp_times = ainp.times.rescale(plotting_time_unit)
  543. mask = (ainp_times > fplon) & (ainp_times < fploff)
  544. force_av_01 = np.mean(ainp.rescale(behav_signal_unit).magnitude[mask])
  545. # setting layout of grip force and object displacement plot
  546. ax5a.set_title('grip force and object displacement', fontdict_titles)
  547. ax5a.yaxis.set_label_position("left")
  548. ax5a.tick_params(direction='in', length=3, labelsize='xx-small',
  549. labelleft=False, labelright=True)
  550. ax5a.set_ylabel('zscore', fontdict_axis)
  551. ax5a.legend(
  552. ai_legend, ai_legend_txt,
  553. bbox_to_anchor=(0.65, .85, 0.25, 0.1), loc=2, handlelength=1.1,
  554. ncol=len(labels), borderaxespad=0., handletextpad=0.4,
  555. prop={'size': 'xx-small'})
  556. # plotting load/pull force of LF and HF trial
  557. force_times = ainp_trialz.times.rescale(plotting_time_unit)
  558. mask = (force_times > fplon_trz) & (force_times < fploff_trz)
  559. force_av_02 = np.mean(ainp_trialz.rescale(behav_signal_unit).magnitude[mask])
  560. bar_width = [0.4, 0.4]
  561. color = [trialx_color, trialz_color]
  562. ax5b.bar([0, 0.6], [force_av_01, force_av_02], bar_width, color=color)
  563. ax5b.set_title('load/pull force', fontdict_titles)
  564. ax5b.set_ylabel(behav_signal_unit.units.dimensionality.latex, fontdict_axis)
  565. ax5b.set_xticks([0, 0.6])
  566. ax5b.set_xticklabels([trialx_trty, trialz_trty], fontdict=fontdict_axis)
  567. ax5b.yaxis.set_label_position("right")
  568. ax5b.tick_params(direction='in', length=3, labelsize='xx-small',
  569. labelleft=False, labelright=True)
  570. # =============================================================================
  571. # PLOT EVENTS across ax5a and add time bar
  572. # =============================================================================
  573. # find trial relevant events
  574. startidx = np.where(event.array_annotations['trial_event_labels'] == 'SR')[0][0]
  575. stopidx = np.where(event.array_annotations['trial_event_labels'] == 'OBB')[0][0]
  576. xticks = []
  577. xticklabels = []
  578. for ev_id, ev in enumerate(event[startidx:stopidx]):
  579. ev_labels = event.array_annotations['trial_event_labels'][startidx:stopidx + 1]
  580. if ev_labels[ev_id] in ['RW-ON']:
  581. ax5a.axvline(ev.rescale(plotting_time_unit), color='k', zorder=0.5)
  582. xticks.append(ev.rescale(plotting_time_unit))
  583. xticklabels.append(ev_labels[ev_id])
  584. elif ev_labels[ev_id] in ['OT', 'OR', 'DO', 'OBB', 'FSRplat-ON',
  585. 'FSRplat-OFF', 'HEplat-ON']:
  586. ev_color = 'k'
  587. xticks.append(ev.rescale(plotting_time_unit))
  588. xticklabels.append(ev_labels[ev_id])
  589. ax5a.axvline(
  590. ev.rescale(plotting_time_unit), color='k', ls='-.', zorder=0.5)
  591. elif ev_labels[ev_id] == 'HEplat-OFF':
  592. ev_color = 'k'
  593. ax5a.axvline(
  594. ev.rescale(plotting_time_unit), color='k', ls='-.', zorder=0.5)
  595. ax5a.set_xticks(xticks)
  596. ax5a.set_xticklabels(xticklabels, fontdict=fontdict_axis, rotation=90)
  597. ax5a.tick_params(axis='x', direction='out', length=3, labelsize='xx-small',
  598. labeltop=False, top=False)
  599. ax5a.set_ylim([-2.0, 2.0])
  600. timebar_ypos = ax5a.get_ylim()[0] + np.diff(ax5a.get_ylim())[0] / 10
  601. timebar_labeloffset = np.diff(ax5a.get_ylim())[0] * 0.02
  602. timebar_xmax = xticks[xticklabels.index('RW-ON')] - 0.1 * pq.s
  603. timebar_xmin = timebar_xmax - 0.25 * pq.s
  604. ax5a.plot([timebar_xmin, timebar_xmax], [timebar_ypos, timebar_ypos], '-',
  605. linewidth=3, color='k')
  606. ax5a.text(timebar_xmin + 0.125 * pq.s, timebar_ypos + timebar_labeloffset,
  607. '250 ms', ha='center', va='bottom', size='xx-small', color='k')
  608. # add time window of ax5a to ax4
  609. ax4.axvspan(ax5a.get_xlim()[0], ax5a.get_xlim()[1], facecolor=[0.9, 0.9, 0.9],
  610. zorder=-0.1, ec=None)
  611. # =============================================================================
  612. # SAVE FIGURE
  613. # =============================================================================
  614. fname = 'data_overview_1_%s' % monkey
  615. for file_format in ['eps', 'png', 'pdf']:
  616. fig.savefig(fname + '.%s' % file_format, dpi=400, format=file_format)