Characterization, model parameters, and fitting

Heavy-end characterization

torch_flash.characterization

Heavy-end characterization independent of thermodynamic model.

CubicFractionProperties dataclass

CubicFractionProperties(
    critical_temperature: Tensor,
    critical_pressure: Tensor,
    acentric_factor: Tensor,
    m: Tensor,
)

Critical properties and acentric factors for cubic-EoS adapters.

Attributes:

Name	Type	Description
critical_temperature	Tensor	Cut critical temperatures in K.
critical_pressure	Tensor	Cut critical pressures in Pa.
acentric_factor	Tensor	Dimensionless cut acentric factors.
m	Tensor	Dimensionless alpha-function parameter associated with the requested SRK or PR correlation.

LumpedDistribution dataclass

LumpedDistribution(
    names: tuple[str, ...],
    carbon_number_bounds: tuple[tuple[int, int], ...],
    mole_fractions: Tensor,
    molar_masses: Tensor,
    densities: Tensor | None,
    properties: dict[str, Tensor],
)

Contiguous pseudo-components produced by a lumping rule.

Attributes:

Name	Type	Description
names	tuple[str, ...]	Pseudo-component labels in increasing carbon-number order.
carbon_number_bounds	tuple[tuple[int, int], ...]	Inclusive lower/upper carbon-number bounds for each lump.
mole_fractions	Tensor	Lump mole fractions.
molar_masses	Tensor	Mole-averaged lump molar masses in kg/mol.
densities	Tensor | None	Ideal-volume-preserving lump densities in kg/m3, or None.
properties	dict[str, Tensor]	Additional named, mass-weighted property tensors.

PseudoComponentCut dataclass

PseudoComponentCut(
    name: str,
    mole_fraction: float,
    normal_boiling_temperature: float,
    specific_gravity: float,
    molar_mass: float,
)

One measured or estimated heavy-end cut in SI units.

Attributes:

Name	Type	Description
name	str	Nonempty pseudo-component identifier.
mole_fraction	float	Nonnegative fraction of the complete fluid composition.
normal_boiling_temperature	float	Normal boiling temperature in K.
specific_gravity	float	Dimensionless liquid specific gravity under the source convention.
molar_mass	float	Molar mass in kg/mol.

SCNDistribution dataclass

SCNDistribution(
    carbon_numbers: Tensor,
    mole_fractions: Tensor,
    molar_masses: Tensor,
    densities: Tensor | None = None,
)

Discrete single-carbon-number representation of a plus fraction.

Attributes:

Name	Type	Description
carbon_numbers	Tensor	Strictly increasing one-dimensional integer-like carbon-number labels.
mole_fractions	Tensor	Nonnegative cut mole fractions; their sum is the represented plus fraction, not generally one.
molar_masses	Tensor	Positive cut molar masses in kg/mol.
densities	Tensor | None	Optional positive cut mass densities in kg/m3.

total_mole_fraction property

total_mole_fraction: Tensor

Return the total feed fraction represented by all SCN cuts.

Returns:

Type	Description
Tensor	Scalar sum of mole_fractions.

average_molar_mass property

average_molar_mass: Tensor

Return the mole-average molar mass.

Returns:

Type	Description
Tensor	Plus-fraction average in kg/mol.

bulk_density property

bulk_density: Tensor | None

Return ideal-volume-mixed bulk mass density.

Returns:

Type	Description
Tensor or None	Density in kg/m3, or None when cut densities are unavailable.

to

to(
    *,
    dtype: dtype | None = None,
    device: device | str | None = None,
) -> SCNDistribution

Return the distribution on a common dtype and device.

Parameters:

Name	Type	Description	Default
dtype	dtype | None	Target floating dtype for fractions, masses, and densities.	None
device	device | str | None	Target device for every tensor.	None

Returns:

Type	Description
SCNDistribution	New immutable distribution. Carbon numbers retain integer dtype.

Source code in src/torch_flash/characterization/types.py

def to(
    self,
    *,
    dtype: torch.dtype | None = None,
    device: torch.device | str | None = None,
) -> SCNDistribution:
    """Return the distribution on a common dtype and device.

    Parameters
    ----------
    dtype:
        Target floating dtype for fractions, masses, and densities.
    device:
        Target device for every tensor.

    Returns
    -------
    SCNDistribution
        New immutable distribution. Carbon numbers retain integer dtype.
    """
    return SCNDistribution(
        self.carbon_numbers.to(device=device),
        self.mole_fractions.to(dtype=dtype, device=device),
        self.molar_masses.to(dtype=dtype, device=device),
        None if self.densities is None else self.densities.to(dtype=dtype, device=device),
    )

pedersen_cubic_properties

pedersen_cubic_properties(
    distribution: SCNDistribution,
    eos: CubicCharacterization,
    parameter_set: ParameterSource = "characterization.pedersen-2024",
) -> CubicFractionProperties

Map characterized SCN cuts to SRK or PR properties.

The default implements Pedersen et al. (2024), Eqs. 5.1-5.5 and Table 5.3. The C200 heavy-aromatic parameter set implements Pedersen, Milter, and Sorensen, SPE Journal 9 (2004), Eqs. 7-14 and Table 4, doi:10.2118/88364-PA. Its alpha parameter follows the published inverse-molar-mass continuation above 1000 g/mol. The correlation coefficients are EoS-specific; the input distribution and density split remain model-neutral.

Parameters:

Name	Type	Description	Default
distribution	SCNDistribution	SCN distribution with molar masses in kg/mol and densities in kg/m3.	required
eos	CubicCharacterization	"SRK" or "PR"; each selects a distinct coefficient table.	required
parameter_set	ParameterSource	Characterization parameter identifier, YAML path, or explicit record.	'characterization.pedersen-2024'

Returns:

Type	Description
CubicFractionProperties	Cut properties with the same shape, dtype, and device as the distribution.

Raises:

Type	Description
ValueError	If the EoS name is invalid, densities are absent, or the fitted acentric-factor inversion has no real solution.
ParameterDatabaseError	If the selected document is not a compatible characterization set.

Source code in src/torch_flash/characterization/cubic.py

def pedersen_cubic_properties(
    distribution: SCNDistribution,
    eos: CubicCharacterization,
    parameter_set: ParameterSource = "characterization.pedersen-2024",
) -> CubicFractionProperties:
    """Map characterized SCN cuts to SRK or PR properties.

    The default implements Pedersen et al. (2024), Eqs. 5.1-5.5 and
    Table 5.3. The C200 heavy-aromatic parameter set implements Pedersen,
    Milter, and Sorensen, *SPE Journal* 9 (2004), Eqs. 7-14 and Table 4,
    doi:10.2118/88364-PA. Its alpha parameter follows the published
    inverse-molar-mass continuation above 1000 g/mol. The correlation
    coefficients are EoS-specific; the input distribution and density split
    remain model-neutral.

    Parameters
    ----------
    distribution:
        SCN distribution with molar masses in kg/mol and densities in kg/m3.
    eos:
        ``"SRK"`` or ``"PR"``; each selects a distinct coefficient table.
    parameter_set:
        Characterization parameter identifier, YAML path, or explicit record.

    Returns
    -------
    CubicFractionProperties
        Cut properties with the same shape, dtype, and device as the
        distribution.

    Raises
    ------
    ValueError
        If the EoS name is invalid, densities are absent, or the fitted
        acentric-factor inversion has no real solution.
    ParameterDatabaseError
        If the selected document is not a compatible characterization set.
    """
    if eos not in ("SRK", "PR"):
        raise ValueError("eos must be 'SRK' or 'PR'")
    if distribution.densities is None:
        raise ValueError("Pedersen cubic properties require characterized SCN densities")
    loaded = load_model_parameters(parameter_set)
    if loaded.model_kind != "characterization":
        raise ParameterDatabaseError(f"{loaded.identifier!r} is not a characterization set")
    tables = loaded.parameters.get("cubic_properties")
    if not isinstance(tables, Mapping):
        raise ParameterDatabaseError(f"{loaded.identifier!r} requires cubic_properties")
    record = tables.get(eos)
    if not isinstance(record, Mapping):
        raise ParameterDatabaseError(f"{loaded.identifier!r} has no {eos} property table")
    c1, c2, c3, c4 = _coefficients(record, "critical_temperature", 4, loaded.identifier)
    d1, d2, d3, d4, d5 = _coefficients(record, "log_critical_pressure", 5, loaded.identifier)
    e1, e2, e3, e4 = _coefficients(record, "m", 4, loaded.identifier)
    w0, w1, w2 = _coefficients(record, "m_to_acentric_factor", 3, loaded.identifier)
    molar_mass = 1.0e3 * distribution.molar_masses
    density = 1.0e-3 * distribution.densities
    critical_temperature = c1 * density + c2 * torch.log(molar_mass) + c3 * molar_mass
    critical_temperature = critical_temperature + c4 / molar_mass
    log_pressure_atm = d1 + d2 * density**d5 + d3 / molar_mass + d4 / molar_mass.square()
    critical_pressure = 101_325.0 * torch.exp(log_pressure_atm)
    m = e1 + e2 * molar_mass + e3 * density + e4 * molar_mass.square()
    threshold_value = record.get("inverse_mass_m_threshold")
    if threshold_value is not None:
        if (
            not isinstance(threshold_value, int | float)
            or threshold_value <= 0.0
            or not torch.isfinite(molar_mass.new_tensor(float(threshold_value)))
        ):
            raise ParameterDatabaseError(
                f"{loaded.identifier!r} inverse_mass_m_threshold must be positive and finite"
            )
        threshold = molar_mass.new_tensor(float(threshold_value))
        above_threshold = molar_mass > threshold
        if bool(above_threshold.any()):
            if distribution.carbon_numbers.numel() < 2:
                raise ValueError(
                    "inverse-mass alpha continuation requires at least two logarithmic-density cuts"
                )
            split = loaded.parameters.get("plus_split")
            relation = split.get("molecular_weight") if isinstance(split, Mapping) else None
            if not isinstance(relation, Mapping):
                raise ParameterDatabaseError(
                    f"{loaded.identifier!r} requires plus_split.molecular_weight"
                )
            slope_value = relation.get("slope")
            intercept_value = relation.get("intercept")
            if not isinstance(slope_value, int | float) or not isinstance(
                intercept_value, int | float
            ):
                raise ParameterDatabaseError(
                    f"{loaded.identifier!r} molecular-weight relation must be numeric"
                )
            threshold_carbon_number = (threshold - float(intercept_value)) / float(slope_value)
            carbon_number = distribution.carbon_numbers.to(density)
            log_span = torch.log(carbon_number[-1]) - torch.log(carbon_number[0])
            if not bool(torch.isfinite(log_span) & (torch.abs(log_span) > 0.0)):
                raise ValueError(
                    "inverse-mass alpha continuation requires distinct carbon-number cuts"
                )
            density_slope = (density[-1] - density[0]) / log_span
            threshold_density = density[0] + density_slope * (
                torch.log(threshold_carbon_number) - torch.log(carbon_number[0])
            )
            threshold_m = e1 + e2 * threshold + e3 * threshold_density + e4 * threshold.square()
            m = torch.where(
                above_threshold,
                threshold * threshold_m / molar_mass,
                m,
            )
    discriminant = w1 * w1 - 4.0 * w2 * (w0 - m)
    if bool((discriminant < 0.0).any()):
        raise ValueError("Pedersen cubic correlation produced no real acentric factor")
    root_first = (-w1 + torch.sqrt(discriminant)) / (2.0 * w2)
    root_second = (-w1 - torch.sqrt(discriminant)) / (2.0 * w2)
    acentric = torch.where(
        torch.abs(root_first) <= torch.abs(root_second),
        root_first,
        root_second,
    )
    return CubicFractionProperties(
        critical_temperature,
        critical_pressure,
        acentric,
        m,
    )

pedersen_density_split

pedersen_density_split(
    distribution: SCNDistribution,
    plus_density: float | Tensor,
    *,
    anchor_density: float | Tensor | None = None,
    anchor_carbon_number: int | None = None,
    parameter_set: ParameterSource = "characterization.pedersen-2024",
) -> SCNDistribution

Assign rho_N = C + D ln(CN) while matching bulk plus density.

Density inputs and results are SI (kg/m3). If no measured anchor is provided, Pedersen's suggested C6 density ratio is applied to the carbon number immediately preceding the plus fraction.

Parameters:

Name	Type	Description	Default
distribution	SCNDistribution	SCN distribution whose moles and molar masses are retained.	required
plus_density	float | Tensor	Positive scalar target bulk density in kg/m3.	required
anchor_density	float | Tensor | None	Optional positive measured density in kg/m3 at anchor_carbon_number.	None
anchor_carbon_number	int | None	Positive anchor label. By default, the number immediately before the first SCN is used.	None
parameter_set	ParameterSource	Pedersen characterization parameter source.	'characterization.pedersen-2024'

Returns:

Type	Description
SCNDistribution	New distribution with per-cut densities matching the bulk ideal-volume balance.

Raises:

Type	Description
ValueError	If target/anchor values or the anchor label are invalid.
ConvergenceError	If a positive density trend cannot satisfy the volume balance.
ParameterDatabaseError	If the selected document does not define the density rule.

Source code in src/torch_flash/characterization/distributions.py

def pedersen_density_split(
    distribution: SCNDistribution,
    plus_density: float | Tensor,
    *,
    anchor_density: float | Tensor | None = None,
    anchor_carbon_number: int | None = None,
    parameter_set: ParameterSource = "characterization.pedersen-2024",
) -> SCNDistribution:
    """Assign ``rho_N = C + D ln(CN)`` while matching bulk plus density.

    Density inputs and results are SI (kg/m3). If no measured anchor is
    provided, Pedersen's suggested C6 density ratio is applied to the carbon
    number immediately preceding the plus fraction.

    Parameters
    ----------
    distribution:
        SCN distribution whose moles and molar masses are retained.
    plus_density:
        Positive scalar target bulk density in kg/m3.
    anchor_density:
        Optional positive measured density in kg/m3 at
        ``anchor_carbon_number``.
    anchor_carbon_number:
        Positive anchor label. By default, the number immediately before the
        first SCN is used.
    parameter_set:
        Pedersen characterization parameter source.

    Returns
    -------
    SCNDistribution
        New distribution with per-cut densities matching the bulk
        ideal-volume balance.

    Raises
    ------
    ValueError
        If target/anchor values or the anchor label are invalid.
    ConvergenceError
        If a positive density trend cannot satisfy the volume balance.
    ParameterDatabaseError
        If the selected document does not define the density rule.
    """
    parameters, source = _characterization_parameters(parameter_set)
    split = parameters.get("plus_split")
    if not isinstance(split, Mapping):
        raise ParameterDatabaseError(f"{source} requires a plus_split mapping")
    if split.get("density_log_carbon_number") is not True:
        raise ParameterDatabaseError(f"{source} does not define the logarithmic density rule")
    target = torch.as_tensor(
        plus_density,
        dtype=distribution.mole_fractions.dtype,
        device=distribution.mole_fractions.device,
    )
    if target.ndim or not bool(torch.isfinite(target) & (target > 0.0)):
        raise ValueError("plus density must be a finite positive scalar")
    if anchor_density is None:
        ratio = _positive_number(split, "default_anchor_density_ratio", source)
        anchor = ratio * target
    else:
        anchor = torch.as_tensor(
            anchor_density,
            dtype=target.dtype,
            device=target.device,
        )
    if anchor.ndim or not bool(torch.isfinite(anchor) & (anchor > 0.0)):
        raise ValueError("anchor density must be a finite positive scalar")
    if anchor_carbon_number is None:
        anchor_carbon_number = int(distribution.carbon_numbers[0]) - 1
    if anchor_carbon_number < 1:
        raise ValueError("anchor carbon number must be positive")

    log_delta = torch.log(distribution.carbon_numbers.to(target)) - torch.log(
        target.new_tensor(float(anchor_carbon_number))
    )
    mass = distribution.mole_fractions * distribution.molar_masses

    def residual(density_slope: Tensor) -> Tensor:
        densities = anchor + density_slope * log_delta
        bulk = mass.sum() / torch.sum(mass / densities)
        return bulk - target

    density_slope = (target - anchor) / torch.clamp(log_delta.mean(), min=1.0e-6)
    lower = -0.95 * anchor / torch.clamp(log_delta.max(), min=1.0e-6)
    upper = 10.0 * target
    for _ in range(30):
        value = residual(density_slope)
        derivative = torch.func.grad(residual)(density_slope)
        density_slope = torch.clamp(density_slope - value / derivative, lower, upper)
    densities = anchor + density_slope * log_delta
    if (
        bool((densities <= 0.0).any())
        or float(torch.abs(residual(density_slope) / target).detach()) > 1.0e-10
    ):
        raise ConvergenceError("Pedersen density split did not satisfy its volume balance")
    return SCNDistribution(
        distribution.carbon_numbers,
        distribution.mole_fractions,
        distribution.molar_masses,
        densities,
    )

pedersen_logarithmic_split

pedersen_logarithmic_split(
    plus_mole_fraction: float | Tensor,
    plus_molar_mass: float | Tensor,
    *,
    first_carbon_number: int = 7,
    max_carbon_number: int | None = None,
    parameter_set: ParameterSource = "characterization.pedersen-2024",
    dtype: dtype | None = None,
    device: device | str | None = None,
) -> SCNDistribution

Split a plus fraction with Pedersen's logarithmic molar distribution.

The finite distribution satisfies both plus-fraction mole and molar-mass balances (Pedersen et al., 2024, Eqs. 5.10-5.12). Molecular weights use Eq. 5.22. Measured extended compositions should be preferred whenever available, as emphasized by the source.

Parameters:

Name	Type	Description	Default
plus_mole_fraction	float | Tensor	Positive scalar mole fraction represented by the plus cut.	required
plus_molar_mass	float | Tensor	Positive scalar average molar mass in kg/mol.	required
first_carbon_number	int	Inclusive SCN label range. The default maximum comes from the parameter set.	7
max_carbon_number	int	Inclusive SCN label range. The default maximum comes from the parameter set.	7
parameter_set	ParameterSource	Pedersen characterization parameter source.	'characterization.pedersen-2024'
dtype	dtype | None	Output tensor options; omitted values use runtime configuration.	None
device	dtype | None	Output tensor options; omitted values use runtime configuration.	None

Returns:

Type	Description
SCNDistribution	Finite distribution preserving total plus moles and average molar mass.

Raises:

Type	Description
ValueError	If scalar inputs or carbon bounds are invalid, or the requested mean mass lies outside the finite SCN range.
ConvergenceError	If the moment-constrained solve does not close the mass balance.
ParameterDatabaseError	If the selected parameter document lacks required coefficients.

Source code in src/torch_flash/characterization/distributions.py

def pedersen_logarithmic_split(
    plus_mole_fraction: float | Tensor,
    plus_molar_mass: float | Tensor,
    *,
    first_carbon_number: int = 7,
    max_carbon_number: int | None = None,
    parameter_set: ParameterSource = "characterization.pedersen-2024",
    dtype: torch.dtype | None = None,
    device: torch.device | str | None = None,
) -> SCNDistribution:
    """Split a plus fraction with Pedersen's logarithmic molar distribution.

    The finite distribution satisfies both plus-fraction mole and molar-mass
    balances (Pedersen et al., 2024, Eqs. 5.10-5.12). Molecular weights use
    Eq. 5.22. Measured extended compositions should be preferred whenever
    available, as emphasized by the source.

    Parameters
    ----------
    plus_mole_fraction:
        Positive scalar mole fraction represented by the plus cut.
    plus_molar_mass:
        Positive scalar average molar mass in kg/mol.
    first_carbon_number, max_carbon_number:
        Inclusive SCN label range. The default maximum comes from the
        parameter set.
    parameter_set:
        Pedersen characterization parameter source.
    dtype, device:
        Output tensor options; omitted values use runtime configuration.

    Returns
    -------
    SCNDistribution
        Finite distribution preserving total plus moles and average molar mass.

    Raises
    ------
    ValueError
        If scalar inputs or carbon bounds are invalid, or the requested mean
        mass lies outside the finite SCN range.
    ConvergenceError
        If the moment-constrained solve does not close the mass balance.
    ParameterDatabaseError
        If the selected parameter document lacks required coefficients.
    """
    dtype, device = resolve_tensor_options(dtype, device)
    parameters, source = _characterization_parameters(parameter_set)
    split = parameters.get("plus_split")
    if not isinstance(split, Mapping):
        raise ParameterDatabaseError(f"{source} requires a plus_split mapping")
    if max_carbon_number is None:
        raw_max = split.get("default_max_carbon_number")
        if not isinstance(raw_max, int):
            raise ParameterDatabaseError(f"{source} default_max_carbon_number must be an integer")
        max_carbon_number = raw_max
    if (
        not isinstance(first_carbon_number, int)
        or not isinstance(max_carbon_number, int)
        or first_carbon_number < 1
        or max_carbon_number < first_carbon_number
    ):
        raise ValueError("carbon-number bounds must be ordered positive integers")
    relation = split.get("molecular_weight")
    if not isinstance(relation, Mapping):
        raise ParameterDatabaseError(f"{source} requires a molecular_weight mapping")
    slope = _positive_number(relation, "slope", source)
    intercept_value = relation.get("intercept")
    if not isinstance(intercept_value, int | float) or not isfinite(intercept_value):
        raise ParameterDatabaseError(f"{source} molecular-weight intercept must be finite")
    intercept = float(intercept_value)

    fraction = torch.as_tensor(plus_mole_fraction, dtype=dtype, device=device)
    target_mass = torch.as_tensor(plus_molar_mass, dtype=dtype, device=device)
    if fraction.ndim or target_mass.ndim:
        raise ValueError("plus mole fraction and molar mass must be scalar")
    if not bool(torch.isfinite(fraction) & (fraction > 0.0)):
        raise ValueError("plus mole fraction must be finite and positive")
    if not bool(torch.isfinite(target_mass) & (target_mass > 0.0)):
        raise ValueError("plus molar mass must be finite and positive")
    carbon_numbers = torch.arange(
        first_carbon_number,
        max_carbon_number + 1,
        dtype=dtype,
        device=device,
    )
    molar_masses = (slope * carbon_numbers + intercept) * 1.0e-3
    if not bool(
        (target_mass >= molar_masses[0] - 1.0e-12) & (target_mass <= molar_masses[-1] + 1.0e-12)
    ):
        raise ValueError("plus molar mass lies outside the selected finite SCN range")

    # Fixed Newton iterations retain an autodiff path through the moment
    # constraint. d<E[M]>/d lambda is cov(M, carbon number).
    log_slope = target_mass.new_zeros(())
    centered_carbon = carbon_numbers - carbon_numbers[0]
    for _ in range(40):
        weights = torch.softmax(log_slope * centered_carbon, dim=0)
        mean_mass = torch.sum(weights * molar_masses)
        mean_carbon = torch.sum(weights * centered_carbon)
        covariance = torch.sum(
            weights * (molar_masses - mean_mass) * (centered_carbon - mean_carbon)
        )
        step = torch.clamp((mean_mass - target_mass) / covariance, -2.0, 2.0)
        log_slope = torch.clamp(log_slope - step, -50.0, 50.0)
    weights = torch.softmax(log_slope * centered_carbon, dim=0)
    mole_fractions = fraction * weights
    mean_mass = torch.sum(weights * molar_masses)
    mass_residual = torch.abs(mean_mass - target_mass)
    mass_scale = torch.maximum(torch.abs(mean_mass), torch.abs(target_mass))
    # The final moment is accumulated in the requested dtype. Permit several
    # ulps of reduction-order variation while retaining the historical
    # float64 absolute gate.
    precision_limit = 8.0 * torch.finfo(target_mass.dtype).eps * float(mass_scale.detach())
    if float(mass_residual.detach()) > max(1.0e-9, precision_limit):
        raise ConvergenceError("Pedersen logarithmic split did not satisfy its mass balance")
    return SCNDistribution(
        carbon_numbers.to(torch.int64),
        mole_fractions,
        molar_masses,
    )

whitson_gamma_split

whitson_gamma_split(
    plus_mole_fraction: float | Tensor,
    plus_molar_mass: float | Tensor,
    *,
    first_carbon_number: int = 7,
    max_carbon_number: int = 80,
    shape: float | Tensor | None = None,
    minimum_molar_mass: float | Tensor | None = None,
    parameter_set: ParameterSource = "characterization.whitson-2000",
    dtype: dtype | None = None,
    device: device | str | None = None,
) -> SCNDistribution

Discretize Whitson's shifted gamma distribution into SCN-like bins.

The last requested bin contains the complete tail to infinite molecular weight, so total moles and average molar mass are preserved exactly apart from floating-point roundoff. shape=1 gives the exponential special case discussed by both Whitson and Pedersen.

Parameters:

Name	Type	Description	Default
plus_mole_fraction	float | Tensor	Positive scalar mole fraction represented by the plus cut.	required
plus_molar_mass	float | Tensor	Positive scalar average molar mass in kg/mol.	required
first_carbon_number	int	Inclusive nominal bin labels. The final bin contains the complete remaining gamma tail.	7
max_carbon_number	int	Inclusive nominal bin labels. The final bin contains the complete remaining gamma tail.	7
shape	float | Tensor | None	Positive gamma shape parameter. None uses the source default.	None
minimum_molar_mass	float | Tensor | None	Positive shifted-gamma lower bound in kg/mol, below the average.	None
parameter_set	ParameterSource	Whitson characterization parameter source.	'characterization.whitson-2000'
dtype	dtype | None	Output tensor options; omitted values use runtime configuration.	None
device	dtype | None	Output tensor options; omitted values use runtime configuration.	None

Returns:

Type	Description
SCNDistribution	Discrete distribution preserving total moles and mean molar mass.

Raises:

Type	Description
ValueError	If scalar values, bounds, shape, or minimum mass are invalid.
ParameterDatabaseError	If the selected document lacks required gamma parameters.

Source code in src/torch_flash/characterization/distributions.py

def whitson_gamma_split(
    plus_mole_fraction: float | Tensor,
    plus_molar_mass: float | Tensor,
    *,
    first_carbon_number: int = 7,
    max_carbon_number: int = 80,
    shape: float | Tensor | None = None,
    minimum_molar_mass: float | Tensor | None = None,
    parameter_set: ParameterSource = "characterization.whitson-2000",
    dtype: torch.dtype | None = None,
    device: torch.device | str | None = None,
) -> SCNDistribution:
    """Discretize Whitson's shifted gamma distribution into SCN-like bins.

    The last requested bin contains the complete tail to infinite molecular
    weight, so total moles and average molar mass are preserved exactly apart
    from floating-point roundoff. ``shape=1`` gives the exponential special
    case discussed by both Whitson and Pedersen.

    Parameters
    ----------
    plus_mole_fraction:
        Positive scalar mole fraction represented by the plus cut.
    plus_molar_mass:
        Positive scalar average molar mass in kg/mol.
    first_carbon_number, max_carbon_number:
        Inclusive nominal bin labels. The final bin contains the complete
        remaining gamma tail.
    shape:
        Positive gamma shape parameter. ``None`` uses the source default.
    minimum_molar_mass:
        Positive shifted-gamma lower bound in kg/mol, below the average.
    parameter_set:
        Whitson characterization parameter source.
    dtype, device:
        Output tensor options; omitted values use runtime configuration.

    Returns
    -------
    SCNDistribution
        Discrete distribution preserving total moles and mean molar mass.

    Raises
    ------
    ValueError
        If scalar values, bounds, shape, or minimum mass are invalid.
    ParameterDatabaseError
        If the selected document lacks required gamma parameters.
    """
    dtype, device = resolve_tensor_options(dtype, device)
    parameters, source = _characterization_parameters(parameter_set)
    gamma = parameters.get("gamma_distribution")
    if not isinstance(gamma, Mapping):
        raise ParameterDatabaseError(f"{source} requires a gamma_distribution mapping")
    if shape is None:
        shape = _positive_number(gamma, "default_shape", source)
    alpha = torch.as_tensor(shape, dtype=dtype, device=device)
    fraction = torch.as_tensor(plus_mole_fraction, dtype=dtype, device=device)
    average_mass = torch.as_tensor(plus_molar_mass, dtype=dtype, device=device)
    if any(value.ndim for value in (alpha, fraction, average_mass)):
        raise ValueError("gamma shape and plus-fraction inputs must be scalar")
    if not bool(
        torch.isfinite(alpha)
        & (alpha > 0.0)
        & torch.isfinite(fraction)
        & (fraction > 0.0)
        & torch.isfinite(average_mass)
        & (average_mass > 0.0)
    ):
        raise ValueError("gamma shape, mole fraction, and molar mass must be finite and positive")
    if first_carbon_number < 1 or max_carbon_number < first_carbon_number:
        raise ValueError("carbon-number bounds must be ordered positive integers")
    if minimum_molar_mass is None:
        relation = gamma.get("recommended_minimum_molecular_weight_relation")
        if not isinstance(relation, Mapping):
            raise ParameterDatabaseError(
                f"{source} requires recommended minimum-molecular-weight coefficients"
            )
        scale = _positive_number(relation, "scale", source)
        multiplier = _positive_number(relation, "multiplier", source)
        exponent = _positive_number(relation, "exponent", source)
        eta = 1.0e-3 * scale * (1.0 - 1.0 / (1.0 + multiplier / alpha**exponent))
    else:
        eta = torch.as_tensor(minimum_molar_mass, dtype=dtype, device=device)
    if eta.ndim or not bool(torch.isfinite(eta) & (eta > 0.0) & (eta < average_mass)):
        raise ValueError("minimum molar mass must be finite, positive, and below the average")
    beta = (average_mass - eta) / alpha

    carbon_numbers = torch.arange(
        first_carbon_number,
        max_carbon_number + 1,
        dtype=dtype,
        device=device,
    )
    boundary_increment = _positive_number(
        gamma,
        "molecular_weight_boundary_increment",
        source,
    )
    # Whitson and Brule (2000), Table 5.4 and program GAMSPL, place
    # molecular-weight boundaries at eta, eta + 14, eta + 28, ... g/mol.
    # The carbon numbers are therefore nominal labels for consecutive bins;
    # they do not define midpoints through the n-paraffin relation.
    offsets = torch.arange(
        carbon_numbers.numel(),
        dtype=dtype,
        device=device,
    )
    lower = eta + boundary_increment * offsets
    upper = torch.cat(
        (
            lower[1:],
            torch.full_like(eta.reshape(1), torch.inf),
        )
    )
    y_lower = torch.clamp((lower - eta) / beta, min=0.0)
    y_upper = torch.clamp((upper - eta) / beta, min=0.0)
    p_lower = torch.special.gammainc(alpha, y_lower)
    p_upper = torch.special.gammainc(alpha, y_upper)
    probabilities = p_upper - p_lower
    first_moment_probability = torch.special.gammainc(alpha + 1.0, y_upper)
    first_moment_probability = first_moment_probability - torch.special.gammainc(
        alpha + 1.0, y_lower
    )
    bin_masses = eta + alpha * beta * first_moment_probability / probabilities
    return SCNDistribution(
        carbon_numbers.to(torch.int64),
        fraction * probabilities,
        bin_masses,
    )

equal_weight_lump

equal_weight_lump(
    distribution: SCNDistribution,
    groups: int,
    *,
    properties: Mapping[str, Tensor] | None = None,
) -> LumpedDistribution

Lump contiguous SCN cuts into approximately equal-weight groups.

Critical properties supplied through properties use Pedersen's weight-average rule (2024, Eqs. 5.29-5.31). Molar mass is mole averaged; density, when present, preserves ideal mixed volume.

Parameters:

Name	Type	Description	Default
distribution	SCNDistribution	Ordered SCN distribution to combine.	required
groups	int	Number of nonempty contiguous lumps, from one through the SCN count.	required
properties	Mapping[str, Tensor] | None	Optional named per-SCN tensors to aggregate by mass weighting.	None

Returns:

Type	Description
LumpedDistribution	Contiguous pseudo-components preserving total moles and mass.

Raises:

Type	Description
ValueError	If groups is invalid or a supplied property name/shape is invalid.

Source code in src/torch_flash/characterization/lumping.py

def equal_weight_lump(
    distribution: SCNDistribution,
    groups: int,
    *,
    properties: Mapping[str, Tensor] | None = None,
) -> LumpedDistribution:
    """Lump contiguous SCN cuts into approximately equal-weight groups.

    Critical properties supplied through ``properties`` use Pedersen's
    weight-average rule (2024, Eqs. 5.29-5.31). Molar mass is mole averaged;
    density, when present, preserves ideal mixed volume.

    Parameters
    ----------
    distribution:
        Ordered SCN distribution to combine.
    groups:
        Number of nonempty contiguous lumps, from one through the SCN count.
    properties:
        Optional named per-SCN tensors to aggregate by mass weighting.

    Returns
    -------
    LumpedDistribution
        Contiguous pseudo-components preserving total moles and mass.

    Raises
    ------
    ValueError
        If ``groups`` is invalid or a supplied property name/shape is invalid.
    """
    count = distribution.carbon_numbers.numel()
    if not isinstance(groups, int) or groups < 1 or groups > count:
        raise ValueError("groups must be an integer from one through the SCN count")
    supplied = {} if properties is None else dict(properties)
    for name, values in supplied.items():
        if not isinstance(name, str) or not name:
            raise ValueError("lumped property names must be non-empty strings")
        if values.shape != distribution.mole_fractions.shape:
            raise ValueError(f"lumped property {name!r} must match the SCN distribution")

    mass = distribution.mole_fractions * distribution.molar_masses
    cumulative = torch.cumsum(mass, dim=0)
    targets = (
        mass.sum()
        * torch.arange(
            1,
            groups,
            dtype=mass.dtype,
            device=mass.device,
        )
        / groups
    )
    boundaries = [0]
    for target in targets:
        candidate = int(torch.searchsorted(cumulative, target).detach()) + 1
        minimum = boundaries[-1] + 1
        maximum = count - (groups - len(boundaries))
        boundaries.append(min(max(candidate, minimum), maximum))
    boundaries.append(count)

    names: list[str] = []
    bounds: list[tuple[int, int]] = []
    fractions: list[Tensor] = []
    molar_masses: list[Tensor] = []
    densities: list[Tensor] = []
    lumped_properties: dict[str, list[Tensor]] = {name: [] for name in supplied}
    for start, stop in pairwise(boundaries):
        current_fraction = distribution.mole_fractions[start:stop]
        current_mass = mass[start:stop]
        total_fraction = current_fraction.sum()
        total_mass = current_mass.sum()
        lower = int(distribution.carbon_numbers[start])
        upper = int(distribution.carbon_numbers[stop - 1])
        names.append(f"C{lower}" if lower == upper else f"C{lower}-C{upper}")
        bounds.append((lower, upper))
        fractions.append(total_fraction)
        molar_masses.append(total_mass / total_fraction)
        if distribution.densities is not None:
            densities.append(
                total_mass / torch.sum(current_mass / distribution.densities[start:stop])
            )
        for name, values in supplied.items():
            lumped_properties[name].append(
                torch.sum(current_mass * values[start:stop]) / total_mass
            )
    return LumpedDistribution(
        tuple(names),
        tuple(bounds),
        torch.stack(fractions),
        torch.stack(molar_masses),
        None if distribution.densities is None else torch.stack(densities),
        {name: torch.stack(values) for name, values in lumped_properties.items()},
    )

Interaction and group-contribution parameters

torch_flash.parameters

Published thermodynamic parameter sets.

ModelParameterSet dataclass

ModelParameterSet(
    identifier: str,
    model_kind: str,
    model: str,
    version: str,
    parameters: ParameterDocument,
    units: ParameterDocument = dict(),
    references: tuple[ParameterDocument, ...] = (),
    description: str = "",
    source: str = "<api>",
)

Validated model parameter set independent of tensor dtype and device.

Parameters can be supplied directly through this dataclass, read from a custom YAML path, or selected from the bundled database by identifier. The stored mappings are recursively read-only. Use :meth:as_dict when a mutable deep copy is needed.

Attributes:

Name	Type	Description
identifier	str	Stable lowercase parameter-set identifier.
model_kind	str	Scientific domain such as "cubic", "activity", or "multiparameter".
model	str	Model/form name declared by the document.
version	str	Parameter-set version or revision.
parameters	ParameterDocument	Recursively read-only numerical/scientific payload.
units	ParameterDocument	Recursively read-only units mapping.
references	tuple[ParameterDocument, ...]	Source records for equations, coefficients, or fits.
description	str	Human-readable parameter-set scope.
source	str	Bundled resource label, custom path, or "<api>".

from_document classmethod

from_document(
    document: Mapping[str, Any], *, source: str = "<api>"
) -> ModelParameterSet

Validate a model-parameter YAML document represented as a mapping.

Parameters:

Name	Type	Description	Default
document	Mapping[str, Any]	Parsed YAML mapping following the torch-flash model-parameter schema.	required
source	str	Label used in diagnostics and stored provenance.	'<api>'

Returns:

Type	Description
ModelParameterSet	Immutable validated parameter document.

Raises:

Type	Description
ParameterDatabaseError	If required headers, identifiers, mappings, units, or references are invalid.

Source code in src/torch_flash/database.py

@classmethod
def from_document(
    cls,
    document: Mapping[str, Any],
    *,
    source: str = "<api>",
) -> ModelParameterSet:
    """Validate a model-parameter YAML document represented as a mapping.

    Parameters
    ----------
    document:
        Parsed YAML mapping following the torch-flash model-parameter
        schema.
    source:
        Label used in diagnostics and stored provenance.

    Returns
    -------
    ModelParameterSet
        Immutable validated parameter document.

    Raises
    ------
    ParameterDatabaseError
        If required headers, identifiers, mappings, units, or references
        are invalid.
    """
    _validate_header(document, _MODEL_FORMAT, source)
    parameters = document.get("parameters")
    if not isinstance(parameters, Mapping):
        raise ParameterDatabaseError(f"{source} requires a 'parameters' mapping")
    units = document.get("units", {})
    references = document.get("references", [])
    if not isinstance(units, Mapping):
        raise ParameterDatabaseError(f"{source} requires 'units' to be a mapping")
    if not isinstance(references, list):
        raise ParameterDatabaseError(f"{source} requires 'references' to be a list")
    description = document.get("description", "")
    if not isinstance(description, str):
        raise ParameterDatabaseError(f"{source} requires 'description' to be a string")
    return cls(
        identifier=_require_string(document, "id", source),
        model_kind=_require_string(document, "model_kind", source),
        model=_require_string(document, "model", source),
        version=_require_string(document, "version", source),
        parameters=parameters,
        units=units,
        references=tuple(cast(list[ParameterDocument], references)),
        description=description,
        source=source,
    )

as_dict

as_dict() -> dict[str, Any]

Return a mutable deep copy of the numerical parameter payload.

Returns:

Type	Description
dict	Deeply mutable copy. Mutating it cannot alter the cached parameter set or any existing model tensors.

Source code in src/torch_flash/database.py

def as_dict(self) -> dict[str, Any]:
    """Return a mutable deep copy of the numerical parameter payload.

    Returns
    -------
    dict
        Deeply mutable copy. Mutating it cannot alter the cached parameter
        set or any existing model tensors.
    """
    return cast(dict[str, Any], _thaw(self.parameters))

CubicInteractionParameters dataclass

CubicInteractionParameters(
    kij: Tensor, lij: Tensor, parameter_set: str
)

Validated cubic attraction and covolume interaction matrices.

Attributes:

Name	Type	Description
kij	Tensor	Symmetric dimensionless attraction-interaction matrix.
lij	Tensor	Symmetric dimensionless cross-covolume-interaction matrix.
parameter_set	str	Identifier of the parameter set from which the matrices were selected.

Notes

Both matrices have zero diagonals and axes ordered exactly as the :class:~torch_flash.components.ComponentSet supplied to the loader.

PPR78GroupContributionParameters dataclass

PPR78GroupContributionParameters(
    group_names: tuple[str, ...],
    group_fractions: Tensor,
    group_a: Tensor,
    group_b: Tensor,
    reference_temperature: float,
    parameter_set: str,
)

Selected PPR78 or E-PPR78 decompositions and interactions.

Attributes:

Name	Type	Description
group_names	tuple[str, ...]	Ordered group identifiers.
group_fractions	Tensor	Normalized component group fractions with shape (ncomponents, ngroups).
group_a, group_b		Symmetric universal interaction matrices in Pa with shape (ngroups, ngroups).
reference_temperature	float	Positive interaction reference temperature in K.
parameter_set	str	Versioned parameter-set identifier.

available_parameter_sets

available_parameter_sets(
    *, model_kind: str | None = None
) -> tuple[str, ...]

Return bundled parameter-set identifiers, optionally filtered by kind.

Parameters:

Name	Type	Description	Default
model_kind	str | None	Exact model-kind filter, for example "cubic" or "activity". None returns every bundled identifier.	None

Returns:

Type	Description
tuple[str, ...]	Sorted immutable tuple of matching identifiers. An unknown kind yields an empty tuple.

Source code in src/torch_flash/database.py

def available_parameter_sets(*, model_kind: str | None = None) -> tuple[str, ...]:
    """Return bundled parameter-set identifiers, optionally filtered by kind.

    Parameters
    ----------
    model_kind:
        Exact model-kind filter, for example ``"cubic"`` or ``"activity"``.
        ``None`` returns every bundled identifier.

    Returns
    -------
    tuple[str, ...]
        Sorted immutable tuple of matching identifiers. An unknown kind yields
        an empty tuple.
    """
    entries, _ = _index()
    identifiers = tuple(sorted(entries))
    if model_kind is None:
        return identifiers
    if model_kind == "multifluid":
        warnings.warn(
            "model_kind='multifluid' is deprecated; use 'multiparameter'",
            DeprecationWarning,
            stacklevel=2,
        )
        model_kind = "multiparameter"
    return tuple(
        identifier
        for identifier in identifiers
        if _load_builtin(identifier).model_kind == model_kind
    )

clear_parameter_caches

clear_parameter_caches() -> None

Clear parsed model-parameter YAML caches.

Use this after intentionally modifying a custom file in a long-running process. Existing model instances and tensors are unaffected.

Returns:

Type	Description
None	The next parameter load reparses its source.

Source code in src/torch_flash/database.py

def clear_parameter_caches() -> None:
    """Clear parsed model-parameter YAML caches.

    Use this after intentionally modifying a custom file in a long-running
    process. Existing model instances and tensors are unaffected.

    Returns
    -------
    None
        The next parameter load reparses its source.
    """
    _load_builtin.cache_clear()
    _load_path.cache_clear()
    _index.cache_clear()

load_model_parameters

load_model_parameters(
    source: ParameterSource,
) -> ModelParameterSet

Load a bundled identifier, custom YAML path, or explicit parameter set.

Parsed bundled and custom files are cached by identifier or resolved path. Call :func:clear_parameter_caches after intentionally modifying a custom file in a long-running process.

Parameters:

Name	Type	Description	Default
source	ParameterSource	Bundled identifier or alias, custom YAML path, or an explicit :class:ModelParameterSet.	required

Returns:

Type	Description
ModelParameterSet	Immutable validated parameter set.

Raises:

Type	Description
KeyError	If a string is neither a bundled identifier/alias nor a YAML path.
FileNotFoundError	If a custom YAML path does not exist.
ParameterDatabaseError	If the document does not satisfy the schema.

Source code in src/torch_flash/database.py

def load_model_parameters(source: ParameterSource) -> ModelParameterSet:
    """Load a bundled identifier, custom YAML path, or explicit parameter set.

    Parsed bundled and custom files are cached by identifier or resolved path.
    Call :func:`clear_parameter_caches` after intentionally modifying a custom
    file in a long-running process.

    Parameters
    ----------
    source:
        Bundled identifier or alias, custom YAML path, or an explicit
        :class:`ModelParameterSet`.

    Returns
    -------
    ModelParameterSet
        Immutable validated parameter set.

    Raises
    ------
    KeyError
        If a string is neither a bundled identifier/alias nor a YAML path.
    FileNotFoundError
        If a custom YAML path does not exist.
    ParameterDatabaseError
        If the document does not satisfy the schema.
    """
    if isinstance(source, ModelParameterSet):
        return source
    if isinstance(source, Path):
        return _load_path(str(source.expanduser().absolute()))
    _, aliases = _index()
    normalized = source.strip().lower()
    if normalized in aliases:
        if normalized.startswith("multifluid."):
            warnings.warn(
                "the 'multifluid.' parameter prefix is deprecated; use 'multiparameter.'",
                DeprecationWarning,
                stacklevel=2,
            )
        return _load_builtin(aliases[normalized])
    candidate = Path(source).expanduser()
    if candidate.suffix.lower() in (".yaml", ".yml"):
        return _load_path(str(candidate.absolute()))
    available = ", ".join(available_parameter_sets())
    raise KeyError(f"unknown model parameter set {source!r}; available: {available}")

cubic_interaction_parameters

cubic_interaction_parameters(
    components: ComponentSet, source: ParameterSource
) -> CubicInteractionParameters

Load general van der Waals one-fluid kij and lij matrices.

Parameters:

Name	Type	Description	Default
components	ComponentSet	Ordered component set defining output matrix axes, dtype, and device.	required
source	ParameterSource	Bundled identifier, YAML path, mapping, or loaded binary-interaction parameter set.	required

Returns:

Type	Description
CubicInteractionParameters	Symmetric matrices on the component set's dtype and device.

Raises:

Type	Description
KeyError	If the parameter inventory does not support a requested component.
ParameterDatabaseError	If model identity, units, defaults, component order, or pair records are invalid.

Notes

The YAML payload declares component_order, optional numeric defaults, and unordered first|second pair records. Missing pair fields use the declared defaults, which are zero when omitted.

Source code in src/torch_flash/parameters/cubic_interactions.py

def cubic_interaction_parameters(
    components: ComponentSet,
    source: ParameterSource,
) -> CubicInteractionParameters:
    """Load general van der Waals one-fluid ``kij`` and ``lij`` matrices.

    Parameters
    ----------
    components
        Ordered component set defining output matrix axes, dtype, and device.
    source
        Bundled identifier, YAML path, mapping, or loaded binary-interaction
        parameter set.

    Returns
    -------
    CubicInteractionParameters
        Symmetric matrices on the component set's dtype and device.

    Raises
    ------
    KeyError
        If the parameter inventory does not support a requested component.
    ParameterDatabaseError
        If model identity, units, defaults, component order, or pair records
        are invalid.

    Notes
    -----
    The YAML payload declares ``component_order``, optional numeric
    ``defaults``, and unordered ``first|second`` pair records. Missing pair
    fields use the declared defaults, which are zero when omitted.
    """
    loaded = load_model_parameters(source)
    if loaded.model_kind != "binary_interaction":
        raise ParameterDatabaseError(
            f"{loaded.identifier!r} is {loaded.model_kind!r}, not 'binary_interaction'"
        )
    if loaded.model != "cubic-vdw-one-fluid":
        raise ParameterDatabaseError(f"{loaded.identifier!r} model must be 'cubic-vdw-one-fluid'")
    if loaded.units.get("kij") != "dimensionless" or loaded.units.get("lij") != "dimensionless":
        raise ParameterDatabaseError(
            f"{loaded.identifier!r} must declare dimensionless kij and lij units"
        )
    parameters = loaded.parameters
    order = _component_order(parameters.get("component_order"), loaded.identifier)
    unsupported = sorted(set(components.names) - set(order))
    if unsupported:
        raise KeyError(
            f"{loaded.identifier} has no interaction parameters for: {', '.join(unsupported)}"
        )
    defaults = parameters.get("defaults", {})
    if not isinstance(defaults, Mapping):
        raise ParameterDatabaseError(f"{loaded.identifier} defaults must be a mapping")
    default_kij = _finite_default(defaults, "kij", loaded.identifier)
    default_lij = _finite_default(defaults, "lij", loaded.identifier)
    shape = (components.ncomponents, components.ncomponents)
    kij = torch.full(
        shape,
        default_kij,
        dtype=components.critical_temperature.dtype,
        device=components.critical_temperature.device,
    )
    lij = torch.full_like(kij, default_lij)
    kij.fill_diagonal_(0.0)
    lij.fill_diagonal_(0.0)
    pairs = parameters.get("pairs", {})
    if not isinstance(pairs, Mapping):
        raise ParameterDatabaseError(f"{loaded.identifier} pairs must be a mapping")
    seen: set[frozenset[str]] = set()
    selected = {name: index for index, name in enumerate(components.names)}
    allowed = set(order)
    for key, record in pairs.items():
        if not isinstance(key, str) or key.count("|") != 1:
            raise ParameterDatabaseError(
                f"{loaded.identifier} pair keys must have the form 'first|second'"
            )
        first, second = key.split("|")
        if first == second:
            raise ParameterDatabaseError(f"{loaded.identifier} cannot contain self-pairs")
        if first not in allowed or second not in allowed:
            raise ParameterDatabaseError(
                f"{loaded.identifier} pair {key!r} names a component outside component_order"
            )
        unordered = frozenset((first, second))
        if unordered in seen:
            raise ParameterDatabaseError(
                f"{loaded.identifier} contains duplicate pair {first}|{second}"
            )
        seen.add(unordered)
        if not isinstance(record, Mapping):
            raise ParameterDatabaseError(f"{loaded.identifier} pair {key!r} must be a mapping")
        values = {}
        for parameter, default in (("kij", default_kij), ("lij", default_lij)):
            value = record.get(parameter, default)
            if not isinstance(value, int | float) or not math.isfinite(value):
                raise ParameterDatabaseError(
                    f"{loaded.identifier} pair {key!r} {parameter} must be finite and numeric"
                )
            values[parameter] = float(value)
        if first not in selected or second not in selected:
            continue
        first_index = selected[first]
        second_index = selected[second]
        kij[first_index, second_index] = kij[second_index, first_index] = values["kij"]
        lij[first_index, second_index] = lij[second_index, first_index] = values["lij"]
    return CubicInteractionParameters(kij, lij, loaded.identifier)

ppr78_group_contribution_parameters

ppr78_group_contribution_parameters(
    components: ComponentSet,
    source: ParameterSource = DEFAULT_PPR78_GROUP_CONTRIBUTION,
    *,
    group_counts: Mapping[str, Mapping[str, float]]
    | Tensor
    | None = None,
) -> PPR78GroupContributionParameters

Load PPR78 or E-PPR78 group parameters for components.

Parameters:

Name	Type	Description	Default
components	ComponentSet	Ordered component set defining selected decompositions, dtype, and device.	required
source	ParameterSource	Bundled identifier, custom YAML path, mapping, or loaded parameter set.	DEFAULT_PPR78_GROUP_CONTRIBUTION
group_counts	Mapping[str, Mapping[str, float]] | Tensor | None	Optional component-by-group counts as a mapping or tensor. Values are normalized into group fractions and override only source component decompositions.	None

Returns:

Type	Description
PPR78GroupContributionParameters	Selected fractions and universal interaction tensors. Parameter sets with explicitly unavailable pairs are reduced to the groups active in components after availability has been checked.

Raises:

Type	Description
KeyError	If a requested component lacks a decomposition.
ValueError	If explicit counts have invalid shapes or values.
ParameterDatabaseError	If source identity, units, group inventory, interactions, or the availability of an active group pair is invalid.

Notes

source may be the default 2004 saturated-hydrocarbon PPR78 set, :data:PPR78_HYDROGEN_WATER_GROUP_CONTRIBUTION, a custom YAML path, or an in-memory :class:~torch_flash.database.ModelParameterSet. The global E-PPR78 parameterization is selected with :data:DEFAULT_EPPR78_GROUP_CONTRIBUTION. The H2/N2/H2O set is the exact active-group submatrix of the 2013 hydrogen and water extensions, not a refit. Explicit group_counts override only the component decompositions; the source still supplies the group inventory and interaction matrices.

Source code in src/torch_flash/parameters/group_contribution.py

def ppr78_group_contribution_parameters(
    components: ComponentSet,
    source: ParameterSource = DEFAULT_PPR78_GROUP_CONTRIBUTION,
    *,
    group_counts: Mapping[str, Mapping[str, float]] | Tensor | None = None,
) -> PPR78GroupContributionParameters:
    """Load PPR78 or E-PPR78 group parameters for ``components``.

    Parameters
    ----------
    components
        Ordered component set defining selected decompositions, dtype, and
        device.
    source
        Bundled identifier, custom YAML path, mapping, or loaded parameter set.
    group_counts
        Optional component-by-group counts as a mapping or tensor. Values are
        normalized into group fractions and override only source component
        decompositions.

    Returns
    -------
    PPR78GroupContributionParameters
        Selected fractions and universal interaction tensors. Parameter sets
        with explicitly unavailable pairs are reduced to the groups active in
        ``components`` after availability has been checked.

    Raises
    ------
    KeyError
        If a requested component lacks a decomposition.
    ValueError
        If explicit counts have invalid shapes or values.
    ParameterDatabaseError
        If source identity, units, group inventory, interactions, or the
        availability of an active group pair is invalid.

    Notes
    -----
    ``source`` may be the default 2004 saturated-hydrocarbon PPR78 set,
    :data:`PPR78_HYDROGEN_WATER_GROUP_CONTRIBUTION`, a custom YAML path, or an
    in-memory :class:`~torch_flash.database.ModelParameterSet`. The global
    E-PPR78 parameterization is selected with
    :data:`DEFAULT_EPPR78_GROUP_CONTRIBUTION`. The
    H2/N2/H2O set is the exact active-group submatrix of the 2013 hydrogen and
    water extensions, not a refit. Explicit ``group_counts`` override only
    the component decompositions; the source still supplies the group
    inventory and interaction matrices.
    """
    loaded = load_model_parameters(source)
    if loaded.model_kind != "group_contribution":
        raise ParameterDatabaseError(
            f"{loaded.identifier!r} is {loaded.model_kind!r}, not 'group_contribution'"
        )
    if loaded.model not in {"PPR78", "E-PPR78"}:
        raise ParameterDatabaseError(f"{loaded.identifier!r} model must be 'PPR78' or 'E-PPR78'")
    expected_units = {
        "A": "Pa",
        "B": "Pa",
        "group_fraction": "dimensionless",
        "reference_temperature": "K",
    }
    if any(loaded.units.get(key) != unit for key, unit in expected_units.items()):
        raise ParameterDatabaseError(
            f"{loaded.identifier!r} must declare PPR78/E-PPR78 A, B, fraction, "
            "and temperature units"
        )
    parameters = loaded.parameters
    groups = _group_names(parameters.get("groups"), loaded.identifier)
    reference_temperature = _finite(
        parameters.get("reference_temperature"),
        f"{loaded.identifier} reference_temperature",
    )
    if reference_temperature <= 0.0:
        raise ParameterDatabaseError(f"{loaded.identifier} reference_temperature must be positive")
    group_a, group_b, unavailable = _interaction_matrices(
        parameters.get("interactions"),
        parameters.get("unavailable_interactions"),
        groups,
        loaded.identifier,
        dtype=components.critical_temperature.dtype,
        device=components.critical_temperature.device,
    )
    decompositions = parameters.get("component_groups") if group_counts is None else group_counts
    fractions = _fraction_matrix(
        components,
        groups,
        decompositions,
        loaded.identifier if group_counts is None else "<api group_counts>",
    )
    groups, fractions, group_a, group_b = _select_available_active_groups(
        groups,
        fractions,
        group_a,
        group_b,
        unavailable,
        loaded.identifier,
    )
    return PPR78GroupContributionParameters(
        groups,
        fractions,
        group_a,
        group_b,
        reference_temperature,
        loaded.identifier,
    )

ppr78_mixing

ppr78_mixing(
    components: ComponentSet,
    source: ParameterSource = DEFAULT_PPR78_GROUP_CONTRIBUTION,
    *,
    group_counts: Mapping[str, Mapping[str, float]]
    | Tensor
    | None = None,
    trainable: bool = False,
) -> PPR78Mixing

Construct a differentiable PPR78 or E-PPR78 mixing rule.

Parameters:

Name	Type	Description	Default
components	ComponentSet	Ordered component set.	required
source	ParameterSource	PPR78 or E-PPR78 parameter source.	DEFAULT_PPR78_GROUP_CONTRIBUTION
group_counts	Mapping[str, Mapping[str, float]] | Tensor | None	Optional explicit component group counts overriding stored decompositions. Counts are normalized to molecular group fractions.	None
trainable	bool	Register the independent, active off-diagonal A/B group interactions as PyTorch parameters.	False

Returns:

Type	Description
PPR78Mixing	Mixing-rule module on the component set's dtype and device.

Raises:

Type	Description
KeyError	If a requested component lacks a stored decomposition.
ValueError	If explicit group counts or tensor values are invalid.
ParameterDatabaseError	If the source inventory is malformed or an interaction between two active E-PPR78 groups is unavailable.

Notes

The implemented correlation is Eq. (5) of Jaubert and Mutelet, Fluid Phase Equilibria 224 (2004) 285--304, doi:10.1016/j.fluid.2004.06.059. E-PPR78 uses the same equation with the global interactions in Jaubert, Qian, Lasala, and Privat, Fluid Phase Equilibria 560 (2022) 113456, doi:10.1016/j.fluid.2022.113456, Table S4. For a source with unavailable pairs, only the groups active in components are retained after every active pair has been validated. All tensors honor the component set's dtype and device.

Source code in src/torch_flash/parameters/group_contribution.py

def ppr78_mixing(
    components: ComponentSet,
    source: ParameterSource = DEFAULT_PPR78_GROUP_CONTRIBUTION,
    *,
    group_counts: Mapping[str, Mapping[str, float]] | Tensor | None = None,
    trainable: bool = False,
) -> PPR78Mixing:
    """Construct a differentiable PPR78 or E-PPR78 mixing rule.

    Parameters
    ----------
    components
        Ordered component set.
    source
        PPR78 or E-PPR78 parameter source.
    group_counts
        Optional explicit component group counts overriding stored
        decompositions. Counts are normalized to molecular group fractions.
    trainable
        Register the independent, active off-diagonal A/B group interactions
        as PyTorch parameters.

    Returns
    -------
    PPR78Mixing
        Mixing-rule module on the component set's dtype and device.

    Raises
    ------
    KeyError
        If a requested component lacks a stored decomposition.
    ValueError
        If explicit group counts or tensor values are invalid.
    ParameterDatabaseError
        If the source inventory is malformed or an interaction between two
        active E-PPR78 groups is unavailable.

    Notes
    -----
    The implemented correlation is Eq. (5) of Jaubert and Mutelet,
    *Fluid Phase Equilibria* 224 (2004) 285--304,
    doi:10.1016/j.fluid.2004.06.059. E-PPR78 uses the same equation with the
    global interactions in Jaubert, Qian, Lasala, and Privat,
    *Fluid Phase Equilibria* 560 (2022) 113456,
    doi:10.1016/j.fluid.2022.113456, Table S4. For a source with unavailable
    pairs, only the groups active in ``components`` are retained after every
    active pair has been validated. All tensors honor the component set's
    dtype and device.
    """
    parameters = ppr78_group_contribution_parameters(
        components,
        source,
        group_counts=group_counts,
    )
    return PPR78Mixing(
        parameters.group_fractions,
        parameters.group_a,
        parameters.group_b,
        reference_temperature=parameters.reference_temperature,
        trainable=trainable,
        parameter_set=parameters.parameter_set,
    )

binary_interaction

binary_interaction(
    components: ComponentSet,
    parameter_set: ParameterSource,
    eos: PetroleumEOS = "PR",
    *,
    aggregate_hydrocarbons: Iterable[str] = (),
) -> Tensor

Select a petroleum cubic-EOS binary-interaction matrix.

Parameters:

Name	Type	Description	Default
components	ComponentSet	Ordered component set defining matrix axes, dtype, and device.	required
parameter_set	ParameterSource	Bundled identifier, YAML path, mapping, or loaded petroleum binary-interaction parameter set.	required
eos	PetroleumEOS	Select the source's Peng--Robinson ("PR") or Soave--Redlich--Kwong ("SRK") table.	'PR'
aggregate_hydrocarbons	Iterable[str]	Component names in components that represent characterized hydrocarbon cuts. Each uses the source's declared aggregate hydrocarbon row; hydrocarbon-hydrocarbon interactions remain zero.	()

Returns:

Type	Description
Tensor	Symmetric dimensionless kij matrix with zero diagonal and shape (ncomponents, ncomponents).

Raises:

Type	Description
ValueError	If eos is unsupported or aggregate_hydrocarbons is a string instead of an iterable of strings, contains non-string items, or declares a nonhydrocarbon.
KeyError	If a requested component is outside the parameter inventory.
ParameterDatabaseError	If the source model identity or table structure is invalid.

Source code in src/torch_flash/parameters/petroleum.py

def binary_interaction(
    components: ComponentSet,
    parameter_set: ParameterSource,
    eos: PetroleumEOS = "PR",
    *,
    aggregate_hydrocarbons: Iterable[str] = (),
) -> Tensor:
    """Select a petroleum cubic-EOS binary-interaction matrix.

    Parameters
    ----------
    components
        Ordered component set defining matrix axes, dtype, and device.
    parameter_set
        Bundled identifier, YAML path, mapping, or loaded petroleum
        binary-interaction parameter set.
    eos
        Select the source's Peng--Robinson (``"PR"``) or
        Soave--Redlich--Kwong (``"SRK"``) table.
    aggregate_hydrocarbons
        Component names in ``components`` that represent characterized
        hydrocarbon cuts. Each uses the source's declared aggregate
        hydrocarbon row; hydrocarbon-hydrocarbon interactions remain zero.

    Returns
    -------
    Tensor
        Symmetric dimensionless ``kij`` matrix with zero diagonal and shape
        ``(ncomponents, ncomponents)``.

    Raises
    ------
    ValueError
        If ``eos`` is unsupported or ``aggregate_hydrocarbons`` is a string
        instead of an iterable of strings, contains non-string items, or
        declares a nonhydrocarbon.
    KeyError
        If a requested component is outside the parameter inventory.
    ParameterDatabaseError
        If the source model identity or table structure is invalid.
    """
    if eos not in ("PR", "SRK"):
        raise ValueError("eos must be 'PR' or 'SRK'")
    loaded = load_model_parameters(parameter_set)
    if loaded.model_kind != "binary_interaction":
        raise ParameterDatabaseError(
            f"{loaded.identifier!r} is {loaded.model_kind!r}, not 'binary_interaction'"
        )
    parameters = loaded.parameters
    hydrocarbons = _string_tuple(parameters.get("hydrocarbons"), "hydrocarbons")
    nonhydrocarbons = _string_tuple(parameters.get("nonhydrocarbons"), "nonhydrocarbons")
    if isinstance(aggregate_hydrocarbons, str):
        raise ValueError(
            "aggregate_hydrocarbons must be an iterable of component-name strings, not a string"
        )
    aggregate_items = tuple(aggregate_hydrocarbons)
    if any(not isinstance(name, str) for name in aggregate_items):
        raise ValueError("aggregate_hydrocarbons must contain only strings")
    aggregate = frozenset(aggregate_items)
    if not aggregate.issubset(components.names):
        missing = ", ".join(sorted(aggregate - set(components.names)))
        raise KeyError(f"aggregate hydrocarbon names are absent from components: {missing}")
    if aggregate & set(nonhydrocarbons):
        names = ", ".join(sorted(aggregate & set(nonhydrocarbons)))
        raise ValueError(f"nonhydrocarbons cannot be declared as aggregate hydrocarbons: {names}")
    supported = frozenset((*hydrocarbons, *nonhydrocarbons, *aggregate))
    unsupported = sorted(set(components.names) - supported)
    if unsupported:
        names = ", ".join(unsupported)
        raise KeyError(f"{loaded.identifier} has no parameters for: {names}")
    aggregate_from = parameters.get("aggregate_from")
    if not isinstance(aggregate_from, str) or aggregate_from not in hydrocarbons:
        raise ParameterDatabaseError(
            f"{loaded.identifier!r} requires a valid aggregate_from component"
        )
    aggregate_index = hydrocarbons.index(aggregate_from)
    eos_tables = parameters.get("eos")
    if not isinstance(eos_tables, Mapping):
        raise ParameterDatabaseError(f"{loaded.identifier!r} requires an eos mapping")
    table = eos_tables.get(eos)
    if not isinstance(table, Mapping):
        raise ParameterDatabaseError(f"{loaded.identifier!r} has no {eos} table")
    hydrocarbon_rows = table.get("hydrocarbon_rows")
    pair_values = table.get("nonhydrocarbon_pairs")
    if not isinstance(hydrocarbon_rows, Mapping) or not isinstance(pair_values, Mapping):
        raise ParameterDatabaseError(f"{loaded.identifier!r} has incomplete {eos} tables")

    kij = torch.zeros(
        (components.ncomponents, components.ncomponents),
        dtype=components.critical_temperature.dtype,
        device=components.critical_temperature.device,
    )
    for row, first in enumerate(components.names):
        for column in range(row):
            second = components.names[column]
            first_hydrocarbon = first in hydrocarbons or first in aggregate
            second_hydrocarbon = second in hydrocarbons or second in aggregate
            if first_hydrocarbon and second_hydrocarbon:
                value = 0.0
            elif first in nonhydrocarbons and second in nonhydrocarbons:
                key = f"{first}|{second}"
                reverse_key = f"{second}|{first}"
                raw_value = pair_values.get(key, pair_values.get(reverse_key))
                if not isinstance(raw_value, int | float):
                    raise ParameterDatabaseError(
                        f"{loaded.identifier!r} lacks pair {first}|{second}"
                    )
                value = float(raw_value)
            else:
                nonhydrocarbon = first if first in nonhydrocarbons else second
                hydrocarbon = second if first in nonhydrocarbons else first
                values = hydrocarbon_rows.get(nonhydrocarbon)
                if not isinstance(values, Sequence) or isinstance(values, str):
                    raise ParameterDatabaseError(
                        f"{loaded.identifier!r} lacks row {nonhydrocarbon!r}"
                    )
                hydrocarbon_index = (
                    aggregate_index
                    if hydrocarbon in aggregate
                    else min(hydrocarbons.index(hydrocarbon), aggregate_index)
                )
                raw_value = values[hydrocarbon_index]
                if not isinstance(raw_value, int | float):
                    raise ParameterDatabaseError("binary-interaction values must be numeric")
                value = float(raw_value)
            kij[row, column] = value
            kij[column, row] = value
    return kij

pedersen_binary_interaction

pedersen_binary_interaction(
    components: ComponentSet,
    eos: PetroleumEOS = "PR",
    *,
    aggregate_hydrocarbons: Iterable[str] = (),
) -> Tensor

Return Pedersen et al. (2024), Table 4.2, kij values.

Parameters:

Name	Type	Description	Default
components	ComponentSet	Ordered component set defining output axes.	required
eos	PetroleumEOS	Source table, "PR" or "SRK".	'PR'
aggregate_hydrocarbons	Iterable[str]	Explicit names of characterized hydrocarbon cuts that use the published C7+ interaction row.	()

Returns:

Type	Description
Tensor	Symmetric dimensionless interaction matrix.

Notes

This is a distinct parameterization from Whitson's Table A-3. The source's aggregate C7+ entries are used for n-heptane through n-decane.

Source code in src/torch_flash/parameters/petroleum.py

def pedersen_binary_interaction(
    components: ComponentSet,
    eos: PetroleumEOS = "PR",
    *,
    aggregate_hydrocarbons: Iterable[str] = (),
) -> Tensor:
    """Return Pedersen et al. (2024), Table 4.2, ``kij`` values.

    Parameters
    ----------
    components
        Ordered component set defining output axes.
    eos
        Source table, ``"PR"`` or ``"SRK"``.
    aggregate_hydrocarbons
        Explicit names of characterized hydrocarbon cuts that use the
        published C7+ interaction row.

    Returns
    -------
    Tensor
        Symmetric dimensionless interaction matrix.

    Notes
    -----
    This is a distinct parameterization from Whitson's Table A-3. The source's
    aggregate C7+ entries are used for n-heptane through n-decane.
    """
    return binary_interaction(
        components,
        "binary-interaction.pedersen-2024",
        eos,
        aggregate_hydrocarbons=aggregate_hydrocarbons,
    )

whitson_binary_interaction

whitson_binary_interaction(
    components: ComponentSet, eos: PetroleumEOS = "PR"
) -> Tensor

Return Whitson and Brule (2000), Table A-3, kij values.

Parameters:

Name	Type	Description	Default
components	ComponentSet	Ordered component set defining output axes.	required
eos	PetroleumEOS	Source table, "PR" or "SRK".	'PR'

Returns:

Type	Description
Tensor	Symmetric dimensionless interaction matrix.

Notes

The source's C7+ values are applied from n-heptane through n-decane. The printed H2S/C7+ value is not silently extrapolated by carbon number.

Source code in src/torch_flash/parameters/petroleum.py

def whitson_binary_interaction(
    components: ComponentSet,
    eos: PetroleumEOS = "PR",
) -> Tensor:
    """Return Whitson and Brule (2000), Table A-3, ``kij`` values.

    Parameters
    ----------
    components
        Ordered component set defining output axes.
    eos
        Source table, ``"PR"`` or ``"SRK"``.

    Returns
    -------
    Tensor
        Symmetric dimensionless interaction matrix.

    Notes
    -----
    The source's C7+ values are applied from n-heptane through n-decane. The
    printed H2S/C7+ value is not silently extrapolated by carbon number.
    """
    return binary_interaction(components, "binary-interaction.whitson-2000", eos)

Parameter fitting

Calibration utilities preserve PyTorch gradients but do not by themselves establish identifiability or independent validation.

torch_flash.fitting

Autodifferentiable thermodynamic parameter estimation.

FitResult dataclass

FitResult(
    losses: tuple[float, ...],
    converged: bool,
    iterations: int,
    final_loss: float,
    stopping_reason: Literal[
        "tolerance", "patience", "iteration-limit"
    ],
)

Optimization history and final convergence diagnostics.

Attributes:

Name	Type	Description
losses	tuple[float, ...]	Scalar objective at the initial point followed by each accepted optimizer iterate.
converged	bool	Whether successive objective values met the relative stopping test.
iterations	int	Number of optimizer iterations executed.
final_loss	float	Minimum objective value retained in the supplied parameter tensors.
stopping_reason	Literal['tolerance', 'patience', 'iteration-limit']	"tolerance", "patience", or "iteration-limit".

Notes

Optimizer convergence does not establish parameter identifiability or model validation. Inspect sensitivities, parameter correlations, and independent holdout behavior separately.

CubicInteractionFitSystem dataclass

CubicInteractionFitSystem(
    components: ComponentSet,
    component_indices: tuple[int, ...],
    states: tuple[PhaseTransitionState, ...],
)

One component subset and transition-state group in a cubic fit.

Attributes:

Name	Type	Description
components	ComponentSet	Ordered component set used to construct this system's cubic model.
component_indices	tuple[int, ...]	Indices mapping components into the master interaction matrix.
states	tuple[PhaseTransitionState, ...]	Phase-transition states used by the joint fit and sensitivity audit.

CubicPhaseTransitionFitResult dataclass

CubicPhaseTransitionFitResult(
    parameters: Tensor,
    parameter_history: tuple[float, ...],
    selected_iteration: int,
    selected_loss: float,
    optimizer_converged: bool,
    iterations: int,
    optimizer_stopping_reason: Literal[
        "tolerance", "patience", "iteration-limit"
    ],
    seed_evaluations: tuple[
        tuple[PhaseTransitionEvaluation, ...], ...
    ],
    sensitivity_matrix: Tensor,
    sensitivity_singular_values: Tensor,
    sensitivity_rank: int,
    sensitivity_condition_number: float,
    calibration_objective: CalibrationObjective,
    optimizer: OptimizerKind,
    sensitivity_kind: Literal[
        "log-transition-pressure", "observed-state-fugacity"
    ],
    fitted_states: tuple[
        tuple[PhaseTransitionState, ...], ...
    ],
)

Bounded cubic-interaction fit and local identifiability diagnostics.

Attributes:

Name	Type	Description
parameters	Tensor	Selected physical interaction vector in the declared pair order.
parameter_history	tuple[float, ...]	Scalar objective at the initial point followed by each accepted optimizer iterate.
selected_iteration	int	Zero-based optimizer-step index of the retained minimum-loss point. Zero identifies the initial point.
selected_loss	float	Loss at parameters.
optimizer_converged	bool	Whether the strict successive-loss stopping test was met.
iterations	int	Number of Adam iterations executed.
optimizer_stopping_reason	Literal['tolerance', 'patience', 'iteration-limit']	"tolerance", "patience", or "iteration-limit".
seed_evaluations	tuple[tuple[PhaseTransitionEvaluation, ...], ...]	Source-parameter branch evaluations in system/state order.
sensitivity_matrix	Tensor	Jacobian of fitted log pressure ratios with respect to physical interaction parameters. Nonfinite rows are excluded.
sensitivity_singular_values	Tensor	Singular values of sensitivity_matrix.
sensitivity_rank	int	Numerical matrix rank using dtype-aware SVD tolerance.
sensitivity_condition_number	float	Largest/smallest singular-value ratio, or infinity when rank deficient.
calibration_objective	CalibrationObjective	Residual formulation minimized during calibration.
optimizer	OptimizerKind	PyTorch optimizer used for parameter updates.
sensitivity_kind	Literal['log-transition-pressure', 'observed-state-fugacity']	Scientific meaning of the rows in sensitivity_matrix.
fitted_states	tuple[tuple[PhaseTransitionState, ...], ...]	Phase-transition states in input system/state order with the selected incipient branch compositions installed as explicit initializers.

build_cubic_interaction_models

build_cubic_interaction_models(
    constructor: Callable[..., StateModel],
    systems: Sequence[CubicInteractionFitSystem],
    parameters: Tensor,
    *,
    kij_pairs: Sequence[tuple[int, int]],
    lij_pairs: Sequence[tuple[int, int]] = (),
) -> tuple[StateModel, ...]

Construct matched cubic systems from one shared interaction vector.

Parameters:

Name	Type	Description	Default
constructor	Callable[..., StateModel]	Cubic model constructor accepting components, kij, and lij.	required
systems	Sequence[CubicInteractionFitSystem]	Component subsets and their master-component index mappings.	required
parameters	Tensor	One-dimensional tensor containing all kij_pairs values followed by all lij_pairs values.	required
kij_pairs	Sequence[tuple[int, int]]	Unique off-diagonal master-component index pairs. Returned models use symmetric zero-diagonal matrices.	required
lij_pairs	Sequence[tuple[int, int]]	Unique off-diagonal master-component index pairs. Returned models use symmetric zero-diagonal matrices.	required

Returns:

Type	Description
tuple	One constructed model per system in input order. Parameter gradients, dtype, and device are preserved.

Raises:

Type	Description
ValueError	If parameters, system mappings, or interaction pairs are invalid.

Source code in src/torch_flash/fitting.py

def build_cubic_interaction_models(
    constructor: Callable[..., StateModel],
    systems: Sequence[CubicInteractionFitSystem],
    parameters: Tensor,
    *,
    kij_pairs: Sequence[tuple[int, int]],
    lij_pairs: Sequence[tuple[int, int]] = (),
) -> tuple[StateModel, ...]:
    """Construct matched cubic systems from one shared interaction vector.

    Parameters
    ----------
    constructor
        Cubic model constructor accepting ``components``, ``kij``, and ``lij``.
    systems
        Component subsets and their master-component index mappings.
    parameters
        One-dimensional tensor containing all ``kij_pairs`` values followed
        by all ``lij_pairs`` values.
    kij_pairs, lij_pairs
        Unique off-diagonal master-component index pairs. Returned models use
        symmetric zero-diagonal matrices.

    Returns
    -------
    tuple
        One constructed model per system in input order. Parameter gradients,
        dtype, and device are preserved.

    Raises
    ------
    ValueError
        If parameters, system mappings, or interaction pairs are invalid.
    """
    if parameters.ndim != 1 or parameters.numel() != len(kij_pairs) + len(lij_pairs):
        raise ValueError("interaction parameter vector does not match declared pairs")
    if not systems:
        raise ValueError("at least one cubic interaction fit system is required")
    component_count = max(max(system.component_indices) for system in systems) + 1
    all_pairs = (*kij_pairs, *lij_pairs)
    if any(
        first < 0
        or second < 0
        or first >= component_count
        or second >= component_count
        or first == second
        for first, second in all_pairs
    ):
        raise ValueError("interaction pairs must be distinct master-component indices")
    if len(set(kij_pairs)) != len(kij_pairs) or len(set(lij_pairs)) != len(lij_pairs):
        raise ValueError("interaction pairs must be unique within each interaction matrix")
    for system in systems:
        if len(system.component_indices) != len(system.components.names) or len(
            set(system.component_indices)
        ) != len(system.component_indices):
            raise ValueError("system component indices must uniquely match its component set")

    kij = _interaction_matrix(parameters[: len(kij_pairs)], kij_pairs, component_count)
    lij = _interaction_matrix(parameters[len(kij_pairs) :], lij_pairs, component_count)
    models = []
    for system in systems:
        indices = torch.tensor(
            system.component_indices,
            dtype=torch.long,
            device=parameters.device,
        )
        models.append(
            constructor(
                system.components,
                kij=kij[indices][:, indices],
                lij=lij[indices][:, indices],
            )
        )
    return tuple(models)

fit_cubic_phase_transition_interactions

fit_cubic_phase_transition_interactions(
    constructor: Callable[..., StateModel],
    systems: Sequence[CubicInteractionFitSystem],
    initial_parameters: Tensor,
    lower_bounds: Tensor,
    upper_bounds: Tensor,
    *,
    kij_pairs: Sequence[tuple[int, int]],
    lij_pairs: Sequence[tuple[int, int]] = (),
    objective: CalibrationObjective = "transition-pressure",
    optimizer: OptimizerKind = "adam",
    learning_rate: float = 0.005,
    weight_decay: float = 0.0,
    momentum: float = 0.0,
    max_iterations: int = 30,
    stopping_tolerance: float = 1e-09,
    no_improvement_patience: int | None = 6,
    fugacity_tolerance: float = 1e-07,
    solver_iterations: int = 24,
    minimum_phase_separation: float = 0.002,
    fugacity_penalty: float = 0.05,
    phase_merge_penalty: float = 5.0,
    parameter_prior_weight: float = 0.0,
) -> CubicPhaseTransitionFitResult

Fit shared conventional-cubic interactions to phase-transition states.

Parameters:

Name	Type	Description	Default
constructor	Callable[..., StateModel]	Cubic model constructor accepting component, kij, and lij arguments.	required
systems	Sequence[CubicInteractionFitSystem]	Component subsets, master-index mappings, and two-phase transition states. Two-phase boundaries and binary three-phase invariants are included in one joint objective.	required
initial_parameters	Tensor	Interior physical interaction vector in kij_pairs then lij_pairs order.	required
lower_bounds	Tensor	Finite physical bounds with the same shape as initial_parameters.	required
upper_bounds	Tensor	Finite physical bounds with the same shape as initial_parameters.	required
kij_pairs	Sequence[tuple[int, int]]	Shared master-component pair identities.	required
lij_pairs	Sequence[tuple[int, int]]	Shared master-component pair identities.	required
objective	CalibrationObjective	"transition-pressure" solves a local transition pressure at every trial parameter vector and minimizes squared log-pressure ratios. "observed-state-fugacity" simultaneously optimizes the shared interactions and one physically separated incipient composition per two-phase observation at its measured temperature and pressure.	'transition-pressure'
optimizer	OptimizerKind	PyTorch optimizer: "adam", "adamw", "sgd", or strong-Wolfe "lbfgs".	'adam'
learning_rate	float	Optimizer learning rate in unconstrained coordinates.	0.005
weight_decay	float	Nonnegative optimizer weight decay. Keep zero unless a regularization prior on the transformed variables is scientifically justified.	0.0
momentum	float	SGD momentum. Ignored by Adam and AdamW.	0.0
max_iterations	int	Maximum optimizer iterations.	30
stopping_tolerance	float	Relative successive-loss stopping threshold.	1e-09
no_improvement_patience	int | None	Stop after this many evaluated losses without a new minimum. None disables this criterion. The best evaluated iterate is retained in either case.	6
fugacity_tolerance	float	Local transition-pressure solver controls. They apply to the "transition-pressure" objective and the post-fit pressure sensitivity audit.	1e-07
solver_iterations	float	Local transition-pressure solver controls. They apply to the "transition-pressure" objective and the post-fit pressure sensitivity audit.	1e-07
minimum_phase_separation	float	Minimum selected-component mole-fraction difference. The simultaneous observed-state objective enforces it by construction, rather than allowing the algebraic homogeneous solution.	0.002
fugacity_penalty	float	Weight on mean squared dimensionless fugacity residual.	0.05
phase_merge_penalty	float	Weight on squared phase-separation shortfall.	5.0
parameter_prior_weight	float	Nonnegative quadratic prior weight on the bounded physical interaction displacement from initial_parameters, scaled by each half-range. Zero disables the prior.	0.0

Returns:

Type	Description
CubicPhaseTransitionFitResult	Minimum-loss bounded parameters, complete optimization history, source-seed branch results, and a local sensitivity SVD.

Raises:

Type	Description
ValueError	If bounds, parameters, pair identities, states, objective, optimizer, or numerical controls are invalid, or the selected objective cannot represent a requested state.

Notes

Branches are identified once at initial_parameters using :func:torch_flash.envelope.evaluate_phase_transition_state; the selected or fallback physically separated incipient composition initializes the same branch during autodiff optimization. The simultaneous observed-state formulation evaluates every experimental state in every optimizer step, avoids differentiating through a failed nested pressure solve, and keeps the latent incipient compositions strictly separated from the parent.

The minimum-loss evaluated iterate is retained even when a later first-order step increases the objective. SVD rank and condition number diagnose only the reported local residual Jacobian; they do not establish a unique global fit.

Source code in src/torch_flash/fitting.py

def fit_cubic_phase_transition_interactions(
    constructor: Callable[..., StateModel],
    systems: Sequence[CubicInteractionFitSystem],
    initial_parameters: Tensor,
    lower_bounds: Tensor,
    upper_bounds: Tensor,
    *,
    kij_pairs: Sequence[tuple[int, int]],
    lij_pairs: Sequence[tuple[int, int]] = (),
    objective: CalibrationObjective = "transition-pressure",
    optimizer: OptimizerKind = "adam",
    learning_rate: float = 5.0e-3,
    weight_decay: float = 0.0,
    momentum: float = 0.0,
    max_iterations: int = 30,
    stopping_tolerance: float = 1.0e-9,
    no_improvement_patience: int | None = 6,
    fugacity_tolerance: float = 1.0e-7,
    solver_iterations: int = 24,
    minimum_phase_separation: float = 2.0e-3,
    fugacity_penalty: float = 0.05,
    phase_merge_penalty: float = 5.0,
    parameter_prior_weight: float = 0.0,
) -> CubicPhaseTransitionFitResult:
    """Fit shared conventional-cubic interactions to phase-transition states.

    Parameters
    ----------
    constructor
        Cubic model constructor accepting component, ``kij``, and ``lij``
        arguments.
    systems
        Component subsets, master-index mappings, and two-phase transition
        states. Two-phase boundaries and binary three-phase invariants are
        included in one joint objective.
    initial_parameters
        Interior physical interaction vector in ``kij_pairs`` then
        ``lij_pairs`` order.
    lower_bounds, upper_bounds
        Finite physical bounds with the same shape as ``initial_parameters``.
    kij_pairs, lij_pairs
        Shared master-component pair identities.
    objective
        ``"transition-pressure"`` solves a local transition pressure at every
        trial parameter vector and minimizes squared log-pressure ratios.
        ``"observed-state-fugacity"`` simultaneously optimizes the shared
        interactions and one physically separated incipient composition per
        two-phase observation at its measured temperature and pressure.
    optimizer
        PyTorch optimizer: ``"adam"``, ``"adamw"``, ``"sgd"``, or
        strong-Wolfe ``"lbfgs"``.
    learning_rate
        Optimizer learning rate in unconstrained coordinates.
    weight_decay
        Nonnegative optimizer weight decay. Keep zero unless a regularization
        prior on the transformed variables is scientifically justified.
    momentum
        SGD momentum. Ignored by Adam and AdamW.
    max_iterations
        Maximum optimizer iterations.
    stopping_tolerance
        Relative successive-loss stopping threshold.
    no_improvement_patience
        Stop after this many evaluated losses without a new minimum. ``None``
        disables this criterion. The best evaluated iterate is retained in
        either case.
    fugacity_tolerance, solver_iterations
        Local transition-pressure solver controls. They apply to the
        ``"transition-pressure"`` objective and the post-fit pressure
        sensitivity audit.
    minimum_phase_separation
        Minimum selected-component mole-fraction difference. The simultaneous
        observed-state objective enforces it by construction, rather than
        allowing the algebraic homogeneous solution.
    fugacity_penalty
        Weight on mean squared dimensionless fugacity residual.
    phase_merge_penalty
        Weight on squared phase-separation shortfall.
    parameter_prior_weight
        Nonnegative quadratic prior weight on the bounded physical interaction
        displacement from ``initial_parameters``, scaled by each half-range.
        Zero disables the prior.

    Returns
    -------
    CubicPhaseTransitionFitResult
        Minimum-loss bounded parameters, complete optimization history,
        source-seed branch results, and a local sensitivity SVD.

    Raises
    ------
    ValueError
        If bounds, parameters, pair identities, states, objective, optimizer,
        or numerical controls are invalid, or the selected objective cannot
        represent a requested state.

    Notes
    -----
    Branches are identified once at ``initial_parameters`` using
    :func:`torch_flash.envelope.evaluate_phase_transition_state`; the selected
    or fallback physically separated incipient composition initializes the
    same branch during autodiff optimization. The simultaneous observed-state
    formulation evaluates every experimental state in every optimizer step,
    avoids differentiating through a failed nested pressure solve, and keeps
    the latent incipient compositions strictly separated from the parent.

    The minimum-loss evaluated iterate is retained even when a later
    first-order step increases the objective. SVD rank and condition number
    diagnose only the reported local residual Jacobian; they do not establish
    a unique global fit.
    """
    if (
        initial_parameters.ndim != 1
        or lower_bounds.shape != initial_parameters.shape
        or upper_bounds.shape != initial_parameters.shape
    ):
        raise ValueError("initial cubic interactions and bounds must be equally shaped vectors")
    if not bool(
        torch.isfinite(initial_parameters).all()
        & torch.isfinite(lower_bounds).all()
        & torch.isfinite(upper_bounds).all()
        & (lower_bounds < initial_parameters).all()
        & (initial_parameters < upper_bounds).all()
    ):
        raise ValueError("initial cubic interactions must lie strictly inside finite bounds")
    if (
        learning_rate <= 0.0
        or max_iterations <= 0
        or stopping_tolerance <= 0.0
        or (no_improvement_patience is not None and no_improvement_patience <= 0)
        or fugacity_tolerance <= 0.0
        or solver_iterations <= 0
        or minimum_phase_separation < 0.0
        or fugacity_penalty < 0.0
        or phase_merge_penalty < 0.0
        or parameter_prior_weight < 0.0
        or weight_decay < 0.0
        or momentum < 0.0
    ):
        raise ValueError("cubic phase-transition fitting controls are invalid")
    if objective not in ("transition-pressure", "observed-state-fugacity"):
        raise ValueError("unknown cubic phase-transition calibration objective")
    if optimizer not in ("adam", "adamw", "sgd", "lbfgs"):
        raise ValueError("unknown PyTorch fitting optimizer")
    if objective != "transition-pressure" and any(
        state.boundary_kind == "liquid-liquid-vapor"
        for system in systems
        for state in system.states
    ):
        raise ValueError(
            "simultaneous observed-state calibration currently accepts only two-phase states"
        )

    midpoint = 0.5 * (lower_bounds + upper_bounds)
    half_width = 0.5 * (upper_bounds - lower_bounds)

    def to_physical(raw: Tensor) -> Tensor:
        return midpoint + half_width * torch.tanh(raw)

    normalized = torch.clamp(
        (initial_parameters - midpoint) / half_width,
        -1.0 + 1.0e-10,
        1.0 - 1.0e-10,
    )
    raw = nn.Parameter(torch.atanh(normalized).detach().clone())
    initial_models = build_cubic_interaction_models(
        constructor,
        systems,
        initial_parameters,
        kij_pairs=kij_pairs,
        lij_pairs=lij_pairs,
    )
    seed_evaluations = tuple(
        evaluate_phase_transition_states(
            model,
            system.states,
            exhaustive_two_phase_starts=True,
        )
        for model, system in zip(initial_models, systems, strict=True)
    )
    if any(
        state.boundary_kind == "liquid-liquid-vapor" and not evaluation.converged
        for system, system_evaluations in zip(
            systems,
            seed_evaluations,
            strict=True,
        )
        for state, evaluation in zip(
            system.states,
            system_evaluations,
            strict=True,
        )
    ):
        raise ValueError(
            "three-phase fit states require residual-converged, non-degenerate seed branches"
        )
    fixed_systems = tuple(
        CubicInteractionFitSystem(
            system.components,
            system.component_indices,
            tuple(
                replace(
                    state,
                    initial_incipient_compositions=(
                        ()
                        if state.boundary_kind == "liquid-liquid-vapor"
                        else (
                            _separated_two_phase_seed(
                                state,
                                evaluation,
                                minimum_phase_separation=minimum_phase_separation,
                            ),
                        )
                    ),
                    initial_three_phase_compositions=(
                        (evaluation.phase_compositions,)
                        if state.boundary_kind == "liquid-liquid-vapor"
                        else ()
                    ),
                )
                for state, evaluation in zip(
                    system.states,
                    system_evaluations,
                    strict=True,
                )
            ),
        )
        for system, system_evaluations in zip(systems, seed_evaluations, strict=True)
    )
    fit_batches = _build_cubic_transition_batches(fixed_systems)
    observed_batches = (
        _build_cubic_observed_state_batches(
            fixed_systems,
            minimum_phase_separation=minimum_phase_separation,
        )
        if objective != "transition-pressure"
        else ()
    )
    raw_compositions = tuple(
        nn.Parameter(
            _encode_separated_composition(
                batch,
                minimum_phase_separation=minimum_phase_separation,
            )
        )
        for batch in observed_batches
    )

    def loss_at(raw_parameters: Tensor, compositions: Sequence[Tensor]) -> Tensor:
        physical_parameters = to_physical(raw_parameters)
        if objective == "observed-state-fugacity":
            loss = _cubic_observed_state_fugacity_loss(
                constructor,
                fixed_systems,
                observed_batches,
                physical_parameters,
                compositions,
                kij_pairs=kij_pairs,
                lij_pairs=lij_pairs,
                minimum_phase_separation=minimum_phase_separation,
            )
        else:
            loss = _cubic_transition_loss(
                constructor,
                fixed_systems,
                fit_batches,
                physical_parameters,
                kij_pairs=kij_pairs,
                lij_pairs=lij_pairs,
                fugacity_tolerance=fugacity_tolerance,
                solver_iterations=solver_iterations,
                minimum_phase_separation=minimum_phase_separation,
                fugacity_penalty=fugacity_penalty,
                phase_merge_penalty=phase_merge_penalty,
            )
        parameter_prior = ((physical_parameters - initial_parameters) / half_width).square().mean()
        return loss + parameter_prior_weight * parameter_prior

    optimization = fit_parameters(
        (raw, *raw_compositions),
        lambda: loss_at(raw, raw_compositions),
        optimizer=optimizer,
        learning_rate=learning_rate,
        weight_decay=weight_decay,
        momentum=momentum,
        max_iterations=max_iterations,
        tolerance=stopping_tolerance,
        no_improvement_patience=no_improvement_patience,
    )
    selected_index = min(
        range(len(optimization.losses)),
        key=optimization.losses.__getitem__,
    )
    selected_parameters = to_physical(raw).detach()
    selected_raw_compositions = tuple(
        raw_composition.detach().clone() for raw_composition in raw_compositions
    )
    if objective != "transition-pressure":
        decoded_compositions = tuple(
            _decode_separated_composition(
                batch,
                raw_composition,
                minimum_phase_separation=minimum_phase_separation,
            ).detach()
            for batch, raw_composition in zip(
                observed_batches,
                selected_raw_compositions,
                strict=True,
            )
        )
        decoded_by_group = {
            (batch.system_index, batch.phase_kinds): compositions
            for batch, compositions in zip(
                observed_batches,
                decoded_compositions,
                strict=True,
            )
        }
        group_positions: dict[tuple[int, tuple[PhaseKind, PhaseKind]], int] = {}
        fitted_states_rows = []
        for system_index, system in enumerate(fixed_systems):
            system_states = []
            for state in system.states:
                phase_kinds: tuple[PhaseKind, PhaseKind] = (
                    ("liquid", "liquid")
                    if state.boundary_kind == "liquid-liquid"
                    else ("liquid", "vapor")
                )
                group_key = (system_index, phase_kinds)
                group_position = group_positions.get(group_key, 0)
                composition = decoded_by_group[group_key][group_position]
                group_positions[group_key] = group_position + 1
                system_states.append(
                    replace(
                        state,
                        initial_incipient_compositions=(composition,),
                    )
                )
            fitted_states = tuple(system_states)
            fitted_states_rows.append(fitted_states)
        result_fitted_states = tuple(fitted_states_rows)
    else:
        result_fitted_states = tuple(system.states for system in fixed_systems)

    sensitivity_parameters = selected_parameters.clone().requires_grad_(True)
    if objective != "transition-pressure":
        residual_rows = torch.cat(
            tuple(
                (residual / residual.shape[-1] ** 0.5).reshape(-1)
                for residual in _cubic_observed_state_fugacity_residuals(
                    constructor,
                    fixed_systems,
                    observed_batches,
                    sensitivity_parameters,
                    selected_raw_compositions,
                    kij_pairs=kij_pairs,
                    lij_pairs=lij_pairs,
                    minimum_phase_separation=minimum_phase_separation,
                )
            )
        )
        sensitivity_outputs = residual_rows
        sensitivity_kind: Literal[
            "log-transition-pressure",
            "observed-state-fugacity",
        ] = "observed-state-fugacity"
    else:
        sensitivity_models = build_cubic_interaction_models(
            constructor,
            fixed_systems,
            sensitivity_parameters,
            kij_pairs=kij_pairs,
            lij_pairs=lij_pairs,
        )
        sensitivity_evaluations = _evaluate_cubic_transition_batches(
            sensitivity_models,
            fit_batches,
            tolerance=fugacity_tolerance,
            max_iterations=solver_iterations,
            minimum_phase_separation=minimum_phase_separation,
        )
        sensitivity_outputs = torch.cat(
            tuple(
                torch.log(evaluation.pressure / batch.reference_pressure)
                for batch, evaluation in zip(
                    fit_batches,
                    sensitivity_evaluations,
                    strict=True,
                )
            )
        )
        sensitivity_kind = "log-transition-pressure"
    finite_sensitivity_outputs = torch.isfinite(sensitivity_outputs.detach())
    if bool(finite_sensitivity_outputs.any()):
        finite_outputs = sensitivity_outputs[finite_sensitivity_outputs]
        sensitivity = torch.autograd.grad(
            finite_outputs,
            sensitivity_parameters,
            grad_outputs=torch.eye(
                finite_outputs.numel(),
                dtype=finite_outputs.dtype,
                device=finite_outputs.device,
            ),
            is_grads_batched=True,
        )[0].detach()
        sensitivity = sensitivity[torch.isfinite(sensitivity).all(dim=-1)]
    else:
        sensitivity = selected_parameters.new_empty((0, selected_parameters.numel()))
    singular_values = (
        torch.linalg.svdvals(sensitivity)
        if sensitivity.numel()
        else selected_parameters.new_empty((0,))
    )
    if singular_values.numel():
        threshold = singular_values[0] * max(sensitivity.shape) * torch.finfo(sensitivity.dtype).eps
        rank = int((singular_values > threshold).sum())
        condition = (
            float(singular_values[0] / singular_values[-1])
            if rank == selected_parameters.numel()
            else float("inf")
        )
    else:
        rank = 0
        condition = float("inf")
    return CubicPhaseTransitionFitResult(
        selected_parameters,
        optimization.losses,
        selected_index,
        optimization.losses[selected_index],
        optimization.converged,
        optimization.iterations,
        optimization.stopping_reason,
        seed_evaluations,
        sensitivity,
        singular_values,
        rank,
        condition,
        objective,
        optimizer,
        sensitivity_kind,
        result_fitted_states,
    )

least_squares_loss

least_squares_loss(
    prediction: Tensor,
    observation: Tensor,
    *,
    scale: Tensor | float = 1.0,
    weights: Tensor | None = None,
) -> Tensor

Return a dimensionless weighted mean-square residual.

Parameters:

Name	Type	Description	Default
prediction	Tensor	Broadcast-compatible model predictions and observations in matching physical units.	required
observation	Tensor	Broadcast-compatible model predictions and observations in matching physical units.	required
scale	Tensor | float	Positive residual scale in the same units. Scalar or broadcastable tensor.	1.0
weights	Tensor | None	Optional nonnegative statistical weights, broadcastable to the residual shape.	None

Returns:

Type	Description
Tensor	Scalar mean of the squared, scaled residuals; weights enter through their square roots.

Source code in src/torch_flash/fitting.py

def least_squares_loss(
    prediction: Tensor,
    observation: Tensor,
    *,
    scale: Tensor | float = 1.0,
    weights: Tensor | None = None,
) -> Tensor:
    """Return a dimensionless weighted mean-square residual.

    Parameters
    ----------
    prediction, observation
        Broadcast-compatible model predictions and observations in matching
        physical units.
    scale
        Positive residual scale in the same units. Scalar or broadcastable
        tensor.
    weights
        Optional nonnegative statistical weights, broadcastable to the
        residual shape.

    Returns
    -------
    Tensor
        Scalar mean of the squared, scaled residuals; weights enter through
        their square roots.
    """
    residual = (prediction - observation) / scale
    if weights is not None:
        residual = residual * torch.sqrt(weights)
    return torch.mean(residual.square())

phase_equilibrium_residual

phase_equilibrium_residual(
    model: StateModel,
    temperature: Tensor,
    pressure: Tensor,
    phase1_composition: Tensor,
    phase2_composition: Tensor,
    *,
    phase_kinds: tuple[PhaseKind, PhaseKind] = (
        "liquid",
        "vapor",
    ),
) -> Tensor

Return component log-fugacity equalities for measured phase pairs.

Parameters:

Name	Type	Description	Default
model	StateModel	Homogeneous-state fugacity model.	required
temperature	Tensor	Scalar or batched temperatures in K and pressures in Pa.	required
pressure	Tensor	Scalar or batched temperatures in K and pressures in Pa.	required
phase1_composition	Tensor	Equally shaped, strictly positive phase mole fractions.	required
phase2_composition	Tensor	Equally shaped, strictly positive phase mole fractions.	required
phase_kinds	tuple[PhaseKind, PhaseKind]	Root request corresponding to each measured phase.	('liquid', 'vapor')

Returns:

Type	Description
Tensor	Dimensionless component residual ln(x1_i phi1_i) - ln(x2_i phi2_i).

Raises:

Type	Description
ValueError	If phase compositions differ in shape or are nonpositive/nonfinite.

Notes

Inputs may contain independent batched states. The residual is dimensionless and is zero when every component has equal fugacity in the two requested phases. Strictly positive compositions are required because the thermodynamic equality is evaluated in logarithmic form.

Source code in src/torch_flash/fitting.py

def phase_equilibrium_residual(
    model: StateModel,
    temperature: Tensor,
    pressure: Tensor,
    phase1_composition: Tensor,
    phase2_composition: Tensor,
    *,
    phase_kinds: tuple[PhaseKind, PhaseKind] = ("liquid", "vapor"),
) -> Tensor:
    """Return component log-fugacity equalities for measured phase pairs.

    Parameters
    ----------
    model
        Homogeneous-state fugacity model.
    temperature, pressure
        Scalar or batched temperatures in K and pressures in Pa.
    phase1_composition, phase2_composition
        Equally shaped, strictly positive phase mole fractions.
    phase_kinds
        Root request corresponding to each measured phase.

    Returns
    -------
    Tensor
        Dimensionless component residual
        ``ln(x1_i phi1_i) - ln(x2_i phi2_i)``.

    Raises
    ------
    ValueError
        If phase compositions differ in shape or are nonpositive/nonfinite.

    Notes
    -----
    Inputs may contain independent batched states.  The residual is
    dimensionless and is zero when every component has equal fugacity in the
    two requested phases.  Strictly positive compositions are required
    because the thermodynamic equality is evaluated in logarithmic form.
    """
    phase1 = normalize_composition(phase1_composition)
    phase2 = normalize_composition(phase2_composition)
    if phase1.shape != phase2.shape:
        raise ValueError("phase-equilibrium compositions must have equal shapes")
    if not bool(
        torch.isfinite(phase1).all()
        & torch.isfinite(phase2).all()
        & (phase1 > 0.0).all()
        & (phase2 > 0.0).all()
    ):
        raise ValueError("phase-equilibrium compositions must be finite and positive")
    return (
        torch.log(phase1)
        + model.log_fugacity_coefficients(
            temperature,
            pressure,
            phase1,
            phase_kinds[0],
        )
        - torch.log(phase2)
        - model.log_fugacity_coefficients(
            temperature,
            pressure,
            phase2,
            phase_kinds[1],
        )
    )

fit_parameters

fit_parameters(
    parameters: Iterable[Parameter],
    closure: Callable[[], Tensor],
    *,
    optimizer: OptimizerKind = "adam",
    learning_rate: float = 0.05,
    weight_decay: float = 0.0,
    momentum: float = 0.0,
    max_iterations: int = 500,
    tolerance: float = 1e-10,
    no_improvement_patience: int | None = None,
) -> FitResult

Fit arbitrary PyTorch thermodynamic parameters with a PyTorch optimizer.

Parameters:

Name	Type	Description	Default
parameters	Iterable[Parameter]	Trainable tensors passed to the selected optimizer.	required
closure	Callable[[], Tensor]	Zero-argument function that recomputes one finite differentiable scalar loss.	required
optimizer	OptimizerKind	"adam", "adamw", full-batch "sgd", or limited-memory "lbfgs" with a strong-Wolfe line search.	'adam'
learning_rate	float	Optimizer learning rate.	0.05
weight_decay	float	Nonnegative optimizer weight decay. A nonzero value is a regularizer, not a numerical acceleration, and requires a defensible prior.	0.0
momentum	float	SGD momentum in [0, 1). Ignored by Adam and AdamW.	0.0
max_iterations	int	Maximum optimizer steps.	500
tolerance	float	Relative successive-loss stopping threshold.	1e-10
no_improvement_patience	int | None	Optional number of consecutive evaluated losses without a new minimum before stopping. The returned converged flag remains False for this performance termination.	None

Returns:

Type	Description
FitResult	Initial and accepted-iterate loss history plus explicit stopping diagnostics. The supplied tensors are restored to the minimum-loss accepted point.

Raises:

Type	Description
ValueError	If no parameters are supplied or the closure returns a nonfinite or nonscalar loss.

Notes

The closure must recompute and return a scalar differentiable loss. Bounds can be imposed by parameterizing physical values through sigmoid/softplus transforms in the model.

Source code in src/torch_flash/fitting.py

def fit_parameters(
    parameters: Iterable[nn.Parameter],
    closure: Callable[[], Tensor],
    *,
    optimizer: OptimizerKind = "adam",
    learning_rate: float = 0.05,
    weight_decay: float = 0.0,
    momentum: float = 0.0,
    max_iterations: int = 500,
    tolerance: float = 1.0e-10,
    no_improvement_patience: int | None = None,
) -> FitResult:
    """Fit arbitrary PyTorch thermodynamic parameters with a PyTorch optimizer.

    Parameters
    ----------
    parameters
        Trainable tensors passed to the selected optimizer.
    closure
        Zero-argument function that recomputes one finite differentiable
        scalar loss.
    optimizer
        ``"adam"``, ``"adamw"``, full-batch ``"sgd"``, or limited-memory
        ``"lbfgs"`` with a strong-Wolfe line search.
    learning_rate
        Optimizer learning rate.
    weight_decay
        Nonnegative optimizer weight decay. A nonzero value is a regularizer,
        not a numerical acceleration, and requires a defensible prior.
    momentum
        SGD momentum in ``[0, 1)``. Ignored by Adam and AdamW.
    max_iterations
        Maximum optimizer steps.
    tolerance
        Relative successive-loss stopping threshold.
    no_improvement_patience
        Optional number of consecutive evaluated losses without a new
        minimum before stopping. The returned ``converged`` flag remains
        ``False`` for this performance termination.

    Returns
    -------
    FitResult
        Initial and accepted-iterate loss history plus explicit stopping
        diagnostics. The supplied tensors are restored to the minimum-loss
        accepted point.

    Raises
    ------
    ValueError
        If no parameters are supplied or the closure returns a nonfinite or
        nonscalar loss.

    Notes
    -----
    The closure must recompute and return a scalar differentiable loss. Bounds
    can be imposed by parameterizing physical values through sigmoid/softplus
    transforms in the model.
    """
    trainable = tuple(parameters)
    if not trainable:
        raise ValueError("at least one trainable parameter is required")
    if (
        optimizer not in ("adam", "adamw", "sgd", "lbfgs")
        or learning_rate <= 0.0
        or weight_decay < 0.0
        or not 0.0 <= momentum < 1.0
        or max_iterations <= 0
        or tolerance < 0.0
    ):
        raise ValueError("fitting controls must be nonnegative with positive rate and iterations")
    if no_improvement_patience is not None and no_improvement_patience <= 0:
        raise ValueError("no-improvement patience must be positive when provided")
    if optimizer == "lbfgs" and (weight_decay != 0.0 or momentum != 0.0):
        raise ValueError("LBFGS does not accept weight decay or momentum")
    torch_optimizer: torch.optim.Optimizer
    if optimizer == "adam":
        torch_optimizer = torch.optim.Adam(
            trainable,
            lr=learning_rate,
            weight_decay=weight_decay,
        )
    elif optimizer == "adamw":
        torch_optimizer = torch.optim.AdamW(
            trainable,
            lr=learning_rate,
            weight_decay=weight_decay,
        )
    elif optimizer == "sgd":
        torch_optimizer = torch.optim.SGD(
            trainable,
            lr=learning_rate,
            momentum=momentum,
            weight_decay=weight_decay,
        )
    else:
        torch_optimizer = torch.optim.LBFGS(
            trainable,
            lr=learning_rate,
            max_iter=1,
            max_eval=10,
            history_size=10,
            line_search_fn="strong_wolfe",
        )
    initial_loss = closure()
    if initial_loss.ndim != 0 or not bool(torch.isfinite(initial_loss)):
        raise ValueError("fitting closure must return one finite scalar loss")
    initial = float(initial_loss.detach())
    history = [initial]
    converged = False
    previous = initial
    best = initial
    best_parameters = tuple(parameter.detach().clone() for parameter in trainable)
    evaluations_without_improvement = 0
    stopping_reason: Literal["tolerance", "patience", "iteration-limit"] = "iteration-limit"
    for _iteration in range(1, max_iterations + 1):
        if optimizer == "lbfgs":

            def line_search_closure() -> Tensor:
                torch_optimizer.zero_grad()
                trial_loss = closure()
                if trial_loss.ndim != 0 or not bool(torch.isfinite(trial_loss)):
                    raise ValueError("fitting closure must return one finite scalar loss")
                trial_loss.backward()
                return trial_loss

            cast(torch.optim.LBFGS, torch_optimizer).step(line_search_closure)
            loss = closure()
        else:
            torch_optimizer.zero_grad()
            loss = closure()
            if loss.ndim != 0 or not bool(torch.isfinite(loss)):
                raise ValueError("fitting closure must return one finite scalar loss")
            loss.backward()
            cast(torch.optim.Optimizer, torch_optimizer).step()
            loss = closure()
        if loss.ndim != 0 or not bool(torch.isfinite(loss)):
            raise ValueError("fitting closure must return one finite scalar loss")
        current = float(loss.detach())
        history.append(current)
        if current < best:
            best = current
            best_parameters = tuple(parameter.detach().clone() for parameter in trainable)
            evaluations_without_improvement = 0
        else:
            evaluations_without_improvement += 1
        if abs(float(previous) - current) <= tolerance * max(1.0, abs(current)):
            converged = True
            stopping_reason = "tolerance"
            break
        if (
            no_improvement_patience is not None
            and evaluations_without_improvement >= no_improvement_patience
        ):
            stopping_reason = "patience"
            break
        previous = current
    with torch.no_grad():
        for parameter, selected in zip(trainable, best_parameters, strict=True):
            parameter.copy_(selected)
    return FitResult(
        tuple(history),
        converged,
        _iteration,
        best,
        stopping_reason,
    )