From e9c872fc0e255fedbe92fd85510fd430781886d9 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Fri, 3 May 2024 14:54:27 +0200 Subject: [PATCH 01/35] Improved MUcv_gui() Resizable window, optimised RAM memory usage, allows to specify a custom csv field delimiter. --- openhdemg/library/muap.py | 128 ++++++++++++++++++++++++-------------- 1 file changed, 82 insertions(+), 46 deletions(-) diff --git a/openhdemg/library/muap.py b/openhdemg/library/muap.py index b0a2067..ceccb19 100644 --- a/openhdemg/library/muap.py +++ b/openhdemg/library/muap.py @@ -25,7 +25,6 @@ import tkinter as tk from tkinter import ttk from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg -import pyperclip def diff(sorted_rawemg): @@ -1587,6 +1586,10 @@ class MUcv_gui(): The STA is calculated over all the firings. muaps_timewindow : int, default 50 Timewindow to compute ST MUAPs in milliseconds. + figsize : list, default [20, 15] + Size of the initial MUAPs figure in centimeters [width, height]. + csv_separator : str, default "\t" + The field delimiter used to create the .csv copied to the clipboard. See also -------- @@ -1621,6 +1624,7 @@ def __init__( n_firings=[0, 50], muaps_timewindow=50, figsize=[25, 20], + csv_separator="\t", ): # On start, compute the necessary information self.emgfile = emgfile @@ -1633,6 +1637,7 @@ def __init__( ) self.sta_xcc = xcc_sta(self.st) self.figsize = figsize + self.csv_separator = csv_separator # After that, set up the GUI self.root = tk.Tk() @@ -1647,20 +1652,33 @@ def __init__( ) self.root.iconbitmap(iconpath) - # Create main frame, assign structure and minimum spacing - self.frm = ttk.Frame(self.root, padding=15) - # Assign grid structure - self.frm.grid() + # Create outer frames, assign structure and minimum spacing + # Left + left_frm = tk.Frame(self.root, padx=2, pady=2) + left_frm.pack(side=tk.LEFT, expand=True, fill="both") + # Right + right_frm = tk.Frame(self.root, padx=4, pady=4) + right_frm.pack(side=tk.TOP, anchor="nw", expand=True, fill="y") + + # Create inner frames + # Top left + top_left_frm = tk.Frame(left_frm, padx=2, pady=2) + top_left_frm.pack(side=tk.TOP, anchor="nw", fill="x") + # Bottom left + self.bottom_left_frm = tk.Frame(left_frm, padx=2, pady=2) + self.bottom_left_frm.pack( + side=tk.TOP, anchor="nw", expand=True, fill="both", + ) # Label MU number combobox - munumber_label = ttk.Label(self.frm, text="MU number", width=15) + munumber_label = ttk.Label(top_left_frm, text="MU number", width=15) munumber_label.grid(row=0, column=0, columnspan=1, sticky=tk.W) # Create a combobox to change MU self.all_mus = list(range(emgfile["NUMBER_OF_MUS"])) self.selectmu_cb = ttk.Combobox( - self.frm, + top_left_frm, textvariable=tk.StringVar(), values=self.all_mus, state='readonly', @@ -1676,91 +1694,101 @@ def __init__( lambda event: self.gui_plot(), ) - # Add 2 empty columns - emp0 = ttk.Label(self.frm, text="", width=15) - emp0.grid(row=0, column=1, columnspan=1, sticky=tk.W) - emp1 = ttk.Label(self.frm, text="", width=15) - emp1.grid(row=0, column=2, columnspan=1, sticky=tk.W) + # Add empty column + emp = ttk.Label(top_left_frm, text="", width=15) + emp.grid(row=0, column=1, columnspan=1, sticky=tk.W) # Create the widgets to calculate CV # Label and combobox to select the matrix column - col_label = ttk.Label(self.frm, text="Column", width=15) - col_label.grid(row=0, column=3, columnspan=1, sticky=tk.W) + col_label = ttk.Label(top_left_frm, text="Column", width=15) + col_label.grid(row=0, column=2, columnspan=1, sticky=tk.W) self.columns = list(self.st[0].keys()) self.col_cb = ttk.Combobox( - self.frm, + top_left_frm, textvariable=tk.StringVar(), values=self.columns, state='readonly', width=15, ) - self.col_cb.grid(row=1, column=3, columnspan=1, sticky=tk.W) + self.col_cb.grid(row=1, column=2, columnspan=1, sticky=tk.W) self.col_cb.current(0) # Label and combobox to select the matrix channels self.rows = list(range(len(list(self.st[0][self.columns[0]].columns)))) - start_label = ttk.Label(self.frm, text="From row", width=15) - start_label.grid(row=0, column=4, columnspan=1, sticky=tk.W) + start_label = ttk.Label(top_left_frm, text="From row", width=15) + start_label.grid(row=0, column=3, columnspan=1, sticky=tk.W) self.start_cb = ttk.Combobox( - self.frm, + top_left_frm, textvariable=tk.StringVar(), values=self.rows, state='readonly', width=15, ) - self.start_cb.grid(row=1, column=4, columnspan=1, sticky=tk.W) + self.start_cb.grid(row=1, column=3, columnspan=1, sticky=tk.W) self.start_cb.current(0) - self.stop_label = ttk.Label(self.frm, text="To row", width=15) - self.stop_label.grid(row=0, column=5, columnspan=1, sticky=tk.W) + self.stop_label = ttk.Label(top_left_frm, text="To row", width=15) + self.stop_label.grid(row=0, column=4, columnspan=1, sticky=tk.W) self.stop_cb = ttk.Combobox( - self.frm, + top_left_frm, textvariable=tk.StringVar(), values=self.rows, state='readonly', width=15, ) - self.stop_cb.grid(row=1, column=5, columnspan=1, sticky=tk.W) + self.stop_cb.grid(row=1, column=4, columnspan=1, sticky=tk.W) self.stop_cb.current(max(self.rows)) # Button to start CV estimation self.ied = emgfile["IED"] self.fsamp = emgfile["FSAMP"] button_est = ttk.Button( - self.frm, + top_left_frm, text="Estimate", command=self.compute_cv, width=15, ) - button_est.grid(row=1, column=6, columnspan=1, sticky="we") + button_est.grid(row=1, column=5, columnspan=1, sticky="we") - # Add empty column - self.emp2 = ttk.Label(self.frm, text="", width=5) - self.emp2.grid(row=0, column=7, columnspan=1, sticky=tk.W) + # Configure column weights + for c in range(6): + if c == 1: + top_left_frm.columnconfigure(c, weight=20) + else: + top_left_frm.columnconfigure(c, weight=1) - # Add text frame to show the results (only CV and RMS) - self.res_df = pd.DataFrame( - data=0, - index=self.all_mus, - columns=["CV", "RMS", "XCC", "Column", "From_Row", "To_Row"], - ) - self.textbox = tk.Text(self.frm, width=20) - self.textbox.grid(row=2, column=8, sticky="ns") - self.textbox.insert('1.0', self.res_df.loc[:, ["CV", "RMS"]].to_string()) + # Align the right_frm + alignment_label = ttk.Label(right_frm, text="", width=15) + alignment_label.pack(side=tk.TOP, fill="x") # Create a button to copy the dataframe to clipboard copy_btn = ttk.Button( - self.frm, + right_frm, text="Copy results", command=self.copy_to_clipboard, - width=20, + width=15, + ) + copy_btn.pack(side=tk.TOP, fill="x", pady=(0, 5)) + + # Add text frame to show the results (only CV and RMS) + self.res_df = pd.DataFrame( + data=0.00, + index=self.all_mus, + columns=["CV", "RMS", "XCC", "Column", "From_Row", "To_Row"], + ) + self.textbox = tk.Text(right_frm, width=25) + self.textbox.pack(side=tk.TOP, expand=True, fill="y") + self.textbox.insert( + '1.0', + self.res_df.loc[:, ["CV", "RMS", "XCC"]].to_string( + float_format="{:.2f}".format + ), ) - copy_btn.grid(row=1, column=8, columnspan=1, sticky="we") # Plot MU 0 while opening the GUI, # this will move the GUI in the background ??. @@ -1788,16 +1816,22 @@ def gui_plot(self): figsize=self.figsize, ) + # If canvas already exists, destroy it + if hasattr(self, 'canvas'): + self.canvas.get_tk_widget().destroy() + # Place the figure in the GUI - canvas = FigureCanvasTkAgg(fig, master=self.frm) - canvas.draw() - canvas.get_tk_widget().grid(row=2, column=0, columnspan=7, sticky="we") + self.canvas = FigureCanvasTkAgg(fig, master=self.bottom_left_frm) + self.canvas.draw_idle() # Await resizing + self.canvas.get_tk_widget().pack( + expand=True, fill="both", padx=0, pady=0, + ) plt.close() def copy_to_clipboard(self): # Copy the dataframe to clipboard in csv format. - pyperclip.copy(self.res_df.to_csv(index=False, sep='\t')) + self.res_df.to_clipboard(excel=True, sep=self.csv_separator) # Define functions for cv estimation def compute_cv(self): @@ -1839,5 +1873,7 @@ def compute_cv(self): self.textbox.replace( '1.0', 'end', - self.res_df.loc[:, ["CV", "RMS"]].round(3).to_string(), + self.res_df.loc[:, ["CV", "RMS", "XCC"]].to_string( + float_format="{:.2f}".format + ), ) From 3fdf2a32d8d4ebaecdbe0f40e9f627b24459bbb6 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Fri, 3 May 2024 14:55:17 +0200 Subject: [PATCH 02/35] Dependency management Removed required dependency pyperclip --- requirements.txt | Bin 366 -> 332 bytes setup.py | 3 +-- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index 9a01b8d606ec3a603c5068c29a0756efd0b62f51..592a3f4a6409389a1a08fb69df942cf21eb1fe70 100644 GIT binary patch delta 10 RcmaFIbcSidx5*NWVE`KN1V;b> delta 41 scmX@Z^p0u5H~9jFN`?Z4RE8pkWQH7uOd#792n`wZ7%YI;XyWHc0QiUsiU0rr diff --git a/setup.py b/setup.py index 853a955..a5dff5f 100644 --- a/setup.py +++ b/setup.py @@ -17,11 +17,10 @@ "openpyxl==3.1.2", "pandas==2.0.3", "pandastable==0.13.1", - "pyperclip==1.8.2", - # "pre-commit==3.7.0", "scipy==1.11.3", "seaborn==0.13.0", "joblib==1.3.2", + # "pre-commit==3.7.0", ] PACKAGES = [ From a69ae6b66b894f810a92b653f787db8483fca0cd Mon Sep 17 00:00:00 2001 From: JABeauchamp <103655220+JABeauchamp@users.noreply.github.com> Date: Sat, 4 May 2024 13:14:41 -0500 Subject: [PATCH 03/35] Initial SVR fits and delta-F --- CODEOWNERS | 5 + openhdemg/library/__init__.py | 2 + openhdemg/library/pic.py | 247 ++++++++++++++++++++++++++++++++++ openhdemg/library/tools.py | 147 ++++++++++++++++++++ requirements.txt | Bin 332 -> 362 bytes 5 files changed, 401 insertions(+) create mode 100644 openhdemg/library/pic.py diff --git a/CODEOWNERS b/CODEOWNERS index bba5f83..c69ef4d 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -7,3 +7,8 @@ # When there are changes in the GUI folder, @GiacomoValliPhD # and @PaulRitsche will be requested for review. /openhdemg/gui/ @GiacomoValliPhD @PaulRitsche + +# When there are changes in the pic module, @GiacomoValliPhD +# and @JABeauchamp will be requested for review. +/openhdemg/library/pic @GiacomoValliPhD @JABeauchamp +/openhdemg/library/analysis/comput_svr @GiacomoValliPhD @JABeauchamp diff --git a/openhdemg/library/__init__.py b/openhdemg/library/__init__.py index 1df624e..58a4fc1 100644 --- a/openhdemg/library/__init__.py +++ b/openhdemg/library/__init__.py @@ -7,6 +7,7 @@ "otbelectrodes", "muap", "info", + "pic" ] from openhdemg.library.openfiles import ( @@ -30,3 +31,4 @@ from openhdemg.library.electrodes import * from openhdemg.library.muap import * from openhdemg.library.info import * +from openhdemg.library.pic import * diff --git a/openhdemg/library/pic.py b/openhdemg/library/pic.py new file mode 100644 index 0000000..ef7aa90 --- /dev/null +++ b/openhdemg/library/pic.py @@ -0,0 +1,247 @@ +""" +This module contains all the functions used to quantify and analyze MU persistent inward currents. +Currently includes delta F. +""" + +import pandas as pd +import numpy as np +from openhdemg.library.tools import compute_svr +from itertools import combinations + + +def compute_deltaf( + emgfile, + smoothfits, + average_method="test_unit_average", + normalization="False", + recruitment_difference_cutoff=1.0, + corr_cutoff=0.7, + controlunitmodulation_cutoff=0.5, + clean="True" + ): + + """ + Conducts a paired motor unit analysis, quantifying delta F bettwen the supplied collection of motor units. + Origional framework for deltaF provided in Gorassini et. al., 2002 : https://journals.physiology.org/doi/full/10.1152/jn.00024.2001 + James (Drew) Beauchamp + - + Parameters + ---------- + emgfile : dict + The dictionary containing the emgfile. + + smoothfits : list of arrays + Smoothed discharge rate estimates. + Each array: motor unit discharge rate x samples aligned in time; instances of non-firing = NaN, please + Your choice of smoothing. See compute_svr gen_svr for example. + + average_method : str, default "test_unit_average" + The method for test MU deltaF value. More to be added TODO + + "test_unit_average" + The average across all possible control units + + normalization : str, {"False", "ctrl_max_desc"}, default "False" + "ctrl_max_desc" + Whether to normalize deltaF values to control unit descending range during test unit firing + See Skarabot et. al., 2023 : https://www.biorxiv.org/content/10.1101/2023.10.16.562612v1 + + recruitment_difference_cutoff : float, default 1 + An exlusion criteria corresponding to the necessary difference between control and test MU recruitement in seconds + corr_cutoff : float, (0,1), default 0.7 + An exclusion criteria corresponding to the correlation between control and test unit discharge rate. + controlunitmodulation_cutoff : float, default 0.5 + An exclusion criteria corresponding to the necessary modulation of control unit discharge rate during test unit firing in Hz. + clean : str, default "True" + To remove values that do not meet exclusion criteria + - + Returns + ------- + dF : pd.DataFrame + A pd.DataFrame containing deltaF values and corresponding MU number + + Warnings + -------- + TODO + + Example + ----- + import openhdemg.library as emg + import numpy as np + import matplotlib.pyplot as plt + + # Load the sample file + emgfile = emg.emg_from_samplefile() + + # Sort MUs based on recruitment order + emgfile = emg.sort_mus(emgfile=emgfile) + + # quantify svr fits + svrfits = emg.tools.compute_svr(emgfile) + + # quantify delta F + dF = emg.pic.compute_deltaf(emgfile=emgfile, smoothfits=svrfits["gensvr"]) + + dF + MU dF + 0 0 NaN + 1 1 NaN + 2 2 NaN + 3 3 1.838382 + 4 4 2.709522 + + # for all possible combinations, not test unit average, MU in this case is pairs (reporter,test) + dF2 = emg.pic.compute_deltaf(emgfile=emgfile, smoothfits=svrfits["gensvr"],average_method='all') + + dF2 + MU dF + 0 (0, 1) NaN + 1 (0, 2) NaN + 2 (0, 3) 2.127461 + 3 (0, 4) NaN + 4 (1, 2) NaN + 5 (1, 3) 1.549303 + 6 (1, 4) NaN + 7 (2, 3) NaN + 8 (2, 4) NaN + 9 (3, 4) 2.709522 + + + """ + # svrfits = compute_SVR(emgfile) # use smooth svr fits + + dfret_ret = [] + mucombo_ret = np.empty(0, int) + if emgfile["NUMBER_OF_MUS"] < 2: # if less than 2 MUs, can not quantify deltaF + dfret_ret = np.nan + mucombo_ret = np.nan*np.ones(1, 2) + else: + combs = combinations(range(emgfile["NUMBER_OF_MUS"]), 2) + # combinations of MUs + # TODO if units are nonconsecutive + + # init + r_ret = [] + dfret = [] + testmu = [] + ctrl_mod = [] + mucombo = [] + rcrt_diff = [] + controlmu = [] + for mucomb in list(combs): # for all possible combinations of MUs + mu1_id, mu2_id = mucomb[0], mucomb[1] # extract possible MU combinations (a unique MU pair) + mucombo.append((mu1_id, mu2_id)) # track current MU combination + + # first MU firings, recruitment, and decrecruitment + mu1_times = np.where(emgfile["BINARY_MUS_FIRING"][mu1_id] == 1)[0] + mu1_rcrt, mu1_drcrt = mu1_times[1], mu1_times[-1] # skip first since idr is defined on second + + # second MU firings, recruitment, and decrecruitment + mu2_times = np.where(emgfile["BINARY_MUS_FIRING"][mu2_id]==1)[0] + mu2_rcrt, mu2_drcrt = mu2_times[1], mu2_times[-1] # skip first since idr is defined on second + + # region of MU overlap + muoverlap = range(max(mu1_rcrt, mu2_rcrt), min(mu1_drcrt, mu2_drcrt)) + + # if MUs do not overlapt by more than two or more samples + if len(muoverlap) < 2: + dfret = np.append(dfret, np.nan) + r_ret = np.append(r_ret, np.nan) + rcrt_diff = np.append(rcrt_diff, np.nan) + ctrl_mod = np.append(ctrl_mod, np.nan) + continue # TODO test + + # corr between units - not always necessary, can be set to 0 when desired + r = pd.DataFrame(zip(smoothfits[mu1_id][muoverlap], smoothfits[mu2_id][muoverlap])).corr() + r_ret = np.append(r_ret, r[0][1]) + + # recruitment diff, necessary to ensure PICs are activated in control unit + rcrt_diff = np.append(rcrt_diff, np.abs(mu1_rcrt-mu2_rcrt)/emgfile["FSAMP"]) + if mu1_rcrt < mu2_rcrt: + controlU = 1 # MU 1 is control unit, 2 is test unit + + if mu1_drcrt < mu2_drcrt: # if control (reporter) unit is not on for entirety of test unit, set last firing to control unit + mu2_drcrt = mu1_drcrt + # this may understimate PICs, other methods can be employed + # delta F: change in control MU discharge rate between test unit recruitment and derecruitment + df = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu2_drcrt] + + # control unit discharge rate modulation while test unit is firing + ctrl_mod = np.append(ctrl_mod, np.nanmax(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)])-np.nanmin(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)])) + + if normalization == "False": + dfret = np.append(dfret, df) + elif normalization == "ctrl_max_desc": # normalize deltaF values to control unit descending range during test unit firing + k = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu1_drcrt] + dfret = np.append(dfret, df/k) + + elif mu1_rcrt > mu2_rcrt: + controlU = 2 # MU 2 is control unit, 1 is test unit + if mu1_drcrt > mu2_drcrt: # if control (reporter) unit is not on for entirety of test unit, set last firing to control unit + mu1_drcrt = mu2_drcrt + # this may understimate PICs, other methods can be employed + # delta F: change in control MU discharge rate between test unit recruitment and derecruitment + df = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu1_drcrt] + + # control unit discharge rate modulation while test unit is firing + ctrl_mod = np.append(ctrl_mod, np.nanmax(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)])-np.nanmin(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)])) + + if normalization == "False": + dfret = np.append(dfret, df) + elif normalization == "ctrl_max_desc": # normalize deltaF values to control unit descending range during test unit firing + k = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu2_drcrt] + dfret = np.append(dfret, df/k) + + elif mu1_rcrt == mu2_rcrt: + if mu1_drcrt > mu2_drcrt: + controlU = 1 # MU 1 is control unit, 2 is test unit + # delta F: change in control MU discharge rate between test unit recruitment and derecruitment + df = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu2_drcrt] + + # control unit discharge rate modulation while test unit is firing + ctrl_mod = np.append(ctrl_mod, np.nanmax(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)])-np.nanmin(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)])) + + if normalization == "False": + dfret = np.append(dfret, df) + elif normalization == "ctrl_max_desc": # normalize deltaF values to control unit descending range during test unit firing + k = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu1_drcrt] + dfret = np.append(dfret, df/k) + else: + controlU = 2 # MU 2 is control unit, 1 is test unit + # delta F: change in control MU discharge rate between test unit recruitment and derecruitment + df = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu1_drcrt] + + # control unit discharge rate modulation while test unit is firing + ctrl_mod = np.append(ctrl_mod, np.nanmax(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)])-np.nanmin(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)])) + + if normalization == "False": + dfret = np.append(dfret, df) + elif normalization == "ctrl_max_desc": # normalize deltaF values to control unit descending range during test unit firing + k = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu2_drcrt] + dfret = np.append(dfret, df/k) + + # collect which MUs were control vs test + controlmu.append(mucombo[-1][controlU-1]) + testmu.append(mucombo[-1][1-controlU//2]) + + if clean == "True": # remove values that dont meet exclusion criteria + rcrt_diff_bin = rcrt_diff > recruitment_difference_cutoff + corr_bin = r_ret > corr_cutoff + ctrl_mod_bin = ctrl_mod > controlunitmodulation_cutoff + clns = np.asarray([rcrt_diff_bin & corr_bin & ctrl_mod_bin]) + dfret[~clns[0]] = np.nan + + if average_method == "test_unit_average": # average across all control units + for ii in range(emgfile["NUMBER_OF_MUS"]): + clean_indices = [index for (index, item) in enumerate(testmu) if item == ii] + if np.isnan(dfret[clean_indices]).all(): + dfret_ret = np.append(dfret_ret, np.nan) + else: + dfret_ret = np.append(dfret_ret, np.nanmean(dfret[clean_indices])) + mucombo_ret = np.append(mucombo_ret, int(ii)) + else: # return all values and corresponding combinations + dfret_ret = dfret + mucombo_ret = mucombo + + dF = pd.DataFrame({'MU': mucombo_ret, 'dF': dfret_ret}) + return dF diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index 541dc43..7ed3ede 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -13,6 +13,9 @@ import warnings from openhdemg.library.mathtools import compute_sil +from sklearn.svm import SVR +from scipy.stats import iqr + def showselect(emgfile, how="ref_signal", title="", titlesize=12, nclic=2): """ @@ -1195,3 +1198,147 @@ def compute_rfd( rfd = rfd * conversion_val return rfd + + +def compute_svr( + emgfile, + gammain=1/1.6, + regparam=1/0.370, + endpointweights_numpulses=5, + endpointweights_magnitude=5, + discontfiring_dur=1.0 + ): + """ + Fit MU Discharge rates with Support Vector Regression, nonlinear regression + Provides smooth and continous estimates of discharge rate usefull for quantification and visualization + Suggested hyperparameters and framework from Beauchamp et. al., 2022 + https://doi.org/10.1088/1741-2552/ac4594 + James (Drew) Beauchamp + + Parameters + ---------- + emgfile : dict + The dictionary containing the emgfile. + gammain : float, default 1/1.6 + The kernel coefficient. + regparam : float, default 1/0.370 + The regularization parameter, must be positive. + endpointweights_numpulses : int, default 5 + Number of discharge instances at the start and end of MU firing to apply a weighting coefficient. + endpointweights_magnitude : int, default 5 + The scaling factor applied to the number of pulses provided by endpointweights_numpulses. + The scaling is applied to the regularization parameter, per sample. + Larger values force the classifier to put more emphasis on the number of discharge instances at the start and end of firing provided by endpointweights_numpulses + discontfiring_dur : int, default 1 + Duration of time in seconds that defines an instnance of discontinuous firing. SVR fits will not be returned at points of discontinuity. + + Returns + ------- + svrfits : pd.DataFrame + A pd.DataFrame containing the smooth/continous MU discharge rates and corresponding time vetors + + Warnings + -------- + TODO + + Example + ----- + import openhdemg.library as emg + import numpy as np + import matplotlib.pyplot as plt + + # Load the sample file + emgfile = emg.emg_from_samplefile() + + # Sort MUs based on recruitment order + emgfile = emg.sort_mus(emgfile=emgfile) + + # quantify svr fits + svrfits = emg.tools.compute_svr(emgfile) + + # quick plot showing results + scl = 12 # offset MUs for viz + idr = emg.compute_idr(emgfile) + for ii in range(len(svrfits["svrfit"])): + xtmp = np.transpose([idr[ii].timesec[1:]]) + ytmp = idr[ii].idr[1:].to_numpy() + plt.scatter(xtmp, ytmp+scl*(ii), color='darkorange') + plt.plot(svrfits["svrtime"][ii], svrfits["svrfit"][ii]+scl*(ii), color='cornflowerblue') + plt.show() + """ + + # TODO input checking and edge cases + idr = compute_idr(emgfile) # calc IDR + + svrfit_acm = [] + svrtime_acm = [] + gensvr_acm = [] + for mu in range(len(idr)): # for all MUs + # train the model on the data + xtmp = np.transpose([idr[mu].timesec[1:]]) # time vector, removing first element + ytmp = idr[mu].idr[1:].to_numpy() # discharge rates, removing first element, since DR has been assigned to second pulse + xdiff = idr[mu].diff_mupulses[1:].values # time between discharges, will use for discontinuity calc + mup = np.array(idr[mu].mupulses[1:]) # motor unit pulses, samples + + # Defining weight vector. A scaling applied to the regularization parameter, per sample. + smpwht = np.ones(len(ytmp)) + smpwht[0:endpointweights_numpulses-1] = endpointweights_magnitude + smpwht[(len(ytmp)-(endpointweights_numpulses-1)):len(ytmp)] = endpointweights_magnitude + + # create an SVR model with a gausian kernel and supplied hyperparams + # origional hyperparameters from Beauchamp et. al., 2022 https://doi.org/10.1088/1741-2552/ac4594 + svr = SVR(kernel='rbf', gamma=gammain, C=np.abs(regparam), epsilon=iqr(ytmp)/11) + svr.fit(xtmp, ytmp, sample_weight=smpwht) + + # Defining prediction vector + # TODO need to add custom range + predind = np.arange(mup[0], mup[-1]+1) # from the second firing to the end of firing, in samples + predtime = (predind/emgfile["FSAMP"]).reshape(-1, 1) # in time (s) + gen_svr = np.nan*np.ones(emgfile["EMG_LENGTH"]) # initialize nan vector for tracking fits aligned in time. usefull for later quant metrics + + # check for discontinous firing + bkpnt = mup[np.where((xdiff > (discontfiring_dur * emgfile["FSAMP"])))[0]] + + # make predictions on the data + if bkpnt.size > 0: # if there is a point of discontinuity + if bkpnt[0] == mup[0]: + smoothfit = np.nan*np.ones(1) + newtm = np.nan*np.ones(1) + bkpnt = bkpnt[1:] + else: + tmptm = predtime[0:np.where((bkpnt[0] >= predind[0:-1]) & (bkpnt[0] < predind[1:]))[0][0]] # break up time vector for first continous range of firing + smoothfit = svr.predict(tmptm) # predict with svr model + newtm = tmptm # track new time vector + + tmpind = predind[0:np.where((bkpnt[0] >= predind[0:-1]) & ( bkpnt[0]= predind[0:-1]) & (bkpnt[ii] < predind[1:]))[0][0] # current index of discontinuity + nextind = np.where((bkpnt[ii+1] >= predind[0:-1]) & (bkpnt[ii+1] <= predind[1:]))[0][0] # next index of discontinuity + + curmup = np.where(mup == bkpnt[ii])[0][0] # MU firing before discontinuity + curind_nmup = np.where((mup[curmup+1] >= predind[0:-1]) & (mup[curmup+1] <= predind[1:]))[0][0] # MU firing after discontinuity + + if curind_nmup >= nextind: # if the next discontinuity is the next MU firing, nan fill + # edge case NEED TO CHECK THE GREATER THAN CASE>> WHY TODO + smoothfit = np.append(smoothfit, np.nan*np.ones(1)) + newtm = np.append(newtm, np.nan*np.ones(1)) + else: # fit next continuous region of firing + smoothfit = np.append(smoothfit, np.nan*np.ones(len(predtime[curind:curind_nmup])-2)) + smoothfit = np.append(smoothfit, svr.predict(predtime[curind_nmup:nextind])) + newtm = np.append(newtm, np.nan*np.ones(len(predtime[curind:curind_nmup])-2)) + newtm = np.append(newtm, predtime[curind_nmup:nextind]) + gen_svr[predind[curind_nmup:nextind]] = svr.predict(predtime[curind_nmup:nextind]) + else: + smoothfit = svr.predict(predtime) + newtm = predtime + gen_svr[predind] = smoothfit + + svrfit_acm.append(smoothfit.copy()) + svrtime_acm.append(newtm.copy()) + gensvr_acm.append(gen_svr.copy()) # append fits, new time vect, time aligned fits + + svrfits = ({"svrfit": svrfit_acm, "svrtime": svrtime_acm, "gensvr": gensvr_acm}) + + return svrfits diff --git a/requirements.txt b/requirements.txt index 592a3f4a6409389a1a08fb69df942cf21eb1fe70..b7f523cf87884e551ecdf411858b09dd3a465b18 100644 GIT binary patch delta 38 qcmX@Z^onVN52G9xLoq`(Lk>eKLn1>FLmq=I5E?S*F_-|cAp-!imIpil delta 7 OcmaFGbcShz4 Date: Wed, 29 May 2024 11:26:15 +0200 Subject: [PATCH 04/35] Added Tracking_gui() class --- openhdemg/library/muap.py | 466 +++++++++++++++++++++++++++++++++++--- 1 file changed, 430 insertions(+), 36 deletions(-) diff --git a/openhdemg/library/muap.py b/openhdemg/library/muap.py index ceccb19..2323fb3 100644 --- a/openhdemg/library/muap.py +++ b/openhdemg/library/muap.py @@ -12,7 +12,7 @@ mle_cv_est, ) from openhdemg.library.electrodes import sort_rawemg -from openhdemg.library.plotemg import plot_muaps, plot_muaps_for_cv +from openhdemg.library.plotemg import plot_idr, plot_muaps, plot_muaps_for_cv from scipy import signal import matplotlib.pyplot as plt from functools import reduce @@ -21,7 +21,6 @@ from joblib import Parallel, delayed import copy import os -import warnings import tkinter as tk from tkinter import ttk from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg @@ -687,7 +686,7 @@ def align_by_xcorr(sta_mu1, sta_mu2, finalduration=0.5): no_nan_sta2 = df2.dropna(axis=1, inplace=False) # Compute 2dxcorr to identify a common lag/delay - normxcorr_df, normxcorr_max = norm_twod_xcorr( + normxcorr_df, _ = norm_twod_xcorr( no_nan_sta1, no_nan_sta2, mode="same" ) @@ -697,6 +696,7 @@ def align_by_xcorr(sta_mu1, sta_mu2, finalduration=0.5): ) normxcorr_df = normxcorr_df.set_index(corr_lags) lag = normxcorr_df.idxmax().median() # First signal compared to second + # TODO Can mean provide better results? allow choice? # Be sure that the lag/delay does not exceed values suitable for the final # expected duration. @@ -736,6 +736,8 @@ def align_by_xcorr(sta_mu1, sta_mu2, finalduration=0.5): # TODO update examples for code="None" +# TODO change parallel processing, allow to specify if parallel +# TODO add GUI # This function exploits parallel processing: # - align and xcorr are processed in parallel @@ -756,7 +758,7 @@ def tracking( exclude_belowthreshold=True, filter=True, show=False, -): +): # TODO add gui example """ Track MUs across two files comparing the MUAPs' shape and distribution. @@ -859,15 +861,16 @@ def tracking( - sta : computes the STA of every MUs. - norm_twod_xcorr : normalised 2-dimensional cross-correlation of STAs of two MUs. + - Tracking_gui : GUI for the visual inspection and manual selection of the + tracking results (directly callable from the tracking function). - remove_duplicates_between : remove duplicated MUs across two different files based on STA. Notes ----- - Parallel processing can improve performances by 5-10 times compared to + Parallel processing can significantly improve performances compared to serial processing. In this function, parallel processing has been - implemented for the tasks involving 2-dimensional cross-correlation, and - plotting. This might change in future releases. + implemented for the tasks involving 2-dimensional cross-correlation. Examples -------- @@ -1019,7 +1022,7 @@ def parallel(mu_file1): # Loop all the MUs of file 1 df1.dropna(axis=1, inplace=True) df2, _ = unpack_sta(aligned_sta2) df2.dropna(axis=1, inplace=True) - normxcorr_df, normxcorr_max = norm_twod_xcorr( + _, normxcorr_max = norm_twod_xcorr( df1, df2, mode="full" ) @@ -1042,13 +1045,12 @@ def parallel(mu_file1): # Loop all the MUs of file 1 res = Parallel(n_jobs=-1)( delayed(parallel)(mu_file1) for mu_file1 in range(emgfile1["NUMBER_OF_MUS"]) - ) + ) # TODO n_jobs -1 remove joblib and allow user to select number of cores t1 = time.time() print(f"\nTime of tracking parallel processing: {round(t1-t0, 2)} Sec\n") # Convert res to pd.DataFrame - tracking_res = [] for pos, i in enumerate(res): if pos == 0: tracking_res = pd.DataFrame(i) @@ -1095,45 +1097,436 @@ def parallel(mu_file1): # Loop all the MUs of file 1 pd.reset_option("display.max_rows") # Plot the MUs pairs - if show: - def parallel(ind): # Function for the parallel execution of plotting - if tracking_res["XCC"].loc[ind] >= threshold: - # Align STA - if not isinstance(custom_muaps, list): - aligned_sta1, aligned_sta2 = align_by_xcorr( - sta_emgfile1[tracking_res["MU_file1"].loc[ind]], - sta_emgfile2[tracking_res["MU_file2"].loc[ind]], - finalduration=0.5, - ) - else: - aligned_sta1 = sta_emgfile1[tracking_res["MU_file1"].loc[ind]] - aligned_sta2 = sta_emgfile2[tracking_res["MU_file2"].loc[ind]] - - title = "MUAPs from MU '{}' in file 1 and MU '{}' in file 2, XCC = {}".format( - tracking_res["MU_file1"].loc[ind], - tracking_res["MU_file2"].loc[ind], - round(tracking_res["XCC"].loc[ind], 2), - ) - plot_muaps( - sta_dict=[aligned_sta1, aligned_sta2], - title=title, - showimmediately=False + """ if show: + def parallel(ind): # Function for the parallel execution of plotting """ # TODO + """ for ind in tracking_res.index: + if tracking_res["XCC"].loc[ind] >= threshold: + # Align STA + if not isinstance(custom_muaps, list): + aligned_sta1, aligned_sta2 = align_by_xcorr( + sta_emgfile1[tracking_res["MU_file1"].loc[ind]], + sta_emgfile2[tracking_res["MU_file2"].loc[ind]], + finalduration=0.5, ) + else: + aligned_sta1 = sta_emgfile1[tracking_res["MU_file1"].loc[ind]] + aligned_sta2 = sta_emgfile2[tracking_res["MU_file2"].loc[ind]] - plt.show() + title = "MUAPs from MU '{}' in file 1 and MU '{}' in file 2, XCC = {}".format( + tracking_res["MU_file1"].loc[ind], + tracking_res["MU_file2"].loc[ind], + round(tracking_res["XCC"].loc[ind], 2), + ) + plot_muaps( + sta_dict=[aligned_sta1, aligned_sta2], + title=title, + showimmediately=False + ) - # Check that the number of plots does not exceed the number of cores + if ind == tracking_res.index[-1]: + plt.show(block=True) + else: + plt.show(block=False) """ + + """ # Check that the number of plots does not exceed the number of cores num_cores = os.cpu_count() if len(tracking_res.index) > num_cores: # If yes, raise a warning warnings.warn("\n\nThere are more plots to show than available cores\n") # Parallel execution of plotting - Parallel(n_jobs=-1)(delayed(parallel)(ind) for ind in tracking_res.index) + Parallel(n_jobs=-1)(delayed(parallel)(ind) for ind in tracking_res.index) """ return tracking_res +class Tracking_gui(): # TODO return updated results, do not run if no pairs are detected when called from tracking, add delete excluded pairs button + """ + GUI for the visual inspection and manual selection of the tracking results. + + Parameters + ---------- + emgfile1 : dict + The dictionary containing the first emgfile. + emgfile2 : dict + The dictionary containing the second emgfile. + tracking_res : pd.DataFrame + The results of the tracking including the MU from file 1, + MU from file 2 and the normalised cross-correlation value (XCC). + This is obtained with the function `tracking()`. + sta_emgfile1 : dict + dict containing a dict of STA (pd.DataFrame) for every MUs from + emgfile1. This is obtained with the function `sta()`. + sta_emgfile2 : dict + dict containing a dict of STA (pd.DataFrame) for every MUs from + emgfile2. This is obtained with the function `sta()`. + align_muaps : bool, default True + Whether to align the MUAPs before plotting. If true, the visualised + MUAPs time window will be 1/2 of the original (because the maximum + allowed shift during MUAPs alignment is 50%). + addrefsig : bool, default True + If True, the REF_SIGNAL is plotted in front of the IDR with a + separated y-axes. + csv_separator : str, default "\t" + The field delimiter used to create the .csv copied to the clipboard. + + Methods + ------- + get_results() + Returns the results of the tracking including the MU from file 1, + MU from file 2, the normalised cross-correlation value (XCC) and the + inclusion or exclusion of the MUs pair. after the GUI is closed. + + See also + -------- + - tracking : Track MUs across two files comparing the MUAPs' shape and + distribution. + + Examples + -------- + Track MUs between two files. + + >>> import openhdemg.library as emg + >>> emgfile_1 = emg.askopenfile(filesource="OPENHDEMG") + >>> emgfile_2 = emg.askopenfile(filesource="OPENHDEMG") + >>> tracking_res = emg.tracking(emgfile_1, emgfile_2) + + Obtained required variables for Tracking_gui(). Pay attention to use the + same derivation specified during tracking and to use an appropriate MUAPs + timewindow, according to the align_muaps option in Tracking_gui(). + + >>> sorted_rawemg_1 = emg.sort_rawemg(emgfile_1) + >>> sorted_rawemg_1 = emg.diff(sorted_rawemg_1) + >>> sta_dict_1 = emg.sta(emgfile_1, sorted_rawemg_1, timewindow=100) + >>> sorted_rawemg_2 = emg.sort_rawemg(emgfile_2) + >>> sorted_rawemg_2 = emg.diff(sorted_rawemg_2) + >>> sta_dict_2 = emg.sta(emgfile_2, sorted_rawemg_2, timewindow=100) + + Inspect the tracking results with a convenient GUI. + + >>> tracking = emg.Tracking_gui( + ... emgfile1=emgfile_1, + ... emgfile2=emgfile_2, + ... tracking_res=tracking_res, + ... sta_emgfile1=sta_dict_1, + ... sta_emgfile2=sta_dict_2, + ... align_muaps=True, + ... addrefsig=True, + ... ) + + Return the updated results for further use in the code. + + >>> updated_results = tracking.get_results() + MU_file1 MU_file2 XCC Inclusion + 0 0 2 0.897364 Excluded + 1 1 0 0.947486 Included + 2 2 1 0.923901 Included + 3 4 8 0.893922 Included + """ + def __init__( + self, + emgfile1, + emgfile2, + tracking_res, + sta_emgfile1, + sta_emgfile2, + align_muaps=True, + addrefsig=True, + csv_separator="\t", + ): + # Define needed variables + self.emgfile1 = emgfile1 + self.emgfile2 = emgfile2 + self.tracking_res = tracking_res + self.sta_emgfile1 = sta_emgfile1 + self.sta_emgfile2 = sta_emgfile2 + self.align_muaps = align_muaps + self.addrefsig = addrefsig + self.csv_separator = csv_separator + + # Add included/excluded label to self.tracking_res + self.tracking_res = self.tracking_res.assign(Inclusion="Included") + + # Define GUI structure + # After that, set up the GUI + self.root = tk.Tk() + self.root.title('MUAPs tracking') + root_path = os.path.dirname(os.path.abspath(__file__)) + iconpath = os.path.join( + root_path, + "..", + "gui", + "gui_files", + "Icon_transp.ico" + ) + self.root.iconbitmap(iconpath) + + # Create outer frames, assign structure and minimum spacing + # Top + top_frm = tk.Frame(self.root, padx=10) + top_frm.grid( + row=0, column=0, columnspan=3, sticky=tk.NSEW, pady=(10, 8), + ) + # Central - Left + self.central_left_frm = tk.Frame( + self.root, padx=1, pady=1, background="gray", + ) + self.central_left_frm.grid( + row=1, rowspan=2, column=0, columnspan=2, sticky=tk.NSEW, + ) + # Central - Right Top + self.central_right_top_frm = tk.Frame( + self.root, padx=0, pady=1, background="gray", + ) + self.central_right_top_frm.grid( + row=1, column=2, columnspan=1, sticky=tk.NSEW, + ) + # Central - Right Bottom + self.central_right_bottom_frm = tk.Frame( + self.root, padx=0, pady=1, background="gray", + ) + self.central_right_bottom_frm.grid( + row=2, column=2, columnspan=1, sticky=tk.NSEW, + ) + # Bottom + bottom_frm = tk.Frame(self.root, padx=1, pady=0, background="gray") + bottom_frm.grid(row=3, column=0, columnspan=3, sticky=tk.NSEW) + + # Assign rows and columns weight to avoid empty spaces + self.root.rowconfigure(1, weight=1) + self.root.rowconfigure(2, weight=1) + + self.root.columnconfigure(0, weight=1) + self.root.columnconfigure(1, weight=1) + self.root.columnconfigure(2, weight=1) + + # Label MU pair selection and create combobox to change MUs pair + # File 1 + mu_pair_label = ttk.Label(top_frm, text="Pair of MUs to visualise:") + mu_pair_label.pack(side=tk.LEFT) + + mu_pairs = list(self.tracking_res.index) + self.select_mu_pair_cb = ttk.Combobox( + top_frm, + values=mu_pairs, + state='readonly', + width=5, + ) + self.select_mu_pair_cb.pack(side=tk.LEFT, padx=(2, 15)) + self.select_mu_pair_cb.current(0) + # gui_plot() takes one positional argument (self), but the bind() + # method is passing two arguments: the event object and the function + # itself. Use lambda to avoid the error. + self.select_mu_pair_cb.bind( + '<>', + lambda event: self.gui_plot(), + ) + + # Button to copy the dataframe to clipboard + copy_btn = ttk.Button( + top_frm, + text="Copy results", + command=self.copy_to_clipboard, + ) + copy_btn.pack(side=tk.RIGHT, padx=(20, 0)) + + # Button to include/exclude MUAPs + btn_include_muaps = ttk.Button( + top_frm, + text="Include/Exclude", + command=self.include_exclude, + ) + btn_include_muaps.pack(side=tk.RIGHT, padx=(0, 20)) + + # Inclusion label + self.inclusion_label = ttk.Label( + top_frm, text="INCLUDED", foreground="green", + ) + self.inclusion_label.pack(side=tk.RIGHT, padx=(20, 20)) + + # Text frame to show the tracking results + self.textbox = tk.Text(bottom_frm, height=10, relief="flat") + self.textbox.pack( + expand=True, side=tk.BOTTOM, fill=tk.X) + self.textbox.insert( + '1.0', + self.tracking_res.to_string(float_format="{:.2f}".format), + ) + + # Highlight the specific row + self.highlight_row(0) + + # Lower the window to prevent flashing + self.root.lower() + + # Plot the first MU pair + self.gui_plot() + + # Bring back the GUI in in the foreground + self.root.lift() + + # Start mainloop + self.root.mainloop() + + def highlight_row(self, row_number): + # Highlight a row in self.textbox + + # Define a tag for highlighting + self.textbox.tag_configure("highlight", background="yellow") + + # Update rownumber for the table format + row_number += 2 + + # Calculate the start and end index of the row in the text widget + start_index = f"{row_number}.0" + end_index = f"{row_number}.end" + + # Add the highlight tag to the row + self.textbox.tag_add("highlight", start_index, end_index) + + def gui_plot(self): + # Display the MUAPs and IDR. + + # Update the textbox highlight + pair = int(self.select_mu_pair_cb.get()) + self.textbox.tag_delete("highlight") + self.highlight_row(pair) + + # Update the inclusion_label + if self.tracking_res.loc[pair, "Inclusion"] == "Excluded": + self.inclusion_label.config( + text="EXCLUDED", + foreground="red" + ) + else: + self.inclusion_label.config( + text="INCLUDED", + foreground="green" + ) + + # Get MUs number + mu1 = int(self.tracking_res["MU_file1"].loc[pair]) + mu2 = int(self.tracking_res["MU_file2"].loc[pair]) + + # MUAPs figure + if self.align_muaps: + aligned_sta1, aligned_sta2 = align_by_xcorr( + self.sta_emgfile1[mu1], + self.sta_emgfile2[mu2], + finalduration=0.5, + ) + muaps_fig = plot_muaps( + sta_dict=[aligned_sta1, aligned_sta2], + title="", + figsize=[5, 5], + showimmediately=False, + tight_layout=False, + ) + else: + muaps_fig = plot_muaps( + sta_dict=[self.sta_emgfile1[mu1], self.sta_emgfile2[mu2]], + title="", + figsize=[5, 5], + showimmediately=False, + tight_layout=False, + ) + + # If canvas already exist, destroy it + if hasattr(self, 'muaps_canvas'): + self.muaps_canvas.get_tk_widget().destroy() + + # Place the MUAPs figure in the canvas + self.muaps_canvas = FigureCanvasTkAgg( + muaps_fig, master=self.central_left_frm, + ) + self.muaps_canvas.draw_idle() # Await resizing + self.muaps_canvas.get_tk_widget().pack( + expand=True, fill="both", padx=0, pady=0, + ) + plt.close(muaps_fig) + + # Display IDR mu1 figure + idr_mu1_fig = plot_idr( + emgfile=self.emgfile1, + munumber=mu1, + addrefsig=self.addrefsig, + figsize=[3, 3], + showimmediately=False, + tight_layout=False, + ) + + if hasattr(self, 'idr_mu1_canvas'): + self.idr_mu1_canvas.get_tk_widget().destroy() + + self.idr_mu1_canvas = FigureCanvasTkAgg( + idr_mu1_fig, master=self.central_right_top_frm, + ) + self.idr_mu1_canvas.draw_idle() # Await resizing + self.idr_mu1_canvas.get_tk_widget().pack( + expand=True, fill="both", padx=0, pady=0, + ) + plt.close(idr_mu1_fig) + + # Display IDR mu2 figure + idr_mu2_fig = plot_idr( + emgfile=self.emgfile2, + munumber=mu2, + addrefsig=self.addrefsig, + figsize=[3, 3], + showimmediately=False, + tight_layout=False, + ) # TODO plot orange as MUAPs + + if hasattr(self, 'idr_mu2_canvas'): + self.idr_mu2_canvas.get_tk_widget().destroy() + + self.idr_mu2_canvas = FigureCanvasTkAgg( + idr_mu2_fig, master=self.central_right_bottom_frm, + ) + self.idr_mu2_canvas.draw_idle() # Await resizing + self.idr_mu2_canvas.get_tk_widget().pack( + expand=True, fill="both", padx=0, pady=0, + ) + plt.close(idr_mu2_fig) + + def include_exclude(self): + # Include or exclude the current MU pair + + # Check if the current pair is included or excluded and update label + pair = int(self.select_mu_pair_cb.get()) + if self.tracking_res.loc[pair, "Inclusion"] == "Included": + self.tracking_res.loc[pair, "Inclusion"] = "Excluded" + self.inclusion_label.config( + text="EXCLUDED", + foreground="red" + ) + else: + self.tracking_res.loc[pair, "Inclusion"] = "Included" + self.inclusion_label.config( + text="INCLUDED", + foreground="green" + ) + + # Update table + self.textbox.replace( + '1.0', + 'end', + self.tracking_res.to_string(float_format="{:.2f}".format), + ) + self.textbox.tag_delete("highlight") + self.highlight_row(pair) + + def copy_to_clipboard(self): + # Copy the dataframe to clipboard in csv format. + + self.tracking_res.to_clipboard(excel=True, sep=self.csv_separator) + + def get_results(self): + # Get the edited tracking_res + + return self.tracking_res + + def remove_duplicates_between( emgfile1, emgfile2, @@ -1561,6 +1954,7 @@ def estimate_cv_via_mle(emgfile, signal): return cv +# TODO add function to return the results class MUcv_gui(): """ Graphical user interface for the estimation of MUs conduction velocity. From 3d4a1e1e877d999d4da40509e3f9bc9db1c6c157 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Wed, 29 May 2024 11:43:05 +0200 Subject: [PATCH 05/35] Fit the necessities of Tracking_GUI() --- openhdemg/library/plotemg.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/openhdemg/library/plotemg.py b/openhdemg/library/plotemg.py index b8d1219..ac8a357 100644 --- a/openhdemg/library/plotemg.py +++ b/openhdemg/library/plotemg.py @@ -148,6 +148,7 @@ def plot_emgsig( for count, thisChannel in enumerate(channels): # Normalise the series norm_raw = min_max_scaling(emgfile["RAW_SIGNAL"][thisChannel]) + # TODO option to scale all the df or the single channel # Add 1 to the previous channel to avoid overlapping norm_raw = norm_raw + count plt.plot(x_axis, norm_raw) @@ -851,8 +852,6 @@ def plot_idr( ) # Check if we have a single MU or a list of MUs to plot. - # In the first case, the use of plt.figure has been preferred to - # plt.subplots for the implementation of the MUs cleaning. if isinstance(munumber, int): ax1 = sns.scatterplot( x=idr[munumber]["timesec" if timeinseconds else "mupulses"], @@ -923,6 +922,7 @@ def plot_muaps( title="MUAPs from STA", figsize=[20, 15], showimmediately=True, + tight_layout=False, ): """ Plot MUAPs obtained from STA from one or multiple MUs. @@ -942,6 +942,11 @@ def plot_muaps( If True (default), plt.show() is called and the figure showed to the user. It is useful to set it to False when calling the function from the GUI. + tight_layout : bool, default False + If True (default), the plt.tight_layout() is called and the figure's + layout is improved. + It is useful to set it to False when calling the function from the GUI + or if the final layout is not correct. Returns ------- @@ -1088,7 +1093,7 @@ def plot_muaps( "Unacceptable number of rows and columns to plot" ) - showgoodlayout(tight_layout=False, despined=True) + showgoodlayout(tight_layout=tight_layout, despined=True) if showimmediately: plt.show() @@ -1328,6 +1333,7 @@ def plot_muaps_for_cv( title="MUAPs for CV", figsize=[20, 15], showimmediately=True, + tight_layout=False, ): """ Visualise MUAPs on which to calculate MUs CV. @@ -1351,6 +1357,11 @@ def plot_muaps_for_cv( If True (default), plt.show() is called and the figure showed to the user. It is useful to set it to False when calling the function from the GUI. + tight_layout : bool, default False + If True (default), the plt.tight_layout() is called and the figure's + layout is improved. + It is useful to set it to False when calling the function from the GUI + or if the final layout is not correct. Returns ------- From 06ad933c66186a33031fb237db06f585846321d7 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Wed, 29 May 2024 16:05:34 +0200 Subject: [PATCH 06/35] Update tracking() Possibility to directly call Tracking_gui() and other improvements --- openhdemg/library/muap.py | 212 +++++++++++++++++++++++++------------- 1 file changed, 141 insertions(+), 71 deletions(-) diff --git a/openhdemg/library/muap.py b/openhdemg/library/muap.py index 2323fb3..0295350 100644 --- a/openhdemg/library/muap.py +++ b/openhdemg/library/muap.py @@ -8,7 +8,7 @@ from openhdemg.library.mathtools import ( norm_twod_xcorr, norm_xcorr, - find_teta, + find_mle_teta, mle_cv_est, ) from openhdemg.library.electrodes import sort_rawemg @@ -696,7 +696,6 @@ def align_by_xcorr(sta_mu1, sta_mu2, finalduration=0.5): ) normxcorr_df = normxcorr_df.set_index(corr_lags) lag = normxcorr_df.idxmax().median() # First signal compared to second - # TODO Can mean provide better results? allow choice? # Be sure that the lag/delay does not exceed values suitable for the final # expected duration. @@ -736,12 +735,7 @@ def align_by_xcorr(sta_mu1, sta_mu2, finalduration=0.5): # TODO update examples for code="None" -# TODO change parallel processing, allow to specify if parallel -# TODO add GUI - -# This function exploits parallel processing: -# - align and xcorr are processed in parallel -# - plotting is processed in parallel +# This function exploits parallel processing for MUAPs alignment and xcorr def tracking( emgfile1, emgfile2, @@ -757,11 +751,18 @@ def tracking( custom_muaps=None, exclude_belowthreshold=True, filter=True, + multiprocessing=True, show=False, -): # TODO add gui example + gui=True, + gui_addrefsig=True, + gui_csv_separator="\t", +): """ Track MUs across two files comparing the MUAPs' shape and distribution. + It is also possible to use a convenient GUI for the inspection of the + obtained MU pairs. + Parameters ---------- emgfile1 : dict @@ -842,14 +843,31 @@ def tracking( filter : bool, default True If true, when the same MU has a match of XCC > threshold with multiple MUs, only the match with the highest XCC is returned. + multiprocessing : bool, default True + If True (default) parallel processing will be used to reduce execution + time. show : bool, default False - Whether to plot the STA of pairs of MUs with XCC above threshold. + Whether to plot the STA of pairs of MUs with XCC above threshold. Set + to False (default) when gui=True to avoid postponing the GUI execution. + If show=True and gui=True, the GUI will be executed after closing all + the figures. + gui : bool, default True + If True (default) a GUI for the visual inspection and manual selection + of the tracking results will be called. + gui_addrefsig : bool, default True + If True, the REF_SIGNAL is plotted in front of the IDR with a + separated y-axes. This is used only when gui=True. + gui_csv_separator : str, default "\t" + The field delimiter used by the GUI to create the .csv copied to the + clipboard. This is used only when gui=True. Returns ------- tracking_res : pd.DataFrame The results of the tracking including the MU from file 1, MU from file 2 and the normalised cross-correlation value (XCC). + If gui=True, an additional column indicating the inclusion/exclusion + of the MUs pairs will also be present in tracking_res. Warns ----- @@ -874,6 +892,26 @@ def tracking( Examples -------- + Track MUs between two OPENHDEMG (.json) files and inspect the results with + a convenient GUI. + + >>> import openhdemg.library as emg + >>> emgfile1 = emg.askopenfile(filesource="OPENHDEMG") + >>> emgfile2 = emg.askopenfile(filesource="OPENHDEMG") + >>> tracking_res = emg.tracking( + ... emgfile1=emgfile1, + ... emgfile2=emgfile2, + ... firings="all", + ... derivation="sd", + ... timewindow=50, + ... threshold=0.8, + ... matrixcode="GR08MM1305", + ... orientation=180, + ... filter=True, + ... show=False, + ... gui=True, + >>> ) + Track MUs between two OTB files and show the filtered results. >>> import openhdemg.library as emg @@ -893,6 +931,7 @@ def tracking( ... exclude_belowthreshold=True, ... filter=True, ... show=False, + ... gui=False, ... ) MU_file1 MU_file2 XCC 0 0 3 0.820068 @@ -934,6 +973,7 @@ def tracking( ... exclude_belowthreshold=True, ... filter=True, ... show=False, + ... gui=False, ... ) """ @@ -997,7 +1037,7 @@ def tracking( else: raise ValueError("custom_muaps is not a list of two dictionaries") - print("\nTracking started") + print("\nTracking started:") # Tracking function to run in parallel def parallel(mu_file1): # Loop all the MUs of file 1 @@ -1039,16 +1079,27 @@ def parallel(mu_file1): # Loop all the MUs of file 1 return res - # Start parallel execution - # Measure running time - t0 = time.time() - - res = Parallel(n_jobs=-1)( - delayed(parallel)(mu_file1) for mu_file1 in range(emgfile1["NUMBER_OF_MUS"]) - ) # TODO n_jobs -1 remove joblib and allow user to select number of cores + if multiprocessing: + # Start parallel execution + res = Parallel(n_jobs=-1, verbose=1)( + delayed(parallel)(mu_file1) for mu_file1 in range(emgfile1["NUMBER_OF_MUS"]) + ) + print("\n") - t1 = time.time() - print(f"\nTime of tracking parallel processing: {round(t1-t0, 2)} Sec\n") + else: + # Start serial execution + t0 = time.time() + + res = [] + for pos, mu_file1 in enumerate(range(emgfile1["NUMBER_OF_MUS"])): + res.append(parallel(mu_file1)) + # Show progress + t1 = time.time() + print( + f"Done {pos+1} out of {emgfile1['NUMBER_OF_MUS']} | " + + f"Elapsed time: {round(t1-t0, 2)} Sec", + ) + print("\n") # Convert res to pd.DataFrame for pos, i in enumerate(res): @@ -1061,32 +1112,39 @@ def parallel(mu_file1): # Loop all the MUs of file 1 # Filter the results if filter: # Sort file by MUs in file 1 and XCC to have first the highest XCC - sorted_res = tracking_res.sort_values( - by=["MU_file1", "XCC"], ascending=False + tracking_res = tracking_res.sort_values( + by=["MU_file1", "XCC"], ascending=False, ) + # Get unique MUs from file 1 - unique = sorted_res["MU_file1"].unique() + unique = tracking_res["MU_file1"].unique() - res_unique = pd.DataFrame(columns=sorted_res.columns) + res_unique = pd.DataFrame(columns=tracking_res.columns) # Get the combo uf unique MUs from file 1 with MUs from file 2 for pos, mu1 in enumerate(unique): - this_res = sorted_res[sorted_res["MU_file1"] == mu1] + this_res = tracking_res[tracking_res["MU_file1"] == mu1] # Fill the result df with the first row (highest XCC) res_unique.loc[pos, :] = this_res.iloc[0, :] # Now repeat the task with MUs from file 2 - sorted_res = res_unique.sort_values( + tracking_res = res_unique.sort_values( by=["MU_file2", "XCC"], ascending=False ) - unique = sorted_res["MU_file2"].unique() - res_unique = pd.DataFrame(columns=sorted_res.columns) + unique = tracking_res["MU_file2"].unique() + res_unique = pd.DataFrame(columns=tracking_res.columns) for pos, mu2 in enumerate(unique): - this_res = sorted_res[sorted_res["MU_file2"] == mu2] + this_res = tracking_res[tracking_res["MU_file2"] == mu2] res_unique.loc[pos, :] = this_res.iloc[0, :] tracking_res = res_unique.sort_values(by=["MU_file1"]) + else: + # Sort file by MUs in file 1 and XCC to have first the highest XCC + tracking_res = tracking_res.sort_values( + by=["MU_file1", "XCC"], ascending=[True, False], + ) + # Print the full results pd.set_option("display.max_rows", None) convert_dict = {"MU_file1": int, "MU_file2": int, "XCC": float} @@ -1097,50 +1155,61 @@ def parallel(mu_file1): # Loop all the MUs of file 1 pd.reset_option("display.max_rows") # Plot the MUs pairs - """ if show: - def parallel(ind): # Function for the parallel execution of plotting """ # TODO - """ for ind in tracking_res.index: - if tracking_res["XCC"].loc[ind] >= threshold: - # Align STA - if not isinstance(custom_muaps, list): - aligned_sta1, aligned_sta2 = align_by_xcorr( - sta_emgfile1[tracking_res["MU_file1"].loc[ind]], - sta_emgfile2[tracking_res["MU_file2"].loc[ind]], - finalduration=0.5, - ) - else: - aligned_sta1 = sta_emgfile1[tracking_res["MU_file1"].loc[ind]] - aligned_sta2 = sta_emgfile2[tracking_res["MU_file2"].loc[ind]] - - title = "MUAPs from MU '{}' in file 1 and MU '{}' in file 2, XCC = {}".format( - tracking_res["MU_file1"].loc[ind], - tracking_res["MU_file2"].loc[ind], - round(tracking_res["XCC"].loc[ind], 2), - ) - plot_muaps( - sta_dict=[aligned_sta1, aligned_sta2], - title=title, - showimmediately=False - ) + if show: + for ind in tracking_res.index: + if tracking_res["XCC"].loc[ind] >= threshold: + # Align STA + if not isinstance(custom_muaps, list): + aligned_sta1, aligned_sta2 = align_by_xcorr( + sta_emgfile1[tracking_res["MU_file1"].loc[ind]], + sta_emgfile2[tracking_res["MU_file2"].loc[ind]], + finalduration=0.5, + ) + else: + aligned_sta1 = sta_emgfile1[tracking_res["MU_file1"].loc[ind]] + aligned_sta2 = sta_emgfile2[tracking_res["MU_file2"].loc[ind]] - if ind == tracking_res.index[-1]: - plt.show(block=True) - else: - plt.show(block=False) """ + title = "MUAPs from MU '{}' in file 1 and MU '{}' in file 2, XCC = {}".format( + tracking_res["MU_file1"].loc[ind], + tracking_res["MU_file2"].loc[ind], + round(tracking_res["XCC"].loc[ind], 2), + ) + plot_muaps( + sta_dict=[aligned_sta1, aligned_sta2], + title=title, + showimmediately=False + ) - """ # Check that the number of plots does not exceed the number of cores - num_cores = os.cpu_count() - if len(tracking_res.index) > num_cores: - # If yes, raise a warning - warnings.warn("\n\nThere are more plots to show than available cores\n") + if ind == tracking_res.index[-1]: + plt.show(block=True) + else: + plt.show(block=False) + + # Call the GUI and return the tracking_res + if gui: + # Check if no pairs have been detected + if tracking_res.empty: + return tracking_res + + tracking_gui = Tracking_gui( + emgfile1=emgfile1, + emgfile2=emgfile2, + tracking_res=tracking_res, + sta_emgfile1=sta_emgfile1, + sta_emgfile2=sta_emgfile2, + align_muaps=False if isinstance(custom_muaps, list) else True, + addrefsig=gui_addrefsig, + csv_separator=gui_csv_separator, + ) - # Parallel execution of plotting - Parallel(n_jobs=-1)(delayed(parallel)(ind) for ind in tracking_res.index) """ + # Get and return updated tracking_res + return tracking_gui.get_results() - return tracking_res + else: + return tracking_res -class Tracking_gui(): # TODO return updated results, do not run if no pairs are detected when called from tracking, add delete excluded pairs button +class Tracking_gui(): # TODO add delete excluded pairs button """ GUI for the visual inspection and manual selection of the tracking results. @@ -1189,7 +1258,7 @@ class Tracking_gui(): # TODO return updated results, do not run if no pairs are >>> import openhdemg.library as emg >>> emgfile_1 = emg.askopenfile(filesource="OPENHDEMG") >>> emgfile_2 = emg.askopenfile(filesource="OPENHDEMG") - >>> tracking_res = emg.tracking(emgfile_1, emgfile_2) + >>> tracking_res = emg.tracking(emgfile_1, emgfile_2, timewindow=50) Obtained required variables for Tracking_gui(). Pay attention to use the same derivation specified during tracking and to use an appropriate MUAPs @@ -1527,6 +1596,7 @@ def get_results(self): return self.tracking_res +# TODO add GUI? def remove_duplicates_between( emgfile1, emgfile2, @@ -1927,7 +1997,7 @@ def estimate_cv_via_mle(emgfile, signal): sig = signal.to_numpy() sig = sig.T - # Prepare the input 1D signals for find_teta + # Prepare the input 1D signals for find_mle_teta if np.shape(sig)[0] > 3: sig1 = sig[1, :] sig2 = sig[2, :] @@ -1935,7 +2005,7 @@ def estimate_cv_via_mle(emgfile, signal): sig1 = sig[0, :] sig2 = sig[1, :] - teta = find_teta( + teta = find_mle_teta( sig1=sig1, sig2=sig2, ied=ied, @@ -2035,7 +2105,7 @@ def __init__( # After that, set up the GUI self.root = tk.Tk() - self.root.title('MUs cv estimation') + self.root.title('MUs CV estimation') root_path = os.path.dirname(os.path.abspath(__file__)) iconpath = os.path.join( root_path, From df24cda8f545884a66015a79ecbf97ef1c8efe66 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Wed, 29 May 2024 16:06:14 +0200 Subject: [PATCH 07/35] Replace find_teta with find_mle_teta --- openhdemg/library/mathtools.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/openhdemg/library/mathtools.py b/openhdemg/library/mathtools.py index f08b7aa..de6c0cb 100644 --- a/openhdemg/library/mathtools.py +++ b/openhdemg/library/mathtools.py @@ -119,19 +119,19 @@ def norm_xcorr(sig1, sig2, out="both"): def norm_twod_xcorr(df1, df2, mode="full"): """ - Normalised 2-dimensional cross-correlation of STAs of two MUS. + Normalised 2-dimensional cross-correlation of 2. - Any pre-processing of the RAW_SIGNAL (i.e., normal, differential or double - differential) can be passed as long as the two inputs have same shape. + The two inputs must have same shape. + When this function is used to cross-correlate MUAPs obtained via STA, + df1 and df2 should contain the unpacked STA of the first and second MU, + respectively, without np.nan columns. Parameters ---------- df1 : pd.DataFrame - A pd.DataFrame containing the STA of the first MU without np.nan - column. + A pd.DataFrame containing the first 2-dimensional signal. df2 : pd.DataFrame - A pd.DataFrame containing the STA of the second MU without np.nan - column. + A pd.DataFrame containing the second 2-dimensional signal. mode : str {"full", "valid", "same"}, default "full" A string indicating the size of the output: @@ -692,7 +692,7 @@ def mle_cv_est(sig, initial_teta, ied, fsamp): Examples -------- - Refer to the examples of find_teta to obtain sig and initial_teta. + Refer to the examples of find_mle_teta to obtain sig and initial_teta. """ # Calculate ied in meters @@ -737,7 +737,7 @@ def mle_cv_est(sig, initial_teta, ied, fsamp): return cv, teta -def find_teta(sig1, sig2, ied, fsamp): +def find_mle_teta(sig1, sig2, ied, fsamp): """ Find the starting value for teta. @@ -811,7 +811,7 @@ def find_teta(sig1, sig2, ied, fsamp): Third, estimate teta. - >>> initial_teta = emg.find_teta( + >>> initial_teta = emg.find_mle_teta( ... sig1=sig1, ... sig2=sig2, ... ied=emgfile["IED"], From fc85a9ba2296657108225a37d9bcdf99ac50e6c9 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Thu, 30 May 2024 10:08:29 +0200 Subject: [PATCH 08/35] Add tracking gui in main GUI --- .../gui/gui_modules/advanced_analyses.py | 51 ++++++++++++------- openhdemg/gui/settings.py | 4 ++ 2 files changed, 37 insertions(+), 18 deletions(-) diff --git a/openhdemg/gui/gui_modules/advanced_analyses.py b/openhdemg/gui/gui_modules/advanced_analyses.py index a4e95cc..82ed6fb 100644 --- a/openhdemg/gui/gui_modules/advanced_analyses.py +++ b/openhdemg/gui/gui_modules/advanced_analyses.py @@ -343,7 +343,7 @@ def advanced_analysis(self): self.head.rowconfigure(row, weight=1) # Specify Signal - signal_value = ("OTB", "DEMUSE", "OPENHDEMG", "CUSTOMCSV") + signal_value = ("OPENHDEMG", "DEMUSE", "OTB", "CUSTOMCSV") signal_entry = ctk.CTkComboBox( self.head, width=150, @@ -471,6 +471,15 @@ def advanced_analysis(self): for child in self.head.winfo_children(): child.grid_configure(padx=5, pady=5) + # Add description for the show_checkbox for Motor Unit Tracking + if self.advanced_method.get() == "Motor Unit Tracking": + description_label = ctk.CTkLabel( + self.head, + text="Leave Show off if using \n the tracking GUI", + font=("Segoe UI", 10.5), + ) + description_label.grid(column=0, row=14, pady=(0, 4)) + # Add Which widget and update the track button # to match functionalities required for duplicate removal if self.advanced_method.get() == "Duplicate Removal": @@ -561,6 +570,7 @@ def advanced_analysis(self): n_firings=self.parent.settings.MUcv_gui__n_firings, muaps_timewindow=self.parent.settings.MUcv_gui__muaps_timewindow, figsize=self.parent.settings.MUcv_gui__figsize, + csv_separator=self.parent.settings.MUcv_gui__csv_separator, ) except AttributeError as e: @@ -804,26 +814,31 @@ def track_mus(self): exclude_belowthreshold=self.exclude_thres.get(), filter=self.filter_adv.get(), show=self.show_adv.get(), + gui=self.parent.settings.tracking__gui, + gui_addrefsig=self.parent.settings.tracking__gui_addrefsig, + gui_csv_separator=self.parent.settings.tracking__gui_csv_separator, ) - # Add result terminal - track_terminal = ttk.LabelFrame( - self.head, text="MUs Tracking Result", height=100, - relief="ridge", - ) - track_terminal.grid( - column=2, - row=0, - columnspan=2, - rowspan=14, - pady=8, - padx=10, - sticky=(N, S, W, E), - ) + # Add result terminal, only if tracking__gui=False to avoid + # uninterrupted processes. + if not self.parent.settings.tracking__gui: + track_terminal = ttk.LabelFrame( + self.head, text="MUs Tracking Result", height=100, + relief="ridge", + ) + track_terminal.grid( + column=2, + row=0, + columnspan=2, + rowspan=14, + pady=8, + padx=10, + sticky=(N, S, W, E), + ) - # Add table containing results to the label frame - track_table = Table(track_terminal, dataframe=tracking_res) - track_table.show() + # Add table containing results to the label frame + track_table = Table(track_terminal, dataframe=tracking_res) + track_table.show() except AttributeError as e: show_error_dialog( diff --git a/openhdemg/gui/settings.py b/openhdemg/gui/settings.py index fb335ea..b14f761 100644 --- a/openhdemg/gui/settings.py +++ b/openhdemg/gui/settings.py @@ -86,6 +86,9 @@ # in tracking() tracking__firings = "all" tracking__derivation = "sd" +tracking__gui = True +tracking__gui_addrefsig = True +tracking__gui_csv_separator = "\t" # in remove_duplicates_between() remove_duplicates_between__firings = "all" @@ -95,6 +98,7 @@ MUcv_gui__n_firings = [0, 50] MUcv_gui__muaps_timewindow = 50 MUcv_gui__figsize = [25, 20] +MUcv_gui__csv_separator = "\t" # --------------------------------- electrodes -------------------------------- # This custom sorting order is valid for all the GUI windows, although the From 0ac6c3db843d316ace984641a62aacfb1a079b2c Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Thu, 30 May 2024 17:36:54 +0200 Subject: [PATCH 09/35] Create test_anlysis.py --- openhdemg/tests/unit/test_anlysis.py | 398 +++++++++++++++++++++++++++ 1 file changed, 398 insertions(+) create mode 100644 openhdemg/tests/unit/test_anlysis.py diff --git a/openhdemg/tests/unit/test_anlysis.py b/openhdemg/tests/unit/test_anlysis.py new file mode 100644 index 0000000..38f9abd --- /dev/null +++ b/openhdemg/tests/unit/test_anlysis.py @@ -0,0 +1,398 @@ +""" +To run the tests using unittest, execute from the openhdemg/tests directory: + python -m unittest discover + +First, you should dowload all the files necessary for the testing and store +them inside openhdemg/tests/fixtures. The files are available at: +https://drive.google.com/drive/folders/1suCZSils8rSCs2E3_K25vRCbN3AFDI7F?usp=sharing + +IMPORTANT: Do not alter the content of the dowloaded folder! + +WARNING!!! Since the library's functions perform complex tasks and return +complex data structures, these tests can verify that no critical errors occur, +but the accuracy of each function must be assessed independently upon creation, +or at each revision of the code. + +WARNING!!! - UNTESTED FUNCTIONS: none +""" + + +import unittest +from openhdemg.tests.unit.functions_for_unit_test import get_directories as getd + +from openhdemg.library.openfiles import emg_from_samplefile + +from openhdemg.library.analysis import ( + compute_thresholds, compute_dr, basic_mus_properties, compute_covisi, + compute_drvariability, +) + +import numpy as np + + +class TestAnalysis(unittest.TestCase): + """ + Test the functions/classes in the analysis module. + """ + + def test_compute_thresholds(self): + """ + Test the compute_thresholds function with the samplefile. + """ + + # Load the decomposed samplefile + emgfile = emg_from_samplefile() + + # Default parameters + res = compute_thresholds( + emgfile=emgfile, + event_="rt_dert", + type_="abs_rel", + n_firings=1, + mvc=1234, + ) + + self.assertAlmostEqual( + res["abs_RT"][0], 86.824, places=2, + ) + self.assertAlmostEqual( + res["abs_DERT"][1], 220.965, places=2, + ) + self.assertAlmostEqual( + res["rel_RT"][2], 12.491, places=2, + ) + self.assertAlmostEqual( + res["rel_DERT"][3], 7.373, places=2, + ) + + # Change n_firings + res = compute_thresholds( + emgfile=emgfile, + event_="rt_dert", + type_="abs_rel", + n_firings=5, + mvc=1234, + ) + + self.assertAlmostEqual( + res["abs_RT"][0], 170.833, places=2, + ) + self.assertAlmostEqual( + res["abs_DERT"][1], 244.561, places=2, + ) + self.assertAlmostEqual( + res["rel_RT"][2], 14.677, places=2, + ) + self.assertAlmostEqual( + res["rel_DERT"][3], 9.313, places=2, + ) + + # Change event_ + res = compute_thresholds( + emgfile=emgfile, + event_="rt", + type_="rel", + n_firings=5, + mvc=1234, + ) + + self.assertAlmostEqual( + res["rel_RT"][0], 13.843, places=2, + ) + + def test_compute_dr(self): + """ + Test the compute_dr function with the samplefile. + """ + + # Load the decomposed samplefile + emgfile = emg_from_samplefile() + + # Ramps duration + t_ramps = 10 * emgfile["FSAMP"] + + # Default parameters + res = compute_dr( + emgfile=emgfile, + n_firings_RecDerec=4, + n_firings_steady=10, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + ) + + self.assertAlmostEqual( + res["DR_rec"][0], 3.341, places=2, + ) + self.assertAlmostEqual( + res["DR_derec"][1], 4.662, places=2, + ) + self.assertAlmostEqual( + res["DR_start_steady"][2], 8.793, places=2, + ) + self.assertAlmostEqual( + res["DR_end_steady"][3], 10.718, places=2, + ) + self.assertAlmostEqual( + res["DR_all_steady"][4], 10.693, places=2, + ) + self.assertAlmostEqual( + res["DR_all"][3], 10.693, places=2, + ) + + # Change n_firings_RecDerec + res = compute_dr( + emgfile=emgfile, + n_firings_RecDerec=10, + n_firings_steady=10, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + ) + + self.assertAlmostEqual( + res["DR_rec"][0], 7.031, places=2, + ) + self.assertAlmostEqual( + res["DR_derec"][1], 5.596, places=2, + ) + + res = compute_dr( + emgfile=emgfile, + n_firings_RecDerec=1, + n_firings_steady=10, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + ) + + self.assertTrue(np.isnan(res["DR_rec"][0])) + self.assertTrue(np.isnan(res["DR_derec"][0])) + + # Change n_firings_steady + res = compute_dr( + emgfile=emgfile, + n_firings_RecDerec=4, + n_firings_steady=36, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + ) + + self.assertAlmostEqual( + res["DR_start_steady"][0], 9.971, places=2, + ) + self.assertAlmostEqual( + res["DR_end_steady"][1], 6.806, places=2, + ) + + res = compute_dr( + emgfile=emgfile, + n_firings_RecDerec=1, + n_firings_steady=1, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + ) + + self.assertAlmostEqual( + res["DR_start_steady"][0], 7.474, places=2, + ) + self.assertAlmostEqual( + res["DR_end_steady"][1], 7.062, places=2, + ) + + def test_basic_mus_properties(self): + """ + Test the basic_mus_properties function with the samplefile. + """ + + # Load the decomposed samplefile + emgfile = emg_from_samplefile() + + # Ramps duration + t_ramps = 10 * emgfile["FSAMP"] + + # Default parameters + res = basic_mus_properties( + emgfile=emgfile, + n_firings_rt_dert=1, + n_firings_RecDerec=4, + n_firings_steady=10, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + accuracy="default", + ignore_negative_ipts=False, + constrain_pulses=[True, 3], + mvc=1234, + ) + + self.assertAlmostEqual( + res["MVC"][0], 1234.0, places=0, + ) + self.assertAlmostEqual( + res["Accuracy"][1], 0.955, places=2, + ) + self.assertAlmostEqual( + res["avg_Accuracy"][0], 0.914, places=2, + ) + self.assertAlmostEqual( + res["COV_steady"][0], 1.316, places=2, + ) + + # Change accuracy estimation + res = basic_mus_properties( + emgfile=emgfile, + n_firings_rt_dert=1, + n_firings_RecDerec=4, + n_firings_steady=10, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + accuracy="SIL_PNR", + ignore_negative_ipts=True, + constrain_pulses=[True, 3], + mvc=1234, + ) + + self.assertAlmostEqual( + res["SIL"][1], 0.830, places=2, + ) + self.assertAlmostEqual( + res["avg_SIL"][0], 0.676, places=2, + ) + self.assertAlmostEqual( + res["PNR"][4], 28.469, places=2, + ) + self.assertAlmostEqual( + res["avg_PNR"][0], 29.113, places=2, + ) + + res = basic_mus_properties( + emgfile=emgfile, + n_firings_rt_dert=1, + n_firings_RecDerec=4, + n_firings_steady=10, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + accuracy="PNR", + ignore_negative_ipts=True, + constrain_pulses=[False, 3], + mvc=1234, + ) + + self.assertAlmostEqual( + res["PNR"][4], 28.079, places=2, + ) + self.assertAlmostEqual( + res["avg_PNR"][0], 29.895, places=2, + ) + + def test_compute_covisi(self): + """ + Test the compute_dr function with the samplefile. + """ + + # Load the decomposed samplefile + emgfile = emg_from_samplefile() + + # Ramps duration + t_ramps = 10 * emgfile["FSAMP"] + + # Default parameters + res = compute_covisi( + emgfile=emgfile, + n_firings_RecDerec=4, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + single_mu_number=-1, + ) + + self.assertAlmostEqual( + res["COVisi_rec"][0], 67.000, places=2, + ) + self.assertAlmostEqual( + res["COVisi_derec"][1], 24.007, places=2, + ) + self.assertAlmostEqual( + res["COVisi_steady"][2], 8.700, places=2, + ) + self.assertAlmostEqual( + res["COVisi_all"][3], 19.104, places=2, + ) + + # Change n_firings_RecDerec and event_ + res = compute_covisi( + emgfile=emgfile, + n_firings_RecDerec=1, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec", + single_mu_number=-1, + ) + + self.assertTrue(np.isnan(res["COVisi_rec"][0])) + self.assertTrue(np.isnan(res["COVisi_derec"][0])) + + # Change single_mu_number + res = compute_covisi( + emgfile=emgfile, + n_firings_RecDerec=4, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + single_mu_number=3, + ) + + self.assertAlmostEqual( + res["COVisi_all"][0], 19.104, places=2, + ) + + def test_compute_drvariability(self): + """ + Test the compute_drvariability function with the samplefile. + """ + + # Load the decomposed samplefile + emgfile = emg_from_samplefile() + + # Ramps duration + t_ramps = 10 * emgfile["FSAMP"] + + # Default parameters + res = compute_drvariability( + emgfile=emgfile, + n_firings_RecDerec=4, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + ) + + self.assertAlmostEqual( + res["DRvar_rec"][0], 109.254, places=2, + ) + self.assertAlmostEqual( + res["DRvar_derec"][1], 21.662, places=2, + ) + self.assertAlmostEqual( + res["DRvar_steady"][2], 8.840, places=2, + ) + self.assertAlmostEqual( + res["DRvar_all"][3], 12.803, places=2, + ) + + # Change n_firings_RecDerec + res = compute_drvariability( + emgfile=emgfile, + n_firings_RecDerec=1, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + ) + + self.assertTrue(np.isnan(res["DRvar_rec"][0])) + self.assertTrue(np.isnan(res["DRvar_derec"][1])) + + +if __name__ == '__main__': + unittest.main() From d2236b66e926ba2aac288236e339ca7c11fedf60 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Mon, 3 Jun 2024 11:50:09 +0200 Subject: [PATCH 10/35] Update docs --- docs/api_analysis.md | 5 ++++- docs/api_muap.md | 7 +++++++ docs/api_pic.md | 3 +++ mkdocs.yml | 1 + 4 files changed, 15 insertions(+), 1 deletion(-) create mode 100644 docs/api_pic.md diff --git a/docs/api_analysis.md b/docs/api_analysis.md index efe84d9..6371c3a 100644 --- a/docs/api_analysis.md +++ b/docs/api_analysis.md @@ -1,6 +1,9 @@ Description ----------- -This module contains all the functions used to analyse the MUs properties when not involving the MUs action potential shape. +This module contains all the functions used to analyse the MUs properties apart from those involving: + +- The MUs action potential shape => Available in the [muap module](api_muap.md) +- DeltaF / PICs estimation => available in the [pic module](api_pic.md)
diff --git a/docs/api_muap.md b/docs/api_muap.md index b7e8fb3..5a1707d 100644 --- a/docs/api_muap.md +++ b/docs/api_muap.md @@ -68,6 +68,13 @@ This module contains functions to produce and analyse MUs anction potentials
+::: openhdemg.library.muap.Tracking_gui + options: + show_root_full_path: False + show_root_heading: True + +
+ ::: openhdemg.library.muap.remove_duplicates_between options: show_root_full_path: False diff --git a/docs/api_pic.md b/docs/api_pic.md new file mode 100644 index 0000000..38d69ce --- /dev/null +++ b/docs/api_pic.md @@ -0,0 +1,3 @@ +Description +----------- +This \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index b537729..750bd91 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -83,6 +83,7 @@ nav: - openfiles: api_openfiles.md - plotemg: api_plotemg.md - analysis: api_analysis.md + - pic: api_pic.md - tools: api_tools.md - mathtools: api_mathtools.md - muap: api_muap.md From d241a023837085194d45e25e67c4aedadb8e0fc9 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Mon, 3 Jun 2024 11:50:15 +0200 Subject: [PATCH 11/35] Update backup_settings.py --- openhdemg/gui/backup_settings.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/openhdemg/gui/backup_settings.py b/openhdemg/gui/backup_settings.py index fb335ea..b14f761 100644 --- a/openhdemg/gui/backup_settings.py +++ b/openhdemg/gui/backup_settings.py @@ -86,6 +86,9 @@ # in tracking() tracking__firings = "all" tracking__derivation = "sd" +tracking__gui = True +tracking__gui_addrefsig = True +tracking__gui_csv_separator = "\t" # in remove_duplicates_between() remove_duplicates_between__firings = "all" @@ -95,6 +98,7 @@ MUcv_gui__n_firings = [0, 50] MUcv_gui__muaps_timewindow = 50 MUcv_gui__figsize = [25, 20] +MUcv_gui__csv_separator = "\t" # --------------------------------- electrodes -------------------------------- # This custom sorting order is valid for all the GUI windows, although the From 026e2953993654dbb582930894ff6f9bb0b90e10 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Mon, 3 Jun 2024 15:57:10 +0200 Subject: [PATCH 12/35] Format compute_svr and compute_deltaf Format functions for PEP 8 agreement and documentation building. --- openhdemg/library/pic.py | 432 ++++++++++++++++++++++--------------- openhdemg/library/tools.py | 242 +++++++++++++-------- 2 files changed, 407 insertions(+), 267 deletions(-) diff --git a/openhdemg/library/pic.py b/openhdemg/library/pic.py index ef7aa90..feb665f 100644 --- a/openhdemg/library/pic.py +++ b/openhdemg/library/pic.py @@ -1,88 +1,102 @@ """ -This module contains all the functions used to quantify and analyze MU persistent inward currents. +This module contains all the functions used to quantify and analyze MU +persistent inward currents. + Currently includes delta F. """ import pandas as pd import numpy as np -from openhdemg.library.tools import compute_svr from itertools import combinations def compute_deltaf( - emgfile, - smoothfits, - average_method="test_unit_average", - normalization="False", - recruitment_difference_cutoff=1.0, - corr_cutoff=0.7, - controlunitmodulation_cutoff=0.5, - clean="True" - ): - + emgfile, + smoothfits, + average_method="test_unit_average", + normalisation="False", + recruitment_difference_cutoff=1.0, + corr_cutoff=0.7, + controlunitmodulation_cutoff=0.5, + clean="True", +): """ - Conducts a paired motor unit analysis, quantifying delta F bettwen the supplied collection of motor units. - Origional framework for deltaF provided in Gorassini et. al., 2002 : https://journals.physiology.org/doi/full/10.1152/jn.00024.2001 - James (Drew) Beauchamp - - + Quantify delta F via paired motor unit analysis. + + Conducts a paired motor unit analysis, quantifying delta F between the + supplied collection of motor units. Origional framework for deltaF provided + in Gorassini et. al., 2002: + https://journals.physiology.org/doi/full/10.1152/jn.00024.2001 + + Author: James (Drew) Beauchamp + Parameters ---------- emgfile : dict The dictionary containing the emgfile. - smoothfits : list of arrays Smoothed discharge rate estimates. - Each array: motor unit discharge rate x samples aligned in time; instances of non-firing = NaN, please + Each array: motor unit discharge rate x samples aligned in time; + instances of non-firing = NaN, please #TODO maybe something is missing after please? Your choice of smoothing. See compute_svr gen_svr for example. + average_method : str {"test_unit_average", "all"}, default "test_unit_average" + The method for test MU deltaF value. More to be added. TODO # TODOs should be moved outside documentation or they will be displayed in the website - average_method : str, default "test_unit_average" - The method for test MU deltaF value. More to be added TODO - - "test_unit_average" - The average across all possible control units + ``test_unit_average`` + The average across all possible control units. - normalization : str, {"False", "ctrl_max_desc"}, default "False" - "ctrl_max_desc" - Whether to normalize deltaF values to control unit descending range during test unit firing - See Skarabot et. al., 2023 : https://www.biorxiv.org/content/10.1101/2023.10.16.562612v1 + ``all`` TODO the user is not aware of what other average methods can be used + abc + normalisation : str {"False", "ctrl_max_desc"}, default "False" + The method for deltaF nomalization. + ``ctrl_max_desc`` + Whether to normalise deltaF values to control unit descending + range during test unit firing. See Skarabot et. al., 2023: + https://www.biorxiv.org/content/10.1101/2023.10.16.562612v1 recruitment_difference_cutoff : float, default 1 - An exlusion criteria corresponding to the necessary difference between control and test MU recruitement in seconds - corr_cutoff : float, (0,1), default 0.7 - An exclusion criteria corresponding to the correlation between control and test unit discharge rate. + An exlusion criteria corresponding to the necessary difference between + control and test MU recruitement in seconds. + corr_cutoff : float (0 to 1), default 0.7 + An exclusion criteria corresponding to the correlation between control + and test unit discharge rate. controlunitmodulation_cutoff : float, default 0.5 - An exclusion criteria corresponding to the necessary modulation of control unit discharge rate during test unit firing in Hz. + An exclusion criteria corresponding to the necessary modulation of + control unit discharge rate during test unit firing in Hz. clean : str, default "True" To remove values that do not meet exclusion criteria - - + Returns ------- - dF : pd.DataFrame - A pd.DataFrame containing deltaF values and corresponding MU number + delta_f : pd.DataFrame + A pd.DataFrame containing deltaF values and corresponding MU number. + TODO here we should explain that the resulting df will be different + depending on average_method. In particular, if average_method="all", + delta_f[MU][row] will contain a tuple representing... instead of a MU + number. + + See also + -------- + - compute_svr : fit MU discharge rates with Support Vector Regression, + nonlinear regression. Warnings -------- - TODO - - Example - ----- - import openhdemg.library as emg - import numpy as np - import matplotlib.pyplot as plt + TODO consider removing section from docstring while not implemented - # Load the sample file - emgfile = emg.emg_from_samplefile() - - # Sort MUs based on recruitment order - emgfile = emg.sort_mus(emgfile=emgfile) - - # quantify svr fits - svrfits = emg.tools.compute_svr(emgfile) - - # quantify delta F - dF = emg.pic.compute_deltaf(emgfile=emgfile, smoothfits=svrfits["gensvr"]) - - dF + Examples + -------- + Quantify delta F using svr fits. + + >>> import openhdemg.library as emg + >>> import numpy as np + >>> emgfile = emg.emg_from_samplefile() + >>> emgfile = emg.sort_mus(emgfile=emgfile) + >>> svrfits = emg.compute_svr(emgfile) + >>> delta_f = emg.compute_deltaf( + ... emgfile=emgfile, smoothfits=svrfits["gensvr"], + ... ) + delta_f MU dF 0 0 NaN 1 1 NaN @@ -90,10 +104,13 @@ def compute_deltaf( 3 3 1.838382 4 4 2.709522 - # for all possible combinations, not test unit average, MU in this case is pairs (reporter,test) - dF2 = emg.pic.compute_deltaf(emgfile=emgfile, smoothfits=svrfits["gensvr"],average_method='all') + For all possible combinations, not test unit average, MU in this case is + pairs (reporter, test). - dF2 + >>> delta_f_2 = emg.compute_deltaf( + ... emgfile=emgfile, smoothfits=svrfits["gensvr"], average_method='all', + ... ) + delta_f_2 MU dF 0 (0, 1) NaN 1 (0, 2) NaN @@ -105,143 +122,204 @@ def compute_deltaf( 7 (2, 3) NaN 8 (2, 4) NaN 9 (3, 4) 2.709522 - - """ - # svrfits = compute_SVR(emgfile) # use smooth svr fits dfret_ret = [] mucombo_ret = np.empty(0, int) - if emgfile["NUMBER_OF_MUS"] < 2: # if less than 2 MUs, can not quantify deltaF + + # If less than 2 MUs, can not quantify deltaF + if emgfile["NUMBER_OF_MUS"] < 2: dfret_ret = np.nan - mucombo_ret = np.nan*np.ones(1, 2) - else: - combs = combinations(range(emgfile["NUMBER_OF_MUS"]), 2) - # combinations of MUs - # TODO if units are nonconsecutive - - # init - r_ret = [] - dfret = [] - testmu = [] - ctrl_mod = [] - mucombo = [] - rcrt_diff = [] - controlmu = [] - for mucomb in list(combs): # for all possible combinations of MUs - mu1_id, mu2_id = mucomb[0], mucomb[1] # extract possible MU combinations (a unique MU pair) - mucombo.append((mu1_id, mu2_id)) # track current MU combination - - # first MU firings, recruitment, and decrecruitment - mu1_times = np.where(emgfile["BINARY_MUS_FIRING"][mu1_id] == 1)[0] - mu1_rcrt, mu1_drcrt = mu1_times[1], mu1_times[-1] # skip first since idr is defined on second - - # second MU firings, recruitment, and decrecruitment - mu2_times = np.where(emgfile["BINARY_MUS_FIRING"][mu2_id]==1)[0] - mu2_rcrt, mu2_drcrt = mu2_times[1], mu2_times[-1] # skip first since idr is defined on second - - # region of MU overlap - muoverlap = range(max(mu1_rcrt, mu2_rcrt), min(mu1_drcrt, mu2_drcrt)) - - # if MUs do not overlapt by more than two or more samples - if len(muoverlap) < 2: - dfret = np.append(dfret, np.nan) - r_ret = np.append(r_ret, np.nan) - rcrt_diff = np.append(rcrt_diff, np.nan) - ctrl_mod = np.append(ctrl_mod, np.nan) - continue # TODO test - - # corr between units - not always necessary, can be set to 0 when desired - r = pd.DataFrame(zip(smoothfits[mu1_id][muoverlap], smoothfits[mu2_id][muoverlap])).corr() - r_ret = np.append(r_ret, r[0][1]) - - # recruitment diff, necessary to ensure PICs are activated in control unit - rcrt_diff = np.append(rcrt_diff, np.abs(mu1_rcrt-mu2_rcrt)/emgfile["FSAMP"]) - if mu1_rcrt < mu2_rcrt: + mucombo_ret = np.nan*np.ones(1, 2) # TODO this line is not working, please replace with nan or what appropriate + + delta_f = pd.DataFrame({'MU': mucombo_ret, 'dF': dfret_ret}) + + return delta_f + + # If more than 2 MUs, quantify deltaF. + # Combinations of MUs. + combs = combinations(range(emgfile["NUMBER_OF_MUS"]), 2) + # TODO if units are nonconsecutive + + # init + r_ret = [] + dfret = [] + testmu = [] + ctrl_mod = [] + mucombo = [] + rcrt_diff = [] + controlmu = [] + for mucomb in list(combs): # For all possible combinations of MUs + # Extract possible MU combinations (a unique MU pair) + mu1_id, mu2_id = mucomb[0], mucomb[1] + # Track current MU combination + mucombo.append((mu1_id, mu2_id)) + + # First MU firings, recruitment, and decrecruitment + mu1_times = np.where(emgfile["BINARY_MUS_FIRING"][mu1_id] == 1)[0] + mu1_rcrt, mu1_drcrt = mu1_times[1], mu1_times[-1] + # Skip first since idr is defined on second + + # Second MU firings, recruitment, and decrecruitment + mu2_times = np.where(emgfile["BINARY_MUS_FIRING"][mu2_id] == 1)[0] + mu2_rcrt, mu2_drcrt = mu2_times[1], mu2_times[-1] + # Skip first since idr is defined on second + + # Region of MU overlap + muoverlap = range( + max(mu1_rcrt, mu2_rcrt), min(mu1_drcrt, mu2_drcrt), + ) + + # If MUs do not overlapt by more than two or more samples + if len(muoverlap) < 2: + dfret = np.append(dfret, np.nan) + r_ret = np.append(r_ret, np.nan) + rcrt_diff = np.append(rcrt_diff, np.nan) + ctrl_mod = np.append(ctrl_mod, np.nan) + continue # TODO test + + # Corr between units - not always necessary, can be set to 0 when + # desired. + r = pd.DataFrame( + zip( + smoothfits[mu1_id][muoverlap], + smoothfits[mu2_id][muoverlap], + ) + ).corr() + r_ret = np.append(r_ret, r[0][1]) + + # Recruitment diff, necessary to ensure PICs are activated in + # control unit. + rcrt_diff = np.append( + rcrt_diff, np.abs(mu1_rcrt-mu2_rcrt)/emgfile["FSAMP"], + ) + if mu1_rcrt < mu2_rcrt: + controlU = 1 # MU 1 is control unit, 2 is test unit + + # If control (reporter) unit is not on for entirety of test + # unit, set last firing to control unit. + if mu1_drcrt < mu2_drcrt: + mu2_drcrt = mu1_drcrt + # This may understimate PICs, other methods can be employed + # delta F: change in control MU discharge rate between test + # unit recruitment and derecruitment. + df = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu2_drcrt] + + # Control unit discharge rate modulation while test unit is + # firing. + ctrl_mod = np.append( + ctrl_mod, + np.nanmax(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)]) + - np.nanmin(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)]), + ) + + if normalisation == "False": + dfret = np.append(dfret, df) + elif normalisation == "ctrl_max_desc": + # Normalise deltaF values to control unit descending range + # during test unit firing. + k = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu1_drcrt] + dfret = np.append(dfret, df/k) + + elif mu1_rcrt > mu2_rcrt: + controlU = 2 # MU 2 is control unit, 1 is test unit + if mu1_drcrt > mu2_drcrt: + # If control (reporter) unit is not on for entirety of + # test unit, set last firing to control unit. + mu1_drcrt = mu2_drcrt + # This may understimate PICs, other methods can be employed. + # delta F: change in control MU discharge rate between test + # unit recruitment and derecruitment. + df = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu1_drcrt] + + # Control unit discharge rate modulation while test unit is + # firing. + ctrl_mod = np.append( + ctrl_mod, + np.nanmax(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)]) + - np.nanmin(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)]), + ) + + if normalisation == "False": + dfret = np.append(dfret, df) + elif normalisation == "ctrl_max_desc": + # Normalise deltaF values to control unit descending range + # during test unit firing. + k = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu2_drcrt] + dfret = np.append(dfret, df/k) + + elif mu1_rcrt == mu2_rcrt: + if mu1_drcrt > mu2_drcrt: controlU = 1 # MU 1 is control unit, 2 is test unit - - if mu1_drcrt < mu2_drcrt: # if control (reporter) unit is not on for entirety of test unit, set last firing to control unit - mu2_drcrt = mu1_drcrt - # this may understimate PICs, other methods can be employed - # delta F: change in control MU discharge rate between test unit recruitment and derecruitment + # delta F: change in control MU discharge rate between + # test unit recruitment and derecruitment. df = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu2_drcrt] - # control unit discharge rate modulation while test unit is firing - ctrl_mod = np.append(ctrl_mod, np.nanmax(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)])-np.nanmin(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)])) + # Control unit discharge rate modulation while test unit is + # firing. + ctrl_mod = np.append( + ctrl_mod, + np.nanmax(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)]) + - np.nanmin(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)]), + ) - if normalization == "False": + if normalisation == "False": dfret = np.append(dfret, df) - elif normalization == "ctrl_max_desc": # normalize deltaF values to control unit descending range during test unit firing + elif normalisation == "ctrl_max_desc": + # Normalise deltaF values to control unit descending + # range during test unit firing. k = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu1_drcrt] dfret = np.append(dfret, df/k) - - elif mu1_rcrt > mu2_rcrt: + else: controlU = 2 # MU 2 is control unit, 1 is test unit - if mu1_drcrt > mu2_drcrt: # if control (reporter) unit is not on for entirety of test unit, set last firing to control unit - mu1_drcrt = mu2_drcrt - # this may understimate PICs, other methods can be employed - # delta F: change in control MU discharge rate between test unit recruitment and derecruitment + # delta F: change in control MU discharge rate between + # test unit recruitment and derecruitment. df = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu1_drcrt] - # control unit discharge rate modulation while test unit is firing - ctrl_mod = np.append(ctrl_mod, np.nanmax(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)])-np.nanmin(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)])) + # Control unit discharge rate modulation while test unit is + # firing. + ctrl_mod = np.append( + ctrl_mod, + np.nanmax(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)]) + - np.nanmin(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)]), + ) - if normalization == "False": + if normalisation == "False": dfret = np.append(dfret, df) - elif normalization == "ctrl_max_desc": # normalize deltaF values to control unit descending range during test unit firing + elif normalisation == "ctrl_max_desc": + # Normalise deltaF values to control unit descending + # range during test unit firing. k = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu2_drcrt] dfret = np.append(dfret, df/k) - elif mu1_rcrt == mu2_rcrt: - if mu1_drcrt > mu2_drcrt: - controlU = 1 # MU 1 is control unit, 2 is test unit - # delta F: change in control MU discharge rate between test unit recruitment and derecruitment - df = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu2_drcrt] - - # control unit discharge rate modulation while test unit is firing - ctrl_mod = np.append(ctrl_mod, np.nanmax(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)])-np.nanmin(smoothfits[mu1_id][range(mu2_rcrt, mu2_drcrt)])) - - if normalization == "False": - dfret = np.append(dfret, df) - elif normalization == "ctrl_max_desc": # normalize deltaF values to control unit descending range during test unit firing - k = smoothfits[mu1_id][mu2_rcrt]-smoothfits[mu1_id][mu1_drcrt] - dfret = np.append(dfret, df/k) - else: - controlU = 2 # MU 2 is control unit, 1 is test unit - # delta F: change in control MU discharge rate between test unit recruitment and derecruitment - df = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu1_drcrt] - - # control unit discharge rate modulation while test unit is firing - ctrl_mod = np.append(ctrl_mod, np.nanmax(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)])-np.nanmin(smoothfits[mu2_id][range(mu1_rcrt, mu1_drcrt)])) - - if normalization == "False": - dfret = np.append(dfret, df) - elif normalization == "ctrl_max_desc": # normalize deltaF values to control unit descending range during test unit firing - k = smoothfits[mu2_id][mu1_rcrt]-smoothfits[mu2_id][mu2_drcrt] - dfret = np.append(dfret, df/k) - - # collect which MUs were control vs test - controlmu.append(mucombo[-1][controlU-1]) - testmu.append(mucombo[-1][1-controlU//2]) - - if clean == "True": # remove values that dont meet exclusion criteria - rcrt_diff_bin = rcrt_diff > recruitment_difference_cutoff - corr_bin = r_ret > corr_cutoff - ctrl_mod_bin = ctrl_mod > controlunitmodulation_cutoff - clns = np.asarray([rcrt_diff_bin & corr_bin & ctrl_mod_bin]) - dfret[~clns[0]] = np.nan - - if average_method == "test_unit_average": # average across all control units - for ii in range(emgfile["NUMBER_OF_MUS"]): - clean_indices = [index for (index, item) in enumerate(testmu) if item == ii] - if np.isnan(dfret[clean_indices]).all(): - dfret_ret = np.append(dfret_ret, np.nan) - else: - dfret_ret = np.append(dfret_ret, np.nanmean(dfret[clean_indices])) - mucombo_ret = np.append(mucombo_ret, int(ii)) - else: # return all values and corresponding combinations - dfret_ret = dfret - mucombo_ret = mucombo - - dF = pd.DataFrame({'MU': mucombo_ret, 'dF': dfret_ret}) - return dF + # Collect which MUs were control vs test + controlmu.append(mucombo[-1][controlU-1]) + testmu.append(mucombo[-1][1-controlU//2]) + + if clean == "True": # Remove values that dont meet exclusion criteria + rcrt_diff_bin = rcrt_diff > recruitment_difference_cutoff + corr_bin = r_ret > corr_cutoff + ctrl_mod_bin = ctrl_mod > controlunitmodulation_cutoff + clns = np.asarray([rcrt_diff_bin & corr_bin & ctrl_mod_bin]) + dfret[~clns[0]] = np.nan + + if average_method == "test_unit_average": + # Average across all control units + for ii in range(emgfile["NUMBER_OF_MUS"]): + clean_indices = [ + index for (index, item) in enumerate(testmu) if item == ii + ] + if np.isnan(dfret[clean_indices]).all(): + dfret_ret = np.append(dfret_ret, np.nan) + else: + dfret_ret = np.append( + dfret_ret, np.nanmean(dfret[clean_indices]), + ) + mucombo_ret = np.append(mucombo_ret, int(ii)) + else: # Return all values and corresponding combinations + dfret_ret = dfret + mucombo_ret = mucombo + + delta_f = pd.DataFrame({'MU': mucombo_ret, 'dF': dfret_ret}) + + return delta_f diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index 7ed3ede..5ae2312 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -10,12 +10,11 @@ import numpy as np import matplotlib.pyplot as plt from scipy import signal +from scipy.stats import iqr +from sklearn.svm import SVR import warnings from openhdemg.library.mathtools import compute_sil -from sklearn.svm import SVR -from scipy.stats import iqr - def showselect(emgfile, how="ref_signal", title="", titlesize=12, nclic=2): """ @@ -147,7 +146,7 @@ def create_binary_firings(emg_length, number_of_mus, mupulses): if not isinstance(mupulses, list): raise ValueError("mupulses is not a list of ndarrays") - # Initialize a pd.DataFrame with zeros + # Initialise a pd.DataFrame with zeros binary_MUs_firing = pd.DataFrame( np.zeros((emg_length, number_of_mus), dtype=int) ) @@ -1201,19 +1200,23 @@ def compute_rfd( def compute_svr( - emgfile, - gammain=1/1.6, - regparam=1/0.370, - endpointweights_numpulses=5, - endpointweights_magnitude=5, - discontfiring_dur=1.0 - ): + emgfile, + gammain=1/1.6, + regparam=1/0.370, + endpointweights_numpulses=5, + endpointweights_magnitude=5, + discontfiring_dur=1.0, +): """ - Fit MU Discharge rates with Support Vector Regression, nonlinear regression - Provides smooth and continous estimates of discharge rate usefull for quantification and visualization - Suggested hyperparameters and framework from Beauchamp et. al., 2022 + Fit MU discharge rates with Support Vector Regression, nonlinear + regression. + + Provides smooth and continous estimates of discharge rate useful for + quantification and visualisation. Suggested hyperparameters and framework + from Beauchamp et. al., 2022 https://doi.org/10.1088/1741-2552/ac4594 - James (Drew) Beauchamp + + Author: James (Drew) Beauchamp Parameters ---------- @@ -1224,121 +1227,180 @@ def compute_svr( regparam : float, default 1/0.370 The regularization parameter, must be positive. endpointweights_numpulses : int, default 5 - Number of discharge instances at the start and end of MU firing to apply a weighting coefficient. + Number of discharge instances at the start and end of MU firing to + apply a weighting coefficient. endpointweights_magnitude : int, default 5 - The scaling factor applied to the number of pulses provided by endpointweights_numpulses. + The scaling factor applied to the number of pulses provided by + endpointweights_numpulses. The scaling is applied to the regularization parameter, per sample. - Larger values force the classifier to put more emphasis on the number of discharge instances at the start and end of firing provided by endpointweights_numpulses + Larger values force the classifier to put more emphasis on the number + of discharge instances at the start and end of firing provided by + endpointweights_numpulses. discontfiring_dur : int, default 1 - Duration of time in seconds that defines an instnance of discontinuous firing. SVR fits will not be returned at points of discontinuity. + Duration of time in seconds that defines an instnance of discontinuous + firing. SVR fits will not be returned at points of discontinuity. Returns ------- svrfits : pd.DataFrame - A pd.DataFrame containing the smooth/continous MU discharge rates and corresponding time vetors + A pd.DataFrame containing the smooth/continous MU discharge rates and + corresponding time vetors. + + See also + -------- + - compute_deltaf : quantify delta F via paired motor unit analysis. Warnings -------- - TODO + TODO consider removing section from docstring while not implemented - Example - ----- - import openhdemg.library as emg - import numpy as np - import matplotlib.pyplot as plt - - # Load the sample file - emgfile = emg.emg_from_samplefile() - - # Sort MUs based on recruitment order - emgfile = emg.sort_mus(emgfile=emgfile) - - # quantify svr fits - svrfits = emg.tools.compute_svr(emgfile) - - # quick plot showing results - scl = 12 # offset MUs for viz - idr = emg.compute_idr(emgfile) - for ii in range(len(svrfits["svrfit"])): - xtmp = np.transpose([idr[ii].timesec[1:]]) - ytmp = idr[ii].idr[1:].to_numpy() - plt.scatter(xtmp, ytmp+scl*(ii), color='darkorange') - plt.plot(svrfits["svrtime"][ii], svrfits["svrfit"][ii]+scl*(ii), color='cornflowerblue') - plt.show() + Examples + -------- + Quantify svr fits. + + >>> import openhdemg.library as emg + >>> import numpy as np + >>> import matplotlib.pyplot as plt + >>> emgfile = emg.emg_from_samplefile() + >>> emgfile = emg.sort_mus(emgfile=emgfile) + >>> svrfits = emg.compute_svr(emgfile) + + Quick plot showing the results. + + >>> scl = 12 # offset MUs for viz + >>> idr = emg.compute_idr(emgfile) + >>> for ii in range(len(svrfits["svrfit"])): + ... xtmp = np.transpose([idr[ii].timesec[1:]]) + ... ytmp = idr[ii].idr[1:].to_numpy() + ... plt.scatter(xtmp, ytmp+scl*(ii), color='darkorange') + ... plt.plot( + ... svrfits["svrtime"][ii], + ... svrfits["svrfit"][ii]+scl*(ii), + ... color='cornflowerblue', + ... ) + >>> plt.show() """ # TODO input checking and edge cases - idr = compute_idr(emgfile) # calc IDR + idr = compute_idr(emgfile) # Calc IDR svrfit_acm = [] svrtime_acm = [] gensvr_acm = [] - for mu in range(len(idr)): # for all MUs - # train the model on the data - xtmp = np.transpose([idr[mu].timesec[1:]]) # time vector, removing first element - ytmp = idr[mu].idr[1:].to_numpy() # discharge rates, removing first element, since DR has been assigned to second pulse - xdiff = idr[mu].diff_mupulses[1:].values # time between discharges, will use for discontinuity calc - mup = np.array(idr[mu].mupulses[1:]) # motor unit pulses, samples - - # Defining weight vector. A scaling applied to the regularization parameter, per sample. + for mu in range(len(idr)): # For all MUs + # Train the model on the data. + # Time vector, removing first element. + xtmp = np.transpose([idr[mu].timesec[1:]]) + # Discharge rates, removing first element, since DR has been assigned + # to second pulse. + ytmp = idr[mu].idr[1:].to_numpy() + # Time between discharges, will use for discontinuity calc + xdiff = idr[mu].diff_mupulses[1:].values + # Motor unit pulses, samples + mup = np.array(idr[mu].mupulses[1:]) + + # Defining weight vector. A scaling applied to the regularization + # parameter, per sample. smpwht = np.ones(len(ytmp)) smpwht[0:endpointweights_numpulses-1] = endpointweights_magnitude smpwht[(len(ytmp)-(endpointweights_numpulses-1)):len(ytmp)] = endpointweights_magnitude - # create an SVR model with a gausian kernel and supplied hyperparams - # origional hyperparameters from Beauchamp et. al., 2022 https://doi.org/10.1088/1741-2552/ac4594 - svr = SVR(kernel='rbf', gamma=gammain, C=np.abs(regparam), epsilon=iqr(ytmp)/11) + # Create an SVR model with a gausian kernel and supplied hyperparams. + # Origional hyperparameters from Beauchamp et. al., 2022: + # https://doi.org/10.1088/1741-2552/ac4594 + svr = SVR( + kernel='rbf', gamma=gammain, C=np.abs(regparam), + epsilon=iqr(ytmp)/11, + ) svr.fit(xtmp, ytmp, sample_weight=smpwht) # Defining prediction vector - # TODO need to add custom range - predind = np.arange(mup[0], mup[-1]+1) # from the second firing to the end of firing, in samples - predtime = (predind/emgfile["FSAMP"]).reshape(-1, 1) # in time (s) - gen_svr = np.nan*np.ones(emgfile["EMG_LENGTH"]) # initialize nan vector for tracking fits aligned in time. usefull for later quant metrics - - # check for discontinous firing - bkpnt = mup[np.where((xdiff > (discontfiring_dur * emgfile["FSAMP"])))[0]] - - # make predictions on the data - if bkpnt.size > 0: # if there is a point of discontinuity + # TODO need to add custom range. + # From the second firing to the end of firing, in samples. + predind = np.arange(mup[0], mup[-1]+1) + predtime = (predind/emgfile["FSAMP"]).reshape(-1, 1) # In time (s) + # Initialise nan vector for tracking fits aligned in time. Usefull for + # later quant metrics. + gen_svr = np.nan*np.ones(emgfile["EMG_LENGTH"]) + + # Check for discontinous firing + bkpnt = mup[ + np.where((xdiff > (discontfiring_dur * emgfile["FSAMP"])))[0] + ] + + # Make predictions on the data + if bkpnt.size > 0: # If there is a point of discontinuity if bkpnt[0] == mup[0]: smoothfit = np.nan*np.ones(1) newtm = np.nan*np.ones(1) bkpnt = bkpnt[1:] else: - tmptm = predtime[0:np.where((bkpnt[0] >= predind[0:-1]) & (bkpnt[0] < predind[1:]))[0][0]] # break up time vector for first continous range of firing - smoothfit = svr.predict(tmptm) # predict with svr model - newtm = tmptm # track new time vector - - tmpind = predind[0:np.where((bkpnt[0] >= predind[0:-1]) & ( bkpnt[0]= predind[0:-1]) & (bkpnt[ii] < predind[1:]))[0][0] # current index of discontinuity - nextind = np.where((bkpnt[ii+1] >= predind[0:-1]) & (bkpnt[ii+1] <= predind[1:]))[0][0] # next index of discontinuity - - curmup = np.where(mup == bkpnt[ii])[0][0] # MU firing before discontinuity - curind_nmup = np.where((mup[curmup+1] >= predind[0:-1]) & (mup[curmup+1] <= predind[1:]))[0][0] # MU firing after discontinuity - - if curind_nmup >= nextind: # if the next discontinuity is the next MU firing, nan fill - # edge case NEED TO CHECK THE GREATER THAN CASE>> WHY TODO + tmptm = predtime[ + 0: np.where( + (bkpnt[0] >= predind[0:-1]) & (bkpnt[0] < predind[1:]) + )[0][0] + ] # Break up time vector for first continous range of firing + smoothfit = svr.predict(tmptm) # Predict with svr model + newtm = tmptm # Track new time vector + + tmpind = predind[ + 0: np.where( + (bkpnt[0] >= predind[0:-1]) & (bkpnt[0] < predind[1:]) + )[0][0] + ] # Sample vector of first continous range of firing + # Fill corresponding sample indices with svr fit + gen_svr[tmpind.astype(np.int64)] = smoothfit + # Add last firing as discontinuity + bkpnt = np.append(bkpnt, mup[-1]) + for ii in range(bkpnt.size-1): # All instances of discontinuity + curind = np.where( + (bkpnt[ii] >= predind[0:-1]) & (bkpnt[ii] < predind[1:]) + )[0][0] # Current index of discontinuity + nextind = np.where( + (bkpnt[ii+1] >= predind[0:-1]) & (bkpnt[ii+1] <= predind[1:]) + )[0][0] # Next index of discontinuity + + # MU firing before discontinuity + curmup = np.where(mup == bkpnt[ii])[0][0] + curind_nmup = np.where( + (mup[curmup+1] >= predind[0:-1]) & (mup[curmup+1] <= predind[1:]) + )[0][0] # MU firing after discontinuity + + # If the next discontinuity is the next MU firing, nan fill + if curind_nmup >= nextind: + # Edge case NEED TO CHECK THE GREATER THAN CASE>> WHY TODO smoothfit = np.append(smoothfit, np.nan*np.ones(1)) newtm = np.append(newtm, np.nan*np.ones(1)) - else: # fit next continuous region of firing - smoothfit = np.append(smoothfit, np.nan*np.ones(len(predtime[curind:curind_nmup])-2)) - smoothfit = np.append(smoothfit, svr.predict(predtime[curind_nmup:nextind])) - newtm = np.append(newtm, np.nan*np.ones(len(predtime[curind:curind_nmup])-2)) + else: # Fit next continuous region of firing + smoothfit = np.append( + smoothfit, + np.nan*np.ones(len(predtime[curind:curind_nmup])-2), + ) + smoothfit = np.append( + smoothfit, svr.predict(predtime[curind_nmup:nextind]), + ) + newtm = np.append( + newtm, + np.nan*np.ones(len(predtime[curind:curind_nmup])-2), + ) newtm = np.append(newtm, predtime[curind_nmup:nextind]) - gen_svr[predind[curind_nmup:nextind]] = svr.predict(predtime[curind_nmup:nextind]) + gen_svr[predind[curind_nmup:nextind]] = svr.predict( + predtime[curind_nmup:nextind] + ) else: smoothfit = svr.predict(predtime) newtm = predtime gen_svr[predind] = smoothfit + # Append fits, new time vect, time aligned fits svrfit_acm.append(smoothfit.copy()) svrtime_acm.append(newtm.copy()) - gensvr_acm.append(gen_svr.copy()) # append fits, new time vect, time aligned fits + gensvr_acm.append(gen_svr.copy()) - svrfits = ({"svrfit": svrfit_acm, "svrtime": svrtime_acm, "gensvr": gensvr_acm}) + svrfits = { + "svrfit": svrfit_acm, + "svrtime": svrtime_acm, + "gensvr": gensvr_acm, + } return svrfits From ecfbb4521049e084ee539d352939d7a3d4c00087 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Mon, 3 Jun 2024 15:58:08 +0200 Subject: [PATCH 13/35] Add compute_svr and compute_deltaf --- docs/api_pic.md | 13 ++++++++++++- docs/api_tools.md | 7 +++++++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/docs/api_pic.md b/docs/api_pic.md index 38d69ce..d2bdfd8 100644 --- a/docs/api_pic.md +++ b/docs/api_pic.md @@ -1,3 +1,14 @@ Description ----------- -This \ No newline at end of file +This module contains all the functions used to quantify and analyze MU persistent inward currents. + +Currently includes delta F. + +
+ +::: openhdemg.library.pic.compute_deltaf + options: + show_root_full_path: False + show_root_heading: True + +
\ No newline at end of file diff --git a/docs/api_tools.md b/docs/api_tools.md index f65b7bd..b421ad0 100644 --- a/docs/api_tools.md +++ b/docs/api_tools.md @@ -104,3 +104,10 @@ shortcuts necessary to operate with the HD-EMG recordings. show_root_heading: True
+ +::: openhdemg.library.tools.compute_svr + options: + show_root_full_path: False + show_root_heading: True + +
From 3c2383d7d5ef4fa968ac8cd5032ecb185cba44c2 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 4 Jun 2024 08:13:13 +0200 Subject: [PATCH 14/35] Update pic.py --- openhdemg/library/pic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openhdemg/library/pic.py b/openhdemg/library/pic.py index feb665f..7f22db8 100644 --- a/openhdemg/library/pic.py +++ b/openhdemg/library/pic.py @@ -18,7 +18,7 @@ def compute_deltaf( recruitment_difference_cutoff=1.0, corr_cutoff=0.7, controlunitmodulation_cutoff=0.5, - clean="True", + clean="True", # TODO should this be a boolean True? ): """ Quantify delta F via paired motor unit analysis. From 6c0b1ecabd2bcc6b6004b6e7bbc0c34a2378451c Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 4 Jun 2024 11:06:07 +0200 Subject: [PATCH 15/35] Add PICs in GUI --- openhdemg/gui/backup_settings.py | 44 ++- .../gui/gui_modules/advanced_analyses.py | 273 ++++++++++++++++-- openhdemg/gui/settings.py | 44 ++- 3 files changed, 301 insertions(+), 60 deletions(-) diff --git a/openhdemg/gui/backup_settings.py b/openhdemg/gui/backup_settings.py index b14f761..d70adff 100644 --- a/openhdemg/gui/backup_settings.py +++ b/openhdemg/gui/backup_settings.py @@ -19,30 +19,31 @@ import numpy as np + # --------------------------------- openfiles --------------------------------- -# in emg_from_demuse() +# In emg_from_demuse() emg_from_demuse__ignore_negative_ipts = False -# in emg_from_otb() +# In emg_from_otb() emg_from_otb__ext_factor = 8 emg_from_otb__refsig = [True, "fullsampled"] emg_from_otb__extras = None emg_from_otb__ignore_negative_ipts = False -# in refsig_from_otb() +# In refsig_from_otb() refsig_from_otb__refsig = "fullsampled" refsig_from_otb__extras = None -# in emg_from_delsys() +# In emg_from_delsys() emg_from_delsys__emg_sensor_name = "Galileo sensor" emg_from_delsys__refsig_sensor_name = "Trigno Load Cell" emg_from_delsys__filename_from = "mus_directory" -# in refsig_from_delsys() +# In refsig_from_delsys() refsig_from_delsys__refsig_sensor_name = "Trigno Load Cell" -# in emg_from_customcsv() +# In emg_from_customcsv() emg_from_customcsv__ref_signal = "REF_SIGNAL" emg_from_customcsv__raw_signal = "RAW_SIGNAL" emg_from_customcsv__ipts = "IPTS" @@ -53,21 +54,21 @@ emg_from_customcsv__fsamp = 2048 emg_from_customcsv__ied = 8 -# in refsig_from_customcsv() +# In refsig_from_customcsv() refsig_from_customcsv__ref_signal = "REF_SIGNAL" refsig_from_customcsv__extras = "EXTRAS" refsig_from_customcsv__fsamp = 2048 -# in save_json_emgfile() +# In save_json_emgfile() save_json_emgfile__compresslevel = 4 # ---------------------------------- analysis --------------------------------- -# in compute_thresholds() +# In compute_thresholds() compute_thresholds__n_firings = 1 -# in basic_mus_properties() +# In basic_mus_properties() basic_mus_properties__n_firings_rt_dert = 1 basic_mus_properties__accuracy = "default" basic_mus_properties__ignore_negative_ipts = False @@ -76,30 +77,45 @@ # ----------------------------------- tools ----------------------------------- -# in resize_emgfile() +# In resize_emgfile() resize_emgfile__how = "ref_signal" resize_emgfile__accuracy = "recalculate" resize_emgfile__ignore_negative_ipts = False +# In compute_svr() +compute_svr__gammain = 1/1.6 +compute_svr__regparam = 1/0.370 +compute_svr__endpointweights_numpulses = 5 +compute_svr__endpointweights_magnitude = 5 +compute_svr__discontfiring_dur = 1.0 + + +# ------------------------------------ pic ------------------------------------ +# In compute_deltaf() +compute_deltaf__recruitment_difference_cutoff = 1.0 +compute_deltaf__corr_cutoff = 0.7 +compute_deltaf__controlunitmodulation_cutoff = 0.5 + # ------------------------------------ muap ----------------------------------- -# in tracking() +# In tracking() tracking__firings = "all" tracking__derivation = "sd" tracking__gui = True tracking__gui_addrefsig = True tracking__gui_csv_separator = "\t" -# in remove_duplicates_between() +# In remove_duplicates_between() remove_duplicates_between__firings = "all" remove_duplicates_between__derivation = "sd" -# in MUcv_gui() +# In MUcv_gui() MUcv_gui__n_firings = [0, 50] MUcv_gui__muaps_timewindow = 50 MUcv_gui__figsize = [25, 20] MUcv_gui__csv_separator = "\t" + # --------------------------------- electrodes -------------------------------- # This custom sorting order is valid for all the GUI windows, although the # documentation is accessible in the api of the electrodes module. diff --git a/openhdemg/gui/gui_modules/advanced_analyses.py b/openhdemg/gui/gui_modules/advanced_analyses.py index 82ed6fb..12eea2f 100644 --- a/openhdemg/gui/gui_modules/advanced_analyses.py +++ b/openhdemg/gui/gui_modules/advanced_analyses.py @@ -101,6 +101,10 @@ def __init__(self, parent): self.filter_adv = BooleanVar() self.show_adv = BooleanVar() self.which_adv = StringVar() + self.average_method_adv = StringVar() + self.smoothing_method_adv = StringVar() + self.normalisation_adv = StringVar() + self.clean_adv = StringVar() # TODO this may become a bool self.emgfile1 = {} self.emgfile2 = {} self.extension_factor_adv = StringVar() @@ -109,27 +113,6 @@ def __init__(self, parent): self.parent = parent self.parent.load_settings() - # Disable config for DELSYS files - try: - if self.parent.resdict["SOURCE"] == "DELSYS": - show_error_dialog( - parent=self, - error=None, - solution=str( - "Advanced Tools for Delsys are only accessible from " - + "the library." - ), - ) - return - except AttributeError: - pass # Allow to open advanced tools even if no file is loaded - except Exception as e: - show_error_dialog( - parent=self, error=e, - solution=str("An unexpected error occurred."), - ) - return - # Open window self.a_window = ctk.CTkToplevel(fg_color="LightBlue4") self.a_window.title("Advanced Tools Window") @@ -170,14 +153,16 @@ def __init__(self, parent): "Motor Unit Tracking", "Duplicate Removal", "Conduction Velocity", + "Persistent Inward Currents", ) self.advanced_method = StringVar() adv_box = ctk.CTkComboBox( self.a_window, - width=200, + width=260, variable=self.advanced_method, values=adv_box_values, state="readonly", + command=self.enable_disable_a_window_elements ) adv_box.grid(row=2, column=1, sticky=(W, E)) self.advanced_method.set("Motor Unit Tracking") @@ -188,14 +173,14 @@ def __init__(self, parent): font=("Segoe UI", 18, "bold"), ).grid(row=3, column=0, sticky=(W, E)) self.mat_orientation_adv = StringVar() - orientation = ctk.CTkComboBox( + self.orientation_combobox = ctk.CTkComboBox( self.a_window, - width=100, + width=260, variable=self.mat_orientation_adv, values=("0", "180"), state="readonly", ) - orientation.grid(row=3, column=1, sticky=(W, E)) + self.orientation_combobox.grid(row=3, column=1, sticky=(W, E)) self.mat_orientation_adv.set("180") # Matrix code @@ -210,14 +195,14 @@ def __init__(self, parent): "GR04MM1305", "GR10MM0808", ) - matrix_code = ctk.CTkComboBox( + self.matrix_code_combobox = ctk.CTkComboBox( self.a_window, - width=150, + width=260, variable=self.mat_code_adv, values=matrix_code_vals, state="readonly", ) - matrix_code.grid(row=4, column=1, sticky=(W, E)) + self.matrix_code_combobox.grid(row=4, column=1, sticky=(W, E)) self.mat_code_adv.set("GR08MM1305") # Trace variable for updating window @@ -240,9 +225,38 @@ def __init__(self, parent): for child in self.a_window.winfo_children(): child.grid_configure(padx=5, pady=5) + # Check emgfile source and adjust available functionalities + try: + if self.parent.resdict["SOURCE"] == "DELSYS": + new_values = ("Persistent Inward Currents", ) + adv_box.configure(values=new_values) + self.advanced_method.set("Persistent Inward Currents") + self.enable_disable_a_window_elements(None) + except AttributeError: + pass # Allow to open advanced tools even if no file is loaded + except Exception as e: + show_error_dialog( + parent=self, error=e, + solution=str("An unexpected error occurred."), + ) + return + # ----------------------------------------------------------------------------------------------- # Advanced Analysis functionalities + def enable_disable_a_window_elements(self, _): + """ + Enable or disable comboboxes in a_window depending on the value of + self.advanced_method. + """ + + if self.advanced_method.get() == "Persistent Inward Currents": + self.orientation_combobox.configure(state="disabled") + self.matrix_code_combobox.configure(state="disabled") + else: + self.orientation_combobox.configure(state="readonly") + self.matrix_code_combobox.configure(state="readonly") + def on_matrix_none_adv(self, *args): """ Handle changes in the matrix code selection in the AdvancedAnalysis @@ -303,7 +317,8 @@ def advanced_analysis(self): depending on the advanced analysis option chosen by the user. It dynamically creates GUI elements like dropdowns, buttons, and checkboxes specific to the selected analysis tool, such as 'Motor Unit - Tracking', 'Conduction Velocity', or 'Duplicate Removal'. + Tracking', 'Conduction Velocity', 'Duplicate Removal' or 'Persistent + Inward Currents'. Raises ------ @@ -512,6 +527,9 @@ def advanced_analysis(self): self.head.destroy() self.a_window.destroy() + # Reload settings + self.parent.load_settings() + if self.mat_code_adv.get() == "None": # Get rows and columns and turn into list @@ -595,8 +613,121 @@ def advanced_analysis(self): ) self.head.destroy() + # TODO better to do specific windows rather than modifying a default + if self.advanced_method.get() == "Persistent Inward Currents": + # Create appropriate widgets + for widget in self.head.winfo_children(): + widget.destroy() + + r = 0 + # Smoothing Method label + smoothing_method_label = ctk.CTkLabel( + self.head, text="Smoothing Method:", + font=("Segoe UI", 18, "bold"), + ) + smoothing_method_label.grid(column=0, row=r) + # Combobox for Smoothing Method + smoothing_method_combobox = ctk.CTkComboBox( + self.head, + values=("Support Vector Regression",), + variable=self.smoothing_method_adv, + state="disabled", + width=260, + ) + smoothing_method_combobox.grid(column=1, row=r) + self.smoothing_method_adv.set("Support Vector Regression") + + r += 1 + # Average Method label + average_method_label = ctk.CTkLabel( + self.head, text="Average Method:", + font=("Segoe UI", 18, "bold"), + ) + average_method_label.grid(column=0, row=r) + # Combobox for Average Method + average_method_combobox = ctk.CTkComboBox( + self.head, + values=("test_unit_average", "all"), + variable=self.average_method_adv, + state="readonly", + width=260, + ) + average_method_combobox.grid(column=1, row=r) + self.average_method_adv.set("test_unit_average") + + r += 1 + # Normalisation label + normalisation_label = ctk.CTkLabel( + self.head, text="Normalisation:", + font=("Segoe UI", 18, "bold"), + ) + normalisation_label.grid(column=0, row=r) + # Combobox for normalisation_label + normalisation_combobox = ctk.CTkComboBox( + self.head, + values=("False", "ctrl_max_desc"), + variable=self.normalisation_adv, + state="readonly", + width=260, + ) + normalisation_combobox.grid(column=1, row=r) + self.normalisation_adv.set("False") + + r += 1 + # clean label # TODO this may become a bool + clean_label = ctk.CTkLabel( + self.head, text="Clean:", + font=("Segoe UI", 18, "bold"), + ) + clean_label.grid(column=0, row=r) + # Combobox for normalisation_label + clean_combobox = ctk.CTkComboBox( + self.head, + values=("True", "False"), + variable=self.clean_adv, + state="readonly", + width=260, + ) + clean_combobox.grid(column=1, row=r) + self.clean_adv.set("True") + + r += 1 + text = ( + "Sort MUs based on recruitment order before PICs estimation." + + "\nThis can be done in the main window." + ) + description_label = ctk.CTkLabel( + self.head, + text=text, + font=("Segoe UI", 10.5), + ) + description_label.grid(column=1, row=r, pady=(5, 0)) + + r += 1 + # Add button to execute MU tracking + compute_pic_button = ctk.CTkButton( + self.head, + text="Compute PIC", + command=self.compute_pic, + ) + compute_pic_button.grid( + column=0, row=14, columnspan=2, sticky=(W, E), + ) + + # Add padding + for child in self.head.winfo_children(): + child.grid_configure(padx=5, pady=5) + + # Configure weights + for row in range(r): + self.head.rowconfigure(row, weight=1) + for column in range(4): + self.head.columnconfigure(column, weight=1) + # Destroy first window to avoid too many pop-ups - self.a_window.withdraw() # TODO check for reference as destroy will evoke tcl error upon continuing, whereas withdraw creates issue when closing. + self.a_window.withdraw() + # TODO check for reference as destroy will evoke tcl error upon + # continuing, whereas withdraw creates issue when closing. ### Define function for advanced analysis tools def open_emgfile1(self): @@ -778,7 +909,7 @@ def track_mus(self): openhdemg.tracking() """ - # Reload settings for tracking + # Reload settings for tracking and duplicates removal self.parent.load_settings() try: @@ -823,7 +954,7 @@ def track_mus(self): # uninterrupted processes. if not self.parent.settings.tracking__gui: track_terminal = ttk.LabelFrame( - self.head, text="MUs Tracking Result", height=100, + self.head, text="MUs Tracking Results", height=100, relief="ridge", ) track_terminal.grid( @@ -958,3 +1089,81 @@ def remove_duplicates_between(self): + "\n - custom_sorting_order in settings" ), ) + + def compute_pic(self): + try: + # Reload settings + self.parent.load_settings() + + # Check that at least 1-2 MU is present, compute PICs and show results in the terminal + if self.parent.resdict["NUMBER_OF_MUS"] < 2: + show_error_dialog( + parent=self, + error=None, + solution=str( + "Not enough MUs for PICs estimation." + ), + ) + return + + # Get smoothed DR + print("\n-----------------------------\nDR smoothing in progress") + if self.smoothing_method_adv.get() == "Support Vector Regression": + svrfits = openhdemg.compute_svr( + emgfile=self.parent.resdict, + gammain=self.parent.settings.compute_svr__gammain, + regparam=self.parent.settings.compute_svr__regparam, + endpointweights_numpulses=self.parent.settings.compute_svr__endpointweights_numpulses, + endpointweights_magnitude=self.parent.settings.compute_svr__endpointweights_magnitude, + discontfiring_dur=self.parent.settings.compute_svr__discontfiring_dur, + ) + else: + pass # TODO to implement + + print( + "-----------------------------\nDeltaF estimation in progress" + ) + delta_f = openhdemg.compute_deltaf( + emgfile=self.parent.resdict, + smoothfits=svrfits["gensvr"], + average_method=self.average_method_adv.get(), + normalisation=self.normalisation_adv.get(), + recruitment_difference_cutoff=self.parent.settings.compute_deltaf__recruitment_difference_cutoff, + corr_cutoff=self.parent.settings.compute_deltaf__corr_cutoff, + controlunitmodulation_cutoff=self.parent.settings.compute_deltaf__controlunitmodulation_cutoff, + clean=self.clean_adv.get(), + ) + print( + "-----------------------------\nDeltaF estimation finished" + + "\n-----------------------------\n" + ) + + # Add results terminal + track_terminal = ttk.LabelFrame( + self.head, text="MUs PICs Results", height=100, width=100, + relief="ridge", + ) + track_terminal.grid( + column=2, + row=0, + columnspan=2, + rowspan=15, + pady=8, + padx=10, + sticky=(N, S, W, E), + ) + + # Add table containing results to the label frame + track_table = Table(track_terminal, dataframe=delta_f) + track_table.show() + + except AttributeError as e: + show_error_dialog( + parent=self, + error=e, + solution=str( + "Please make sure to load a file in the main window " + + "for Persistent Inward Currents estimation." + ), + ) + self.head.destroy() diff --git a/openhdemg/gui/settings.py b/openhdemg/gui/settings.py index b14f761..d70adff 100644 --- a/openhdemg/gui/settings.py +++ b/openhdemg/gui/settings.py @@ -19,30 +19,31 @@ import numpy as np + # --------------------------------- openfiles --------------------------------- -# in emg_from_demuse() +# In emg_from_demuse() emg_from_demuse__ignore_negative_ipts = False -# in emg_from_otb() +# In emg_from_otb() emg_from_otb__ext_factor = 8 emg_from_otb__refsig = [True, "fullsampled"] emg_from_otb__extras = None emg_from_otb__ignore_negative_ipts = False -# in refsig_from_otb() +# In refsig_from_otb() refsig_from_otb__refsig = "fullsampled" refsig_from_otb__extras = None -# in emg_from_delsys() +# In emg_from_delsys() emg_from_delsys__emg_sensor_name = "Galileo sensor" emg_from_delsys__refsig_sensor_name = "Trigno Load Cell" emg_from_delsys__filename_from = "mus_directory" -# in refsig_from_delsys() +# In refsig_from_delsys() refsig_from_delsys__refsig_sensor_name = "Trigno Load Cell" -# in emg_from_customcsv() +# In emg_from_customcsv() emg_from_customcsv__ref_signal = "REF_SIGNAL" emg_from_customcsv__raw_signal = "RAW_SIGNAL" emg_from_customcsv__ipts = "IPTS" @@ -53,21 +54,21 @@ emg_from_customcsv__fsamp = 2048 emg_from_customcsv__ied = 8 -# in refsig_from_customcsv() +# In refsig_from_customcsv() refsig_from_customcsv__ref_signal = "REF_SIGNAL" refsig_from_customcsv__extras = "EXTRAS" refsig_from_customcsv__fsamp = 2048 -# in save_json_emgfile() +# In save_json_emgfile() save_json_emgfile__compresslevel = 4 # ---------------------------------- analysis --------------------------------- -# in compute_thresholds() +# In compute_thresholds() compute_thresholds__n_firings = 1 -# in basic_mus_properties() +# In basic_mus_properties() basic_mus_properties__n_firings_rt_dert = 1 basic_mus_properties__accuracy = "default" basic_mus_properties__ignore_negative_ipts = False @@ -76,30 +77,45 @@ # ----------------------------------- tools ----------------------------------- -# in resize_emgfile() +# In resize_emgfile() resize_emgfile__how = "ref_signal" resize_emgfile__accuracy = "recalculate" resize_emgfile__ignore_negative_ipts = False +# In compute_svr() +compute_svr__gammain = 1/1.6 +compute_svr__regparam = 1/0.370 +compute_svr__endpointweights_numpulses = 5 +compute_svr__endpointweights_magnitude = 5 +compute_svr__discontfiring_dur = 1.0 + + +# ------------------------------------ pic ------------------------------------ +# In compute_deltaf() +compute_deltaf__recruitment_difference_cutoff = 1.0 +compute_deltaf__corr_cutoff = 0.7 +compute_deltaf__controlunitmodulation_cutoff = 0.5 + # ------------------------------------ muap ----------------------------------- -# in tracking() +# In tracking() tracking__firings = "all" tracking__derivation = "sd" tracking__gui = True tracking__gui_addrefsig = True tracking__gui_csv_separator = "\t" -# in remove_duplicates_between() +# In remove_duplicates_between() remove_duplicates_between__firings = "all" remove_duplicates_between__derivation = "sd" -# in MUcv_gui() +# In MUcv_gui() MUcv_gui__n_firings = [0, 50] MUcv_gui__muaps_timewindow = 50 MUcv_gui__figsize = [25, 20] MUcv_gui__csv_separator = "\t" + # --------------------------------- electrodes -------------------------------- # This custom sorting order is valid for all the GUI windows, although the # documentation is accessible in the api of the electrodes module. From 3435c7c93781034b6eba862eeba874aa2f3ade31 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 4 Jun 2024 15:47:32 +0200 Subject: [PATCH 16/35] Update remove_duplicates_between() --- openhdemg/gui/gui_modules/advanced_analyses.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openhdemg/gui/gui_modules/advanced_analyses.py b/openhdemg/gui/gui_modules/advanced_analyses.py index 12eea2f..71afc60 100644 --- a/openhdemg/gui/gui_modules/advanced_analyses.py +++ b/openhdemg/gui/gui_modules/advanced_analyses.py @@ -1051,6 +1051,7 @@ def remove_duplicates_between(self): custom_sorting_order=self.parent.settings.custom_sorting_order, filter=self.filter_adv.get(), show=self.show_adv.get(), + gui=False, which=self.which_adv.get(), ) From 1e46396960fd19b24b7851b4db4deabd3ba474c0 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 4 Jun 2024 15:48:10 +0200 Subject: [PATCH 17/35] Add gui in remove_duplicates_between() --- openhdemg/library/muap.py | 77 ++++++++++++++++++++++++++++++++++----- 1 file changed, 68 insertions(+), 9 deletions(-) diff --git a/openhdemg/library/muap.py b/openhdemg/library/muap.py index 0295350..d13febb 100644 --- a/openhdemg/library/muap.py +++ b/openhdemg/library/muap.py @@ -1431,8 +1431,10 @@ def __init__( # Plot the first MU pair self.gui_plot() - # Bring back the GUI in in the foreground + # Bring back the GUI in the foreground self.root.lift() + self.root.attributes('-topmost', True) + self.root.after_idle(self.root.attributes, '-topmost', False) # Start mainloop self.root.mainloop() @@ -1596,7 +1598,6 @@ def get_results(self): return self.tracking_res -# TODO add GUI? def remove_duplicates_between( emgfile1, emgfile2, @@ -1611,7 +1612,11 @@ def remove_duplicates_between( custom_sorting_order=None, custom_muaps=None, filter=True, + multiprocessing=True, show=False, + gui=True, + gui_addrefsig=True, + gui_csv_separator="\t", which="munumber", ): """ @@ -1695,8 +1700,20 @@ def remove_duplicates_between( filter : bool, default True If true, when the same MU has a match of XCC > threshold with multiple MUs, only the match with the highest XCC is returned. + multiprocessing : bool, default True + If True (default) parallel processing will be used to reduce execution + time. show : bool, default False Whether to plot the STA of pairs of MUs with XCC above threshold. + gui : bool, default True + If True (default) a GUI for the visual inspection and manual selection + of the tracking results will be called. + gui_addrefsig : bool, default True + If True, the REF_SIGNAL is plotted in front of the IDR with a + separated y-axes. This is used only when gui=True. + gui_csv_separator : str, default "\t" + The field delimiter used by the GUI to create the .csv copied to the + clipboard. This is used only when gui=True. which : str {"munumber", "accuracy"}, default "munumber" How to remove the duplicated MUs. @@ -1713,6 +1730,8 @@ def remove_duplicates_between( tracking_res : pd.DataFrame The results of the tracking including the MU from file 1, MU from file 2 and the normalised cross-correlation value (XCC). + If gui=True, an additional column indicating the inclusion/exclusion + of the MUs pairs will also be present in tracking_res. See also -------- @@ -1723,7 +1742,29 @@ def remove_duplicates_between( Examples -------- - Remove duplicated MUs between two OTB files and save the emgfiles + Remove duplicated MUs between two OPENHDEMG files and inspect the tracking + outcome via a convenient GUI. Then Save the emgfiles without duplicates. + Of the 2 duplicated MUs, the one with the lowest accuracy is removed. + + >>> import openhdemg.library as emg + >>> emgfile1 = emg.askopenfile(filesource="OPENHDEMG") + >>> emgfile2 = emg.askopenfile(filesource="OPENHDEMG") + >>> emgfile1, emgfile2, tracking_res = emg.remove_duplicates_between( + ... emgfile1, + ... emgfile2, + ... firings="all", + ... derivation="mono", + ... timewindow=50, + ... threshold=0.9, + ... matrixcode="GR08MM1305", + ... orientation=180, + ... gui=True, + ... which="accuracy", + ... ) + >>> emg.asksavefile(emgfile1) + >>> emg.asksavefile(emgfile2) + + Remove duplicated MUs between two OTB files and directly save the emgfiles without duplicates. The duplicates are removed from the file with more MUs. @@ -1743,16 +1784,18 @@ def remove_duplicates_between( ... n_cols=None, ... filter=True, ... show=False, + ... gui=False, ... which="munumber", ... ) >>> emg.asksavefile(emgfile1) >>> emg.asksavefile(emgfile2) Remove duplicated MUs between two files where channels are sorted with a - custom order and save the emgfiles without duplicates. Of the 2 duplicated - MUs, the one with the lowest accuracy is removed. + custom order and directly save the emgfiles without duplicates. Of the 2 + duplicated MUs, the one with the lowest accuracy is removed. >>> import openhdemg.library as emg + >>> import numpy as np >>> emgfile1 = emg.askopenfile(filesource="CUSTOMCSV") >>> emgfile2 = emg.askopenfile(filesource="CUSTOMCSV") >>> custom_sorting_order = [ @@ -1776,6 +1819,7 @@ def remove_duplicates_between( ... custom_sorting_order=custom_sorting_order, ... filter=True, ... show=False, + ... gui=False, ... which="accuracy", ... ) >>> emg.asksavefile(emgfile1) @@ -1802,14 +1846,27 @@ def remove_duplicates_between( custom_muaps=custom_muaps, exclude_belowthreshold=True, filter=filter, + multiprocessing=multiprocessing, show=show, + gui=gui, + gui_addrefsig=gui_addrefsig, + gui_csv_separator=gui_csv_separator, ) + # If the tracking gui has been used to include/exclude pairs, use the + # included pairs only. + if gui: + tracking_res_cleaned = tracking_res[ + tracking_res["Inclusion"] == "Included" + ] + else: + tracking_res_cleaned = tracking_res + # Identify how to remove MUs if which == "munumber": if emgfile1["NUMBER_OF_MUS"] >= emgfile2["NUMBER_OF_MUS"]: # Remove MUs from emgfile1 - mus_to_remove = list(tracking_res["MU_file1"]) + mus_to_remove = list(tracking_res_cleaned["MU_file1"]) emgfile1 = delete_mus( emgfile=emgfile1, munumber=mus_to_remove, if_single_mu="remove" ) @@ -1818,7 +1875,7 @@ def remove_duplicates_between( else: # Remove MUs from emgfile2 - mus_to_remove = list(tracking_res["MU_file2"]) + mus_to_remove = list(tracking_res_cleaned["MU_file2"]) emgfile2 = delete_mus( emgfile=emgfile2, munumber=mus_to_remove, if_single_mu="remove" @@ -1831,7 +1888,7 @@ def remove_duplicates_between( # on ACCURACY value. to_remove1 = [] to_remove2 = [] - for i, row in tracking_res.iterrows(): + for i, row in tracking_res_cleaned.iterrows(): acc1 = emgfile1["ACCURACY"].loc[int(row["MU_file1"])] acc2 = emgfile2["ACCURACY"].loc[int(row["MU_file2"])] @@ -2258,8 +2315,10 @@ def __init__( # this will move the GUI in the background ??. self.gui_plot() - # Bring back the GUI in in the foreground + # Bring back the GUI in the foreground self.root.lift() + self.root.attributes('-topmost', True) + self.root.after_idle(self.root.attributes, '-topmost', False) # Start the main loop self.root.mainloop() From f1d8bdd2e6f6ff242fef1ee8c1116b57a71308cb Mon Sep 17 00:00:00 2001 From: Paul Date: Sat, 8 Jun 2024 20:37:58 +0200 Subject: [PATCH 18/35] Fixed iconbitmap path adavanced_analysis.py --- openhdemg/gui/gui_modules/advanced_analyses.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openhdemg/gui/gui_modules/advanced_analyses.py b/openhdemg/gui/gui_modules/advanced_analyses.py index 82ed6fb..49db8cd 100644 --- a/openhdemg/gui/gui_modules/advanced_analyses.py +++ b/openhdemg/gui/gui_modules/advanced_analyses.py @@ -137,7 +137,8 @@ def __init__(self, parent): # Set window icon head_path = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) iconpath = head_path + "/gui_files/Icon_transp.ico" - self.a_window.iconbitmap(default=iconpath) + self.a_window.iconbitmap(iconpath) + if platform.startswith("win"): self.a_window.after(200, lambda: self.a_window.iconbitmap(iconpath)) From 758dc55612032fbee5c9e119ea1cd0f238f3bbfe Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 11 Jun 2024 10:31:58 +0200 Subject: [PATCH 19/35] Update setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index a5dff5f..fdbb904 100644 --- a/setup.py +++ b/setup.py @@ -19,7 +19,7 @@ "pandastable==0.13.1", "scipy==1.11.3", "seaborn==0.13.0", - "joblib==1.3.2", + "scikit-learn==1.5.0", # "pre-commit==3.7.0", ] From 61d59d2356c21e388d84ccb449db1730bfaeda0c Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 11 Jun 2024 10:32:18 +0200 Subject: [PATCH 20/35] Minor spelling fix --- openhdemg/library/mathtools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openhdemg/library/mathtools.py b/openhdemg/library/mathtools.py index de6c0cb..dac3022 100644 --- a/openhdemg/library/mathtools.py +++ b/openhdemg/library/mathtools.py @@ -18,7 +18,7 @@ def min_max_scaling(series_or_df): """ Min-max scaling of pd.series or pd.dataframes. - Min-max feature scaling is often simply referred to as normalization, + Min-max feature scaling is often simply referred to as normalisation, which rescales the dataset feature to a range of 0 - 1. It's calculated by subtracting the feature's minimum value from the value and then dividing it by the difference between the maximum and minimum From f5966438ae1ee312d4c09db5d0cc55922e103e82 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 11 Jun 2024 10:37:56 +0200 Subject: [PATCH 21/35] Fix remasining iconbitmap on MAC --- openhdemg/gui/gui_modules/advanced_analyses.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openhdemg/gui/gui_modules/advanced_analyses.py b/openhdemg/gui/gui_modules/advanced_analyses.py index de8417b..2dab858 100644 --- a/openhdemg/gui/gui_modules/advanced_analyses.py +++ b/openhdemg/gui/gui_modules/advanced_analyses.py @@ -343,7 +343,7 @@ def advanced_analysis(self): # Set window icon head_path = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) iconpath = head_path + "/gui_files/Icon_transp.ico" - self.head.iconbitmap(default=iconpath) + self.head.iconbitmap(iconpath) if platform.startswith("win"): self.head.after(200, lambda: self.head.iconbitmap(iconpath)) From 31a54333c4aeff5319b6b5f828eb3413ade2d31b Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 11 Jun 2024 14:36:44 +0200 Subject: [PATCH 22/35] Added idr_range in analysis functions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The functions “compute_dr”, “compute_drvariability”, “compute_covisi” and “basic_mus_properties” now allow to exclude firings outside a custom range of frequencies. --- openhdemg/gui/backup_settings.py | 4 + openhdemg/gui/gui_modules/mu_properties.py | 11 +++ openhdemg/gui/settings.py | 4 + openhdemg/library/analysis.py | 103 +++++++++++++++++++++ 4 files changed, 122 insertions(+) diff --git a/openhdemg/gui/backup_settings.py b/openhdemg/gui/backup_settings.py index d70adff..b477037 100644 --- a/openhdemg/gui/backup_settings.py +++ b/openhdemg/gui/backup_settings.py @@ -68,11 +68,15 @@ # In compute_thresholds() compute_thresholds__n_firings = 1 +# In compute_dr() +compute_dr__idr_range = None + # In basic_mus_properties() basic_mus_properties__n_firings_rt_dert = 1 basic_mus_properties__accuracy = "default" basic_mus_properties__ignore_negative_ipts = False basic_mus_properties__constrain_pulses = [True, 3] +basic_mus_properties__idr_range = None # ----------------------------------- tools ----------------------------------- diff --git a/openhdemg/gui/gui_modules/mu_properties.py b/openhdemg/gui/gui_modules/mu_properties.py index ea75c24..7022c63 100644 --- a/openhdemg/gui/gui_modules/mu_properties.py +++ b/openhdemg/gui/gui_modules/mu_properties.py @@ -318,6 +318,9 @@ def compute_mu_threshold(self): compute_thresholds in library. """ + # Update settings + self.parent.load_settings() + try: # Compute thresholds self.parent.mu_thresholds = openhdemg.compute_thresholds( @@ -374,6 +377,9 @@ def compute_mu_dr(self): compute_dr in library. """ + # Update settings + self.parent.load_settings() + try: # Compute discharge rates self.parent.mus_dr = openhdemg.compute_dr( @@ -381,6 +387,7 @@ def compute_mu_dr(self): n_firings_RecDerec=int(self.firings_rec.get()), n_firings_steady=int(self.firings_ste.get()), event_=self.dr_event.get(), + idr_range=self.parent.settings.compute_dr__idr_range, ) # Display results self.parent.display_results(self.parent.mus_dr) @@ -433,6 +440,9 @@ def basic_mus_properties(self): basic_mus_properties in library. """ + # Update settings + self.parent.load_settings() + try: # Calculate properties self.parent.mu_prop_df = openhdemg.basic_mus_properties( @@ -440,6 +450,7 @@ def basic_mus_properties(self): n_firings_rt_dert=self.parent.settings.basic_mus_properties__n_firings_rt_dert, n_firings_RecDerec=int(self.b_firings_rec.get()), n_firings_steady=int(self.b_firings_ste.get()), + idr_range=self.parent.settings.basic_mus_properties__idr_range, accuracy=self.parent.settings.basic_mus_properties__accuracy, ignore_negative_ipts=self.parent.settings.basic_mus_properties__ignore_negative_ipts, constrain_pulses=self.parent.settings.basic_mus_properties__constrain_pulses, diff --git a/openhdemg/gui/settings.py b/openhdemg/gui/settings.py index d70adff..b477037 100644 --- a/openhdemg/gui/settings.py +++ b/openhdemg/gui/settings.py @@ -68,11 +68,15 @@ # In compute_thresholds() compute_thresholds__n_firings = 1 +# In compute_dr() +compute_dr__idr_range = None + # In basic_mus_properties() basic_mus_properties__n_firings_rt_dert = 1 basic_mus_properties__accuracy = "default" basic_mus_properties__ignore_negative_ipts = False basic_mus_properties__constrain_pulses = [True, 3] +basic_mus_properties__idr_range = None # ----------------------------------- tools ----------------------------------- diff --git a/openhdemg/library/analysis.py b/openhdemg/library/analysis.py index d8f297f..d2159b5 100644 --- a/openhdemg/library/analysis.py +++ b/openhdemg/library/analysis.py @@ -196,6 +196,7 @@ def compute_dr( start_steady=-1, end_steady=-1, event_="rec_derec_steady", + idr_range=None, ): """ Calculate the discharge rate (DR). @@ -233,6 +234,12 @@ def compute_dr( ``steady`` DR is calculated during the steady-state phase. + idr_range : None or list, default None + If idr_range is a list [lower_limit, upper_limit], only firings with an + instantaneous discharge rate (IDR) within the limits are used for DR + calculation. lower_limit and upper_limit should be in pulses per + second. See examples section. + If idr_range is None, all the firings are used for DR calculation. Returns ------- @@ -300,6 +307,21 @@ def compute_dr( 1 14.440561 10.019572 11.822081 11.683134 2 7.293547 5.846093 7.589531 8.055731 3 13.289651 9.694317 11.613640 11.109796 + + The firings used for DR calculation can be filtered to avoid too closed + firings (e.g., doublets) or intervals of inactivity. + + >>> emgfile = emg.emg_from_samplefile() + >>> idr = emg.compute_dr( + ... emgfile, start_steady=15000, end_steady=51000, idr_range=[1, 50], + ... ) + >>> idr + DR_rec DR_derec ... DR_end_steady DR_all_steady DR_all + 0 3.341579 4.606835 ... 8.506270 7.895837 7.657269 + 1 5.701081 4.662196 ... 6.271989 6.916566 6.814687 + 2 5.699017 3.691367 ... 7.221309 8.183408 7.949294 + 3 7.548770 5.449581 ... 10.399834 11.164994 10.693076 + 4 8.344515 5.333535 ... 9.694317 10.750855 10.543011 """ # Check that all the inputs are correct @@ -324,6 +346,23 @@ def compute_dr( idr = compute_idr(emgfile=emgfile) + # Filter firings outside the idr_range, if required + if idr_range is not None: + if not isinstance(idr_range, list): + raise ValueError( + "idr_range can be None or a list of 2 numbers." + + f"A{type(idr_range)} was passed instead." + ) + else: + if len(idr_range) != 2: + raise ValueError( + "idr_range can be None or a list of 2 numbers." + + f"The list contains {len(idr_range)} numbers instead." + ) + for mu in idr.keys(): + idr[mu]["idr"] = idr[mu]["idr"][idr[mu]["idr"] > idr_range[0]] + idr[mu]["idr"] = idr[mu]["idr"][idr[mu]["idr"] < idr_range[1]] + # Check if we need to manually select the area for the steady-state phase title = ( "Select the start/end area of the steady-state by hovering the mouse" + @@ -457,6 +496,7 @@ def basic_mus_properties( n_firings_steady=10, start_steady=-1, end_steady=-1, + idr_range=None, accuracy="default", ignore_negative_ipts=False, constrain_pulses=[True, 3], @@ -496,6 +536,13 @@ def basic_mus_properties( The start and end point (in samples) of the steady-state phase. If < 0 (default), the user will need to manually select the start and end of the steady-state phase. + idr_range : None or list, default None + If idr_range is a list [lower_limit, upper_limit], only firings with an + instantaneous discharge rate (IDR) within the limits are used for DR + and coefficient of variation of interspike interval (COVisi) + calculation. lower_limit and upper_limit should be in pulses per + second. See compute_dr() examples section. + If idr_range is None, all the firings are used for DR calculation. accuracy : str {"default", "SIL", "PNR", "SIL_PNR"}, default "default" The accuracy measure to return. @@ -732,6 +779,7 @@ def basic_mus_properties( n_firings_steady=n_firings_steady, start_steady=start_steady, end_steady=end_steady, + idr_range=idr_range, ) exportable_df = pd.concat([exportable_df, mus_dr], axis=1) @@ -742,6 +790,7 @@ def basic_mus_properties( start_steady=start_steady, end_steady=end_steady, event_="steady", + idr_range=idr_range, ) exportable_df = pd.concat([exportable_df, covisi], axis=1) @@ -763,6 +812,7 @@ def compute_covisi( start_steady=-1, end_steady=-1, event_="rec_derec_steady", + idr_range=None, single_mu_number=-1, ): """ @@ -800,6 +850,12 @@ def compute_covisi( ``steady`` covisi is calculated during the steady-state phase. + idr_range : None or list, default None + If idr_range is a list [lower_limit, upper_limit], only firings with an + instantaneous discharge rate (IDR) within the limits are used for + COVisi calculation. lower_limit and upper_limit should be in pulses per + second. See compute_dr() examples section. + If idr_range is None, all the firings are used for DR calculation. single_mu_number : int, default -1 Number of the specific MU to compute the COVisi. If single_mu_number >= 0, only the COVisi of the entire contraction @@ -884,6 +940,29 @@ def compute_covisi( # We use the idr pd.DataFrame to calculate the COVisi idr = compute_idr(emgfile=emgfile) + # Filter firings outside the idr_range, if required + if idr_range is not None: + if not isinstance(idr_range, list): + raise ValueError( + "idr_range can be None or a list of 2 numbers." + + f"A{type(idr_range)} was passed instead." + ) + else: + if len(idr_range) != 2: + raise ValueError( + "idr_range can be None or a list of 2 numbers." + + f"The list contains {len(idr_range)} numbers instead." + ) + idr_range[0] = emgfile["FSAMP"] / idr_range[0] + idr_range[1] = emgfile["FSAMP"] / idr_range[1] + for mu in idr.keys(): + idr[mu]["diff_mupulses"] = idr[mu]["diff_mupulses"][ + idr[mu]["diff_mupulses"] < idr_range[0] + ] + idr[mu]["diff_mupulses"] = idr[mu]["diff_mupulses"][ + idr[mu]["diff_mupulses"] > idr_range[1] + ] + # Check if we need to analyse all the MUs or a single MU if single_mu_number < 0: # Check if we need to manually select the area for the steady-state @@ -981,6 +1060,7 @@ def compute_drvariability( start_steady=-1, end_steady=-1, event_="rec_derec_steady", + idr_range=None, ): """ Calculate the DR variability. @@ -1018,6 +1098,12 @@ def compute_drvariability( ``steady`` DR variability is calculated during the steady-state phase. + idr_range : None or list, default None + If idr_range is a list [lower_limit, upper_limit], only firings with an + instantaneous discharge rate (IDR) within the limits are used for DR + variability calculation. lower_limit and upper_limit should be in + pulses per second. See compute_dr() examples section. + If idr_range is None, all the firings are used for DR calculation. Returns ------- @@ -1090,6 +1176,23 @@ def compute_drvariability( # We use the idr pd.DataFrame to calculate the COVisi idr = compute_idr(emgfile=emgfile) + # Filter firings outside the idr_range, if required + if idr_range is not None: + if not isinstance(idr_range, list): + raise ValueError( + "idr_range can be None or a list of 2 numbers." + + f"A{type(idr_range)} was passed instead." + ) + else: + if len(idr_range) != 2: + raise ValueError( + "idr_range can be None or a list of 2 numbers." + + f"The list contains {len(idr_range)} numbers instead." + ) + for mu in idr.keys(): + idr[mu]["idr"] = idr[mu]["idr"][idr[mu]["idr"] > idr_range[0]] + idr[mu]["idr"] = idr[mu]["idr"][idr[mu]["idr"] < idr_range[1]] + # Check if we need to manually select the area for the steady-state phase if event_ == "rec_derec_steady" or event_ == "steady": title = ( From cce8e5b9c32ce22ffdd8d4b76810796bf089b9bf Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Wed, 12 Jun 2024 16:31:11 +0200 Subject: [PATCH 23/35] Add plot_smoothed_dr() --- docs/api_plotemg.md | 7 + openhdemg/library/analysis.py | 12 +- openhdemg/library/plotemg.py | 243 ++++++++++++++++++++++++++++++++++ openhdemg/library/tools.py | 21 ++- 4 files changed, 265 insertions(+), 18 deletions(-) diff --git a/docs/api_plotemg.md b/docs/api_plotemg.md index dc815d4..8da6229 100644 --- a/docs/api_plotemg.md +++ b/docs/api_plotemg.md @@ -53,6 +53,13 @@ This module contains all the functions used to visualise the content of the impo
+::: openhdemg.library.plotemg.plot_smoothed_dr + options: + show_root_full_path: False + show_root_heading: True + +
+ ::: openhdemg.library.plotemg.plot_muaps options: show_root_full_path: False diff --git a/openhdemg/library/analysis.py b/openhdemg/library/analysis.py index d2159b5..25c0ae4 100644 --- a/openhdemg/library/analysis.py +++ b/openhdemg/library/analysis.py @@ -350,13 +350,13 @@ def compute_dr( if idr_range is not None: if not isinstance(idr_range, list): raise ValueError( - "idr_range can be None or a list of 2 numbers." + + "idr_range can be None or a list of 2 numbers. " + f"A{type(idr_range)} was passed instead." ) else: if len(idr_range) != 2: raise ValueError( - "idr_range can be None or a list of 2 numbers." + + "idr_range can be None or a list of 2 numbers. " + f"The list contains {len(idr_range)} numbers instead." ) for mu in idr.keys(): @@ -944,13 +944,13 @@ def compute_covisi( if idr_range is not None: if not isinstance(idr_range, list): raise ValueError( - "idr_range can be None or a list of 2 numbers." + + "idr_range can be None or a list of 2 numbers. " + f"A{type(idr_range)} was passed instead." ) else: if len(idr_range) != 2: raise ValueError( - "idr_range can be None or a list of 2 numbers." + + "idr_range can be None or a list of 2 numbers. " + f"The list contains {len(idr_range)} numbers instead." ) idr_range[0] = emgfile["FSAMP"] / idr_range[0] @@ -1180,13 +1180,13 @@ def compute_drvariability( if idr_range is not None: if not isinstance(idr_range, list): raise ValueError( - "idr_range can be None or a list of 2 numbers." + + "idr_range can be None or a list of 2 numbers. " + f"A{type(idr_range)} was passed instead." ) else: if len(idr_range) != 2: raise ValueError( - "idr_range can be None or a list of 2 numbers." + + "idr_range can be None or a list of 2 numbers. " + f"The list contains {len(idr_range)} numbers instead." ) for mu in idr.keys(): diff --git a/openhdemg/library/plotemg.py b/openhdemg/library/plotemg.py index ac8a357..fae3371 100644 --- a/openhdemg/library/plotemg.py +++ b/openhdemg/library/plotemg.py @@ -6,6 +6,7 @@ import matplotlib.pyplot as plt import seaborn as sns import pandas as pd +import numpy as np from openhdemg.library.tools import compute_idr from openhdemg.library.mathtools import min_max_scaling @@ -884,6 +885,8 @@ def plot_idr( ax1.set_xlabel("Time (Sec)" if timeinseconds else "Samples") elif len(munumber) == 1: + # TODO remove this part and check len before. + # See plot_smoothed_dr as a reference ax1 = sns.scatterplot( x=idr[munumber[0]]["timesec" if timeinseconds else "mupulses"], y=idr[munumber[0]]["idr"], @@ -917,6 +920,246 @@ def plot_idr( return fig +def plot_smoothed_dr( + emgfile, + smoothfits, + munumber="all", + addidr=True, + stack=True, + addrefsig=True, + timeinseconds=True, + figsize=[20, 15], + showimmediately=True, + tight_layout=True, +): + """ + Plot the instantaneous discharge rate (IDR). + + Parameters + ---------- + emgfile : dict + The dictionary containing the emgfile. + smoothfits : pd.DataFrame + Smoothed discharge rate estimates aligned in time. Columns should + contain the smoothed discharge rate for each MU. + munumber : str, int or list, default "all" + + ``all`` + IDR of all the MUs is plotted. + + Otherwise, a single MU (int) or multiple MUs (list of int) can be + specified. + The list can be passed as a manually-written list or with: + munumber=[*range(0, 12)]. + We need the " * " operator to unpack the results of range into a list. + munumber is expected to be with base 0 (i.e., the first MU in the file + is the number 0). + addidr : bool, default True + Whether to show also the IDR behind the smoothed DR line. + stack : bool, default True + Whether to stack multiple MUs. If False, all the MUs smoothed DR will + be plotted on the same line. + addrefsig : bool, default True + If True, the REF_SIGNAL is plotted in front of the MUs IDR with a + separated y-axes. + timeinseconds : bool, default True + Whether to show the time on the x-axes in seconds (True) + or in samples (False). + figsize : list, default [20, 15] + Size of the figure in centimeters [width, height]. + showimmediately : bool, default True + If True (default), plt.show() is called and the figure showed to the + user. + It is useful to set it to False when calling the function from the GUI. + tight_layout : bool, default True + If True (default), the plt.tight_layout() is called and the figure's + layout is improved. + It is useful to set it to False when calling the function from the GUI. + + Returns + ------- + fig : pyplot `~.figure.Figure` + + See also + -------- + - compute_svr : Fit MU discharge rates with Support Vector Regression, + nonlinear regression. + + Notes + ----- + munumber = "all" corresponds to: + munumber = [*range(0, emgfile["NUMBER_OF_MUS"])] + + Examples + -------- + Smooth MUs DR using Support Vector Regression. + + >>> import openhdemg.library as emg + >>> import pandas as pd + >>> emgfile = emg.emg_from_samplefile() + >>> emgfile = emg.sort_mus(emgfile=emgfile) + >>> svrfits = emg.compute_svr(emgfile) + + Plot the stacked smoothed DR of all the MUs and overlay the IDR and the + reference signal. + + >>> smoothfits = pd.DataFrame(svrfits["gensvr"]).transpose() + >>> emg.plot_smoothed_dr( + >>> emgfile, + >>> smoothfits=smoothfits, + >>> munumber="all", + >>> addidr=True, + >>> stack=True, + >>> addrefsig=True, + >>> ) + """ + + # Compute the IDR + idr = compute_idr(emgfile=emgfile) + + if addrefsig: + if not isinstance(emgfile["REF_SIGNAL"], pd.DataFrame): + raise TypeError( + "REF_SIGNAL is probably absent or it is not contained in a dataframe" + ) + + # Check if all the MUs have to be plotted + if isinstance(munumber, list): + if len(munumber) == 1: # Manage exception of single MU + munumber = munumber[0] + elif isinstance(munumber, str): + if emgfile["NUMBER_OF_MUS"] == 1: # Manage exception of single MU + munumber = 0 + else: + munumber = [*range(0, emgfile["NUMBER_OF_MUS"])] + + # Check smoothfits type + if not isinstance(smoothfits, pd.DataFrame): + raise TypeError("smoothfits must be a pd.DataFrame") + + # Use the subplot function to allow for the use of twinx() + fig, ax1 = plt.subplots( + figsize=(figsize[0] / 2.54, figsize[1] / 2.54), + num="Smoothed DR", + ) + + # Check if we have a single MU or a list of MUs to plot. + if isinstance(munumber, int): + # Plot IDR + if addidr: + sns.scatterplot( + x=idr[munumber]["timesec" if timeinseconds else "mupulses"], + y=idr[munumber]["idr"], + ax=ax1, + alpha=0.4, + ) # sns.scatterplot breaks if x or y are nan? + + # Plot smoothed DR + if timeinseconds: + x_smooth = np.arange(len(smoothfits[munumber])) / emgfile["FSAMP"] + else: + x_smooth = np.arange(len(smoothfits[munumber])) + ax1.plot(x_smooth, smoothfits[munumber], linewidth=2) + + # Labels + ax1.set_ylabel( + "MU {} (PPS)".format(munumber) + ) # Useful because if the MU is empty it won't show the channel number + ax1.set_xlabel("Time (Sec)" if timeinseconds else "Samples") + + elif isinstance(munumber, list) and stack: + for count, thisMU in enumerate(munumber): + # Plot IDR + if addidr: + # Normalise the IDR series and add 1 to the previous MUs to + # avoid overlapping of the MUs. + norm_idr = min_max_scaling(idr[thisMU]["idr"]) + count + sns.scatterplot( + x=idr[thisMU]["timesec" if timeinseconds else "mupulses"], + y=norm_idr, + ax=ax1, + alpha=0.4, + ) + + # Plot smoothed DR + if timeinseconds: + x_smooth = np.arange( + len(smoothfits[thisMU]) + ) / emgfile["FSAMP"] + else: + x_smooth = np.arange(len(smoothfits[thisMU])) + + if addidr: + # Scale factor to match smoothed DR and IDR when stacked + delta = ( + (smoothfits[thisMU].max() - smoothfits[thisMU].min()) / + (idr[thisMU]["idr"].max() - idr[thisMU]["idr"].min()) + ) + delta_min = ( + (smoothfits[thisMU].min() - idr[thisMU]["idr"].min()) / + (idr[thisMU]["idr"].max() - idr[thisMU]["idr"].min()) + ) + norm_fit = min_max_scaling(smoothfits[thisMU]) + norm_fit = norm_fit * delta + delta_min + count + ax1.plot(x_smooth, norm_fit, linewidth=2) + else: + norm_fit = min_max_scaling(smoothfits[thisMU]) + count + ax1.plot(x_smooth, norm_fit, linewidth=2) + + # Ensure correct and complete ticks on the left y axis + ax1.set_yticks([*range(len(munumber))]) + ax1.set_yticklabels([str(mu) for mu in munumber]) + # Set axes labels + ax1.set_ylabel("Motor units") + ax1.set_xlabel("Time (Sec)" if timeinseconds else "Samples") + + elif isinstance(munumber, list) and not stack: + for count, thisMU in enumerate(munumber): + # Plot IDR + if addidr: + sns.scatterplot( + x=idr[thisMU]["timesec" if timeinseconds else "mupulses"], + y=idr[thisMU]["idr"], + ax=ax1, + alpha=0.4, + ) + + # Plot smoothed DR + if timeinseconds: + x_smooth = np.arange( + len(smoothfits[thisMU]) + ) / emgfile["FSAMP"] + else: + x_smooth = np.arange(len(smoothfits[thisMU])) + ax1.plot(x_smooth, smoothfits[thisMU], linewidth=2) + + # Set axes labels + ax1.set_ylabel("Discharge rate (PPS)") + ax1.set_xlabel("Time (Sec)" if timeinseconds else "Samples") + + else: + raise TypeError( + "While calling the plot_idr function, you should pass an integer, a list or 'all' to munumber" + ) + + if addrefsig: + ax2 = ax1.twinx() + # Plot the ref signal + xref = ( + emgfile["REF_SIGNAL"].index / emgfile["FSAMP"] + if timeinseconds + else emgfile["REF_SIGNAL"].index + ) + sns.lineplot(y=emgfile["REF_SIGNAL"][0], x=xref, color="0.4", ax=ax2) + ax2.set_ylabel("MVC") + + showgoodlayout(tight_layout, despined="2yaxes" if addrefsig else False) + if showimmediately: + plt.show() + + return fig + + def plot_muaps( sta_dict, title="MUAPs from STA", diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index 5ae2312..193ac9a 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -1267,18 +1267,15 @@ def compute_svr( Quick plot showing the results. - >>> scl = 12 # offset MUs for viz - >>> idr = emg.compute_idr(emgfile) - >>> for ii in range(len(svrfits["svrfit"])): - ... xtmp = np.transpose([idr[ii].timesec[1:]]) - ... ytmp = idr[ii].idr[1:].to_numpy() - ... plt.scatter(xtmp, ytmp+scl*(ii), color='darkorange') - ... plt.plot( - ... svrfits["svrtime"][ii], - ... svrfits["svrfit"][ii]+scl*(ii), - ... color='cornflowerblue', - ... ) - >>> plt.show() + >>> smoothfits = pd.DataFrame(svrfits["gensvr"]).transpose() + >>> emg.plot_smoothed_dr( + >>> emgfile, + >>> smoothfits=smoothfits, + >>> munumber="all", + >>> addidr=False, + >>> stack=True, + >>> addrefsig=True, + >>> ) """ # TODO input checking and edge cases From 08987f4d9fc0fb3af879c0214822324802a56ad2 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Wed, 12 Jun 2024 17:18:00 +0200 Subject: [PATCH 24/35] Docstring error --- openhdemg/library/plotemg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openhdemg/library/plotemg.py b/openhdemg/library/plotemg.py index fae3371..fedc1b9 100644 --- a/openhdemg/library/plotemg.py +++ b/openhdemg/library/plotemg.py @@ -933,7 +933,7 @@ def plot_smoothed_dr( tight_layout=True, ): """ - Plot the instantaneous discharge rate (IDR). + Plot the smoothed discharge rate. Parameters ---------- From 79ba14c843bc1d12fff281a2566519c97ff74a43 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Wed, 12 Jun 2024 17:18:32 +0200 Subject: [PATCH 25/35] Update test_analysis for idr_range --- openhdemg/tests/unit/test_anlysis.py | 70 ++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/openhdemg/tests/unit/test_anlysis.py b/openhdemg/tests/unit/test_anlysis.py index 38f9abd..4c7bd62 100644 --- a/openhdemg/tests/unit/test_anlysis.py +++ b/openhdemg/tests/unit/test_anlysis.py @@ -202,6 +202,22 @@ def test_compute_dr(self): res["DR_end_steady"][1], 7.062, places=2, ) + # Change idr_range + res = compute_dr( + emgfile=emgfile, + n_firings_RecDerec=4, + n_firings_steady=10, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + idr_range=[7, 10], + ) + + self.assertTrue(np.isnan(res["DR_end_steady"][0])) + self.assertAlmostEqual( + res["DR_all"][1], 7.644, places=2, + ) + def test_basic_mus_properties(self): """ Test the basic_mus_properties function with the samplefile. @@ -287,6 +303,29 @@ def test_basic_mus_properties(self): res["avg_PNR"][0], 29.895, places=2, ) + # Change idr_range + res = basic_mus_properties( + emgfile=emgfile, + n_firings_rt_dert=1, + n_firings_RecDerec=4, + n_firings_steady=10, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + idr_range=[7, 10], + accuracy="default", + ignore_negative_ipts=False, + constrain_pulses=[True, 3], + mvc=1234, + ) + + self.assertTrue(np.isnan(res["DR_end_steady"][0])) + self.assertAlmostEqual( + res["DR_all"][1], 7.644, places=2, + ) + self.assertAlmostEqual( + res["COVisi_steady"][2], 7.677, places=2, + ) + def test_compute_covisi(self): """ Test the compute_dr function with the samplefile. @@ -348,6 +387,22 @@ def test_compute_covisi(self): res["COVisi_all"][0], 19.104, places=2, ) + # Change idr_range + res = compute_covisi( + emgfile=emgfile, + n_firings_RecDerec=4, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + idr_range=[7, 10], + single_mu_number=-1, + ) + + self.assertTrue(np.isnan(res["COVisi_rec"][0])) + self.assertAlmostEqual( + res["COVisi_steady"][2], 7.677, places=2, + ) + def test_compute_drvariability(self): """ Test the compute_drvariability function with the samplefile. @@ -393,6 +448,21 @@ def test_compute_drvariability(self): self.assertTrue(np.isnan(res["DRvar_rec"][0])) self.assertTrue(np.isnan(res["DRvar_derec"][1])) + # Change idr_range + res = compute_drvariability( + emgfile=emgfile, + n_firings_RecDerec=4, + start_steady=0 + t_ramps, + end_steady=emgfile["EMG_LENGTH"] - t_ramps, + event_="rec_derec_steady", + idr_range=[7, 10], + ) + + self.assertTrue(np.isnan(res["DRvar_rec"][0])) + self.assertAlmostEqual( + res["DRvar_all"][1], 6.466, places=2, + ) + if __name__ == '__main__': unittest.main() From 4642e71d97ff4ea0abebeac6cd33d903ae82f8ba Mon Sep 17 00:00:00 2001 From: JABeauchamp <103655220+JABeauchamp@users.noreply.github.com> Date: Wed, 12 Jun 2024 12:08:33 -0500 Subject: [PATCH 26/35] Adding PIC unit test, fixing todos --- openhdemg/library/pic.py | 31 +++++----- openhdemg/library/tools.py | 4 -- openhdemg/tests/unit/test_pic.py | 99 ++++++++++++++++++++++++++++++++ 3 files changed, 114 insertions(+), 20 deletions(-) create mode 100644 openhdemg/tests/unit/test_pic.py diff --git a/openhdemg/library/pic.py b/openhdemg/library/pic.py index 7f22db8..c03296b 100644 --- a/openhdemg/library/pic.py +++ b/openhdemg/library/pic.py @@ -18,7 +18,7 @@ def compute_deltaf( recruitment_difference_cutoff=1.0, corr_cutoff=0.7, controlunitmodulation_cutoff=0.5, - clean="True", # TODO should this be a boolean True? + clean=True, ): """ Quantify delta F via paired motor unit analysis. @@ -37,16 +37,17 @@ def compute_deltaf( smoothfits : list of arrays Smoothed discharge rate estimates. Each array: motor unit discharge rate x samples aligned in time; - instances of non-firing = NaN, please #TODO maybe something is missing after please? + instances of non-firing = NaN Your choice of smoothing. See compute_svr gen_svr for example. average_method : str {"test_unit_average", "all"}, default "test_unit_average" - The method for test MU deltaF value. More to be added. TODO # TODOs should be moved outside documentation or they will be displayed in the website + The method for test MU deltaF value. More to be added. ``test_unit_average`` The average across all possible control units. - ``all`` TODO the user is not aware of what other average methods can be used - abc + ``all`` + This returns all possible MU pairs + normalisation : str {"False", "ctrl_max_desc"}, default "False" The method for deltaF nomalization. @@ -54,6 +55,7 @@ def compute_deltaf( Whether to normalise deltaF values to control unit descending range during test unit firing. See Skarabot et. al., 2023: https://www.biorxiv.org/content/10.1101/2023.10.16.562612v1 + recruitment_difference_cutoff : float, default 1 An exlusion criteria corresponding to the necessary difference between control and test MU recruitement in seconds. @@ -63,27 +65,24 @@ def compute_deltaf( controlunitmodulation_cutoff : float, default 0.5 An exclusion criteria corresponding to the necessary modulation of control unit discharge rate during test unit firing in Hz. - clean : str, default "True" + clean : bool, default True To remove values that do not meet exclusion criteria Returns ------- delta_f : pd.DataFrame A pd.DataFrame containing deltaF values and corresponding MU number. - TODO here we should explain that the resulting df will be different - depending on average_method. In particular, if average_method="all", - delta_f[MU][row] will contain a tuple representing... instead of a MU - number. + The resulting df will be different depending on average_method. + In particular, if average_method="all", delta_f[MU][row] will + contain a tuple representing the indices of the two motor units + for each given pair (reporter, test) and their corresponding + deltaF value. See also -------- - compute_svr : fit MU discharge rates with Support Vector Regression, nonlinear regression. - Warnings - -------- - TODO consider removing section from docstring while not implemented - Examples -------- Quantify delta F using svr fits. @@ -130,7 +129,7 @@ def compute_deltaf( # If less than 2 MUs, can not quantify deltaF if emgfile["NUMBER_OF_MUS"] < 2: dfret_ret = np.nan - mucombo_ret = np.nan*np.ones(1, 2) # TODO this line is not working, please replace with nan or what appropriate + mucombo_ret = np.nan*np.ones([1, 2]) delta_f = pd.DataFrame({'MU': mucombo_ret, 'dF': dfret_ret}) @@ -296,7 +295,7 @@ def compute_deltaf( controlmu.append(mucombo[-1][controlU-1]) testmu.append(mucombo[-1][1-controlU//2]) - if clean == "True": # Remove values that dont meet exclusion criteria + if clean: # Remove values that dont meet exclusion criteria rcrt_diff_bin = rcrt_diff > recruitment_difference_cutoff corr_bin = r_ret > corr_cutoff ctrl_mod_bin = ctrl_mod > controlunitmodulation_cutoff diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index 5ae2312..b93b885 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -1250,10 +1250,6 @@ def compute_svr( -------- - compute_deltaf : quantify delta F via paired motor unit analysis. - Warnings - -------- - TODO consider removing section from docstring while not implemented - Examples -------- Quantify svr fits. diff --git a/openhdemg/tests/unit/test_pic.py b/openhdemg/tests/unit/test_pic.py new file mode 100644 index 0000000..736a5b5 --- /dev/null +++ b/openhdemg/tests/unit/test_pic.py @@ -0,0 +1,99 @@ +""" +To run the tests using unittest, execute from the openhdemg/tests directory: + python -m unittest discover + +First, you should dowload all the files necessary for the testing and store +them inside openhdemg/tests/fixtures. The files are available at: +https://drive.google.com/drive/folders/1suCZSils8rSCs2E3_K25vRCbN3AFDI7F?usp=sharing + +IMPORTANT: Do not alter the content of the dowloaded folder! + +WARNING!!! Since the library's functions perform complex tasks and return +complex data structures, these tests can verify that no critical errors occur, +but the accuracy of each function must be assessed independently upon creation, +or at each revision of the code. + +WARNING!!! - UNTESTED FUNCTIONS: none +""" + + +import unittest +from openhdemg.tests.unit.functions_for_unit_test import get_directories as getd + +from openhdemg.library.openfiles import emg_from_samplefile + +from openhdemg.library.tools import ( + compute_svr, +) + +from openhdemg.library.pic import ( + compute_deltaf, +) + +import numpy as np + + +class TestPIC(unittest.TestCase): + """ + Test the functions/classes in the pic module. + """ + + def test_compute_deltaf(self): + """ + Test the compute_deltaf function with the samplefile. + """ + + # Load the decomposed samplefile + emgfile = emg_from_samplefile() + + # Generate smooth discharge rate estimates + svrfits = compute_svr(emgfile) + + # Default parameters, test_unit_average + res = compute_deltaf( + emgfile, + smoothfits=svrfits["gensvr"], + average_method="test_unit_average", + normalisation="False", + recruitment_difference_cutoff=1.0, + corr_cutoff=0.7, + controlunitmodulation_cutoff=0.5, + clean=True, + ) + + expected_values = [np.nan, 2.709522, 1.838382, np.nan, np.nan] + + for ii in range(len(expected_values)): + if np.isnan(expected_values[ii]): + self.assertTrue(np.isnan(res["dF"][ii])) + else: + self.assertAlmostEqual( + res["dF"][ii], expected_values[ii], places=2, + ) + + # Compute for all MU pairs + res = compute_deltaf( + emgfile, + smoothfits=svrfits["gensvr"], + average_method="all", + normalisation="False", + recruitment_difference_cutoff=1.0, + corr_cutoff=0.7, + controlunitmodulation_cutoff=0.5, + clean=True, + ) + + expected_values = [np.nan, np.nan, np.nan, np.nan, 2.709522, np.nan, + np.nan, 2.127461, 1.549303, np.nan] + + for ii in range(len(expected_values)): + if np.isnan(expected_values[ii]): + self.assertTrue(np.isnan(res["dF"][ii])) + else: + self.assertAlmostEqual( + res["dF"][ii], expected_values[ii], places=2, + ) + + +if __name__ == '__main__': + unittest.main() From d525bac3dacbb0cb68a9da9986d16f4279822a1c Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Thu, 13 Jun 2024 12:04:33 +0200 Subject: [PATCH 27/35] Update docs --- docs/gui_advanced.md | 27 ++++++++++-- docs/gui_settings.md | 60 ++++++++++++++++++++++++++- docs/md_graphics/gui/pics_window.png | Bin 0 -> 17212 bytes docs/overrides/main.html | 10 ++--- docs/what's-new.md | 44 ++++++++++++++++++++ 5 files changed, 130 insertions(+), 11 deletions(-) create mode 100644 docs/md_graphics/gui/pics_window.png diff --git a/docs/gui_advanced.md b/docs/gui_advanced.md index c2999b7..4a974ed 100644 --- a/docs/gui_advanced.md +++ b/docs/gui_advanced.md @@ -8,13 +8,14 @@ Please note, the `Advanced Tools` might not be available for all the files, as s ## Start a Specific Tool -So far, we have included three advanced analyses in the *openhdemg* GUI. +So far, we have included four advanced analyses in the *openhdemg* GUI. - `Motor Unit Tracking` - `Duplicate Removal` - `Conduction Velocity Estimation` +- `Persistent Inward Currents` -For all of those, the specification of a `Matrix Code` and a `Matrix Orientation` is required. +For some of these analyses, the specification of a `Matrix Code` and a `Matrix Orientation` is required. The `Matrix Code` must be specified according to the one you used during acquisition. So far, the implemented codes are: @@ -33,6 +34,8 @@ Rows, Columns: 13, 5 ``` means that your File has 65 channels organised over 13 rows and 5 columns. +The use of `None` is suggested anytime your loaded file contains already sorted EMG channels and you want to avoid further sorting. + If you selected one of the built-in sorting orders (e.g., `GR08MM1305`, `GR04MM1305`, `GR10MM0808`), you need to specify also the `Orientation` in row two and column four in the left side of the `Plot Window`. The `Orientaion` must match the one of your matrix during acquisition. You can find a reference image for the `Orientation` at the bottom in the right side of the `Plot Window`. `Orientation` is ignored when `Matrix Code` is `None` or `Custom order`. Once you specified these parameter, you can click the `Advaned Analysis` button to start your analysis. @@ -66,7 +69,7 @@ When you want to track MUs across two different files, you need to select the `M 7. The `Show` checkbox indicates whether to plot the spike triggered average of pairs of MUs with cross-correlation above `Threshold`. -8. By clicking the `Track` button, you can start the analysis. The tracking results will be displayed in the `MUs Tracking Resul` output in the right side of the `MUs Tracking Window`. +8. By clicking the `Track` button, you can start the analysis. The tracking results will be displayed in the `MUs Tracking Result` output in the right side of the `MUs Tracking Window`. ## Duplicate Removal @@ -111,6 +114,24 @@ Prior to calculation of the `Conduction Velocity` you need to load a file in the We are now at the end of describing the advanced functions included in the *openhdemg* GUI. If you want to take a look at more basic stuff, check out the [basic](gui_basics.md). +## Persistent Inward Currents + +When you want to estimate Persistent Inward Currents, you need to select the `Persistent Inward Currents` option. In this case, `Matrix Code` and `Matrix Orentation` will be locked, as these are not needed for the analysis. Once you clicked the `Advanced Analysis` button, the `Persistent Inward Currents Window` will pop-up. + +![persistent_inward_currents](md_graphics/gui/pics_window.png) + +Please note that it is suggested to [sort the MUs](gui_basics.md/#motor-unit-sorting) based on recruitment order before performing the PICs estimation. This can be done in the main GUI window. + +1. In row 1 of this window you can select the smoothing tecnique to adopt (at the time of writing, only the `Support Vector Regression` is available, but more are under development). + +2. Select the `Average Method` to use for test MU deltaF value. + +3. Select the method for deltaF `Normalisation`. + +4. Select the `Clean` option if you want to remove values that do not meet exclusion criteria. + +5. By clicking the `Compute PIC` button, you can start the analysis. The results will be displayed in the right side of the `Persistent Inward Currents Window`. + ## More questions? We hope that this tutorial was useful. If you need any additional information, do not hesitate to read the answers or ask a question in the [*openhdemg* discussion section](https://github.com/GiacomoValliPhD/openhdemg/discussions){:target="_blank"}. If you are not familiar with GitHub discussions, please read this [post](https://github.com/GiacomoValliPhD/openhdemg/discussions/42){:target="_blank"}. This will allow the *openhdemg* community to answer your questions. diff --git a/docs/gui_settings.md b/docs/gui_settings.md index ec04841..17a872a 100644 --- a/docs/gui_settings.md +++ b/docs/gui_settings.md @@ -8,6 +8,10 @@ Accessing these GUI settings is simple: just click on the Gear icon located at t Upon clicking the settings icon, a Python file will open in your default text editor. Here, you'll be able to modify the settings according to your requirements. +Always remember to **save the Python settings file after changing the variables**, or the changes will not be effective. + +If the Python file does not open, please read the [Troubleshooting section](#troubleshooting). + ## Modify GUI settings The GUI settings file (named *settings.py*) is organised into distinct topic sections. @@ -18,7 +22,7 @@ All the variables that can be modified are labelled as: In this way, the values that each variable can assume can be discovered from the API section of this website by navigating into the topic, then into the Name of the function and then looking at the specific variable. -The variables that can be modified from the GUI settings, as well as their values, can differ between different *openhdemg* releases. Therefore, the user is always encouraged to check the specifics APIs, and not to rely on this guide, which only serves didactical purposes. +**The variables that can be modified from the GUI settings, as well as their values, can differ between different *openhdemg* releases. Therefore, the user is always encouraged to check the specifics APIs, and not to rely on this guide, which only serves didactical purposes.** ### openfiles @@ -89,6 +93,22 @@ resize_emgfile__ignore_negative_ipts = False The corresponding APIs for the function *resize_emgfile* can be accessed [here](api_tools.md/#openhdemg.library.tools.resize_emgfile). +### pic + +The section named *pic* controls the Persistent Inward Currents estimation. + +The API section for this topic is [here](api_pic.md). + +For example, in the following code snippet we can adjust how to estimate Delta F. + +```Python +compute_deltaf__recruitment_difference_cutoff = 1.0 +compute_deltaf__corr_cutoff = 0.7 +compute_deltaf__controlunitmodulation_cutoff = 0.5 +``` + +The corresponding APIs for the function *compute_deltaf* can be accessed [here](api_pic.md/#openhdemg.library.pic.compute_deltaf). + ### muap The section named *muap* controls all the functionalities that require MU action potentials. Currently, it controls the behaviour of functionalities such as MU tracking, duplicate removal and conduction velocity estimation, among others. @@ -189,6 +209,44 @@ Additionally, you cannot remove the unused variables from the settings file! You If you accidentally modify some variables and the GUI stops working properly, you can restore the original settings by copying and pasting the content of the *backup_settings.py* file. This will be visible in the file explorer of your editor next to the *settings.py* file. +## Troubleshooting + +### The settings don't show up + +If clicking the Gear icon doesn't open the *settings.py* file, it might be because your operating system doesn't recognize the .py file extension. This can happen if you've never opened a Python file before. To solve this: + +1. Double-click on any Python file (do this from outside Visual Studio Code, or it will take care of the process and your operating system will not associate the Python file to a specific software). + +2. A window should prompt you to choose which software to use to open the file. + +3. Select Visual Studio Code and set it as the default application for .py files. + +4. After doing this, try restarting the *openhdemg* GUI and Visual Studio Code. The Gear icon should now function correctly. + +If the issue persists, please continue reading. + +### Locate the settings.py file + +If the *settings.py* file does not open after clicking on the Gear icon, you can still manually navigate to this file and change the settings. The *settings.py* file can be accessed navigating the file explorer of your editor (usually on the left side of Visual Studio Code). + +In your file explorer, navigate as follows: + +1. Click on your virtual environment folder + +2. Click on *Lib* + +3. Click on *openhdemg* + +4. Click on *gui* + +5. Here you will find the *settings.py* file. Double click it and edit the settings as needed. + +6. Save the file + +### Changes are not effective + +Always remember to save the Python settings file after changing the variables, or the changes will not be effective. It is not necessary to restart the GUI when changes are made to the settings.py file. + ## More questions? We hope that this tutorial was useful. If you need any additional information, do not hesitate to read the answers or ask a question in the [*openhdemg* discussion section](https://github.com/GiacomoValliPhD/openhdemg/discussions){:target="_blank"}. If you are not familiar with GitHub discussions, please read this [post](https://github.com/GiacomoValliPhD/openhdemg/discussions/42){:target="_blank"}. This will allow the *openhdemg* community to answer your questions. \ No newline at end of file diff --git a/docs/md_graphics/gui/pics_window.png b/docs/md_graphics/gui/pics_window.png new file mode 100644 index 0000000000000000000000000000000000000000..ff932d7898d02f3067ea884a0f2948a16fcbc6c6 GIT binary patch literal 17212 zcmbun1yo(hvOj!qm*5&Kkl+^F-7UB~1U- zAP{W5m711|mV!L5DcFw5*bHo9&g5a|0MJ1oejyJBV^bS*7jhGGODlT;%A=+>N^&bR z0ZL5{1r`MdG4r=p(q2yHs$Pm}rd~FtJZ6+ag2?#wZ zKaFW@0(Nx~pri!o6t@V_^=v-{gLoLwZ`073jI)W7lcuLn4* zc{-RgzcF_PyE>VgOSqZayHNd>#FJcjm90F?ZMDR$?9A<*0qF`*zUJWmFI4USpt5nY z{6!T5+k%}`fo)|jK*`C$!OroU^dvuC89QT3b4@EVm$!e7_!mLm+}`qu@~_c;Q=T}q zb$wc5b$csdPHcaTru7$1fRcxkne{J{_CJ#Ty&k;3c~k(K0W7lr?)4W*ModiE32b3y z3(!QK4h#7kDN*uQz;@)O{oNE;J6@oxm4&Ccu?x_MjfL$sBMT=ZE4vyC4=)=#FAEnv z3mY%XA0qkP8aNVW#xBPHSNo@pP0kPWP*C8NwsLj>J9+-T^~Y{kHFx~u>W@oXtKZw4 zoc#A<^BS8z$w+|G!3k{UYHDuwn^j=AKd8=N3m11|Cvy=?V4Z<^Sy)&B()1)(e%ip~ ztW2ydOq_r3;r`Yf80~+a8}px>F+Ux=f3e2T{NI@8|J~p}PAbsvk9&Xy08|O{zjX<4 z@o#NpZV#v+CqUPr`wScb`b^FWSe5^}8X(ZC1XCj;qrPE!T4X2pH{s#Y zZ;aa7+dJBa$%l-FjIz`9^j3-Y@j5}Up4Xnms06DE9PpBuzS{U^c!m51?J0W7YS0Cb zpehwvYezO1n0%;sfxPAd(5n#-P#>Pm0DKK3U{-=cWJLH9g*=OywDSSVtLmQd4tPgI1xJCKvi;Y=!GxYxVQ1u@XQM>lV1@K5)i!NF5Gzei4Im1FRfq; zkKSX8dK({i1D`szpb#`Ax`+TDh#B$$bAsi|moMLgJ*DG25-36_aw#GxQbgWSgm+|j zAW7S>5U^=M=fT0~QQ(4JRgjmBzFXxs!r=?T6bslG|H$?#DWN}|h0M5z#E9UA+1$wJ zItY_}$J|H@gsSvVjL6C%1`iMA3m@w13{pQ=H{xDvTmjTSl7qCiGYEv<{qzf!LXS=i z0+EAc#6{FR-XCQ9E9+{eZnw?GW^2I*iGc z$43ng&g_p|dHC+1OuTC3-QL^Id^q67rAuA<#5f=`#S2&WLrvN{P$V6bSziTOx0(g2 zyGx#CHV`o=C3YfWPG~vlJjYLtygTlU_5t&1+5NCV^T(j)_VR7AvaT*lAUUE*o;JcN zhW{_I<|W|5)A?w*#iI!XBIdLWuSel;yB-ElcP}6Wl82`t_0x-Bip!IGQX>Y2BfWsY zieR$7+ejftq@*cBsaa}iTJdo|6YyCG_0CV%%{o@TTud!`p*>o@@+zj-A80(7{C#Y1 z@%|RwW?jjDL@DFOZa^uUc5l{3m9Sr^wn+YA3-Wk{A4?2o4orJ$9*cR(>;CvqzqPy3 zzhbY)7j+AA@Ml`P&7Fp6_m!Sj1h-xok{CYFT>4?biZHN?zCa-1dYDOVyJIyTQ9bQi zY`n1Xb9roppXEQ>I@jb&^j9IJ!bb|!N?<|7E=V<>lQ@NumX`i?cQzUo71c5@Y;k59 zC&0R!$?xMGyIT+Q(0790rZe)aXTk7^)31Ewxp zQYwdqRU&j*;f!F|c?#T_v^>_f%NYz13Mx`mM>%ZvWjg7(!ZQ9Lj0uCDOXuiGv#a_V zIZlqcsP&ZaLF&lgO@3WUDI}rrw)T73BU9(|+ABC%5kyD&{lK^Mq6Z|>$BZ|7_pTTf zEm!M{^A1E!gn}1JmO0#aYD@ufrDw}R{J*6&6YBqCbAvsq+0$O}aijN9uk`KuVfuwA z6jq8*ajp$t&nMDzC3rS1o2jAU;o(}++p^1Cq0kv>zmrB-CF|I1R;Nz$FZ@tUm(L%5 z`u3fH_=~p}t@zw}ri#3AetqSQyQ(KV$~`&QMe;O=?>}`MGM_GuLX;wXYacNsC93@6 zU}=EjbF*Ojy6)w~VrGoZ;l0Bn6y&r4;l4n+BSaAi$+I?+lPQAtF6Eng`Uk`Maetb*NMb?6` ztS;!IjxssH5xE9oU^0DL$n;$8wH8#2G;IhIC8=vX&kGGgXOW0{cgL=5?BjthNd*?q zd%v^1jw6Da{MFGFJU^c7#|b^?fZMAX{x6re5R3Z9ITdb}!MirThFm`68_Uy#A6%3_ z((1aQCaO5E_;zrmHd`mzxn`p2j)-T4H!09Y!WN6r*YD#(g#hM+E;dxdj7gnJ_l}FP zIOe}fG`*zcx`idf8;t`bWP*-lqE0nbe^cQiA;sdUkRp)P`fy0q&-Kuiy2qEOH7L&% zkh?W4{ruP9Bu4+E@t|wQ#-4w@=dXkf;sve|+kqP;W7qUs6gd9+iHEG_N%r$C>{iyI z$5R$mNx4JSuoJ&2A#HS_lqi5B-RMW#P16~QJoVLoDdJYYOr_eW_>svsY4x77oJnk;wH!19kXR^pVbV;451ePv6zN*Zg;J(&;eA z-!r{2gj{ZWY^Qz8=DuFu+V;<+cZG8JD}c6G(&8_K$Y$el4C0h_!8QKwaVO$_N1GJJT69S&5BZFB>|-ReHws5u(nOdjGYk6`+X5h{LTovvd^CX za@MrfM*DrIAn-*cOW|{QIog9cOKj_7rpq?QqHL-w&23Qw$0L(xw{O{MFVB$e^-`+; zqtM4o;DqCXd>=6wuE4AG^uv-DTKbYBA zAfAP~?M7-#yLa{$Y2{&O1wx#yILlK{!9}6oLR({^|CI*{uoY7hUq`+5>nV6P@q01Z zsHHJsML-|G^yGlyc=k8L5#J2x?rzvbO`ub{@Tn{&>@)1nrF044Z}2B)J?ouH^mdfT zNdhnzVj#H&4#pd&0Kncz{!_QN+n{Fxx?D1Tsb6T94KT$?(p#>A>;|4Sqe(!-ljG~K zN;A6A_E72y;lN^J^v|f@*9rQ6HM#dXf}jDKGh%J&u_j=88lp=tr+|dj_2yhwT+BEW zr+w|TNJckeO8!Uag^0C`LBw0FdJD7w1uU6nQi||TD1X^R${bXp(f0xi?6_bqla5pU z`AiGuY%(`cAAQLinfuyG_q|&?;em=?j6=OTuByB+!j@sJ;UTrSswwN*1t`#)VxKGa z&Ff;66@{VnT^{eoK_&>dO3sl!$6KG4Z*=M!WQ^LBFqykN`A5}Jkp*wA50QtfGfs+S znwHSzp$EUMA0tsh%06A623Hzo?BIA0!h;;5XD!`iO^!qAnth!6(D&)hJ0bcFs6>X~ zDp9(X{Kf3lP?u;a zPD#(O0t4dpSx1s^>9B?`>AO4CqK!p2>4gO^goK3J3hsW6 z86G6emusd?EE!t&Kj$$bB{-)eBhdo6#=wn}UNgzQmh!KaLiV0f2OmX9^n%X^Vl!Vh zXcCrKOVC_7WyY+p+{qgreOsPGWKK-gIHHe<;G=vqjzyYOL(s@kk!aC~Y&olo?yf#7 z&M`#&tzZegZ$kVgsZornRtLP2E-JIDHaa(Vhr}K0)?szv5Lt{`7^5Ax8s=`XLiz|90q6{fdJC#!Kc7 ziIvNK+OILNC)lb2Jt*&U-(2TJLf`Ia@Xg45s0fN}-2S+?^iIv-I(qd0NqO@u0g)S* z;B}5|Eq!UONPWBUFc$jwy#$1PLh(>ypu)D^fINd%Uozv}w-*wL&>(DvJtn7R676X^ zLJq4pw19wdLQ0CMdt;w7Vq8Jm=Rekq8mMTi*hm4>9m3b(%r=Aa{47 zWdf@|?{pp;d!%OOm4b|~YgVtv6O5`}OF5yxnf3W%VB>~q64sVn)?}Wke&MT{MfIbN zfVLt;f?^F4#|v}KvplxU+T9}0AQ*4oMB;_i2h|5cFshl8Bcuh%FW4o4*vxnoXMBY} zB+%dqHOVFL#~FvF3p}o4UDw3Z>lahwdY2V7=_u)bIu&SgKk8G8;r|k^NwPVee0p{^ zaqJa!VXb77D z1w~|0WvSK+zL3Ged-%rVZMo~#bNT)i>P?hElG}h;#RKp8!LZlh*A<(v8#w=QY|@)c zi!=);hjSfIHq-PTC+qYQP*6o%H_78_xoYv(QfHQF_=ku3NlMG9ebrkL%pupeWe_4P^ zKZj&BlV*qtf}}rvy!|yfDg5pZ+T(DJhsB}xn|0dJr5S(FEsxceu`SZW* zDwKQC7rbd`t~3hcqbDwqeM?nfC#yPO{ zB*DsP%{ABt8~r%_0m3vSD}D&DshYW$%YK3bF)ANljj)6@_)==q$q{oZbz0P<(eO14 zfx1AxA2o8gtFCmMAst&JqZLjPR#Ajwgcv2yN>GN(Z zy)b!~YgZtbHA(7JTx?^coRllXZ4{!Wy~_bAt~RN!=lt}^$sa8$(nU(}CChcZu5P&= zo~SZu`{++jXV_|Z^4g7`5sAeeIm9qC<@z)>X)?Pb@vy}@Y-RdQoTlC*-NVm}ja{su zVpn^S2_)2OJE^KEx)rN=z6vv7z~)OYCqoCPk;U7N%Dh$@G5MKl9+fnRY6fj^HSgXU zv+gxx(#)qC)Pth$w)6ZXH|9vx4OM!9E@Y}Su2?#V9n3x2xjEtX51U#T(%cq$H z2Z_JR8VTZoO}fo&N3Wxt@dCilIp1{h;0&7Kp9MKDBjrr%AtCl~`*Pj2l*-6^A%nQ> z4u^z=v7&O9G}RgJ=5Q%p^g$mq?`hp{ z&$GI0NiwPAmmhX~kNn^UTgy4+EaI>8m6W_q_|O{r>%h#NXj_EJar}#f+l(So|FJ!J zLkJS`C}p)JahaKOH48BTUWHyq*)>X2;2@(%^}C%6kSHhSYPuVI(@8EuR{aa;8fdQf zc+M{n=zra{mvvjcY2!XS* zx}%`vEAV!6hGMMG`*N{A9l2zc{>8x6F8JJtYiD$MxkfjqM<1on9YyqgK&{Wr51pMH zOK3>Fb5#~#6*7IKbOMH6rIA%GP^f#|hcsPEzlTb@+UV(ocx~IzZA?$6)+iq`t-kC) z_k=qt)@)1geoo>DBA!&5xpl?sq(u{Dv02n}ttDD3b{f`=r>nXS$&U2ECg)krCV>U@ zq5rZSM@Ep$8_4=e03CXh6&8$3NT-M{$_YwLq;4-;U_Bz0LlGHQ%qXj0fj8vp>`oyI z$bCfSKTJ?WH=M+1iZVyQX~s}NWuM#+AxV3?-QKK(QItwjzRWti-29RER}K29 z9R0?nhK*+Flz|H5WQDEPtiHK)YWL`~2H;T|ej?P>ZPNMmaH%G^Zcy3J45^r%_nq1} zmGB^GFIou~z(7B$p@$)!9DN#@#el`NR;QsvK;vZLS;yz=l2RYs0w?p@!WU{Oan~X5 zBe|0R$m`a3m(wxN=(Yr9?1b?|OmdNatb(&bulyBJq~8jC+_8AletYY3!e6dmGuHbs%Xf70YMPP-9lAaeu4IE+rE0o$*(je|Xw~4m=#7a(B!rX@KKMQy zb=f?pX|-(ZgizxJ;|lhp+}gpB+&{Vvg0F*A4R2F4p@l(?FF^6l6$FJ|(OhnCVSC|6 zb`{)s+$BBXs{2kGeaNK09?v%EW=jMnf2NN-&y}VQ-yEHc*_&ocTB=yoUsVhlq}TXd z$0YayBHB%>68nmM1hxuQnvq8%dwRY9nRA?SM7<_Ge4-}C;n_2$I}$ajV68GP9H+LI z0{K+G@`W@{gGAP)t^Fcs+QhNJT&c7|%WY%@sT)&4;wn{2&%iHJk!r5|q!%k6Mh(yH z7JtP@_ea+A!9F9aP=Z}}DSs*}4un-Il^k~WjkSeJkPN{VJg_?=CAts0GnWU<4V@UJ zEWEwa$nq`4QyAv+%VO*BplM{)?v)qy$6ij{{E1rFDy$Dbe*=gA6v#;%AO(_}K>rn4 z{a3W{dod%WSDk(09fOr>a|?MmH;eaCSbPy4E%1QgFDY8XFYL7pD>=Haf>ePh4ymYo zl&nIl?ctrrtth5(Riz<(jo*@vwzX_LLFz|MN+M+T0uXwdkixudnFSVC8O9qqV2WLh z95x`^#Ng|1)v#4W6#v3JyGwf147Ihb2%_qVyCZ{B>wXNr_TIUxgb_ASNs(==8C{;W^yT+s<$}%jHwvfBBYvZfTUBt?+56m-fDG!xvG&etSc6H$OF) zp7&-cys-|Cu_jm^E_kd?t2Sz$swy%a?9@Jj-e5R7KZ|K7P~%fI*KLT%a}a%MqpT%Y z)7p}i-6jR#pN?C!kbTCiRj%G}Ukks>RB~FY2oQ^ljh-@h*Rs>s(oTy)Nd5I`wdRjxqMqiy+ggaKP&c@n zO4?$dnk!qAZojE9Gnc)Njv?HoZ!nU1v?`0SqMm8vw=_NXyy?nQsL zu?9XRNf-n+lPRH90fisUKd6-32Q?NK?|dz6x1KEaon-HPuJ96hDn|0{?o7i$ z{gT2(2o$_sS<^xK0ljdIQ1slYUip|WT#RJ?2>MupdpVKaE>ZsXL1o`(tQLX|GgeJ@ zSWkg$T_4CnhI_*o`t@6G zOzO1FmP{V@QaN1{OSti%dX!S>!*iHbkXZms<=p8&RP9*ReR<5=k9y?? zPE9d`voyDD-}y4=t9L9(hioP({qY<>(MLHlVwzi;vo4mhb~!EuG*^^MWWj-C%{k$x z8tc(|ljj*+{7hh9D=^Npfs|8g^WK!(G_e&<-LP|DETSjVrp&!6WAz*Ne4~vwel46X zcwO>g^%q3#05E|W{yEZj^Y>Hy4|wN2~3Sgo#w z#q*Yt87DDk9zsdOV8kV2h=A+Q3{t0%%HzBWdN z;+6~87lrk%q5M`ocNTmzCCzi$mzkTtyEV-G$oz5=>vojfsqF7OoA1`5;3&IPsX1Tv zO8MyGPK4ci$_W)U4qDhv5sKwkL_cnB7*XNivzIvBsbDv$|-d)I2^!QWR*w@k_4&t#DYCg`_E7DS*iO6Z!;q_VpU ziq1!P1CmLs;#qql#}c4GF}t0T5ab1dWbJ4T_hyA9t$fap!WAl|S4UHH_%Ex8NJz!z zd>V<^vp;(;aI>t^p1jI(YJ-<#>LFu=-q!2&F=fRliC;fyKt_u_As!Lm;+(ifup}HqwK1pdJ^lKSjF zi`};o)H;ViO7U~yU-Om4$MLIFD2u%})qy>agpx0Z;vh~Bv8U}&_OF>4zI}T3Yi!Km zFYZSmcMpm_BZBs%SK+YkqMQVJuhBqH?)}#Vsfql-t;_e~XB!Nm_S$va4MV~-tF90M zjN3CQmR)K%|EA3Hb27Tvm?~{;`P4V1>ejU{jZ4yI2(nOn2U6YK;^m5vAJIHk9L`Qk zk?df;-_0*EblZwgt7y0h8j%pf5%|wY&f|BeUJaCTmEnO}J!hpO-WVL*w}=!3H$}kn zQyjn!7Mdc0-nf}tMiL=t>7CZY7~i27=<0LLe7vBFhSym0>R}Vy4&1~`NnaSUqa^=wFuC&*ZtOzBBfY4k zF~Z(zuw{t@^0~w-cY9kHUUPpgc+?EP2r8h9MhgL^ot7CilnJTKVTrGWH9&M)UNfU& ztPh6@$Pmm@*5v}rY@qpdTpZ4GZBd>+hkGGsZUwR+Jsb+aSX*$H5A%NLG(AI^N%l%C z-^nnJCYso%_Vuo5y??V}LUPkLFoFO-KSCg%ZE6aN|72H2A5c|Tz5LyZ%{3GCCEXMU zhifP+NSn?T3f`qxK%GN!)H*t@SoqbI4aijC(>ty=uaHk z6IfYm@B6Ho^$x1kV~6&qt@z@z{3%jqVC2)_j!7#vz=TS`#8onztVPgyxl+UW26Agw ztH+bAf=fxo7SZR<4-PA1RX(56@P-ETuSMM$*Q3q76Mp2@bN_jGX%*4w{aBFArVVaG zdpN~ax%>%n_ZS~em{*cb1=ToETefgnz}4>vd`3-?jbxEzoad=j#X6fRuQpZ4>>$do z=;0nW6UcxI7_0lzzZ*C6^Q(|wLt%wX(Pk%W>#kiZGF;}m-OHs2C__~!<<@15FkXFA z1vF4YVX@4;Ta3N&dyMaE{dF`X5LV7QFTyMfAPxdt~|HEfD zz9IDp_l1Z^ktX^6+a1U6Gmz%pyYDUdsQ+=Wx%aw@h;E9P#IO;)CoAN>4m$S1ziG4n zuA>-|pD<${WOz-57buaY={EFw*S*VPDAgXXt+$VmebP`R+eIkAVo15zp5*&$eZx~- z<{#fE2lpbq$n#p7rnEQEakZo970+|&s8g}!QxmoIRCMp2Q(5fGW~t0@b35J5+}0T1 zZrEizg5FSa;DA6ReNVLjnKa|m)LoH@nnwjM!}bH*;e*IEV72E{IWyLV7Gm`9@;?Yq z6s)87VaWyJQbPcw1hnQD6yT|DB2uXZezn2(Xf`}a^3NB-a>#c{J710!)fm$T{n4@heWBNms^QYtgmt^%Rp2|E7dXm+~ZW zx^cBJFQyeS(TW5TGp~l?p}F-QUR%*uWy&m83Q9Bf(?3Jon=XT@l&Mx;;g`dt$NYDq zi2+ys;IfMX6g?oa-=F#*c> zvA8+KN_BOs>S8>{l;&w|N(wwSv1t(!B@v#bm5tNaT>a$O@J}&2;(3i@hzj}uIyB*k z7R~GW$uNXtYH*5|>bHlSxw`?=LnI=6Bd3^pd7n!?{wMeB z->$H68exg>g9pUr`fFOMu0h8OxSAaq0hmgC`+BG+Atd=SgbBEf#)qssacDJ0*gh*% zYNdO{{ji5ujs+S>o+D7>s?#^F)sbARJj15?U}Zg)3KdXYhL&QjzmURO&WU0gy50!o zsjnNwi2_Od$LYj#}I?|;eS&^4p zg`jqEsjww*)RTLe^)v6mqQ!i%bV=F2H-7Y+*nObQ%O^CP2@N{!G8=tkRy+)h|8VRc zCuEpq-$~Zwp*Pd04N5l#JBGwQyTy`i({i<*IS_KU0UI?SZ}z9+8R$ZRwCAoL7|j3+ z7?hi?cVeAdw|V2Ls5D1!r-AeN%BvuGt~Axu$m2CU6@>;dnXiM$RgcXIPxd0@-%h7) zaEx|;w9`y^R7(v(wf6HAx@K!_6y{7H>7TJ%G`zgy0*)l07kEN@I+R&%jmQlFW8V3l z9#Y=Dm;u_vdusd1AS{V*f=V+!-?>iA4c#n0e>fVcuL+C=aWmdMMvwz2=hSWFjTiIf zNMMvMVa2ddui$#M)l(wG3M(mRoNrId+QwXv522dyYLImzqOm>(1}0|!S_qGz{Lmw- zu}acT{GDFD#?eZtYk5jcEGp=&8a}6&z9}(GdC{}8wGjNhWluST)$Do~B~q^*#xmw+ zE&S{biMgxHFRC8Q3l9tfQ@B?*y(0@By1Fb;sYn_Lu5bqjV(l*?Ar@$%2=Rlu#cUS z>Ovooth#-}dz(4$vM5k0Vl%^|1X@1@%MpEDOU19O28KI%BL-Wqe zvTjQ?kOaQ|ug)LCconqO&O6kFO(c2rtf7died+4CBXkxKQO-AP>Jrryg0J1q4Qviw zKS^yHJEwFSSk;JrPDJz&6XNW#FlW`4E!QZX$_XwuFfg@cBcY$Lc=C#b>#CzSlEdME zP3hHT4YisT7LzQ&hcoy}%@eem{$9U~Ng&$4WA6iC^GbGo)u$w?F*Om_+*KVPVA_@a zx%iZ_N)rK}-^gLIMz#GsVhr5HDweHtFft3Rd^q#fz$s(8hn>dwR#9~o%>ecpE#C(6 zp=FL6-bn^b2AvReGVLG;aI)GQ7=5bp?^-zh*M4>FnnA?~OJ%Cvz?`KJ`a zw*zVi^eda_EmV~|$7Ycd4DGW;b8;G;eRTUTR*;#u;zv&>yreGti!M-nYC3eVuc6Rb-j(dS1xvKq znG0f;Oh1+pUPW@?6T|}P{?^>f1b9KHigR>IbsTFRF;!JPh&XR|=dC3tIeI$ASUr|@ z*a6e4eScE#%Ax0>`+;ZKZM@|u5>mcqvzggQlbz}W*PUqUNVb#sMOB6WM+;_h@w_HA zO7oSo@y*hwo4Tc12)&qk!pQDzzI;}rFN#%?O7fAD!9H3lYR2|#+VaRd2D?F5yLc~2 z&+5$^8Icax(p#=gnbHrDDmk_)$hz?LVwiGq!SnO;N2>MT7SK#GRL;*_K=a#=g3Ga( z;;;<+X&nydj|+@#n;43|3ZveCufz*~9(5tvtWO^9ak_qg`$x3`Py_)JAml%0X=){g zQ{~klxb`SRM{;h?BK~@|kOWS50;qVfxjE0=nXox=3nTy0&c*H?#08*p;fe1SlJrpJ zFiYdX7Kh8N#6MnR0CxF&w(w^L3hXn?$o~ta`(L2F?0L$MmU}e)qlui~7El7iVF0`l z(bGBsK;vYk|17ouSPQJr1NgvF)Xf#rtMa>tIk^lW?=&%PL$v^q$Lh&sa#jX_qt`eV z(oiH!{+WLwe@{?mN_eQN%`|BCMMC$&vp}uJzwuHjn_5>mC$}9&3ObD}Q5*O0u6aYq zu2$kYEnJ&PvnR7j=_Q7FV)3cJbYk3Vn4gxRUs;`=ZGCfOm(m8TkQD82cUvH18C+Ol z(N_9`kI%sGnYfRZ$SLo{jCYkg=dwYWb@f29_T8YZZW4Om+#W4O_*741e@iu9THRCH z7*djX>)FiR-D+s!Xbb>xow7mw4j-!l3vNZ=oD%I>TYQ9LU|s?#V@L+TITU(&G#a<6 z7Mk$SS`??&N3Hx&AZdmt3~^4MNkvWB!laTw|1hDhD`qg2|!Kii#dOqKhBD0rSIbe zzkcvuuRVINc&G_H)%Uz31SoHc{>kyCP+Rqg0TycOw^{q01WhDYQQKEVBr@a`+Hdbp z^PS81_%<2CV@b}MGpi*zfrfR5wxtH0oLVaLOW(SVXY{}9lKO0qaFj)N=V@n@>rwM3 znBx)M#v*iqb-ZbOgR-cR)skdQrGPo_@&-+1)+fvYn_GgsCnelrk5-AzxiVAxZJH(r zTyK`&LH4I%8aryd%or@0C3?v2vWJ2efxb?FsYKu27Gp%fOFz`Phmm@bD{wqjp|)~1TpY^kon0;)$3pNCxv>8$1e=yw$n z$`%7fXu>8;&zed*^#=f=5*9RU(qx?3AsS8pl7nq%noaZFO@}gcNL78QEc_GW3--7T zM+mSG-i#B6$VH!~EM<`UlFg4}fKYns<@ES4w8QdRH^-l$w657a@36smJ=fj>^92u8 zqkZj1049y*NP$l(o-pH9(bh7>nY6GX-a|GjD=y=5IYpZUqT^0q04EpMel|uQAV0PK z4V_7y*G6!~{&LdnvyU!2vih=kFD46*kBZPfQ!1%4DDi+NVyjhw%QRcrwRPMPHo83$oodI0T!GIJ>$ZG1XgswIe z$aBXp4%e*YG6npM?4xT)sSKZcajV@C!4Fw>HQb%DRv{nHPuVb_fG6c!4ra}_tGWSP zVMFj(>*2Nj&QvI%d>j0n^7LPX6a1MXsso||2KGq59CYE^yOmx*hlB@&+_Mq<6)kE? z9s&srlYWsuYSyF;B7RnBfx8#r0>fYanNIE#_kYf-{Gla-YyHo-{Lj{V)bB#Z#Diy>E@tE{blc`ynXcB# zF(jR6ZLC^Pv*mM|JJeYwW^Vn$ub4Ri}&J?3$Cv8Of@`GoJk%s03dTxzVOA zfwhlv-+ldNLyZ@t?L9-zEq?k*)sho*-GQ3hl0?=C%PM&}F8?*;7Hzp2LJ8=sR+zc{!?&3&F5d)eR4JwN zhI#Aw4Qf@{mdEJ7v}M=MIz5M!tUk{62n^f$_4K|PSlcLkTR7Tt5z@LGr8$1B*j5>o z;Du(~l%|Ord8)4n@kIYsxSt&D{>6JZ;0mWf2R`F?*btFfC0Q^FMMJ0YY2fEy4adZ3 z?+S?a_mxcs@(&H$oaF&$vhv;IXwd9q4Ca@)!%Ty#JA_O#L?O{=OW6h#sykc6ClS2> zL{yTkC}-o%^<2sE>mz)De9pnhYfdVa9@r+OA)nXxgUtKLbBqhAws3(y!-e8i?h^q1 z7^c)rMQ@P^60YWfb0h|j~)uRV!1o%U&68q!aBk5*qV-#El@m-%9wR0zh=HJKsB z_#B=57xs^ypXYT#l%55m+?SZDYnICWiOD3ikDY7p-_2m+EciB`%DCT_q#WM2E{Jh47Y&D`6`uw?1b|-09bnOUYM|BwR3g0YS-|(-IC_# z3%Mk5vOI?^_C8B7#3v*rahSGJAxI|UDJcO1q$Fhsaz_iv5;LTXcIwuirNB-!H`7qd z{b7lpAX!EU4P^2R9-zqe{KBo~jj_vCx*If>W}M<6E|b{ymoKrRE-91}K2$vR&YpY+ zIQ}x}X3jV;P?&wK>BG9m)K6h*K#O@-A?48xVq^La%A8=I7OZ;C(|pLHejxRdLKe<; z;_tn#C||S>*P3~wGZ`4yYqC7wJbm;$uHLu8zU1d)v^bHH-k153iv}pS1lFeM+~&Dm*y?y@l@_(}22PUa_q^jPog@x0T;um;(mG10hJDOM)ET(YvAG9MlLC$ZTt` zbZE(ZjdW6X!TMe8mI2c`k}(gL#HMcNL8rUbvU`qi?+-s!AB*;dR)6Mh?3^7A>-v*` z`x|++4{ON1wa%J6Ue|H3Kw*vgYUI&&$Dq?CX@@0^Rc3C6ve&CQh8slBLN#0?J;7Ww z;A{o9v_1;`3tH`Q0CMQz4xSFXHoO|bPgT|bjLCIZTV%ErfWdB$F0+<5UQOD@C_n~o zhOMQ2>E;{@{qj(QGa3y^B;1HW#q+CCicChpEuioUYiWY{L`MsF?3y3Sj1slv zxNxEwzpGPN_Ox|)2)5C%0k#)6Y{t_8^iHvC^gJR5i=5XR-Rd}hQ$*>9ckGmR4(k@r z7TFEzZ8rMjG;|DW|7P25bMq|^{Msz)9*ZFGh~6HY-k;nBS#wRGLy)_zVuvZLPcMO! zMRgFd-YKwg=+Pl(t03rXt+t(uWKeF+wn%KMln>U5^N`R6;iXa7`F-rzQO4D=oun6m zknj=U+>H25XPh3;n4c^GUjok$`oFF}mN-}dEN{Z}(^ttJJDccxyeG*U$??Yi@j}q2 zaW(Eb%rjJAmEztLvv|0V26pjeg;|nRw;H#9)v%okV$zvO*an}5@dDCa(*B5z+MI*6 z>Uzl!Yn;=4<@sHNz;ZU9r(pEiFCKW6*>k~YvP?6p$aTEuDddEAQSHI(><_5$ysK@n zd6wuvEj|PP`Os6DJ`lD1uDJi_z40IKcl|AX`m;>`pE^O2N&iY~a?kRAjEmYq7QyVU zUf+4tYf(tus~>1ID+JR@hZlo~BT%L0Dn3+*fm)ZST8v~5;1+NL zndiZMT3x}5Hp3W6V2&Qx(kxcnpufGb7VJN(*fkh7P%qHzyU?{zFQE5uO2#+m(e!22 zrMY>DLwII1?=&q0=D^1sAR8Xy=+_ciglb^|Qx&T8XJFeO4-2$0n}R?k)z$g_R|9oaXv_~Pt~+Gu9yfKz`+c^HS>8GX0XudvaY4! zE53nWU-b3a#fI>KxEsv7Y{ioQbuENVQcLHueiUX4gHOvHY2o^m;!5CxJ`fq0jPC*k z%pj2-+JNB$6-vPUIz7VD)ntyY_2quv!Gggdss=qW#`Wc0lp$w2lyjH;i9Y|riPN>n z%>BchwEc=Y^J|>0JL9HYxz23DKoGFvzuIFLL6N!3qaG<~_YV`L0!3eZSZ%m}4VJi{ z=)w3a-7hB?rskT%4G*JHGy8@YWsPRrc(l&m*>@cVg$>9Fc28{UQA@c+*>|jzB5UlC zZT&1R71$<HQ2xSzKqwJZkZ{@4nGW~_}OnEC-6uJ?@|$tRPoHp`BOhAYaB2Q{q6 zgQ4z^lh?P^Y-tzUBfhI8oAYf}9(fGO7M^+oZf$0dzp~(+9lMJ3j+9f(*Lw~_9m-Cy zS-=9-NP$#!J+Eu+S;lx~MqY{flt&Ztv0* zjLXSvo^fCGhgmpA_wV~R5^PcEa^t%bp+EP;@HiCK{@DZjjO1yo6Z{oc9=?F9P${UM zM~SF!KPw4^#rnI@3~;F!f*&#Un%y_nuEgQKjdHyA65CveC#MYFNI9nDn>o5?2Bn&( z4ca+4b~Vh~Lf8ngKkr2G9CFT?zeD&+32akwV4E_3WjXIs(kl0`e6u~+R_2v9SZ}e7 z3j;&NFu>3F)iqZStf1@T@hmwdh1Iwg-+l;Lk$@dLBELH_%n1_AwK!4`>fUX0ko|bd zq<5_$dErtuqfs#8PBiqOF_ka&Hl}|L%5f{X=^g_%%?h_^)QzdHuJ}8!O3*-UmygjG zFo%!UxN_}5k26}We-=#i-A;Mv)zB+@kEj)%Q4)cv<27btS#&gx$p<@hN#?e>W2Q*r zts0*GeRH(YnbT=C|0|MZ@1wrFX6M>Dg6lLc;$`e&`Xqf&^KIUP?a0dED3d7x%j0oa zD26*139VFD6oON`VKw=f`OgdrL_jqLlcoEAolP?J*d7iIr8N9As8+AmVpFoZ#BRZNuey@9p(-m`dFO^vdsHk<9PNpo!T#zX*z%^EC-ir#gqw zy@h(v_CPiMhya(vMhh?^T5i*|H$2GsJgAh1`Eg5HY^v#fZRl2g@W^P;lN^}-w7)bA zy!R`9a)#={43ztGVok)t4RR)gH3w@4Z8v0(3=B_pdEg%sXLrnY!$U*gXY8Iu!2u7+ zl3G826%hw+4OQG*uQ&;f?M*y3GzwRQx5xq7nf@Jg!;gd|znnK~QU>+?PGpZwk6(dj z)Pq-ITTGu`p-ow+2ZOymyr%{qm%!zZ8^S96IbU+^yf?;rO?I9XY<1D@=HtWnLIuB4 z=;!paG=mrOJ;2LXBEc93Bu{T^{rxrHKLbQp(DU~k50@64KXHWn1 OLPkPSyh79{@c#kx@L;_F literal 0 HcmV?d00001 diff --git a/docs/overrides/main.html b/docs/overrides/main.html index 43eddac..d22b522 100644 --- a/docs/overrides/main.html +++ b/docs/overrides/main.html @@ -7,13 +7,9 @@ {% block announce %}
-

2 X 🎉 == Double News!

-

The 4th openhdemg beta release is now available! =>

- See what's new +

🚀 Release V0.1.0

+

The first openhdemg stable version is here! =>

+ See the new functionalities
-

We are thrilled to announce the publication of our Tutorial article! =>

- Read it Now -
-
{% endblock %} \ No newline at end of file diff --git a/docs/what's-new.md b/docs/what's-new.md index 4224d8e..4ccd32a 100644 --- a/docs/what's-new.md +++ b/docs/what's-new.md @@ -1,3 +1,47 @@ +## :octicons-tag-24: 0.1.0 +:octicons-clock-24: June 2024 + +This release is aimed at expanding the functionalities of the library with new and improved analysis tools. It also marks the exit from the beta phase! + +### Backward Compatibility + +This version is fully backward compatible with v0.1.0-beta.3 and 4. + +### Major Achievements +- **More functionalities** + +### Major Changes + +- **New modules**: + + - With this release, we introduce the module `pic`. A module dedicated to Persistent Inward Currents estimation. It now allows to estimate PICs via Delta F and will be enriched with new functionalities in the near future. + +- **New classes**: + + - The class `Tracking_gui` provides a convenient interface for the visual inspection of the tracking results, improving the flexibility and accuracy of MUs tracking based on MUAPs shape. + +- **Updated classes**: + + - The class `MUcv_gui` has been optimised to reduce RAM memory usage and can now be expanded to full-screen. Furthermore, it now accepts custom separators for copying the results and pasting them into Excel (or any other text/tabular document). + +- **New functions**: + + - The function `compute_svr` allows to fit MU discharge rates with Support Vector Regression, nonlinear regression. + - The function `compute_deltaf` allows to quantify delta F via paired motor unit analysis. + - The function `plot_smoothed_dr` allows to visualise the smoothed discharge rate, with or without IDR and with or without stacking MUs. + +- **Updated functions**: + + - The functions `tracking()` and `remove_duplicates_between()` now allow to directly show the `Tracking_gui` after completing the tracking procedure for additional inspection of the results. They also allow to select whether to use parallel processing to improve execution speed. + - The functions `plot_muaps` and `plot_muaps_for_cv` now allow to use a tight layout for the output figure. + - The functions `compute_dr`, `compute_drvariability`, `compute_covisi` and `basic_mus_properties` now allow to exclude firings outside a custom range of frequencies. + +- **GUI updates**: + + - The openhdemg GUI allows access to the `Tracking_gui` during MUs tracking (not after duplicates removal) and to perform PICs estimation via the advanced tools window. It also allows tuning the behaviour of these functionalities through the settings window. + +
+ ## :octicons-tag-24: 0.1.0-beta.4 :octicons-clock-24: April 2024 From 3ec7e96aa6eeef91c8c7c77c2b095b9724a324f5 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Thu, 13 Jun 2024 13:56:37 +0200 Subject: [PATCH 28/35] Final adjustments for release 0.1.0 --- CODEOWNERS | 4 ++-- openhdemg/__init__.py | 2 +- openhdemg/library/pic.py | 3 --- openhdemg/library/tools.py | 3 +-- requirements.txt | Bin 362 -> 344 bytes 5 files changed, 4 insertions(+), 8 deletions(-) diff --git a/CODEOWNERS b/CODEOWNERS index c69ef4d..857800c 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -10,5 +10,5 @@ # When there are changes in the pic module, @GiacomoValliPhD # and @JABeauchamp will be requested for review. -/openhdemg/library/pic @GiacomoValliPhD @JABeauchamp -/openhdemg/library/analysis/comput_svr @GiacomoValliPhD @JABeauchamp +/openhdemg/library/pic.py @GiacomoValliPhD @JABeauchamp +/openhdemg/library/tools.py @GiacomoValliPhD @JABeauchamp diff --git a/openhdemg/__init__.py b/openhdemg/__init__.py index 578cf8c..a05eb9a 100644 --- a/openhdemg/__init__.py +++ b/openhdemg/__init__.py @@ -1,3 +1,3 @@ __all__ = ["__version__"] -__version__ = "0.1.0-beta.4" +__version__ = "0.1.0" diff --git a/openhdemg/library/pic.py b/openhdemg/library/pic.py index c03296b..adf74f4 100644 --- a/openhdemg/library/pic.py +++ b/openhdemg/library/pic.py @@ -47,7 +47,6 @@ def compute_deltaf( ``all`` This returns all possible MU pairs - normalisation : str {"False", "ctrl_max_desc"}, default "False" The method for deltaF nomalization. @@ -55,7 +54,6 @@ def compute_deltaf( Whether to normalise deltaF values to control unit descending range during test unit firing. See Skarabot et. al., 2023: https://www.biorxiv.org/content/10.1101/2023.10.16.562612v1 - recruitment_difference_cutoff : float, default 1 An exlusion criteria corresponding to the necessary difference between control and test MU recruitement in seconds. @@ -88,7 +86,6 @@ def compute_deltaf( Quantify delta F using svr fits. >>> import openhdemg.library as emg - >>> import numpy as np >>> emgfile = emg.emg_from_samplefile() >>> emgfile = emg.sort_mus(emgfile=emgfile) >>> svrfits = emg.compute_svr(emgfile) diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index 688dc72..c2dff8c 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -1255,8 +1255,7 @@ def compute_svr( Quantify svr fits. >>> import openhdemg.library as emg - >>> import numpy as np - >>> import matplotlib.pyplot as plt + >>> import pandas as pd >>> emgfile = emg.emg_from_samplefile() >>> emgfile = emg.sort_mus(emgfile=emgfile) >>> svrfits = emg.compute_svr(emgfile) diff --git a/requirements.txt b/requirements.txt index b7f523cf87884e551ecdf411858b09dd3a465b18..5703586c4834fb944fe0a669ab0ec4ae310be1a4 100644 GIT binary patch delta 32 mcmaFGbc1Pw5hHIgLo!1qLpBhXFz8MWWE5sIWzb_VU;qG*u?GAA delta 64 zcmcb?^onVN5u-&GLq0eGkhBFtLk2wtV<0wS;9@9d$Og)!G9)q-G33EzO@Ok7 F3;?o~3nl;n From 90c634ef454487645042e925ba0608341c2d2282 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Thu, 13 Jun 2024 15:36:38 +0200 Subject: [PATCH 29/35] Update tools.py Minor fixes --- openhdemg/library/tools.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index c2dff8c..5b2bf7e 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -294,8 +294,8 @@ def resize_emgfile( # Double check that start_, end_ are within the real range. if start_ < 0: start_ = 0 - if end_ > emgfile["REF_SIGNAL"].shape[0]: - end_ = emgfile["REF_SIGNAL"].shape[0] + if end_ > emgfile["RAW_SIGNAL"].shape[0]: + end_ = emgfile["RAW_SIGNAL"].shape[0] # Create the object to store the resized emgfile. rs_emgfile = copy.deepcopy(emgfile) From 3b26ad7b46d76d9f8af1f107b0fa64bd0dd34fe3 Mon Sep 17 00:00:00 2001 From: Paul Date: Sat, 15 Jun 2024 18:26:43 +0200 Subject: [PATCH 30/35] Fixed: opening of advanced analysis window --- openhdemg/gui/gui_modules/advanced_analyses.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openhdemg/gui/gui_modules/advanced_analyses.py b/openhdemg/gui/gui_modules/advanced_analyses.py index 49db8cd..d7cfa28 100644 --- a/openhdemg/gui/gui_modules/advanced_analyses.py +++ b/openhdemg/gui/gui_modules/advanced_analyses.py @@ -328,7 +328,7 @@ def advanced_analysis(self): # Set window icon head_path = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) iconpath = head_path + "/gui_files/Icon_transp.ico" - self.head.iconbitmap(default=iconpath) + self.head.iconbitmap(iconpath) if platform.startswith("win"): self.head.after(200, lambda: self.head.iconbitmap(iconpath)) From 5ec71b5131026b0ed0d61303b063234730a37e16 Mon Sep 17 00:00:00 2001 From: JABeauchamp <103655220+JABeauchamp@users.noreply.github.com> Date: Sat, 15 Jun 2024 14:34:34 -0500 Subject: [PATCH 31/35] fix bug w/ only last mup discontinuous --- openhdemg/library/tools.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index b93b885..05836cd 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -1244,7 +1244,7 @@ def compute_svr( ------- svrfits : pd.DataFrame A pd.DataFrame containing the smooth/continous MU discharge rates and - corresponding time vetors. + corresponding time vectors. See also -------- @@ -1323,6 +1323,7 @@ def compute_svr( bkpnt = mup[ np.where((xdiff > (discontfiring_dur * emgfile["FSAMP"])))[0] ] + bkpnt = bkpnt[np.where(bkpnt!=mup[-1])] # Make predictions on the data if bkpnt.size > 0: # If there is a point of discontinuity From 03d037a1a3c9ab73ad932dd3077b16344d0f87ed Mon Sep 17 00:00:00 2001 From: JABeauchamp <103655220+JABeauchamp@users.noreply.github.com> Date: Sat, 15 Jun 2024 14:35:13 -0500 Subject: [PATCH 32/35] formatting --- openhdemg/library/tools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index 05836cd..58a48f1 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -1323,7 +1323,7 @@ def compute_svr( bkpnt = mup[ np.where((xdiff > (discontfiring_dur * emgfile["FSAMP"])))[0] ] - bkpnt = bkpnt[np.where(bkpnt!=mup[-1])] + bkpnt = bkpnt[np.where(bkpnt != mup[-1])] # Make predictions on the data if bkpnt.size > 0: # If there is a point of discontinuity From e85649d42a51c82d8cf4f5b20222bc5d6ea32fc8 Mon Sep 17 00:00:00 2001 From: JABeauchamp <103655220+JABeauchamp@users.noreply.github.com> Date: Mon, 17 Jun 2024 17:34:32 -0500 Subject: [PATCH 33/35] fix edge conditions for SVR fits --- openhdemg/library/tools.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index 58a48f1..c9574d1 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -1327,10 +1327,25 @@ def compute_svr( # Make predictions on the data if bkpnt.size > 0: # If there is a point of discontinuity - if bkpnt[0] == mup[0]: + if bkpnt[0] == mup[0]: # When first firing is discontinuity smoothfit = np.nan*np.ones(1) newtm = np.nan*np.ones(1) bkpnt = bkpnt[1:] + tmptm = predtime[ + 0: np.where( + (bkpnt[0] >= predind[0:-1]) & (bkpnt[0] < predind[1:]) + )[0][0] + ] # Break up time vector for first continous range of firing + smoothfit = np.append(smoothfit, svr.predict(tmptm)) + # Predict with svr model + newtm = np.append(newtm, tmptm) # Track new time vector + tmpind = predind[ + 0: np.where( + (bkpnt[0] >= predind[0:-1]) & (bkpnt[0] < predind[1:]) + )[0][0] + ] # Sample vector of first continous range of firing + # Fill corresponding sample indices with svr fit + gen_svr[tmpind.astype(np.int64)] = svr.predict(tmptm) else: tmptm = predtime[ 0: np.where( From 6e632976fdf348c2e88ac245235afeb29856be0b Mon Sep 17 00:00:00 2001 From: JABeauchamp <103655220+JABeauchamp@users.noreply.github.com> Date: Mon, 17 Jun 2024 19:04:48 -0500 Subject: [PATCH 34/35] Discontinuity bug fix --- openhdemg/library/tools.py | 21 +++------------------ 1 file changed, 3 insertions(+), 18 deletions(-) diff --git a/openhdemg/library/tools.py b/openhdemg/library/tools.py index c9574d1..7c929c8 100644 --- a/openhdemg/library/tools.py +++ b/openhdemg/library/tools.py @@ -1293,7 +1293,7 @@ def compute_svr( # Time between discharges, will use for discontinuity calc xdiff = idr[mu].diff_mupulses[1:].values # Motor unit pulses, samples - mup = np.array(idr[mu].mupulses[1:]) + mup = np.array(idr[mu].mupulses) # Defining weight vector. A scaling applied to the regularization # parameter, per sample. @@ -1313,7 +1313,7 @@ def compute_svr( # Defining prediction vector # TODO need to add custom range. # From the second firing to the end of firing, in samples. - predind = np.arange(mup[0], mup[-1]+1) + predind = np.arange(mup[1], mup[-1]+1) predtime = (predind/emgfile["FSAMP"]).reshape(-1, 1) # In time (s) # Initialise nan vector for tracking fits aligned in time. Usefull for # later quant metrics. @@ -1330,22 +1330,7 @@ def compute_svr( if bkpnt[0] == mup[0]: # When first firing is discontinuity smoothfit = np.nan*np.ones(1) newtm = np.nan*np.ones(1) - bkpnt = bkpnt[1:] - tmptm = predtime[ - 0: np.where( - (bkpnt[0] >= predind[0:-1]) & (bkpnt[0] < predind[1:]) - )[0][0] - ] # Break up time vector for first continous range of firing - smoothfit = np.append(smoothfit, svr.predict(tmptm)) - # Predict with svr model - newtm = np.append(newtm, tmptm) # Track new time vector - tmpind = predind[ - 0: np.where( - (bkpnt[0] >= predind[0:-1]) & (bkpnt[0] < predind[1:]) - )[0][0] - ] # Sample vector of first continous range of firing - # Fill corresponding sample indices with svr fit - gen_svr[tmpind.astype(np.int64)] = svr.predict(tmptm) + bkpnt = bkpnt[1:] else: tmptm = predtime[ 0: np.where( From 274d898207475c055c648f2bd03eda37f0176754 Mon Sep 17 00:00:00 2001 From: Giacomo Valli_PhD <81100252+GiacomoValliPhD@users.noreply.github.com> Date: Tue, 18 Jun 2024 13:38:18 +0200 Subject: [PATCH 35/35] Release 0.1.0 --- docs/isek_jek_tutorials.md | 2 +- setup.py | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/docs/isek_jek_tutorials.md b/docs/isek_jek_tutorials.md index c69dec6..7095818 100644 --- a/docs/isek_jek_tutorials.md +++ b/docs/isek_jek_tutorials.md @@ -26,7 +26,7 @@ You can download the sample files and the sample scripts [here](https://drive.go
-## 2023 ISEK Workshop +## 2024 ISEK Workshop **Workshop: Simplified analysis of motor unit properties with *openhdemg*.** diff --git a/setup.py b/setup.py index fdbb904..105d373 100644 --- a/setup.py +++ b/setup.py @@ -1,13 +1,12 @@ -# Copyright (C) 2023 Giacomo Valli, Paul Ritsche +# Copyright (C) 2022 - 2024. The openhdemg community. -# For the openhdemg version # To read the content of the README or description file from pathlib import Path # To install required dependencies from setuptools import setup -import openhdemg as emg +import openhdemg INSTALL_REQUIRES = [ "customtkinter==5.2.1", @@ -68,7 +67,7 @@ "Source Code": "https://github.com/GiacomoValliPhD/openhdemg", "Bug Tracker": "https://github.com/GiacomoValliPhD/openhdemg/issues", }, - version=emg.__version__, + version=openhdemg.__version__, install_requires=INSTALL_REQUIRES, include_package_data=True, packages=PACKAGES,