ChemistryLab.jl

Appearance

Utilities ​

Module ​

ChemistryLab.ChemistryLab Module
julia
ChemistryLab

Top-level module for parsing, representing and manipulating chemical formulas, species, stoichiometric matrices and ThermoFun / PHREEQC-like data (Parkhurst and Appelo, 2013).

Overview

  • Canonical containers: AtomGroup, Formula.

  • Species representations: Species, CemSpecies.

  • Parsers and IO helpers for ThermoFun / PHREEQC data (substances, reactions).

  • Stoichiometric matrix construction and reaction helpers.

  • Utilities for units, colored and Unicode terminal output.

Examples

julia
julia> using ChemistryLab

julia> f = Formula("H2O")
Formula{Int64}
    formula: H2O ◆ H₂O
composition: H => 2, O => 1
     charge: 0

julia> f[:H]
2

julia> H2O = Species(f; name="Vapour", symbol="H₂O", aggregate_state=AS_GAS, class=SC_GASFLUID)
Species{Int64}
           name: Vapour
         symbol: H₂O
        formula: H2O ◆ H₂O
          atoms: H => 2, O => 1
         charge: 0
aggregate_state: AS_GAS
          class: SC_GASFLUID
     properties: M = 0.0180149999937744 kg mol⁻¹

julia> CO2 = Species(
           Dict(:C=>1, :O=>2);
           name="Carbon dioxide",
           symbol="CO₂",
           aggregate_state=AS_GAS,
           class=SC_GASFLUID,
       ) # definition from Dict
Species{Int64}
           name: Carbon dioxide
         symbol: CO₂
        formula: CO2 ◆ CO₂
          atoms: O => 2, C => 1
         charge: 0
aggregate_state: AS_GAS
          class: SC_GASFLUID
     properties: M = 0.04400899998479143 kg mol⁻¹

julia> O2 = Species("O2"; name="Dioxygen", symbol="O₂", aggregate_state=AS_GAS, class=SC_GASFLUID) # definition from String
Species{Int64}
           name: Dioxygen
         symbol: O₂
        formula: O2 ◆ O₂
          atoms: O => 2
         charge: 0
aggregate_state: AS_GAS
          class: SC_GASFLUID
     properties: M = 0.03199799998894218 kg mol⁻¹

julia> C3H8 = Species(
           "C₃H₈"; name="Propane", symbol="C₃H₈", aggregate_state=AS_GAS, class=SC_GASFLUID
       ) # definition from Unicode String
Species{Int64}
           name: Propane
         symbol: C₃H₈
        formula: C₃H₈ ◆ C3H8
          atoms: C => 3, H => 8
         charge: 0
aggregate_state: AS_GAS
          class: SC_GASFLUID
     properties: M = 0.04409699998476102 kg mol⁻¹

julia> r = Reaction([C3H8, O2, CO2, H2O])
  equation: C₃H₈ + 5O₂ = 4H₂O + 3CO₂
 reactants: C₃H₈ => 1, O₂ => 5
  products: H₂O => 4, CO₂ => 3
    charge: 0

Physical constants ​

The gas constant, the Faraday constant and the electric constant, taken from DynamicQuantities rather than written down, each in two forms: the bare SI number for arithmetic and the dimensioned quantity for anything that must carry its unit. See Where the numbers come from for how to reach these, and everything else the package already knows, from a script.

ChemistryLab.AVOGADRO Constant
julia
AVOGADRO

The Avogadro constant as a plain Float64 in mol⁻¹, derived from AVOGADRO_Q.

ChemistryLab.AVOGADRO_Q Constant
julia
AVOGADRO_Q

The Avogadro constant with its dimensions, 6.02214076e23 mol⁻¹, from DynamicQuantities.Constants.N_A (exact by the 2019 SI redefinition).

Surface chemistry states site densities per square nanometer and the rest of this package works in moles per square meter; this is the only bridge between them, so it is taken from the library rather than written down.

ChemistryLab.CALORIE Constant
julia
CALORIE

The thermochemical calorie as a plain Float64 in joules, derived from CALORIE_Q: 4.184.

ChemistryLab.CALORIE_Q Constant
julia
CALORIE_Q

The thermochemical calorie with its dimensions, 4.184 J by definition: the unit of the energies of the older databases (kcal in PHREEQC files, cal in SUPCRT). DynamicQuantities does not define it.

ChemistryLab.FARADAY Constant
julia
FARADAY

The Faraday constant as a plain Float64 in C/mol, derived from FARADAY_Q. Used by the Nernst relation between pe and Eh.

ChemistryLab.FARADAY_Q Constant
julia
FARADAY_Q

The Faraday constant with its dimensions, 96485.332… C/mol, from DynamicQuantities.Constants.F (CODATA, exact by the 2019 SI redefinition).

ChemistryLab.P_STANDARD Constant
julia
P_STANDARD

The standard-state pressure as a plain Float64 in pascals, derived from P_STANDARD_Q: 1.0e5.

ChemistryLab.P_STANDARD_Q Constant
julia
P_STANDARD_Q

The standard-state pressure with its dimensions, 1 bar: the pressure at which the standard properties of the databases read by this package are tabulated, and the P° of an ideal gas's activity xᵢ P/P°. It is a convention, IUPAC's since 1982, and not a measured constant, which is why it does not come from DynamicQuantities.Constants.

ChemistryLab.R_GAS Constant
julia
R_GAS

The molar gas constant as a plain Float64 in J/(mol·K), derived from R_GAS_Q rather than written again.

This is the form the inner loops want: an activity model is evaluated once per species per Newton iteration, on Float64 and on ForwardDiff.Dual, and a dimensional Quantity there costs an allocation each time for a factor whose units are known at the call site anyway.

ChemistryLab.R_GAS_Q Constant
julia
R_GAS_Q

The molar gas constant with its dimensions, 8.314462618… J/(mol·K), taken from DynamicQuantities.Constants.R (CODATA).

Use this wherever the surrounding quantities carry units. For the plain number used inside an activity model or a rate law, see R_GAS.

ChemistryLab.T_STANDARD Constant
julia
T_STANDARD

The reference temperature as a plain Float64 in kelvins, derived from T_STANDARD_Q: 298.15.

ChemistryLab.T_STANDARD_Q Constant
julia
T_STANDARD_Q

The reference temperature with its dimensions, 298.15 K: the temperature at which the databases read by this package tabulate their standard properties and their equilibrium constants. It is a convention, as P_STANDARD_Q is.

ChemistryLab.T_ZERO_CELSIUS Constant
julia
T_ZERO_CELSIUS

The zero of the Celsius scale as a plain Float64 in kelvins, derived from T_ZERO_CELSIUS_Q: 273.15.

ChemistryLab.T_ZERO_CELSIUS_Q Constant
julia
T_ZERO_CELSIUS_Q

The zero of the Celsius scale with its dimensions, 273.15 K, by the definition of the scale. DynamicQuantities has no affine units, and a temperature in degrees Celsius is converted to kelvins by adding it.

ChemistryLab.VACUUM_PERMITTIVITY Constant
julia
VACUUM_PERMITTIVITY

The electric constant as a plain Float64 in F/m, derived from VACUUM_PERMITTIVITY_Q. Used by the Gouy-Chapman relation between a surface charge density and the potential it raises, in DiffuseLayer.

ChemistryLab.VACUUM_PERMITTIVITY_Q Constant
julia
VACUUM_PERMITTIVITY_Q

The electric constant with its dimensions, 8.8541878128e-12 F/m, from DynamicQuantities.Constants.eps_0 (CODATA).

ChemistryLab.RT_over_F Method
julia
RT_over_F(T)

RT/F in volts at temperature T in kelvin — the natural scale of an electrochemical potential, and the factor in Nernst's relation. At 25 °C, RT ln(10)/F is the familiar 0.05916 V per unit of pe.

Numerical floors ​

The small amounts the solvers and the activity models work with, each defined once with the reason for its value: the default ϵ of the solvers, the floor of the activities, the floor of the certificate and the exponent below which an amount is taken as zero.

ChemistryLab._ACTIVITY_FLOOR Constant
julia
_ACTIVITY_FLOOR

The molality, in mol per kilogram of water, below which an activity model reads a species at this molality: 1e-30. It keeps log n finite, and nothing else. The amount it stands for is that molality times the water of the state (_activity_floor_amount), so that the floor follows the size of the system: until 0.41.0 it was 1e-30 mol whatever the system, the same amount for a kilogram of paste and for a pore solution of a few grams.

It is lower than ϵ = 1e-16, the floor that bounds the interior-point back ends and fills the absent products of a cold state, which is too high a floor for the activities of a cement paste. With a few grams of water per 100 g of binder, H⁺ at pH 13.5 to 14 is 1e-16 mol: floored there, its activity no longer followed its amount, a solve left it at 3e-100 mol, and a pH read from the amount came out 0.09 high.

Other codes floor lower, on one recurring argument. Leal et al. (2017) recall that in an interior-point minimization an unstable species ends at an amount of the order of the barrier parameter, and that Leal et al. (2016) recommended it below 1e-25 so that such a species holds less than one molecule (1/N_A = 1.66e-24 mol) in a system of one mole. GEMS3K (Kulik et al., 2013) takes that same molecule as the least amount it considers (lowPosNum, 1.66e-24 mol) and eliminates a solution species below DcMin = 1e-30 mol (ms_multi.h). PHREEQC sets a molality to zero belowMIN_LM = −30in log (global_structures.h), the form taken here. Reaktoro bounds every amount byepsilon = 1e-16mol, a bound of its solver from which it takes its barrier parameter (EquilibriumOptions.hpp), the role ourϵhas. Soϵ keeps its other two roles, and the activity models floor at this molality.

ChemistryLab._AMOUNT_FLOOR Constant
julia
_AMOUNT_FLOOR

The default of the keyword ϵ, 1e-16 mol: the smallest amount the interior-point back ends are bounded by, the amount of an absent product in a cold state, and the size of the regularizations that keep a square root or a ratio finite at zero (sqrt(I + ϵ)). It is no longer the floor of the activities; see _ACTIVITY_FLOOR.

ChemistryLab._CERTIFICATE_FLOOR Constant
julia
_CERTIFICATE_FLOOR

The default floor of optimality_certificate, 1e-25 mol: an amount above it is interior, and held to the stationarity equality; below it a pure phase is held to its inequality and a member of a present phase to the one-sided form of the equality (OptimaSolver's kkt_certificate, whose default it is too).

ChemistryLab._LOG_UNDERFLOW Constant
julia
_LOG_UNDERFLOW

The exponent below which an amount computed as exp(w) is taken as the truncation exp(-700) ≈ 1e-304: the bound OptimaSolver's dual Newton clamps its log amounts to (W_FLOOR), and the symmetric bound under which exp does not overflow.

ChemistryLab._activity_floor Method
julia
_activity_floor(p) -> Real

The floor an activity model applies to the amounts of p: p.ϵa, set by _build_params to min(ϵ, _activity_floor_amount(state)), or p.ϵ for a parameter tuple built by hand without it.

Miscellaneous helpers ​

ChemistryLab.COL_CHARGE Constant

Terminal color for charge display (cyan bold).

ChemistryLab.COL_PAR Constant

Terminal color for parentheses (magenta bold).

ChemistryLab.COL_STOICH_EXT Constant

Terminal color for external stoichiometric coefficients (yellow bold).

ChemistryLab.COL_STOICH_INT Constant

Terminal color for internal stoichiometric coefficients (red bold).

ChemistryLab._TPMemo Type
julia
_TPMemo{K, V}

A memo of a pure function of the temperature and the pressure, asked for at the same arguments by every species and every state an evaluation builds (the density of water, a Newton solve of its equation of state, and what is derived from it). On plain numbers it keeps a table of every call, emptied past _MEMO_MAX entries; on dual numbers it keeps the last call, matched by identity of its arguments (===: values and partials). Both under a lock. What it returns is what the function computes, bit for bit.

ChemistryLab._adim Method
julia
_adim(x) -> Real

x as a bare number, checking that it carries no dimension.

Why not extend the Base math functions to quantities

Extending them to DynamicQuantities.Quantity by stripping the unit first,

julia
Base.log(x::Quantity) = log(ustrip(x))

would be type piracy — the function and the type both belong elsewhere — and so apply to every package loaded beside this one, whether or not it asked. And it would not merely remove a dimension check; it would make the result depend on an invisible normalization, because ustrip returns the value in SI base units:

julia
log(2u"m")       = 0.693      instead of raising
log(1u"mol/L")   = 6.908      the same quantity written two ways,
log(1u"mmol/L")  = 0.0        two different answers

What this does instead

Dispatch rather than a branch, so Float64, ForwardDiff.Dual and Symbolics.Num all take the first method untouched — none of them carries a dimension — and only a genuine Quantity pays for the test.

A dimensionless Quantity passes: a ratio whose units cancel is a number, and saying so explicitly is the point. Anything else raises.

ChemistryLab._ensure_unit Method
julia
_ensure_unit(unit, x::Real) -> Quantity
_ensure_unit(unit, x) -> Quantity

Ensure a value is a Quantity with the given unit.

  • Plain Real → interpreted as SI, wrapped: x * unit.

  • Quantity → converted to unit via safe_uconvert.

ChemistryLab._primal Method
julia
_primal(x::Real) -> Real

Extract the primal (non-derivative) value of x for use in control-flow decisions.

For plain Real numbers this is a no-op. For forward-mode AD types such as ForwardDiff.Dual (which are subtypes of Real) this returns x itself, but since those types overload comparison operators to compare only the primal value, writing _primal(x) > threshold makes the intent explicit and ensures that branching decisions are always taken on physical values rather than derivatives.

Adding a method _primal(x::YourDualType) = x.value enables support for non-standard AD frameworks that do not overload comparison operators.

Examples

julia
julia> _primal(3.14)
3.14
ChemistryLab.force_uconvert Method
julia
force_uconvert(qout::UnionAbstractQuantity, q) -> AbstractQuantity

Strip units from q (converting to qout dimensions when possible) and re-attach the unit of qout. Unlike safe_uconvert, this always returns a quantity with the units of qout, even when the input is dimensionless.

ChemistryLab.print_title Method
julia
print_title(title; crayon=:none, indent="", style=:none)

Print a formatted title with optional styling and indentation.

Arguments

  • title: title string to display.

  • crayon: Crayon for colored output (default :none).

  • indent: left indentation string (default "").

  • style: formatting style - :none, :underline, or :box (default :none).

Examples

julia
julia> print_title("Test Title"; style=:none)
Test Title

julia> print_title("Section"; style=:underline, indent="  ")
  Section
  ───────
ChemistryLab.root_type Method
julia
root_type(T) -> Type

Get the root type wrapper from a potentially parameterized type.

ChemistryLab.safe_uconvert Method
julia
safe_uconvert(qout::UnionAbstractQuantity{<:Any, <:AbstractSymbolicDimensions}, q::UnionAbstractQuantity{<:Any, <:Dimensions}) -> UnionAbstractQuantity
safe_uconvert(qout::UnionAbstractQuantity{<:Any, <:AbstractSymbolicDimensions}, q) -> Any

Convert q to the unit type represented by qout using uconvert when q is a quantity with compatible dimensions. If q is a plain numeric value, return q unchanged.

Arguments

  • qout : unitful quantity target (e.g., u"m", u"kg").

  • q : a quantity or a plain numeric value.

Returns

  • Converted quantity (of type qout) or the original q when q has no units.

Examples

julia
julia> safe_uconvert(us"m", 5u"cm")
0.05 m

julia> safe_uconvert(us"m", 3.0)
3.0
ChemistryLab.safe_ustrip Method
julia
safe_ustrip(unit::UnionAbstractQuantity, q::UnionAbstractQuantity) -> Number
safe_ustrip(unit::UnionAbstractQuantity, q) -> Number

Safely remove units from q, converting to unit when their dimensions match. If q is a plain numeric value (no units) it is returned unchanged.

When the dimensions of unit and q differ (e.g. unit is dimensionless but q carries a physical dimension) the conversion is skipped and q is stripped in its native unit. This prevents a DimensionError in contexts such as ThermoFactory where a dimensionless placeholder unit u"1" may be passed alongside a quantity with physical dimensions.

Examples

julia
julia> safe_ustrip(u"m", 5u"m")
5.0

julia> safe_ustrip(u"cm", 1u"m")
100.0

julia> safe_ustrip(u"m", 3.0)
3.0

Sub/superscripts ​

ChemistryLab.dict_all_normal_to_sub Constant

Extended mapping from normal-line characters (digits, letters, operators) to Unicode subscripts.

ChemistryLab.dict_frac_unicode Constant

Mapping from Rational fractions to their Unicode vulgar fraction strings.

ChemistryLab.dict_normal_to_sub Constant

Mapping from normal-line digits and signs to their Unicode subscript equivalents.

ChemistryLab.dict_normal_to_super Constant

Mapping from normal-line characters to their Unicode superscript equivalents.

ChemistryLab.dict_sub_to_normal Constant

Mapping from Unicode subscript digits to their normal-line equivalents.

ChemistryLab.dict_super_to_normal Constant

Mapping from Unicode superscript characters to their normal-line equivalents.

ChemistryLab.dict_unicode_frac Constant

Mapping from Unicode vulgar fraction characters to Rational values.

ChemistryLab.all_normal_to_sub Method
julia
all_normal_to_sub(s::AbstractString) -> String

Convert all normal characters (including letters and operators) to subscripts.

ChemistryLab.from_subscriptnumber Method
julia
from_subscriptnumber(s::String) -> Int

Parse a Unicode subscript number string to an integer.

Examples

julia
julia> from_subscriptnumber("₄₂")
42

julia> from_subscriptnumber("₋₃")
-3
ChemistryLab.from_superscriptnumber Method
julia
from_superscriptnumber(s::String) -> Int

Parse a Unicode superscript number string to an integer.

Examples

julia
julia> from_superscriptnumber("⁴²")
42

julia> from_superscriptnumber("⁻²")
-2
ChemistryLab.issubscript Method
julia
issubscript(c::Char) -> Bool

Return whether c is a numeric subscript.

Examples

julia
julia> issubscript('₂')
true

julia> issubscript('2')
false
ChemistryLab.issuperscript Method
julia
issuperscript(c::Char) -> Bool

Return whether c is a numeric superscript or ⁺/⁻.

Examples

julia
julia> issuperscript('²')
true

julia> issuperscript('a')
false
ChemistryLab.normal_to_sub Method
julia
normal_to_sub(s::AbstractString) -> String

Convert all normal characters in s to numeric subscripts.

Examples

julia
julia> normal_to_sub("H2O")
"H₂O"
ChemistryLab.normal_to_super Method
julia
normal_to_super(s::AbstractString) -> String

Convert all normal characters or +/- in s to numeric superscripts.

Examples

julia
julia> normal_to_super("2+")
"²⁺"
ChemistryLab.sub_to_normal Method
julia
sub_to_normal(s::AbstractString) -> String

Convert all numeric subscripts in s to normal line.

Examples

julia
julia> sub_to_normal("H₂O")
"H2O"
ChemistryLab.subscriptnumber Method
julia
subscriptnumber(i::Integer) -> String

Convert an integer to its Unicode subscript representation.

Examples

julia
julia> subscriptnumber(42)
"₄₂"

julia> subscriptnumber(-3)
"₋₃"
ChemistryLab.super_to_normal Method
julia
super_to_normal(s::AbstractString) -> String

Convert all numeric superscripts or ⁺/⁻ in s to normal line.

Examples

julia
julia> super_to_normal("Ca²⁺")
"Ca2+"
ChemistryLab.superscriptnumber Method
julia
superscriptnumber(i::Integer) -> String

Convert an integer to its Unicode superscript representation.

Examples

julia
julia> superscriptnumber(42)
"⁴²"

julia> superscriptnumber(-2)
"⁻²"

Provenance ​

A number that says how it was obtained and from where, so the claim survives out of a data file and into a table, a figure or a fitted result. See Where the numbers come from for when to reach for it.

ChemistryLab.ProvenanceKind Type
julia
ProvenanceKind

How a number came to be, ordered from the weakest claim to the strongest.

  • PROV_PLACEHOLDER — a value standing in for one nobody has yet. It exists so a calculation can be demonstrated, and it is not evidence of anything.

  • PROV_ESTIMATED — borrowed from an analog or inferred from a correlation. Somebody reasoned, nobody measured.

  • PROV_FITTED — identified from data by an optimization. As sound as the data, the model and the identifiability of the parameter, which is three separate questions.

  • PROV_PUBLISHED — taken from a named source that reports it as a value.

  • PROV_MEASURED — determined experimentally in the work at hand.

  • PROV_UNSTATED — nothing is claimed. Deliberately the weakest, so a number that forgot to say where it came from never strengthens a result.

The order is what weakest uses, and the ordering is the point: a result is only as sound as its least sound input.

ChemistryLab.Traced Type
julia
struct Traced{T}

A value carrying how it was obtained and from where.

Traced is not a Real. That is deliberate: a number that flowed silently into arithmetic would arrive at the far end with its provenance gone, and the whole point is that losing it should be an act rather than an accident. Unwrap it with value at the boundary where it enters a calculation, which is where a reader can see the claim being dropped.

Fields

  • value: the number.

  • kind: a ProvenanceKind.

  • source: free text — a database and version, a paper, the analog a value was borrowed from, the dataset it was fitted to.

  • uncertainty: what the source says about how well it is known, in the same unit as value, or nothing when it says nothing. Published sorption compilations report it per entry — ClaySor 2023 (Marinich et al., 2025) writes error: 0.17 beside a log K — and a compilation where half the constants are known to 0.1 and half to 0.5 is not one number's worth of information.

Example

julia
julia> using ChemistryLab

julia> ε = Traced(0.05, PROV_ESTIMATED, "analog: SO4-2");

julia> ChemistryLab.value(ε), provenance(ε)
(0.05, PROV_ESTIMATED)

julia> is_evidence(ε)
false

value and source are not exported, and reaching for them through the module is not an oversight. Both names are taken: Symbolics and SciMLBase each export a value, so exporting one here would make the bare name ambiguous for anyone loading either beside this package — and a doc page that did would fail on use rather than on import.

See also: weakest, provenance_report.

ChemistryLab.is_evidence Method
julia
is_evidence(t) -> Bool

Whether t is a claim about the world somebody stands behind — PROV_MEASURED or PROV_PUBLISHED — as opposed to something estimated, fitted, standing in, or unstated.

A fitted value is deliberately not evidence here. It may be excellent, and whether it is depends on the data, the model and the parameter's identifiability, which is a judgement this predicate has no way to make.

ChemistryLab.provenance Method
julia
provenance(t::Traced) -> ProvenanceKind

How t was obtained.

ChemistryLab.provenance_report Method
julia
provenance_report(values) -> NamedTuple

How a collection of Traced values divides by provenance: a count per kind, the weakest present, and whether every one of them is evidence.

Meant to be printed beside a result. A table of coefficients where nine are published and one is a placeholder is a different object from one where all ten are published, and the difference does not show in any of the numbers.

values may be any iterable of Traced, including the values of a Dict.

ChemistryLab.source Method
julia
source(t::Traced) -> String

Where t came from, as free text.

ChemistryLab.uncertainty Method
julia
uncertainty(t::Traced)

What the source says about how well t is known, or nothing when it says nothing — which is not the same as saying it is exact.

ChemistryLab.value Method
julia
value(t::Traced)

The number, with its provenance dropped — explicitly, which is the only way it should be dropped.

ChemistryLab.weakest Method
julia
weakest(ts...) -> ProvenanceKind

The least sound provenance among its arguments — what a result computed from them can honestly claim.

Nothing here propagates it automatically, because Traced is not a number and does not do arithmetic. This is the function to call where a derived quantity is assembled, so that the claim is written down rather than inherited by accident.

Identifiability ​

Which parameters a measurement can actually determine, and which it only appears to. See Where the numbers come from for how this closes onto Traced: a parameter the data did not constrain comes back as a placeholder rather than as a fitted value.

ChemistryLab.Identifiability Type
julia
struct Identifiability

What a measurement can and cannot determine about a parameter vector, at one point.

Fields

  • J: the log-sensitivity matrix ∂y/∂log θ.

  • S, U, V: its singular value decomposition. V[:, k] is the parameter combination the k-th singular value belongs to.

  • condition: S[1]/S[end], or Inf when the problem is exactly rank-deficient — see below.

  • rank: how many directions the data constrain — see identifiable_rank: with a noise level, those whose standard error in log θ is below tol; without one, those before the largest gap in the spectrum that exceeds gap. It is a count of directions and not of parameters; null_participation is what says which parameters those directions leave undetermined.

  • correlation: the parameter correlation matrix from (JᵀJ)⁻¹.

  • stderr: approximate relative standard errors, σ√diag((JᵀJ)⁻¹), or nothing when there is no noise level σ.

  • rmse: the residual root-mean-square, or nothing.

  • noise: the noise level σ the rank and the errors were read with: the one given, or else the residual standard deviation of the fit; nothing when there is neither.

  • names: parameter names, for reading the output.

These are linearized errors at one point. They say which numbers in a fit deserve to be quoted; they are not confidence intervals, and a model that is nonlinear in its parameters — which is most of them — will have a likelihood that is not the ellipse this describes.

Fewer observations than parameters

The problem is then exactly rank-deficient, and that case is reported rather than smoothed over: condition is Inf, and null_participation accounts for the directions the data cannot see at all — not only the ones they see badly. correlation and stderr, on the other hand, come from a pseudo-inverse of a singular JᵀJ, so they are finite where the truth is not, and understate the error on any parameter with weight in that null space. The participation is what to read there; it is also what marks those parameters PROV_PLACEHOLDER in as_traced.

ChemistryLab.as_traced Method
julia
as_traced(id::Identifiability, θ; source, null_threshold = 0.1)
    -> Vector{Traced}

The fitted parameters, each carrying PROV_FITTED, the dataset it was fitted to, and its linearized standard error as uncertainty.

This is where an identification meets Traced: a number that came out of an optimization is not the same kind of claim as one somebody measured, and is_evidence returns false for it on purpose.

A parameter with more than null_threshold of its weight in the directions the data do not constrain is marked PROV_PLACEHOLDER instead — a value the fit had to leave somewhere is not a value the fit determined. That test is null_participation and not the parameter's position relative to the rank: the rank counts directions, the position is the caller's packing order, and confusing the two flags whichever of a trading-off pair was listed second. Pass null_threshold = Inf to report them all as fitted and take the responsibility.

ChemistryLab.identifiability Method
julia
identifiability(forward, θ; observed = nothing, noise = nothing,
                names = nothing, gap = 5.0, tol = 1.0) -> Identifiability

How much of θ the output of forward determines.

What a measurement determines depends on its noise, and the noise level σ is what turns the linearized errors into numbers. It is noise when given: the standard deviation of the instrument, in the units of y. Otherwise, with observed, it is the residual standard deviation of the fit, √(Σr² / (n − p)) over the n − p degrees of freedom that p parameters leave to n observations, and never less than √eps times the root-mean-square of the computed curve: a synthetic observation fitted exactly leaves a residual at the rounding of the arithmetic, which says nothing about a measurement and would count as determined a direction that moves the curve by less than its last digits. With neither, there is no σ: stderr comes back nothing rather than pretending, and the rank is read off the gap in the spectrum (identifiable_rank). The sensitivity is log_sensitivity's, exact; relstep is accepted and ignored, as there.

Reading it

A condition of a few means every direction is constrained. Several orders of magnitude means the data determine a combination of parameters, and V[:, end] names which one — the entries of that column say which parameters trade off against each other.

correlation is the sharper instrument when two parameters are collinear: hydration_calibration.jl found k₁ and n₁ correlated at −0.96, which says the data see a product and not its factors, and therefore that only one of the two can be fitted. Which one to keep is a modeling judgement, not a statistical one — there the rate constant was kept and the shape exponent fixed, because a rate constant is the quantity a different clinker plausibly changes.

See also: identifiable_rank, as_traced.

ChemistryLab.identifiable_rank Method
julia
identifiable_rank(S; noise = nothing, tol = 1.0, gap = 5.0) -> Int

How many directions a singular spectrum S of ∂y/∂log θ constrains.

With a noise level σ (noise, a standard deviation in the units of y), a fit determines the combination of log θ along the k-th singular direction to a standard error of σ / S[k], and the direction counts when that error is below tol: one by default, a factor e on the combination. The ratio of two quantities in the units of y has none. This is the definition the rank is meant to carry, and it takes three things: the sensitivity, the noise and a precision.

Without a noise level, the rank is read off the spectrum alone: the number of singular values before its largest ratio exceeding gap. A threshold would have units; a gap does not.

Why the noise decides when there is one

A gap separates singular values, not what the data see from what they do not: whether a direction is seen depends on the noise, which the spectrum does not carry. Two cases of this repository part the two rules:

  • scripts/hydration_calibration.jl, six rate parameters against one heat curve. The spectrum is [424, 108, 63.6, 6.46, 1.98, 0.152], and the residual of 25.9 J/g over thirty instants is a standard deviation of 29.0 J/g on the twenty-four degrees of freedom the six parameters leave. The standard errors of the six directions are then 0.07, 0.27, 0.46, 4.5, 15 and 190 in log θ: three combinations are determined, and the other three are known to within factors of 89 and more. The largest ratio, 13, is between the fifth singular value and the sixth, so the gap answers five.

  • Three separated peaks of a synthetic thermogram (test/thermogram.jl), a midpoint of 400 to 950 K and a width of 12 to 20 K each, fitted exactly. Every direction is determined. A log-sensitivity carries the size of its parameter, so the three midpoints come out 15 to 33 times above the three widths, and the largest ratio, 5.8, falls between the two groups: the gap answers three.

What the gap answers is the count that holds over the widest range of noise levels: the k-th direction is in for σ < tol S[k] and the next one out for σ > tol S[k+1], a range of S[k]/S[k+1] in σ. That is the best guess when the noise is unknown, and an assumption about the measurement, which the measurement may not meet: the heat curve above does not.

The default gap of 5 was set on the first spectrum when it was measured by central differences, [420, 100, 60, 6.3, 1.4, 0.20], whose largest ratio, 9.5, was between the third and the fourth. Exact derivatives moved the smallest singular value by a quarter, enough to move the largest ratio, and the rule with it.

Which gap, when there are several

Without a noise level, the cut is at the largest qualifying ratio, not the first one. A rate law measured against its own shrinking-core exponent has a ratio of 6 between its first two singular values and a third at the rounding of the arithmetic: only the sum of two exponents is visible. Cutting at the first ratio would say that the amplitude alone is determined, when the amplitude and the sum are, which is two. A factor of six is ordinary conditioning; a singular value at the rounding is a structure.

Returns length(S) when no ratio exceeds gap, which is the honest answer for a flat spectrum: everything is constrained, or nothing distinguishes what is not.

ChemistryLab.log_sensitivity Method
julia
log_sensitivity(forward, θ) -> Matrix

∂y/∂log θⱼ = θⱼ ∂y/∂θⱼ, with y = forward(θ), exact, by forward-mode differentiation.

Differentiating with respect to the logarithm is what makes the columns comparable: a rate constant and a dimensionless exponent have no common unit, and a matrix mixing ∂y/∂k with ∂y/∂n has a singular spectrum that says more about the units than about the data.

forward must accept ForwardDiff dual numbers and return a vector. Every forward model of this package does, a kinetic run included: an equilibrium inside it is differentiated by the implicit-function theorem at its certified answer. The cost is one evaluation of forward on dual numbers per chunk of parameters (up to twelve at a time), and the derivatives carry no step: an exact degeneracy between two parameters shows as a singular value at the rounding of the computation, not as a condition number of a few hundred.

The logarithm is undefined at zero, and so is this: rescale or shift a parameter that can vanish first.

relstep is deprecated, accepted and ignored: the derivatives are exact.

ChemistryLab.nparameters Method
julia
nparameters(id::Identifiability) -> Int

How many parameters were differentiated against.

Not length(id.S): the singular spectrum is as long as the SMALLER of the two dimensions, so with fewer observations than parameters it is shorter than the parameter vector. Every count of parameters reads this, and every count of directions reads length(id.S) — conflating the two is what made an under-determined fit report a shorter answer than it was asked about.

ChemistryLab.null_participation Method
julia
null_participation(id::Identifiability) -> Vector{Float64}

For each parameter, how much of it lies in the directions the data do not constrain: Σ_{k > rank} V[i,k]², which is between 0 and 1.

Why this and not "beyond the rank"

A rank of 2 out of 3 says the data constrain two directions in parameter space. It says nothing about which parameters are determined, because the directions are the singular vectors and the parameter order is whatever the caller packed. Reading the rank as "the first two parameters are fine" confuses an index with a subspace, and for a model where two parameters trade off it picks the wrong one — the one that happened to be listed second.

The participation is the honest version. In a model where a and c enter only as their product, it comes out [0.5, 0, 0.5]: neither a nor c is determined on its own, b is, and no ordering was consulted to say so.