include lmfit

src/pychemelt/thermal_oligomer.py

@@ -692,11 +692,10 @@ def fit_thermal_unfolding_three_state_global(
                 raise ValueError('CpTh must be large enough for fitting. If you do not wish to fit the Cp values, omit this parameter.')
 
             # Parameters for T1, T2, will be added later via gridsearch
-            p0 = [0, 250, 0, 250, self.Cp0]
+            p0 = [0, 200, 0, 200, self.Cp0]
 
         else:
-            p0 = [0, 250, 0, 250, 0]
-
+            p0 = [0, 200, 0, 200, 0]
 
         params_names = [
             'Tm1 (°C)',
@@ -787,7 +786,7 @@ def fit_thermal_unfolding_three_state_global(
             tm1_lower = 15
             tm1_upper = self.user_max_temp + 20
 
-            tm2_lower = 15
+            tm2_lower = 30
             tm2_upper = self.user_max_temp + 20
 
         low_bounds[0] = tm1_lower
@@ -810,10 +809,19 @@ def fit_thermal_unfolding_three_state_global(
             p0[3] = adjust_value_to_interval(p0[3], dh2_lower, dh2_upper, 1)
 
         else:
-            dh1_lower = 10
+            # Set dh1_lower 30 for monomer, 50 for dimer, 70 for trimer and 90 for tetramer, and dh1_upper to 500 for all
+            lower_value = 30
+            if "Dimer" in self.model:
+                lower_value = 50
+            elif "Trimer" in self.model:
+                lower_value = 70
+            elif "Tetramer" in self.model:
+                lower_value = 90
+
+            dh1_lower = lower_value
             dh1_upper = 500
 
-            dh2_lower = 10
+            dh2_lower = lower_value
             dh2_upper = 500
 
         low_bounds[1] = dh1_lower
@@ -837,9 +845,9 @@ def fit_thermal_unfolding_three_state_global(
 
         else:
 
-            cp_lower, cp_upper = 0.1, CpTh - 0.3
+            cp_lower, cp_upper = 0.1, CpTh - 0.4
 
-            low_bounds[4] = cp_lower
+            low_bounds[4]  = cp_lower
             high_bounds[4] = cp_upper
 
             # Verify that the Cp initial guess is within the user-defined limits
@@ -851,13 +859,11 @@ def fit_thermal_unfolding_three_state_global(
 
         signal_fx = map_three_state_model_to_signal_fx(self.model)
 
-        has_nmeric_intermediate = not "monomeric_intermediate" in self.model.lower()
-
         if t1_init != 0:
-            p0[0], low_bounds[0], high_bounds[0] = t1_init, np.max([t1_init - 15, 0]), t1_init + 20
+            p0[0], low_bounds[0], high_bounds[0] = t1_init, np.max([t1_init - 15, 20]), t1_init + 20
 
         if t2_init != 0:
-            p0[2], low_bounds[2], high_bounds[2] = t2_init, np.max([t2_init - 15, 0]), t2_init + 20
+            p0[2], low_bounds[2], high_bounds[2] = t2_init, np.max([t2_init - 15, 20]), t2_init + 20
 
         kwargs = {
             'oligomer_concentrations': self.oligomer_concentrations_expanded,
@@ -872,23 +878,17 @@ def fit_thermal_unfolding_three_state_global(
 
         fit_fx = fit_oligomer_unfolding_three_states_single_slopes
 
-        step = 6
-        num_rows = len(self.oligomer_concentrations)
-
-        if num_rows > 3:
-            step += 2
-        if num_rows > 4:
-            step += 2
+        step = 7
 
         if t1_init == 0 or t2_init == 0:
 
             if self.limited_tm:
                 test_T1s = np.arange(np.max([tm1_lower, 20]), tm1_upper, step)
-                test_T2s = np.arange(np.max([tm2_lower, 20]) + step, tm2_upper, step)
+                test_T2s = np.arange(np.max([tm2_lower, 35]) + step, tm2_upper, step)
 
             else:
-                test_T1s = np.arange(np.max([self.global_min_temp + 10, 20]), self.global_max_temp - 20*has_nmeric_intermediate, step)
-                test_T2s = np.arange(np.max([self.global_min_temp + 10, 20]) + step, self.global_max_temp + 5, step)
+                test_T1s = np.arange(np.max([self.global_min_temp + 10, 20]), self.global_max_temp, step)
+                test_T2s = np.arange(np.max([self.global_min_temp + 20, 20]) + step, self.global_max_temp + step, step)
 
             if t1_init != 0:
                 test_T1s = np.array([t1_init])
@@ -900,6 +900,9 @@ def fit_thermal_unfolding_three_state_global(
 
             df = pd.DataFrame(combinations, columns=['t1', 't2'])
 
+            # Remove all combinations where t1 is 30 degrees higher than t2
+            df = df[~(df['t1'] - df['t2'] >= 30)]
+
             df_tm = pd.DataFrame(np.zeros_like(combinations), columns=['t1', 't2'], dtype=float)
             df_dh = pd.DataFrame(np.zeros_like(combinations), columns=['dh1', 'dh2'], dtype=float)
 
@@ -910,10 +913,13 @@ def fit_thermal_unfolding_three_state_global(
             kwargs['list_of_signals'] = self.signal_lst_expanded_subset
 
             for index, row in df.iterrows():
+
                 kwargs['t1'] = row['t1']
                 kwargs['t2'] = row['t2']
 
-                fit_params, cov, pred, result, minimizer = fit_oligomer_unfolding_three_states_single_slopes(**kwargs)
+                fit_params, cov, pred, result, minimizer = fit_oligomer_unfolding_three_states_single_slopes(
+                    **kwargs,
+                    max_nfev=6000) # We need to limit max_nfev for a faster initial grid search
 
                 #using the fitted parameters as a base for fitting
                 df_tm.iloc[index, 0] = fit_params[0]
@@ -933,16 +939,18 @@ def fit_thermal_unfolding_three_state_global(
             dh1_init, dh2_init = df_dh['dh1'][idx], df_dh['dh2'][idx]
             p0[1], p0[3] = dh1_init, dh2_init
 
-            low_bounds[0], low_bounds[2] = t1_init - 15, t2_init - 15
+            low_bounds[0], low_bounds[2] = t1_init - 14, t2_init - 14
             high_bounds[0], high_bounds[2] = t1_init + 18, t2_init + 18
 
+            low_bounds[1], low_bounds[3] = np.max([dh1_init - 100, 30]), np.max([dh2_init - 100, 30])
+            high_bounds[1], high_bounds[3] = dh1_init + 150, dh2_init + 150
+
             kwargs['initial_parameters'] = p0
             kwargs['low_bounds'] = low_bounds
             kwargs['high_bounds'] = high_bounds
             kwargs['t1'] = None
             kwargs['t2'] = None
 
-
         # Now use the whole dataset
         kwargs['list_of_temperatures'] = self.temp_lst_expanded
         kwargs['list_of_signals'] = self.signal_lst_expanded
@@ -966,7 +974,7 @@ def fit_thermal_unfolding_three_state_global(
             result=result,
             minimizer=minimizer,
             fit_m_value=False,
-            three_state_model=True,
+            three_state_model=True
         )
 
         rel_errors = relative_errors(global_fit_params, cov)

src/pychemelt/utils/fitting.py

@@ -989,6 +989,8 @@ def fit_oligomer_unfolding_three_states_single_slopes(
     dh2=None,
     CpTh_value=None,
     method="least_squares",
+    max_nfev=None
+
 ):
     """
     Vectorized and optimized version of global thermal unfolding fitting of oligomers.
@@ -1041,17 +1043,17 @@ def fit_oligomer_unfolding_three_states_single_slopes(
         high_bounds[0] = temperature_to_kelvin(high_bounds[0])
     else:
         initial_parameters[0] = temperature_to_kelvin(t1)
-        low_bounds[0] = initial_parameters[0] - 12
-        high_bounds[0] = initial_parameters[0] + 18
+        low_bounds[0] = initial_parameters[0] - 8
+        high_bounds[0] = initial_parameters[0] + 12
 
     if t2 is None:
         initial_parameters[2] = temperature_to_kelvin(initial_parameters[2])
         low_bounds[2] = temperature_to_kelvin(low_bounds[2])
         high_bounds[2] = temperature_to_kelvin(high_bounds[2])
     else:
         initial_parameters[2] = temperature_to_kelvin(t2)
-        low_bounds[2] = initial_parameters[2] - 12
-        high_bounds[2] = initial_parameters[2] + 18
+        low_bounds[2] = initial_parameters[2] - 8
+        high_bounds[2] = initial_parameters[2] + 12
 
     if dh1 is not None:
         initial_parameters[1] = dh1
@@ -1228,7 +1230,7 @@ def model(pars):
     def residuals(pars):
         return model(pars) - y_all
 
-    minimizer = lmfit.Minimizer(residuals, params_lmfit, calc_covar=True)
+    minimizer = lmfit.Minimizer(residuals, params_lmfit, calc_covar=True,max_nfev=max_nfev)
     result = minimizer.minimize(method=method)
 
     global_fit_params = np.array([result.params[name].value for name in param_names], dtype=float)

tests/test_fitting_therm_oligo.py

@@ -21,8 +21,8 @@
 RNG_SEED = 2
 TEMP_START = 30.0
 TEMP_STOP = 90.0
-N_TEMPS = 80
-CONCS = np.arange(10, 100, 10)*1e-6
+N_TEMPS = 70
+CONCS = np.array([1,3,9,81])*1e-6
 
 # Model / ground-truth parameters
 DHm_VAL = 250
@@ -72,7 +72,7 @@ def test_fit_monomer_unfolding_single_slopes_constant():
         temp_list.append(temp_range)
 
     p0 = [Tm_VAL, DHm_VAL, CP0_VAL] + [INTERCEPT_N] *len(concs) + [INTERCEPT_U] * len(concs)
-    low_bounds = [TEMP_START, TEMP_START, 1]   + [1e-5]*(2*len(concs))
+    low_bounds = [TEMP_START, DHm_VAL-50, 1]   + [1e-5]*(2*len(concs))
     high_bounds = [TEMP_STOP, DHm_VAL+100, 5] + [1e3]*(2*len(concs))
 
     kwargs = {
@@ -84,7 +84,7 @@ def test_fit_monomer_unfolding_single_slopes_constant():
         'baseline_unfolded_fx':constant_baseline, 
     }
 
-    global_fit_params, cov, predicted_lst, _, _ = fit_oligomer_unfolding_single_slopes(
+    global_fit_params, cov, predicted_lst, result, minimizer = fit_oligomer_unfolding_single_slopes(
         initial_parameters=p0,
         low_bounds=low_bounds,
         high_bounds=high_bounds,
@@ -93,7 +93,12 @@ def test_fit_monomer_unfolding_single_slopes_constant():
 
     expected = [Tm_VAL, DHm_VAL, CP0_VAL]
 
-    np.testing.assert_allclose(global_fit_params[:3], expected, rtol=0.1, atol=0)
+    error = np.abs(np.sqrt(np.diag(cov)))
+    params_minus_error = (global_fit_params - error)*0.98
+    params_plus_error = (global_fit_params + error)*1.02
+
+    for i in range(3):
+        assert params_minus_error[i] <= expected[i] <= params_plus_error[i], f"Parameter {i}: {expected[i]} not in [{params_minus_error[i]}, {params_plus_error[i]}]"
 
     # Fit with fixed Tm
     p0_tm = p0.copy()
@@ -114,7 +119,12 @@ def test_fit_monomer_unfolding_single_slopes_constant():
 
     expected = [DHm_VAL, CP0_VAL]
 
-    np.testing.assert_allclose(global_fit_params[:2], expected, rtol=0.1, atol=0)
+    error = np.abs(np.sqrt(np.diag(cov)))
+    params_minus_error = (global_fit_params - error)*0.98
+    params_plus_error = (global_fit_params + error)*1.02
+
+    for i in range(2):
+        assert params_minus_error[i] <= expected[i] <= params_plus_error[i], f"Parameter {i}: {expected[i]} not in [{params_minus_error[i]}, {params_plus_error[i]}]"
 
     # End of - Fit with fixed Tm
 
@@ -138,7 +148,13 @@ def test_fit_monomer_unfolding_single_slopes_constant():
 
     expected = [Tm_VAL, CP0_VAL]
 
-    np.testing.assert_allclose(global_fit_params[:2], expected, rtol=0.1, atol=0)
+    error = np.abs(np.sqrt(np.diag(cov)))
+    params_minus_error = (global_fit_params - error)*0.98
+    params_plus_error = (global_fit_params + error)*1.02
+
+    for i in range(2):
+        assert params_minus_error[i] <= expected[i] <= params_plus_error[i], f"Parameter {i}: {expected[i]} not in [{params_minus_error[i]}, {params_plus_error[i]}]"
+
 
     # End of - Fit with fixed dH
 
@@ -787,8 +803,13 @@ def test_fit_monomer_unfolding_many_signals_constant():
     )
 
     expected = [Tm_VAL, DHm_VAL, CP0_VAL]
+    error = np.abs(np.sqrt(np.diag(cov)))
+    params_minus_error = (global_fit_params - error)*0.99
+    params_plus_error = (global_fit_params + error)*1.01
+
+    for i in range(3):
+        assert params_minus_error[i] <= expected[i] <= params_plus_error[i], f"Parameter {i}: {expected[i]} not in [{params_minus_error[i]}, {params_plus_error[i]}]"
 
-    np.testing.assert_allclose(global_fit_params[:3], expected, rtol=0.1, atol=0)
 
     # Fit scale factor
 
@@ -813,7 +834,13 @@ def test_fit_monomer_unfolding_many_signals_constant():
 
     expected = [Tm_VAL, DHm_VAL, CP0_VAL]
 
-    np.testing.assert_allclose(global_fit_params[:3], expected, rtol=0.1, atol=0)
+    expected = [Tm_VAL, DHm_VAL, CP0_VAL]
+    error = np.abs(np.sqrt(np.diag(cov)))
+    params_minus_error = (global_fit_params - error*2) # Two sigma 
+    params_plus_error = (global_fit_params + error*2) 
+
+    for i in range(3):
+        assert params_minus_error[i] <= expected[i] <= params_plus_error[i], f"Parameter {i}: {expected[i]} not in [{params_minus_error[i]}, {params_plus_error[i]}]"
 
     # Fit cp value set
 
@@ -831,7 +858,14 @@ def test_fit_monomer_unfolding_many_signals_constant():
 
     expected = [Tm_VAL, DHm_VAL]
 
-    np.testing.assert_allclose(global_fit_params[:2], expected, rtol=0.1, atol=0)
+    error = np.abs(np.sqrt(np.diag(cov)))
+    params_minus_error = (global_fit_params - error)*0.99
+    params_plus_error = (global_fit_params + error)*1.01
+
+    for i in range(2):
+        assert params_minus_error[i] <= expected[i] <= params_plus_error[i], f"Parameter {i}: {expected[i]} not in [{params_minus_error[i]}, {params_plus_error[i]}]"
+
+
 
 def test_fit_dimer_unfolding_many_signals_constant():
     signal_fx = map_two_state_model_to_signal_fx("Dimer")

tests/test_fitting_therm_oligo_2.py

@@ -23,7 +23,7 @@
 TEMP_START = 20.0
 TEMP_STOP = 90.0
 N_TEMPS = 80
-CONCS = np.arange(10, 100, 10)*1e-6
+CONCS = np.array([10,20,30,40,50])*1e-6
 
 # Model / ground-truth parameters
 DHm_VAL = 250
@@ -453,6 +453,7 @@ def test_fit_monomer_unfolding_many_signals_exponential():
     np.testing.assert_allclose(global_fit_params[:2], [Tm_VAL,DHm_VAL], rtol=0.1, atol=1e-2)
 
 def test_fit_dimer_unfolding_many_signals_exponential():
+
     signal_fx = map_two_state_model_to_signal_fx("Dimer")
 
     signal_list = []
@@ -472,8 +473,8 @@ def test_fit_dimer_unfolding_many_signals_exponential():
     p0 += [C_N_VAL, C_U_VAL]
     p0 += [ALPHA_N_VAL,ALPHA_U_VAL]
 
-    low_bounds  = [0 for _ in p0]
-    high_bounds = [1e3 for _ in p0]
+    low_bounds  = [Tm_VAL-10,DHm_VAL-50] + [1e-3 for _ in p0[2:]]
+    high_bounds = [Tm_VAL+10,DHm_VAL+50] + [500 for _ in p0[2:]]
 
     kwargs = {
         'list_of_temperatures' : temp_list,
@@ -517,7 +518,7 @@ def test_fit_trimer_unfolding_many_signals_exponential():
     p0 += [ALPHA_N_VAL,ALPHA_U_VAL]
 
     low_bounds = [-0.1 for _ in p0]
-    high_bounds = [1e3 for _ in p0]
+    high_bounds = [500 for _ in p0]
 
 
     kwargs = {
@@ -561,7 +562,7 @@ def test_fit_tetramer_unfolding_many_signals_exponential():
     p0 += [ALPHA_N_VAL,ALPHA_U_VAL]
 
     low_bounds = [-0.1 for _ in p0]
-    high_bounds = [1e3 for _ in p0]
+    high_bounds = [500 for _ in p0]
 
 
     kwargs = {

tests/test_fitting_therm_oligo_Cp_value.py

@@ -27,7 +27,7 @@
 TEMP_START = 30.0
 TEMP_STOP = 90.0
 N_TEMPS = 80
-CONCS = np.arange(10, 100, 10)*1e-6
+CONCS = np.array([2,6,24,72])*1e-6
 
 # Two state
 # Model / ground-truth parameters