[cjkunselman18]

This pull request uses Jansson derivatives to calculate and output the gradient of each contribution to the likelihood function for use in gradient-based optimization.

[bocklund]

some changes I want to make and I think it's about ready to merge:

1. we copy of the Database objects in context creation to prevent mutating the variables in the global state to prevent weird behavior or race conditions with the lasting Database object. PyCalphad's Database.__deepcopy__ shallow copies symbols and we might need a more detailed fix/check for correct behavior.
2. for performance, since gradients are a performance regression in function call time, i'd like to make the gradients behavior opt-in so that gradients are only computed if requested by the calling optimizer. this also has the benefit of allowing us to keep the test changes small and not thread through a bunch of API changes to handle gradient outputs
3. double check whether sign error in gradients for SciPy optimizers
4. investigate use of adjoint method to reduce computational time of gradients. would require more implementation work and maybe some pycalphad changes to support, but should make it significantly cheaper to evaluate gradients. it could also solve the problem with parameters-as-conditions and hitting the hard NumPy limit of too many dimensions in an array (which we hit because PyCalphad encodes each parameter-condition as a condition array dimension)

[bocklund]

I have been thinking on this awhile and made several prototypes. I think I finally have something worth sharing. I've pushed a bunch of changes. The original intent of being able to compute gradients, of course, stays.

Changes

What I've done is to:

• Allow parameters to stay as symbols rather than as conditions. This fixes a issues with things like having hard limits on number of parameters (because PyCalphad was making them array dimensions and NumPy is limited to 64 dimensional arrays). We don't need to copy or modify databases anymore. There's now no global state so parallelization should work fine.
• Make gradients opt-in. This way we don't have to change the tests and don't have to pay for gradients if we don't want to use them, BUT
• (big change) we now use a different method to get gradients, which makes the gradients go from 5-20x slower than likelihood only to only about 10% slower (1.1x slower) than likelihood only. The previous revision did one JanssonDerivative calculation per (output, parameter) pair, rebuilding the solver state each time. Now each converged equilibrium's tangent (dμ, dNP, dy) with respect to every parameter comes from one least-squares solve of PyCalphad's equilibrium matrix. Each parameter's right-hand side is the column it would have as a free state variable. This is espei/_parameter_sensitivities.py, in NumPy, with no PyCalphad changes.
• fixed a bug introduced in this PR where ZPF data from a subsystem without some of the fitted parameters failed. Binary ZPF data in a ternary fit raised "Conditions not yet supported" in starting_point, because each ZPF dataset's phase records only know the parameters of its subsystem, while every fitted parameter was passed as a condition. Master handles this case; a regression test is added
• (existing bug) _FORM fixed-configuration data ignored excluded_model_contributions for the phase; only the reference states excluded them. shift_reference_state builds HMR/SMR/CPMR before the contributions were zeroed.
• (existing bug) ZPFResidual and EquilibriumPropertyResidual read the parameter vector in sorted-name order whatever the order of symbols_to_fit, while the activity and fixed-configuration residuals followed symbols_to_fit. All four now follow symbols_to_fit. setup_context sorts, so ESPEI's own fits were not affected.

How to use it

Currently there's no YAML API. You need to use the Python API.

Here's an example to run:

# Fit the Cu-Mg tutorial's generated parameters with L-BFGS-B and likelihood gradients.
# Run from ESPEI's docs/tutorials/cu-mg-example directory. Needs PyCalphad develop .
import time
import timeit
import numpy as np
from pycalphad import Database, binplot, variables as v
from espei.datasets import apply_tags, load_datasets, recursive_glob
from espei.error_functions.context import setup_context
from espei.optimizers.opt_scipy import SciPyOptimizer
from espei.plot import dataplot
from espei.utils import database_symbols_to_fit, unpack_piecewise

# The tutorial's datasets and tags (run_mcmc.yaml) and generated parameters
datasets = load_datasets(recursive_glob("input-data"))
apply_tags(datasets, {"dft": {"excluded_model_contributions": ["idmix", "mag"]},
                      "estimated-entropy": {"excluded_model_contributions": ["idmix", "mag"]}})
dbf = Database("Cu-Mg-generated.tdb")
symbols = database_symbols_to_fit(dbf)

# Time each residual's likelihood with and without its gradient
ctx = setup_context(dbf, datasets, symbols)
x = np.array([unpack_piecewise(dbf.symbols[s]) for s in symbols])
for residual in ctx["residual_objs"]:
    residual.get_likelihood_and_gradient(x)  # the first call compiles the parameter derivatives
    t_likelihood = min(timeit.repeat(lambda: residual.get_likelihood(x), number=1, repeat=3))
    t_gradient = min(timeit.repeat(lambda: residual.get_likelihood_and_gradient(x), number=1, repeat=3))
    print(f"{type(residual).__name__}: likelihood {t_likelihood:.3f} s, with gradient {t_gradient:.3f} s")

# Maximize the likelihood with its gradient
opt = SciPyOptimizer(dbf)
t0 = time.perf_counter()
fit_dbf = opt.fit(symbols, datasets, method="L-BFGS-B", jac=True, maxiter=50)
print(f"L-BFGS-B: {time.perf_counter() - t0:.1f} s, {opt.minimize_result.nfev} evaluations, lnL {-opt.minimize_result.fun:.1f}")
fit_dbf.to_file("Cu-Mg-gradient.tdb", if_exists="overwrite")

comps = ["CU", "MG", "VA"]
phases = list(fit_dbf.phases.keys())
conds = {v.P: 101325, v.T: (500, 1500, 10), v.X("MG"): (0, 1, 0.02)}
ax = binplot(fit_dbf, comps, phases, conds)
dataplot(comps, phases, conds, datasets, ax=ax)
ax.figure.savefig("Cu-Mg-gradient.png", bbox_inches="tight")

However, it's worth noting that we have had better results in our papers [1,2] if we scale the parameters to order unity:

# Fit the Cu-Mg tutorial's generated parameters with L-BFGS-B in parameters scaled to order unity.
# Run from ESPEI's docs/tutorials/cu-mg-example directory. Needs PyCalphad develop
import time
import numpy as np
from scipy.optimize import minimize
from pycalphad import Database, binplot, variables as v
from espei.datasets import apply_tags, load_datasets, recursive_glob
from espei.error_functions.context import setup_context
from espei.plot import dataplot
from espei.utils import database_symbols_to_fit, unpack_piecewise

# The tutorial's datasets and tags (run_mcmc.yaml) and generated parameters
datasets = load_datasets(recursive_glob("input-data"))
apply_tags(datasets, {"dft": {"excluded_model_contributions": ["idmix", "mag"]},
                      "estimated-entropy": {"excluded_model_contributions": ["idmix", "mag"]}})
dbf = Database("Cu-Mg-generated.tdb")
ctx = setup_context(dbf, datasets, database_symbols_to_fit(dbf))
symbols = ctx["symbols_to_fit"]

# Optimize z = theta / |theta_0|, so every variable starts at +-1
theta0 = np.array([unpack_piecewise(dbf.symbols[s]) for s in symbols])
scale = np.where(theta0 != 0, np.abs(theta0), 1.0)

def objective(z):
    """Negative log-likelihood and its gradient with respect to the scaled parameters"""
    lnlike, gradient = 0.0, np.zeros(z.size)
    for residual in ctx["residual_objs"]:
        likelihood, residual_gradient = residual.get_likelihood_and_gradient(z * scale)
        lnlike += likelihood
        gradient += residual_gradient
    return -lnlike, -gradient * scale  # dlnL/dz = |theta_0| dlnL/dtheta

t0 = time.perf_counter()
result = minimize(objective, theta0 / scale, jac=True, method="L-BFGS-B", options={"maxiter": 50})
print(f"L-BFGS-B: {time.perf_counter() - t0:.1f} s, {result.nfev} evaluations, lnL {-result.fun:.1f}")

# Fitted symbols are parameters to the residuals, so updating the Database afterwards is safe
dbf.symbols.update(dict(zip(symbols, result.x * scale)))
dbf.to_file("Cu-Mg-gradient-scaled.tdb", if_exists="overwrite")

comps = ["CU", "MG", "VA"]
phases = list(dbf.phases.keys())
conds = {v.P: 101325, v.T: (500, 1500, 10), v.X("MG"): (0, 1, 0.02)}
ax = binplot(dbf, comps, phases, conds)
dataplot(comps, phases, conds, datasets, ax=ax)
ax.figure.savefig("Cu-Mg-gradient-scaled.png", bbox_inches="tight")

it takes about a minute to produce:

[1] C. Kunselman, B. Bocklund, R. Otis, R. Arróyave, Gradient-informed Hamiltonian Monte Carlo for multicomponent CALPHAD model optimization and uncertainty quantification, Acta Materialia (2026) 122412. https://doi.org/10.1016/j.actamat.2026.122412.
[2] C. Kunselman, B. Bocklund, R. Otis, R. Arróyave, Analytical gradient-based optimization of CALPHAD model parameters, Acta Materialia 300 (2025) 121453. https://doi.org/10.1016/j.actamat.2025.121453.

Performance

Times are per evaluation, medians, and the central-difference column is 2·N_θ likelihoods:

	Likelihood	Likelihood + gradient	Previous revision	Central differences
Cu–Mg tutorial (22 parameters, all four data types)	1.01 s	1.12 s (1.11×)	24.5 s	44.5 s
Cr–Fe–Ti–V, 10 parameters	1.78 s	1.94 s (1.09×)	fails (item 1)	35.6 s
Cr–Fe–Ti–V, 200 parameters	2.59 s	2.83 s (1.09×)	over its limit	1036 s

Per equilibrium, the tangent costs 0.2–0.5 ms at 10–200 parameters, against a 10–14 ms solve.

What's left to do:

I might also try to take a pass at simplifying the NumPy backend. Ultimately, we could move this tangent approach back into PyCalphad and it would be a bit faster, but

1. changes would be more invasive (about 2-3x more lines of code to put it in PyCaplhad's Cython core),
2. this is a bit experimental and I don't want to commit to those invasive changes yet,
3. we should be able to migrate cleanly if/when this get ported to PyCalphad

These require the latest PyCalphad develop branch, but no other changes on PyCalphad. It will be compatible with PyCalphad 0.12 when that is released.