Merge pull request #21 from PPeitsch/19-fix-kissinger-module

19 fix kissinger module

__init__.py

@@ -6,7 +6,7 @@
 approaches, and result visualization for thermal analysis data.
 """
 
-__version__ = "0.1.0"
+__version__ = "0.2.1"
 
 from . import data_preprocessing
 from . import data_import
@@ -19,10 +19,9 @@
 from . import parallel_processing
 from . import utility_functions
 
-# Convenience imports for commonly used functions
 from .data_preprocessing import smooth_data
 from .data_import import tga_importer, dsc_importer
-from .result_visualization import plot_activation_energy
+from .model_fitting_methods import kissinger_method
 
 __all__ = [
     "data_preprocessing",
@@ -38,6 +37,5 @@
     "smooth_data",
     "tga_importer",
     "dsc_importer",
-    "synthetic_data",
-    "plot_activation_energy",
+    "kissinger_method",
 ]

model_fitting_methods/kissinger.py

@@ -1,17 +1,33 @@
 """Kissinger method for non-isothermal kinetics analysis."""
 
 import numpy as np
-from scipy import stats
 from scipy.optimize import fsolve
 from typing import Tuple
 import logging
+import statsmodels.api as sm
 
 logging.basicConfig(level=logging.INFO)
 logger = logging.getLogger(__name__)
 
 R = 8.314  # Gas constant in J/(mol·K)
 
 
+def kissinger_nonlinear_eq(t: float, e_a: float, a: float, b: float) -> float:
+    """
+    Kissinger nonlinear equation for peak temperature calculation.
+
+    Args:
+        t (float): Temperature in K.
+        e_a (float): Activation energy in J/mol.
+        a (float): Pre-exponential factor in min^-1.
+        b (float): Heating rate in K/min.
+
+    Returns:
+        float: Value of the Kissinger equation at the given temperature.
+    """
+    return (e_a * b) / (R * t ** 2) - a * np.exp(-e_a / (R * t))
+
+
 def calculate_t_p(e_a: float, a: float, beta: np.ndarray) -> np.ndarray:
     """
     Calculate peak temperatures using the Kissinger equation.
@@ -24,25 +40,22 @@ def calculate_t_p(e_a: float, a: float, beta: np.ndarray) -> np.ndarray:
     Returns:
         np.ndarray: Calculated peak temperatures in K.
     """
-
-    def kissinger_equation(t, b):
-        return (e_a * b) / (R * t ** 2) - a * np.exp(-e_a / (R * t))
-
     t_p = np.zeros_like(beta)
     for i, b in enumerate(beta):
-        t_p[i] = fsolve(kissinger_equation, x0=e_a / (R * 20), args=(b,))[0]
-
+        try:
+            t_p[i] = fsolve(kissinger_nonlinear_eq, x0=e_a / (R * 20), args=(e_a, a, b))[0]
+        except (RuntimeError, ValueError) as e:
+            logger.warning(f"Failed to converge for heating rate {b}: {e}. Using initial guess.")
+            t_p[i] = e_a / (R * 20)
     return t_p
 
 
-def kissinger_equation(t_p: np.ndarray, e_a: float, a: float, beta: np.ndarray) -> np.ndarray:
+def kissinger_equation(t_p: np.ndarray, beta: np.ndarray) -> np.ndarray:
     """
     Kissinger equation for non-isothermal kinetics.
 
     Args:
         t_p (np.ndarray): Peak temperatures in K.
-        e_a (float): Activation energy in J/mol.
-        a (float): Pre-exponential factor in min^-1.
         beta (np.ndarray): Heating rates in K/min.
 
     Returns:
@@ -73,19 +86,31 @@ def kissinger_method(t_p: np.ndarray, beta: np.ndarray) -> Tuple[float, float, f
         raise ValueError("Peak temperatures and heating rates must be positive")
 
     x = 1 / t_p
-    y = np.log(beta / t_p ** 2)
+    y = kissinger_equation(t_p, beta)
 
-    slope, intercept, r_value, _, stderr = stats.linregress(x, y)
+    X = sm.add_constant(x)
+    model = sm.OLS(y, X).fit()
+
+    slope = model.params[1]
+    intercept = model.params[0]
+    se_slope = model.bse[1]
+    se_intercept = model.bse[0]
 
     e_a = -R * slope
+    se_e_a = R * se_slope
     a = np.exp(intercept + np.log(e_a / R))
+    se_ln_a = np.sqrt(se_intercept ** 2 + (se_e_a / e_a) ** 2)
 
-    se_e_a = R * stderr
-    se_ln_a = np.sqrt((stderr / slope) ** 2 + (se_e_a / e_a) ** 2)
-
-    r_squared = r_value ** 2
+    r_squared = model.rsquared
 
     logger.info(f"Kissinger analysis completed. E_a = {e_a / 1000:.2f} ± {se_e_a / 1000:.2f} kJ/mol, "
                 f"A = {a:.2e} min^-1, R^2 = {r_squared:.4f}")
 
+    adjusted_r_squared = model.rsquared_adj
+    f_statistic = model.fvalue
+    p_values = model.pvalues
+
+    logger.info(f"Adjusted R^2 = {adjusted_r_squared:.4f}, F-statistic = {f_statistic:.2f}")
+    logger.info(f"P-values: {p_values}")
+
     return e_a, a, se_e_a, se_ln_a, r_squared

requirements.txt

@@ -0,0 +1,10 @@
+numpy~=1.24.3
+matplotlib~=3.7.5
+setuptools~=68.2.0
+pandas~=2.0.3
+scipy~=1.10.1
+statsmodels~=0.14.1
+scikit-learn>=0.24.0
+seaborn>=0.11.0
+pytest>=6.2.0
+pylint>=2.7.0

setup.py

@@ -2,20 +2,22 @@
 
 setup(
     name="pkynetics",
-    version="0.1.0",
+    version="0.2.1",
     packages=find_packages(),
     install_requires=[
-        "numpy>=1.20.0",
-        "pandas>=1.2.0",
-        "scipy>=1.6.0",
+        "numpy~=1.24.3",
+        "matplotlib~=3.7.5",
+        "pandas~=2.0.3",
+        "scipy~=1.10.1",
+        "statsmodels~=0.14.1",
         "scikit-learn>=0.24.0",
-        "matplotlib>=3.3.0",
         "seaborn>=0.11.0",
-        "tensorflow>=2.5.0",
-        "torch>=1.8.0",
     ],
     extras_require={
-        "dev": ["pytest>=6.2.0", "pylint>=2.7.0"],
+        "dev": [
+            "pytest>=6.2.0",
+            "pylint>=2.7.0",
+        ],
     },
     entry_points={
         "console_scripts": [