#!/usr/bin/env python3
'''
Copyright 2021 The Regents of the University of Colorado
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
Author: Christian Rickert <christian.rickert@cuanschutz.edu>
Title: ParamAP (parametrization of sinoatrial myocyte action potentials)
DOI: https://doi.org/10.5281/zenodo.823742
URL: https://github.com/christianrickert/ParamAP/
'''
# imports
# runtime
import fnmatch
import functools
import gc
import math
import os
import sys
# numpy
import numpy as np
# scipy
import scipy.signal as sp_sig
import scipy.spatial as sp_spat
import scipy.stats as sp_stat
# matplotlib
import matplotlib.backends.backend_pdf as mpbp
import matplotlib.pyplot as mpp
# variables
APEXT = 0.5 # margin of ap extension (%)
FFRAME = 0.5 # time interval for filtering (ms)
POLYNOM = 2 # polynomial order used for filtering
# functions
def askboolean(dlabel="custom boolean", dval=True):
"""Returns a boolean provided by the user."""
if dval: # True
dstr = "Y/n"
else: # False
dstr = "y/N"
while True:
uchoice = input(dlabel + " [" + dstr + "]: ") or dstr
if uchoice in ["Y/n", "Y", "y", "Yes", "yes"]:
print("True\n")
return True # break
elif uchoice in ["y/N", "N", "n", "No", "no"]:
print("False\n")
return False # break
else:
continue
def askext(dlabel="custom extension", dext='atf'):
"""Returns a file extention provided by the user."""
while True:
uext = str(input("Enter " + dlabel + " [" + dext + "]: ")).lower() or dext
if uext not in ["dat", "log", "pdf"] and len(uext) == 3:
print(uext + "\n")
return uext # break
else:
print("Invalid file extension!\n")
continue
def askunit(dlabel="custom unit", daxis='', dunit=''):
"""Returns a unit provided by the user."""
while True:
uunit = input("Enter " + dlabel + " [" + dunit + "]: ") or dunit
if daxis in ["x", "X"]:
if uunit in ["ms", "s"]:
print(uunit + "\n")
return uunit # break
else:
print("Invalid unit for X-axis!\n")
continue
elif daxis in ["y", "Y"]:
if uunit in ["mV", "V"]:
print(uunit + "\n")
return uunit # break
else:
print("Invalid unit for Y-axis!\n")
continue
def askvalue(dlabel="custom value", dval=1.0, dunit="", dtype="float"):
"""Returns a value provided by the user."""
while True:
try:
uval = float(input("Enter " + dlabel + " [" + str(dval) + "]" + dunit + ": ") or dval)
break
except ValueError:
print("Non-numerical input!\n")
continue
if dtype == "float": # default
pass
elif dtype == "int":
uval = int(round(uval))
print(str(uval) + "\n")
return uval
def getfiles(path='/home/user/', pattern='*'):
"""Returns all files in path matching the pattern."""
abspath = os.path.abspath(path)
for fileobject in os.listdir(abspath):
filename = os.path.join(abspath, fileobject)
if os.path.isfile(filename) and fnmatch.fnmatchcase(fileobject, pattern):
yield os.path.join(abspath, filename)
def getneighbors(origin_i=np.empty(0), vicinity=np.empty(0), origin_x=np.empty(0), origin_y=np.empty(0), hwidth=float("inf"), fheight=0.0, limit=None, within=float("inf"), bads=False):
"""Returns all nearest-neighbors in ascending (i.e. increasing distance) order."""
neighbors = np.zeros(0)
badorigins = np.zeros(0)
vicinity_kdt = sp_spat.KDTree(list(zip(vicinity, np.zeros(vicinity.size)))) # KDTree for the nearest neighbors search
for origin in origin_i:
neighbor_left, neighbor_right = False, False
for position in vicinity_kdt.query([origin, 0.0], k=limit, distance_upper_bound=within)[1]: # return nearest neighbors in ascending order
if not neighbor_left or not neighbor_right:
neighbor = vicinity[position]
if (abs(origin_x[origin]-origin_x[neighbor]) <= hwidth) and (abs(origin_y[origin]-origin_y[neighbor]) >= fheight): # relative criteria for minima left and right of maximum
if not neighbor_left and (neighbor < origin): # criteria for minimum left of maximum only
neighbors = np.append(neighbors, neighbor)
neighbor_left = True
elif not neighbor_right and (neighbor > origin): # criteria for minimum right of maximum only
neighbors = np.append(neighbors, neighbor)
neighbor_right = True
else: # odd origins with missing neighbors
badorigins = np.append(badorigins, np.argwhere(origin == origin_i))
neighbors = np.sort(np.unique(neighbors)) # unique elements only
if neighbors.size <= 1: # missing neighbor
if neighbor_left:
neighbors = np.append(neighbors, 0.0) # append neighbor_right
if neighbor_right:
neighbors = np.insert(neighbors, 0, 0.0) # insert neighbor_left
badorigins = np.sort(np.unique(badorigins))
return (neighbors.astype(int), badorigins.astype(int)) if bads else (neighbors.astype(int))
def getrunavg(xdata=np.empty(0), xinterval=FFRAME, porder=POLYNOM):
"""Returns the running average count based on a given time interval."""
tmprun = int(round(xinterval/(xdata[1]-xdata[0])))
while tmprun <= porder: # prevents filtering
tmprun += 1
return (tmprun) if tmprun % 2 else (tmprun + 1) # odd number
def getpartitions(pstart=0, pstop=100, pint=5, pmin=10):
"""Returns a partition list in percent to segment an interval."""
plist = []
for part_l in list(range(int(pstart), int(pstop)+int(pint), int(pint))):
for part_r in list(range(int(pstart), int(pstop)+int(pint), int(pint))):
if part_r > part_l and part_r-part_l >= int(pmin): # no duplication or empty partitions, minimum size
plist.append([part_l, part_r]) # return statement removes the outmost list
return plist
def getbestlinearfit(xaxis=np.empty(0), yaxis=np.empty(0), xmin=0.0, xmax=1.0, pstart=0, pstop=100, pint=1, pmin=10):
"""Returns the best linear fit from segments of an interval."""
bst_r = 0 # regression coefficient
seg_i = np.argwhere((xaxis >= xmin) & (xaxis <= xmax)).ravel() # analyzing partial segment only
seg_t = xaxis[seg_i[-1]]-xaxis[seg_i[0]] # full interval from partial segment
seg_m, seg_n, seg_r = 0.0, 0.0, 0.0
for partition in getpartitions(pstart, pstop, pint, pmin):
seg_i = np.argwhere((xaxis >= (avgfmin_x[0]+(seg_t*partition[0]/100))) & (xaxis <= (avgfmin_x[0]+(seg_t*partition[1]/100)))).ravel() # 'ravel()' required for 'sp_stat.linregress()'
seg_x = xaxis[seg_i]
seg_y = yaxis[seg_i]
seg_m, seg_n, seg_r = sp_stat.linregress(seg_x, seg_y)[0:3] # tuple unpacking and linear regression of partial ap segment
if math.pow(seg_r, 2.0) >= math.pow(bst_r, 2.0):
bst_m, bst_n, bst_r = seg_m, seg_n, seg_r
bst_i, bst_x, bst_y = seg_i, seg_x, seg_y
# print(partition[0], " - ", partition[1], " : ", str(partition[1]-partition[0]), " ~ ", str(math.pow(bst_r, 2.0))) # shows progress, but is slow!
return (bst_i, bst_x, bst_y, bst_m, bst_n, bst_r)
def mpp_setup(title='Plot title', xlabel='Time [ ]', ylabel='Voltage [ ]'):
"""Provides a title and axes labels to a Matplotlib plot."""
mpp.title(title)
mpp.xlabel(xlabel)
mpp.ylabel(ylabel)
def readfile(inputfile='name'):
""" read an atf file into a numpy array """
defux = ["ms", "s"]
defuy = ["mV", "V"]
inpunits = False
try:
with open(inputfile, 'r') as in_data:
full_record = [] # python list
for _ in range(0, 2):
full_record.append(in_data.readline().strip().split()) # first and second record
full_record.append([in_data.readline().strip() for line in range(0, int(full_record[1][0]))]) # optional record
full_record.append([in_data.readline().strip().split('\t')]) # title record
full_record.append(np.genfromtxt(in_data, dtype='float', delimiter='\t', unpack=True)) # data record, use: 'np.loadtxt()' for limited memory usage, but lower speed
except FileNotFoundError:
print("Error! File not found.")
read_result = None
raise
else:
inpuxy = full_record[3][0]
inpux = str(inpuxy[0]).split()[-1][1:-2]
inpuy = str(inpuxy[1]).split()[-1][1:-2]
if inpux in defux and inpuy in defuy:
inpunits = True
inp_xy = full_record[4]
read_result = inp_xy, inpunits, inpux, inpuy
return read_result
def toggle_mpp_imode(activate=True):
"""Toggle matplotlib interactive mode on or off."""
inter_active = mpp.isinteractive()
if activate and not inter_active:
mpp.ion() # turn interactive mode on
elif not activate and inter_active:
mpp.ioff() # turn interactive mode off
# main routine
AUTHOR = "Copyright 2021 The Regents of the University of Colorado"
SEPARBOLD = 79*'='
SEPARNORM = 79*'-'
SOFTWARE = "ParamAP"
VERSION = "version 1.3," # (2021-05-29)
WORKDIR = SOFTWARE # working directory for parameterization
print('{0:^79}'.format(SEPARBOLD) + os.linesep)
GREETER = '{0:<{w0}}{1:<{w1}}{2:<{w2}}'.format(SOFTWARE, VERSION, AUTHOR, w0=len(SOFTWARE)+1, w1=len(VERSION)+1, w2=len(AUTHOR)+1)
INTERMEDIATELINE = '{0:}'.format("University of Colorado, Anschutz Medical Campus")
DISCLAIMER = "ParamAP is distributed in the hope that it will be useful, but it comes without\nany guarantee or warranty. This program is free software; you can redistribute\nit and/or modify it under the terms of the GNU General Public License:"
URL = "https://www.gnu.org/licenses/gpl-2.0.en.html"
print('{0:^79}'.format(GREETER) + os.linesep)
print('{0:^79}'.format(INTERMEDIATELINE) + os.linesep)
print('{0:^79}'.format(DISCLAIMER) + os.linesep)
print('{0:^79}'.format(URL) + os.linesep)
print('{0:^79}'.format(SEPARBOLD) + os.linesep)
# customize use case
AUTORUN = askboolean("Use automatic mode?", False)
SERIES = askboolean("Run time series analysis?", False)
APMODE = askboolean("Analyze action potentials?", True)
print('{0:^79}'.format(SEPARNORM))
# set up working directory
WORKPATH = os.path.abspath(WORKDIR)
if not os.path.exists(WORKPATH):
os.mkdir(WORKPATH)
print("FOLDER:\t" + WORKPATH + "\n")
FILE = 0 # file
EXTENSION = askext(dlabel="custom file type", dext='atf') # file extension used to filter files in working directory
if SERIES:
AVG_FRAME = askvalue(dlabel="analysis frame time", dval=5000.0, dunit=' ms') # time interval for series analysis (ms)
ATFFILES = getfiles(path=WORKDIR, pattern=("*." + EXTENSION))
for ATFFILE in ATFFILES: # iterate through files
name = os.path.splitext(os.path.split(ATFFILE)[1])[0]
print('{0:^79}'.format(SEPARNORM))
print("FILE:\t" + str(name) + os.linesep)
AP_AMP = 50.0 # minimum acceptable ap amplitude (mV)
AP_HWD = 250.0 # maximum acceptable ap half width (ms)
AP_MAX = 50.0 # maximum acceptable ap value (mV)
AP_MIN = -10.0 # minimum acceptable ap value (mV)
MDP_MAX = -50.0 # maximum acceptable mdp value (mV)
MDP_MIN = -90.0 # minimum acceptable mdp value (mV)
WM_DER = 1.0 # window multiplier for derivative filtering
WM_MAX = 4.0 # window multiplier for maximum detection
WM_MIN = 16.0 # window multiplier for minimum detection
# read file raw data
sys.stdout.write(">> READING... ")
sys.stdout.flush()
RAW_XY, UNITS, UNIT_X, UNIT_Y = readfile(ATFFILE)
if not UNITS: # missing or incomplete units from header
print("\n")
UNIT_X = askunit(dlabel="X-axis unit", daxis="X", dunit=UNIT_X)
UNIT_Y = askunit(dlabel="Y-axis unit", daxis="Y", dunit=UNIT_Y)
sys.stdout.write(1*"\t")
TOMS = 1000.0 if UNIT_X == "s" else 1.0
RAW_XY[0] *= TOMS # full X-axis, UNIT_X = "ms"
raw_x = RAW_XY[0] # partial X-axis for time series analysis
TOMV = 1000.0 if UNIT_Y == "V" else 1.0
RAW_XY[1] *= TOMV # full Y-axis, UNIT_Y = "mV"
raw_y = RAW_XY[1] # partial Y-axis for time series analysis
runavg = getrunavg(RAW_XY[0]) # used for filtering and peak detection
ipg_t = RAW_XY[0][1]-RAW_XY[0][0] # time increment for interpolation grid
if not APMODE: # avoid noise artifacts in beat detection mode
runavg = 10.0*runavg+1
WM_MAX *= 1.5
WM_MIN = WM_MAX
avg_start = RAW_XY[0][0] # interval start for averaging
avg_stop = RAW_XY[0][-1] # interval stop for averaging
sys.stdout.write(8*"\t" + " [OK]\n")
sys.stdout.flush()
while True: # repeat data analysis for current file
STARTPDF = True # overwrite existing file
SEGMENT = 0.0
while True: # time series analysis
try:
# create raw data plot
sys.stdout.write(">> PLOTTING... ")
sys.stdout.flush()
mpp_setup(title="Raw data: " + name, xlabel='Time (ms)', ylabel='Voltage (mV)')
mpp.plot(raw_x, raw_y, '0.75') # raw data (grey line)
if STARTPDF:
pdf_file = mpbp.PdfPages(os.path.join(WORKDIR, name + ".pdf"), keep_empty=False) # multi-pdf file
STARTPDF = False # append existing file
mpp.tight_layout() # avoid label overlaps
if SEGMENT == 0.0:
mpp.savefig(pdf_file, format='pdf', dpi=600) # save before .show()!
sys.stdout.write(8*"\t" + " [OK]\n")
sys.stdout.flush()
if not AUTORUN:
toggle_mpp_imode(activate=True)
mpp.show()
# set parameters for averaging
sys.stdout.write(">> SETTING... ")
sys.stdout.flush()
if not AUTORUN:
print("\n")
if SEGMENT == 0.0: # initialize values
avg_start = askvalue(dlabel="analysis start time", dval=avg_start, dunit=' ms')
avg_stop = askvalue(dlabel="analysis stop time", dval=avg_stop, dunit=' ms')
AP_MAX = askvalue(dlabel="upper limit for maxima", dval=AP_MAX, dunit=' mV')
AP_MIN = askvalue(dlabel="lower limit for maxima", dval=AP_MIN, dunit=' mV')
MDP_MAX = askvalue(dlabel="upper limit for minima", dval=MDP_MAX, dunit=' mV')
MDP_MIN = askvalue(dlabel="lower limit for minima", dval=MDP_MIN, dunit=' mV')
if APMODE:
AP_HWD = askvalue(dlabel="maximum peak half width", dval=AP_HWD, dunit=' ms')
AP_AMP = askvalue(dlabel="minimum peak amplitude", dval=AP_AMP, dunit=' mV')
runavg = askvalue(dlabel="running average window size", dval=runavg, dunit='', dtype='int')
WM_DER = askvalue(dlabel="window multiplier for derivative", dval=WM_DER, dunit='')
WM_MAX = askvalue(dlabel="window multiplier for maxima", dval=WM_MAX, dunit='')
WM_MIN = askvalue(dlabel="window multiplier for minima", dval=WM_MIN, dunit='')
mpp.clf() # clear canvas
if SEGMENT == 0.0: # set first frame
tmp_start = avg_start + (SEGMENT*AVG_FRAME if SERIES else 0.0)
tmp_stop = (tmp_start + AVG_FRAME) if SERIES else avg_stop
raw_i = np.argwhere((RAW_XY[0] >= tmp_start) & (RAW_XY[0] <= tmp_stop)).ravel()
raw_x = RAW_XY[0][raw_i[0]:raw_i[-1]+1]
raw_y = RAW_XY[1][raw_i[0]:raw_i[-1]+1]
sys.stdout.write(("" if AUTORUN else 1*"\t") + 8*"\t" + " [OK]\n")
sys.stdout.flush()
# filter noise of raw data with Savitzky-Golay
sys.stdout.write(">> FILTERING... ")
sys.stdout.flush()
rawf_y = sp_sig.savgol_filter(raw_y, runavg, POLYNOM, mode='nearest')
sys.stdout.write(7*"\t" + " [OK]\n")
sys.stdout.flush()
# detect extrema in filtered raw data
sys.stdout.write(">> SEARCHING... ")
sys.stdout.flush()
if AUTORUN: # use unrestricted dataset (slower)
# detect maxima in filtered raw data
tmpavg = int(round(WM_MAX*runavg)) if int(round(WM_MAX*runavg)) % 2 else int(round(WM_MAX*runavg))+1
rawfmax_iii = np.asarray(sp_sig.argrelmax(rawf_y, order=tmpavg)).ravel() # unfiltered maxima
rawfmax_x = raw_x[rawfmax_iii]
rawfmax_y = rawf_y[rawfmax_iii]
# detect minima in filtered raw data
tmpavg = int(round(WM_MIN*runavg)) if int(round(WM_MIN*runavg)) % 2 else int(round(WM_MIN*runavg))+1
rawfmin_iii = np.asarray(sp_sig.argrelmin(rawf_y, order=tmpavg)).ravel() # unfiltered minima
rawfmin_x = raw_x[rawfmin_iii]
rawfmin_y = rawf_y[rawfmin_iii]
sys.stdout.write(7*"\t" + " [OK]\n")
sys.stdout.flush()
else: # use restricted dataset (faster)
# detect maxima in filtered raw data
tmpmax_x = raw_x[np.intersect1d(np.argwhere(rawf_y >= AP_MIN), np.argwhere(rawf_y <= AP_MAX))]
tmpmax_y = rawf_y[np.intersect1d(np.argwhere(rawf_y >= AP_MIN), np.argwhere(rawf_y <= AP_MAX))]
tmpavg = int(round(WM_MAX*runavg)) if int(round(WM_MAX*runavg)) % 2 else int(round(WM_MAX*runavg))+1
rawfmax_iii = np.asarray(sp_sig.argrelmax(tmpmax_y, order=tmpavg)).ravel() # unfiltered maxima
rawfmax_ii = np.asarray(np.where(np.in1d(raw_x.ravel(), np.intersect1d(raw_x, tmpmax_x[rawfmax_iii]).ravel()).reshape(raw_x.shape))).ravel() # back to full dataset
rawfmax_x = raw_x[rawfmax_ii]
rawfmax_y = rawf_y[rawfmax_ii]
# detect minima in filtered raw data
tmpmin_x = raw_x[np.intersect1d(np.argwhere(rawf_y >= MDP_MIN), np.argwhere(rawf_y <= MDP_MAX))]
tmpmin_y = rawf_y[np.intersect1d(np.argwhere(rawf_y >= MDP_MIN), np.argwhere(rawf_y <= MDP_MAX))]
tmpavg = int(round(WM_MIN*runavg)) if int(round(WM_MIN*runavg)) % 2 else int(round(WM_MIN*runavg))+1
rawfmin_iii = np.asarray(sp_sig.argrelmin(tmpmin_y, order=tmpavg)).ravel() # unfiltered minima
rawfmin_ii = np.asarray(np.where(np.in1d(raw_x.ravel(), np.intersect1d(raw_x, tmpmin_x[rawfmin_iii]).ravel()).reshape(raw_x.shape))).ravel()
rawfmin_x = raw_x[rawfmin_ii]
rawfmin_y = rawf_y[rawfmin_ii]
sys.stdout.write(7*"\t" + " [OK]\n")
sys.stdout.flush()
# analyze and reduce extrema in filtered raw data
sys.stdout.write(">> REDUCING... ")
sys.stdout.flush()
rawfmax_m = np.mean(rawfmax_y) # rough estimate due to assignment errors
rawfmin_m = np.mean(rawfmin_y)
rawfmaxmin_m = (rawfmax_m + rawfmin_m) / 2.0 # center between unreduced maxima and minima within limits (may differ from average of AVGMAX and AVGMIN)
if AUTORUN: # estimate range for reduction of extrema
# reduce maxima from unrestricted dataset
rawfmax_ii = np.argwhere(rawfmax_y >= rawfmaxmin_m).ravel() # use center to discriminate between maxima and minima
rawfmax_x = rawfmax_x[rawfmax_ii]
rawfmax_y = rawfmax_y[rawfmax_ii]
rawfmax_std = np.std(rawfmax_y, ddof=1) # standard deviation from the (estimated) arithmetic mean
AP_MAX = np.mean(rawfmax_y) + 4.0 * rawfmax_std # 99% confidence interval
AP_MIN = np.mean(rawfmax_y) - 4.0 * rawfmax_std
rawfmax_ii = functools.reduce(np.intersect1d, (rawfmax_iii, np.argwhere(rawf_y >= AP_MIN), np.argwhere(rawf_y <= AP_MAX)))
rawfmax_x = raw_x[rawfmax_ii]
rawfmax_y = rawf_y[rawfmax_ii]
# reduce minima from unrestricted dataset
rawfmin_ii = np.argwhere(rawfmin_y <= rawfmaxmin_m)
rawfmin_x = rawfmin_x[rawfmin_ii].ravel()
rawfmin_y = rawfmin_y[rawfmin_ii].ravel()
rawfmin_std = np.std(rawfmin_y, ddof=1)
MDP_MAX = np.mean(rawfmin_y) + 4.0 * rawfmin_std
MDP_MIN = np.mean(rawfmin_y) - 4.0 * rawfmin_std
rawfmin_ii = functools.reduce(np.intersect1d, (rawfmin_iii, np.argwhere(rawf_y >= MDP_MIN), np.argwhere(rawf_y <= MDP_MAX)))
rawfmin_x = raw_x[rawfmin_ii]
rawfmin_y = rawf_y[rawfmin_ii]
if APMODE: # check extrema for consistency - reduce maxima
badmax_ii = np.zeros(0)
badmin_ii = np.zeros(0)
rawfmin_i, badmax_ii = getneighbors(rawfmax_ii, rawfmin_ii, raw_x, rawf_y, AP_HWD, AP_AMP, bads=True)
rawfmax_i = np.delete(rawfmax_ii, badmax_ii)
rawfmin_i = rawfmin_i.astype(int) # casting required for indexing
# check extrema for boundary violations - reduce maxima and minima
while True: # rough check, assignment happens later
if rawfmax_i[0] < rawfmin_i[0]: # starts with a maximum
rawfmax_i = rawfmax_i[1:]
continue
elif rawfmin_i[1] < rawfmax_i[0]: # starts with two minima
rawfmin_i = rawfmin_i[1:]
continue
elif rawfmax_i[-1] > rawfmin_i[-1]: # ends with a maximum
rawfmax_i = rawfmax_i[0:-1]
continue
elif rawfmin_i[-2] > rawfmax_i[-1]: # ends with two minima
rawfmin_i = rawfmin_i[0:-1]
continue
else:
break
rawfmax_x = raw_x[rawfmax_i] # filtered and extracted maxima
rawfmax_y = rawf_y[rawfmax_i]
# assign minima to corresponding maxima - reduce minima
minmaxmin = np.asarray([3*[0] for i in range(rawfmax_i.size)]) # [[min_left_index, max_index, min_right_index], ...]
rawfmin_kdt = sp_spat.KDTree(list(zip(rawfmin_i, np.zeros(rawfmin_i.size))))
i = 0 # index
for max_i in rawfmax_i:
MIN_LEFT, MIN_RIGHT = False, False
minmaxmin[i][1] = max_i
for order_i in rawfmin_kdt.query([max_i, 0.0], k=None)[1]:
min_i = rawfmin_i[order_i]
if not MIN_LEFT and (min_i < max_i):
minmaxmin[i][0] = min_i
MIN_LEFT = True
elif not MIN_RIGHT and (min_i > max_i):
minmaxmin[i][2] = min_i
MIN_RIGHT = True
i += 1
rawfmin_i = np.unique(minmaxmin[:, [0, 2]].ravel())
rawfmin_x = raw_x[rawfmin_i] # filtered and extracted minima
rawfmin_y = rawf_y[rawfmin_i]
# find largest distance between left minima and maxima
ipg_hwl, ipg_tmp = 0.0, 0.0
for min_l, max_c in minmaxmin[:, [0, 1]]:
ipg_tmp = raw_x[max_c] - raw_x[min_l]
if ipg_tmp > ipg_hwl:
ipg_hwl = ipg_tmp
# find largest distance between right minima and maxima
ipg_hwr, ipg_tmp = 0.0, 0.0
for max_c, min_r in minmaxmin[:, [1, 2]]:
ipg_tmp = raw_x[min_r] - raw_x[max_c]
if ipg_tmp > ipg_hwr:
ipg_hwr = ipg_tmp
else: # beating rate
rawfmax_x = raw_x[rawfmax_ii] # pre-filtered maxima
rawfmax_y = rawf_y[rawfmax_ii]
rawfmin_x = raw_x[rawfmin_ii] # pre-filtered minima
rawfmin_y = rawf_y[rawfmin_ii]
rawfmax_m = np.mean(rawfmax_y) # refined estimate due to exlusion (ap_mode)
rawfmin_m = np.mean(rawfmin_y)
if rawfmax_y.size == 0: # no APs detected
raise Warning
else: # two or more APs
frate = 60000.0*(rawfmax_y.size/(rawfmax_x[-1]-rawfmax_x[0])) if rawfmax_y.size > 1 else float('nan') # AP firing rate (FR) [1/min]
sys.stdout.write(8*"\t" + " [OK]\n")
sys.stdout.flush()
# create extrema plot
sys.stdout.write(">> PLOTTING... ")
sys.stdout.flush()
mpp_setup(title="Extrema: " + name, xlabel='Time (ms)', ylabel='Voltage (mV)')
mpp.plot([raw_x[0], raw_x[-1]], [0.0, 0.0], '0.85') # X-Axis (grey line)
mpp.plot([raw_x[0], raw_x[-1]], [rawfmaxmin_m, rawfmaxmin_m], 'k--') # center between unfiltered maxima and unfiltered minima, i.e. not between AVGMAX and AVGMIN (black dashed line)
mpp.plot(raw_x, raw_y, '0.50', raw_x, rawf_y, 'r') # raw data and averaged data (grey, red line)
mpp.plot([raw_x[0], raw_x[-1]], [AP_MAX, AP_MAX], 'b') # upper limit for maxima (blue dotted line)
mpp.plot([raw_x[0], raw_x[-1]], [AP_MIN, AP_MIN], 'b:') # lower limit for maxima (blue dotted line)
mpp.plot([rawfmax_x, rawfmax_x], [AP_MIN, AP_MAX], 'b') # accepted maxima (blue line)
mpp.plot([raw_x[0], raw_x[-1]], [MDP_MIN, MDP_MIN], 'g') # lower limit for minima (green line)
mpp.plot([raw_x[0], raw_x[-1]], [MDP_MAX, MDP_MAX], 'g:') # upper limit for minima (green dotted line)
mpp.plot([rawfmin_x, rawfmin_x], [MDP_MIN, MDP_MAX], 'g') # accepted minima (green line)
mpp.plot([rawfmax_x[0], rawfmax_x[-1]], [rawfmax_m, rawfmax_m], 'k') # average of maxima, time interval used for firing rate count (black line)
mpp.plot([rawfmin_x[0], rawfmin_x[-1]], [rawfmin_m, rawfmin_m], 'k') # average of minima (black line)
mpp.plot(raw_x[rawfmax_ii], rawf_y[rawfmax_ii], 'bo') # maxima (blue dots)
mpp.plot(raw_x[rawfmin_ii], rawf_y[rawfmin_ii], 'go') # minima (green dots)
mpp.figtext(0.12, 0.91, "{0:<s} {1:<.4G}".format("AVGMAX (mV):", rawfmax_m), ha='left', va='center')
mpp.figtext(0.12, 0.88, "{0:<s} {1:<.4G}".format("FR (AP/min):", frate), ha='left', va='center')
mpp.figtext(0.12, 0.85, "{0:<s} {1:<.4G}".format("AVGMIN (mV):", rawfmin_m), ha='left', va='center')
mpp.tight_layout()
sys.stdout.write(8*"\t" + " [OK]\n")
sys.stdout.flush()
sys.stdout.write(">> SAVING... ")
sys.stdout.flush()
mpp.savefig(pdf_file, format='pdf', dpi=600)
sys.stdout.write(8*"\t" + " [OK]\n")
sys.stdout.flush()
if not AUTORUN:
toggle_mpp_imode(activate=False)
mpp.show()
mpp.clf()
if APMODE:
# slice raw data segments by minima and align by maxima
sys.stdout.write(">> AVERAGING... ")
sys.stdout.flush()
# align ap segments by maxima, extend and average ap segments
ipg_max = float((1.0+APEXT)*ipg_hwr)
ipg_min = -float((1.0+APEXT)*ipg_hwl)
avg_x = np.arange(ipg_min, ipg_max, ipg_t, dtype='float64') # interpolation grid
avgxsize = avg_x.size
avg_y = np.zeros(avgxsize, dtype='float64') # ap average array
mpp.subplot2grid((4, 1), (0, 0), rowspan=3) # upper subplot
TIMESTAMP = "[" + str(round(tmp_start, 2)) + "ms-" + str(round(tmp_stop, 2)) + "ms]"
mpp_setup(title='Analysis: ' + name + ' ' + TIMESTAMP, xlabel='Time (ms)', ylabel='Voltage (mV)')
mpp.plot([avg_x[0], avg_x[-1]], [0.0, 0.0], '0.85') # X-axis
N = 0 # current maximum
for min_l, max_c, min_r in minmaxmin: # slicing of ap segments, extend ap parts if possible
minext_l = int(min_l - APEXT*(max_c - min_l)) # use int for index slicing
minext_r = int(min_r + APEXT*(min_r - max_c))
# prepare ap SEGMENT
tmp_x = np.asarray(raw_x[:] - raw_x[max_c]) # align by maximum
tmp_y = np.interp(avg_x, tmp_x, raw_y[:])
# average ap segments
if N == 0: # first average
avg_y = np.copy(tmp_y)
else: # all other averages
i = 0 # array index
NW = (1.0/(N+1.0)) # new data weight
pw = (N/(N+1.0)) # previous data weight
for y in np.nditer(avg_y, op_flags=['readwrite']):
y[...] = pw*y + NW*tmp_y[i] # integrate raw data into averaged data
i += 1
N += 1
mpp.plot(avg_x, tmp_y, '0.75') # plot aligned raw data segments
sys.stdout.write("\t\t\t\t\t\t\t [OK]\n")
sys.stdout.flush()
# analyze AP parameters with given criteria
sys.stdout.write(">> ANALYZING... ")
sys.stdout.flush()
# filter noise of averaged data with Savitzky-Golay
avgf_y = sp_sig.savgol_filter(avg_y, runavg, POLYNOM, mode='nearest')
# detect "Peak potential: Maximum potential of AP" (PP) (mV)
avgfmax_i = np.argwhere(avg_x == 0.0) # data point for maximum centered
if not avgfmax_i.size > 0: # data point for maximum left or right of center
tmpavg = int(round(WM_MAX*runavg)) if int(round(WM_MAX*runavg)) % 2 else int(round(WM_MAX*runavg))+1
avgfmax_ii = np.asarray(sp_sig.argrelmax(avgf_y, order=tmpavg)).ravel() # find all maxima
avgfmax_i = avgfmax_ii[np.argmin(np.abs(avg_x[avgfmax_ii] - 0.0))] # return the maximum closest to X = 0.0
avgfmax_x = avg_x[avgfmax_i]
avgfmax_y = avgf_y[avgfmax_i]
pp_y = float(avgfmax_y)
pp = pp_y
# detect and reduce (several) minima in filtered average data,
tmpavg = int(round(WM_MIN*runavg)) if int(round(WM_MIN*runavg)) % 2 else int(round(WM_MIN*runavg))+1
avgfmin_ii = np.asarray(sp_sig.argrelmin(avgf_y, order=tmpavg)).ravel() # find all minima
avgfmin_i = getneighbors(np.asarray([avgfmax_i]), avgfmin_ii, avg_x, avgf_y, AP_HWD, AP_AMP)
avgfmin_x = avg_x[avgfmin_i]
avgfmin_y = avgf_y[avgfmin_i]
# determine "Maximum diastolic potential 1: Minimum potential preceding PP" (MDP1) (mV)
mdp1_i = avgfmin_i[0]
mdp1_x = avg_x[mdp1_i]
mdp1_y = avgf_y[mdp1_i]
mdp1 = mdp1_y
# determine "Maximum diastolic potential 2: Minimum potential following PP" (MDP2) (mV)
mdp2_i = avgfmin_i[-1]
mdp2_x = avg_x[mdp2_i]
mdp2_y = avgf_y[mdp2_i]
mdp2 = mdp2_y
# determine "Cycle length: Time interval MDP1-MDP2" (CL) (ms)
cl = float(mdp2_x - mdp1_x)
# determine "Action potential amplitude: Potential difference of PP minus MDP2" (APA) (mV)
apa = pp - mdp2
# determine "AP duration 30: Time interval at 30% of maximum repolarization" (APD30) (ms)
apd30_l = (pp - 0.30*apa) # threshold value
apd30_i = functools.reduce(np.intersect1d, (np.argwhere(avgf_y > apd30_l), np.argwhere(avg_x >= mdp1_x), np.argwhere(avg_x <= mdp2_x)))
apd30_x = (avg_x[apd30_i[0]-1], avg_x[apd30_i[-1]+1]) # equal to or smaller than apd30_l
apd30_y = (avgf_y[apd30_i[0]-1], avgf_y[apd30_i[-1]+1])
apd30 = float(apd30_x[-1] - apd30_x[0])
# determine "AP duration 50: Time interval at 50% of maximum repolarization" (APD50) (ms)
apd50_l = (pp - 0.50*apa) # threshold value
apd50_i = functools.reduce(np.intersect1d, (np.argwhere(avgf_y > apd50_l), np.argwhere(avg_x >= mdp1_x), np.argwhere(avg_x <= mdp2_x)))
apd50_x = (avg_x[apd50_i[0]-1], avg_x[apd50_i[-1]+1]) # equal to or smaller than apd50_l
apd50_y = (avgf_y[apd50_i[0]-1], avgf_y[apd50_i[-1]+1])
apd50 = float(apd50_x[-1] - apd50_x[0])
# determine "AP duration 90: Time interval at 90% of maximum repolarization" (APD90) (ms)
apd90_l = pp - 0.90*apa
apd90_i = functools.reduce(np.intersect1d, (np.argwhere(avgf_y > apd90_l), np.argwhere(avg_x >= mdp1_x), np.argwhere(avg_x <= mdp2_x)))
apd90_x = (avg_x[apd90_i[0]-1], avg_x[apd90_i[-1]+1]) # equal to or smaller than apd90_l
apd90_y = (avgf_y[apd90_i[0]-1], avgf_y[apd90_i[-1]+1])
apd90 = float(apd90_x[-1] - apd90_x[0])
# calculate derivative of averaged data (mV/ms)
avgfg_y = np.ediff1d(avgf_y) # dY/1, differences between values
avgfg_y = np.insert(avgfg_y, 0, avgfg_y[0]) # preserve array size
avgfg_y = avgfg_y / ipg_t # dY/dX, differences per increment
# filter derivative of averaged data
tmpavg = int(round(WM_DER*runavg)) if int(round(WM_DER*runavg)) % 2 else int(round(WM_DER*runavg))+1
avgfgf_y = sp_sig.savgol_filter(avgfg_y, tmpavg, POLYNOM, mode='nearest')
# determine "Maximum upstroke velocity: Maximum of derivative between MDP1 and PP" (MUV) (mV/ms)
tmpavg = int(round(WM_MAX*runavg)) if int(round(WM_MAX*runavg)) % 2 else int(round(WM_MAX*runavg))+1
avgfgfmax_ii = functools.reduce(np.intersect1d, (sp_sig.argrelmax(avgfgf_y, order=tmpavg), np.argwhere(avg_x >= mdp1_x), np.argwhere(avg_x <= avgfmax_x)))
avgfgfmax_i = getneighbors(np.asarray([avgfmax_i]), avgfgfmax_ii, avg_x, avgfgf_y)[0] # avoid errors from large ap part extensions
avgfgfmax_x = avg_x[avgfgfmax_i]
avgfgfmax_y = avgfgf_y[avgfgfmax_i]
muv = float(avgfgfmax_y)
# determine "Maximum repolarization rate: Minimum of derivative between PP and MDP2" (MRR) (mV/ms)
tmpavg = int(round(WM_MIN*runavg)) if int(round(WM_MIN*runavg)) % 2 else int(round(WM_MIN*runavg))+1
avgfgfmin_ii = functools.reduce(np.intersect1d, (sp_sig.argrelmin(avgfgf_y, order=tmpavg), np.argwhere(avg_x >= avgfmax_x), np.argwhere(avg_x <= mdp2_x)))
avgfgfmin_i = getneighbors(np.asarray([apd90_i[-1]+1]), avgfgfmin_ii, avg_x, avgfgf_y)[0] # mrr or trr
# determine "Transient repolarization rate: Second minimum of derivative between PP and MDP2 after PP, if distinct from MRR" (TRR) (mV/ms)
avgfgfmin_i = np.append(avgfgfmin_i, getneighbors(np.asarray([avgfgfmax_i]), avgfgfmin_ii, avg_x, avgfgf_y)[1]) # trr only
avgfgfmin_x = avg_x[avgfgfmin_i]
avgfgfmin_y = avgfgf_y[avgfgfmin_i]
mrr = float(avgfgf_y[avgfgfmin_i][0])
if avgfgfmin_i[0] == avgfgfmin_i[1]: # no trr
WITH_TRR = False
trr = float("nan")
else:
WITH_TRR = True
trr = float(avgfgf_y[avgfgfmin_i][1]) # True
# approximate diastolic duration in filtered derivative
da_i, da_x, da_y, da_m, da_n, da_r = getbestlinearfit(avg_x, avgfgf_y, mdp1_x, apd90_x[0], 10, 90, 1, 40) # get a baseline for the derivative before exceeding the threshold
# determine "Threshold potential: Potential separating DD and APD." (THR) (mV)
thr_i = functools.reduce(np.intersect1d, (np.argwhere(avgfgf_y >= ((da_m*avg_x + da_n) + 0.5)), np.argwhere(avg_x >= avg_x[da_i[-1]]), np.argwhere(avg_x <= apd50_x[0])))[0].astype(int) # determine baseline-corrected threshold level
thr_x = avg_x[thr_i]
thr_y = avgf_y[thr_i]
thr = float(thr_y)
# determine "Late AP duration"
lapd_l = thr
lapd_i = functools.reduce(np.intersect1d, (np.argwhere(avgf_y > lapd_l), np.argwhere(avg_x >= avgfmax_x), np.argwhere(avg_x <= mdp2_x)))
lapd_x = (0.5*(avg_x[lapd_i[-1]]+avg_x[lapd_i[-1]+1])) # equal to or smaller than lapd_l
lapd_y = (0.5*(avgf_y[lapd_i[-1]]+avgf_y[lapd_i[-1]+1]))
lapd = float(mdp2_x - lapd_x)
# determine "Early diastolic duration: Time from MDP1 to end of linear fit for DDR" (EDD) (ms)
edd_i, edd_x, edd_y, edd_m, edd_n, edd_r = getbestlinearfit(avg_x, avgf_y, mdp1_x, thr_x, 10, 50, 1, 20) # fit EDD within the threshold level determined earlier
edd = float(edd_x[-1]-mdp1_x)
# determine "Diastolic depolarization rate: Potential change rate at end of EDD" (DDR) (mV/ms)
ddr = float(edd_m) # or: np.mean(avgfgf_y[edd_i])
# determine "Diastolic duration: EDD plus LDD" (DD) (ms)
dd = float(thr_x - mdp1_x)
# determine "Late diastolic duration: Time from end of linear fit for DDR to THR" (LDD) (ms)
ldd = float(thr_x - edd_x[-1])
# determine "Action potential duration: Time between THR and MDP2" (APD) (ms)
apd = float(mdp2_x - thr_x)
sys.stdout.write("\t\t\t\t\t\t\t [OK]\n")
sys.stdout.flush()
# create analysis plot
sys.stdout.write(">> PLOTTING... ") # the X-axis and the individual segments are already plotted during averaging
sys.stdout.flush()
mpp.plot([mdp1_x, thr_x], [mdp1_y, mdp1_y], 'k-.') # DD (black dashed/dotted line)
mpp.plot([thr_x, mdp2_x], [mdp2_y, mdp2_y], 'k') # APD (black line)
mpp.plot([lapd_x, mdp2_x], [lapd_y, lapd_y], 'r-.') # LAPD (black line)
mpp.plot([apd50_x[0], apd50_x[1]], [apd50_y[1], apd50_y[1]], 'k') # APD50 (black line)
mpp.plot([apd90_x[0], apd90_x[1]], [apd90_y[1], apd90_y[1]], 'k') # APD90 (black line)
mpp.plot([mdp1_x, mdp1_x], [mdp1_y, 0.0], 'k:') # MDP1 indicator (black dotted line)
mpp.plot([mdp2_x, mdp2_x], [mdp2_y, 0.0], 'k:') # MDP2 indicator (black dotted line)
mpp.plot([avgfgfmax_x, avgfgfmax_x], [mdp2_y, avgf_y[avgfgfmax_i]], 'k:') # MUV indicator (black dotted line)
mpp.plot([avgfgfmin_x[0], avgfgfmin_x[0]], [mdp2_y, avgf_y[avgfgfmin_i[0]]], 'k:') # MRR indicator (black dotted line)
mpp.plot([edd_x[-1], edd_x[-1]], [mdp2_y, 0.0], 'k:') # EDD/LDD separator (black dashed line)
mpp.plot([thr_x, thr_x], [thr_y, 0.0], 'k:') # DD/APD upper separator (black dotted line)
mpp.plot([thr_x, thr_x], [mdp2_y, thr_y], 'k:') # DD/APD lower separator (black dotted line)
mpp.plot(avg_x, avg_y, 'k', avg_x, avgf_y, 'r') # averaged data and filtered averaged data (black, red lines)
mpp.plot(avg_x[edd_i], avgf_y[edd_i], 'g') # best linear fit segment for DDR (green line)
mpp.plot(avg_x, (edd_m*avg_x + edd_n), 'k--') # DDR (black dashed line)
mpp.plot([edd_x[-1]], [edd_y[-1]], 'ko') # EDD-LDD separator (black dot)
mpp.plot(lapd_x, lapd_y, 'ro') # LAPD (black dot)
mpp.plot([apd50_x[1]], [apd50_y[1]], 'ko') # APD50 (black dots)
mpp.plot(apd90_x[1], apd90_y[1], 'ko') # APD90 (black dots)
mpp.plot(thr_x, avgf_y[thr_i], 'ro') # THR (red dot)
mpp.plot(avgfgfmax_x, avgf_y[avgfgfmax_i], 'wo') # MUV (white dot)
mpp.plot(avgfgfmin_x[0], avgf_y[avgfgfmin_i[0]], 'wo') # MRR (white dot)
if WITH_TRR:
mpp.plot([avgfgfmin_x[1], avgfgfmin_x[1]], [mdp2_y, avgf_y[avgfgfmin_i[1]]], 'k:') # trr indicator (black dotted line)
mpp.plot(avgfgfmin_x[1], avgf_y[avgfgfmin_i[1]], 'wo') # trr (white dot)
mpp.plot(avgfmax_x, pp_y, 'bo') # PP (blue dot)
mpp.plot(avgfmin_x, avgfmin_y, 'go') # MDP1, MDP2 (green dots)
mpp.figtext(0.12, 0.91, "{0:<s} {1:<.4G}".format("APs (#):", rawfmax_y.size), ha='left', va='center')
mpp.figtext(0.12, 0.88, "{0:<s} {1:<.4G}".format("FR (AP/min):", frate), ha='left', va='center')
mpp.figtext(0.12, 0.85, "{0:<s} {1:<.4G}".format("CL (ms):", cl), ha='left', va='center')
mpp.figtext(0.12, 0.82, "{0:<s} {1:<.4G}".format("DD (ms):", dd), ha='left', va='center')
mpp.figtext(0.12, 0.79, "{0:<s} {1:<.4G}".format("EDD (ms):", edd), ha='left', va='center')
mpp.figtext(0.12, 0.76, "{0:<s} {1:<.4G}".format("LDD (ms):", ldd), ha='left', va='center')
mpp.figtext(0.12, 0.73, "{0:<s} {1:<.4G}".format("APD (ms):", apd), ha='left', va='center')
mpp.figtext(0.12, 0.70, "{0:<s} {1:<.4G}".format("LAPD (ms):", lapd), ha='left', va='center')
mpp.figtext(0.12, 0.67, "{0:<s} {1:<.4G}".format("APD50 (ms):", apd50), ha='left', va='center')
mpp.figtext(0.12, 0.64, "{0:<s} {1:<.4G}".format("APD90 (ms):", apd90), ha='left', va='center')
mpp.figtext(0.12, 0.61, "{0:<s} {1:<.4G}".format("MDP1 (mV):", mdp1), ha='left', va='center')
mpp.figtext(0.12, 0.58, "{0:<s} {1:<.4G}".format("MDP2 (mV):", mdp2), ha='left', va='center')
mpp.figtext(0.12, 0.55, "{0:<s} {1:<.4G}".format("THR (mV):", thr), ha='left', va='center')
mpp.figtext(0.12, 0.52, "{0:<s} {1:<.4G}".format("PP (mV):", pp), ha='left', va='center')
mpp.figtext(0.12, 0.49, "{0:<s} {1:<.4G}".format("APA (mV):", apa), ha='left', va='center')
mpp.figtext(0.12, 0.46, "{0:<s} {1:<.4G}".format("DDR (mV/ms):", ddr), ha='left', va='center')
mpp.figtext(0.12, 0.43, "{0:<s} {1:<.4G}".format("MUV (mV/ms):", muv), ha='left', va='center')
mpp.figtext(0.12, 0.40, "{0:<s} {1:<.4G}".format("TRR (mV/ms):", trr), ha='left', va='center')
mpp.figtext(0.12, 0.37, "{0:<s} {1:<.4G}".format("MRR (mV/ms):", mrr), ha='left', va='center')
mpp.subplot2grid((4, 1), (3, 0)) # lower subplot
mpp_setup(title="", xlabel='Time (ms)', ylabel='(mV/ms)')
mpp.plot([avg_x[0], avg_x[-1]], [0.0, 0.0], '0.85') # x axis
mpp.plot([avgfgfmin_x[0], avgfgfmin_x[0]], [avgfgfmin_y[0], avgfgfmax_y], 'k:') # MRR indicator (black dotted line)
if WITH_TRR:
mpp.plot([avgfgfmin_x[1], avgfgfmin_x[1]], [avgfgfmin_y[1], avgfgfmax_y], 'k:') # TRR indicator (black dotted line)
mpp.plot([thr_x, thr_x], [avgfgf_y[thr_i], avgfgfmax_y], 'k:') # THR indicator (black dotted line)
mpp.plot(avg_x, avgfg_y, 'c', avg_x, avgfgf_y, 'm') # derivative and filtered derivative
mpp.plot(avg_x[da_i], avgfgf_y[da_i], 'g') # best linear fit SEGMENT for THR (green line)
mpp.plot(avg_x, (da_m*avg_x + da_n), 'k--') # best linear fit for THR (black dashed line)
mpp.plot(thr_x, avgfgf_y[thr_i], 'ro') # THR (red dot)
mpp.plot(avgfgfmax_x, avgfgfmax_y, 'bo') # derivative maximum (blue dot)
mpp.plot(avgfgfmin_x, avgfgfmin_y, 'go') # derivative minima (green dots)
sys.stdout.write(8*"\t" + " [OK]\n")
sys.stdout.flush()
# data summary
sys.stdout.write(">> SAVING... ")
sys.stdout.flush()
avg_file = os.path.join(WORKDIR, name + "_" + TIMESTAMP + "_avg.dat")
UHEADER = "" +\
"Analysis start time: " + 4*"\t" + str(tmp_start) + " ms\n" + \
"Analysis stop time:" + 4*"\t" + str(tmp_stop) + " ms\n" + \
"Upper limit for maxima:" + 3*"\t" + str(AP_MAX) + " mV\n" + \
"Lower limit for maxima:" + 3*"\t" + str(AP_MIN) + " mV\n" + \
"Upper limit for minima:" + 3*"\t" + str(MDP_MAX) + " mV\n" + \
"Lower limit for minima:" + 3*"\t" + str(MDP_MIN) + " mV\n" + \
"Maximum peak half width:" + 3*"\t" + str(AP_HWD) + " ms\n" + \
"Minimum peak amplitude:" + 3*"\t" + str(AP_AMP) + " mV\n" + \
"Running average window size:" + 2*"\t" + str(runavg) + "\n" + \
"Window multiplier for derivative:" + "\t" + str(WM_DER) + "\n" + \
"Window multiplier for maxima:" + 2*"\t" + str(WM_MAX) + "\n" + \
"Window multiplier for minima:" + 2*"\t" + str(WM_MIN) + "\n" + \
"Time (ms)" + "\t" + "Averaged signal (mV)" + "\t" + "Filtered average (mV)"
np.savetxt(avg_file, np.column_stack((avg_x, avg_y, avgf_y)), fmt='%e', delimiter='\t', header=UHEADER)
mpp.tight_layout()
mpp.savefig(pdf_file, format='pdf', dpi=600)
sum_file = os.path.join(WORKDIR, "ParamAP.log")
NEWFILE = not bool(os.path.exists(sum_file))
with open(sum_file, 'a') as targetfile: # append file
if NEWFILE: # write header
targetfile.write("{0:s}\t{1:s}\t{2:s}\t{3:s}\t{4:s}\t{5:s}\t{6:s}\t{7:s}\t{8:s}\t{9:s}\t{10:s}\t{11:s}\t{12:s}\t{13:s}\t{14:s}\t{15:s}\t{16:s}\t{17:s}\t{18:s}\t{19:s}\t{20:s}\t{21:s}".format(
"File ( )", "Start (ms)", "Stop (ms)", "APs (#)", "FR (AP/min)", "CL (ms)", "DD (ms)", "EDD (ms)", "LDD (ms)", "APD (ms)", "LAPD (ms)", "APD50 (ms)", "APD90 (ms)", "MDP1 (mV)", "MDP2 (mV)", "THR (mV)", "PP (mV)", "APA (mV)", "DDR (mV/ms)", "MUV (mV/ms)", "TRR (mV/ms)", "MRR (mV/ms)") + "\n")
targetfile.write("{0:s}\t{1:4G}\t{2:4G}\t{3:4G}\t{4:4G}\t{5:4G}\t{6:4G}\t{7:4G}\t{8:4G}\t{9:4G}\t{10:4G}\t{11:4G}\t{12:4G}\t{13:4G}\t{14:4G}\t{15:4G}\t{16:4G}\t{17:4G}\t{18:4G}\t{19:4G}\t{20:4G}\t{21:4G}".format(
name, tmp_start, tmp_stop, rawfmax_y.size, frate, cl, dd, edd, ldd, apd, lapd, apd50, apd90, mdp1, mdp2, thr, pp, apa, ddr, muv, trr, mrr) + "\n")
targetfile.flush()
sys.stdout.write(8*"\t" + " [OK]\n")
sys.stdout.flush()
if not AUTORUN:
mpp.show()
except IndexError as ierr: # check running average and window multiplier
sys.stdout.write("\n" + 9*"\t" + " [ER]")
print("\r ## Run failed. Detection of extrema or threshold failed.")
except PermissionError as perr: # file already opened or storage read-only
sys.stdout.write("\n" + 9*"\t" + " [ER]")
print("\r ## Run failed. File access denied by system.")
except Warning as werr: # increase averaging window time
sys.stdout.write("\n" + 9*"\t" + " [ER]")
print("\r ## Run failed. Identification of action potentials failed.")
except Exception as uerr: # unknown
sys.stdout.write("\n" + 9*"\t" + " [UN]")
print("\r ## Run failed. Error was: {0}".format(uerr) + ".")
except KeyboardInterrupt as kerr: # user canceled this file
sys.stdout.write("\n" + 9*"\t" + " [KO]")
print("\r ## Run skipped. Canceled by user.")
if SERIES: # check for next frame
if tmp_stop + AVG_FRAME <= avg_stop:
SEGMENT += 1.0
tmp_start = avg_start + SEGMENT*AVG_FRAME # prepare next frame for preview
tmp_stop = tmp_start + AVG_FRAME
raw_i = np.argwhere((RAW_XY[0] >= tmp_start) & (RAW_XY[0] <= tmp_stop)).ravel()
raw_x = RAW_XY[0][raw_i[0]:raw_i[-1]+1]
raw_y = RAW_XY[1][raw_i[0]:raw_i[-1]+1]
print()
print("RUN:\t" + str(int(SEGMENT + 1)) + "/" + str(math.floor((avg_stop-avg_start)/AVG_FRAME)))
print()
else: # not enough data left in file
break
else: # no time series analysis
break
if not AUTORUN: # check for next file
print()
NEXTFILE = askboolean("Continue with next file?", True)
if NEXTFILE:
break
else: # re-run current file
raw_x = RAW_XY[0] # recover original rawdata
raw_y = RAW_XY[1]
continue
else: # autorun
break
# housekeeping after each file
FILE += 1
sys.stdout.write(">> CLEANING... ")
sys.stdout.flush()
pdf_file.close() # close multi-pdf file and remove if empty
mpp.clf() # clear canvas
gc.collect() # start garbage collection to prevent memory fragmentation
sys.stdout.write(8*"\t" + " [OK]\n")
sys.stdout.flush()
# print summary
print('{0:^79}'.format(SEPARBOLD))
SUMMARY = "End of run: " + str(FILE) + str(" files" if FILE != 1 else " file") + " processed."
print('{0:^79}'.format(SUMMARY))
print('{0:^79}'.format(SEPARBOLD) + os.linesep)
WAIT = input("Press ENTER to end this program.")