Merge pull request SarthakJariwala#19 from SarthakJariwala/development

Improvements to Lifetime Analysis Code

Author: SarthakJariwala

PythonGUI_apps/DataBrowser.py

@@ -85,8 +85,10 @@ def load_app(self):
 
 
 def run():
+    app = QtGui.QApplication(sys.argv)#.instance()
+    app.setStyle("Fusion")
     win = MainWindow()
-    QtGui.QApplication.instance().exec_()
+    sys.exit(app.exec_())
     return
 
 run()

PythonGUI_apps/Lifetime_analysis/Fit_functions.py

@@ -9,10 +9,10 @@
 from scipy.optimize import differential_evolution
 from scipy.special import gamma
 
-def stretch_exp_fit(TRPL, t, Tc = (0,1e5), Beta = (0,1), A = (0,1e6)):
+def stretch_exp_fit(TRPL, t, Tc = (0,1e5), Beta = (0,1), A = (0,1e6), noise=(0,1e6)):
 
-    def exp_stretch(t, tc, beta, a):
-        return (a * np.exp(-((1.0 / tc) * t) ** beta)) 
+    def exp_stretch(t, tc, beta, a, noise):
+        return ((a * np.exp(-((1.0 / tc) * t) ** beta)) + noise) 
 
     def avg_tau_from_exp_stretch(tc, beta):
         return (tc / beta) * gamma(1.0 / beta)
@@ -25,13 +25,14 @@ def residuals(params):#params are the parameters to be adjusted by differential
             tc = params[0]
             beta = params[1]
             a = params[2]
+            noise = params[3]
 
-            PL_sim = exp_stretch(t,tc,beta,a)
+            PL_sim = exp_stretch(t,tc,beta,a, noise)
 
-            Resid= np.sqrt(np.sum(((PL_sim-TRPL)**2)/(np.sqrt(PL_sim)**2)))
+            Resid= np.sum(((PL_sim-TRPL)**2)/(np.sqrt(TRPL)**2))
             return Resid #returns the difference between the PL data and simulated data
 
-        bounds = [Tc, Beta, A] 
+        bounds = [Tc, Beta, A, noise] 
 
         result = differential_evolution(residuals, bounds)
         return result.x
@@ -41,20 +42,21 @@ def residuals(params):#params are the parameters to be adjusted by differential
     tc = p[0]
     beta = p[1]
     a = p[2]
+    noise = p[3]
 
-    PL_fit = exp_stretch(t,tc,beta,a)
+    PL_fit = exp_stretch(t,tc,beta,a, noise)
 
     avg_tau = avg_tau_from_exp_stretch(tc,beta)
 
-    return tc, beta, a, avg_tau, PL_fit
+    return tc, beta, a, avg_tau, PL_fit, noise
 
-def double_exp_fit(TRPL, t, tau1_bounds=(0,100), a1_bounds=(0,1e5), tau2_bounds=(0,100), a2_bounds=(0,1e5)):
+def double_exp_fit(TRPL, t, tau1_bounds=(0,1000), a1_bounds=(0,1e6), tau2_bounds=(0,10000), a2_bounds=(0,1e5), noise=(0,1e6)):
 
     def single_exp(t, tau, a):
-        return (a * np.exp(-((1.0 / tau)*t) ))
+        return (a * np.exp(-((1.0 / tau)*t)))
 
-    def double_exp(t, tau1, a1, tau2, a2):
-        return ((single_exp(t, tau1, a1)) + (single_exp(t, tau2, a2)))
+    def double_exp(t, tau1, a1, tau2, a2, noise):
+        return ((single_exp(t, tau1, a1)) + (single_exp(t, tau2, a2)) + noise)
 
     def avg_tau_from_double_exp(tau1, a1, tau2, a2):
         return (((tau1*a1) + (tau2*a2))/(a1+a2))
@@ -68,13 +70,14 @@ def residuals(params):#params are the parameters to be adjusted by differential
             a1 = params[1]
             tau2 = params[2]
             a2 = params[3]
+            noise = params[4]
 
-            PL_sim = double_exp(t,tau1, a1, tau2, a2)
+            PL_sim = double_exp(t,tau1, a1, tau2, a2, noise)
 
-            Resid= np.sqrt(np.sum(((PL_sim-TRPL)**2)/(np.sqrt(PL_sim)**2)))
+            Resid= np.sum(((PL_sim-TRPL)**2)/(np.sqrt(TRPL)**2))
             return Resid #returns the difference between the PL data and simulated data
 
-        bounds = [tau1_bounds, a1_bounds, tau2_bounds, a2_bounds] 
+        bounds = [tau1_bounds, a1_bounds, tau2_bounds, a2_bounds, noise] 
 
         result = differential_evolution(residuals, bounds)
         return result.x
@@ -85,17 +88,18 @@ def residuals(params):#params are the parameters to be adjusted by differential
     a1 = p[1]
     tau2 = p[2]
     a2 = p[3]
+    noise = p[4]
 
-    PL_fit = double_exp(t, tau1, a1, tau2, a2)
+    PL_fit = double_exp(t, tau1, a1, tau2, a2, noise)
 
     avg_tau = avg_tau_from_double_exp(tau1, a1, tau2, a2)
 
-    return tau1, a1, tau2, a2, avg_tau, PL_fit
+    return tau1, a1, tau2, a2, avg_tau, PL_fit, noise
 
-def single_exp_fit(TRPL, t, tau_bounds=(0,1000), a_bounds=(0,1e5)):
+def single_exp_fit(TRPL, t, tau_bounds=(0,10000), a_bounds=(0,1e6), noise=(0,1e6)):
 
-    def single_exp(t, tau, a):
-        return (a * np.exp(-((1.0 / tau)*t) ))
+    def single_exp(t, tau, a, noise):
+        return (a * np.exp(-((1.0 / tau)*t) ) + noise)
 
     def Diff_Ev_Fit_singleExp(TRPL):
         TRPL = TRPL
@@ -104,13 +108,14 @@ def residuals(params):#params are the parameters to be adjusted by differential
             #Variable Rates
             tau = params[0]
             a = params[1]
+            noise = params[2]
 
-            PL_sim = single_exp(t, tau, a)
+            PL_sim = single_exp(t, tau, a, noise)
 
-            Resid= np.sqrt(np.sum(((PL_sim-TRPL)**2)/(np.sqrt(PL_sim)**2)))
+            Resid= np.sum(((PL_sim-TRPL)**2)/(np.sqrt(TRPL)**2))
             return Resid #returns the difference between the PL data and simulated data
 
-        bounds = [tau_bounds, a_bounds] 
+        bounds = [tau_bounds, a_bounds, noise] 
 
         result = differential_evolution(residuals, bounds)
         return result.x
@@ -119,8 +124,9 @@ def residuals(params):#params are the parameters to be adjusted by differential
 
     tau = p[0]
     a = p[1]
+    noise = p[2]
 
-    PL_fit = single_exp(t, tau, a)
+    PL_fit = single_exp(t, tau, a, noise)
 
-    return tau, a, PL_fit
+    return tau, a, PL_fit, noise
 

PythonGUI_apps/Lifetime_analysis/Lifetime_plot_fit.py

@@ -25,12 +25,12 @@
 # local module imports
 try:
     from Lifetime_analysis.Fit_functions import stretch_exp_fit, double_exp_fit, single_exp_fit
-    from Lifetime_analysis.picoharp_phd import read_picoharp_phd
+    from Lifetime_analysis.read_ph_phd import read_picoharp_phd
     from Lifetime_analysis.Fit_functions_with_irf import fit_exp_stretch_diffev, fit_exp_stretch_fmin_tnc, fit_multi_exp_diffev, fit_multi_exp_fmin_tnc
 except:
     from Fit_functions import stretch_exp_fit, double_exp_fit, single_exp_fit
     from Fit_functions_with_irf import fit_exp_stretch_diffev, fit_exp_stretch_fmin_tnc, fit_multi_exp_diffev, fit_multi_exp_fmin_tnc
-    from picoharp_phd import read_picoharp_phd
+    from read_ph_phd import read_picoharp_phd
 
 """Recylce params for plotting"""
 plt.rc('xtick', labelsize = 20)
@@ -95,10 +95,12 @@ def __init__(self):
         self.file = None
         self.out = None # output file after fitting
         self.data_list = []
-        self.fit_lifetime_called = False
+        self.fit_lifetime_called_w_irf = False
+        self.fit_lifetime_called_wo_irf = False
         self.x_mem = [] # containers for adding x data to memory
         self.y_mem = [] # containers for adding y data to memory
         self.best_fit_mem = [] # containers for adding best fit data to memory
+        self.best_fit_mem_x = [] # containers for adding best fit data to memory
         self.legend = [] # containers for adding legend to memory
 
         #variables accounting for data received from FLIM analysis
@@ -127,6 +129,10 @@ def open_with_skip_rows_window(self):
         skip_rows = self.skip_rows_window.ui.skip_rows_spinBox.value()
         if ".txt" in self.filename[0]:
             self.file = np.loadtxt(self.filename[0], skiprows=skip_rows)
+            
+            if self.file.ndim == 1: # if there is only one trace, reshape to 2D
+                self.file = self.file.reshape(self.file.shape[0], 1)
+    
         elif ".csv" in self.filename[0]:
             self.file = np.genfromtxt(self.filename[0], skip_header=skip_rows, delimiter=",")
 
@@ -147,6 +153,10 @@ def open_irf_with_skip_rows_window(self):
         irf_skip_rows = self.irf_skip_rows_window.ui.skip_rows_spinBox.value()
         if ".txt" in self.irf_filename[0]:
             self.irf_file = np.loadtxt(self.irf_filename[0], skiprows=irf_skip_rows)
+            
+            if self.irf_file.ndim == 1: # if there is only one trace, reshape to 2d array
+                self.irf_file = self.irf_file.reshape(self.irf_file.shape[0], 1)
+                        
         elif ".csv" in self.irf_filename[0]:
             self.irf_file = np.genfrontxt(self.irf_filename[0], skip_header=irf_skip_rows, delimiter=",")
 
@@ -258,7 +268,8 @@ def plot(self):
                 x,y = self.acquire_settings() #get data
 
             self.ui.plot.plot(x, y, clear=self.ui.clear_plot_checkBox.isChecked(), pen=pg.mkPen(self.plot_color))
-            self.fit_lifetime_called = False
+            self.fit_lifetime_called_w_irf = False
+            self.fit_lifetime_called_wo_irf = False
 
             try:
                 self.ui.Result_textBrowser.setText("Integral Counts :\n" "{:.2E}".format(
@@ -303,7 +314,7 @@ def fit_and_plot(self):
                 self.ui.plot.plot(t, y, clear=self.ui.clear_plot_checkBox.isChecked(), pen=pg.mkPen(self.plot_color))
 
                 if fit_func == "Stretched Exponential": #stretch exponential tab
-                    tc, beta, a, avg_tau, PL_fit = stretch_exp_fit(TRPL_interp, t)
+                    tc, beta, a, avg_tau, PL_fit, noise = stretch_exp_fit(TRPL_interp, t)
                     self.out = np.empty((len(t), 3))
                     self.out[:,0] = t #time
                     self.out[:,1] = TRPL_interp #Raw PL 
@@ -313,11 +324,12 @@ def fit_and_plot(self):
                                                        "\nFit Method: " + "diff_ev" + #TODO : change when diff_ev and fmin_tnc implemented for non-irf
                                                        "\nAverage Lifetime = " + str(avg_tau)+ " ns"
                                                        "\nCharacteristic Tau = " + str(tc)+" ns"
-                                                       "\nBeta = "+str(beta))
+                                                       "\nBeta = "+str(beta)+
+                                                       "\nNoise = "+ str(noise))
                     self.ui.average_lifetime_spinBox.setValue(avg_tau)
 
                 elif fit_func == "Double Exponential": #double exponential tab
-                    tau1, a1, tau2, a2, avg_tau, PL_fit = double_exp_fit(TRPL_interp, t)
+                    tau1, a1, tau2, a2, avg_tau, PL_fit, noise = double_exp_fit(TRPL_interp, t)
                     self.out = np.empty((len(t), 3))
                     self.out[:,0] = t #time
                     self.out[:,1] = TRPL_interp #Raw PL 
@@ -329,11 +341,12 @@ def fit_and_plot(self):
                                                        "\nTau 1 = " + str(tau1)+" ns"
                                                        "\nA 1 = " + str(a1)+
                                                        "\nTau 2 = " + str(tau2)+" ns"
-                                                       "\nA 2 = " + str(a2))
+                                                       "\nA 2 = " + str(a2)+
+                                                       "\nNoise = "+ str(noise))
                     #TODO - once tau_avg implemented, set average lifetime spinbox to tau_avg value
 
                 elif fit_func == "Single Exponential": #single exponential tab
-                    tau, a, PL_fit = single_exp_fit(TRPL_interp, t)
+                    tau, a, PL_fit, noise = single_exp_fit(TRPL_interp, t)
                     self.out = np.empty((len(t), 3))
                     self.out[:,0] = t #time
                     self.out[:,1] = TRPL_interp #Raw PL 
@@ -342,12 +355,14 @@ def fit_and_plot(self):
                     self.ui.Result_textBrowser.setText("Fit Results:\n\nFit Function: Single Exponential"
                                                        "\nFit Method: " + "diff_ev" +
                                                        "\nLifetime = " + str(tau)+ " ns"
-                                                       "\nA = " + str(a))
+                                                       "\nA = " + str(a)+
+                                                       "\nNoise = "+ str(noise))
                     self.ui.average_lifetime_spinBox.setValue(tau)
 
                 #add fit params to data_list
                 self.data_list.append("Data Channel: " + str(self.ui.Data_channel_spinBox.value()) + "\n" + self.ui.Result_textBrowser.toPlainText())
-                self.fit_lifetime_called = True
+                self.fit_lifetime_called_wo_irf = True
+                self.fit_lifetime_called_w_irf = False
 
                 self.ui.plot.setLabel('left', 'Intensity', units='a.u.')
                 self.ui.plot.setLabel('bottom', 'Time (ns)')
@@ -475,7 +490,8 @@ def fit_and_plot_with_irf(self):
 
                 #add fit params to data_list
                 self.data_list.append("Data Channel: " + str(self.ui.Data_channel_spinBox.value()) + "\n" + self.ui.Result_textBrowser.toPlainText())
-                self.fit_lifetime_called = True
+                self.fit_lifetime_called_w_irf = True
+                self.fit_lifetime_called_wo_irf = False
         except Exception as e:
             self.ui.Result_textBrowser.append(format(e))
 
@@ -548,12 +564,19 @@ def clean_up_after_fig_export(self):
         self.y_mem = []
         self.legend = []
         self.best_fit_mem = []
+        self.best_fit_mem_x = []
 
     def add_trace_to_mem(self):
         try:
-            if self.fit_lifetime_called == True:
+            if self.fit_lifetime_called_w_irf == True:
                 self.x_mem.append(self.out[:,0])
                 self.y_mem.append(self.out[:,1])
+                self.best_fit_mem_x.append(self.out[:,0])
+                self.best_fit_mem.append(self.out[:,2])
+            elif self.fit_lifetime_called_wo_irf == True:
+                self.x_mem.append(self.acquire_settings()[0])
+                self.y_mem.append(self.acquire_settings()[1])
+                self.best_fit_mem_x.append(self.out[:,0])
                 self.best_fit_mem.append(self.out[:,2])
             else:
                 self.x_mem.append(self.acquire_settings()[0])
@@ -578,8 +601,8 @@ def pub_ready_plot_export(self):
                 plt.tick_params(direction='out', length=8, width=3.5)
                 for i in range(len(self.x_mem)):
                     plt.plot(self.x_mem[i], self.y_mem[i], label=str(self.legend[i]))
-                    if self.fit_lifetime_called == True:
-                        plt.plot(self.x_mem[i], self.best_fit_mem[i],'k--')
+                    if self.fit_lifetime_called_w_irf == True or self.fit_lifetime_called_wo_irf == True:
+                        plt.plot(self.best_fit_mem_x[i], self.best_fit_mem[i],'k--')
 
                 plt.yscale('log')
                 plt.xlabel("Time (ns)", fontsize=20, fontweight='bold')