from .frontend_classes import Dataset, RePro, Stimulus
from .util import BoltzmannFit, unzip_if_needed, gaussian_kernel, zero_crossings
import numpy as np
import nixio as nix
from scipy.stats import circstd
# from scipy.optimize import curve_fit
import os
import subprocess
from tqdm import tqdm
from IPython import embed
class BaselineData:
"""
Class representing the Baseline data that has been recorded within a given Dataset.
"""
def __init__(self, dataset: Dataset):
self.__spike_data = []
self.__eod_data = []
self.__eod_times = []
self.__dataset = dataset
self.__repros = None
self.__cell = dataset.cells[0] # Beware: Assumption that there is only a single cell
self._get_data()
def _get_data(self):
if not self.__dataset:
return
self.__repros = RePro.find("BaselineActivity", cell_id=self.__cell.id)
for i in tqdm(range(len(self.__repros)), desc="loading data"):
r = self.__repros[i]
sd = self.__read_spike_data(r)
if sd is not None and len(sd) > 1:
self.__spike_data.append(sd)
else:
continue
self.__eod_data.append(self.__read_eod_data(r, self.__spike_data[-1][-1]))
def valid(self):
# fixme implement me!
pass
def __read_spike_data(self, r: RePro):
if self.__dataset.has_nix:
return self.__read_spike_data_from_nix(r)
else:
return self.__read_spike_data_from_directory(r)
def __read_eod_data(self, r: RePro, duration):
if self.__dataset.has_nix:
return self.__read_eod_data_from_nix(r, duration)
else:
return self.__read_eod_data_from_directory(r, duration)
def __get_serial_correlation(self, times, max_lags=50):
if times is None or len(times) < max_lags:
return None
isis = np.diff(times)
unbiased = isis - np.mean(isis, 0)
norm = sum(unbiased ** 2)
a_corr = np.correlate(unbiased, unbiased, "same") / norm
a_corr = a_corr[int(len(a_corr) / 2):]
return a_corr[:max_lags]
def serial_correlation(self, max_lags=50):
"""
Returns the serial correlation of the interspike intervals.
@param max_lags: The number of lags to take into account
@return: the serial correlation as a function of the lag
"""
scs = []
for sd in self.__spike_data:
if sd is None or len(sd) < 100:
continue
corr = self.__get_serial_correlation(sd, max_lags=max_lags)
if corr is not None:
scs.append(corr)
return scs
def circular_std(self):
circular_stds = []
for i in range(self.size):
phases = self.__spike_phases(index=i)
circular_stds.append(circstd(phases))
return circular_stds
def eod_frequency(self):
eod_frequencies = []
for i in range(self.size):
eod, time = self.eod(i)
xings = zero_crossings(eod, time, interpolate=False)
eod_frequencies.append(len(xings)/xings[-1])
return 0.0 if len(eod_frequencies) < 1 else np.mean(eod_frequencies)
def __spike_phases(self, index=0): # fixme buffer this stuff
e_times = self.eod_times(index=index)
eod_period = np.mean(np.diff(e_times))
phases = np.zeros(len(self.spikes(index)))
for i, st in enumerate(self.spikes(index)):
last_eod_index = np.where(e_times <= st)[0]
if len(last_eod_index) == 0:
continue
phases[i] = (st - e_times[last_eod_index[-1]]) / eod_period * 2 * np.pi
return phases
def eod_times(self, index=0, interpolate=True):
if index >= self.size:
return None
if len(self.__eod_times) < len(self.__eod_data):
eod, time = self.eod(index)
times = zero_crossings(eod, time, interpolate=interpolate)
else:
times = self.__eod_times[index]
return times
@property
def dataset(self):
return self.__dataset
@property
def cell(self):
cells = self.__dataset.cells
return cells if len(cells) > 1 else cells[0]
@property
def subject(self):
subjects = self.__dataset.subjects
return subjects if len(subjects) > 1 else subjects[0]
def spikes(self, index: int=0):
"""Get the spike times of the spikes recorded in the given baseline recording.
Args:
index (int, optional): If the baseline activity has been recorded several times, the index can be given. Defaults to 0.
Returns:
: [description]
"""
return self.__spike_data[index] if len(self.__spike_data) >= index else None
def membrane_voltage(self, index: int=0):
if index >= self.size:
raise IndexError("Index %i out of bounds for size %i!" % (index, self.size))
if not self.__dataset.has_nix:
print('Sorry, this is not supported for non-nixed datasets. Implement it at '
'fishbook.reproclasses.BaselineData.membrane_voltage and send a pull request!')
return None, None
else:
rp = self.__repros[index]
data_source = os.path.join(self.__dataset.data_source, self.__dataset.id + ".nix")
f = nix.File.open(data_source, nix.FileMode.ReadOnly)
b = f.blocks[0]
t = b.tags[rp.id]
if not t:
print("Tag not found!")
try:
data = t.retrieve_data("V-1")[:]
time = np.asarray(t.references["V-1"].dimensions[0].axis(len(data)))
except:
data = np.empty()
time = np.empty()
f.close()
return time, data
def eod(self, index: int=0):
eod = self.__eod_data[index] if len(self.__eod_data) >= index else None
time = np.arange(len(eod)) / self.__dataset.samplerate
return eod, time
@property
def burst_index(self):
bi = []
for i, sd in enumerate(self.__spike_data):
if len(sd) < 2:
continue
et = self.eod_times(index=i)
eod_period = np.mean(np.diff(et))
isis = np.diff(sd)
bi.append(np.sum(isis < (1.5 * eod_period))/len(isis))
return bi
@property
def coefficient_of_variation(self):
cvs = []
for d in self.__spike_data:
isis = np.diff(d)
cvs.append(np.std(isis)/np.mean(isis))
return cvs
@property
def vector_strength(self):
vss = []
spike_phases = []
for i, sd in enumerate(self.__spike_data):
phases = self.__spike_phases(i)
ms_sin_alpha = np.mean(np.sin(phases)) ** 2
ms_cos_alpha = np.mean(np.cos(phases)) ** 2
vs = np.sqrt(ms_cos_alpha + ms_sin_alpha)
vss.append(vs)
spike_phases.append(phases)
return vss, spike_phases
@property
def size(self):
return len(self.__spike_data)
def __str__(self):
return "Baseline data of cell %s " % self.__cell.id
def __read_eod_data_from_nix(self, r: RePro, duration) -> np.ndarray:
data_source = os.path.join(self.__dataset.data_source, self.__dataset.id + ".nix")
if not os.path.exists(data_source):
print("Data not found! Trying from directory")
return self.__read_eod_data_from_directory(r, duration)
f = nix.File.open(data_source, nix.FileMode.ReadOnly)
b = f.blocks[0]
t = b.tags[r.id]
if not t:
print("Tag not found!")
try:
data = t.retrieve_data("EOD")[:]
except:
data = np.empty()
f.close()
return data
def __read_eod_data_from_directory(self, r: RePro, duration) -> np.ndarray:
sr = self.__dataset.samplerate
unzip_if_needed(self.__dataset.data_source, "trace-2.raw")
eod = np.fromfile(self.__dataset.data_source + "/trace-2.raw", np.float32)
eod = eod[:int(duration * sr)]
return eod
def __read_spike_data_from_nix(self, r: RePro) -> np.ndarray:
data_source = os.path.join(self.__dataset.data_source, self.__dataset.id + ".nix")
if not os.path.exists(data_source):
print("Data not found! Trying from directory")
return self.__read_spike_data_from_directory(r)
f = nix.File.open(data_source, nix.FileMode.ReadOnly)
b = f.blocks[0]
t = b.tags[r.id]
if not t:
print("Tag not found!")
try:
data = t.retrieve_data("Spikes-1")[:]
except:
data = None
f.close()
if len(data) < 100:
data = None
return data
def __read_spike_data_from_directory(self, r) -> np.ndarray:
data = []
data_source = os.path.join(self.__dataset.data_source, "basespikes1.dat")
if os.path.exists(data_source):
found_run = False
with open(data_source, 'r') as f:
l = f.readline()
while l:
if "index" in l:
index = int(l.strip("#").strip().split(":")[-1])
found_run = index == r.run
if l.startswith("#Key") and found_run:
data = self.__do_read(f)
break
l = f.readline()
if len(data) < 100:
return None
return np.asarray(data)
@staticmethod
def __do_read(f) -> np.ndarray:
data = []
f.readline()
unit = f.readline().strip("#").strip()
scale = 0.001 if unit == "ms" else 1
l = f.readline()
while l and "#" not in l and len(l.strip()) > 0:
data.append(float(l.strip())*scale)
l = f.readline()
return np.asarray(data)
class FIData:
"""
Class representing the data recorded with the relacs FI-Curve repro. The instance will load the data upon
construction which may take a while.
FI Data offers convenient access to the spike and local EOD data as well as offers conveince methods to get the
firing rate and also to fit a Boltzmann function to the the FI curve.
"""
def __init__(self, dataset: Dataset):
"""
Constructor.
:param dataset: The dataset entity for which the fi curve repro data should be loaded.
"""
self.__spike_data = []
self.__contrasts = []
self.__eod_data = []
self.__eod_times = []
self.__dataset = dataset
self.__repros = None
self.__cell = dataset.cells[0] # Beware: Assumption that there is only a single cell
self._get_data()
pass
def _get_data(self):
if not self.__dataset:
return
self.__repros = RePro.find("FICurve", cell_id=self.__cell.id)
for r in self.__repros:
sd, c, eods, time = self.__read_spike_data(r)
if sd is not None and len(sd) > 1:
self.__spike_data.extend(sd)
if eods:
self.__eod_data.extend(eods)
self.__contrasts.extend(c)
if time:
self.__eod_times.extend(time)
else:
continue
def __read_spike_data(self, repro: RePro):
"""
:param repro:
:return: spike data and the respective contrasts
"""
if self.__dataset.has_nix:
return self.__read_spikes_from_nix(repro)
else:
return self.__read_spikes_from_directory(repro)
def __do_read_spike_data_from_nix(self, mtag: nix.pycore.MultiTag, stimulus: Stimulus, repro: RePro):
r_settings = repro.settings.split("\n")
s_settings = stimulus.settings.split("\n")
delay = 0.0
contrast = 0.0
for s in r_settings:
if "delay:" in s:
delay = float(s.split(":")[-1])
break
for s in s_settings:
if "Contrast:" in s and "PreContrast" not in s and "\t\t" not in s and "+-" not in s:
contrast = float(s.split(":")[-1])
break
start_time = stimulus.start_time - delay
end_time = stimulus.start_time + stimulus.duration
eod_da = mtag.references["LocalEOD-1"]
start_index_eod = eod_da.dimensions[0].index_of(start_time)
end_index_eod = eod_da.dimensions[0].index_of(end_time)
local_eod = eod_da[start_index_eod:end_index_eod]
spikes = self.__all_spikes[(self.__all_spikes >= start_time) & (self.__all_spikes < end_time)] - start_time - delay
time = np.asarray(eod_da.dimensions[0].axis(end_index_eod - start_index_eod)) - delay
return spikes, local_eod, time, contrast
def __read_spikes_from_nix(self, repro: RePro):
spikes = []
eods = []
time = []
contrasts = []
stimuli = Stimulus.find(cell_id=repro.cell_id, repro_id=repro.id)
if len(stimuli) == 0:
return spikes, contrasts, eods, time
data_source = os.path.join(self.__dataset.data_source, self.__dataset.id + ".nix")
if not os.path.exists(data_source):
print("Data not found! Trying from directory")
return self.__read_spikes_from_directory(repro)
f = nix.File.open(data_source, nix.FileMode.ReadOnly)
b = f.blocks[0]
self.__all_spikes = b.data_arrays["Spikes-1"][:]
mt = None
for i in tqdm(range(len(stimuli)), desc="Loading data"):
s = stimuli[i]
if not mt or mt.id != s.multi_tag_id:
mt = b.multi_tags[s.multi_tag_id]
sp, eod, t, c = self.__do_read_spike_data_from_nix(mt, s, repro)
spikes.append(sp)
eods.append(eod)
time.append(t)
contrasts.append(c)
f.close()
return spikes, contrasts, eods, time
def __do_read_data_block(self, f, l):
spikes = []
while len(l.strip()) > 0 and "#" not in l:
spikes.append(float(l.strip())/1000)
l = f.readline()
return spikes, l
def __read_spikes_from_directory(self, repro: RePro):
spikes = []
contrasts = []
times = []
print("Warning! Exact reconstruction of stimulus order not possible for old relacs files!")
data_source = os.path.join(self.__dataset.data_source, "fispikes1.dat")
delay = 0.0
pause = 0.0
for s in repro.settings.split(", "):
if "pause" in s:
s = s.split(":")[-1].strip()
pause = float(s[:-2])/1000 if "ms" in s else float(s[:-1])
if "delay" in s:
s = s.split(":")[-1].strip()
delay = float(s[:-2])/1000 if "ms" in s else float(s[:-1])
t_start = -delay
t_end = pause + delay
time = np.arange(t_start, t_end, 1./repro.dataset.samplerate)
if os.path.exists(data_source):
with open(data_source, 'r') as f:
line = f.readline()
fish_intensity = None
stim_intensity = None
while line:
line = line.strip().lower()
if "index" in line:
fish_intensity = None
stim_intensity = None
if "intensity = " in line:
if "true intensity = " in line:
fish_intensity = float(line.split("=")[-1].strip()[:-2])
elif "pre" not in line:
stim_intensity = float(line.split("=")[-1].strip()[:-2])
if len(line) > 0 and "#" not in line: # data line
sp, line = self.__do_read_data_block(f, line)
spikes.append(sp)
times.append(time)
contrasts.append((stim_intensity/fish_intensity-1)*100)
continue
line = f.readline()
return spikes, contrasts, None, times
@property
def size(self) -> int:
"""
The number of recorded trials
:return: An integer with the number of trials.
"""
return len(self.__spike_data)
def spikes(self, index=-1):
"""
The spike times recorded in the specified trial(s)
:param index: the index of the trial. Default of -1 indicates that all data should be returned.
:return:
"""
if 0 <= index < self.size:
return self.__spike_data[index]
else:
return self.__spike_data
def eod(self, index=-1):
"""
The local eod (including the stimulus) measurement of the selected trial(s).
:param index: the index of the trial. Default of -1 indicates that all data should be returned.
:return: Either two vectors representing time and the local eod or two lists of such vectors
"""
if len(self.__eod_data) == 0:
print("EOD data not available for old-style relacs data.")
return None, None
if 0 <= index < self.size:
return self.__eod_times[index], self.__eod_data[index]
else:
return self.__eod_times, self.__eod_data
def contrast(self, index=-1):
"""
The stimulus contrast used in the respective trial(s).
:param index: the index of the trial. Default of -1 indicates that all data should be returned.
:return: Either a single scalar representing the contrast, or a list of such scalars, one entry for each trial.
"""
if 0 <= index < self.size:
return self.__contrasts[index]
else:
return self.__contrasts
def time_axis(self, index=-1):
"""
Get the time axis of a single trial or a list of time-vectors for all trials.
:param index: the index of the trial. Default of -1 indicates that all data should be returned.
:return: Either a single vector representing time, or a list of such vectors, one for each trial.
"""
if 0 <= index < self.size:
return self.__eod_times[index]
else:
return self.__eod_times
def rate(self, index=0, kernel_width=0.005):
"""
Returns the firing rate for a single trial.
:param index: The index of the trial. 0 <= index < size
:param kernel_width: The width of the gaussian kernel in seconds
:return: tuple of time and rate
"""
t = self.time_axis(index)
dt = np.mean(np.diff(t))
sp = self.spikes(index)
binary = np.zeros(t.shape)
spike_indices = ((sp - t[0]) / dt).astype(int)
binary[spike_indices[(spike_indices >= 0) & (spike_indices < len(binary))]] = 1
g = gaussian_kernel(kernel_width, dt)
rate = np.convolve(binary, g, mode='same')
return t, rate
def boltzmann_fit(self, start_time=0.01, end_time=0.05, kernel_width=0.005):
"""
Extracts the average firing rate within a time window from the averaged across trial firing rate.
The analysis time window is specified by the start_time and end_time parameters. Firing rate is estimated by
convolution with a Gaussian kernel of a given width. All parameters are given in 's'.
:param start_time: the start of the analysis window.
:param end_time: the end of the analysis window.
:param kernel_width: standard deviation of the Gaussian kernel used for firing rate estimation.
:return: object of type BoltzmannFit
"""
contrasts = np.zeros(self.size)
rates = np.zeros(self.size)
for i in range(self.size):
contrasts[i] = np.round(self.contrast(i))
t, r = self.rate(i, kernel_width)
rates[i] = np.mean(r[(t >= start_time) & (t < end_time)])
boltzmann_fit = BoltzmannFit(contrasts, rates)
return boltzmann_fit
class FileStimulusData:
def __init__(self, dataset: Dataset):
"""
Constructor.
:param dataset: The dataset entity for which the filestimulus repro data should be loaded.
"""
self.__spike_data = []
self.__contrasts = []
self.__stimuli = []
self.__dataset = dataset
self.__repros = None
self.__cell = dataset.cells[0] # Beware: Assumption that there is only a single cell
self._get_data()
def _get_data(self):
if not self.__dataset:
return
self.__repros = RePro.find("FileStimulus", cell_id=self.__cell.id)
for r in self.__repros:
sd, c, stims = self.__read_spike_data_from_nix(r) if self.__dataset.has_nix else self.__read_spike_data_from_directory(r)
if sd is not None and len(sd) > 1:
self.__spike_data.extend(sd)
self.__contrasts.extend(c)
self.__stimuli.extend(stims)
else:
continue
def __do_read_spike_data_from_nix(self, mt: nix.pycore.MultiTag, stimulus: Stimulus, repro: RePro):
spikes = None
contrast = 0.0
stim_file = ""
r_settings = repro.settings.split("\n")
s_settings = stimulus.settings.split("\n")
delay = 0.0
for s in r_settings:
if "delay:" in s:
delay = float(s.split(":")[-1])
break
start_time = stimulus.start_time - delay
end_time = stimulus.start_time + mt.extents[stimulus.index]
contrast = 0.0 # this is a quick fix!!!
embed()
for s in s_settings:
if "Contrast:" in s and "PreContrast" not in s and "\t\t" not in s and "+-" not in s:
contrast = float(s.split(":")[-1])
break
return spikes, contrast, stim_file
local_eod = eod_da[start_index_eod:end_index_eod]
spikes = self.__all_spikes[(self.__all_spikes >= start_time) & (self.__all_spikes < end_time)] - start_time - delay
time = np.asarray(eod_da.dimensions[0].axis(end_index_eod - start_index_eod)) - delay
return spikes, local_eod, time, contrast
return spikes, contrast, stim_file
def __read_spike_data_from_nix(self, repro: RePro):
spikes = []
contrasts = []
stim_files = []
stimuli = Stimulus.find(cell_id=repro.cell_id, repro_id=repro.id)
if len(stimuli) == 0:
return spikes, contrasts, stim_files
data_source = os.path.join(self.__dataset.data_source, self.__dataset.id + ".nix")
if not os.path.exists(data_source):
print("Data not found! Trying from directory")
return self.__read_spike_data_from_directory(repro)
f = nix.File.open(data_source, nix.FileMode.ReadOnly)
b = f.blocks[0]
self.__all_spikes = b.data_arrays["Spikes-1"][:]
mt = None
for i in tqdm(range(len(stimuli)), desc="Loading data"):
s = stimuli[i]
if not mt or mt.id != s.multi_tag_id:
mt = b.multi_tags[s.multi_tag_id]
sp, c, stim = self.__do_read_spike_data_from_nix(mt, s, repro)
spikes.append(sp)
contrasts.append(c)
stim_files.append(stim)
f.close()
return spikes, contrasts, stim_files
def __read_spike_data_from_directory(self, repro: RePro):
print("not yet my friend!")
spikes = []
contrast = 0.0
stim = None
return spikes, contrast, stim
def read_stimulus(self, index=0):
pass
if __name__ == "__main__":
# dataset = Dataset(dataset_id='2011-06-14-ag')
dataset = Dataset(dataset_id="2018-09-13-ac-invivo-1")
# dataset = Dataset(dataset_id='2013-04-18-ac')
fi_curve = FileStimulusData(dataset)
embed()