Enhance fixed-point methods and optimization methods

.github/workflows/test.yml

@@ -10,12 +10,14 @@ jobs:
       fail-fast: false
       matrix:
         resolution: ["lowest-direct", "highest"]
-        python-version: ["3.10", "3.11", "3.12", "3.13"]
+        python-version: ["3.11", "3.12", "3.13", "3.14"]
         exclude:
           - python-version: "3.12"
             resolution: "lowest-direct"
           - python-version: "3.13"
             resolution: "lowest-direct" #numba 0.60 requires python < 3.13
+          - python-version: "3.14"
+            resolution: "lowest-direct" #numba 0.60 requires python < 3.13
     steps:
       - name: Checkout code
         uses: actions/checkout@v4

.vscode/extensions.json

@@ -1,6 +1,7 @@
 {
     "recommendations": [
         "ms-python.python",
-        "charliermarsh.ruff"
+        "charliermarsh.ruff",
+        "astral-sh.ty"
     ]
 }

.vscode/settings.json

@@ -7,10 +7,14 @@
     ],
     "python.testing.unittestEnabled": false,
     "python.testing.pytestEnabled": true,
+    "python.languageServer": "Pylance",
+    "python.analysis.diagnosticMode": "openFilesOnly",
+    "ty.disableLanguageServices": true,
     "[python]": {
         "editor.rulers": [
             90
         ],
+        "editor.inlayHints.enabled": "off",
         "editor.formatOnSave": true,
         "editor.defaultFormatter": "charliermarsh.ruff",
         "editor.codeActionsOnSave": {

README.md

@@ -5,7 +5,7 @@
 ![PyPI - Downloads](https://img.shields.io/pypi/dm/polykin?label=PyPI%20downloads)
 [![Ask DeepWiki](https://deepwiki.com/badge.svg)](https://deepwiki.com/HugoMVale/polykin)
 
-PolyKin is an open-source polymerization kinetics library for Python. It is still at an early
+PolyKin is an open-source Python library for polymerization kinetics and related chemical engineering calculations. It is still at an early
 development stage, but the following modules can already be used:
 
 - Activity coefficient models
@@ -53,8 +53,8 @@ development stage, but the following modules can already be used:
   - [x] Termination composite model
 - Math
   - [x] Joint confidence regions
-  - [x] Root and fixed-point solvers
-- [ ] Models
+  - [x] Fixed-point solvers (scalar and vector)
+  - [x] Root solvers (scalar end vector)
 - Physical property correlations
   - [x] Antoine
   - [x] DIPPR

docs/installation.md

@@ -1,6 +1,6 @@
 # Installation
 
-`PolyKin` requires Python >= 3.10, because it makes use of recent type hint syntax. It further relies on a number of mature and well-maintained mathematical/scientific libraries:
+`PolyKin` requires Python >= 3.11. It further relies on a number of mature and well-maintained mathematical/scientific libraries:
 `matplotlib`, `mpmath`, `numba`, `numpy`, `pydantic`, `scipy`, etc.
 
 ## Install from PyPI

docs/reference/math/derivatives/_SUMMARY.md

@@ -1,6 +1,7 @@
 * [](index.md)
 * [derivative_centered](derivative_centered.md)
 * [derivative_complex](derivative_complex.md)
+* [gradient_forward](gradient_forward.md)
 * [hessian_forward](hessian_forward.md)
 * [hessian2_centered](hessian2_centered.md)
 * [jacobian_forward](jacobian_forward.md)

docs/reference/math/derivatives/gradient_forward.md

@@ -0,0 +1,6 @@
+# polykin.math.derivatives
+
+::: polykin.math.derivatives
+    options:
+        members:
+            - gradient_forward

docs/reference/math/fixpoint/_SUMMARY.md

@@ -1,3 +1,6 @@
 * [](index.md)
 * [fixpoint_anderson](fixpoint_anderson.md)
+* [fixpoint_damped](fixpoint_damped.md)
+* [fixpoint_dem](fixpoint_dem.md)
+* [fixpoint_steffensen](fixpoint_steffensen.md)
 * [fixpoint_wegstein](fixpoint_wegstein.md)

docs/reference/math/fixpoint/fixpoint_damped.md

@@ -0,0 +1,6 @@
+# polykin.math.fixpoint
+
+::: polykin.math.fixpoint
+    options:
+        members:
+            - fixpoint_damped

docs/reference/math/fixpoint/fixpoint_dem.md

@@ -0,0 +1,6 @@
+# polykin.math.fixpoint
+
+::: polykin.math.fixpoint
+    options:
+        members:
+            - fixpoint_dem

docs/reference/math/fixpoint/fixpoint_steffensen.md

@@ -0,0 +1,6 @@
+# polykin.math.fixpoint
+
+::: polykin.math.fixpoint
+    options:
+        members:
+            - fixpoint_steffensen

docs/reference/math/fixpoint/index.md

@@ -1 +1,21 @@
 # Fixed-Point Solvers (polykin.math.fixpoint)
+
+## Scalar Problems
+
+[fixpoint_steffensen] is the only solver in this group specialized for scalar fixed-point problems $g(x)=x$. It accelerates direct substitution through Aitken's delta-squared process, often converging significantly faster than plain fixed-point iteration when the fixed-point sequence is smooth and locally contractive.
+
+## Multivariable Problems
+
+[fixpoint_damped] is the most basic option for systems $\mathbf{g}(\mathbf{x})$. It applies direct substitution with optional damping, making each iteration inexpensive and easy to interpret. This is a reasonable first choice when the mapping is already contractive or only mildly unstable, especially if function evaluations are cheap.
+
+[fixpoint_wegstein] improves on direct substitution by estimating a componentwise secant slope and converting it into a bounded acceleration factor. In practice, it often offers a favorable balance between robustness and speed for weakly to moderately coupled systems. It is usually more effective than simple damping, but because the acceleration is applied componentwise, performance can degrade when the variables are strongly coupled.
+
+[fixpoint_dem] targets problems whose slow convergence is dominated by a single eigenvalue of the fixed-point iteration. By estimating that dominant mode and periodically applying a promotion step, it can reduce iteration counts for recycle-type problems and other slowly convergent systems. It is more specialized than Wegstein's method and is most useful when there is a clear dominant convergence mode.
+
+[fixpoint_anderson] is generally the most powerful solver in this family for difficult coupled systems. It combines several previous iterates through a least-squares acceleration step, often delivering substantially better convergence than direct substitution-based methods. The additional linear algebra introduces some overhead, but this cost is usually justified when evaluating $\mathbf{g}(\mathbf{x})$ is expensive or when simpler accelerators stall.
+
+[fixpoint_steffensen]: fixpoint_steffensen.md
+[fixpoint_damped]: fixpoint_damped.md
+[fixpoint_wegstein]: fixpoint_wegstein.md
+[fixpoint_dem]: fixpoint_dem.md
+[fixpoint_anderson]: fixpoint_anderson.md

docs/reference/math/optimization/OptimumResult.md

@@ -0,0 +1,6 @@
+# polykin.math.optimization
+
+::: polykin.math.optimization
+    options:
+        members:
+            - OptimumResult

docs/reference/math/optimization/VectorOptimumResult.md

@@ -0,0 +1,6 @@
+# polykin.math.optimization
+
+::: polykin.math.optimization
+    options:
+        members:
+            - VectorOptimumResult

docs/reference/math/optimization/_SUMMARY.md

@@ -0,0 +1,7 @@
+* [](index.md)
+* [fmin_brent](fmin_brent.md)
+* [fmin_nelder_mead](fmin_nelder_mead.md)
+* [fmin_qnewton](fmin_qnewton.md)
+* [fmin_secant](fmin_secant.md)
+* [OptimumResult](OptimumResult.md)
+* [VectorOptimumResult](VectorOptimumResult.md)

docs/reference/math/optimization/fmin_brent.md

@@ -0,0 +1,6 @@
+# polykin.math.optimization
+
+::: polykin.math.optimization
+    options:
+        members:
+            - fmin_brent

docs/reference/math/optimization/fmin_nelder_mead.md

@@ -0,0 +1,6 @@
+# polykin.math.optimization
+
+::: polykin.math.optimization
+    options:
+        members:
+            - fmin_nelder_mead

docs/reference/math/optimization/fmin_qnewton.md

@@ -0,0 +1,6 @@
+# polykin.math.optimization
+
+::: polykin.math.optimization
+    options:
+        members:
+            - fmin_qnewton

docs/reference/math/optimization/fmin_secant.md

@@ -0,0 +1,6 @@
+# polykin.math.optimization
+
+::: polykin.math.optimization
+    options:
+        members:
+            - fmin_secant

docs/reference/math/optimization/index.md

@@ -0,0 +1,22 @@
+# Optimization Solvers (polykin.math.optimization)
+
+## Univariate Optimization
+
+These solvers cover the typical trade-off between local speed and global robustness for minimization of scalar objective functions $f(x)$. They differ in the use of derivatives and whether a bounded interval is required.
+
+[fmin_secant] is a fast local method for smooth problems when two reasonable initial guesses are available. It uses a secant approximation of the first derivative, evaluated via centered finite differences or complex-step differentiation. It typically converges quickly near a well-behaved minimum but is not robust to poor initial guesses, weak curvature, or nonsmooth objectives.
+
+[fmin_brent] is a robust bounded method for cases where an interval containing a minimum is known. It combines golden-section search with inverse parabolic interpolation, achieving reliable convergence while exploiting smoothness when present. It does not require derivatives and is the preferred choice for safe one-dimensional minimization within a bracketed interval.
+
+## Multivariate Optimization
+
+These solvers provide different trade-offs between robustness, efficiency, and derivative requirements for unconstrained minimization of multivariate objective functions, $f(\mathbf{x})$.
+
+[fmin_nelder_mead] is a derivative-free method for unconstrained minimization of multivariate functions. It evolves a simplex of trial points through reflection, expansion, contraction, and shrink operations, making it useful when derivatives are unavailable, unreliable, or too expensive to compute. In exchange for this flexibility, convergence is usually slower than with gradient-based methods, and performance may degrade for poorly scaled or high-dimensional problems.
+
+[fmin_qnewton] is a fast local method for smooth unconstrained problems when gradient information is available analytically or can be estimated accurately. It combines Newton-type curvature information with optional BFGS Hessian updates and supports either a line-search or dogleg globalization strategy to improve robustness from remote initial guesses. This is generally the preferred solver when the objective is differentiable and a reasonably good starting point is available. As with other quasi-Newton methods, however, convergence remains fundamentally local.
+
+[fmin_secant]: fmin_secant.md
+[fmin_brent]: fmin_brent.md
+[fmin_nelder_mead]: fmin_nelder_mead.md
+[fmin_qnewton]: fmin_qnewton.md

docs/reference/math/roots/index.md

@@ -1 +1,24 @@
 # Root Solvers (polykin.math.roots)
+
+## Single Equation
+
+These three solvers span the usual tradeoff between local speed and global robustness for real-valued equations $f(x)=0$. They differ mainly in how much information they require and how sensitive they are to the initial guess.
+
+[root_newton] is the fastest option for smooth problems with a good initial guess. It uses a Newton–Raphson update and estimates derivatives via complex-step differentiation, so no user-supplied derivative is needed. Convergence is typically very rapid, but the function must accept complex inputs, and—as with any Newton method—it may fail or converge to the wrong root if initialized poorly.
+
+[root_secant] removes the need for derivatives by approximating the slope from two initial guesses. It retains fast local convergence on well-behaved problems and works when complex-step differentiation is not applicable. It is generally less robust than Brent’s method and less decisive than Newton near difficult roots, but provides a useful middle ground between simplicity and speed.
+
+[root_brent] is the safest choice when a bracketing interval is available. By combining bisection, secant steps, and inverse quadratic interpolation, it maintains the reliability of bracketed methods while recovering much of the speed of open methods. The only requirement is an interval where the function changes sign; when this is available, Brent’s method is typically the best choice for robust production use.
+
+## System of Equations
+
+[rootvec_qnewton] targets systems of nonlinear equations $\boldsymbol{F}(\boldsymbol{x})=\boldsymbol{0}$, which are inherently more challenging than scalar problems. There is no multidimensional analogue of bracketing, and thus no method combining global guarantees with high efficiency comparable to Brent’s method.
+
+Instead, such solvers rely on local linearization. Classical Newton methods use the full Jacobian, but forming and factorizing it can be expensive for large systems. `rootvec_qnewton` supports two approaches: recomputing the full Jacobian at each iteration (either user-supplied or via finite differences), or updating an approximate Jacobian using Broyden’s rank-1 formula. In this sense, Broyden’s method plays a role analogous to the secant method in the scalar case.
+
+Full Jacobian updates improve robustness at higher cost, while Broyden updates reduce per-iteration cost at the expense of reliability, especially with poor initial guesses or strong nonlinearity. As with other quasi-Newton methods, convergence is local, so good initial guesses remain important.
+
+[root_newton]: root_newton.md
+[root_secant]: root_secant.md
+[root_brent]: root_brent.md
+[rootvec_qnewton]: rootvec_qnewton.md

docs/tutorials/_SUMMARY.md

@@ -6,6 +6,7 @@
 * [Distributions](distributions.ipynb)
 * [Equations of State](equations_of_state.ipynb)
 * [Kinetic Coefficients](kinetic_coefficients.ipynb)
+* [Optimization](optimization.ipynb)
 * [Physical Property Equations](physprop_equations.ipynb)
 * [Thermal Conductivity](thermal_conductivity.ipynb)
 * [Viscosity](viscosity.ipynb)

mkdocs.yml

@@ -178,6 +178,7 @@ nav:
           - Derivatives: reference/math/derivatives/
           - Miscellaneous: reference/math/misc/
           - Fixed-Point Solvers: reference/math/fixpoint/
+          - Optimization Solvers: reference/math/optimization/
           - Root Solvers: reference/math/roots/
           - Special Functions: reference/math/special/
       - Properties:

pyproject.toml

@@ -13,7 +13,7 @@ classifiers = [
     "Topic :: Scientific/Engineering :: Chemistry",
 ]
 
-requires-python = ">=3.10.0"
+requires-python = ">=3.11.0"
 
 dependencies = [
     "matplotlib>=3.8",
@@ -34,6 +34,7 @@ dev = [
     "pytest-cov>=6.0.0",
     "pytest-xdist>=3.6.1",
     "ruff>=0.14.9",
+    "ty>=0.0.33",
 ]
 docs = [
     "black>=26.1.0",
@@ -58,7 +59,7 @@ addopts = ["-v", "-s", "-n=2"]
 
 [tool.ruff]
 line-length = 90
-target-version = "py310"
+target-version = "py311"
 
 [tool.ruff.lint]
 select = ["E", "F", "W", "I", "UP", "D", "NPY201"]
@@ -72,6 +73,15 @@ pydocstyle = { convention = 'numpy' }
 [tool.ruff.format]
 line-ending = "auto"
 
+[tool.ty.environment]
+root = ["./src"]
+
+[tool.ty.src]
+include = ["src", "tests"]
+
+[tool.ty.terminal]
+output-format = "concise"
+
 [build-system]
 requires = ["hatchling"]
 build-backend = "hatchling.build"

src/polykin/copolymerization/binary.py

@@ -8,8 +8,8 @@
 from numba import jit
 from scipy.integrate import solve_ivp
 
+from polykin.math.machine import eps
 from polykin.utils.exceptions import ODESolverError
-from polykin.utils.math import eps
 from polykin.utils.typing import FloatArray, FloatArrayLike, FloatMatrix, FloatVectorLike
 
 __all__ = [

src/polykin/copolymerization/fitting.py

@@ -471,7 +471,7 @@ def sse(r: tuple[float, float], atol=1e-4, rtol=1e-4) -> float:
     npar = 2
     if ndata > npar:
         s_sq = sse_opt / (ndata - npar)
-        H = hessian2_centered(sse, r_opt, h=1e-4)
+        H = hessian2_centered(sse, r_opt, epsf=1e-6)
         Hinv = np.linalg.inv(H)
         cov = 2 * Hinv * s_sq
     else:

src/polykin/copolymerization/multicomponent.py

@@ -8,8 +8,8 @@
 from numpy import exp
 from scipy.integrate import solve_ivp
 
+from polykin.math.machine import eps
 from polykin.utils.exceptions import ODESolverError, ShapeError
-from polykin.utils.math import eps
 from polykin.utils.typing import (
     FloatArray,
     FloatArrayLike,

src/polykin/copolymerization/penultimate.py

@@ -9,7 +9,7 @@
 import numpy as np
 
 from polykin.kinetics import Arrhenius
-from polykin.utils.math import eps
+from polykin.math.machine import eps
 from polykin.utils.tools import check_bounds
 from polykin.utils.typing import FloatArray, FloatArrayLike, IntArrayLike
 

src/polykin/copolymerization/terminal.py

@@ -15,8 +15,9 @@
 from scipy.integrate import solve_ivp
 
 from polykin.kinetics import Arrhenius
-from polykin.math import root_brent
-from polykin.utils.math import convert_FloatOrVectorLike_to_FloatVector, eps
+from polykin.math.machine import eps
+from polykin.math.roots import root_brent
+from polykin.utils.math import convert_FloatOrVectorLike_to_FloatVector
 from polykin.utils.tools import check_bounds, check_in_set, custom_repr
 from polykin.utils.typing import (
     FloatArray,

src/polykin/distributions/gold.py

@@ -7,7 +7,7 @@
 
 from polykin.distributions.analyticaldistributions import poisson
 from polykin.math import root_brent
-from polykin.utils.math import eps
+from polykin.math.machine import eps
 from polykin.utils.typing import FloatArray, IntArrayLike
 
 

src/polykin/hmt/diffusion.py

@@ -5,8 +5,8 @@
 from numpy import cos, exp, pi, sin, sqrt
 from scipy.special import erfc
 
+from polykin.math.machine import eps
 from polykin.math.special import ierfc, roots_xcotx, roots_xtanx
-from polykin.utils.math import eps
 
 __all__ = [
     "profile_constc_semiinf",
@@ -299,9 +299,7 @@ def uptake_constc_sheet(Fo: float) -> float:
     else:
         # Solution for normal times
         B = -(Fo * pi**2) / 4
-        S = sum(
-            1 / (2 * n + 1) ** 2 * exp(B * (2 * n + 1) ** 2) for n in range(0, N - 1)
-        )
+        S = sum(1 / (2 * n + 1) ** 2 * exp(B * (2 * n + 1) ** 2) for n in range(0, N - 1))
         result = 1 - (8 / pi**2) * S
 
     return result
@@ -433,9 +431,7 @@ def uptake_convection_sheet(Fo: float, Bi: float) -> float:
         N = 4  # Number of terms in series expansion (optimal value)
         x = roots_xtanx(Bi, N)
         b2 = x**2
-        S = sum(
-            exp(-b2[n] * Fo) / (b2[n] * (b2[n] + Bi * (Bi + 1))) for n in range(0, N)
-        )
+        S = sum(exp(-b2[n] * Fo) / (b2[n] * (b2[n] + Bi * (Bi + 1))) for n in range(0, N))
         return 1 - (2 * Bi**2) * S
     else:
         return 0
@@ -501,9 +497,7 @@ def uptake_convection_sphere(Fo: float, Bi: float) -> float:
         N = 4  # Number of terms in series expansion (optimal value)
         x = roots_xcotx(Bi, N)
         b2 = x**2
-        S = sum(
-            exp(-b2[n] * Fo) / (b2[n] * (b2[n] + Bi * (Bi - 1))) for n in range(0, N)
-        )
+        S = sum(exp(-b2[n] * Fo) / (b2[n] * (b2[n] + Bi * (Bi - 1))) for n in range(0, N))
         return 1 - (6 * Bi**2) * S
     else:
         return 0