First attempt at Pull request from Alex B

PYMEcs/Analysis/MINFLUX.py

@@ -1,47 +1,13 @@
-from scipy.stats import binned_statistic
-from PYMEcs.IO.MINFLUX import get_stddev_property
-import numpy as np
-import matplotlib.pyplot as plt
-from PYMEcs.pyme_warnings import warn
+import os
 
+import matplotlib.pyplot as plt
+import numpy as np
 import pandas as pd
-from PYMEcs.misc.matplotlib import boxswarmplot
-
-def plot_stats(ds,ax,errdict,sdmax=None,swarmsize=3,siteKey='siteID'):
-    df = site_stats(ds,errdict,siteKey=siteKey)
-    kwargs = dict(swarmsize=swarmsize,width=0.2,annotate_means=True,annotate_medians=True,swarmalpha=0.4)
-    boxswarmplot(df,ax=ax,**kwargs)
-    ax.set_ylim(0,sdmax)
-    ax.set_ylabel('precision [nm]')
-    return df
-
-def site_stats(ds,sitedict,siteKey='siteID'):
-    uids, idx = np.unique(ds[siteKey],return_index=True)
-    sitestats = {}
-    for key in sitedict:
-        prop = sitedict[key]
-        sitestats[key] = ds[prop][idx]
-    df = pd.DataFrame.from_dict(sitestats)
-    return df[df > 0].dropna() # this drops rows with zeroes; these should not occur but apparently do; probably a bug somewhere
+from scipy.stats import binned_statistic
 
-def plotsitestats(p,origamiErrorLimit=10,figsize=None,swarmsize=3,siteKey='siteID',fignum=None):
-    uids = np.unique(p[siteKey])
-    fig, axs = plt.subplots(2, 2, figsize=figsize,num=fignum)
-    plot_stats(p,axs[0, 0],dict(xd_sd_corr='error_x',x_sd='error_x_nc',x_sd_trace='error_x_ori'),
-                sdmax=origamiErrorLimit,swarmsize=swarmsize,siteKey=siteKey)
-    plot_stats(p,axs[0, 1],dict(yd_sd_corr='error_y',y_sd='error_y_nc',y_sd_trace='error_y_ori'),
-                sdmax=origamiErrorLimit,swarmsize=swarmsize,siteKey=siteKey)
-    if p.mdh.get('MINFLUX.Is3D'):
-        plot_stats(p,axs[1, 0],dict(zd_sd_corr='error_z',z_sd='error_z_nc',z_sd_trace='error_z_ori'),
-                    sdmax=origamiErrorLimit,swarmsize=swarmsize,siteKey=siteKey)
+from PYMEcs.IO.MINFLUX import get_stddev_property
+from PYMEcs.pyme_warnings import warn
 
-    all_axes = dict(xd_sd_corr='error_x',yd_sd_corr='error_y')
-    if p.mdh.get('MINFLUX.Is3D'):
-        all_axes['zd_sd_corr'] = 'error_z'
-    df_allaxes = plot_stats(p,axs[1, 1],all_axes,sdmax=origamiErrorLimit,swarmsize=swarmsize,siteKey=siteKey)
-    fig.suptitle('Site stats for %d sites' % df_allaxes.shape[0])
-    plt.tight_layout()
-    return df_allaxes
 
 def propcheck_density_stats(ds,warning=True):
     for prop in ['clst_area','clst_vol','clst_density','clst_stdz']:
@@ -90,6 +56,9 @@ def plot_density_stats(ds,objectID='dbscanClumpID',scatter=False):
         ax1[1].scattered_boxplot(sz,labels=['Stddev Z'],showmeans=True)
     plt.tight_layout()
 
+from PYMEcs.misc.matplotlib import boxswarmplot
+
+
 def plot_density_stats_sns(ds,objectID='dbscanClumpID'):
     if not propcheck_density_stats(ds):
         return
@@ -113,16 +82,46 @@ def plot_density_stats_sns(ds,objectID='dbscanClumpID'):
     fig.suptitle('Density stats for %d clusters' % dens.size)
     plt.tight_layout()
 
+    # --- Save the values to csv --- (Alex B addition)
+    df = pd.DataFrame({
+        "Metric": ["cluster density"],
+        "Mean": [np.mean(dens)],
+        "Median": [np.median(dens)],
+        "Unit": ["#/um^2"]})
+    # --- Show the head of df in the console --- (Alex B addition)
+    print(df.head())
+
+    # --- Save as a csv file --- (Alex B addition)
+    # Save the cluster size file
+    #fn = os.path.splitext(ds.mdh.get('MINFLUX.Filename'))[0]
+    dirpath, filename = os.path.split(ds.filename)
+    fn = ds.mdh.get('MINFLUX.TimeStamp')
+    df.to_csv(os.path.join(dirpath, fn + "_clusterDensity.csv"), 
+              index=False, header=True)
+
+    # save data CSV file with dbscanClustered
+    df_dbscanClustered = ds.dataSources['dbscanClustered'].to_pandas()
+    df_dbscanClustered.to_csv(
+        os.path.join(dirpath, fn + "_dbscanClustered_data.csv"), index=False, header=True)
+    # Used as a reminder for LocRate csv saving
+    # print(
+    #   f'\ncsv name is: {fn + "_clusterDensity.csv"}\nIf you did not load a session, csv file and figures will be saved on the desktop')
+    # By default the file is saved on the Desktop, if a session file is used, it is saved in the same directory as the session file.
+
     return dens
 
 def plot_stats_minflux(deltas, durations, tdintrace, efo_or_dtovertime, times,
-                       showTimeAverages=False, dsKey=None, areaString=None):
+                       showTimeAverages=False, dsKey=None, areaString=None, timestamp=None, dirPath=None):
     from scipy.stats import iqr
-
+
+    # --- Create the figure (plot with 2x2 subplots) ---
     fig, (ax1, ax2) = plt.subplots(2, 2)
+
+    # --- Compute some statistics for TBT plot (top-right = ax1[0]) ---
     dtmedian = np.median(deltas)
     dtmean = np.mean(deltas)
     dtiqr = iqr(deltas,rng=(10, 90)) # we are going for the 10 to 90 % range
+
     h = ax1[0].hist(deltas,bins=40,range=(0,dtmean + 2*dtiqr))
     ax1[0].plot([dtmedian,dtmedian],[0,h[0].max()])
     # this is time between one dye molecule and the next dye molecule being seen
@@ -131,14 +130,16 @@ def plot_stats_minflux(deltas, durations, tdintrace, efo_or_dtovertime, times,
              verticalalignment='bottom', transform=ax1[0].transAxes)
     ax1[0].text(0.95, 0.7, '  mean %.2f s' % dtmean, horizontalalignment='right',
              verticalalignment='bottom', transform=ax1[0].transAxes)
-    if not areaString is None:
+    if areaString is not None:
         ax1[0].text(0.95, 0.6, areaString, horizontalalignment='right',
              verticalalignment='bottom', transform=ax1[0].transAxes)
 
+    # --- Compute some statistics for trace duration plot (top-left = ax1[1]) ---
     durmedian = np.median(durations)
     durmean = np.mean(durations)
     duriqr = iqr(durations,rng=(10, 90))
     h = ax1[1].hist(durations,bins=40,range=(0,durmean + 2*duriqr))
+
     ax1[1].plot([durmedian,durmedian],[0,h[0].max()])
     ax1[1].set_xlabel('duration of "traces" [s]')
     ax1[1].text(0.95, 0.8, 'median %.0f ms' % (1e3*durmedian), horizontalalignment='right',
@@ -147,11 +148,13 @@ def plot_stats_minflux(deltas, durations, tdintrace, efo_or_dtovertime, times,
              verticalalignment='bottom', transform=ax1[1].transAxes)
     # ax1[1].set_xlim(0,durmean + 2*duriqr) # superfluous since we are using the range keyword in hist
 
+    # --- Compute some statistics for time between localisations in same trace (bottom-left = ax2[0]) ---
     tdintrace_ms = 1e3*tdintrace
     tdmedian = np.median(tdintrace_ms)
     tdmean = np.mean(tdintrace_ms)
     tdiqr = iqr(tdintrace_ms,rng=(10, 90))
     h = ax2[0].hist(tdintrace_ms,bins=50,range=(0,tdmean + 2*tdiqr))
+
     ax2[0].plot([tdmedian,tdmedian],[0,h[0].max()])
     # these are times between repeated localisations of the same dye molecule
     ax2[0].set_xlabel('time between localisations in same trace [ms]')
@@ -161,6 +164,7 @@ def plot_stats_minflux(deltas, durations, tdintrace, efo_or_dtovertime, times,
              verticalalignment='bottom', transform=ax2[0].transAxes)
     # ax2[0].set_xlim(0,tdmean + 2*tdiqr) # superfluous since we are using the range keyword in hist
 
+    # --- Compute the TBT running time average (bottom-right = ax2[1]) ---
     if showTimeAverages:
         ax2[1].plot(times,efo_or_dtovertime)
         ax2[1].set_xlabel('TBT running time average [s]')
@@ -175,7 +179,30 @@ def plot_stats_minflux(deltas, durations, tdintrace, efo_or_dtovertime, times,
         plt.suptitle('Location rate analysis from datasource %s' % dsKey)
     plt.tight_layout()
 
+    # --- Calculate dimensions and area of the image ---
+    dimension = areaString.split(' ')[1]
+    area = float(dimension.split('x')[0])*float(dimension.split('x')[1])
+
+    # --- Save medians values to csv --- (Alex B addition)
+    df = pd.DataFrame({
+        "Metric": ["TBT (Time Between Traces)", "ROI dimension", "ROI area", "Trace Duration", "Time Between Localizations"],
+        "Median": [dtmedian.round(2), dimension, area, (durmedian * 1e3).round(0), tdmedian.round(0)],  # Convert seconds -> milliseconds where needed
+        "Mean": [dtmean.round(2), dimension, area, (durmean * 1e3).round(0), tdmean.round(0)],  # Convert seconds -> milliseconds where needed
+        "IQR": [dtiqr.round(2), dimension, area, (duriqr * 1e3).round(0), tdiqr.round(0)],  # Convert seconds -> milliseconds where needed
+        "Unit": ["s","um^2","um^2", "ms", "ms"]    })
+
+    # --- Show the head of df in the console --- (Alex B addition)
+    print(df.head())
+
+    # --- Save as a csv file --- (Alex B addition)
+    # Save the LocRate file
+    df.to_csv(os.path.join(dirPath, timestamp + "_locRate.csv"), index=False, header=True)
 
+    # --- Save the figure --- (Alex B addition)
+    plt.savefig(os.path.join(dirPath, timestamp + '_locRate.png'),
+                dpi=300, bbox_inches='tight')
+    print(f'Figure and *.csv file saved to {os.path.join(dirPath, timestamp + "_locRate.png")}')
+
 # this function assumes a pandas dataframe
 # the pandas frame should generally be generated via the function minflux_npy2pyme from PYMEcs.IO.MINFLUX
 def analyse_locrate_pdframe(datain,use_invalid=False,showTimeAverages=True):
@@ -214,7 +241,10 @@ def analyse_locrate_pdframe(datain,use_invalid=False,showTimeAverages=True):
     tdiff = data['tim'].values[1:]-data['tim'].values[:-1]
     tdsmall = tdiff[tdiff <= 0.1]
     tdmedian = np.median(tdsmall)
-
+
+    dirpath, filename = os.path.split(datain.filename)
+    timeStamp = datain.mdh.get('MINFLUX.TimeStamp')
+
     start_times = data.loc[startindex,'tim'][:-1].values # we use those for binning the deltas, we discard final time to match size of deltas
 
     if has_invalid and use_invalid:
@@ -226,13 +256,16 @@ def analyse_locrate_pdframe(datain,use_invalid=False,showTimeAverages=True):
     if showTimeAverages:
         delta_averages, bin_edges, binnumber = binned_statistic(start_times,deltas,statistic='mean', bins=50)
         delta_av_times = 0.5*(bin_edges[:-1] + bin_edges[1:]) # bin centres
-        plot_stats_minflux(deltas, durations_proper, tdiff, tdmedian, delta_averages, delta_av_times, showTimeAverages=True)
+        plot_stats_minflux(deltas, durations_proper, tdiff, tdmedian, 
+                           delta_averages, delta_av_times, showTimeAverages=True, 
+                           timestamp=timeStamp, dirPath=dirpath)
     else:
-        plot_stats_minflux(deltas, durations_proper, tdiff, tdmedian, data['efo'], None)
+        plot_stats_minflux(deltas, durations_proper, tdiff, tdmedian, data['efo'], None, 
+                           timestamp=timeStamp, dirPath=dirpath)
 
 
 # similar version but now using a pipeline
-def analyse_locrate(data,datasource='Localizations',showTimeAverages=True, plot=True):
+def analyse_locrate(data,datasource='Localizations',showTimeAverages=True, plot=True, timestamp=None):
     curds = data.selectedDataSourceKey
     data.selectDataSource(datasource)
     bins = np.arange(int(data['clumpIndex'].max())+1) + 0.5
@@ -258,14 +291,17 @@ def analyse_locrate(data,datasource='Localizations',showTimeAverages=True, plot=
     leny_um = 1e-3*(data['y'].max()-data['y'].min())
     area_string = 'area %.1fx%.1f um^2' % (lenx_um,leny_um)
     data.selectDataSource(curds)
-
+
+    dirpath, filename = os.path.split(data.filename)
+    timeStamp = data.mdh.get('MINFLUX.TimeStamp')
+
     if plot:
         if showTimeAverages:
             delta_averages, bin_edges, binnumber = binned_statistic(starts[:-1],deltas,statistic='mean', bins=50)
             delta_av_times = 0.5*(bin_edges[:-1] + bin_edges[1:]) # bin centres
             plot_stats_minflux(deltas, durations_proper, tdintrace, delta_averages, delta_av_times,
-                               showTimeAverages=True, dsKey = datasource, areaString=area_string)
+                               showTimeAverages=True, dsKey = datasource, areaString=area_string, timestamp=timeStamp, dirPath=dirpath)
         else:
-            plot_stats_minflux(deltas, durations_proper, tdintrace, data['efo'], None, dsKey = datasource, areaString=area_string)
+            plot_stats_minflux(deltas, durations_proper, tdintrace, data['efo'], None, dsKey = datasource, areaString=area_string, timestamp=timeStamp, dirPath=dirpath)
 
     return (starts,ends,deltas,durations_proper,tdintrace)

PYMEcs/Analysis/NPC.py

@@ -1,5 +1,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
+
 from PYMEcs.pyme_warnings import warn
 
 piover4 = np.pi/4.0
@@ -31,6 +32,8 @@
 #     return (crot.T[:,0],crot.T[:,1])
 
 from circle_fit import taubinSVD
+
+
 def fitcirc(x,y,sigma=None):
     pcs = np.vstack((x,y)).T
     xc, yc, r, sigma = taubinSVD(pcs)
@@ -222,8 +225,9 @@ def segnotorad(sec):
     else:
         return Nlabeled
 
-from scipy.special import binom
 from scipy.optimize import curve_fit
+from scipy.special import binom
+
 
 def pn(k,n,p):
     return (binom(n,k)*(np.power(p,k)*np.power((1-p),(n-k))))
@@ -328,9 +332,10 @@ def to3vecs(x,y,z):
 def xyzfrom3vec(v):
     return (v[:,0],v[:,1],v[:,2])
 
-from scipy.spatial.transform import Rotation as R
 from scipy.interpolate import RegularGridInterpolator
 from scipy.signal import fftconvolve
+from scipy.spatial.transform import Rotation as R
+
 
 def fpinterpolate(fp3d,x,y,z,method='linear', bounds_error=True, fill_value=np.nan):
     # V[i,j,k] = 100*x[i] + 10*y[j] + z[k]
@@ -617,7 +622,7 @@ def normalize_points(self,zclip=None,mode='mean'):
             raise RuntimeError("unknown mode '%s', should be mean or median" % mode)
         npts = self.points - self.offset
         nt = self.t
-        if not zclip is None:
+        if zclip is not None:
             zgood = (npts[:,2] > -zclip)*(npts[:,2] < zclip)
             npts = npts[zgood,:]
             nt = nt[zgood]
@@ -881,7 +886,6 @@ def measure_labeleff(self,nthresh=1,do_plot=False,printpars=False,refit=False):
             self.measurements.append([nt,nb])
 
     def plot_labeleff(self,thresh=None):
-        from PYMEcs.misc.utils import get_timestamp_from_filename
         if len(self.measurements) < 10:
             raise RuntimeError("not enough measurements, need at least 10, got %d" %
                                len(self.measurements))
@@ -895,6 +899,26 @@ def plot_labeleff(self,thresh=None):
         plt.figure()
         plotcdf_npc3d(nlab,timestamp=get_timestamp_from_filename(self.filename),thresh=thresh)
 
+        # --- Alex B addiiton for auto-saving of LE-figure ---
+        # Function is called in NPCcalcLM.py (PYMEcs>experimental>OnAnalyse3DNPCsByID_auto_save)
+
+    def plot_labeleff_for_auto_save(self,thresh=None):
+        if len(self.measurements) < 10:
+            raise RuntimeError("not enough measurements, need at least 10, got %d" %
+                               len(self.measurements))
+        meas = np.array(self.measurements)
+        nlab = meas.sum(axis=1)
+        # fill with trailing zeros if we have a known number of NPCs but have fewer measurements
+        # the "missing NPCs" typically represent NPCs with no events
+        if int(self.known_number) > 0 and nlab.shape[0] < int(self.known_number):
+            nlab = np.pad(nlab,((0,int(self.known_number)-nlab.shape[0])))
+
+        LEfig = plt.figure()
+        plotcdf_npc3d(nlab,timestamp=get_timestamp_from_filename(self.filename),thresh=thresh)
+        return LEfig
+
+        # --- Fin Alex B addition ---
+
     def diam(self):
         diams = []
         for npc in self.npcs: