Updated VST Analysis #8

- Implemented solution by E.Kosari using means between peaks and troughs
- New method is enabled via 'Stick-Slip-Detection'
- Fitting for VST tests now uses percentiles instead of raw datapoints
- Customization with new options
- Added dot animation while waiting for tasks

requirements.txt

@@ -9,3 +9,4 @@ traitlets
 uncertainties
 uncertainties.egg
 pyqt5
+terminal-animation
\ No newline at end of file

rstevaluation/rstevaluation/analysis.py

 """ Advanced data analysis functions for ringshear tester data """
 import numpy as np
 import scipy.optimize as spopt
-import scipy.stats as stats
+import scipy.signal as spsig
+import scipy.stats as spstats
 from uncertainties.core import ufloat
 
 
@@ -25,8 +26,8 @@ def rst_analmut(x, y):
 
     # M = M[np.nonzero(np.isfinite(M))]
     # C = C[np.nonzero(np.isfinite(C))]
-    statsM = stats.norm.fit([d.n for d in M])
-    statsC = stats.norm.fit([d.n for d in C])
+    statsM = spstats.norm.fit([d.n for d in M])
+    statsC = spstats.norm.fit([d.n for d in C])
     uncM = ufloat(statsM[0], 2*statsM[1])
     uncC = ufloat(statsC[0], 2*statsC[1])
     fric_mut = (uncM, uncC)
@@ -80,7 +81,10 @@ def fitfric(x, a, b):
 
 
 # %%=====================VELOCITY STEPPING ANALYSIS============================
-def vst_analysis(data):
+def vst_analysis(data, config):
+    """
+    Simply fits the apparent friction with the current loading velocity.
+    """
     mu_app = data['shearstress']/data['normalstress']  # Apparent friction
     logvel = np.log10(data['velocity'])
 
@@ -91,11 +95,75 @@ def vst_analysis(data):
     onlyfinite = np.isfinite(logvel)
     mu_app = mu_app[onlyfinite]
     logvel = logvel[onlyfinite]
+    # Remove values where logvel < -4
+    mu_app = mu_app[logvel >= -4]
+    logvel = logvel[logvel >= -4]
+    pfit, pcov, x, y = vst_fit(
+        fitfric, logvel, mu_app, config
+    )
+    perr = 2*np.sqrt(np.diag(pcov))
+    name_fit = (
+        r'y = ({:.4f}$\pm${:.5f})log10(x) +' +
+        '\n\t\t' +
+        r'({:.4f}$\pm${:.5f})').format(pfit[0],
+                                       perr[0],
+                                       pfit[1],
+                                       perr[1])
+    return pfit, perr, name_fit, logvel, mu_app, x, y
 
-    pfit, pcov = spopt.curve_fit(
-        fitfric,
-        logvel,
-        mu_app
+
+# %%=====================VELOCITY STEPPING ANALYSIS============================
+def vst_analysis_alternative(data, config):
+    """
+    Alternative method for fitting apparent friction provided by E. Kosari when
+    the data shows strong stick-slip.
+    First the stick-slip events are detected with peak detection. Then the
+    peaks are fit with a log-linear equation.
+    """
+    mu_app = data['shearstress']/data['normalstress']  # Apparent friction
+    logvel = np.log10(data['velocity'])
+
+    peaks = spsig.find_peaks(
+        mu_app,
+        prominence=config.getfloat('parameters', 'peak_prominence')
+    )[0].tolist()
+    troughs = spsig.find_peaks(
+        -mu_app,
+        prominence=config.getfloat('parameters', 'peak_prominence')
+    )[0].tolist()
+
+    # If the first trough is before a peak delete it
+    while troughs[0] < peaks[0]:
+        troughs.pop(0)
+    # Remove peaks if they are not followed by a trough:
+    while peaks[-1] > troughs[-1]:
+        peaks.pop(-1)
+
+    if len(peaks) != len(troughs):
+        troughs = [
+            troughs[int(np.argwhere(p < np.array(troughs))[0])]
+            for p in peaks
+        ]
+
+    # Calculate means for fitting
+    mu_app_fit = np.array(
+        [(mu_app[p] + mu_app[t])/2 for p, t in zip(peaks, troughs)]
+    )
+    logvel_fit = logvel[peaks]
+
+    # Remove nonfinite values for fitting
+    onlyfinite = np.isfinite(mu_app_fit)
+    mu_app_fit = mu_app_fit[onlyfinite]
+    logvel_fit = logvel_fit[onlyfinite]
+    onlyfinite = np.isfinite(logvel_fit)
+    mu_app_fit = mu_app_fit[onlyfinite]
+    logvel_fit = logvel_fit[onlyfinite]
+    # Remove values where logvel < -4
+    mu_app_fit = mu_app_fit[logvel_fit >= -4]
+    logvel_fit = logvel_fit[logvel_fit >= -4]
+    # Do fitting
+    pfit, pcov, x, y = vst_fit(
+        fitfric, logvel_fit, mu_app_fit, config
     )
     perr = 2*np.sqrt(np.diag(pcov))
     name_fit = (
@@ -105,10 +173,36 @@ def vst_analysis(data):
                                        perr[0],
                                        pfit[1],
                                        perr[1])
-    return pfit, perr, name_fit
+    return pfit, perr, name_fit, logvel_fit, mu_app_fit, x, y
 
 
+def vst_fit(fitfunc, logvel, mu_app, config):
+    """
+    Takes a special subsample per velocity step.
+    """
+    vel_round = np.around(logvel, config.getint('parameters', 'vel_accuracy'))
+    vel_un = np.unique(vel_round)
+    # sample_counts = [len(np.argwhere(vel_round == v)) for v in vel_un]
+    # min_samples = np.min(sample_counts)
+    x = []
+    y = []
+    prct = np.linspace(1, 99, config.getint('parameters', 'fit_percentiles'))
+    for vu in vel_un:
+        # selection = np.argwhere(vel_round == vu)[:min_samples]
+        # x.extend(vel_round[selection])
+        # y.extend(mu_app[selection])
+        selection = mu_app[np.argwhere(vel_round == vu)]
+        for p in prct:
+            x.append(vu)
+            y.append(np.percentile(selection, p))
+    x = np.squeeze(np.array(x))
+    y = np.squeeze(np.array(y))
+    pfit, pcov = spopt.curve_fit(fitfunc, x, y)
+    return pfit, pcov, x, y
+
 # %%======================BOOTSTRAP LINEAR REGRESSION==========================
+
+
 def fit_bootstrap(datax, datay, function=None,
                   yerr_systematic=0.0, nsigma=2, nsets=100,
                   yerr=None):

rstevaluation/rstevaluation/files.py

@@ -9,6 +9,7 @@ import operator
 import os
 import shutil
 
+import animation
 import nptdms
 import numpy as np
 import uncertainties as unc
@@ -19,9 +20,9 @@ from uncertainties import unumpy as unp
 from rstevaluation import data as rstdat
 
 
+@animation.wait()
 def convert(path, file_in, config):
     """Reads data from source files and returns a dictionary with data"""
-
     if isinstance(config, configparser.ConfigParser):
         var = {
             k: json.loads(config['parameters'][k])
@@ -81,7 +82,6 @@ def convert(path, file_in, config):
             np.sqrt(10**2 + (data['shear_smth_nom']*0.005)**2)
         )
     )
-
     print(file_in+' read')
     return data
 

rstevaluation/rstevaluation/mainGUI.py

@@ -33,7 +33,7 @@ class RST_pick_GUI(tk.Tk):
         self.is_rst = tk.IntVar()
         self.cfg = configparser.ConfigParser()
         self.cfg.optionxform = str
-        self.path_cfg.set('rst_evaluation_')
+        self.path_cfg.set('')
         self.params = dict()
         self.protocol('WM_DELETE_WINDOW', self.closing_menu)
 
@@ -139,6 +139,7 @@ class RST_pick_GUI(tk.Tk):
             text='Stick-Slip Detection',
             variable=self.stick_slip)
         self.opt_stickslip.grid(column=1, row=5, sticky='EW')
+        self.opt_stickslip['state'] = tk.DISABLED
 
         # More Options
         icon = rsticons.icon_options()
@@ -244,6 +245,7 @@ class RST_pick_GUI(tk.Tk):
             self.entry_output.insert(0, f_path + '/')
         else:
             print(event.widget._name)
+        self.switch_pred_load()
 
     def start_processing(self, event):
         """Initiates Processing of data"""
@@ -255,13 +257,18 @@ class RST_pick_GUI(tk.Tk):
         if messagebox.askyesno('Close Program', 'Quit?'):
             sys.exit(0)
 
-    def switch_pred_load(self, event):
+    def switch_pred_load(self, event=None):
         """ Switches the predefined normal load to inactive for VST """
         try:
             if self.is_rst.get():
                 self.opt_pred['state'] = tk.DISABLED
+                self.opt_stickslip['state'] = tk.NORMAL
             else:
                 self.opt_pred['state'] = tk.NORMAL
+                self.opt_stickslip['state'] = tk.DISABLED
+            if not event:
+                self.opt_pred['state'] = tk.DISABLED
+                self.opt_stickslip['state'] = tk.NORMAL
         except AttributeError as _:
             pass
 
@@ -323,6 +330,7 @@ class RST_pick_GUI(tk.Tk):
         self.is_rst.set(self.cfg.getint('options', 'rst'))
         self.rev_pick.set(self.cfg.getint('options', 'rev_pick'))
         self.pred_norm.set(self.cfg.getint('options', 'pred_norm'))
+        self.stick_slip.set(self.cfg.getint('options', 'stick_slip'))
         # Read parameters and store them as a dictionary of StringVars
         # This automatically links them to changes made in the Options dialog.
         for item in self.cfg['parameters'].items():
@@ -343,6 +351,7 @@ class RST_pick_GUI(tk.Tk):
         self.cfg['options']['rst'] = str(self.is_rst.get())
         self.cfg['options']['rev_pick'] = str(self.rev_pick.get())
         self.cfg['options']['pred_norm'] = str(self.pred_norm.get())
+        self.cfg['options']['stick_slip'] = str(self.stick_slip.get())
 
         self.cfg['paths']['path_in'] = self.path_in.get()
         self.cfg['paths']['path_out'] = self.path_out.get()
@@ -370,14 +379,6 @@ class RST_pick_GUI(tk.Tk):
     def refresh_entries_from_cfg(self, path=None, init=False):
         """Refresh entries from config file"""
         if init:
-            try:
-                c = self.cfg.read(self.path_cfg.get())
-                if not(c):
-                    self.create_default_config()
-                    self.cfg.read(self.path_cfg.get())
-                    print('created default config')
-                print('read default config')
-            except Exception:
             self.create_default_config()
             self.cfg.read(self.path_cfg.get())
             print('created default config')
@@ -385,7 +386,7 @@ class RST_pick_GUI(tk.Tk):
             try:
                 c = self.cfg.read(self.path_cfg.get())
                 if not(c):
-                    self.create_default_config(path=path)
+                    self.create_default_config(path=path.get())
                     self.cfg.read(self.path_cfg.get())
                     print('created project config')
                 print('read project config')
@@ -428,6 +429,7 @@ class RST_pick_GUI(tk.Tk):
             'rev_pick': 1,    # if 1: the user can review the picked peaks
             'is_VST': 0,      # if 1: the file is a VST-Test file
             'pred_norm': 1,   # if 1: uses predefined normal stresses
+            'stick_slip': 0   # if 1: uses alternative methods for vst and rst
         }
 
         config['parameters'] = {
@@ -442,7 +444,10 @@ class RST_pick_GUI(tk.Tk):
             'smoothwindow': 5,  # Smoothing Factor for plots
             'pred_norms': [500, 1000, 2000, 4000, 8000, 16000],
             'nsamples': 10000,  # Maximum number of samples to use
-            'downsample_filter': 0,  # Use downsampling filter or not.
+            'downsample_filter': 0,  # Use downsampling filter or not
+            'peak_prominence': 0.09,  # Sensitivity for peak finding
+            'vel_accuracy': 2,  # Rounding accuracy for logarithm of velocity
+            'fit_percentiles': 20,  # Number of percentiles to take for vst fit
         }
 
         config['units'] = {
@@ -458,6 +463,9 @@ class RST_pick_GUI(tk.Tk):
             'pred_norms': 'Pa',
             'nsamples': 'samples',
             'downsample_filter': '(0 = No, 1 = Yes)',
+            'peak_prominence': 'friction difference',
+            'vel_accuracy': 'decimals',
+            'fit_percentiles': '',
         }
         self.path_cfg.set(config['paths']['path_cfg'])
         with open(self.path_cfg.get(), 'w') as f:

rstevaluation/rstevaluation/plots.py

 """ Plotting functions for rst-evaluation """
 import os
 
+import animation
 import matplotlib.gridspec as gridspec
 import matplotlib.pyplot as plt
 import numpy as np
@@ -12,6 +13,7 @@ from rstevaluation import analysis as rst_analysis
 
 
 # ===================LINEAR REGRESSION=========================================
+@animation.wait()
 def plotstd(path, name, strength, fricdata):
     """ Plots the results of the linear regression """
     plt.rcParams['savefig.format'] = 'pdf'
@@ -77,6 +79,7 @@ def plotstd(path, name, strength, fricdata):
 
 
 # %%======================PLOT HISTOGRAMS======================================
+@animation.wait()
 def plothist(path, name, strength, data_mut):
     """ Plots the result of the mutual linear regression as histograms """
     title_mu = [
@@ -185,6 +188,7 @@ def plothist(path, name, strength, data_mut):
 
 
 # %%=======================PLOT  TS============================================
+@animation.wait()
 def plotts(path, name, exp_data, normal_stress):
     """ Plots all experiments into one plot """
     plt.rcParams['figure.figsize'] = (10, 8)
@@ -238,7 +242,8 @@ def plotts(path, name, exp_data, normal_stress):
 
 
 # %%=================PLOT VST ANALYSIS=========================================
-def plotVST(path, data, pfit, perr, name_fit):
+@animation.wait()
+def plotVST(path, data, pfit, perr, name_fit, x, y, xq, yq):
     """ Plots VST analysis data """
     plt.rcParams['savefig.format'] = 'pdf'
     plt.rcParams['font.size'] = 10
@@ -276,9 +281,10 @@ def plotVST(path, data, pfit, perr, name_fit):
     ax2.set_xlim(0, np.nanmax(displ) + np.nanmax(displ)/100)
 
     # shear velocity vs. fricton:
-    ax3.plot(vel, fric, 'k.', ms=2, rasterized=True)
+    ax3.plot(10**x, y, 'k.', ms=2, rasterized=True, label='original data')
+    ax3.plot(10**xq, yq, 'r.', ms=2, rasterized=True, label='fit data')
     ax3.plot(
-        vel, rst_analysis.fitfric(np.log10(vel), *pfit),
+        10**x, rst_analysis.fitfric(x, *pfit),
         'g-',
         label='curve fit: ' + name_fit
     )
@@ -286,13 +292,13 @@ def plotVST(path, data, pfit, perr, name_fit):
     ax3.set_xlabel(r'Shear velocity $v$ [$\mathregular{mm s^{-1}}$]')
     ax3.set_ylabel('Friction')
     ax3.set_xlim(
-        np.nanmin(vel)-(np.nanmin(vel)/5),
-        np.nanmax(vel)+np.nanmax(vel)/5
+        np.nanmin(10**x)-(np.nanmin(10**x)/5),
+        np.nanmax(10**x)+np.nanmax(10**x)/5
     )
     ax3.legend(fontsize=8,
                facecolor='w',
                edgecolor='k',
-               loc='upper right',
+               #    loc='upper right',
                framealpha=0.5)
     fname = os.path.splitext(data['name'])[0]
     fig.suptitle(fname, fontsize=14, y=0.95)

rstevaluation/rstevaluation/processing.py

@@ -5,6 +5,8 @@ import os
 import subprocess
 import sys
 
+from tqdm import tqdm
+
 from rstevaluation import analysis as rstanalysis
 from rstevaluation import files as rstfiles
 from rstevaluation import picking as rstpicking
@@ -26,11 +28,17 @@ def evaluate(config):
         evaluate_data(eval_data, exp_data, config)
     else:
         print('Velocity Stepping Test')
-        for exp in exp_data:
-            pfit, perr, name_fit = rstanalysis.vst_analysis(exp)
+        for exp in tqdm(exp_data, desc='Analysing', leave=False):
+            if config.getint('options', 'stick_slip'):
+                pf, pe, nfit, x, y, xq, yq = rstanalysis.vst_analysis_alternative(
+                    exp, config)
+            else:
+                pf, pe, nfit, x, y, xq, yq = rstanalysis.vst_analysis(
+                    exp, config)
             rstplots.plotVST(
                 config['paths']['path_out'],
-                exp, pfit, perr, name_fit)
+                exp, pf, pe, nfit,
+                x, y, xq, yq)
     with open(config['paths']['path_cfg'], 'w') as cfg_file:
         config.write(cfg_file)
     open_explorer(config['paths']['path_out'])
@@ -53,6 +61,7 @@ def get_data(config):
         rstfiles.check_tdms(file_list, config)
 
     # Data is stored as a list of dictionaries containing the data
+    print('Loading data files')
     exp_data = [
         rstfiles.convert(config['paths']['path_in'], f, config)
         for f in file_list

setup.cfg

@@ -20,6 +20,7 @@ install_requires =
     tqdm
     uncertainties
     pyqt5
+    terminal-animation
 include_package_data = True
 
 [options.packages.find]