"""
Miguel Capllonch Juan
This module performs analysis of the results
"""
import numpy as np
import os
import csv
from neuron import h
import workspace as ws
def remove_previous_results():
""" Remove all previous results, either text or figures """
resultsdir = ws.folders["data/results"]
# Remove all .csv in the results
for item in os.listdir(resultsdir):
if item.endswith(".csv"):
os.remove(os.path.join(resultsdir, item))
# Remove all images
imgdir = os.path.join(resultsdir, 'images')
for item in os.listdir(imgdir):
os.remove(os.path.join(imgdir, item))
def set_time_recording():
""" Set the recording for the time vector of the simulation """
# Record time
ws.tvec = h.Vector()
ws.tvec.record(h._ref_t)
def set_recording(seg, var):
""" Set up a recording vector in hoc for any given segment
in the system """
recording = h.Vector()
# recording.buffer_size(ws.nt)
exec('recording.record(seg._ref_%s)'%var)
ws.recordings.append(recording)
ws.counters['recordings'] += 1
def record_cables():
""" Set up the recording vectors for the desired cables and the
desired variables """
dr = ws.settings["data records"]
rc = dr["record cables"]
ra = dr["record axons"]
if rc == "all":
# Record all types of cables
for cable in ws.cables:
if cable.cabletype == "NAELC":
cable.set_recordings()
elif cable.cabletype == "Axon":
if ra == "all":
# Record all axons
cable.set_recordings()
elif rc == "NAELC":
# Record only NAELC
for cable in ws.cables:
if cable.cabletype == "NAELC":
cable.set_recordings()
elif rc == "axons":
# Record only axons
if ra == "all":
# Record all axons
for cable in ws.cables:
if cable.cabletype == "Axon":
cable.set_recordings()
elif rc == "None":
# Nothing to record
pass
def remove_artefacts(time, recs):
""" Remove the stimulation artefacts in the electrode recordings """
# Start the analysis
# Turn time and recordings into arrays to be able to work with them
time = np.array(time)
recs = np.array(recs)
# Make everything positive
rpos = np.abs(recs)
# Average
ravg = rpos.mean()
# Identify where rpos is higher than its mean: there should be the artefacts
rhigh = np.where(rpos > ravg)
# Identify the separation between pulses
dt = np.gradient(time).mean()
# The 1.0000001 factor is just to make sure we get > and not >=
try:
cuts = np.where(np.gradient(time[rhigh]) > 1.000001 * dt)[0][0::2]
except ValueError:
# This may happen, for instance, if np.gradient can't be
# performed. In case of error, just say there's no cuts
cuts = []
# Split the pulses
segments = []
for cut in cuts:
ttrr = []
ttrr.append(time[rhigh][:cut + 1])
ttrr.append(recs[rhigh][:cut + 1])
segments.append(ttrr)
# Last pulse
ttrr = []
try:
ttrr.append(time[rhigh][cut + 1:])
ttrr.append(recs[rhigh][cut + 1:])
except UnboundLocalError:
# There are no cuts
pass
else:
segments.append(ttrr)
# Identify the jumps at the cuts
# Choose which data are outisde the pulses. Use rhigh as a mask
toutsp = np.ma.masked_where(rpos > ravg, time)
routsp = np.ma.masked_where(rpos > ravg, recs)
# Get the jumps from that
jumps = np.where(routsp.mask[1:] != routsp.mask[:-1])[0]
# Pulse counter
k = 0
# (Average) jump heights of each pulse
pulseheights = []
# New values for the recordings during the pulses
rnews = []
# Iterate over jumps to get the avg. jump height for each pulse
for i, jump in enumerate(jumps):
jheight = abs(recs[jump + 1] - recs[jump])
if i % 2 == 0:
jumpup = jheight
elif len(segments) > 0:
rawrec = segments[k][1]
avgjh = 0.5 * (jumpup + jheight) * np.sign(rawrec)
# Subtract this value from the corresponding segment in recs
rnew = rawrec - avgjh
rnews.append(rnew)
# Add one to the pulse counter
k += 1
else:
rnews = recs.copy()
# Finally, get the entire thing
finalrecs = recs.copy()
if len(segments) > 0:
finalrecs[rhigh] = [item for sublist in rnews for item in sublist]
return finalrecs
def save_vmem(data):
""" Save the membrane potentials """
f_sti = open(os.path.join(ws.folders["data/results"], 'vmemb_stm_site.csv'), 'w')
f_rec = open(os.path.join(ws.folders["data/results"], 'vmemb_rec_site.csv'), 'w')
writers = {
'Sti': csv.writer(f_sti),
'Rec': csv.writer(f_rec)
}
for i, vecs in enumerate(data):
for k, vv in vecs.items():
writers[k].writerow([i] + vv.as_numpy().tolist())
f_sti.close()
f_rec.close()
def save_electrode_recordings():
""" Save the recordings from all electrodes into files """
# Time
tvarr = ws.tvec.as_numpy()
# Iterate over electrodes:
for name, electrode in ws.electrodes.items():
if electrode.role == 'recording':
# Iterate over rings on the electrode
for ring in electrode.rings:
ring_number = ring.ring_number
# Iterate over pads on the ring
for pad, rec in ring.recordings.items():
# Recording ID
rec_id = rec['recording id']
# Fetch the vector with the data from ws
v = ws.recordings[rec_id]
# Save with artefacts
with open(os.path.join(ws.folders["data/results"],
'recordings_E_%s_R%iP%i_withartefacts.csv'\
%(name.replace(' ', '_'), ring_number, pad)), 'w') as f:
fw = csv.writer(f)
fw.writerow(['Time (ms)', 'V (mV)'])
for tt, vv in zip(tvarr, v):
fw.writerow([tt, vv])
# Remove the artefacts
try:
v = remove_artefacts(tvarr, v.as_numpy())
except Exception as e:
message = "ERROR in %s.save_electrode_recordings:\n"%__name__ + type(e).__name__ + ':\n' + str(e)
ws.log(message)
# Don't terminate, just go on
# ws.terminate()
# Save
with open(os.path.join(ws.folders["data/results"],
'recordings_E_%s_R%iP%i_noartefacts.csv'\
%(name.replace(' ', '_'), ring_number, pad)), 'w') as f:
fw = csv.writer(f)
fw.writerow(['Time (ms)', 'V (mV)'])
for tt, vv in zip(tvarr, v):
fw.writerow([tt, vv])
def save_cable_recordings():
""" Save the recordings from cables into files """
# Time
tvarr = ws.tvec.as_numpy()
# Iterate over cables
for cable in ws.cables:
if cable.has_recordings:
# Open file
with open(os.path.join(ws.folders["data/results"], '%s%04i.csv'\
%(cable.cabletype, cable.specific_number)), 'w') as f:
fw = csv.writer(f)
# Write the most basic properties of this cable:
# Position and radius
fw.writerow([cable.x, cable.y, cable.r])
# Write the anatomy of the axon along the z-axis
lp = ["segment positions:"]
for i, seg in enumerate(ws.segments[cable.cable_number]):
lp.append(str(seg))
lp.append(cable.properties['zprofile'][i])
fw.writerow(lp)
# Iterate over recorded variables
for var, recs in cable.recordings.items():
# Write the name of the variable
fw.writerow([var])
# Iterate over recording vectors (or segments)
for rec in recs:
# Get numerical data
data = ws.recordings[rec["number"]].as_numpy()
# Store results
# fw.writerow([ws.sectypes[cable.cable_number][i]] + [x for x in data])
fw.writerow([rec["segment"]] + [x for x in data])
def save_memcurrents(i_mem_, i_mem_ly2):
""" Save the membrane currents into files """
x, y = ws.nvt.x, ws.nvt.y
with open(os.path.join(ws.folders["data/results"], 'membranecurrents.csv'), 'w') as f:
fw = csv.writer(f)
fw.writerow(['Fiber', 'secname', 'x', 'y', 'z', 'i_membrane time series', 'i_membrane on layer 2'])
for i, fiber in i_mem_.items():
for j, segment in enumerate(fiber):
fw.writerow([i, ws.segments[i][j].sec.name(), x[i],
y[i], ws.zprofs[i][j]] + segment.tolist() + \
i_mem_ly2[i][j].tolist())
def store_i_mem(i_t):
""" Store the membrane currents of all the cells.
So far, this only works for axons with 2 extracellular layers """
for i in np.arange(ws.nNAELC, ws.nc, 1):
for j, seg in enumerate(ws.segments[i]):
if 'node' in seg.sec.name():
ws.i_mem_[i][j, i_t] = seg.i_membrane * seg.area()
# Outwards is positive
ws.i_mem_ly2[i][j, i_t] = (seg.vext[0] - seg.vext[1]) * seg.xg[0] * seg.area()
else:
ws.i_mem_[i][j, i_t] = 0.
ws.i_mem_[i][j, i_t] = seg.i_membrane * seg.area()
# Outwards is positive
ws.i_mem_ly2[i][j, i_t] = (seg.vext[0] - seg.vext[1]) * seg.xg[0] * seg.area()
def postprocessing():
""" Do all the post-processing tasks after a simulation """
# Save the recordings from cables into files
save_cable_recordings()
# Save the electrode recordings into files
# if that's what we want
if ws.settings["electrode recordings"]:
# and if the nerve model we're using is compatible with that
if ws.settings["nerve model"] != "simple":
save_electrode_recordings()
