misc: Curve fitting to biexponential function

Removing B in lieu of 1-A gets rid of "bad" results (since B can be 0, and therefore tau2 could have huge/bad values, leading to a huge / nonsensical value of the lifetime) Co-authored-by: Moritz Sallermann <[email protected]>
amritagos · Aug 23, 2024 · 027aaf1 · 027aaf1
1 parent a6bd557
commit 027aaf1
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 8 deletions.
diff --git a/examples/hydrogen_bond_tcf/lineplot.py b/examples/hydrogen_bond_tcf/lineplot.py
@@ -93,11 +93,10 @@ def insert_inset(image_path, axis, rel_height,
 # ---------------------------------------------
 # Curve fitting 
 # initial guess A, tau1, B, tau2
-params, fit_t, fit_ac, lifetime = fit_biexponential(tau_val, tcf_val, [0.5, 1, 1, 2])
+params, fit_t, fit_ac, lifetime = fit_biexponential(tau_val, tcf_val, [0.5, 1, 2])
 print("Lifetime is ", lifetime, "ps \n")
-A, tau1, B, tau2 = params
+A, tau1, tau2 = params
 # If A+B is around 1.01, the fit is bad even if it looks okay
-print(f"Sum A + B = {A+B} and should be 1.0\n")
 print(f"The time constants are {tau1} ps and {tau2} ps\n")
 # ---------------------------------------------
 ax1 = fig.add_subplot(gs[0,0])
@@ -125,6 +124,7 @@ def insert_inset(image_path, axis, rel_height,
 # xtick_vec = np.append(xtick_vec, 2.25)
 # ax1.set_xticks(xtick_vec)
 # ax1.set_ylim([0.0,25])
+# ax1.set_yscale("log")
 ax1.set_xlim([0.0,10])
 
 # ---------------------------------------------------------------------

diff --git a/soluanalysis/misc.py b/soluanalysis/misc.py
@@ -8,7 +8,7 @@
 
 
 def biexponential_model(
-    t: float, A: float, tau1: float, B: float, tau2: float
+    t: float, A: float, tau1: float, tau2: float
 ) -> float:
     """Fit data to a biexponential function (sum of two exponential decays). A and B should sum to 1.0
     C(t) = A*exp(-t/tau1) + B*exp(-t/tau2)
@@ -23,13 +23,13 @@ def biexponential_model(
     Returns:
         float: _description_
     """
-    return A * np.exp(-t / tau1) + B * np.exp(-t / tau2)
+    return A * np.exp(-t / tau1) + (1-A) * np.exp(-t / tau2)
 
 
 def fit_biexponential(
     tau_timeseries: Union[List[float], npt.NDArray],
     tcf_timeseries: Union[List[float], npt.NDArray],
-    initial_guess: Union[List[float], npt.NDArray] = [0.5, 1, 1, 2],
+    initial_guess: Union[List[float], npt.NDArray] = [0.5, 1, 2],
 ) -> Tuple[Tuple[float, float, float], npt.NDArray, npt.NDArray, float]:
     """Fit a biexponential function of the form (A*exp(-t/tau1) + B*exp(-t/tau2)) to the tau values and time autocorrelation function values (using the continuous bond definition).
     Using the relation from Gowers et al. (2015).
@@ -59,9 +59,9 @@ def fit_biexponential(
     fit_ac = biexponential_model(fit_t, *params)
 
     # Extract the parameters
-    A, tau1, B, tau2 = params
+    A, tau1, tau2 = params
 
     # Compute the lifetime
-    lifetime = A * tau1 + B * tau2
+    lifetime = A * tau1 + (1-A) * tau2
 
     return params, fit_t, fit_ac, lifetime