data_overview_1.py 27 KB

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