update: inertia benefit from caching + move from pkg_resources to importlib.metada

.github/workflows/python-package.yml

@@ -33,18 +33,23 @@ jobs:
         with:
           poetry-version: ${{ matrix.poetry-version }}
 
+      - name: Configure Poetry to use matrix Python version
+        run: |
+          poetry env use python${{ matrix.python-version }}
+
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
           if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
           poetry install --extras "cli fast solvents"
 
       - name: Lint with Ruff
-      - run: |
+        run: |
           poetry run ruff check --output-format=github .
 
       - name: Lint with Black
-        uses: psf/black@stable
+        run: |
+          poetry run black --check . 
 
       - name: Test with pytest
         run: |

.vscode/extensions.json

@@ -0,0 +1,5 @@
+{
+    "recommendations": [
+        "kevinrose.vsc-python-indent"
+    ]
+}

overreact/__init__.py

@@ -4,7 +4,7 @@
 
 __docformat__ = "restructuredtext"
 
-import pkg_resources as _pkg_resources
+from importlib.metadata import version
 
 from overreact.api import (
     get_enthalpies,
@@ -48,7 +48,7 @@
     "unparse_reactions",
 ]
 
-__version__ = _pkg_resources.get_distribution(__name__).version
+__version__ = version(__name__)
 __license__ = "MIT"  # I'm too lazy to get it from setup.py...
 
 __headline__ = "📈 Create and analyze chemical microkinetic models built from computational chemistry data."

overreact/_misc.py

@@ -6,15 +6,270 @@
 from __future__ import annotations
 
 import contextlib
+from copy import deepcopy
 from functools import lru_cache as cache
+from functools import wraps
 
 import numpy as np
+from numpy import arange, array, hstack, newaxis, prod
 from scipy.stats import cauchy, norm
 
 import overreact as rx
 from overreact import _constants as constants
 
 
+def _central_diff_weights(Np, ndiv=1):
+    """
+    Return weights for an Np-point central derivative.
+
+    Assumes equally-spaced function points.
+
+    If weights are in the vector w, then
+    derivative is w[0] * f(x-ho*dx) + ... + w[-1] * f(x+h0*dx)
+
+    Extracted directly from Scipy 'finite_differences' module.
+    (https://github.com/scipy/scipy/blob/d1073acbc804b721cfe356969d8461cdd25a7839/scipy/stats/_finite_differences.py)
+
+    Parameters
+    ----------
+    Np : int
+        Number of points for the central derivative.
+    ndiv : int, optional
+        Number of divisions. Default is 1.
+
+    Returns
+    -------
+    w : ndarray
+        Weights for an Np-point central derivative. Its size is `Np`.
+
+    Notes
+    -----
+    Can be inaccurate for a large number of points.
+
+    Examples
+    --------
+    We can calculate a derivative value of a function.
+
+    >>> def f(x):
+    ...     return 2 * x**2 + 3
+    >>> x = 3.0 # derivative point
+    >>> h = 0.1 # differential step
+    >>> Np = 3 # point number for central derivative
+    >>> weights = _central_diff_weights(Np) # weights for first derivative
+    >>> vals = [f(x + (i - Np/2) * h) for i in range(Np)]
+    >>> sum(w * v for (w, v) in zip(weights, vals))/h
+    11.79999999999998
+
+    This value is close to the analytical solution:
+    f'(x) = 4x, so f'(3) = 12
+
+    References
+    ----------
+    .. [1] https://en.wikipedia.org/wiki/Finite_difference
+
+    """
+    if Np < ndiv + 1:
+        msg = "Number of points must be at least the derivative order + 1."
+        raise ValueError(msg)
+    if Np % 2 == 0:
+        msg = "The number of points must be odd."
+        raise ValueError(msg)
+    from scipy import linalg
+
+    ho = Np >> 1
+    x = arange(-ho, ho + 1.0)
+    x = x[:, newaxis]
+    X = x**0.0
+    for k in range(1, Np):
+        X = hstack([X, x**k])
+    return prod(arange(1, ndiv + 1), axis=0) * linalg.inv(X)[ndiv]
+
+
+def _derivative(func, x0, dx=1.0, n=1, args=(), order=3):
+    """
+    Find the nth derivative of a function at a point.
+
+    Given a function, use a central difference formula with spacing `dx` to
+    compute the nth derivative at `x0`.
+
+    Extracted directly from Scipy 'finite_differences' module.
+    (https://github.com/scipy/scipy/blob/d1073acbc804b721cfe356969d8461cdd25a7839/scipy/stats/_finite_differences.py)
+
+    Parameters
+    ----------
+    func : function
+        Input function.
+    x0 : float
+        The point at which the nth derivative is found.
+    dx : float, optional
+        Spacing.
+    n : int, optional
+        Order of the derivative. Default is 1.
+    args : tuple, optional
+        Arguments
+    order : int, optional
+        Number of points to use, must be odd.
+
+    Notes
+    -----
+    Decreasing the step size too small can result in round-off error.
+
+    Examples
+    --------
+    >>> def f(x):
+    ...     return x**3 + x**2
+    >>> _derivative(f, 1.0, dx=1e-6)
+    4.9999999999217337
+    """
+    first_deriv_weight_map = {
+        3: array([-1, 0, 1]) / 2.0,
+        5: array([1, -8, 0, 8, -1]) / 12.0,
+        7: array([-1, 9, -45, 0, 45, -9, 1]) / 60.0,
+        9: array([3, -32, 168, -672, 0, 672, -168, 32, -3]) / 840.0,
+    }
+
+    second_deriv_weight_map = {
+        3: array([1, -2.0, 1]),
+        5: array([-1, 16, -30, 16, -1]) / 12.0,
+        7: array([2, -27, 270, -490, 270, -27, 2]) / 180.0,
+        9: array([-9, 128, -1008, 8064, -14350, 8064, -1008, 128, -9]) / 5040.0,
+    }
+
+    if order < n + 1:
+        msg = "'order' (the number of points used to compute the derivative), must be at least the derivative order 'n' + 1."
+        raise ValueError(msg)
+    if order % 2 == 0:
+        msg = (
+            "'order' (the number of points used to compute the derivative) must be odd."
+        )
+        raise ValueError(msg)
+
+    # pre-computed for n=1 and 2 and low-order for speed.
+    if n == 1:
+        if order == 3:
+            weights = first_deriv_weight_map.get(3)
+        elif n == 1 and order == 5:
+            weights = first_deriv_weight_map.get(5)
+        elif n == 1 and order == 7:
+            weights = first_deriv_weight_map.get(7)
+        elif n == 1 and order == 9:
+            weights = first_deriv_weight_map.get(9)
+        else:
+            weights = _central_diff_weights(order, 1)
+    # TODO(mrauen): I couldn't find a case in overreact where we use the second (or higher) derivatives. Therefore, I think we can delete this piece of code...Or maybe just leave it here for the future implementations (who knows)
+    elif n == 2:
+        if order == 3:
+            weights = second_deriv_weight_map.get(3)
+        elif n == 2 and order == 5:
+            weights = second_deriv_weight_map.get(5)
+        elif n == 2 and order == 7:
+            weights = second_deriv_weight_map.get(7)
+        elif n == 2 and order == 9:
+            weights = second_deriv_weight_map.get(9)
+        else:
+            weights = _central_diff_weights(order, 2)
+    else:
+        weights = _central_diff_weights(order, n)
+
+    val = 0.0
+    ho = order >> 1
+    for k in range(order):
+        val += weights[k] * func(x0 + (k - ho) * dx, *args)
+    return val / prod((dx,) * n, axis=0)
+
+
+def make_hashable(obj):
+    """
+    Given an array, list or set make it immutable by transforming it into a tuple.
+
+    Parameters
+    ----------
+    obj : array
+
+    Returns
+    -------
+    tuple
+
+    Notes
+    -----
+    List comprehension it's key here for list and set, otherwise it will return a tuple with only the first item.
+    """
+    if isinstance(obj, np.ndarray):
+        return (tuple(obj.shape), tuple(obj.ravel()))
+    elif isinstance(obj, (list, set)):
+        return tuple(make_hashable(item) for item in obj)
+    else:
+        return obj
+
+
+def copy_unhashable(maxsize=128, typed=False):
+    """
+    Cache resultant tuples while handling the received unhashable types (array, list, dictionaries).
+
+    Convert unhashable arguments into hashable before passing it to 'lru_cache'. Then, reconstruct the (now) hashable tuple back to return it for the function caller. A copy of the received argument is made in order to prevent errors and side-effects to the original array/list/etc.
+
+    Parameters
+    ----------
+    maxsize : int
+        Cache size limit. Default from functools.lru_cache()
+    typed : bool
+        If set to True, arguments of different types will be cache separately. Default from functools.lru_cache()
+    func : function
+        The function to be wrapped and cached
+
+    Returns
+    -------
+    function
+        A wrapper version of the original function that is cacheable now
+    """
+
+    def decorator(func):
+        @cache(maxsize=maxsize, typed=typed)
+        @wraps(func)
+        def cached_func(*hashable_args, **hashable_kwargs):
+            args = []
+            kwargs = {}
+
+            def convert_back(arg):
+                if isinstance(arg, tuple) and len(arg) == 2:
+                    shape, flat_data = arg
+                    if (
+                        isinstance(shape, tuple)
+                        and all(isinstance(dim, (int, np.integer)) for dim in shape)
+                        and isinstance(flat_data, tuple)
+                    ):
+                        if len(flat_data) == 0 or any(dim <= 0 for dim in shape):
+                            return np.array([])
+                        try:
+                            return np.array(flat_data).reshape(shape)
+                        except ValueError as e:
+                            msg = f"Reshape error: {e} - shape: {shape}, data: {flat_data}"
+                            raise ValueError(msg)
+                return arg
+
+            args = [convert_back(arg) for arg in hashable_args]
+            for k, v in hashable_kwargs.items():
+                kwargs[k] = convert_back(v)
+            args = tuple(args)
+            return func(*args, **kwargs)
+
+        def wrapper(*args, **kwargs):
+            wrapper_hashable_args = []
+            wrapper_hashable_kwargs = {}
+
+            wrapper_hashable_args = [make_hashable(arg) for arg in args]
+            for k, v in kwargs.items():
+                wrapper_hashable_kwargs[k] = make_hashable(v)
+            wrapper_hashable_args = tuple(wrapper_hashable_args)
+            return deepcopy(
+                cached_func(*wrapper_hashable_args, **wrapper_hashable_kwargs),
+            )
+
+        return wrapper
+
+    return decorator
+
+
 def _find_package(package):
     """Check if a package exists without importing it.
 
@@ -739,7 +994,7 @@
     return primes
 
 
-@cache(maxsize=1000000)
+@cache
 def _vdc(n, b=2):
     """Help haltonspace."""
     res, denom = 0, 1

overreact/api.py

@@ -8,12 +8,12 @@
 from __future__ import annotations
 
 __all__ = [
-    "get_k",
-    "get_kappa",
-    "get_freeenergies",
-    "get_entropies",
     "get_enthalpies",
+    "get_entropies",
+    "get_freeenergies",
     "get_internal_energies",
+    "get_k",
+    "get_kappa",
 ]
 
 
@@ -22,11 +22,11 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from scipy.misc import derivative
 
 import overreact as rx
 from overreact import _constants as constants
 from overreact import coords, rates, tunnel
+from overreact._misc import _derivative as derivative
 
 if TYPE_CHECKING:
     from overreact.core import Scheme
@@ -74,7 +74,6 @@ def get_internal_energies(
     for name in compounds:
         logger.info(f"calculate internal energy: {name}")
 
-        # TODO(schneiderfelipe): inertia might benefit from caching
         moments, _, _ = coords.inertia(
             compounds[name].atommasses,
             compounds[name].atomcoords,
@@ -140,7 +139,6 @@ def get_enthalpies(
     for name in compounds:
         logger.info(f"calculate enthalpy: {name}")
 
-        # TODO(schneiderfelipe): inertia might benefit from caching
         moments, _, _ = coords.inertia(
             compounds[name].atommasses,
             compounds[name].atomcoords,
@@ -233,7 +231,6 @@ def get_entropies(
             )
         symmetry_number = coords.symmetry_number(point_group)
 
-        # TODO(schneiderfelipe): inertia might benefit from caching
         moments, _, _ = coords.inertia(
             compounds[name].atommasses,
             compounds[name].atomcoords,