notebook.community

Edit and run 





In [1]:



    


import numpy as np
import scipy as sp
import seaborn as sn
import matplotlib.pylab as pl
%matplotlib inline

import os 
os.chdir('../src/')

from EDFOperator import EDFOperator
from HDFEyeOperator import HDFEyeOperator
from EyeSignalOperator import EyeSignalOperator

os.chdir('../test/')

sn.set(style="ticks")



    











In [2]:



    


edf_file = 'tk_1_2014-10-06_16.44.11.edf'

low_pass_pupil_f, high_pass_pupil_f = 6.0, 0.01
alias = 'test_1'

# initialize the hdfeyeoperator
ho = HDFEyeOperator(os.path.expanduser('~/Downloads/test.h5'))
# insert the edf file contents only when the h5 is not present.
if not os.path.isfile(os.path.expanduser('~/Downloads/test.h5')):
    ho.add_edf_file(edf_file)
    ho.edf_message_data_to_hdf(alias = alias)
    ho.edf_gaze_data_to_hdf(alias = alias, pupil_hp = high_pass_pupil_f, pupil_lp = low_pass_pupil_f)



    


















    






16-10-12_11-45-05 - INFO - EDFOperator - started with tk_1_2014-10-06_16.44.11.edf
16-10-12_11-45-05 - INFO - EDF2ASCOperator - <CommandLineOperator.EDF2ASCOperator object at 0x103e59210> initialized with file tk_1_2014-10-06_16.44.11.edf
16-10-12_11-45-05 - INFO - EDFOperator - reading trials from tk_1_2014-10-06_16.44.11.msg
[['motion_direction', 'cycle_speed', 'grating_background_lum', 'stim_eye_correspondence', 'which_eye_stimulated', 'spatial_freq', 'mask_radius', 'rotation_speed', 'red_green_ratio', 'trial_nr'], ['motion_direction', 'cycle_speed', 'grating_background_lum', 'stim_eye_correspondence', 'which_eye_stimulated', 'spatial_freq', 'mask_radius', 'rotation_speed', 'red_green_ratio', 'trial_nr'], ['motion_direction', 'cycle_speed', 'grating_background_lum', 'stim_eye_correspondence', 'which_eye_stimulated', 'spatial_freq', 'mask_radius', 'rotation_speed', 'red_green_ratio', 'trial_nr'], ['motion_direction', 'cycle_speed', 'grating_background_lum', 'stim_eye_correspondence', 'which_eye_stimulated', 'spatial_freq', 'mask_radius', 'rotation_speed', 'red_green_ratio', 'trial_nr']]
16-10-12_11-45-05 - INFO - EDFOperator - reading key_events from tk_1_2014-10-06_16.44.11.msg
16-10-12_11-45-05 - INFO - EDFOperator - reading eyelink events from tk_1_2014-10-06_16.44.11.msg
16-10-12_11-45-05 - INFO - EDFOperator - reading sounds from tk_1_2014-10-06_16.44.11.msg
16-10-12_11-45-05 - INFO - EDFOperator - cleaning gaze information from tk_1_2014-10-06_16.44.11.gaz
16-10-12_11-45-05 - ERROR - EDFOperator - gaze information already cleaned
16-10-12_11-45-05 - INFO - EDFOperator - identifying recording blocks from tk_1_2014-10-06_16.44.11.gaz
16-10-12_11-45-08 - INFO - EDFOperator - 1 raw blocks discovered
16-10-12_11-45-08 - INFO - EDFOperator - 695256.000000 raw block duration, threshold 30000.000000
16-10-12_11-45-08 - INFO - EDFOperator - 1 correct duration blocks discovered
16-10-12_11-45-11 - INFO - EDFOperator - found data from block 0 of shape (695254, 6)
16-10-12_11-45-11 - INFO - HDFEyeOperator - Adding message data from tk_1_2014-10-06_16.44.11.edf to group  test_1 to /Users/knapen/Downloads/test.h5
16-10-12_11-45-11 - INFO - HDFEyeOperator - Adding gaze data from tk_1_2014-10-06_16.44.11.edf to group  test_1 to /Users/knapen/Downloads/test.h5
16-10-12_11-45-11 - INFO - EyeSignalOperator - Interpolating blinks using interpolate_blinks
16-10-12_11-45-12 - INFO - EyeSignalOperator - Band-pass filtering of pupil signals, hp = 0.010, lp = 6.000
16-10-12_11-45-12 - INFO - EyeSignalOperator - Regressing blinks, saccades and gaze position of pupil signals










    






EyeSignalOperator.py:407: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
EyeSignalOperator.py:435: VisibleDeprecationWarning: boolean index did not match indexed array along dimension 0; dimension is 695254 but corresponding boolean dimension is 695253
EyeSignalOperator.py:478: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
EyeSignalOperator.py:500: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.










    






[(695254,), (695254,), (695254,), (695254,), (695254,)]
16-10-12_11-45-12 - INFO - EyeSignalOperator - Nuisance GLM Results, pearsons R (p) is 0.409 (0.0000)
16-10-12_11-45-12 - INFO - EyeSignalOperator - Nuisance GLM Results, Blink, Saccade beta 504.730, -150.869)
16-10-12_11-45-12 - INFO - EyeSignalOperator - Nuisance GLM Results, Gaze x, y, intercept beta -3.706, 3.052, -2.642 )
16-10-12_11-45-14 - INFO - EyeSignalOperator - Time_frequency_decomposition_pupil, with filterbank [ 0.1         0.06305834  0.03976354  0.02507422  0.01581139  0.0099704
  0.00628717  0.00396458  0.0025    ]
16-10-12_11-45-14 - INFO - HDFEyeOperator - Performed T-F analysis of type lp_butterworth
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.10000
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.06306
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.03976
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.02507
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.01581
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.00997
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.00629
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.00396
16-10-12_11-45-14 - INFO - HDFEyeOperator - Saved T-F analysis 0.00250










    

















    
























In [4]:



    


# downsample for plotting
downsample_rate = 10

# load times per session:
trial_times = ho.read_session_data(alias, 'trials')
trial_phase_times = ho.read_session_data(alias, 'trial_phases')

# check at what timestamps the recording started:
session_start_EL_time = np.array( trial_phase_times[np.array(trial_phase_times['trial_phase_index'] == 1) * np.array(trial_phase_times['trial_phase_trial'] == 0)]['trial_phase_EL_timestamp'] )[0]
session_stop_EL_time = np.array(trial_times['trial_end_EL_timestamp'])[-1]

# and, find some aspects of the recording such as sample rate and recorded eye
sample_rate = ho.sample_rate_during_period([session_start_EL_time, session_stop_EL_time], alias)
eye = ho.eye_during_period([session_start_EL_time, session_stop_EL_time], alias)
if len(eye) > 0:
    eye = ['L','R'][0]



    











In [5]:



    


#
# plot some of the pupil signal's processing steps
#

pupil_raw = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil', requested_eye = eye))
pupil_int = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil_int', requested_eye = eye))

pupil_bp = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil_bp', requested_eye = eye))
pupil_lp = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil_lp', requested_eye = eye))
pupil_hp = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil_hp', requested_eye = eye))

x = sp.signal.decimate(np.arange(len(pupil_raw)) / float(sample_rate), downsample_rate, 1)
pup_raw_dec = sp.signal.decimate(pupil_raw, downsample_rate, 1)
pup_int_dec = sp.signal.decimate(pupil_int, downsample_rate, 1)

pupil_bp_dec = sp.signal.decimate(pupil_bp, downsample_rate, 1)
pupil_lp_dec = sp.signal.decimate(pupil_lp, downsample_rate, 1)
pupil_hp_dec = sp.signal.decimate(pupil_hp, downsample_rate, 1)

# plot interpolated pupil:
fig = pl.figure(figsize = (10,9))
s = fig.add_subplot(311)
pl.plot(x, pup_raw_dec, 'b'); pl.plot(x, pup_int_dec, 'g')
pl.ylabel('pupil size'); pl.xlabel('time (s)')
pl.legend(['raw pupil', 'blink interpolated pupil'])
s.set_title('tk')

ymin = pupil_raw.min(); ymax = pupil_raw.max()
tps = (list(trial_phase_times[trial_phase_times['trial_phase_index'] == 2]['trial_phase_EL_timestamp']) - session_start_EL_time, list(trial_phase_times[trial_phase_times['trial_phase_index'] == 3]['trial_phase_EL_timestamp']) - session_start_EL_time)
for i in range(tps[0].shape[0]):
    pl.axvline(x = tps[0][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'r')
    pl.axvline(x = tps[1][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'k')
s.set_ylim(ymin = pup_int_dec.min()-100, ymax = pup_int_dec.max()+100)
s.set_xlim(xmin = tps[0][0] / float(sample_rate), xmax = tps[1][-1] / float(sample_rate))
sn.despine(offset=10)

s = fig.add_subplot(312)
pl.plot(x, pupil_bp_dec, 'b'); pl.plot(x, pupil_lp_dec, 'g');
pl.ylabel('pupil size'); pl.xlabel('time (s)')
pl.legend(['band_passed', 'lowpass'])
s.set_title('tk')

ymin = pupil_raw.min(); ymax = pupil_raw.max()
tps = (list(trial_phase_times[trial_phase_times['trial_phase_index'] == 2]['trial_phase_EL_timestamp']) - session_start_EL_time, list(trial_phase_times[trial_phase_times['trial_phase_index'] == 3]['trial_phase_EL_timestamp']) - session_start_EL_time)
for i in range(tps[0].shape[0]):
    pl.axvline(x = tps[0][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'r')
    pl.axvline(x = tps[1][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'k')
# s.set_ylim(ymin = pup_int_dec.min()-100, ymax = pup_int_dec.max()+100)
s.set_xlim(xmin = tps[0][0] / float(sample_rate), xmax = tps[1][-1] / float(sample_rate))
sn.despine(offset=10)

s = fig.add_subplot(313)
pl.plot(x, pupil_bp_dec, 'b'); pl.plot(x, pupil_hp_dec, 'g');
pl.ylabel('pupil size'); pl.xlabel('time (s)')
pl.legend(['band_passed', 'highpass'])
s.set_title('tk')

ymin = pupil_raw.min(); ymax = pupil_raw.max()
tps = (list(trial_phase_times[trial_phase_times['trial_phase_index'] == 2]['trial_phase_EL_timestamp']) - session_start_EL_time, list(trial_phase_times[trial_phase_times['trial_phase_index'] == 3]['trial_phase_EL_timestamp']) - session_start_EL_time)
for i in range(tps[0].shape[0]):
    pl.axvline(x = tps[0][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'r')
    pl.axvline(x = tps[1][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'k')
# s.set_ylim(ymin = pup_int_dec.min()-100, ymax = pup_int_dec.max()+100)
s.set_xlim(xmin = tps[0][0] / float(sample_rate), xmax = tps[1][-1] / float(sample_rate))
sn.despine(offset=10)
pl.tight_layout()



    


















    






/Users/knapen/anaconda3/envs/i2.7/lib/python2.7/site-packages/scipy/signal/signaltools.py:3073: FutureWarning:  Note: Decimate's zero_phase keyword argument will default to True in a future release. Until then, decimate defaults to one-way filtering for backwards compatibility. Ideally, always set this argument explicitly.
  "explicitly.", FutureWarning)










    
























In [6]:



    


#
# plot the 'regress blinks' processing step, 
# which estimates the effects of blinks, saccades and foreshortening on the data
# and regresses these effects out.
#

pupil_bp = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil_bp_zscore', requested_eye = eye))
pupil_bp_c = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil_bp_clean_zscore', requested_eye = eye))

x = sp.signal.decimate(np.arange(len(pupil_bp)) / float(sample_rate), downsample_rate, 1)
pupil_bp_dec = sp.signal.decimate(pupil_bp, downsample_rate, 1)
pupil_bp_c_dec = sp.signal.decimate(pupil_bp_c, downsample_rate, 1)


# plot interpolated pupil:
fig = pl.figure(figsize = (10,9))
for trial in range(4):# loop across trials
    s = fig.add_subplot(4,1,trial+1)
    s.set_title('trial %i'%trial)
    pl.plot(x, pupil_bp_dec, 'k', lw = 2.5)
    pl.plot(x, pupil_bp_c_dec, 'r', lw = 1.5)
    pl.ylabel('pupil size'); pl.xlabel('time (s)')
    pl.legend(['uncleaned pupil', 'blink, saccade and foreshortening regressed pupil'])

    ymin = pupil_bp.min(); ymax = pupil_bp.max()
    tps = (list(trial_phase_times[trial_phase_times['trial_phase_index'] == 2]['trial_phase_EL_timestamp']) - session_start_EL_time, list(trial_phase_times[trial_phase_times['trial_phase_index'] == 3]['trial_phase_EL_timestamp']) - session_start_EL_time)
    for i in range(tps[0].shape[0]):
        pl.axvline(x = tps[0][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'g')
        pl.axvline(x = tps[1][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'k')
    # s.set_ylim(ymin = pup_int_dec.min()-100, ymax = pup_int_dec.max()+100)
    s.set_xlim(xmin = tps[0][trial] / float(sample_rate), xmax = tps[-1][trial] / float(sample_rate))
    sn.despine(offset=10)
pl.tight_layout()



    


















    
























In [7]:



    


#
# plot the 'regress blinks' processing step, 
# which estimates the effects of blinks, saccades and foreshortening on the data
# and regresses these effects out.
#

pupil_lp = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil_lp', requested_eye = eye))
pupil_lp_c = np.squeeze(ho.signal_during_period(time_period = [session_start_EL_time, session_stop_EL_time], alias = alias, signal = 'pupil_lp_clean', requested_eye = eye))

x = sp.signal.decimate(np.arange(len(pupil_lp)) / float(sample_rate), downsample_rate, 1)
pupil_lp_dec = sp.signal.decimate(pupil_lp, downsample_rate, 1)
pupil_lp_c_dec = sp.signal.decimate(pupil_lp_c, downsample_rate, 1)


# plot interpolated pupil:
fig = pl.figure(figsize = (10,9))
for trial in range(4):# loop across trials
    s = fig.add_subplot(4,1,trial+1)
    s.set_title('trial %i'%trial)
    pl.plot(x, pupil_lp_dec, 'k', lw = 2.5)
    pl.plot(x, pupil_lp_c_dec, 'r', lw = 1.5)
    pl.ylabel('pupil size'); pl.xlabel('time (s)')
    pl.legend(['uncleaned pupil', 'blink, saccade and foreshortening regressed pupil'], loc= 4)

    ymin = pupil_bp.min(); ymax = pupil_bp.max()
    tps = (list(trial_phase_times[trial_phase_times['trial_phase_index'] == 2]['trial_phase_EL_timestamp']) - session_start_EL_time, list(trial_phase_times[trial_phase_times['trial_phase_index'] == 3]['trial_phase_EL_timestamp']) - session_start_EL_time)
    for i in range(tps[0].shape[0]):
        pl.axvline(x = tps[0][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'g')
        pl.axvline(x = tps[1][i] / float(sample_rate), ymin = ymin, ymax = ymax, color = 'k')
    s.set_ylim(ymin = 5000, ymax = 9000)
    s.set_xlim(xmin = tps[0][trial] / float(sample_rate), xmax = tps[-1][trial] / float(sample_rate))
    sn.despine(offset=10)
pl.tight_layout()



    


















    
























In [ ]:

Content source: tknapen/hedfpy

Similar notebooks:

notebook.community | gallery | about