# Implementation notes

**The below text is LLM-generated, based on the implementation. It has not
yet been reviewed and edited, and should not be trusted without checking
for correctness.**

equiconc solves for equilibrium concentrations using the convex optimization
approach of [Dirks et al. (2007), *SIAM Review* 49(1), 65--88](https://doi.org/10.1137/060651100).
This page summarizes the mathematical formulation.

## Problem setup

Consider a system with \(M\) monomer species and \(J\) total species (monomers +
complexes). Each species \(j\) has a concentration \(c_j\), and each complex is
defined by its stoichiometry: how many copies of each monomer it contains.

The **stoichiometry matrix** \(\mathbf{A} \in \mathbb{R}^{M \times J}\) has
entry \(A_{ij}\) equal to the number of copies of monomer \(i\) in species \(j\).
For monomer species, \(\mathbf{A}\) contains the \(M \times M\) identity matrix.

## Equilibrium conditions

At equilibrium, concentrations satisfy two constraints:

**Mass conservation.** The total amount of each monomer is conserved:

\[
\sum_{j=1}^{J} A_{ij}\, c_j = c_i^0 \quad \text{for } i = 1, \ldots, M
\]

where \(c_i^0\) is the total (initial) concentration of monomer \(i\).

**Chemical equilibrium.** Each complex concentration is related to the free
monomer concentrations by the partition function:

\[
c_j = Q_j \prod_{i=1}^{M} \left(\tilde{c}_i\right)^{A_{ij}}
\]

where \(\tilde{c}_i\) is the free concentration of monomer \(i\), and

\[
Q_j = \exp\!\left(-\frac{\Delta G_j^\circ}{RT}\right)
\]

is the Boltzmann factor for species \(j\) (with \(Q_i = 1\) for monomers,
i.e., \(\Delta G^\circ = 0\)).

## The dual problem

Rather than solving the primal (nonlinear mass conservation equations) directly,
equiconc minimizes the **convex dual** objective:

\[
f(\boldsymbol{\lambda}) = -\boldsymbol{\lambda}^\top \mathbf{c}^0
+ \sum_{j=1}^{J} Q_j \exp\!\left(\sum_{i=1}^{M} A_{ij}\, \lambda_i\right)
\]

where \(\boldsymbol{\lambda} \in \mathbb{R}^M\) are the dual variables
(one per monomer species). The key insight is that the dimension of this
optimization equals the number of monomer species (typically 2--10),
making it fast regardless of how many complexes exist.

### Gradient and Hessian

The gradient is:

\[
\nabla f(\boldsymbol{\lambda}) = -\mathbf{c}^0 + \mathbf{A}\, \mathbf{c}(\boldsymbol{\lambda})
\]

where \(c_j(\boldsymbol{\lambda}) = Q_j \exp(\mathbf{A}^\top \boldsymbol{\lambda})_j\).
At the optimum, \(\nabla f = 0\) recovers the mass conservation equations.

The Hessian is:

\[
\mathbf{H}_f = \mathbf{A}\, \text{diag}(\mathbf{c})\, \mathbf{A}^\top
\]

Since all \(c_j > 0\), the Hessian is positive definite, guaranteeing that
\(f\) is strictly convex and has a unique minimum.

## Trust-region Newton method

equiconc uses a **trust-region Newton method** with **dog-leg steps**:

1. At each iteration, compute the Newton step
   \(\mathbf{p}_N = -\mathbf{H}^{-1} \nabla f\) via Cholesky decomposition
   (always succeeds for the linear objective since \(\mathbf{H}_f \succ 0\)).

2. If \(\|\mathbf{p}_N\| \leq \Delta\) (trust-region radius), take the full
   Newton step.

3. Otherwise, compute the Cauchy (steepest descent) step and interpolate along
   the dog-leg path to stay within the trust region.

4. Update the trust-region radius based on the ratio of actual to predicted
   reduction.

## Objective surface: linear vs. log

equiconc exposes two objective surfaces for the trust-region step,
selected via `SolverOptions.objective`:

- **Linear** (default) minimizes \(f(\boldsymbol\lambda)\) directly. The
  Hessian \(\mathbf H_f = \mathbf A\, \text{diag}(\mathbf c)\, \mathbf A^\top\)
  is positive semi-definite by construction, so plain Cholesky always
  succeeds.
- **Log** minimizes the transformed objective
  \[
  g(\boldsymbol\lambda) = \ln f(\boldsymbol\lambda).
  \]
  Because \(f > 0\) at any feasible iterate (and at the optimum), \(g\)
  shares the minimizer of \(f\) but compresses its exponential dynamic
  range, which can dramatically reduce the iteration count on stiff
  systems (very strong binding, asymmetric \(\mathbf c^0\)). Its gradient
  and Hessian are
  \[
  \nabla g = \frac{\nabla f}{f}, \qquad
  \mathbf H_g = \frac{\mathbf H_f}{f} - \frac{(\nabla f)(\nabla f)^\top}{f^2}.
  \]
  The rank-one subtraction makes \(\mathbf H_g\) potentially **indefinite
  away from the optimum**, so \(g\) is *not* globally convex. equiconc
  compensates with on-the-fly diagonal **modified Cholesky**
  regularization, guards \(f > 0\) at every iterate, and rejects any
  trial step whose quadratic model predicts an ascent.

Both paths share the dual variable \(\boldsymbol\lambda\), the same
primal recovery \(c_j = \exp(\log Q_j + (\mathbf A^\top\boldsymbol\lambda)_j)\),
and the same convergence test on the primal mass-conservation residual,
so results are interchangeable to within tolerance.

### Log-space evaluation

To prevent floating-point overflow when concentrations span many orders of
magnitude, the objective, gradient, and Hessian are evaluated in log space:

\[
\log c_j = \log Q_j + \sum_{i} A_{ij}\, \lambda_i
\]

Values of \(\log c_j\) are clamped to prevent \(\exp(\cdot)\) from exceeding
`f64::MAX`.

## Recovering primal concentrations

Once the dual optimum \(\boldsymbol{\lambda}^*\) is found, primal
concentrations are recovered as:

\[
c_j^* = \exp\!\left(\log Q_j + (\mathbf{A}^\top \boldsymbol{\lambda}^*)_j\right)
\]

## References

- Dirks, R. M., Bois, J. S., Schaeffer, J. M., Winfree, E., & Pierce, N. A.
  (2007). Thermodynamic analysis of interacting nucleic acid strands.
  *SIAM Review*, 49(1), 65--88.
  [doi:10.1137/060651100](https://doi.org/10.1137/060651100)
