API documentation

pupildata.py

Main object-oriented entry point

class pypillometry.pupildata.FakePupilData(pupil, sampling_rate=None, time=None, event_onsets=None, event_labels=None, name=None, sim_params={}, real_baseline=None, real_response_coef=None)

Simulated pupil data for validation purposes.

__init__(pupil, sampling_rate=None, time=None, event_onsets=None, event_labels=None, name=None, sim_params={}, real_baseline=None, real_response_coef=None)

Constructor for artifical pupil data.

plot(plot_range=- inf, inf, interactive=False, baseline=True, response=False, model=True, simulated=True, units='sec')

Make a plot of the pupil data using matplotlib or pypillometry.convenience.plot_pupil_ipy() if interactive=True.

Parameters

  • plot_range (tuple (start,end): plot from start to end (in units of units)) –

  • baseline (plot baseline if estimated) –

  • response (plot response if estimated) –

  • model (plot full model if baseline and response have been estimated) –

  • simulated (plot also the "ground-truth" baseline and response (i.e., the simulated one)?) –

  • interactive (if True, plot with sliders to adjust range) –

  • units (one of "sec"=seconds, "ms"=millisec, "min"=minutes, "h"=hours) –

Return type

None

scale(mean=None, sd=None, inplace=False)

Scale the pupillary signal by subtracting mean and dividing by sd. If these variables are not provided, use the signal’s mean and std.

Parameters

  • mean (mean to subtract from signal) –

  • sd (sd to scale with) –

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

Note

Scaling-parameters are being saved in the scale_params argument.

Return type

None

sub_slice(start=- inf, end=inf, units='sec')

Return a new PupilData object that is a shortened version of the current one (contains all data between start and end in units given by units (one of “ms”, “sec”, “min”, “h”).

unscale(mean=None, sd=None, inplace=False)

Scale back to original values using either values provided as arguments or the values stored in scale_params.

Parameters

  • mean (mean to add from signal) –

  • sd (sd to scale with) –

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

class pypillometry.pupildata.PupilData(pupil, sampling_rate=None, time=None, event_onsets=None, event_labels=None, name=None, keep_orig=True, fill_time_discontinuities=True)

Class representing pupillometric data.

__init__(pupil, sampling_rate=None, time=None, event_onsets=None, event_labels=None, name=None, keep_orig=True, fill_time_discontinuities=True)
Parameters

  • name (Optional[str]) – name of the dataset or None (in which case a random string is selected)

  • time (Union[ndarray, List[float], None]) – timing array or None, in which case the time-array goes from [0,maxT] using sampling_rate (in ms)

  • pupil (Union[ndarray, List[float]]) – pupillary data at times time assumed to be in ms

  • event_onsets (Union[ndarray, List[float], None]) – time-onsets of any events that are to be modelled in the pupil

  • event_labels (Union[ndarray, List[float], None]) – for each event in event_onsets, provide a label

  • sampling_rate (float) – sampling-rate of the pupillary signal in Hz

  • keep_orig (bool) – keep a copy of the original dataset? If True, a copy of the PupilData object as initiated in the constructor is stored in member PupilData.original

  • fill_time_discontinuities (bool) – sometimes, when the eyetracker loses signal, no entry in the EDF is made; when this option is True, such entries will be made and the signal set to 0 there

__len__()

Return number of sampling points in the pupil data.

Return type

int

__repr__()

Return a string-representation of the dataset.

Return type

str

__weakref__

list of weak references to the object (if defined)

add_to_history(event)

Add event to history

apply_history(obj)

Apply history of operations done on self to obj.

obj: PupilData

object of class PupilData to which the operations are to be transferred

copy of the PupilData-object to which the operations in self were applied

Detect blinks in the pupillary signal using several strategies. First, blinks are detected as consecutive sequence of blink_val (f.eks., 0 or NaN). Second, blinks are defined as everything between two crossings of the velocity profile (from negative to positive).

Detected blinks are put into member blinks (matrix 2 x nblinks where start and end are stored as indexes) and member blink_mask which codes for each sampling point whether there is a blink (1) or not (0).

Finally, detected blinks have to be at least min_duration duration (in units).

Parameters

  • min_duration (float) – minimum duration for a sequence of missing numbers to be treated as blink

  • blink_val (float) – “missing value” code

  • winsize (float) – window-size for smoothing for velocity profile (in units)

  • vel_onset (float) – negative velocity that needs to be crossed; arbitrary units that depend on sampling rate etc

  • vel_offset (float) – positive velocity that needs to be exceeded; arbitrary units that depend on sampling rate etc

  • min_onset_len (int) – minimum number of consecutive samples that crossed threshold in the velocity profile to detect as onset (to avoid noise-induced changes)

  • min_offset_len (int) – minimum number of consecutive samples that crossed threshold in the velocity profile to detect as offset (to avoid noise-induced changes)

  • strategies (list of strategies to use) – so far, use a list containing any combination of “zero” and “velocity”

  • units (str) – one of “ms”, “sec”, “min”, “h”

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

Implements the blink-interpolation method by Mahot (2013).

Mahot, 2013: https://figshare.com/articles/A_simple_way_to_reconstruct_pupil_size_during_eye_blinks/688001.

This procedure relies heavily on eye-balling (reconstructing visually convincing signal), so a “plot” option is provided that will plot many diagnostics (see paper linked above) that can help to set good parameter values for winsize, vel_onset, vel_offset and margin.

Parameters

  • winsize (float) – size of the Hanning-window in ms

  • vel_onset (float) – velocity-threshold to detect the onset of the blink

  • vel_offset (float) – velocity-threshold to detect the offset of the blink

  • margin (Tuple[float,float]) – margin that is subtracted/added to onset and offset (in ms)

  • blinkwindow (float) – how much time before and after each blink to include (in ms)

  • interp_type (str) – type of interpolation accepted by scipy.interpolate.interp1d()

  • plot (True, str or None) – if a string, the plot is going to be saved to a multipage PDF file; if None, no plotting is done if True, plot is not saved but produced

  • plot_dim (tuple nrow x ncol) – number of subplots

  • plot_figsize (tuple (width, height)) – dimensions for each figure

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

Interpolation of missing data “in one go”. Detection of blinks happens using Mahot (2013), see blink_onsets_mahot().

Parameters

  • winsize (float) – size of the Hanning-window in ms

  • vel_onset (float) – velocity-threshold to detect the onset of the blink

  • vel_offset (float) – velocity-threshold to detect the offset of the blink

  • margin (Tuple[float,float]) – margin that is subtracted/added to onset and offset (in ms)

  • interp_type (str) – type of interpolation accepted by scipy.interpolate.interp1d()

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

Merge together blinks that are close together. Some subjects blink repeatedly and standard detection/interpolation can result in weird results. This function simply treats repeated blinks as one long blink.

Parameters

  • distance (float) – merge together blinks that are closer together than distance in ms

  • remove_signal (bool) – if True, set all signal values during the “new blinks” to zero so that detect_blinks() will pick them up; interpolation will work either way

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

Plot the detected blinks into separate figures each with nrow x ncol subplots.

Parameters

  • pdf_file (str or None) – if the name of a file is given, the figures are saved into a multi-page PDF file

  • ncol (int) – number of columns for the blink-plots

  • pre_blink (float) – extend plot a certain time before each blink (in ms)

  • post_blink (float) – extend plot a certain time after each blink (in ms)

  • units (str) – units in which the signal is plotted

  • plot_index (bool) – plot a number with the blinks’ index (e.g., for identifying abnormal blinks)

Returns

  • list of plt.Figure objects each with nrow*ncol subplots

  • in Jupyter Notebook, those are displayed inline one after the other

copy(new_name=None)

Make and return a deep-copy of the pupil data.

downsample(fsd, dsfac=False, inplace=False)

Simple downsampling scheme using mean within the downsampling window. See pypillometry.baseline.downsample().

Parameters

  • fsd (float) – new sampling-rate or decimate-factor

  • dsfac (bool) – if False, fsd is the new sampling rate; if True, fsd is the decimate factor

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

drop_original(**kwargs)

Drop original dataset from record (e.g., to save space).

estimate_baseline(method='envelope_iter_bspline_2', inplace=False, **kwargs)

Apply one of the baseline-estimation methods.

Parameters

Note

the results of the estimation is stored in member baseline

estimate_response(npar='free', tmax='free', verbose=50, bounds={'npar': 1, 20, 'tmax': 100, 2000}, inplace=False)

Estimate pupil-response based on event-onsets, see pypillometry.pupil.pupil_response().

npar: float

npar-parameter for the canonical response-function or “free”; in case of “free”, the function optimizes for this parameter

tmax: float

tmax-parameter for the canonical response-function or “free”; in case of “free”, the function optimizes for this parameter

bounds: dict

in case that one or both parameters are estimated, give the lower and upper bounds for the parameters

inplace: bool

if True, make change in-place and return the object if False, make and return copy before making changes

Note

the results of the estimation is stored in members response, response_x (design matrix) and response_pars

classmethod from_file(fname)

Reads a PupilData object from a pickle-file. Use as pypillometry.PupilData.from_file("yourfile.pd").

Parameters

fname (str) – filename

get_erpd(erpd_name, event_select, baseline_win=None, time_win=- 500, 2000)

Extract event-related pupil dilation (ERPD). No attempt is being made to exclude overlaps of the time-windows.

Parameters

  • erpd_name (str) – identifier for the result (e.g., “cue-locked” or “conflict-trials”)

  • baseline_win (tuple (float,float) or None) –

    if None, no baseline-correction is applied if tuple, the mean value in the window in milliseconds (relative to time_win) is

    subtracted from the single-trial ERPDs (baseline-correction)

  • event_select (str or function) – variable describing which events to select and align to - if str: use all events whose label contains the string - if function: apply function to all labels, use those where the function returns True

  • time_win (Tuple[float, float]) – time before and after event to include (in ms)

lowpass_filter(cutoff, order=2, inplace=False)

Lowpass-filter signal using a Butterworth-filter, see pypillometry.baseline.butter_lowpass_filter().

Parameters

  • cutoff (float) – lowpass-filter cutoff

  • order (int) – filter order

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

Return number of detected blinks. Should be run after detect_blinks().

Return type

int

nevents()

Return number of events in pupillometric data.

Return type

int

plot(plot_range=- inf, inf, interactive=False, baseline=True, response=False, model=True, highlight_blinks=True, highlight_interpolated=True, units='sec')

Make a plot of the pupil data using matplotlib or pypillometry.convenience.plot_pupil_ipy() if interactive=True.

Parameters

  • plot_range (tuple (start,end)) – plot from start to end (in units of units)

  • baseline (bool) – plot baseline if estimated

  • response (bool) – plot response if estimated

  • model (bool) – plot full model if baseline and response have been estimated

  • interactive (bool) – if True, plot with sliders to adjust range

  • units (str) – one of “sec”=seconds, “ms”=millisec, “min”=minutes, “h”=hours

Return type

None

plot_segments(overlay=None, pdffile=None, interv=1, figsize=15, 5, ylim=None, **kwargs)

Plot the whole dataset chunked up into segments (usually to a PDF file).

Parameters

  • pdffile (str or None) – file name to store the PDF; if None, no PDF is written

  • interv (float) – duration of each of the segments to be plotted (in minutes)

  • figsize (Tuple[int,int]) – dimensions of the figures

  • kwargs – arguments passed to PupilData.plot()

Returns

figs

Return type

list of matplotlib.Figure objects

print_history()

Pretty-print the history of the current dataset (which manipulations where done on it).

reset_time(t0=0, inplace=False)

Resets time so that the time-array starts at time zero (t0). Resets onsets etc.

Parameters

  • t0 (float) – time at which the PupilData’s time-vector starts

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

scale(mean=None, sd=None, inplace=False)

Scale the pupillary signal by subtracting mean and dividing by sd. If these variables are not provided, use the signal’s mean and std.

Parameters

  • mean (float) – mean to subtract from signal

  • sd (float) – sd to scale with

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

Note

Scaling-parameters are being saved in the scale_params argument.

size_bytes()

Return size of current dataset in bytes.

smooth_window(window='hanning', winsize=11, inplace=False)

Apply smoothing of the signal using a moving window. See smooth_window().

Parameters

  • window (str) –

    (the type of window from ‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’);

    flat window will produce a moving average smoothing.

  • winsize (float) – the length of the window in ms

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

stat_per_event(interval, event_select=None, statfct=<function mean>, return_missing=None)

Return result of applying a statistical function to pupillometric data in a given interval relative to event-onsets. For example, extract mean pupil-size in interval before trial onset.

Parameters

  • event_select (str or function) – variable describing which events to select and align to - if str: use all events whose label contains the string - if function: apply function to all labels, use those where the function returns True

  • interval (tuple (min,max)) – time-window in ms relative to event-onset (0 is event-onset)

  • statfct (function) – function mapping np.array to a single number

  • return_missing (None, "nmiss", "prop") – if None, only an array with the stats per event is return if “nmiss”, returns a tuple (stat, nmiss) where nmiss is the number of missing vales in the timewin if “prop”, return a tuple (stat, prop_miss) where prop_miss is the proportion missing vales in the timewin

Returns

result – number of event-onsets long result array

Return type

np.array

sub_slice(start=- inf, end=inf, units='sec')

Return a new PupilData object that is a shortened version of the current one (contains all data between start and end in units given by units (one of “ms”, “sec”, “min”, “h”).

Parameters

  • start (float) – start for new dataset

  • end (float) – end of new dataset

  • units (str) – time units in which start and end are provided

summary()

Return a summary of the PupilData-object.

Return type

dict

unscale(mean=None, sd=None, inplace=False)

Scale back to original values using either values provided as arguments or the values stored in scale_params.

Parameters

  • mean (float) – mean to add from signal

  • sd (float) – sd to scale with

  • inplace (bool) – if True, make change in-place and return the object if False, make and return copy before making changes

write_file(fname)

Save to file (using pickle).

Parameters

fname (str) – filename

pypillometry.pupildata.create_fake_pupildata(**kwargs)

Return a pyillometry.pupildata.FakePupilData object by buildling it with pypillometry.fakedata.get_dataset().

Parameters

  • ntrials (int) – number of trials

  • isi (float) – inter-stimulus interval in seconds

  • rtdist (tuple (float,float)) – mean and std of a (truncated at zero) normal distribution to generate response times

  • pad (float) – padding before the first and after the last event in seconds

  • fs (float) – sampling rate in Hz

  • baseline_lowpass (float) – cutoff for the lowpass-filter that defines the baseline (highest allowed frequency in the baseline fluctuations)

  • evoked_response_perc (float) – amplitude of the pupil-response as proportion of the baseline

  • response_fluct_sd (float) – How much do the amplitudes of the individual events fluctuate? This is determined by drawing each individual pupil-response to a single event from a (positive) normal distribution with mean as determined by evoked_response_perc and sd response_fluct_sd (in units of evoked_response_perc).

  • prf_npar (tuple (float,float)) – (mean,std) of the npar parameter from pypillometry.pupil.pupil_kernel(). If the std is exactly zero, then the mean is used for all pupil-responses. If the std is positive, npar is taken i.i.d. from ~ normal(mean,std) for each event.

  • prf_tmax (tuple (float,float)) – (mean,std) of the tmax parameter from pypillometry.pupil.pupil_kernel(). If the std is exactly zero, then the mean is used for all pupil-responses. If the std is positive, tmax is taken i.i.d. from ~ normal(mean,std) for each event.

  • prop_spurious_events (float) – Add random events to the pupil signal. prop_spurious_events is expressed as proportion of the number of real events.

  • noise_amp (float) – Amplitude of random gaussian noise that sits on top of the simulated signal. Expressed in units of mean baseline pupil diameter.

pypillometry.pupildata.plotpd(*args, subplots=False, baseline=False)

Plotting for PupilData objects.

Parameters

subplots (plot the different PupilData-objects in the same plot or subplots) –

pypillometry.pupildata.plotpd_ia(*args, figsize=16, 8, baseline=True, events=0)

Interactive plotting for multiple PupilData objects.

Parameters

  • args (PupilData datasets to plot) –

  • figsize (dimensions of the plot) –

  • baseline (plot baselines, too?) –

  • events (plot event-markers? if None, no events are plotted, otherwise events) – is the index of the PupilData object to take the events from

erpd.py

Event-related pupil dilation.

class pypillometry.erpd.ERPD(name, tx, erpd, missing, baselines)

Class representing a event-related pupillary dilation (ERPD) for one subject.

classmethod from_file(fname)

Reads a ERPD object from a pickle-file. Use as pypillometry.ERPD.from_file("yourfile.pd").

Parameters

fname (str) – filename

plot(overlays=None, meanfct=<function mean>, varfct=<function sem>, plot_missing=True)

Plot mean and error-ribbons using varct.

Parameters

  • overlays (single or sequence of ERPDSingleSubject-objects) – the overlays will be added to the same plot

  • meanfct (callable) – mean-function to apply to the single-trial ERPDs for plotting

  • varfct (callable or None) – function to calculate error-bands (e.g., numpy.std() for standard-deviation or scipy.stats.sem() for standard-error) if None, no error bands are plotted

  • plot_missing (bool) – plot percentage interpolated/missing data per time-point?

summary()

Return a summary of the PupilData-object.

Return type

dict

write_file(fname)

Save to file (using pickle).

Parameters

fname (str) – filename

pypillometry.erpd.group_erpd(datasets, erpd_name, event_select, baseline_win=None, time_win=(-500, 2000), subj_meanfct=<function mean>)

Calculate group-level ERPDs by applying subj_meanfct to each subj-level ERPD.

Parameters

  • datasets (list of PupilData objects) – one PupilData object for each subject that should go into the group-level ERPD.

  • erpd_name (str) – identifier for the result (e.g., “cue-locked” or “conflict-trials”)

  • baseline_win (tuple (float,float) or None) –

    if None, no baseline-correction is applied if tuple, the mean value in the window in milliseconds (relative to time_win) is

    subtracted from the single-trial ERPDs (baseline-correction)

  • event_select (str or function) – variable describing which events to select and align to - if str: use all events whose label contains the string - if function: apply function to all labels, use those where the function returns True

  • time_win (Tuple[float, float]) – time before and after event to include (in ms)

  • subj_meanfct (fct) – function to summarize each individual ERPD

Returns

Return type

an ERPD object for the group

pypillometry.erpd.plot_erpds(erpds)

Plot a list of ERPD objects.

io.py

Read/Write data from/to disk.

pypillometry.io.pd_read_pickle(fname)

Read the PupilData-object pdobj from file using pickle.

Parameters

fname (str) – filename or URL to load data from

Returns

pdobj – loaded dataset

Return type

PupilData

pypillometry.io.pd_write_pickle(pdobj, fname)

Store the PupilData-object pdobj in file using pickle.

Parameters

  • pdobj (PupilData) – dataset to save

  • fname (str) – filename to save to

preproc.py

preprocessing functions (blinks, smoothing, missing data…)

Method for finding the on- and offset for each blink (excluding transient). See https://figshare.com/articles/A_simple_way_to_reconstruct_pupil_size_during_eye_blinks/688001.

Parameters

  • sy (np.array) – pupil data

  • blinks (np.array (nblinks x 2)) – blink onset/offset matrix (contiguous zeros)

  • smooth_winsize (int (odd)) – size of the Hanning-window in sampling points

  • vel_onset (float) – velocity-threshold to detect the onset of the blink

  • vel_offset (float) – velocity-threshold to detect the offset of the blink

  • margin (tuple (int,int)) – margin that is subtracted/added to onset and offset (in sampling points)

  • blinkwindow (int) – how much time before and after each blink to include (in sampling points)

Detect blinks as everything between a fast downward and a fast upward-trending PD-changes.

This works similarly to blink_onsets_mahot().

Parameters

  • sy (np.array) – pupil data

  • smooth_winsize (int (odd)) – size of the Hanning-window in sampling points

  • vel_onset (float) – velocity-threshold to detect the onset of the blink

  • vel_offset (float) – velocity-threshold to detect the offset of the blink

  • min_onset_len (int) – minimum number of consecutive samples that cross the threshold to detect onset

  • min_offset_len (int) – minimum number of consecutive samples that cross the threshold to detect offset

Detect blinks as consecutive sequence of blink_val (f.eks., 0 or NaN) of at least min_duration successive values in the signal sy. Detected blinks are put a matrix blinks (nblinks x 2) where start and end are stored as indexes.

Parameters

  • sy (np.array (float)) – signal

  • min_duration (int) – minimum number of consecutive samples for a sequence of missing numbers to be treated as blink

  • blink_val – “missing value” code

Returns

Return type

np.array (nblinks x 2) containing the indices of the start/end of the blinks

pypillometry.preproc.smooth_window(x, window_len=11, window='hanning')

smooth the data using a window with requested size.

This method is based on the convolution of a scaled window with the signal. The signal is prepared by introducing reflected copies of the signal (with the window size) in both ends so that transient parts are minimized in the begining and end part of the output signal.

adapted from SciPy Cookbook: https://scipy-cookbook.readthedocs.io/items/SignalSmooth.html.

Parameters

  • x (the input signal) –

  • window_len (the dimension of the smoothing window; should be an odd integer) –

  • window (the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman') – flat window will produce a moving average smoothing.

Returns

np.array

Return type

the smoothed signal

fakedata.py

Generate artificial pupil-data.

pypillometry.fakedata.generate_pupil_data(event_onsets, fs=1000, pad=5000, baseline_lowpass=0.2, evoked_response_perc=0.02, response_fluct_sd=1, prf_npar=10.35, 0, prf_tmax=917.0, 0, prop_spurious_events=0.2, noise_amp=0.0005)

Generate artificial pupil data as a sum of slow baseline-fluctuations on which event-evoked responses are “riding”.

Parameters

  • event_onsets (list) – list of all events that evoke a response (in seconds)

  • fs (float) – sampling rate in Hz

  • pad (float) – append pad milliseconds of signal after the last event is decayed

  • baseline_lowpass (float) – cutoff for the lowpass-filter that defines the baseline (highest allowed frequency in the baseline fluctuations)

  • evoked_response_perc (float) – amplitude of the pupil-response as proportion of the baseline

  • response_fluct_sd (float) – How much do the amplitudes of the individual events fluctuate? This is determined by drawing each individual pupil-response to a single event from a (positive) normal distribution with mean as determined by evoked_response_perc and sd response_fluct_sd (in units of evoked_response_perc).

  • prf_npar (tuple (float,float)) – (mean,std) of the npar parameter from pypillometry.pupil.pupil_kernel(). If the std is exactly zero, then the mean is used for all pupil-responses. If the std is positive, npar is taken i.i.d. from ~ normal(mean,std) for each event.

  • prf_tmax (tuple (float,float)) – (mean,std) of the tmax parameter from pypillometry.pupil.pupil_kernel(). If the std is exactly zero, then the mean is used for all pupil-responses. If the std is positive, tmax is taken i.i.d. from ~ normal(mean,std) for each event.

  • prop_spurious_events (float) – Add random events to the pupil signal. prop_spurious_events is expressed as proportion of the number of real events.

  • noise_amp (float) – Amplitude of random gaussian noise that sits on top of the simulated signal. Expressed in units of mean baseline pupil diameter.

Returns

  • tx, sy (np.array) – time and simulated pupil-dilation (n)

  • x0 (np.array) – baseline (n)

  • delta_weights (np.array) – pupil-response strengths (len(event_onsets))

pypillometry.fakedata.get_dataset(ntrials=100, isi=2000, rtdist=1000, 500, fs=1000, pad=5000, **kwargs)

Convenience function to run generate_pupil_data() with standard parameters. :param ntrials: number of trials :type ntrials: int :param isi: inter-stimulus interval in milliseconds :type isi: float :param rtdist: mean and std of a (truncated at zero) normal distribution to generate response times :type rtdist: tuple (float,float) :param fs: sampling rate :type fs: float :param pad: padding before the first and after the last event in seconds :type pad: float :param kwargs: arguments for pypillometry.fakedata.generate_pupil_data() :type kwargs: dict

Returns

  • tx, sy (np.array) – time and simulated pupil-dilation (n)

  • baseline (np.array) – baseline (n)

  • event_onsets (np.array) – timing of the simulated event-onsets (stimuli and responses not separated)

  • response_coef (np.array) – pupil-response strengths (len(event_onsets))

baseline.py

Functions for estimating tonic pupillary fluctuations (baseline).

pypillometry.baseline.baseline_envelope(tx, sy, event_onsets, fs=1000, lp=2, prominence_thr=80, interp_method='cubic')

Extract baseline based on the lower envelope of the (filtered) signal.

Steps:

  • filter away noise

  • detect high-prominence throughs in the signal

  • calculate lower envelope based on these peaks

Parameters

  • tx (np.ndarray) – time-vector in seconds

  • sy (np.ndarray) – raw pupil signal

  • event_onsets (list) – onsets of events (stimuli/responses) in seconds

  • fs (float) – sampling rate in Hz

  • lp (float) – low-pass filter cutoff for removing random noise

  • prominence_thr (float in [0,100]) – percentile of the prominence distribution (of the peaks) to use for determining prominent peaks (see scipy.stats.peak_prominences())

  • interp_method (string, one of ["linear", "cubic", "spline"]) –

    “linear” - linear interpolation between the high-prominence peaks “cubic” - cubic interpolation through all high-prominence peaks “spline” - a smoothing spline that is guaranteed to go through all

    high-prominence peaks and smoothes through all the other (lower-prominence) peaks

Returns

base – baseline estimate

Return type

np.array

pypillometry.baseline.baseline_envelope_iter_bspline(tx, sy, event_onsets, fs, fsd=10, lp=2, lam_sig=1, lam_min=1, lam_max=100, verbose=0)

Extract baseline based (re-)estimating the lower envelope using B-splines. See notebook baseline_interp.ipynb for details. The signal is downsampled (to fsd Hz) for speed.

Parameters

  • tx (np.ndarray) – time-vector in seconds

  • sy (np.ndarray) – raw pupil signal

  • event_onsets (list) – onsets of events (stimuli/responses) in milliseconds

  • fs (float) – sampling rate in Hz

  • fsd (float) – downsampled sampling rate (if too slow, decrease)

  • lp (float) – lowpass-filter cutoff (Hz)

  • lam_sig (float) – parameter steering how much the baseline is shaped by the non-peaks of the signal

  • lam_min,lam_max (float) – parameters mapping how much low- and high-prominence peaks influence the baseline

  • verbose (int [0, 100]) – how much information to print (0 nothing, 100 everything)

Returns

(txd,syd,base2,base1) – txd: downsampled time-array syd: downsampled and lowpass-filtered pupil signal base1: is the estimated base after the first round base2: is the final baseline estimate

Return type

tuple

pypillometry.baseline.baseline_pupil_model(tx, sy, event_onsets, fs=1000, lp1=2, lp2=0.2)

Extract baseline based on filtering after removing stim-locked activity.

Steps:

  • filter away noise

  • regress out event-locked activity from the filtered signal using NNLS

  • remove modeled signal from filtered data

  • run another lowpass-filter to get rid of spurious signals

Parameters

  • tx (np.ndarray) – time-vector in seconds

  • sy (np.ndarray) – raw pupil signal

  • event_onsets (list) – onsets of events (stimuli/responses) in seconds

  • fs (float) – sampling rate in Hz

  • lp1 (float) – low-pass filter cutoff for removing random noise

  • lp2 (float) – low-pass filter cutoff for removing spurious peaks from the baseline-signal

Returns

base – baseline estimate

Return type

np.array

pypillometry.baseline.bspline(txd, knots, spline_degree=3)

Re-implementation from https://mc-stan.org/users/documentation/case-studies/splines_in_stan.html. Similar behaviour as R’s bs() function from the splines-library.

Parameters

  • txd (np.array) – time-vector

  • knots (np.array) – location of the knots

  • spline_degree (int) – degree of the spline

Returns

B – matrix of basis functions

Return type

np.array

pypillometry.baseline.butter_lowpass(cutoff, fs, order=5)

Get lowpass-filter coefficients for Butterworth-filter.

Parameters

  • cutoff (float) – lowpass-filter cutoff

  • fs (float) – sampling rate

  • order (int) – filter order

Returns

(b,a) – filter coefficients

Return type

tuple (float,float)

pypillometry.baseline.butter_lowpass_filter(data, cutoff, fs, order=5)

Lowpass-filter signal using a Butterworth-filter.

Parameters

  • data (np.array) – data to lowpass-filter

  • cutoff (float) – lowpass-filter cutoff

  • fs (float) – sampling rate

  • order (int) – filter order

Returns

y – filtered data

Return type

np.array

pypillometry.baseline.downsample(y, R)

Simple downsampling scheme using mean within the downsampling window.

Parameters

  • y (np.array) – signal to downsample

  • R (int) – decimate-factor

Returns

y – downsampled data

Return type

np.array

pupil.py

functions related to pupillary responses.

pypillometry.pupil.plot_prf(npar=10.1, tmax=930, max_duration='estimate', fs=500, **kwargs)

Plot profile of the pupil-response function (PRF) with parameters npar and tmax.

pypillometry.pupil.pupil_build_design_matrix(tx, event_onsets, fs, npar, tmax, max_duration='estimate')

Construct design matrix (nevents x n). Each column has a single pupil-kernel with parameters npar, tmax starting at each event_onset.

Parameters

  • tx (np.array) – in ms

  • event_onsets (np.array) – timing of the events

  • fs (float) – sampling rate (Hz)

  • npar (float) – n in the equation of pypillometry.pupil.pupil_kernel()

  • tmax (float) – t_{max} in pypillometry.pupil.pupil_kernel()

  • max_duration (float or "estimate") – either maximum duration in milliseconds or string “estimate” which causes the function to determine an appropriate duration based on the kernel parameters npar and tmax

Returns

x1 – design matrix

Return type

np.array (nevents x n)

pypillometry.pupil.pupil_get_max_duration(npar, tmax, thr=1e-08, stepsize=1.0)

Get the time when the PRF with parameters \(n\) and \(t_{max}\) is decayed to thr. This gives an indication of how long the duration parameter in pupil_kernel() should be chosen.

Parameters

  • npar,tmax (float) – PRF parameters, see pupil_kernel()

  • thr (float) – desired value to which the PRF is decayed (the lower thr, the longer the duration)

  • stepsize (float) – precision of the maximum duration in ms

Returns

tdur – first value of t so that PRF(t)<thr

Return type

float

pypillometry.pupil.pupil_kernel(duration=4000, fs=1000, npar=10.1, tmax=930.0)

According to Hoeks and Levelt (1993, https://link.springer.com/content/pdf/10.3758%2FBF03204445.pdf), the PRF can be described by the so-called Erlang gamma function $$h_{HL}(t) = t^n e^{-nt/t_{max}}$$ which we normalize to $$h(t)=\frac{1}{h_{max}} h_{HL}(t)$$ where $$h_{max} = \max_t{\left(h_{HL}(t)\right)} = e^{-n}t_{max}^{n}$$ which yields a maximum value of 1 for this function. The function \(h(t)\) is implemented in pp.pupil_kernel().

Parameters

  • duration (float) – in ms; maximum of the time window for which to calculate the PRF [0,duration]

  • fs (float) – sampling rate for resolving the PRF

  • npar (float) – n in the equation above

  • tmax (float) – t_{max} in the equation above

Returns

h – sampled version of h(t) over the interval [0,`duration`] with sampling rate fs

Return type

np.array

pypillometry.pupil.pupil_kernel_t(t, npar, tmax)

According to Hoeks and Levelt (1993, https://link.springer.com/content/pdf/10.3758%2FBF03204445.pdf), the PRF can be described by the so-called Erlang gamma function $$h_{HL}(t) = t^n e^{-nt/t_{max}}$$ which we normalize to $$h(t)=\frac{1}{h_{max}} h_{HL}(t)$$ where $$h_{max} = \max_t{\left(h_{HL}(t)\right)} = e^{-n}t_{max}^{n}$$ which yields a maximum value of 1 for this function. The function \(h(t)\) is implemented in pupil_kernel().

This version of the function evaluates the PRF at inputs t.

Parameters

  • t (float/np.array) – in ms

  • npar (float) – n in the equation above

  • tmax (float) – t_{max} in the equation above

Returns

h – PRF evaluated at t

Return type

np.array

pypillometry.pupil.pupil_response(tx, sy, event_onsets, fs, npar='free', tmax='free', verbose=10, bounds={'npar': 1, 20, 'tmax': 100, 2000}, display_progress=True)

Estimate pupil-response based on event-onsets.

txnp.ndarray

time-vector in milliseconds

synp.ndarray

(baseline-corrected) pupil signal

event_onsetslist

onsets of events (stimuli/responses) in milliseconds

fsfloat

sampling rate in Hz

npar: float

npar-parameter for the canonical response-function or “free”; in case of “free”, the function optimizes for this parameter

tmax: float

tmax-parameter for the canonical response-function or “free”; in case of “free”, the function optimizes for this parameter

bounds: dict

in case that one or both parameters are estimated, give the lower and upper bounds for the parameters

pypillometry.pupil.pupilresponse_nnls(tx, sy, event_onsets, fs, npar=10.1, tmax=930)

Estimate single-event pupil responses based on canonical PRF (pupil_kernel()) using non-negative least-squares (NNLS).

Parameters

  • tx (np.ndarray) – time-vector in milliseconds

  • sy (np.ndarray) – (baseline-corrected) pupil signal

  • event_onsets (list) – onsets of events (stimuli/responses) in seconds

  • fs (float) – sampling rate in Hz

  • npar,tmax (float) – parameters for pypillometry.pupil.pupil_kernel()

Returns

(coef,pred,resid) – coef: purely-positive regression coefficients pred: predicted signal resid: residuals (sy-pred)

Return type

tuple

pypillometry.pupil.stat_event_interval(tx, sy, event_onsets, interval, statfct=<function mean>)

Return result of applying a statistical function to pupillometric data in a given interval relative to event-onsets. For example, extract mean pupil-size in interval before trial onset.

Parameters

  • tx (np.ndarray) – time-vector in milliseconds

  • sy (np.ndarray) – (baseline-corrected) pupil signal

  • event_onsets (list) – onsets of events (stimuli/responses) in seconds

  • interval (tuple (min,max)) – time-window in ms relative to event-onset (0 is event-onset)

  • statfct (function) – function mapping np.array to a single number

Returns

result – number of event-onsets long result array

Return type

np.array

convenience.py

Some convenience functions.

pypillometry.convenience.StanModel_cache(model_code, model_name=None, **kwargs)

Use just as you would stan

pypillometry.convenience.p_asym_laplac(y, mu, sigma, tau)

Asymmetric laplace distribution https://en.wikipedia.org/wiki/Asymmetric_Laplace_distribution; Parametrization as here: https://github.com/stan-dev/stan/issues/2312

tau in [0,1]

pypillometry.convenience.p_asym_laplac_kappa(y, mu, lam, kappa)

Asymmetric laplace distribution https://en.wikipedia.org/wiki/Asymmetric_Laplace_distribution; Wikipedia parametrization.

kappa in [0, infty] where 1 means symmetry

pypillometry.convenience.plot_pupil_ipy(tx, sy, event_onsets=None, overlays=None, overlay_labels=None, blinks=None, interpolated=None, figsize=16, 8, xlab='ms', nsteps=100)

Plotting with interactive adjustment of plotting window. To use this, do

$ pip install ipywidgets $ jupyter nbextension enable –py widgetsnbextension $ jupyter labextension install @jupyter-widgets/jupyterlab-manager

Parameters

  • tx (np.ndarray) – time-vector in seconds

  • sy (np.ndarray) – raw pupil signal

  • event_onsets (list) – onsets of events (stimuli/responses) in seconds

  • overlays (tuple of np.array) – signals to overlay over the plot, given as tuple of arrays of same length as tx

  • overlay_labels (tuple of strings) – labels for the overlays to be displayed in the legend

  • figsize (tuple of int) – dimensions for the plot

  • xlab (str) – label for x-axis

  • nsteps (int) – number of steps for slider

pypillometry.convenience.sizeof_fmt(num, suffix='B')

Convert number of bytes in num into human-readable string representation. Taken from https://stackoverflow.com/a/1094933

pypillometry.convenience.trans_logistic_vec(x, a, b, inverse=False)

vectorized version of trans_logistic()

goes from [a,b] to [-inf,+inf] and back; inverse=False: [a,b] -> [-inf, +inf] inverse=True: [-inf,+inf] -> [a,b] if a or b is +/-infty, a logarithmic/exponential transform is used

Indices and tables