libeq

SolverData

Bases: BaseModel

Source code in src/libeq/data_structure.py

class SolverData(BaseModel):
    model_config = ConfigDict(extra="forbid")

    distribution_opts: DistributionParameters = DistributionParameters()
    titration_opts: SimulationTitrationParameters = SimulationTitrationParameters()
    potentiometry_options: PotentiometryOptions = PotentiometryOptions()

    components: List[str]
    stoichiometry: Np2DArrayInt8
    solid_stoichiometry: Np2DArrayInt8
    log_beta: Np1DArrayFp64
    log_beta_sigma: Np1DArrayFp64 = np.array([])
    log_beta_ref_dbh: Np2DArrayFp64 = np.empty((0, 2))
    log_ks: Np1DArrayFp64 = np.array([])
    log_ks_sigma: Np1DArrayFp64 = np.array([])
    log_ks_ref_dbh: Np2DArrayFp64 = np.empty((0, 2))

    charges: Np1DArrayFp64 = np.array([])

    ionic_strength_dependence: bool = False
    reference_ionic_str_species: Np1DArrayFp64 | float = 0
    reference_ionic_str_solids: Np1DArrayFp64 | float = 0
    dbh_params: Np1DArrayFp64 = np.zeros(8)

    @computed_field
    @cached_property
    def species_charges(self) -> Np1DArrayFp64:
        return (self.stoichiometry * self.charges[:, np.newaxis]).sum(axis=0)

    @computed_field
    @cached_property
    def solid_charges(self) -> Np1DArrayFp64:
        return (self.solid_stoichiometry * self.charges[:, np.newaxis]).sum(axis=0)

    @computed_field(repr=False)
    @cached_property
    def z_star_species(self) -> Np1DArrayFp64:
        return (self.stoichiometry * (self.charges[:, np.newaxis] ** 2)).sum(
            axis=0
        ) - self.species_charges**2

    @computed_field(repr=False)
    @cached_property
    def p_star_species(self) -> Np1DArrayFp64:
        return self.stoichiometry.sum(axis=0) - 1

    @computed_field(repr=False)
    @cached_property
    def z_star_solids(self) -> Np1DArrayFp64:
        return (self.solid_stoichiometry * (self.charges[:, np.newaxis] ** 2)).sum(
            axis=0
        ) - self.solid_charges**2

    @computed_field(repr=False)
    @cached_property
    def p_star_solids(self) -> Np1DArrayFp64:
        return self.solid_stoichiometry.sum(axis=0)

    @computed_field
    @cached_property
    def dbh_values(self) -> Dict[str, Np1DArrayFp64]:
        result = dict()
        for phase, iref, per_species_cde, z, p in zip(
            ("species", "solids"),
            (self.reference_ionic_str_species, self.reference_ionic_str_solids),
            (self.log_beta_ref_dbh, self.log_ks_ref_dbh),
            (self.z_star_species, self.z_star_solids),
            (self.p_star_species, self.p_star_solids),
        ):
            dbh_values = dict()
            dbh_values["azast"] = self.dbh_params[0] * z
            dbh_values["adh"] = self.dbh_params[0]
            dbh_values["bdh"] = self.dbh_params[1]
            dbh_values["cdh"] = self.dbh_params[2] * p + self.dbh_params[3] * z
            dbh_values["ddh"] = self.dbh_params[4] * p + self.dbh_params[5] * z
            dbh_values["edh"] = self.dbh_params[6] * p + self.dbh_params[7] * z
            dbh_values["fib"] = np.sqrt(iref) / (1 + self.dbh_params[1] * np.sqrt(iref))

            not_zero_columns = np.where(np.any(per_species_cde != 0, axis=0))[0]
            for i in not_zero_columns:
                dbh_values["cdh"][i] = per_species_cde[0][i]
                dbh_values["ddh"][i] = per_species_cde[1][i]
                dbh_values["edh"][i] = per_species_cde[2][i]

            result[phase] = dbh_values
        return result

    @computed_field
    @cached_property
    def species_names(self) -> List[str]:
        return self.components + _assemble_species_names(
            self.components, self.stoichiometry
        )

    @computed_field
    @cached_property
    def solids_names(self) -> List[str]:
        return _assemble_species_names(self.components, self.solid_stoichiometry)

    @computed_field
    @cached_property
    def nc(self) -> int:
        return self.stoichiometry.shape[0]

    @computed_field
    @cached_property
    def ns(self) -> int:
        return self.stoichiometry.shape[1]

    @computed_field
    @cached_property
    def nf(self) -> int:
        return self.solid_stoichiometry.shape[1]

    @classmethod
    def load_from_bstac(cls, file_path: str) -> "SolverData":
        data = dict()
        with open(file_path, "r") as file:
            lines = file.readlines()
        parsed_data = parse_BSTAC_file(lines)

        temperature = parsed_data["TEMP"]
        data["stoichiometry"] = np.array(
            [
                [d[key] for key in d if key.startswith("IX")]
                for d in parsed_data["species"]
            ]
        ).T
        data["solid_stoichiometry"] = np.empty(
            (data["stoichiometry"].shape[0], 0), dtype=np.int8
        )
        data["log_beta"] = np.array([d["BLOG"] for d in parsed_data["species"]])

        data["charges"] = np.array(parsed_data.get("charges", []))
        data["components"] = parsed_data["comp_name"]
        data["ionic_strength_dependence"] = parsed_data["ICD"] != 0
        if data["ionic_strength_dependence"]:
            data["reference_ionic_str_species"] = np.array(
                [parsed_data["IREF"] for _ in range(data["stoichiometry"].shape[1])]
            )
            data["reference_ionic_str_solids"] = np.array(
                [
                    parsed_data["IREF"]
                    for _ in range(data["solid_stoichiometry"].shape[1])
                ]
            )
            data["dbh_params"] = [
                parsed_data[i] for i in ["AT", "BT", "c0", "c1", "d0", "d1", "e0", "e1"]
            ]

        titration_options = [
            PotentiometryTitrationsParameters(
                c0=np.array([c["C0"] for c in t["components_concentrations"]]),
                ct=np.array([c["CTT"] for c in t["components_concentrations"]]),
                electro_active_compoment=(
                    t["titration_comp_settings"][1]
                    if t["titration_comp_settings"][1] != 0
                    else len(data["components"]) - 1
                ),
                e0=t["potential_params"][0],
                e0_sigma=t["potential_params"][1],
                slope=(
                    t["potential_params"][4]
                    if t["potential_params"][4] != 0
                    else (temperature + 273.15) / 11.6048 * 2.303
                ),
                v0=t["v_params"][0],
                v0_sigma=t["v_params"][1],
                v_add=np.array(t["volume"]),
                emf=np.array(t["potential"]),
            )
            for t in parsed_data["titrations"]
        ]

        weights = parsed_data.get("MODE", 1)
        match weights:
            case 0:
                weights = "calculated"
            case 1:
                weights = "constants"
            case 2:
                weights = "given"
            case _:
                raise ValueError("Invalid MODE value")

        data["potentiometry_options"] = PotentiometryOptions(
            titrations=titration_options,
            weights=weights,
            px_range=[parsed_data["PHI"], parsed_data["PHF"]],
            beta_flags=[s["KEY"] for s in parsed_data["species"]],
            conc_flags=[],
            pot_flags=[],
        )
        return cls(**data)

    @classmethod
    def load_from_pyes(cls, pyes_data: str | dict) -> "SolverData":
        if isinstance(pyes_data, str):
            with open(pyes_data, "r") as file:
                pyes_data = json.load(file)
        data = dict()
        data["components"] = list(pyes_data["compModel"]["Name"].values())

        data["stoichiometry"] = np.row_stack(
            [
                list(pyes_data["speciesModel"][col].values())
                for col in data["components"]
            ]
        )
        data["log_beta"] = np.array(list(pyes_data["speciesModel"]["LogB"].values()))
        data["log_beta_sigma"] = np.array(
            list(pyes_data["speciesModel"]["Sigma"].values())
        )
        data["log_beta_ref_dbh"] = np.vstack(
            (
                list(pyes_data["speciesModel"]["CGF"].values()),
                list(pyes_data["speciesModel"]["DGF"].values()),
                list(pyes_data["speciesModel"]["EGF"].values()),
            )
        )

        data["solid_stoichiometry"] = np.row_stack(
            [
                list(pyes_data["solidSpeciesModel"][col].values())
                for col in data["components"]
            ]
        )
        data["log_ks"] = np.array(
            list(pyes_data["solidSpeciesModel"]["LogKs"].values())
        )
        data["log_ks_sigma"] = np.array(
            list(pyes_data["solidSpeciesModel"]["Sigma"].values())
        )
        data["log_ks_ref_dbh"] = np.vstack(
            (
                list(pyes_data["solidSpeciesModel"]["CGF"].values()),
                list(pyes_data["solidSpeciesModel"]["DGF"].values()),
                list(pyes_data["solidSpeciesModel"]["EGF"].values()),
            )
        )

        data["charges"] = np.array(list(pyes_data["compModel"]["Charge"].values()))
        data["ionic_strength_dependence"] = pyes_data["imode"] != 0
        data["reference_ionic_str_species"] = np.array(
            list(pyes_data["speciesModel"]["Ref. Ionic Str."].values())
        )
        data["reference_ionic_str_species"] = np.where(
            data["reference_ionic_str_species"] == 0,
            pyes_data["ris"],
            data["reference_ionic_str_species"],
        )

        data["reference_ionic_str_solids"] = np.array(
            list(pyes_data["solidSpeciesModel"]["Ref. Ionic Str."].values())
        )
        data["reference_ionic_str_solids"] = np.where(
            data["reference_ionic_str_solids"] == 0,
            pyes_data["ris"],
            data["reference_ionic_str_solids"],
        )

        data["dbh_params"] = [
            pyes_data[name] for name in ["a", "b", "c0", "c1", "d0", "d1", "e0", "e1"]
        ]

        data["distribution_opts"] = DistributionParameters(
            c0=np.array(list(pyes_data["concModel"]["C0"].values())),
            c0_sigma=np.array(list(pyes_data["concModel"]["Sigma C0"].values())),
            initial_log=pyes_data.get("initialLog"),
            final_log=pyes_data.get("finalLog"),
            log_increments=pyes_data.get("logInc"),
            independent_component=pyes_data.get("ind_comp"),
        )

        data["titration_opts"] = SimulationTitrationParameters(
            c0=np.array(list(pyes_data["concModel"]["C0"].values())),
            c0_sigma=np.array(list(pyes_data["concModel"]["Sigma C0"].values())),
            ct=np.array(list(pyes_data["concModel"]["CT"].values())),
            ct_sigma=np.array(list(pyes_data["concModel"]["Sigma CT"].values())),
            v0=pyes_data.get("v0"),
            v_increment=pyes_data.get("vinc"),
            n_add=pyes_data.get("nop"),
        )
        return cls(**data)

model_config class-attribute instance-attribute

model_config = ConfigDict(extra='forbid')

distribution_opts class-attribute instance-attribute

distribution_opts: DistributionParameters = (
    DistributionParameters()
)

titration_opts class-attribute instance-attribute

titration_opts: SimulationTitrationParameters = (
    SimulationTitrationParameters()
)

potentiometry_options class-attribute instance-attribute

potentiometry_options: PotentiometryOptions = (
    PotentiometryOptions()
)

components instance-attribute

components: List[str]

stoichiometry instance-attribute

stoichiometry: Np2DArrayInt8

solid_stoichiometry instance-attribute

solid_stoichiometry: Np2DArrayInt8

log_beta instance-attribute

log_beta: Np1DArrayFp64

log_beta_sigma class-attribute instance-attribute

log_beta_sigma: Np1DArrayFp64 = array([])

log_beta_ref_dbh class-attribute instance-attribute

log_beta_ref_dbh: Np2DArrayFp64 = empty((0, 2))

log_ks class-attribute instance-attribute

log_ks: Np1DArrayFp64 = array([])

log_ks_sigma class-attribute instance-attribute

log_ks_sigma: Np1DArrayFp64 = array([])

log_ks_ref_dbh class-attribute instance-attribute

log_ks_ref_dbh: Np2DArrayFp64 = empty((0, 2))

charges class-attribute instance-attribute

charges: Np1DArrayFp64 = array([])

ionic_strength_dependence class-attribute instance-attribute

ionic_strength_dependence: bool = False

reference_ionic_str_species class-attribute instance-attribute

reference_ionic_str_species: Np1DArrayFp64 | float = 0

reference_ionic_str_solids class-attribute instance-attribute

reference_ionic_str_solids: Np1DArrayFp64 | float = 0

dbh_params class-attribute instance-attribute

dbh_params: Np1DArrayFp64 = zeros(8)

species_charges cached property

species_charges: Np1DArrayFp64

solid_charges cached property

solid_charges: Np1DArrayFp64

z_star_species cached property

z_star_species: Np1DArrayFp64

p_star_species cached property

p_star_species: Np1DArrayFp64

z_star_solids cached property

z_star_solids: Np1DArrayFp64

p_star_solids cached property

p_star_solids: Np1DArrayFp64

dbh_values cached property

dbh_values: Dict[str, Np1DArrayFp64]

species_names cached property

species_names: List[str]

solids_names cached property

solids_names: List[str]

nc cached property

nc: int

ns cached property

ns: int

nf cached property

nf: int

load_from_bstac classmethod

load_from_bstac(file_path: str) -> SolverData

Source code in src/libeq/data_structure.py

@classmethod
def load_from_bstac(cls, file_path: str) -> "SolverData":
    data = dict()
    with open(file_path, "r") as file:
        lines = file.readlines()
    parsed_data = parse_BSTAC_file(lines)

    temperature = parsed_data["TEMP"]
    data["stoichiometry"] = np.array(
        [
            [d[key] for key in d if key.startswith("IX")]
            for d in parsed_data["species"]
        ]
    ).T
    data["solid_stoichiometry"] = np.empty(
        (data["stoichiometry"].shape[0], 0), dtype=np.int8
    )
    data["log_beta"] = np.array([d["BLOG"] for d in parsed_data["species"]])

    data["charges"] = np.array(parsed_data.get("charges", []))
    data["components"] = parsed_data["comp_name"]
    data["ionic_strength_dependence"] = parsed_data["ICD"] != 0
    if data["ionic_strength_dependence"]:
        data["reference_ionic_str_species"] = np.array(
            [parsed_data["IREF"] for _ in range(data["stoichiometry"].shape[1])]
        )
        data["reference_ionic_str_solids"] = np.array(
            [
                parsed_data["IREF"]
                for _ in range(data["solid_stoichiometry"].shape[1])
            ]
        )
        data["dbh_params"] = [
            parsed_data[i] for i in ["AT", "BT", "c0", "c1", "d0", "d1", "e0", "e1"]
        ]

    titration_options = [
        PotentiometryTitrationsParameters(
            c0=np.array([c["C0"] for c in t["components_concentrations"]]),
            ct=np.array([c["CTT"] for c in t["components_concentrations"]]),
            electro_active_compoment=(
                t["titration_comp_settings"][1]
                if t["titration_comp_settings"][1] != 0
                else len(data["components"]) - 1
            ),
            e0=t["potential_params"][0],
            e0_sigma=t["potential_params"][1],
            slope=(
                t["potential_params"][4]
                if t["potential_params"][4] != 0
                else (temperature + 273.15) / 11.6048 * 2.303
            ),
            v0=t["v_params"][0],
            v0_sigma=t["v_params"][1],
            v_add=np.array(t["volume"]),
            emf=np.array(t["potential"]),
        )
        for t in parsed_data["titrations"]
    ]

    weights = parsed_data.get("MODE", 1)
    match weights:
        case 0:
            weights = "calculated"
        case 1:
            weights = "constants"
        case 2:
            weights = "given"
        case _:
            raise ValueError("Invalid MODE value")

    data["potentiometry_options"] = PotentiometryOptions(
        titrations=titration_options,
        weights=weights,
        px_range=[parsed_data["PHI"], parsed_data["PHF"]],
        beta_flags=[s["KEY"] for s in parsed_data["species"]],
        conc_flags=[],
        pot_flags=[],
    )
    return cls(**data)

load_from_pyes classmethod

load_from_pyes(pyes_data: str | dict) -> SolverData

Source code in src/libeq/data_structure.py

@classmethod
def load_from_pyes(cls, pyes_data: str | dict) -> "SolverData":
    if isinstance(pyes_data, str):
        with open(pyes_data, "r") as file:
            pyes_data = json.load(file)
    data = dict()
    data["components"] = list(pyes_data["compModel"]["Name"].values())

    data["stoichiometry"] = np.row_stack(
        [
            list(pyes_data["speciesModel"][col].values())
            for col in data["components"]
        ]
    )
    data["log_beta"] = np.array(list(pyes_data["speciesModel"]["LogB"].values()))
    data["log_beta_sigma"] = np.array(
        list(pyes_data["speciesModel"]["Sigma"].values())
    )
    data["log_beta_ref_dbh"] = np.vstack(
        (
            list(pyes_data["speciesModel"]["CGF"].values()),
            list(pyes_data["speciesModel"]["DGF"].values()),
            list(pyes_data["speciesModel"]["EGF"].values()),
        )
    )

    data["solid_stoichiometry"] = np.row_stack(
        [
            list(pyes_data["solidSpeciesModel"][col].values())
            for col in data["components"]
        ]
    )
    data["log_ks"] = np.array(
        list(pyes_data["solidSpeciesModel"]["LogKs"].values())
    )
    data["log_ks_sigma"] = np.array(
        list(pyes_data["solidSpeciesModel"]["Sigma"].values())
    )
    data["log_ks_ref_dbh"] = np.vstack(
        (
            list(pyes_data["solidSpeciesModel"]["CGF"].values()),
            list(pyes_data["solidSpeciesModel"]["DGF"].values()),
            list(pyes_data["solidSpeciesModel"]["EGF"].values()),
        )
    )

    data["charges"] = np.array(list(pyes_data["compModel"]["Charge"].values()))
    data["ionic_strength_dependence"] = pyes_data["imode"] != 0
    data["reference_ionic_str_species"] = np.array(
        list(pyes_data["speciesModel"]["Ref. Ionic Str."].values())
    )
    data["reference_ionic_str_species"] = np.where(
        data["reference_ionic_str_species"] == 0,
        pyes_data["ris"],
        data["reference_ionic_str_species"],
    )

    data["reference_ionic_str_solids"] = np.array(
        list(pyes_data["solidSpeciesModel"]["Ref. Ionic Str."].values())
    )
    data["reference_ionic_str_solids"] = np.where(
        data["reference_ionic_str_solids"] == 0,
        pyes_data["ris"],
        data["reference_ionic_str_solids"],
    )

    data["dbh_params"] = [
        pyes_data[name] for name in ["a", "b", "c0", "c1", "d0", "d1", "e0", "e1"]
    ]

    data["distribution_opts"] = DistributionParameters(
        c0=np.array(list(pyes_data["concModel"]["C0"].values())),
        c0_sigma=np.array(list(pyes_data["concModel"]["Sigma C0"].values())),
        initial_log=pyes_data.get("initialLog"),
        final_log=pyes_data.get("finalLog"),
        log_increments=pyes_data.get("logInc"),
        independent_component=pyes_data.get("ind_comp"),
    )

    data["titration_opts"] = SimulationTitrationParameters(
        c0=np.array(list(pyes_data["concModel"]["C0"].values())),
        c0_sigma=np.array(list(pyes_data["concModel"]["Sigma C0"].values())),
        ct=np.array(list(pyes_data["concModel"]["CT"].values())),
        ct_sigma=np.array(list(pyes_data["concModel"]["Sigma CT"].values())),
        v0=pyes_data.get("v0"),
        v_increment=pyes_data.get("vinc"),
        n_add=pyes_data.get("nop"),
    )
    return cls(**data)

species_concentration

species_concentration(
    concentration, log_beta, stoichiometry, full=False
)

Calculate the species concentrations through the mass action law.

\[ S_{i} = \beta_i \prod_{j=1}^{N_c} C_j^{p_{ij}} \]

With \(S_i\) being the concentration of the species \(i\), \(\beta_i\) the equilibrium constant of the species \(i\), \(C_j\) the concentration of the component \(j\), and \(p_{ij}\) the stoichiometric coefficient of the component \(j\) in the species \(i\).

Parameters:

Name	Type	Description	Default
concentration	ndarray	The concentration array of shape (n, c+p), where n is the number of points c is the number of components and p is the number of solid species.	required
log_beta	ndarray	The logarithm of the equilibrium constants with shape (n, s), where s is the number of solid species.	required
stoichiometry	ndarray	The stoichiometry matrix with shape (n, s), where s is the number of soluble species.	required
full	bool	If True, return the concentrations of all species including the original concentrations. If False, return only the concentrations of the new species.	False

Returns:

Type	Description
ndarray	The calculated species concentrations.

Source code in src/libeq/utils.py

def species_concentration(
    concentration,
    log_beta,
    stoichiometry,
    full=False,
):
    r"""
    Calculate the species concentrations through the mass action law.

    $$
    S_{i} = \beta_i \prod_{j=1}^{N_c} C_j^{p_{ij}}
    $$

    With $S_i$ being the concentration of the species $i$, $\beta_i$ the equilibrium constant of the species $i$,
    $C_j$ the concentration of the component $j$, and $p_{ij}$ the stoichiometric coefficient of the component $j$ in the species $i$.

    Parameters
    ----------
    concentration : numpy.ndarray
        The concentration array of shape (n, c+p), where n is the number of points c is the number of components and p is the number of solid species.
    log_beta : numpy.ndarray
        The logarithm of the equilibrium constants with shape (n, s), where s is the number of solid species.
    stoichiometry : numpy.ndarray
        The stoichiometry matrix with shape (n, s), where s is the number of soluble species.
    full : bool, optional
        If True, return the concentrations of all species including the original concentrations.
        If False, return only the concentrations of the new species.

    Returns
    -------
    numpy.ndarray
        The calculated species concentrations.

    """
    nc = stoichiometry.shape[0]
    _c = np.log10(concentration[:, :nc])

    cext = 10 ** (log_beta + _c @ stoichiometry)

    if full:
        p = np.concatenate((concentration, cext), axis=1)
    else:
        p = cext

    return p

uncertanties

uncertanties(
    concentrations,
    stoichiometry,
    solid_stoichiometry,
    log_b,
    log_ks,
    log_beta_sigma,
    log_ks_sigma,
    conc_sigma,
    indepenent_comp: int | None = None,
)

Calculate the uncertainties for the components and species given the input.

Args: concentrations (ndarray): Array of shape (n, m) representing the concentrations of the components and species. stoichiometry (ndarray): Array of shape (n, p) representing the stoichiometric coefficients of the components and species. solid_stoichiometry (ndarray): Array of shape (p, q) representing the stoichiometric coefficients of the solid species. log_b (ndarray): Array of shape (n,) representing the logarithm of the beta values. log_ks (ndarray): Array of shape (n,) representing the logarithm of the equilibrium constants. beta_sigma (float): Standard deviation of the beta values. ks_sigma (float): Standard deviation of the equilibrium constants. conc_sigma (ndarray): Array of shape (n,) representing the standard deviation of the concentrations. indepenent_comp (int | None, optional): Index of the independent component. Defaults to None.

Returns: tuple: A tuple containing the uncertainties for the species and solid species.

Source code in src/libeq/errors.py

def uncertanties(
    concentrations,
    stoichiometry,
    solid_stoichiometry,
    log_b,
    log_ks,
    log_beta_sigma,
    log_ks_sigma,
    conc_sigma,
    indepenent_comp: int | None = None,
):
    """
    Calculate the uncertainties for the components and species given the input.

    Args:
        concentrations (ndarray): Array of shape (n, m) representing the concentrations of the components and species.
        stoichiometry (ndarray): Array of shape (n, p) representing the stoichiometric coefficients of the components and species.
        solid_stoichiometry (ndarray): Array of shape (p, q) representing the stoichiometric coefficients of the solid species.
        log_b (ndarray): Array of shape (n,) representing the logarithm of the beta values.
        log_ks (ndarray): Array of shape (n,) representing the logarithm of the equilibrium constants.
        beta_sigma (float): Standard deviation of the beta values.
        ks_sigma (float): Standard deviation of the equilibrium constants.
        conc_sigma (ndarray): Array of shape (n,) representing the standard deviation of the concentrations.
        indepenent_comp (int | None, optional): Index of the independent component. Defaults to None.

    Returns:
        tuple: A tuple containing the uncertainties for the species and solid species.
    """
    # Get betas from log betas
    nc = stoichiometry.shape[0]
    ns = stoichiometry.shape[1]
    nf = solid_stoichiometry.shape[1]
    num_points = concentrations.shape[0]
    species_sigma_result = np.zeros((num_points, nc + ns))
    solid_sigma_result = np.zeros((num_points, nf))
    distribution = indepenent_comp is not None
    all_beta_sigma = log_beta_sigma * np.log(10) * (10**log_b)
    all_ks_sigma = log_ks_sigma * np.log(10) * (10**log_ks)

    for point in range(num_points):
        beta = 10 ** log_b[point]
        ks = 10 ** log_ks[point]
        beta_sigma = all_beta_sigma[point]
        ks_sigma = all_ks_sigma[point]

        c_free = concentrations[point, :nc]
        c_spec = concentrations[point, nc + nf :]
        c_solid = concentrations[point, nc : nc + nf]
        # c_free, c_spec, c_solid = np.atleast_2d(c_free, c_spec, c_solid)

        saturation_index = _compute_saturation_index(
            c_free[np.newaxis, :], log_ks[point], solid_stoichiometry
        )[0]

        with_solids = any(c_solid > 0)

        to_skip = np.concatenate(([False for _ in range(nc)], c_solid == 0))
        if with_solids:
            nt = nc + nf
        else:
            nt = nc

        # Define dimension of arrays required
        M = np.zeros(shape=(nt, nt))

        der_free_beta = np.zeros(shape=(nc, ns))
        der_free_tot = np.zeros(shape=(nc, nc))
        der_free_ks = np.zeros(shape=(nc, nf))

        der_solid_beta = np.zeros(shape=(nf, ns))
        der_solid_tot = np.zeros(shape=(nf, nc))
        der_solid_ks = np.zeros(shape=(nf, nf))

        b = -stoichiometry * (c_spec / beta)
        d = np.identity(nc)
        f = np.zeros(shape=(nc, nf))

        # Compute common matrix term
        M[:nc, :nc] = (
            (
                np.tile(c_spec, (nc, nc, 1))
                / np.tile(c_free.reshape((nc, 1)), (nc, 1, ns))
            )
            * np.tile(stoichiometry, (nc, 1, 1))
            * np.rot90(np.tile(stoichiometry, (nc, 1, 1)), -1, axes=(0, 1))
        ).sum(axis=-1)
        M[:nc, :nc] += d

        if with_solids:
            M[:nc, nc:nt] = solid_stoichiometry
            M[nc:nt, :nc] = solid_stoichiometry.T * (
                np.tile(saturation_index, (nc, 1)).T / np.tile(c_free, (nt - nc, 1))
            )

            f = np.concatenate((f, np.diag(saturation_index / ks)), axis=0)
            b = np.concatenate(
                (
                    b,
                    [[0 for _ in range(ns)] for _ in range(nf)],
                )
            )
            d = np.concatenate(
                (
                    d,
                    [[0 for _ in range(nc)] for _ in range(nf)],
                )
            )

            der_solid_beta = np.delete(der_solid_beta, c_solid == 0, axis=0)
            der_solid_tot = np.delete(der_solid_tot, c_solid == 0, axis=0)
            der_solid_ks = np.delete(der_solid_ks, c_solid == 0, axis=0)

            M = np.delete(M, to_skip, axis=0)
            M = np.delete(M, to_skip, axis=1)

            b = np.delete(b, to_skip, axis=0)

            d = np.delete(d, to_skip, axis=0)

            f = np.delete(f, to_skip, axis=0)

        if indepenent_comp is not None:
            M = np.delete(M, indepenent_comp, axis=0)
            M = np.delete(M, indepenent_comp, axis=1)

            b = np.delete(b, indepenent_comp, axis=0)
            d = np.delete(d, indepenent_comp, axis=0)
            f = np.delete(f, indepenent_comp, axis=0)

            der_free_beta = np.delete(der_free_beta, indepenent_comp, 0)
            der_free_tot = np.delete(der_free_tot, indepenent_comp, 0)
            der_free_ks = np.delete(der_free_ks, indepenent_comp, 0)

        # Solve the systems of equations
        for i in range(ns):
            solution = np.linalg.solve(M, b[:, i])
            der_free_beta[:, i] = solution[: (nc - 1 if distribution else nc)]
            if with_solids:
                der_solid_beta[:, i] = solution[(nc - 1 if distribution else nc) :]

        for r in range(nc):
            solution = np.linalg.solve(M, d[:, r])
            der_free_tot[:, r] = solution[: (nc - 1 if distribution else nc)]
            if with_solids:
                der_solid_tot[:, r] = solution[(nc - 1 if distribution else nc) :]

        if with_solids:
            for k, skip in enumerate(to_skip[-nf:]):
                if skip:
                    continue
                solution = np.linalg.solve(M, f[:, k])
                der_free_ks[:, k] = solution[: (nc - 1 if distribution else nc)]
                der_solid_ks[:, k] = solution[(nc - 1 if distribution else nc) :]

        if with_solids:
            null_solids_index = np.nonzero(c_solid == 0)[0]
            if null_solids_index.size:
                der_solid_beta = np.insert(der_solid_beta, null_solids_index, 0, axis=0)
                der_solid_tot = np.insert(der_solid_tot, null_solids_index, 0, axis=0)
                der_solid_ks = np.insert(der_solid_ks, null_solids_index, 0, axis=0)

        if distribution:
            der_free_beta = np.insert(der_free_beta, indepenent_comp, 0, axis=0)
            der_free_tot = np.insert(der_free_tot, indepenent_comp, 0, axis=0)
            der_free_ks = np.insert(der_free_ks, indepenent_comp, 0, axis=0)

        # Compute derivatives for the species
        der_spec_beta = (
            np.rot90(np.tile(stoichiometry.T, (ns, 1, 1)), -1)
            * (
                np.stack([np.tile(c_spec, (ns, 1)).T for _ in range(nc)], axis=-1)
                / c_free
            )
            * np.tile(der_free_beta.T, (ns, 1, 1))
        ).sum(axis=-1) + np.diag(c_spec / beta)

        der_spec_tot = (
            np.rot90(np.tile(stoichiometry.T, (nc, 1, 1)), -1)
            * (
                np.stack([np.tile(c_spec, (nc, 1)).T for _ in range(nc)], axis=-1)
                / c_free
            )
            * np.tile(der_free_tot.T, (ns, 1, 1))
        ).sum(axis=-1)

        der_spec_ks = (
            np.rot90(np.tile(stoichiometry.T, (nf, 1, 1)), -1)
            * (
                np.stack([np.tile(c_spec, (nf, 1)).T for _ in range(nc)], axis=-1)
                / c_free
            )
            * np.tile(der_free_ks.T, (ns, 1, 1))
        ).sum(axis=-1)

        # Calculate uncertanity for components and species given the input
        comp_sigma = np.sqrt(
            ((der_free_beta**2) * (beta_sigma**2)).sum(axis=1)
            + ((der_free_tot**2) * (conc_sigma[point] ** 2)).sum(axis=1)
            + ((der_free_ks**2) * (ks_sigma**2)).sum(axis=1)
        )

        species_sigma = np.sqrt(
            ((der_spec_beta**2) * (beta_sigma**2)).sum(axis=1)
            + ((der_spec_tot**2) * (conc_sigma[point] ** 2)).sum(axis=1)
            + ((der_spec_ks**2) * (ks_sigma**2)).sum(axis=1)
        )
        species_sigma = np.concatenate((comp_sigma, species_sigma))

        if with_solids:
            solid_sigma = np.sqrt(
                ((der_solid_beta**2) * (beta_sigma**2)).sum(axis=1)
                + ((der_solid_tot**2) * (conc_sigma[point] ** 2)).sum(axis=1)
                + ((der_solid_ks**2) * (ks_sigma**2)).sum(axis=1)
            )
        else:
            solid_sigma = np.zeros(shape=nf)

        species_sigma_result[point] = species_sigma
        solid_sigma_result[point] = solid_sigma

    return species_sigma_result, solid_sigma_result

PotentiometryOptimizer

PotentiometryOptimizer(data: SolverData, reporter=None)

Source code in src/libeq/optimizers/potentiometry.py

def PotentiometryOptimizer(data: SolverData, reporter=None):
    def f_obj(c):
        """
        Given the concentrations of the components, calculate the objective function value.

        Parameters:
        -------
        x : numpy.ndarray
            The concentrations of the components.

        Returns:
        -------
        emf : numpy.ndarray
            The calculcated potential from components.

        """
        electroactive = fhsel(c)
        calc_remf = np.log(electroactive)
        return np.ravel(calc_remf)

    def free_conc(updated_beta):
        nonlocal _initial_guess
        incoming_beta = updated_beta / 2.303
        gen = ravel(data.log_beta, incoming_beta, beta_flags)
        log_beta = np.fromiter(gen, dtype=float)
        # Solve the system of equations
        c, *_ = solve_equilibrium_equations(
            stoichiometry=stoichiometry,
            solid_stoichiometry=solid_stoichiometry,
            original_log_beta=log_beta,
            original_log_ks=original_log_ks,
            total_concentration=total_concentration,
            outer_fiexd_point_params=outer_fixed_point_params,
            initial_guess=_initial_guess,
            full=True,
        )
        _initial_guess = c[:, : stoichiometry.shape[0]]
        return c

    def jacobian(concentration):
        """
        Calculate the jacobian matrix of the objective function.

        Parameters:
        -------
        x : numpy.ndarray
            The concentrations of the components.

        Returns:
        -------
        jac : numpy.ndarray
            The jacobian matrix of the objective function.

        """
        nc = stoichiometry.shape[0]
        J = np.zeros(shape=(concentration.shape[0], nc, nc))
        diagonals = np.einsum(
            "ij,jk->ijk", concentration[:, nc:], np.eye(concentration.shape[1] - nc)
        )
        # Compute Jacobian for soluble components only
        J = stoichiometry @ diagonals @ stoichiometry.T
        J[:, range(nc), range(nc)] += concentration[:, :nc]

        B = stoichiometry[np.newaxis, ...] * concentration[..., np.newaxis, nc:]
        dcdb = np.squeeze(np.linalg.solve(J, -B))
        return fhsel(dcdb[..., np.flatnonzero(beta_flags)]).T

    def text_reporter(*args):
        print(f"iteration n.{args[0]}")
        print("x", args[1])
        print("dx", args[2])
        print("sigma", args[3])
        print("----------------\n")

    # Load the n titrations with their potential from the data file
    emf = [t.emf for t in data.potentiometry_options.titrations]
    emf0 = (t.e0 for t in data.potentiometry_options.titrations)
    slope = (t.slope for t in data.potentiometry_options.titrations)
    v_add = [t.v_add for t in data.potentiometry_options.titrations]

    ll, ul = data.potentiometry_options.px_range

    reduced_emf = [
        build_reduced_emf(emf_, emf0_, slope_)
        for emf_, emf0_, slope_ in zip(emf, emf0, slope)
    ]
    if ul + ll != 0:
        idx_to_keep = [
            (-red_emf >= ll * 2.303) & (-red_emf <= ul * 2.303)
            for red_emf in reduced_emf
        ]
        reduced_emf = [red_emf[idx] for red_emf, idx in zip(reduced_emf, idx_to_keep)]
        emf = [emf[idx] for emf, idx in zip(emf, idx_to_keep)]
        v_add = [v_add[idx] for v_add, idx in zip(v_add, idx_to_keep)]
    else:
        idx_to_keep = [None for _ in reduced_emf]

    full_emf = np.concatenate(reduced_emf, axis=0).ravel()
    n_exp_points = full_emf.shape[0]

    if data.potentiometry_options.weights == "constants":
        weights = np.ones(n_exp_points)
    elif data.potentiometry_options.weights == "calculated":
        e0_sigma = [t.e0_sigma for t in data.potentiometry_options.titrations]
        v0_sigma = [t.v0_sigma for t in data.potentiometry_options.titrations]

        weights = np.concatenate(
            [
                compute_weights(emf_, v_add_, e0_sigma_, v0_sigma_)
                for emf_, v_add_, e0_sigma_, v0_sigma_ in zip(
                    emf, v_add, e0_sigma, v0_sigma
                )
            ],
            axis=0,
        ).ravel()

    elif data.potentiometry_options.weights == "given":
        raise NotImplementedError("User given weights are not implemented yet.")

    slices = list(accumulate([0] + [s.shape[0] for s in reduced_emf]))
    electro_active_components = [
        t.electro_active_compoment for t in data.potentiometry_options.titrations
    ]
    fhsel = partial(hselect, hindices=electro_active_components, slices=slices[:-1])

    beta_flags = np.array(data.potentiometry_options.beta_flags).astype(int)
    beta_flags = np.where(beta_flags == -1, 0, beta_flags)

    (
        stoichiometry,
        solid_stoichiometry,
        original_log_beta,
        original_log_ks,
        charges,
        independent_component_activity,
    ) = _prepare_common_data(data)

    total_concentration = np.vstack(
        [
            _titration_total_c(t, i)
            for t, i in zip(data.potentiometry_options.titrations, idx_to_keep)
        ]
    )

    original_log_beta = np.tile(original_log_beta, (total_concentration.shape[0], 1))
    original_log_ks = np.tile(original_log_ks, (total_concentration.shape[0], 1))

    outer_fixed_point_params = _assemble_outer_fixed_point_params(
        data, charges, independent_component_activity
    )

    _initial_guess, *_ = solve_equilibrium_equations(
        stoichiometry=stoichiometry,
        solid_stoichiometry=solid_stoichiometry,
        original_log_beta=original_log_beta,
        original_log_ks=original_log_ks,
        total_concentration=total_concentration,
        outer_fiexd_point_params=outer_fixed_point_params,
        initial_guess=None,
        full=False,
    )

    if outer_fixed_point_params["ionic_strength_dependence"] is True:
        print(
            "Ionic strength dependence for potentiometry oprimization is not implemented yet."
        )
        outer_fixed_point_params["ionic_strength_dependence"] = False

    x, concs, return_extra = levenberg_marquardt(
        np.fromiter(unravel(data.log_beta, beta_flags), dtype=float) * 2.303,
        full_emf,
        f_obj,
        free_conc,
        jacobian,
        weights,
        report=reporter,
    )

    b_error, cor_matrix, cov_matrix = fit_final_calcs(
        return_extra["jacobian"], return_extra["residuals"], return_extra["weights"]
    )

    return x, concs, b_error, cor_matrix, cov_matrix, return_extra

EqSolver

EqSolver(
    data: SolverData,
    mode: Literal[
        "titration", "distribution"
    ] = "titration",
)

Solve the equilibrium equations for a given chemical system.

The solver uses a conjuntion of methods to solve the problem at hand. In particular:

• The Positive Continious Fraction Method (PCFM) is used to presolve the equilibrium equations.
• The Newton-Raphson method is used to solve the equilibrium equations.
• The outer fixed point method is used to solve the equilibrium equations for solids.

Parameters:

Name	Type	Description	Default
data	SolverData	The input data containing all the necessary information for the solver.	required
mode	(titration, distribution)	The mode of operation for the solver. Default is "titration".	"titration"

Returns:

Name	Type	Description
result	ndarray	The calculated equilibrium concentrations.
log_beta	ndarray	The logarithm of the stability constants.
log_ks	ndarray	The logarithm of the solubility products.
saturation_index	ndarray	The calculated saturation indices for solid phases.
total_concentration	ndarray	The total concentrations used in the calculations.

Source code in src/libeq/solver/solver.py

def EqSolver(
    data: SolverData, mode: Literal["titration", "distribution"] = "titration"
):
    """
    Solve the equilibrium equations for a given chemical system.

    The solver uses a conjuntion of methods to solve the problem at hand. In particular:

    - The Positive Continious Fraction Method (PCFM) is used to presolve the equilibrium equations.
    - The Newton-Raphson method is used to solve the equilibrium equations.
    - The outer fixed point method is used to solve the equilibrium equations for solids.

    Parameters
    ----------
    data : SolverData
        The input data containing all the necessary information for the solver.
    mode : {"titration", "distribution"}, optional
        The mode of operation for the solver. Default is "titration".

    Returns
    -------
    result : ndarray
        The calculated equilibrium concentrations.
    log_beta : ndarray
        The logarithm of the stability constants.
    log_ks : ndarray
        The logarithm of the solubility products.
    saturation_index : ndarray
        The calculated saturation indices for solid phases.
    total_concentration : ndarray
        The total concentrations used in the calculations.
    """
    # Get the total concentration values depending if mode is titration or distribution
    if mode == "titration":
        (
            stoichiometry,
            solid_stoichiometry,
            original_log_beta,
            original_log_ks,
            charges,
            independent_component_activity,
            total_concentration,
        ) = _prepare_titration_data(data)
    elif mode == "distribution":
        (
            stoichiometry,
            solid_stoichiometry,
            original_log_beta,
            original_log_ks,
            charges,
            independent_component_activity,
            total_concentration,
            independent_component,
            independent_component_concentration,
        ) = _prepare_distribution_data(data)
    else:
        raise ValueError("Invalid work mode")

    outer_fiexd_point_params = _assemble_outer_fixed_point_params(
        data, charges, independent_component_activity
    )

    # Solve the equilibrium equations
    result, log_beta, log_ks, saturation_index, total_concentration = (
        solve_equilibrium_equations(
            stoichiometry=stoichiometry,
            solid_stoichiometry=solid_stoichiometry,
            original_log_beta=original_log_beta,
            original_log_ks=original_log_ks,
            total_concentration=total_concentration,
            outer_fiexd_point_params=outer_fiexd_point_params,
        )
    )

    if mode == "distribution":
        result, log_beta, log_ks, total_concentration = _expand_result(
            result,
            independent_component,
            independent_component_concentration,
            total_concentration,
            log_beta,
            log_ks,
            data.stoichiometry,
            data.solid_stoichiometry,
        )
    return result, log_beta, log_ks, saturation_index, total_concentration