From 5c727981a33bfaeb558c5834f68af6aafdc02248 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Thu, 6 Aug 2026 08:17:24 +0800 Subject: [PATCH 01/10] refactor(core): redox/site-occupancy helpers with instance-aware Fe split --- src/phasetools/core/README.md | 21 ++ src/phasetools/core/base.py | 4 +- src/phasetools/core/phase_properties.py | 330 +++++++++++++++--------- 3 files changed, 224 insertions(+), 131 deletions(-) diff --git a/src/phasetools/core/README.md b/src/phasetools/core/README.md index c87f88e..6e1b37d 100644 --- a/src/phasetools/core/README.md +++ b/src/phasetools/core/README.md @@ -22,3 +22,24 @@ The `core` submodule provides the foundational classes and low-level bridging lo - `get_phase_mg2_number`: Phase-wide $Mg\#$ (Divalent Iron only). - `get_phase_fe_split`: Heuristic splitting of total iron into $\text{Fe}^{2+}$ and $\text{Fe}^{3+}$. - `calculate_kd_fe_mg`: Distribution coefficients between phases. + +### Solvus (multi-instance) handling +MAGEMin reports coexisting solvus limbs as repeated entries in `out.ph` +(e.g. two `dio` clinopyroxenes or two `amp` amphiboles). `phase_frac` +sums them; the per-phase composition helpers take an `instance` argument: + +- integer index (default `0`) — that instance, with a `UserWarning` + noting how many other instances exist (`_phase_indices` raises + `ValueError` for non-integer, non-`'all'` values); +- `'all'` — one value per instance, returned as numpy arrays + (`get_oxide_apfu` → `{oxide: ndarray}`, `extract_end_member`, + `get_phase_mg_number`, `get_phase_mg2_number`, `get_phase_fe_split`). + +`MAGEMinPTGridCalculator.extract_from_grid` / `single_point_calc` / +`generate_2D_grid` accept the same `instance` argument. With +`instance='all'` they emit one column per instance keyed with a +`_0`, `_1`, ... suffix (e.g. `ox_apfu_Na2O_0`); grid points where the +phase has fewer instances are NaN-padded. `mol_frac`/`wt_frac`/ +`vol_frac` always give the summed totals (matching `phase_frac`); with +`instance='all'` additional per-instance columns `mol_frac_0`, +`mol_frac_1`, ... are emitted alongside them. diff --git a/src/phasetools/core/base.py b/src/phasetools/core/base.py index 8a13bde..d81bd75 100644 --- a/src/phasetools/core/base.py +++ b/src/phasetools/core/base.py @@ -75,9 +75,9 @@ def setup_bulk_composition(self, Xoxides, X, sys_in, rm_list=None): except: pass - def _extract_fe_split_from_apfu(self, out, phase): + def _extract_fe_split_from_apfu(self, out, phase, instance=0): """ Internal: Calculate Fe2+ and Fe3+ amounts for a phase using an excess oxygen heuristic. """ from .phase_properties import get_phase_fe_split - return get_phase_fe_split(out, phase) + return get_phase_fe_split(out, phase, instance=instance) diff --git a/src/phasetools/core/phase_properties.py b/src/phasetools/core/phase_properties.py index f54efe8..02a0e1f 100644 --- a/src/phasetools/core/phase_properties.py +++ b/src/phasetools/core/phase_properties.py @@ -1,24 +1,93 @@ import numpy as np +import warnings -def get_oxide_apfu(out, ph, oxides): - """Extract oxide amounts from APFU output for a specific phase.""" - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - oxide_values = np.array(phase_obj.Comp_apfu, dtype=float) - oxide_names = [str(ox) for ox in out.oxides] - oxide_moles_dict = {ox: value for ox, value in zip(oxide_names, oxide_values)} +def _phase_indices(out, phase, instance=0): + """Return the index/indices of ``phase`` in ``out.ph``. - results = {} - for oxide in oxides: - results[oxide] = oxide_moles_dict.get(oxide, 0.0) - except (ValueError, IndexError): - results = {oxide: 0.0 for oxide in oxides} + Parameters + ---------- + out : object + MAGEMin output object with a ``ph`` attribute. + phase : str + Phase name. + instance : int or {'all'}, default=0 + Which instance to resolve. An integer index selects that instance + (negative indices count from the end). ``'all'`` returns every + occurrence. When the phase appears more than once (a solvus -- e.g. + two coexisting clinopyroxenes or amphiboles, reported by MAGEMin as + repeated entries in ``out.ph``), a warning notes how many other + instances exist. + + Returns + ------- + list[int] + Indices of ``phase`` in ``out.ph`` (empty list if the phase is + absent or the requested instance does not exist). + """ + idx = [i for i, p in enumerate(out.ph) if str(p) == phase] + if instance == 'all': + return idx + if not isinstance(instance, (int, np.integer)): + raise ValueError(f"instance must be an integer index or 'all', got {instance!r}") + if not idx: + return [] + if instance < 0: + instance = len(idx) + instance + if instance < 0 or instance >= len(idx): + warnings.warn( + f"phase '{phase}' has {len(idx)} instance(s); requested instance " + f"{instance} does not exist -- treating as absent.", + UserWarning, stacklevel=3) + return [] + if len(idx) > 1: + warnings.warn( + f"phase '{phase}' has {len(idx)} instances (solvus) in this " + f"output; using instance {instance} (there are {len(idx) - 1} " + f"others). Use instance='all' for per-instance arrays.", + UserWarning, stacklevel=3) + return [idx[instance]] - return results +def get_oxide_apfu(out, ph, oxides, instance=0): + """Extract oxide amounts from APFU output for a specific phase. -def get_phase_chemistry(out, ph, oxides, sys_in): + Parameters + ---------- + out : object + MAGEMin output object. + ph : str + Phase name. + oxides : list[str] + Oxides to extract. + instance : int or {'all'}, default=0 + For a phase that appears multiple times (a solvus), an integer + index selects that instance (warning if more than one exists) and + ``'all'`` returns one value per instance as numpy arrays. + + Returns + ------- + dict + ``{oxide: value}`` for an integer index, or ``{oxide: ndarray}`` + (one entry per phase instance) for ``'all'``. + """ + idx = _phase_indices(out, ph, instance) + if not idx: + return {oxide: 0.0 for oxide in oxides} if instance != 'all' \ + else {oxide: np.zeros(0) for oxide in oxides} + oxide_names = [str(ox) for ox in out.oxides] + per = [] + for i in idx: + try: + values = np.array(out.SS_vec[i].Comp_apfu, dtype=float) + per.append({ox: float(values[oxide_names.index(ox)]) + if ox in oxide_names else 0.0 for ox in oxides}) + except (ValueError, IndexError, AttributeError): + per.append({ox: 0.0 for ox in oxides}) + if instance == 'all': + return {ox: np.array([d[ox] for d in per]) for ox in oxides} + return per[0] + +def get_phase_chemistry(out, ph, oxides, sys_in, instance=0): """ Extract oxide concentrations (wt% or mol%) for a specific phase. @@ -32,46 +101,61 @@ def get_phase_chemistry(out, ph, oxides, sys_in): List of oxides to extract. sys_in : str Unit system ('wt' or 'mol'). + instance : int or {'all'}, default=0 + For a phase that appears multiple times (a solvus), an integer + index selects that instance (warning if more than one exists) and + ``'all'`` returns one value per instance as numpy arrays. Returns ------- dict - Oxide concentrations. + Oxide concentrations for an integer index, or per-instance arrays + for ``'all'``. """ - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - - if sys_in.casefold() == 'wt': - # Comp_wt is weight fraction (0-1) for oxides in the phase - values = np.array(phase_obj.Comp_wt, dtype=float) * 100.0 - else: - # Comp is molar fraction (0-1) for oxides in the phase - values = np.array(phase_obj.Comp, dtype=float) * 100.0 - - oxide_names = [str(ox) for ox in out.oxides] - oxide_dict = {ox: value for ox, value in zip(oxide_names, values)} - - results = {} - for oxide in oxides: - results[oxide] = oxide_dict.get(oxide, 0.0) - except (ValueError, IndexError): - results = {oxide: 0.0 for oxide in oxides} + idx = _phase_indices(out, ph, instance) + if not idx: + return {oxide: 0.0 for oxide in oxides} if instance != 'all' \ + else {oxide: np.zeros(0) for oxide in oxides} + oxide_names = [str(ox) for ox in out.oxides] + per = [] + for i in idx: + try: + phase_obj = out.SS_vec[i] + if sys_in.casefold() == 'wt': + # Comp_wt is weight fraction (0-1) for oxides in the phase + values = np.array(phase_obj.Comp_wt, dtype=float) * 100.0 + else: + # Comp is molar fraction (0-1) for oxides in the phase + values = np.array(phase_obj.Comp, dtype=float) * 100.0 + per.append({ox: float(values[oxide_names.index(ox)]) + if ox in oxide_names else 0.0 for ox in oxides}) + except (ValueError, IndexError, AttributeError): + per.append({ox: 0.0 for ox in oxides}) + if instance == 'all': + return {ox: np.array([d[ox] for d in per]) for ox in oxides} + return per[0] - return results +def extract_end_member(phase, MAGEMinOutput, end_member, sys_in, instance=0): + """Extract specific end-member fraction from MAGEMin output. -def extract_end_member(phase, MAGEMinOutput, end_member, sys_in): - """Extract specific end-member fraction from MAGEMin output.""" - try: - phase_ind = MAGEMinOutput.ph.index(phase) - em_index = MAGEMinOutput.SS_vec[phase_ind].emNames.index(end_member) - if sys_in.casefold() == 'wt': - data = MAGEMinOutput.SS_vec[phase_ind].emFrac_wt[em_index] - else: - data = MAGEMinOutput.SS_vec[phase_ind].emFrac[em_index] - except (ValueError, IndexError): - data = 0. - return data + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as a numpy array. + """ + idx = _phase_indices(MAGEMinOutput, phase, instance) + if not idx: + return 0.0 if instance != 'all' else np.zeros(0) + vals = [] + for i in idx: + try: + em_index = MAGEMinOutput.SS_vec[i].emNames.index(end_member) + if sys_in.casefold() == 'wt': + vals.append(float(MAGEMinOutput.SS_vec[i].emFrac_wt[em_index])) + else: + vals.append(float(MAGEMinOutput.SS_vec[i].emFrac[em_index])) + except (ValueError, IndexError, AttributeError): + vals.append(0.0) + return float(vals[0]) if instance != 'all' else np.array(vals) def phase_frac(phase, MAGEMinOutput, sys_in): """ @@ -97,7 +181,7 @@ def phase_frac(phase, MAGEMinOutput, sys_in): except: return 0.0 -def get_phase_mg_number(out, ph): +def get_phase_mg_number(out, ph, instance=0): """ Calculate Mg# (molar Mg / (Mg + Fe_total)) for a specific phase. @@ -107,109 +191,97 @@ def get_phase_mg_number(out, ph): Matches the logic used by MAGEMin's 'ss_MgNum' mode by pulling MgO and FeO directly from the phase's Comp_apfu array. Supports 'FeO', 'Fe' (sb24), and 'Fe2O3' fallback. + + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as a numpy array. """ - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - - oxide_names = [str(ox) for ox in out.oxides] - - # Pull Mg + idx = _phase_indices(out, ph, instance) + if not idx: + return 0.0 if instance != 'all' else np.zeros(0) + oxide_names = [str(ox) for ox in out.oxides] + + def _mg(i): try: - mg_idx = oxide_names.index('MgO') - mg = float(phase_obj.Comp_apfu[mg_idx]) - except ValueError: - mg = 0.0 - - # Pull Fe (total iron atoms) - fe = 0.0 - if 'FeO' in oxide_names: - fe_idx = oxide_names.index('FeO') - fe += float(phase_obj.Comp_apfu[fe_idx]) - # If both are present, we sum them (though unlikely in standard MAGEMin output) - if 'Fe2O3' in oxide_names: - fe2o3_idx = oxide_names.index('Fe2O3') - fe += 2.0 * float(phase_obj.Comp_apfu[fe2o3_idx]) - elif 'Fe' in oxide_names: - fe_idx = oxide_names.index('Fe') - fe += float(phase_obj.Comp_apfu[fe_idx]) - elif 'Fe2O3' in oxide_names: - fe2o3_idx = oxide_names.index('Fe2O3') - fe += 2.0 * float(phase_obj.Comp_apfu[fe2o3_idx]) - else: - # No iron components found - if mg == 0: return 0.0 - return 1.0 # Pure Mg phase - - denominator = mg + fe - if denominator == 0: + phase_obj = out.SS_vec[i] + mg = float(phase_obj.Comp_apfu[oxide_names.index('MgO')]) if 'MgO' in oxide_names else 0.0 + fe = 0.0 + if 'FeO' in oxide_names: + fe += float(phase_obj.Comp_apfu[oxide_names.index('FeO')]) + if 'Fe2O3' in oxide_names: + fe += 2.0 * float(phase_obj.Comp_apfu[oxide_names.index('Fe2O3')]) + elif 'Fe' in oxide_names: + fe += float(phase_obj.Comp_apfu[oxide_names.index('Fe')]) + elif 'Fe2O3' in oxide_names: + fe += 2.0 * float(phase_obj.Comp_apfu[oxide_names.index('Fe2O3')]) + else: + return 1.0 if mg > 0 else 0.0 + denom = mg + fe + return mg / denom if denom else 0.0 + except (ValueError, IndexError, AttributeError): return 0.0 - return mg / denominator - except (ValueError, IndexError, AttributeError): - return 0.0 + vals = [_mg(i) for i in idx] + return float(vals[0]) if instance != 'all' else np.array(vals) -def get_phase_fe_split(out, ph): +def get_phase_fe_split(out, ph, instance=0): """ Calculate Fe2+ and Fe3+ amounts for a phase using an excess oxygen heuristic. Works for both traditional FeO-Fe2O3 bases and MAGEMin's O-basis (ig, mp). + + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as numpy arrays. """ - try: - ox_to_query = ['FeO', 'Fe2O3', 'Fe', 'O'] - apfu = get_oxide_apfu(out, ph, ox_to_query) - - feo_val = apfu.get("FeO", 0.0) - fe2o3_val = apfu.get("Fe2O3", 0.0) - fe_metal_val = apfu.get("Fe", 0.0) - atomic_o = apfu.get("O", 0.0) - - # 1. Calculate Total Fe atoms (Atoms per formula unit) - if atomic_o > 0: - # MAGEMin O-basis (ig, mp) or sb24 basis - # If Fe component is present (sb24), use it; otherwise FeO is total iron. - if fe_metal_val > 0: - total_fe = fe_metal_val - else: - total_fe = feo_val - else: - # Traditional FeO/Fe2O3 basis - total_fe = feo_val + 2.0 * fe2o3_val - - # 2. Calculate Fe3+ atoms using excess oxygen heuristic - # excess_o identifies oxygen atoms added beyond the stoichiometric baseline. - # Works for both 'O as total oxygen' and 'O as excess oxygen' components. - excess_o = max(atomic_o - np.round(atomic_o), 0.0) - fe3 = 2.0 * fe2o3_val + 2.0 * excess_o - - # 3. Divalent iron is the remainder - fe2 = max(total_fe - fe3, 0.0) + apfu = get_oxide_apfu(out, ph, ['FeO', 'Fe2O3', 'Fe', 'O'], instance=instance) - return { - "fe2": fe2, - "fe3": fe3, - } - except: - return {"fe2": 0.0, "fe3": 0.0} + feo = np.asarray(apfu.get("FeO", 0.0), dtype=float) + fe2o3 = np.asarray(apfu.get("Fe2O3", 0.0), dtype=float) + fem = np.asarray(apfu.get("Fe", 0.0), dtype=float) + ato = np.asarray(apfu.get("O", 0.0), dtype=float) + + if ato.size == 0: + empty = {"Fe2": np.zeros(0), "Fe3": np.zeros(0)} + return {"Fe2": 0.0, "Fe3": 0.0} if instance != 'all' else empty -def get_phase_mg2_number(out, ph): + # 1. Calculate Total Fe atoms (Atoms per formula unit) + if np.any(ato > 0): + # MAGEMin O-basis (ig, mp) or sb24 basis + total_fe = fem if np.any(fem > 0) else feo + else: + # Traditional FeO/Fe2O3 basis + total_fe = feo + 2.0 * fe2o3 + + # 2. Calculate Fe3+ atoms using excess oxygen heuristic + excess_o = np.maximum(ato - np.floor(ato), 0.0) + fe3 = 2.0 * fe2o3 + 2.0 * excess_o + + # 3. Divalent iron is the remainder + fe2 = np.maximum(np.asarray(total_fe) - fe3, 0.0) + + if instance != 'all': + return {"Fe2": float(np.squeeze(fe2)), "Fe3": float(np.squeeze(fe3))} + return {"Fe2": fe2, "Fe3": fe3} + +def get_phase_mg2_number(out, ph, instance=0): """ Calculate Mg# (molar Mg / (Mg + Fe2+)) for a specific phase. Uses an excess oxygen heuristic to split total iron into Fe2+ and Fe3+. """ try: - apfu = get_oxide_apfu(out, ph, ['MgO']) + apfu = get_oxide_apfu(out, ph, ['MgO'], instance=instance) mg = apfu.get('MgO', 0.0) - split = get_phase_fe_split(out, ph) - fe2 = split['fe2'] + split = get_phase_fe_split(out, ph, instance=instance) + fe2 = split['Fe2'] denominator = mg + fe2 - if denominator == 0: - return 0.0 + if np.any(np.asarray(denominator) <= 0): + return 0.0 if instance != 'all' else np.zeros(len(mg)) - return mg / denominator + return (mg / denominator) if instance != 'all' else np.asarray(mg / denominator, dtype=float) except: return 0.0 From da59f4ed00b5d51a792a20f5d2dfb37233e9c3cf Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Thu, 6 Aug 2026 08:17:09 +0800 Subject: [PATCH 02/10] feat(utils): add FeOt/FeO/O and non-normalising bulk-rock converters --- src/phasetools/utils/README.md | 2 + src/phasetools/utils/bulk_rock.py | 181 ++++++++++++++++++++++++++++++ 2 files changed, 183 insertions(+) diff --git a/src/phasetools/utils/README.md b/src/phasetools/utils/README.md index 94aa1c2..ccd6c91 100644 --- a/src/phasetools/utils/README.md +++ b/src/phasetools/utils/README.md @@ -7,6 +7,8 @@ The `utils` submodule contains helper functions for chemistry, math, and general ### `bulk_rock.py` - Core stoichiometric engine for managing bulk compositions. - Provides molar mass lookups and unit conversion helpers (e.g., mol% to wt% and vice versa). +- `mol_fractions_to_wt_fractions` / `wt_fractions_to_mol_fractions`: non-normalising converters for single components or subsets (e.g., FeO/Fe2O3/FeOt iron redox conversions). +- `split_feot_to_feo_o` / `express_bulk_in_feo_o_basis`: split total iron (FeOt) into the MAGEMin **FeO + O** redox basis at a target Fe³⁺/FeOt fraction, **conserving FeOt** (FeO column = total Fe; O = f·FeOt/2 per 2FeO + O → Fe₂O₃). Useful for redox sweeps and for converting measured (FeOt-only) bulks into MAGEMin format. - Standardises oxide lists and manages cation-oxide mapping. ### `general.py` diff --git a/src/phasetools/utils/bulk_rock.py b/src/phasetools/utils/bulk_rock.py index 752d683..ea799c0 100644 --- a/src/phasetools/utils/bulk_rock.py +++ b/src/phasetools/utils/bulk_rock.py @@ -44,6 +44,25 @@ def convert_mol_percent_to_wt_percent(mol_percents, components, mass_dict): wt_percents = [(mol * mass_dict[comp] / total_mass) * 100 for comp, mol in zip(components, mol_percents)] return wt_percents +def atomic_frac_to_wt_frac(atomic_frac, mass_dict): + """Convert atomic (molar) site fractions to weight-based site fractions. + + Parameters + ---------- + atomic_frac : dict + Mapping of component names to their atomic fractions (summing to 1.0). + mass_dict : dict + Mapping of component names to their atomic/molecular masses. + + Returns + ------- + dict + Weight-based fractions (summing to 1.0). + """ + total_mass = sum(atomic_frac[k] * mass_dict[k] for k in atomic_frac) + return {k: atomic_frac[k] * mass_dict[k] / total_mass for k in atomic_frac} + + def convert_wt_percent_to_mol_percent(wt_percents, components, mass_dict): """Generic conversion from weight (mass) percent to mole percent.""" total_moles = 0 @@ -52,6 +71,68 @@ def convert_wt_percent_to_mol_percent(wt_percents, components, mass_dict): mol_percents = [((wt / mass_dict[comp]) / total_moles) * 100 for comp, wt in zip(components, wt_percents)] return mol_percents +def mol_fractions_to_wt_fractions(mol, components, mass_dict): + """Convert mole fractions to weight fractions (no normalisation). + + Unlike :func:`convert_mol_percent_to_wt_percent`, the output is not + normalised to sum to 100. This means the function can be applied to a + single component (e.g. a single oxide for an iron redox conversion) or + a subset of a composition, as well as a full composition. + + Parameters + ---------- + mol : float or array_like + Mole fraction(s) of each component. + components : list of str + Component names corresponding to each input value. + mass_dict : dict + Mapping of component names to molecular masses. + + Returns + ------- + float or list or numpy.ndarray + Weight fraction(s) of each component. A scalar input returns a + scalar, a list input returns a list, and any other array-like + input returns a numpy array. + """ + if np.isscalar(mol): + return float(mol) * mass_dict[components[0]] + mol_arr = np.asarray(mol, dtype=float) + masses = np.array([mass_dict[comp] for comp in components], dtype=float) + result = mol_arr * masses + return result.tolist() if isinstance(mol, list) else result + +def wt_fractions_to_mol_fractions(wt, components, mass_dict): + """Convert weight fractions to mole fractions (no normalisation). + + Unlike :func:`convert_wt_percent_to_mol_percent`, the output is not + normalised to sum to 100. This means the function can be applied to a + single component (e.g. a single oxide for an iron redox conversion) or + a subset of a composition, as well as a full composition. + + Parameters + ---------- + wt : float or array_like + Weight fraction(s) of each component. + components : list of str + Component names corresponding to each input value. + mass_dict : dict + Mapping of component names to molecular masses. + + Returns + ------- + float or list or numpy.ndarray + Mole fraction(s) of each component. A scalar input returns a + scalar, a list input returns a list, and any other array-like + input returns a numpy array. + """ + if np.isscalar(wt): + return float(wt) / mass_dict[components[0]] + wt_arr = np.asarray(wt, dtype=float) + masses = np.array([mass_dict[comp] for comp in components], dtype=float) + result = wt_arr / masses + return result.tolist() if isinstance(wt, list) else result + def convert_wt_percent_to_moles(wt_percents, components, mass_dict, total_weight): """Convert weight (mass) percentages to moles.""" moles = [] @@ -76,3 +157,103 @@ def convert_moles_to_mol_percent(moles, components): total = sum(moles_dict.values()) return {comp: (moles_dict[comp] / total) * 100 for comp in components} + +def split_feot_to_feo_o(feot_moles, fe3_frac): + """ + Split total iron (FeOt) into the MAGEMin ``FeO + O`` redox pair at a + target Fe3+/FeOt fraction, conserving the total iron budget. + + Parameters + ---------- + feot_moles : float + Total iron in mole units (atoms of Fe, equivalently the amount of + FeO that would carry all iron as Fe2+). + fe3_frac : float + Target Fe3+/FeOt fraction in ``[0, 1]``. 0 = fully reduced (all + Fe2+), 1 = fully oxidised (all Fe3+). + + Returns + ------- + (feo, o) : tuple[float, float] + ``feo`` is the total-iron column (all Fe expressed as FeO, equal to + ``feot_moles``) and ``o`` is the excess oxygen required to oxidise + the target fraction, per ``2FeO + O -> Fe2O3``. + + Notes + ----- + Molar bookkeeping (each Fe2O3 carries 2 Fe atoms and needs 1 O): + Fe3+ atoms = 2 * O => O = fe3_frac * FeOt / 2 + Fe2+ atoms = FeOt - 2 * O (implicitly held by the FeO column) + """ + fe3_frac = float(fe3_frac) + if not 0.0 <= fe3_frac <= 1.0: + raise ValueError(f"fe3_frac must be in [0, 1], got {fe3_frac!r}") + feot_moles = float(feot_moles) + return feot_moles, fe3_frac * feot_moles / 2.0 + +def express_bulk_in_feo_o_basis(X, Xoxides, fe3_frac, feo_oxide="FeO", + fe2o3_oxide="Fe2O3", o_oxide="O"): + """ + Express a bulk composition in the MAGEMin ``FeO + O`` redox basis at a + target Fe3+/FeOt fraction, conserving total iron. + + The ``FeO`` column is set to total iron (all Fe expressed as FeO, + ``FeOt``) and ``O`` is set to the excess oxygen giving the requested + Fe3+/FeOt partition. Any ``Fe2O3`` component is removed (set to 0) and + an existing ``O`` component is overwritten; missing components are + appended. The returned list preserves the input oxide order. + + Parameters + ---------- + X : array-like + Bulk composition values in MOLE units (mol fractions or mol%). Use + weight-based converters first if the input is in wt%. + Xoxides : list of str + Oxide names corresponding to ``X`` (may include 'FeO', 'Fe2O3', 'O' + in any combination). + fe3_frac : float + Target Fe3+/FeOt fraction in ``[0, 1]``. + feo_oxide, fe2o3_oxide, o_oxide : str + Component names used in ``Xoxides``. + + Returns + ------- + (X_new, Xoxides_new) : tuple[list, list] + Composition in the ``FeO + O`` basis (unnormalised -- pass to + ``convertBulk4MAGEMin`` or normalise afterwards). + + Examples + -------- + >>> X = [50.0, 8.0, 0.5] # SiO2, FeO, O (mol%) + >>> ox = ['SiO2', 'FeO', 'O'] + >>> X2, ox2 = express_bulk_in_feo_o_basis(X, ox, fe3_frac=0.1) + >>> ox2 + ['SiO2', 'FeO', 'O'] + >>> X2[2] == 0.1 * X[1] / 2.0 # O = fe3_frac * FeOt / 2 + True + """ + X = [float(v) for v in X] + Xoxides = list(Xoxides) + feo_i = Xoxides.index(feo_oxide) if feo_oxide in Xoxides else None + fe2o3_i = Xoxides.index(fe2o3_oxide) if fe2o3_oxide in Xoxides else None + o_i = Xoxides.index(o_oxide) if o_oxide in Xoxides else None + + feo_mol = X[feo_i] if feo_i is not None else 0.0 + fe2o3_mol = X[fe2o3_i] if fe2o3_i is not None else 0.0 + feot = feo_mol + 2.0 * fe2o3_mol # total Fe atoms (moles) + + feo_total, o_excess = split_feot_to_feo_o(feot, fe3_frac) + + if feo_i is not None: + X[feo_i] = feo_total + else: + X.append(feo_total) + Xoxides.append(feo_oxide) + if fe2o3_i is not None: + X[fe2o3_i] = 0.0 + if o_i is not None: + X[o_i] = o_excess + else: + X.append(o_excess) + Xoxides.append(o_oxide) + return X, Xoxides From cfa6d7a1b2d4de8cd55f09757bc11401aeb652e7 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:19:16 +0800 Subject: [PATCH 03/10] fix(utils): type hints, 0-d array edge case, duplicate oxide validation --- src/phasetools/utils/bulk_rock.py | 41 ++++++++++++++++++++++++++----- 1 file changed, 35 insertions(+), 6 deletions(-) diff --git a/src/phasetools/utils/bulk_rock.py b/src/phasetools/utils/bulk_rock.py index ea799c0..00da5f2 100644 --- a/src/phasetools/utils/bulk_rock.py +++ b/src/phasetools/utils/bulk_rock.py @@ -44,7 +44,7 @@ def convert_mol_percent_to_wt_percent(mol_percents, components, mass_dict): wt_percents = [(mol * mass_dict[comp] / total_mass) * 100 for comp, mol in zip(components, mol_percents)] return wt_percents -def atomic_frac_to_wt_frac(atomic_frac, mass_dict): +def atomic_frac_to_wt_frac(atomic_frac: dict[str, float], mass_dict: dict[str, float]) -> dict[str, float]: """Convert atomic (molar) site fractions to weight-based site fractions. Parameters @@ -71,7 +71,9 @@ def convert_wt_percent_to_mol_percent(wt_percents, components, mass_dict): mol_percents = [((wt / mass_dict[comp]) / total_moles) * 100 for comp, wt in zip(components, wt_percents)] return mol_percents -def mol_fractions_to_wt_fractions(mol, components, mass_dict): +def mol_fractions_to_wt_fractions( + mol: float | list | np.ndarray, components: list[str], mass_dict: dict[str, float] +) -> float | list | np.ndarray: """Convert mole fractions to weight fractions (no normalisation). Unlike :func:`convert_mol_percent_to_wt_percent`, the output is not @@ -94,15 +96,24 @@ def mol_fractions_to_wt_fractions(mol, components, mass_dict): Weight fraction(s) of each component. A scalar input returns a scalar, a list input returns a list, and any other array-like input returns a numpy array. + + Notes + ----- + Scalar input assumes a single component -- pass ``components=[comp]`` + with the oxide name. """ if np.isscalar(mol): return float(mol) * mass_dict[components[0]] mol_arr = np.asarray(mol, dtype=float) + if mol_arr.ndim == 0: + return float(mol_arr.item() * mass_dict[components[0]]) masses = np.array([mass_dict[comp] for comp in components], dtype=float) result = mol_arr * masses return result.tolist() if isinstance(mol, list) else result -def wt_fractions_to_mol_fractions(wt, components, mass_dict): +def wt_fractions_to_mol_fractions( + wt: float | list | np.ndarray, components: list[str], mass_dict: dict[str, float] +) -> float | list | np.ndarray: """Convert weight fractions to mole fractions (no normalisation). Unlike :func:`convert_wt_percent_to_mol_percent`, the output is not @@ -129,6 +140,8 @@ def wt_fractions_to_mol_fractions(wt, components, mass_dict): if np.isscalar(wt): return float(wt) / mass_dict[components[0]] wt_arr = np.asarray(wt, dtype=float) + if wt_arr.ndim == 0: + return float(wt_arr.item() / mass_dict[components[0]]) masses = np.array([mass_dict[comp] for comp in components], dtype=float) result = wt_arr / masses return result.tolist() if isinstance(wt, list) else result @@ -158,7 +171,7 @@ def convert_moles_to_mol_percent(moles, components): total = sum(moles_dict.values()) return {comp: (moles_dict[comp] / total) * 100 for comp in components} -def split_feot_to_feo_o(feot_moles, fe3_frac): +def split_feot_to_feo_o(feot_moles: float, fe3_frac: float) -> tuple[float, float]: """ Split total iron (FeOt) into the MAGEMin ``FeO + O`` redox pair at a target Fe3+/FeOt fraction, conserving the total iron budget. @@ -191,8 +204,14 @@ def split_feot_to_feo_o(feot_moles, fe3_frac): feot_moles = float(feot_moles) return feot_moles, fe3_frac * feot_moles / 2.0 -def express_bulk_in_feo_o_basis(X, Xoxides, fe3_frac, feo_oxide="FeO", - fe2o3_oxide="Fe2O3", o_oxide="O"): +def express_bulk_in_feo_o_basis( + X: list[float], + Xoxides: list[str], + fe3_frac: float, + feo_oxide: str = "FeO", + fe2o3_oxide: str = "Fe2O3", + o_oxide: str = "O", +) -> tuple[list[float], list[str]]: """ Express a bulk composition in the MAGEMin ``FeO + O`` redox basis at a target Fe3+/FeOt fraction, conserving total iron. @@ -232,6 +251,16 @@ def express_bulk_in_feo_o_basis(X, Xoxides, fe3_frac, feo_oxide="FeO", >>> X2[2] == 0.1 * X[1] / 2.0 # O = fe3_frac * FeOt / 2 True """ + if len(Xoxides) != len(X): + raise ValueError( + f"X ({len(X)} items) and Xoxides ({len(Xoxides)} items) must have the same length" + ) + seen = set() + for ox in Xoxides: + if ox in seen: + raise ValueError(f"Duplicate oxide name in Xoxides: {ox!r}") + seen.add(ox) + X = [float(v) for v in X] Xoxides = list(Xoxides) feo_i = Xoxides.index(feo_oxide) if feo_oxide in Xoxides else None From f7875937bd06a14e7778b056bb593d6b06fa6e8d Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:40:28 +0800 Subject: [PATCH 04/10] refactor(calculators): add juliacall imports, grid mesh, and instance-aware extraction to pt_grid --- src/phasetools/calculators/pt_grid.py | 247 +++++++++++++++++--------- 1 file changed, 163 insertions(+), 84 deletions(-) diff --git a/src/phasetools/calculators/pt_grid.py b/src/phasetools/calculators/pt_grid.py index f5d1f94..8c2fb74 100644 --- a/src/phasetools/calculators/pt_grid.py +++ b/src/phasetools/calculators/pt_grid.py @@ -1,8 +1,9 @@ import numpy as np import sys +from juliacall import Main as jl, convert as jlconvert from ..core.base import MAGEMinBase -from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number -from ..utils.bulk_rock import atomic_mass_dict, convert_mol_percent_to_wt_percent +from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, _phase_indices +from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac from phasetools import MAGEMin_C class MAGEMinPTGridCalculator(MAGEMinBase): @@ -20,17 +21,19 @@ def calculate_grid(self, P, T): P = np.atleast_1d(P) T = np.atleast_1d(T) - if P.shape != T.shape: - if P.ndim == 1 and T.ndim == 1: - P_orig, T_orig = P, T - P, T = np.meshgrid(P_orig, T_orig) - P = P.flatten() - T = T.flatten() - else: - raise ValueError(f"P and T must have the same shape or both be 1D. Got {P.shape} and {T.shape}") + if P.ndim == 1 and T.ndim == 1 and P.shape[0] != T.shape[0]: + P_orig, T_orig = P, T + P, T = np.meshgrid(P_orig, T_orig) + P = P.flatten() + T = T.flatten() + elif P.shape != T.shape: + raise ValueError(f"P and T must have the same shape or both be 1D. Got {P.shape} and {T.shape}") + + P_jl = jlconvert(jl.Vector[jl.Float64], P.astype(float)) + T_jl = jlconvert(jl.Vector[jl.Float64], T.astype(float)) out = MAGEMin_C.multi_point_minimization( - P, T, self.data, X=self.X, Xoxides=self.Xoxides, + P_jl, T_jl, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list ) sys.stdout.flush() @@ -77,38 +80,41 @@ def get_phase_endmembers(self, phase, grid_out=None): if out is None: return [] - # If it's a single output object (from single_point_calc) - if not isinstance(out, (list, np.ndarray)): + # If it's a single output object (from single_point_calc) — not a list/array/Vector + if not isinstance(out, (list, np.ndarray)) and hasattr(out, 'ph'): if phase in out.ph: ph_index = out.ph.index(phase) - return [str(n) for n in out.SS_vec[ph_index].emNames] + if ph_index < len(out.SS_vec): + return [str(n) for n in out.SS_vec[ph_index].emNames] return [] - # If it's a grid + # If it's a grid (list, numpy array, or Julia Vector) for o in out: if phase in o.ph: ph_index = o.ph.index(phase) - return [str(n) for n in o.SS_vec[ph_index].emNames] + if ph_index < len(o.SS_vec): + return [str(n) for n in o.SS_vec[ph_index].emNames] + return [] return [] - def _extract_cations_from_apfu(self, out, phase, cations, sys_in): + def _extract_cations_from_apfu(self, out, phase, cations, sys_in, instance=0): """Internal: Extract cation ratios (e.g., XMg, XFe) for a specific phase.""" ox_to_query = ['MgO', 'MnO', 'CaO', 'FeO', 'Fe2O3', 'Fe', 'O'] - apfu = get_oxide_apfu(out, phase, ox_to_query) - - mg = apfu.get("MgO", 0.0) - mn = apfu.get("MnO", 0.0) - ca = apfu.get("CaO", 0.0) - - feo = apfu.get("FeO", 0.0) - fe2o3 = apfu.get("Fe2O3", 0.0) - fe_metal = apfu.get("Fe", 0.0) - atomic_o = apfu.get("O", 0.0) + apfu = get_oxide_apfu(out, phase, ox_to_query, instance=instance) + + mg = np.asarray(apfu.get("MgO", 0.0), dtype=float) + mn = np.asarray(apfu.get("MnO", 0.0), dtype=float) + ca = np.asarray(apfu.get("CaO", 0.0), dtype=float) + + feo = np.asarray(apfu.get("FeO", 0.0), dtype=float) + fe2o3 = np.asarray(apfu.get("Fe2O3", 0.0), dtype=float) + fe_metal = np.asarray(apfu.get("Fe", 0.0), dtype=float) + atomic_o = np.asarray(apfu.get("O", 0.0), dtype=float) # Calculate total Fe as FeO equivalent (molar atoms) - if atomic_o > 0: + if np.any(atomic_o > 0): # MAGEMin O-basis (ig, mp) or sb24 basis - if fe_metal > 0: + if np.any(fe_metal > 0): fe = fe_metal else: fe = feo @@ -117,21 +123,37 @@ def _extract_cations_from_apfu(self, out, phase, cations, sys_in): fe = feo + 2.0 * fe2o3 total = mg + mn + fe + ca - if total <= 0: - return {c: 0.0 for c in cations} + if instance != 'all': + if total.ndim == 0 and total <= 0: + return {f"cat_{c}": 0.0 for c in cations} + total = np.maximum(total, 1e-12) + vals = {"Mg": mg/total, "Mn": mn/total, "Fe": fe/total, "Ca": ca/total} + if sys_in.casefold() == 'wt': + vals = atomic_frac_to_wt_frac(vals, atomic_mass_dict) + return {f"cat_{c}": float(vals.get(c, 0.0)) for c in cations} - vals = {"Mg": mg/total, "Mn": mn/total, "Fe": fe/total, "Ca": ca/total} - + n = len(total) + total_safe = np.where(total > 0, total, np.nan) + vals = {"Mg": mg/total_safe, "Mn": mn/total_safe, + "Fe": fe/total_safe, "Ca": ca/total_safe} if sys_in.casefold() == 'wt': - keys = list(vals.keys()) - raw_vals = [vals[k] for k in keys] - wt_percents = convert_mol_percent_to_wt_percent(raw_vals, keys, atomic_mass_dict) - vals = {k: v/100.0 for k, v in zip(keys, wt_percents)} + vals = atomic_frac_to_wt_frac(vals, atomic_mass_dict) + return {f"cat_{c}": np.asarray(vals.get(c, 0.0), dtype=float) for c in cations} - return {f"cat_{c}": vals.get(c, 0.0) for c in cations} + def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, grid_out=None, instance=0): + """Extract phase properties from a previously calculated grid. - def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, grid_out=None): - """Extract phase properties from a previously calculated grid.""" + Parameters + ---------- + instance : int or {'all'}, default=0 + For a phase that appears multiple times (a solvus), an integer + index selects that instance and ``'all'`` returns one column + per instance, keyed with a ``_0``, ``_1`` ... suffix (missing + instances at a given point are NaN). ``mol_frac``/``wt_frac``/ + ``vol_frac`` are always the summed totals; with ``'all'`` + additional per-instance columns ``mol_frac_0``, ``mol_frac_1`` + ... are emitted. + """ out = grid_out if grid_out is not None else self.last_grid_out if out is None: raise ValueError("No grid results found. Run calculate_grid first or provide grid_out.") @@ -144,54 +166,92 @@ def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None results = { "mol_frac": np.zeros(P_len), "wt_frac": np.zeros(P_len), "vol_frac": np.zeros(P_len), } - + + all_inst = instance == 'all' + if all_inst: + n_inst = 0 + for o in out: + n_inst = max(n_inst, len(_phase_indices(o, phase, 'all'))) + # per-instance fraction columns alongside the summed totals + for k in range(n_inst): + results[f"mol_frac_{k}"] = np.zeros(P_len) + results[f"wt_frac_{k}"] = np.zeros(P_len) + results[f"vol_frac_{k}"] = np.zeros(P_len) + else: + n_inst = 1 + + def _make_keys(prefix, names): + if all_inst: + return [f"{prefix}{n}_{k}" for n in names for k in range(n_inst)] + return [f"{prefix}{n}" for n in names] + if end_members: - for em in end_members: results[f"em_{em}"] = np.zeros(P_len) + for key in _make_keys("em_", end_members): results[key] = np.zeros(P_len) if oxides: - for ox in oxides: results[f"ox_apfu_{ox}"] = np.zeros(P_len) + for key in _make_keys("ox_apfu_", oxides): results[key] = np.zeros(P_len) if chemistry: - for ox in chemistry: results[f"chem_{ox}"] = np.zeros(P_len) + for key in _make_keys("chem_", chemistry): results[key] = np.zeros(P_len) if cations: - for c in cations: results[f"cat_{c}"] = np.zeros(P_len) + for key in _make_keys("cat_", cations): results[key] = np.zeros(P_len) if mg_number: - results["mg_number"] = np.zeros(P_len) + for key in _make_keys("Mg_number", [""]): results[key] = np.zeros(P_len) if fe_split: - results["fe2"] = np.zeros(P_len) - results["fe3"] = np.zeros(P_len) + for key in _make_keys("Fe2", [""]): results[key] = np.zeros(P_len) + for key in _make_keys("Fe3", [""]): results[key] = np.zeros(P_len) for i in range(P_len): - if phase in out[i].ph: - results["mol_frac"][i] = phase_frac(phase, out[i], 'mol') - results["wt_frac"][i] = phase_frac(phase, out[i], 'wt') - results["vol_frac"][i] = phase_frac(phase, out[i], 'vol') - - if end_members: - for em in end_members: - results[f"em_{em}"][i] = extract_end_member(phase, out[i], em, self.sys_in) - if oxides: - apfu = get_oxide_apfu(out[i], phase, oxides) - for ox in oxides: results[f"ox_apfu_{ox}"][i] = apfu.get(ox, 0.0) - if chemistry: - chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in) - for ox in chemistry: results[f"chem_{ox}"][i] = chem.get(ox, 0.0) - if cations: - cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in) - for c in cations: results[f"cat_{c}"][i] = cat_vals[f"cat_{c}"] - if mg_number: - results["mg_number"][i] = get_phase_mg_number(out[i], phase) - if fe_split: - split = self._extract_fe_split_from_apfu(out[i], phase) - results["fe2"][i] = split["fe2"] - results["fe3"][i] = split["fe3"] + if phase not in out[i].ph: + continue + results["mol_frac"][i] = phase_frac(phase, out[i], 'mol') + results["wt_frac"][i] = phase_frac(phase, out[i], 'wt') + results["vol_frac"][i] = phase_frac(phase, out[i], 'vol') + + if all_inst: + n_idx = _phase_indices(out[i], phase, 'all') + fm = [float(out[i].ph_frac[j]) for j in n_idx] + fw = [float(out[i].ph_frac_wt[j]) for j in n_idx] + fv = [float(out[i].ph_frac_vol[j]) for j in n_idx] + for k in range(n_inst): + results[f"mol_frac_{k}"][i] = fm[k] if k < len(fm) else np.nan + results[f"wt_frac_{k}"][i] = fw[k] if k < len(fw) else np.nan + results[f"vol_frac_{k}"][i] = fv[k] if k < len(fv) else np.nan + + def _store(key, value): + if all_inst: + vals = np.asarray(value, dtype=float) + n = len(vals) + for k in range(n_inst): + results[f"{key}_{k}"][i] = vals[k] if k < n else np.nan + else: + results[key][i] = value + + if end_members: + for em in end_members: + _store(f"em_{em}", extract_end_member(phase, out[i], em, self.sys_in, instance=instance)) + if oxides: + apfu = get_oxide_apfu(out[i], phase, oxides, instance=instance) + for ox in oxides: _store(f"ox_apfu_{ox}", apfu.get(ox, 0.0)) + if chemistry: + chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in, instance=instance) + for ox in chemistry: _store(f"chem_{ox}", chem.get(ox, 0.0)) + if cations: + cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in, instance=instance) + for c in cations: _store(f"cat_{c}", cat_vals[f"cat_{c}"]) + if mg_number: + _store("Mg_number", get_phase_mg_number(out[i], phase, instance=instance)) + if fe_split: + split = self._extract_fe_split_from_apfu(out[i], phase, instance=instance) + _store("Fe2", split["Fe2"]) + _store("Fe3", split["Fe3"]) return results - def generate_2D_grid(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False): + def generate_2D_grid(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, instance=0): """Convenience wrapper.""" self.calculate_grid(P, T) - return self.extract_from_grid(phase, end_members, oxides, chemistry, cations, mg_number, fe_split) + return self.extract_from_grid(phase, end_members, oxides, chemistry, cations, mg_number, fe_split, instance=instance) - def single_point_calc(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False): + def single_point_calc(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, instance=0): """Single-point calculation.""" out = MAGEMin_C.single_point_minimization(P, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list) sys.stdout.flush() @@ -203,20 +263,39 @@ def single_point_calc(self, P, T, phase, end_members=None, oxides=None, chemistr results["wt_frac"] = phase_frac(phase, out, 'wt') results["vol_frac"] = phase_frac(phase, out, 'vol') + all_inst = instance == 'all' + n_inst = len(_phase_indices(out, phase, 'all')) if all_inst else 1 + + if all_inst: + for k, j in enumerate(_phase_indices(out, phase, 'all')): + results[f"mol_frac_{k}"] = float(out.ph_frac[j]) + results[f"wt_frac_{k}"] = float(out.ph_frac_wt[j]) + results[f"vol_frac_{k}"] = float(out.ph_frac_vol[j]) + + def _store(key, value): + if all_inst: + vals = np.asarray(value, dtype=float) + n = len(vals) + for k in range(n_inst): + results[f"{key}_{k}"] = vals[k] if k < n else np.nan + else: + results[key] = value + if end_members: - for em in end_members: results[f"em_{em}"] = extract_end_member(phase, out, em, self.sys_in) + for em in end_members: + _store(f"em_{em}", extract_end_member(phase, out, em, self.sys_in, instance=instance)) if oxides: - apfu = get_oxide_apfu(out, phase, oxides) - for ox in oxides: results[f"ox_apfu_{ox}"] = apfu.get(ox, 0.0) + apfu = get_oxide_apfu(out, phase, oxides, instance=instance) + for ox in oxides: _store(f"ox_apfu_{ox}", apfu.get(ox, 0.0)) if chemistry: - chem = get_phase_chemistry(out, phase, chemistry, self.sys_in) - for ox in chemistry: results[f"chem_{ox}"] = chem.get(ox, 0.0) + chem = get_phase_chemistry(out, phase, chemistry, self.sys_in, instance=instance) + for ox in chemistry: _store(f"chem_{ox}", chem.get(ox, 0.0)) if cations: - cat_vals = self._extract_cations_from_apfu(out, phase, cations, self.sys_in) - for c in cations: results[f"cat_{c}"] = cat_vals[f"cat_{c}"] + cat_vals = self._extract_cations_from_apfu(out, phase, cations, self.sys_in, instance=instance) + for c in cations: _store(f"cat_{c}", cat_vals[f"cat_{c}"]) if fe_split: - split = self._extract_fe_split_from_apfu(out, phase) - results["fe2"] = split["fe2"] - results["fe3"] = split["fe3"] + split = self._extract_fe_split_from_apfu(out, phase, instance=instance) + _store("Fe2", split["Fe2"]) + _store("Fe3", split["Fe3"]) return results, out From 99abc1c260100f7cd3f12fa3a58213414801f6b0 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:40:35 +0800 Subject: [PATCH 05/10] refactor(calculators): return per-instance bundles from extract_from_grid --- src/phasetools/calculators/garnet.py | 12 +- src/phasetools/calculators/pt_grid.py | 241 +++++++++-------- src/phasetools/core/README.md | 45 +++- tests/test_solvus_instances.py | 365 ++++++++++++++++++++++++++ 4 files changed, 525 insertions(+), 138 deletions(-) create mode 100644 tests/test_solvus_instances.py diff --git a/src/phasetools/calculators/garnet.py b/src/phasetools/calculators/garnet.py index 4565236..7e594d3 100644 --- a/src/phasetools/calculators/garnet.py +++ b/src/phasetools/calculators/garnet.py @@ -51,12 +51,18 @@ def _extract_garnet_elements_from_oxides(self, out, sys_in): return Mg, Mn, Fe, Ca def generate_2D_grid_gt_endmembers(self, P, T): - """Compute garnet endmember fractions over a P-T grid.""" + """Compute garnet endmember fractions over a P-T grid. + + Garnet is a single-instance phase, so its single per-instance + bundle is returned directly, preserving the historical key names + (``em_py``, ``em_alm``, ``mol_frac``, ...). If garnet ever + appeared as a solvus, the list of bundles would be returned + instead. + """ self.calculate_grid(P, T) # Automatic discovery of end-members res = self.extract_from_grid("g", end_members='auto') - - return res + return res[0] if len(res) == 1 else res def generate_2D_grid_gt_elements(self, P, T): """Compute garnet element fractions (Mg, Mn, Fe, Ca) over a P-T grid.""" diff --git a/src/phasetools/calculators/pt_grid.py b/src/phasetools/calculators/pt_grid.py index 8c2fb74..7fee7ec 100644 --- a/src/phasetools/calculators/pt_grid.py +++ b/src/phasetools/calculators/pt_grid.py @@ -2,7 +2,7 @@ import sys from juliacall import Main as jl, convert as jlconvert from ..core.base import MAGEMinBase -from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, _phase_indices +from ..core.phase_properties import extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, _phase_indices from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac from phasetools import MAGEMin_C @@ -97,10 +97,13 @@ def get_phase_endmembers(self, phase, grid_out=None): return [] return [] - def _extract_cations_from_apfu(self, out, phase, cations, sys_in, instance=0): - """Internal: Extract cation ratios (e.g., XMg, XFe) for a specific phase.""" + def _extract_cations_from_apfu(self, out, phase, cations, sys_in): + """Internal: per-instance cation ratios (e.g., XMg, XFe). + + Returns ``cat_`` arrays, one entry per instance of the phase. + """ ox_to_query = ['MgO', 'MnO', 'CaO', 'FeO', 'Fe2O3', 'Fe', 'O'] - apfu = get_oxide_apfu(out, phase, ox_to_query, instance=instance) + apfu = get_oxide_apfu(out, phase, ox_to_query, instance='all') mg = np.asarray(apfu.get("MgO", 0.0), dtype=float) mn = np.asarray(apfu.get("MnO", 0.0), dtype=float) @@ -123,36 +126,44 @@ def _extract_cations_from_apfu(self, out, phase, cations, sys_in, instance=0): fe = feo + 2.0 * fe2o3 total = mg + mn + fe + ca - if instance != 'all': - if total.ndim == 0 and total <= 0: - return {f"cat_{c}": 0.0 for c in cations} - total = np.maximum(total, 1e-12) - vals = {"Mg": mg/total, "Mn": mn/total, "Fe": fe/total, "Ca": ca/total} - if sys_in.casefold() == 'wt': - vals = atomic_frac_to_wt_frac(vals, atomic_mass_dict) - return {f"cat_{c}": float(vals.get(c, 0.0)) for c in cations} - - n = len(total) total_safe = np.where(total > 0, total, np.nan) - vals = {"Mg": mg/total_safe, "Mn": mn/total_safe, - "Fe": fe/total_safe, "Ca": ca/total_safe} + with np.errstate(divide='ignore', invalid='ignore'): + vals = {"Mg": mg/total_safe, "Mn": mn/total_safe, + "Fe": fe/total_safe, "Ca": ca/total_safe} + if sys_in.casefold() == 'wt': vals = atomic_frac_to_wt_frac(vals, atomic_mass_dict) - return {f"cat_{c}": np.asarray(vals.get(c, 0.0), dtype=float) for c in cations} - def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, grid_out=None, instance=0): + return {f"cat_{c}": np.asarray(vals[c], dtype=float) for c in cations} + + def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, grid_out=None): """Extract phase properties from a previously calculated grid. - Parameters - ---------- - instance : int or {'all'}, default=0 - For a phase that appears multiple times (a solvus), an integer - index selects that instance and ``'all'`` returns one column - per instance, keyed with a ``_0``, ``_1`` ... suffix (missing - instances at a given point are NaN). ``mol_frac``/``wt_frac``/ - ``vol_frac`` are always the summed totals; with ``'all'`` - additional per-instance columns ``mol_frac_0``, ``mol_frac_1`` - ... are emitted. + MAGEMin reports coexisting solvus limbs as repeated entries in + ``out.ph`` (e.g. two clinopyroxenes ``dio``). Returns a list + with one bundle per instance, indexed ``res[0]``, ``res[1]`` + ...: + + ``res[0]`` -- first (or only) instance; every key is an array + over the grid points (scalars for ``single_point_calc``). + + The key shape is identical for every phase -- a phase with a + single instance simply has a one-element list, always at + ``res[0]``. Sum the per-instance ``mol_frac`` arrays yourself + if you want the total. + + Each bundle key is NaN where the phase (or that instance) is + absent at a grid point. If the phase never appears in the grid, + the returned list is empty. + + Notes + ----- + - Bundle ``res[k]`` is the k-th occurrence of the phase in that + point's ``out.ph``; solvus branch ordering may swap between + grid points, so track both limbs when plotting isopleths. + - ``end_members='auto'`` discovers end-members from the first + occurrence of the phase; solvus limbs normally share the same + solution model and endmember set. """ out = grid_out if grid_out is not None else self.last_grid_out if out is None: @@ -163,139 +174,125 @@ def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None end_members = self.get_phase_endmembers(phase, out) P_len = len(out) - results = { - "mol_frac": np.zeros(P_len), "wt_frac": np.zeros(P_len), "vol_frac": np.zeros(P_len), - } - - all_inst = instance == 'all' - if all_inst: - n_inst = 0 - for o in out: - n_inst = max(n_inst, len(_phase_indices(o, phase, 'all'))) - # per-instance fraction columns alongside the summed totals - for k in range(n_inst): - results[f"mol_frac_{k}"] = np.zeros(P_len) - results[f"wt_frac_{k}"] = np.zeros(P_len) - results[f"vol_frac_{k}"] = np.zeros(P_len) - else: - n_inst = 1 - def _make_keys(prefix, names): - if all_inst: - return [f"{prefix}{n}_{k}" for n in names for k in range(n_inst)] - return [f"{prefix}{n}" for n in names] + # n_inst = max occurrences of the phase across the whole grid + n_inst = 0 + for o in out: + n_inst = max(n_inst, len(_phase_indices(o, phase, 'all'))) + + instances = [{} for _ in range(n_inst)] + + def _precreate(prefix, names): + for k in range(n_inst): + for n in names: + instances[k][f"{prefix}{n}"] = np.full(P_len, np.nan) + for k in range(n_inst): + instances[k]["mol_frac"] = np.full(P_len, np.nan) + instances[k]["wt_frac"] = np.full(P_len, np.nan) + instances[k]["vol_frac"] = np.full(P_len, np.nan) if end_members: - for key in _make_keys("em_", end_members): results[key] = np.zeros(P_len) + _precreate("em_", end_members) if oxides: - for key in _make_keys("ox_apfu_", oxides): results[key] = np.zeros(P_len) + _precreate("ox_apfu_", oxides) if chemistry: - for key in _make_keys("chem_", chemistry): results[key] = np.zeros(P_len) + _precreate("chem_", chemistry) if cations: - for key in _make_keys("cat_", cations): results[key] = np.zeros(P_len) + _precreate("cat_", cations) if mg_number: - for key in _make_keys("Mg_number", [""]): results[key] = np.zeros(P_len) + for k in range(n_inst): + instances[k]["Mg_number"] = np.full(P_len, np.nan) if fe_split: - for key in _make_keys("Fe2", [""]): results[key] = np.zeros(P_len) - for key in _make_keys("Fe3", [""]): results[key] = np.zeros(P_len) + for k in range(n_inst): + instances[k]["Fe2"] = np.full(P_len, np.nan) + instances[k]["Fe3"] = np.full(P_len, np.nan) for i in range(P_len): if phase not in out[i].ph: continue - results["mol_frac"][i] = phase_frac(phase, out[i], 'mol') - results["wt_frac"][i] = phase_frac(phase, out[i], 'wt') - results["vol_frac"][i] = phase_frac(phase, out[i], 'vol') - - if all_inst: - n_idx = _phase_indices(out[i], phase, 'all') - fm = [float(out[i].ph_frac[j]) for j in n_idx] - fw = [float(out[i].ph_frac_wt[j]) for j in n_idx] - fv = [float(out[i].ph_frac_vol[j]) for j in n_idx] + n_idx = _phase_indices(out[i], phase, 'all') + for k, j in enumerate(n_idx): + instances[k]["mol_frac"][i] = float(out[i].ph_frac[j]) + instances[k]["wt_frac"][i] = float(out[i].ph_frac_wt[j]) + instances[k]["vol_frac"][i] = float(out[i].ph_frac_vol[j]) + # instances k >= len(n_idx) stay NaN + + def _fill(key, value): + vals = np.atleast_1d(np.asarray(value, dtype=float)) + n = len(vals) for k in range(n_inst): - results[f"mol_frac_{k}"][i] = fm[k] if k < len(fm) else np.nan - results[f"wt_frac_{k}"][i] = fw[k] if k < len(fw) else np.nan - results[f"vol_frac_{k}"][i] = fv[k] if k < len(fv) else np.nan - - def _store(key, value): - if all_inst: - vals = np.asarray(value, dtype=float) - n = len(vals) - for k in range(n_inst): - results[f"{key}_{k}"][i] = vals[k] if k < n else np.nan - else: - results[key][i] = value + if k < n: + instances[k][key][i] = vals[k] + # else stays NaN if end_members: for em in end_members: - _store(f"em_{em}", extract_end_member(phase, out[i], em, self.sys_in, instance=instance)) + _fill(f"em_{em}", extract_end_member(phase, out[i], em, self.sys_in, instance='all')) if oxides: - apfu = get_oxide_apfu(out[i], phase, oxides, instance=instance) - for ox in oxides: _store(f"ox_apfu_{ox}", apfu.get(ox, 0.0)) + apfu = get_oxide_apfu(out[i], phase, oxides, instance='all') + for ox in oxides: _fill(f"ox_apfu_{ox}", apfu.get(ox, np.zeros(0))) if chemistry: - chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in, instance=instance) - for ox in chemistry: _store(f"chem_{ox}", chem.get(ox, 0.0)) + chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in, instance='all') + for ox in chemistry: _fill(f"chem_{ox}", chem.get(ox, np.zeros(0))) if cations: - cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in, instance=instance) - for c in cations: _store(f"cat_{c}", cat_vals[f"cat_{c}"]) + cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in) + for c in cations: _fill(f"cat_{c}", cat_vals[f"cat_{c}"]) if mg_number: - _store("Mg_number", get_phase_mg_number(out[i], phase, instance=instance)) + _fill("Mg_number", get_phase_mg_number(out[i], phase, instance='all')) if fe_split: - split = self._extract_fe_split_from_apfu(out[i], phase, instance=instance) - _store("Fe2", split["Fe2"]) - _store("Fe3", split["Fe3"]) + split = self._extract_fe_split_from_apfu(out[i], phase, instance='all') + _fill("Fe2", split["Fe2"]) + _fill("Fe3", split["Fe3"]) - return results + return instances - def generate_2D_grid(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, instance=0): + def generate_2D_grid(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False): """Convenience wrapper.""" self.calculate_grid(P, T) - return self.extract_from_grid(phase, end_members, oxides, chemistry, cations, mg_number, fe_split, instance=instance) + return self.extract_from_grid(phase, end_members, oxides, chemistry, cations, mg_number, fe_split) + + def single_point_calc(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False): + """Single-point calculation. - def single_point_calc(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, instance=0): - """Single-point calculation.""" + Returns ``(bundles, out)`` where ``bundles`` is a list with one + bundle per instance of the phase (``bundles[0]``, ``bundles[1]``, + ...), keyed like the grid-level ``extract_from_grid`` with scalar + values. The list is empty when the phase is not present at this + P-T. + """ out = MAGEMin_C.single_point_minimization(P, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list) sys.stdout.flush() - results = {"mol_frac": 0.0, "wt_frac": 0.0, "vol_frac": 0.0, "present": False} + instances = [] if phase in out.ph: - results["present"] = True - results["mol_frac"] = phase_frac(phase, out, 'mol') - results["wt_frac"] = phase_frac(phase, out, 'wt') - results["vol_frac"] = phase_frac(phase, out, 'vol') - - all_inst = instance == 'all' - n_inst = len(_phase_indices(out, phase, 'all')) if all_inst else 1 - - if all_inst: - for k, j in enumerate(_phase_indices(out, phase, 'all')): - results[f"mol_frac_{k}"] = float(out.ph_frac[j]) - results[f"wt_frac_{k}"] = float(out.ph_frac_wt[j]) - results[f"vol_frac_{k}"] = float(out.ph_frac_vol[j]) + n_idx = _phase_indices(out, phase, 'all') + instances = [{"mol_frac": 0.0, "wt_frac": 0.0, "vol_frac": 0.0} + for _ in n_idx] + for k, j in enumerate(n_idx): + instances[k]["mol_frac"] = float(out.ph_frac[j]) + instances[k]["wt_frac"] = float(out.ph_frac_wt[j]) + instances[k]["vol_frac"] = float(out.ph_frac_vol[j]) def _store(key, value): - if all_inst: - vals = np.asarray(value, dtype=float) - n = len(vals) - for k in range(n_inst): - results[f"{key}_{k}"] = vals[k] if k < n else np.nan - else: - results[key] = value + vals = np.atleast_1d(np.asarray(value, dtype=float)) + for k in range(len(vals)): + instances[k][key] = vals[k] if end_members: for em in end_members: - _store(f"em_{em}", extract_end_member(phase, out, em, self.sys_in, instance=instance)) + _store(f"em_{em}", extract_end_member(phase, out, em, self.sys_in, instance='all')) if oxides: - apfu = get_oxide_apfu(out, phase, oxides, instance=instance) - for ox in oxides: _store(f"ox_apfu_{ox}", apfu.get(ox, 0.0)) + apfu = get_oxide_apfu(out, phase, oxides, instance='all') + for ox in oxides: _store(f"ox_apfu_{ox}", apfu.get(ox, np.zeros(0))) if chemistry: - chem = get_phase_chemistry(out, phase, chemistry, self.sys_in, instance=instance) - for ox in chemistry: _store(f"chem_{ox}", chem.get(ox, 0.0)) + chem = get_phase_chemistry(out, phase, chemistry, self.sys_in, instance='all') + for ox in chemistry: _store(f"chem_{ox}", chem.get(ox, np.zeros(0))) if cations: - cat_vals = self._extract_cations_from_apfu(out, phase, cations, self.sys_in, instance=instance) + cat_vals = self._extract_cations_from_apfu(out, phase, cations, self.sys_in) for c in cations: _store(f"cat_{c}", cat_vals[f"cat_{c}"]) if fe_split: - split = self._extract_fe_split_from_apfu(out, phase, instance=instance) + split = self._extract_fe_split_from_apfu(out, phase, instance='all') _store("Fe2", split["Fe2"]) _store("Fe3", split["Fe3"]) - - return results, out + + return instances, out diff --git a/src/phasetools/core/README.md b/src/phasetools/core/README.md index 6e1b37d..2aa41f3 100644 --- a/src/phasetools/core/README.md +++ b/src/phasetools/core/README.md @@ -26,20 +26,39 @@ The `core` submodule provides the foundational classes and low-level bridging lo ### Solvus (multi-instance) handling MAGEMin reports coexisting solvus limbs as repeated entries in `out.ph` (e.g. two `dio` clinopyroxenes or two `amp` amphiboles). `phase_frac` -sums them; the per-phase composition helpers take an `instance` argument: +sums them; the low-level composition helpers in `phase_properties.py` +(`get_oxide_apfu`, `get_phase_chemistry`, `extract_end_member`, +`get_phase_mg_number`, `get_phase_mg2_number`, `get_phase_fe_split`) +take an `instance` argument: - integer index (default `0`) — that instance, with a `UserWarning` - noting how many other instances exist (`_phase_indices` raises - `ValueError` for non-integer, non-`'all'` values); -- `'all'` — one value per instance, returned as numpy arrays - (`get_oxide_apfu` → `{oxide: ndarray}`, `extract_end_member`, - `get_phase_mg_number`, `get_phase_mg2_number`, `get_phase_fe_split`). + noting how many other instances exist; out-of-range → zeros; +- `'all'` — one value per instance, returned as numpy arrays. `MAGEMinPTGridCalculator.extract_from_grid` / `single_point_calc` / -`generate_2D_grid` accept the same `instance` argument. With -`instance='all'` they emit one column per instance keyed with a -`_0`, `_1`, ... suffix (e.g. `ox_apfu_Na2O_0`); grid points where the -phase has fewer instances are NaN-padded. `mol_frac`/`wt_frac`/ -`vol_frac` always give the summed totals (matching `phase_frac`); with -`instance='all'` additional per-instance columns `mol_frac_0`, -`mol_frac_1`, ... are emitted alongside them. +`generate_2D_grid` return a **list of per-instance bundles**, indexed +`res[0]`, `res[1]`, ... : + +```python +res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) +c0, c1 = res # limb 0, limb 1 +c0['mol_frac'] # array over grid points (limb 0) +c0['ox_apfu_Na2O'] # array over grid points (limb 0) +total_dio = c0['mol_frac'] + c1['mol_frac'] # sum the limbs yourself +``` + +- Bundle keys are unsuffixed (`mol_frac`, `wt_frac`, `vol_frac`, + `ox_apfu_*`, `chem_*`, `em_*`, `cat_*`, `Mg_number`, `Fe2`, `Fe3`). +- The shape is the **same for every phase** — a single-instance phase + (e.g. garnet) has a one-element list, always at `res[0]`. +- Each bundle value is an array over the grid points, **NaN** where the + phase (or that limb) is absent; if the phase never appears, the list + is empty. No totals are computed — sum the `mol_frac` bundles if + you want the whole-phase fraction. +- Bundle `res[k]` is the k-th occurrence of the phase in that point's + `out.ph`; solvus branch ordering may swap between grid points, so + track both limbs when plotting isopleths. `end_members='auto'` + discovers end-members from the first occurrence (solvus limbs share + the same solution model). `MAGEMinGarnetCalculator. + generate_2D_grid_gt_endmembers` returns the single-instance bundle + directly, preserving the historical `em_py`, `mol_frac`, ... keys. diff --git a/tests/test_solvus_instances.py b/tests/test_solvus_instances.py new file mode 100644 index 0000000..ac35947 --- /dev/null +++ b/tests/test_solvus_instances.py @@ -0,0 +1,365 @@ +"""Mock-based tests for solvus (multi-instance phase) handling. + +MAGEMin reports coexisting solvus limbs as repeated entries in ``out.ph`` +(e.g. two clinopyroxenes ``dio``, or two amphiboles ``amp``). + +Two layers are tested: + +* the low-level composition helpers (``get_oxide_apfu``, + ``extract_end_member``, ...) take ``instance`` (integer index or + ``'all'`` for per-instance numpy arrays); +* the grid-level ``extract_from_grid`` / ``single_point_calc`` return a + list of per-instance bundles indexed ``res[0]``, ``res[1]`` ... Bundle + keys are unsuffixed (``mol_frac``, ``ox_apfu_Na2O``, ``em_py``, ...), + and each value is an array over the grid points (NaN where the phase + or that limb is absent). The shape is the same for every phase; a + single-instance phase just has a one-element list. + +No live Julia runtime is needed -- ``out`` objects are mocked. +""" + +import unittest +import warnings +import numpy as np +from unittest.mock import MagicMock, patch + +from phasetools.calculators.pt_grid import MAGEMinPTGridCalculator +from phasetools.core.phase_properties import ( + get_oxide_apfu, get_phase_chemistry, extract_end_member, + get_phase_mg_number, get_phase_mg2_number, get_phase_fe_split, + _phase_indices, +) + +# mpe-style oxide ordering +OXIDES = ['H2O', 'SiO2', 'Al2O3', 'CaO', 'MgO', 'FeO', 'K2O', 'Na2O', + 'TiO2', 'MnO', 'O'] + + +def _make_ss(apfu, em_names, em_frac): + """Build a mocked SS_vec entry (one phase instance).""" + p = MagicMock() + p.Comp_apfu = apfu + # molar fractions (0-1) aligned with OXIDES + p.Comp = [0.0, 0.6, 0.03, 0.15, 0.13, 0.02, 0.0, 0.06, 0.0, 0.01, 0.0] + p.Comp_wt = [0.0, 0.55, 0.05, 0.12, 0.10, 0.02, 0.0, 0.05, 0.0, 0.01, 0.0] + p.emNames = em_names + p.emFrac = em_frac + p.emFrac_wt = [f * 0.9 for f in em_frac] + return p + + +def _dio_0(): + # diopside-rich limb: low Na + apfu = [0.0, 1.98, 0.05, 0.74, 0.64, 0.12, 0.0, 0.26, 0.0, 0.01, 6.0] + return _make_ss(apfu, ['di', 'hed', 'om', 'jac'], [0.60, 0.10, 0.25, 0.05]) + + +def _dio_1(): + # omphacite-rich limb: high Na + apfu = [0.0, 1.96, 0.09, 0.59, 0.51, 0.18, 0.0, 0.41, 0.0, 0.01, 6.0] + return _make_ss(apfu, ['di', 'hed', 'om', 'jac'], [0.40, 0.10, 0.40, 0.10]) + + +def _garnet(): + apfu = [0.0, 3.02, 2.0, 0.55, 0.70, 1.70, 0.0, 0.0, 0.0, 0.03, 12.0] + return _make_ss(apfu, ['py', 'alm', 'gr', 'spss'], [0.20, 0.60, 0.16, 0.04]) + + +def _make_out(ph_list, ss_vec, ph_frac=None): + out = MagicMock() + out.ph = ph_list + out.oxides = OXIDES + n = len(ph_list) + if ph_frac is None: + ph_frac = [0.1] * n + out.ph_frac = ph_frac + out.ph_frac_wt = list(ph_frac) + out.ph_frac_vol = list(ph_frac) + out.SS_vec = ss_vec + return out + + +class TestPhaseIndices(unittest.TestCase): + """_phase_indices resolution.""" + + def setUp(self): + self.out = _make_out(['dio', 'q', 'dio', 'g'], [_dio_0(), None, _dio_1(), _garnet()]) + + def test_all_returns_every_occurrence(self): + self.assertEqual(_phase_indices(self.out, 'dio', 'all'), [0, 2]) + + def test_default_zero_warns_with_others(self): + with self.assertWarnsRegex(UserWarning, r"2 instances.*1.*others"): + self.assertEqual(_phase_indices(self.out, 'dio', 0), [0]) + + def test_single_instance_does_not_warn(self): + with warnings.catch_warnings(): + warnings.simplefilter("error") + self.assertEqual(_phase_indices(self.out, 'g', 0), [3]) + + def test_out_of_range_warns_and_empty(self): + with self.assertWarnsRegex(UserWarning, r"instance 2 does not exist"): + self.assertEqual(_phase_indices(self.out, 'dio', 2), []) + + def test_absent_phase_empty(self): + self.assertEqual(_phase_indices(self.out, 'ep', 0), []) + + def test_invalid_instance_type_raises(self): + with self.assertRaises(ValueError): + _phase_indices(self.out, 'dio', 'first') + + +class TestSolvusHelpers(unittest.TestCase): + """Per-instance behaviour of the composition helpers.""" + + def setUp(self): + self.out = _make_out(['dio', 'q', 'dio', 'g'], [_dio_0(), None, _dio_1(), _garnet()]) + + def test_oxide_apfu_default_is_first(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO']) + self.assertAlmostEqual(apfu['Na2O'], 0.26) + self.assertAlmostEqual(apfu['MgO'], 0.64) + + def test_oxide_apfu_second_instance(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO'], instance=1) + self.assertAlmostEqual(apfu['Na2O'], 0.41) + self.assertAlmostEqual(apfu['MgO'], 0.51) + + def test_oxide_apfu_all_returns_arrays(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO'], instance='all') + self.assertTrue(np.allclose(apfu['Na2O'], [0.26, 0.41])) + self.assertTrue(np.allclose(apfu['MgO'], [0.64, 0.51])) + + def test_phase_chemistry_second_instance(self): + chem = get_phase_chemistry(self.out, 'dio', ['Na2O'], 'mol', instance=1) + self.assertAlmostEqual(chem['Na2O'], 0.06 * 100.0) # Comp[7] * 100 + + def test_extract_end_member_per_instance(self): + self.assertAlmostEqual(extract_end_member('dio', self.out, 'di', 'mol'), 0.60) + self.assertAlmostEqual(extract_end_member('dio', self.out, 'di', 'mol', instance=1), 0.40) + vals = extract_end_member('dio', self.out, 'di', 'mol', instance='all') + self.assertTrue(np.allclose(vals, [0.60, 0.40])) + + def test_mg_number_per_instance(self): + # Mg# = Mg/(Mg+Fe); first limb Mg-rich, second more Fe-rich + self.assertGreater(get_phase_mg_number(self.out, 'dio'), + get_phase_mg_number(self.out, 'dio', instance=1)) + vals = get_phase_mg_number(self.out, 'dio', instance='all') + self.assertEqual(vals.shape, (2,)) + self.assertAlmostEqual(vals[0], get_phase_mg_number(self.out, 'dio')) + + def test_mg2_number_per_instance(self): + m0 = get_phase_mg2_number(self.out, 'dio') + m1 = get_phase_mg2_number(self.out, 'dio', instance=1) + self.assertGreater(m0, m1) + vals = get_phase_mg2_number(self.out, 'dio', instance='all') + self.assertEqual(vals.shape, (2,)) + + def test_fe_split_per_instance(self): + # O-basis: excess O = 0, so Fe2 = total Fe, Fe3 = 0 for both limbs + s0 = get_phase_fe_split(self.out, 'dio') + s1 = get_phase_fe_split(self.out, 'dio', instance=1) + self.assertAlmostEqual(s0['Fe2'], 0.12) + self.assertAlmostEqual(s1['Fe2'], 0.18) + all_s = get_phase_fe_split(self.out, 'dio', instance='all') + self.assertTrue(np.allclose(all_s['Fe2'], [0.12, 0.18])) + + def test_out_of_range_returns_zeros(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O'], instance=5) + self.assertEqual(apfu['Na2O'], 0.0) + + +class TestExtractFromGridSolvus(unittest.TestCase): + """Grid extraction: list of per-instance bundles (res[0], res[1], ...).""" + + def _make_calc(self, grid_out): + calc = MAGEMinPTGridCalculator.__new__(MAGEMinPTGridCalculator) + calc.sys_in = 'mol' + calc.last_grid_out = grid_out + calc.rm_list = None + return calc + + def _dio_grid(self): + # point 0: dio solvus (2 limbs); point 1: single dio + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + g1 = _make_out(['q', 'g', 'dio'], + [None, _garnet(), _dio_0()], + ph_frac=[0.5, 0.2, 0.3]) + return [g0, g1] + + def test_uniform_shape_single_and_solvus(self): + """A single-instance phase and a solvus return the same bundle keys.""" + calc = self._make_calc(self._dio_grid()) + solvus = calc.extract_from_grid('dio', oxides=['Na2O', 'MgO'], + mg_number=True, fe_split=True) + single = calc.extract_from_grid('g', oxides=['MgO'], + mg_number=True, fe_split=True) + + # dio: two bundles; garnet: one bundle (always at res[0]) + self.assertEqual(len(solvus), 2) + self.assertEqual(len(single), 1) + + for bundle in solvus + single: + for key in ('mol_frac', 'wt_frac', 'vol_frac', + 'ox_apfu_Na2O' if 'ox_apfu_Na2O' in bundle else 'ox_apfu_MgO', + 'Mg_number', 'Fe2', 'Fe3'): + self.assertIn(key, bundle) + # unsuffixed keys inside each bundle (no _0/_1, no total_*) + self.assertIn('mol_frac', solvus[0]) + self.assertNotIn('mol_frac_0', solvus[0]) + self.assertNotIn('total_mol_frac', solvus[0]) + + def test_limb_bundles(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', oxides=['Na2O']) + c0, c1 = res + # limb fractions (0.3/0.1 at point 0; point 1 has only limb 0) + self.assertTrue(np.allclose(c0['mol_frac'], [0.3, 0.3])) + self.assertAlmostEqual(c0['wt_frac'][0], 0.3) + self.assertAlmostEqual(c0['vol_frac'][0], 0.3) + self.assertAlmostEqual(c1['mol_frac'][0], 0.1) + self.assertAlmostEqual(c0['ox_apfu_Na2O'][0], 0.26) + self.assertAlmostEqual(c1['ox_apfu_Na2O'][0], 0.41) + # user can sum the limbs themselves + total = c0['mol_frac'] + c1['mol_frac'] + self.assertAlmostEqual(total[0], 0.4) + + def test_nan_when_limb_absent(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) + c0, c1 = res + # point 1 has a single dio -> limb-1 values NaN + self.assertTrue(np.isnan(c1['mol_frac'][1])) + self.assertTrue(np.isnan(c1['wt_frac'][1])) + self.assertTrue(np.isnan(c1['ox_apfu_Na2O'][1])) + self.assertTrue(np.isnan(c1['Mg_number'][1])) + # limb 0 is filled + self.assertAlmostEqual(c0['mol_frac'][1], 0.3) + self.assertAlmostEqual(c0['ox_apfu_Na2O'][1], 0.26) + + def test_absent_phase_bundles_nan(self): + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + g2 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc([g0, g2]) + res = calc.extract_from_grid('dio', oxides=['Na2O']) + c0, c1 = res + # point 1: dio absent -> NaN in every bundle + self.assertTrue(np.isnan(c0['mol_frac'][1])) + self.assertTrue(np.isnan(c1['mol_frac'][1])) + self.assertTrue(np.isnan(c0['ox_apfu_Na2O'][1])) + + def test_absent_everywhere_empty(self): + """Phase never stable -> empty list.""" + g0 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + g1 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.5, 0.5]) + calc = self._make_calc([g0, g1]) + res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) + self.assertEqual(res, []) + + def test_single_instance_phase_filled(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('g', oxides=['MgO']) + bundle = res[0] + self.assertAlmostEqual(bundle['mol_frac'][0], 0.4) + self.assertAlmostEqual(bundle['ox_apfu_MgO'][0], 0.70) + + def test_cations_per_instance(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', cations=['Mg', 'Fe']) + c0, c1 = res + self.assertIn('cat_Mg', c0) + self.assertIn('cat_Mg', c1) + # first limb is Mg-richer + self.assertGreater(c0['cat_Mg'][0], c1['cat_Mg'][0]) + + def test_end_members_per_instance(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', end_members=['di', 'om']) + c0, c1 = res + self.assertAlmostEqual(c0['em_di'][0], 0.60) + self.assertAlmostEqual(c1['em_di'][0], 0.40) + self.assertAlmostEqual(c0['em_om'][0], 0.25) + self.assertAlmostEqual(c1['em_om'][0], 0.40) + + +class TestSinglePointCalcSolvus(unittest.TestCase): + """single_point_calc returns the same nested-bundle schema.""" + + def _make_calc(self): + calc = MAGEMinPTGridCalculator.__new__(MAGEMinPTGridCalculator) + calc.sys_in = 'mol' + calc.data = None + calc.X = None + calc.Xoxides = None + calc.rm_list = None + return calc + + def test_solvus_returns_limb_bundles(self): + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g0 + res, out = calc.single_point_calc(10.0, 600.0, 'dio', + oxides=['Na2O', 'MgO']) + self.assertEqual(len(res), 2) + c0, c1 = res + self.assertAlmostEqual(c0['mol_frac'], 0.3) + self.assertAlmostEqual(c1['mol_frac'], 0.1) + self.assertAlmostEqual(c0['wt_frac'], 0.3) + self.assertAlmostEqual(c1['wt_frac'], 0.1) + self.assertAlmostEqual(c0['ox_apfu_Na2O'], 0.26) + self.assertAlmostEqual(c1['ox_apfu_Na2O'], 0.41) + self.assertIs(out, g0) + + def test_single_instance_has_one_bundle(self): + g = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g + res, _ = calc.single_point_calc(10.0, 600.0, 'g', oxides=['MgO']) + self.assertEqual(len(res), 1) + bundle = res[0] + self.assertAlmostEqual(bundle['mol_frac'], 0.4) + self.assertAlmostEqual(bundle['ox_apfu_MgO'], 0.70) + + def test_absent_phase_no_bundles(self): + g = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g + res, _ = calc.single_point_calc(10.0, 600.0, 'dio') + self.assertEqual(res, []) + + +class TestGarnetEndmembersSuffix(unittest.TestCase): + """generate_2D_grid_gt_endmembers returns the single-instance bundle.""" + + def test_returns_bundle_with_historical_keys(self): + from phasetools.calculators.garnet import MAGEMinGarnetCalculator + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + bundle = { + 'mol_frac': np.array([0.10, 0.15]), + 'wt_frac': np.array([0.11, 0.16]), + 'vol_frac': np.array([0.12, 0.17]), + 'em_py': np.array([0.20, 0.30]), + 'em_alm': np.array([0.60, 0.50]), + 'em_gr': np.array([0.16, 0.17]), + 'em_spss': np.array([0.04, 0.03]), + } + with patch.object(calc, 'calculate_grid', return_value=None), \ + patch.object(calc, 'extract_from_grid', return_value=[bundle]): + res = calc.generate_2D_grid_gt_endmembers([10.0], [600.0]) + self.assertIn('em_py', res) + self.assertIn('em_alm', res) + self.assertIn('mol_frac', res) + self.assertTrue(np.allclose(res['em_py'], [0.20, 0.30])) + + +if __name__ == '__main__': + unittest.main() From fe9996ae594a05e746307267483877e441903c95 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:17:05 +0800 Subject: [PATCH 06/10] fix(core): instance-aware type hints, element-wise Fe split, safe scalar conversion --- src/phasetools/core/phase_properties.py | 33 ++++++++++++------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/src/phasetools/core/phase_properties.py b/src/phasetools/core/phase_properties.py index 02a0e1f..4a78141 100644 --- a/src/phasetools/core/phase_properties.py +++ b/src/phasetools/core/phase_properties.py @@ -2,7 +2,7 @@ import warnings -def _phase_indices(out, phase, instance=0): +def _phase_indices(out: object, phase: str, instance: int | str = 0) -> list[int]: """Return the index/indices of ``phase`` in ``out.ph``. Parameters @@ -246,12 +246,12 @@ def get_phase_fe_split(out, ph, instance=0): return {"Fe2": 0.0, "Fe3": 0.0} if instance != 'all' else empty # 1. Calculate Total Fe atoms (Atoms per formula unit) - if np.any(ato > 0): - # MAGEMin O-basis (ig, mp) or sb24 basis - total_fe = fem if np.any(fem > 0) else feo - else: - # Traditional FeO/Fe2O3 basis - total_fe = feo + 2.0 * fe2o3 + # Use np.where so each instance follows its own basis (O-bearing vs traditional). + total_fe = np.where( + ato > 0, + np.where(fem > 0, fem, feo), + feo + 2.0 * fe2o3 + ) # 2. Calculate Fe3+ atoms using excess oxygen heuristic excess_o = np.maximum(ato - np.floor(ato), 0.0) @@ -261,10 +261,10 @@ def get_phase_fe_split(out, ph, instance=0): fe2 = np.maximum(np.asarray(total_fe) - fe3, 0.0) if instance != 'all': - return {"Fe2": float(np.squeeze(fe2)), "Fe3": float(np.squeeze(fe3))} + return {"Fe2": float(fe2.item()), "Fe3": float(fe3.item())} return {"Fe2": fe2, "Fe3": fe3} -def get_phase_mg2_number(out, ph, instance=0): +def get_phase_mg2_number(out: object, ph: str, instance: int | str = 0) -> float | np.ndarray: """ Calculate Mg# (molar Mg / (Mg + Fe2+)) for a specific phase. @@ -273,16 +273,15 @@ def get_phase_mg2_number(out, ph, instance=0): try: apfu = get_oxide_apfu(out, ph, ['MgO'], instance=instance) mg = apfu.get('MgO', 0.0) - + split = get_phase_fe_split(out, ph, instance=instance) fe2 = split['Fe2'] - - denominator = mg + fe2 - if np.any(np.asarray(denominator) <= 0): - return 0.0 if instance != 'all' else np.zeros(len(mg)) - - return (mg / denominator) if instance != 'all' else np.asarray(mg / denominator, dtype=float) - except: + + denominator = np.asarray(mg + fe2, dtype=float) + if instance != 'all': + return float(mg / denominator) if denominator > 0 else 0.0 + return np.where(denominator > 0, mg / denominator, 0.0) + except Exception: return 0.0 def calculate_kd_fe_mg(out, phase1, phase2, use_fe2_only=False): From efde58981e4eaace514ed1752083b45d8213c2b8 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:40:43 +0800 Subject: [PATCH 07/10] feat(calculators): garnet fe_basis, wt fractionation, instance-aware extraction --- src/phasetools/calculators/garnet.py | 144 +++++++++++++++++--- src/phasetools/calculators/phase_search.py | 62 ++++++++- src/phasetools/calculators/pt_estimation.py | 4 +- 3 files changed, 181 insertions(+), 29 deletions(-) diff --git a/src/phasetools/calculators/garnet.py b/src/phasetools/calculators/garnet.py index 7e594d3..5236ab6 100644 --- a/src/phasetools/calculators/garnet.py +++ b/src/phasetools/calculators/garnet.py @@ -2,51 +2,91 @@ import sys from .pt_grid import MAGEMinPTGridCalculator from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu -from ..utils.bulk_rock import atomic_mass_dict, convert_mol_percent_to_wt_percent +from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac from phasetools import MAGEMin_C from juliacall import Main as jl, convert as jlconvert class MAGEMinGarnetCalculator(MAGEMinPTGridCalculator): """High-level wrappers for garnet-focused MAGEMin calculations.""" - def __init__(self, db="ig", dataset=636, verbose=False): + + def __init__(self, db="ig", dataset=636, verbose=False, fe_basis="FeOt"): + """ + Parameters + ---------- + db : str, default="ig" + Thermodynamic database label. + dataset : int, default=636 + Thermodynamic dataset version. + verbose : bool, default=False + Whether to print progress information. + fe_basis : str, default="FeOt" + Iron basis used for garnet X-site fractions: + + * ``'FeOt'`` -- all iron treated as Fe2+ (total iron, + ``FeO + 2*Fe2O3`` APFU) placed on the divalent site. This is + the standard community convention for garnet end-members and + X-site fractions (e.g. Williams & Grambling 1990; Krogh Ravna + 2000) and is consistent with the pyralspite garnet model used + by MAGEMin and Holland-Powell-type databases. + * ``'Fe2+'`` -- only the stoichiometrically estimated ferrous + iron (``_extract_fe_split_from_apfu``) is placed on the + divalent site, excluding Fe3+. Use only when garnet Fe3+ is + known to be significant (oxidised eclogites, skarns) or + measured directly (XANES, Mössbauer). + + Case-insensitive. + """ super().__init__(db, dataset, verbose) + basis = str(fe_basis).casefold() + if basis not in ("feot", "fe2+", "fe2"): + raise ValueError(f"fe_basis must be 'FeOt' or 'Fe2+', got {fe_basis!r}") + self.fe_basis = "fe2+" if basis.startswith("fe2") else "feot" def _extract_garnet_elements_from_oxides(self, out, sys_in): """ Extract garnet Mg-Mn-Fe-Ca cation fractions for the divalent (X) site. - Strictly isolates Fe2+ to ensure X-site fractions sum to 1.0. + + The Fe basis is controlled by ``self.fe_basis``: + + * ``'feot'`` (default) -- all iron treated as Fe2+ (total iron, + ``FeO + 2*Fe2O3`` APFU) placed on the divalent site. The standard + community convention for garnet X-site fractions. + * ``'fe2+'`` -- only the stoichiometrically estimated ferrous iron + is placed on the divalent site, excluding Fe3+. """ if 'g' not in out.ph: return 0.0, 0.0, 0.0, 0.0 - # Isolating divalent iron using the robust Fe-split method - split = self._extract_fe_split_from_apfu(out, 'g') - fe2_moles = split['fe2'] - - elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO']) + if self.fe_basis == 'feot': + # Total iron on the divalent site (all Fe as Fe2+) + elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO', 'FeO', 'Fe2O3']) + fe_moles = elements.get("FeO", 0.0) + 2.0 * elements.get("Fe2O3", 0.0) + else: + # Ferrous iron only, from the excess-oxygen Fe2+/Fe3+ split + elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO']) + fe_moles = self._extract_fe_split_from_apfu(out, 'g')['Fe2'] + mg_moles = elements.get("MgO", 0.0) mn_moles = elements.get("MnO", 0.0) ca_moles = elements.get("CaO", 0.0) # Total atoms in the X-site (should be approx 3.0) - total_x_site_moles = mg_moles + mn_moles + fe2_moles + ca_moles - + total_x_site_moles = mg_moles + mn_moles + fe_moles + ca_moles if total_x_site_moles <= 0: return 0.0, 0.0, 0.0, 0.0 # Normalise to 1.0 (mole fractions of the divalent site) Mg = mg_moles / total_x_site_moles Mn = mn_moles / total_x_site_moles - Fe = fe2_moles / total_x_site_moles + Fe = fe_moles / total_x_site_moles Ca = ca_moles / total_x_site_moles if sys_in.casefold() == 'wt': - wt_percent_list = convert_mol_percent_to_wt_percent( - [Mg, Mn, Fe, Ca], - ["Mg", "Mn", "Fe", "Ca"], + wt_frac = atomic_frac_to_wt_frac( + {"Mg": Mg, "Mn": Mn, "Fe": Fe, "Ca": Ca}, atomic_mass_dict, ) - Mg, Mn, Fe, Ca = [val / 100.0 for val in wt_percent_list] + Mg, Mn, Fe, Ca = wt_frac["Mg"], wt_frac["Mn"], wt_frac["Fe"], wt_frac["Ca"] return Mg, Mn, Fe, Ca @@ -127,7 +167,61 @@ def _gt_single_point_from_jl(self, P, T, X_jl, Xoxides_jl, sys_in, rm_list=None) return gt_frac, gt_wt, gt_vol, Mg, Mn, Fe, Ca, out def gt_along_path(self, P, T, fractionate=False, normalise_start=True): - """Calculate garnet fractions and element chemistry along a P-T path.""" + """Calculate garnet fractions and element chemistry along a P-T path. + + Parameters + ---------- + P : array-like + Pressure values along the path (kbar). + T : array-like + Temperature values along the path (°C). + fractionate : bool, default=False + If True, fractionate garnet from the bulk composition as it grows. + normalise_start : bool, default=True + Controls how the first P-T point is treated: + + * ``True`` — The first P-T point is treated as a nucleation + barrier: garnet volume starts at zero and only new growth is + modelled. The initial garnet fraction (if any) is used as a + baseline for measuring incremental growth but is **not** removed + from the reactive bulk. Use when the path starts outside the + garnet stability field or at the nucleation threshold. + + * ``False`` — The first P-T point has an initial garnet volume + (overstepped nucleation). That fraction is removed from the + bulk at step 0, and subsequent growth is measured relative to + it. Use when the path starts well inside the garnet stability + field. + + Returns + ------- + gt_mol_frac : numpy.ndarray + Garnet molar fraction at each P-T point. + gt_wt_frac : numpy.ndarray + Garnet weight fraction at each P-T point. + gt_vol_frac : numpy.ndarray + Garnet volume fraction at each P-T point. + Mgi : numpy.ndarray + Garnet X-site Mg fraction at each P-T point. + Mni : numpy.ndarray + Garnet X-site Mn fraction at each P-T point. + Fei : numpy.ndarray + Garnet X-site Fe fraction at each P-T point. The Fe basis + (total Fe as Fe2+ by default, or ferrous-only) is set by the + ``fe_basis`` argument passed to :meth:`__init__`. + Cai : numpy.ndarray + Garnet X-site Ca fraction at each P-T point. + X_along_path : numpy.ndarray + Bulk composition after each step's fractionation. Each row is + normalised to sum to 1. Row ``i`` is the bulk **after** step + ``i``'s fractionation has been applied. + + Notes + ----- + Fractionation uses the **current-step** garnet composition (not the + growth-increment composition) — a first-order approximation valid for + small P-T steps. + """ from .phase_search import PhaseFunctions X = self.X @@ -144,6 +238,7 @@ def gt_along_path(self, P, T, fractionate=False, normalise_start=True): X_along_path = np.zeros(shape=(n_points, len(self._Xoxides_py)) ) gt_frac_max_previous = 0. + gt_wt_max_previous = 0. phase_functions = PhaseFunctions(db=self.db, dataset=self.dataset, verbose=self.verbose) if fractionate else None if phase_functions: # Sync standardised state to the helper instance @@ -158,19 +253,28 @@ def gt_along_path(self, P, T, fractionate=False, normalise_start=True): gt_wt_frac[i] = gt_wt gt_vol_frac[i] = gt_vol + # Select the fraction basis consistent with sys_in + if self.sys_in.casefold() == 'wt': + gt_frac_for_fractionation = gt_wt + gt_frac_max_prev_for_fractionation = gt_wt_max_previous + else: + gt_frac_for_fractionation = gt_frac + gt_frac_max_prev_for_fractionation = gt_frac_max_previous + if phase_functions is not None: if i == 0 and not normalise_start: - if gt_frac > 0: - X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=gt_frac) + if gt_frac_for_fractionation > 0: + X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=gt_frac_for_fractionation) X = jlconvert(jl.Vector[jl.Float64], X_py) elif i > 0: - frac_amount = max(gt_frac - gt_frac_max_previous, 0.0) + frac_amount = max(gt_frac_for_fractionation - gt_frac_max_prev_for_fractionation, 0.0) if frac_amount > 0: X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=frac_amount) X = jlconvert(jl.Vector[jl.Float64], X_py) - X_along_path[i] = np.array(X) + X_along_path[i] = np.array(X) / np.sum(X) gt_frac_max_previous = max(gt_frac_max_previous, gt_frac) + gt_wt_max_previous = max(gt_wt_max_previous, gt_wt) Mgi[i] = Mg Mni[i] = Mn diff --git a/src/phasetools/calculators/phase_search.py b/src/phasetools/calculators/phase_search.py index 54fe8a0..cdfa8d9 100644 --- a/src/phasetools/calculators/phase_search.py +++ b/src/phasetools/calculators/phase_search.py @@ -41,24 +41,72 @@ def liquidus_func(T): return result.root def fractionate_phase(self, phase, out, sys_in, frac_amount=None): - """Perform batch fractionation of a phase from the bulk rock composition.""" + """Perform batch fractionation of a phase from the bulk rock composition. + + Removes a specified fraction of a phase from the current bulk composition + and renormalises the remainder. The output is always normalised to sum to 1 + (molar or weight fractions, depending on ``sys_in``). + + MAGEMin orders ``out.ph`` with solution phases first (``SS_vec``), then + pure phases (``PP_vec``). This method handles both indexing ranges + transparently. + + .. note:: + + When a phase appears multiple times in ``out.ph`` (e.g. two + coexisting pyroxenes on a solvus), only the **first** instance is + removed. Use ``phase_frac()`` to obtain the summed fraction before + calling this method if the total solvus fraction is needed. + + Parameters + ---------- + phase : str + Phase name as it appears in ``out.ph`` (e.g. ``'g'``, ``'liq'``). + out : MAGEMinOutput + Raw output from a MAGEMin single-point minimisation. + sys_in : str + ``'mol'`` or ``'wt'`` — determines which bulk composition and + phase composition vectors are used. + frac_amount : float or None, optional + Fraction of the phase to remove. If ``None``, the full phase + fraction from the minimisation result is used. Values of 0 or + less are treated as a no-op and return the current bulk unchanged. + + Returns + ------- + numpy.ndarray + Renormalised bulk composition after fractionation. + """ if sys_in.casefold() == "wt": current_X = out.bulk_wt else: current_X = out.bulk - + if phase in out.ph: phase_ind = out.ph.index(phase) - if sys_in.casefold() == "wt": - ph_comp = np.array(out.SS_vec[phase_ind].Comp_wt) - ph_frac = out.ph_frac_wt[phase_ind] + # MAGEMin orders: solution phases (SS_vec) first, then pure phases (PP_vec) + if phase_ind < out.n_SS: + if sys_in.casefold() == "wt": + ph_comp = np.array(out.SS_vec[phase_ind].Comp_wt) + ph_frac = out.ph_frac_wt[phase_ind] + else: + ph_comp = np.array(out.SS_vec[phase_ind].Comp) + ph_frac = out.ph_frac[phase_ind] else: - ph_comp = np.array(out.SS_vec[phase_ind].Comp) - ph_frac = out.ph_frac[phase_ind] + pp_ind = phase_ind - out.n_SS + if sys_in.casefold() == "wt": + ph_comp = np.array(out.PP_vec[pp_ind].Comp_wt) + ph_frac = out.ph_frac_wt[phase_ind] + else: + ph_comp = np.array(out.PP_vec[pp_ind].Comp) + ph_frac = out.ph_frac[phase_ind] if frac_amount is None: frac_amount = ph_frac + if frac_amount <= 0: + return np.array(current_X, dtype=float) + if frac_amount >= 1.0: warnings.warn(f"fractionate_phase: requested frac_amount={frac_amount} >= 1.0; skipping.") return np.array(current_X, dtype=float) diff --git a/src/phasetools/calculators/pt_estimation.py b/src/phasetools/calculators/pt_estimation.py index ca7fe64..88227d7 100644 --- a/src/phasetools/calculators/pt_estimation.py +++ b/src/phasetools/calculators/pt_estimation.py @@ -26,8 +26,8 @@ def _get_phase_composition(self, out, phase, components, comp_type='element'): # Handle special components Fe2 and Fe3 first if 'Fe2' in components or 'Fe3' in components: split = self._extract_fe_split_from_apfu(out, phase) - element_map['Fe2'] = split['fe2'] - element_map['Fe3'] = split['fe3'] + element_map['Fe2'] = split['Fe2'] + element_map['Fe3'] = split['Fe3'] for ox, stoichiometry in self._stoich_map.items(): for el, mult in stoichiometry.items(): From 01f03de19c0d867e630956151034146bef39e6e7 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:21:45 +0800 Subject: [PATCH 08/10] fix(calculators): add type hints to garnet __init__, document fe2 alias --- src/phasetools/calculators/garnet.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/phasetools/calculators/garnet.py b/src/phasetools/calculators/garnet.py index 5236ab6..8d2a30c 100644 --- a/src/phasetools/calculators/garnet.py +++ b/src/phasetools/calculators/garnet.py @@ -9,7 +9,9 @@ class MAGEMinGarnetCalculator(MAGEMinPTGridCalculator): """High-level wrappers for garnet-focused MAGEMin calculations.""" - def __init__(self, db="ig", dataset=636, verbose=False, fe_basis="FeOt"): + def __init__( + self, db: str = "ig", dataset: int = 636, verbose: bool = False, fe_basis: str = "FeOt" + ) -> None: """ Parameters ---------- @@ -33,6 +35,7 @@ def __init__(self, db="ig", dataset=636, verbose=False, fe_basis="FeOt"): divalent site, excluding Fe3+. Use only when garnet Fe3+ is known to be significant (oxidised eclogites, skarns) or measured directly (XANES, Mössbauer). + * ``'Fe2'`` or ``'fe2'`` -- alias for ``'Fe2+'``. Case-insensitive. """ From b115c3c2e1f0e7d40424f7cf939c7d40a38742cf Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:41:01 +0800 Subject: [PATCH 09/10] test: add Fe-basis and fractionation tests --- tests/test_fe_basis.py | 100 ++++++++++ tests/test_fractionation.py | 363 ++++++++++++++++++++++++++++++++++++ 2 files changed, 463 insertions(+) create mode 100644 tests/test_fe_basis.py create mode 100644 tests/test_fractionation.py diff --git a/tests/test_fe_basis.py b/tests/test_fe_basis.py new file mode 100644 index 0000000..1565e5c --- /dev/null +++ b/tests/test_fe_basis.py @@ -0,0 +1,100 @@ +"""Mock-based tests for the garnet ``fe_basis`` option. + +Verifies that ``MAGEMinGarnetCalculator._extract_garnet_elements_from_oxides`` +honours the ``fe_basis`` flag: + +* ``'FeOt'`` (default) -- total Fe (FeO + 2*Fe2O3 APFU) placed on the + divalent X-site. This is the standard community convention for garnet + end-members / X-site fractions. +* ``'Fe2+'`` -- only the stoichiometrically estimated ferrous iron is + placed on the divalent site, excluding Fe3+. + +No live Julia runtime is needed -- ``get_oxide_apfu`` and the Fe2+/Fe3+ +split are mocked. +""" + +import unittest +from unittest.mock import MagicMock, patch + +from phasetools.calculators.garnet import MAGEMinGarnetCalculator + + +class TestGarnetFeBasis(unittest.TestCase): + """Fe-basis flag on the garnet X-site extraction.""" + + def _make_calc(self, fe_basis): + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + calc.fe_basis = fe_basis # canonical lowercase form set by __init__ + return calc + + def test_feot_uses_total_iron(self): + """FeOt puts FeO + 2*Fe2O3 on the X-site and normalises to 1.""" + calc = self._make_calc('feot') + out = MagicMock() + out.ph = ['g', 'q'] + apfu = {'MgO': 0.6, 'MnO': 0.1, 'CaO': 0.8, 'FeO': 1.2, 'Fe2O3': 0.15} + with patch('phasetools.calculators.garnet.get_oxide_apfu', return_value=apfu): + Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol') + fe_total = 1.2 + 2.0 * 0.15 + total = 0.6 + 0.1 + 0.8 + fe_total + self.assertAlmostEqual(Fe, fe_total / total, places=6) + self.assertAlmostEqual(Mg, 0.6 / total, places=6) + self.assertAlmostEqual(Mn, 0.1 / total, places=6) + self.assertAlmostEqual(Ca, 0.8 / total, places=6) + self.assertAlmostEqual(Mg + Mn + Fe + Ca, 1.0, places=6) + + def test_fe2_uses_split(self): + """Fe2+ uses only the ferrous split, excluding Fe3+.""" + calc = self._make_calc('fe2+') + out = MagicMock() + out.ph = ['g', 'q'] + apfu = {'MgO': 0.6, 'MnO': 0.1, 'CaO': 0.8} + split = {'Fe2': 1.35, 'Fe3': 0.15} + with patch('phasetools.calculators.garnet.get_oxide_apfu', return_value=apfu), \ + patch.object(calc, '_extract_fe_split_from_apfu', return_value=split): + Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol') + total = 0.6 + 0.1 + 0.8 + 1.35 + self.assertAlmostEqual(Fe, 1.35 / total, places=6) + self.assertAlmostEqual(Mg + Mn + Fe + Ca, 1.0, places=6) + self.assertAlmostEqual(Mn, 0.1 / total, places=6) + + def test_absent_garnet_returns_zeros(self): + """No garnet in the assemblage -> all-zero X-site fractions.""" + calc = self._make_calc('feot') + out = MagicMock() + out.ph = ['q', 'dio'] + Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol') + self.assertEqual((Mg, Mn, Fe, Ca), (0.0, 0.0, 0.0, 0.0)) + + def test_invalid_basis_raises(self): + """Unsupported fe_basis values must raise ValueError at construction.""" + with patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__', + return_value=None): + with self.assertRaises(ValueError): + MAGEMinGarnetCalculator(db='ig', fe_basis='Fe3') + + def test_default_is_feot(self): + """The default fe_basis is 'FeOt' (community convention).""" + with patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__', + return_value=None): + calc = MAGEMinGarnetCalculator(db='ig') + self.assertEqual(calc.fe_basis, 'feot') + calc = MAGEMinGarnetCalculator(db='ig', fe_basis='Fe2+') + self.assertEqual(calc.fe_basis, 'fe2+') + calc = MAGEMinGarnetCalculator(db='ig', fe_basis='feot') + self.assertEqual(calc.fe_basis, 'feot') + + +class TestGarnetGeneratorFeBasis(unittest.TestCase): + """GarnetGenerator forwards fe_basis to the garnet calculator.""" + + @patch('phasetools.models.garnet_growth.MAGEMinGarnetCalculator.__init__') + def test_forwards_fe_basis(self, mock_init): + from phasetools.models.garnet_growth import GarnetGenerator + GarnetGenerator(db='mpe', fe_basis='Fe2+') + _, kwargs = mock_init.call_args + self.assertEqual(kwargs.get('fe_basis'), 'Fe2+') + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_fractionation.py b/tests/test_fractionation.py new file mode 100644 index 0000000..a67bc8d --- /dev/null +++ b/tests/test_fractionation.py @@ -0,0 +1,363 @@ +"""Mock-based tests for fractionation correctness fixes. + +These tests verify the fixes described in the implementation plan: +- Fix A: get_retrograde_concentrations crash (garnet_growth.py) +- Fix B: mol/wt unit mismatch in gt_along_path (garnet.py) +- Fix C: pure-phase IndexError in fractionate_phase (phase_search.py) +- Fix D: X_along_path normalise to 1 (garnet.py) +- Fix H1: self.X permanent mutation in run_fractional_stages (magma_ocean.py) + +No live Julia runtime is needed — all MAGEMin calls are mocked. +""" + +import unittest +import numpy as np +from unittest.mock import MagicMock, patch, PropertyMock + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _make_mock_out(phases, n_SS, bulk_mol, bulk_wt, ph_frac_mol, ph_frac_wt, comps_mol, comps_wt): + """Build a minimal mock MAGEMin output object.""" + out = MagicMock() + out.ph = phases + out.n_SS = n_SS + out.bulk = np.array(bulk_mol, dtype=float) + out.bulk_wt = np.array(bulk_wt, dtype=float) + out.ph_frac = list(ph_frac_mol) + out.ph_frac_wt = list(ph_frac_wt) + + ss_vec = [] + pp_vec = [] + for idx, ph in enumerate(phases): + obj = MagicMock() + obj.Comp = np.array(comps_mol[idx], dtype=float) + obj.Comp_wt = np.array(comps_wt[idx], dtype=float) + if idx < n_SS: + ss_vec.append(obj) + else: + pp_vec.append(obj) + + out.SS_vec = ss_vec + out.PP_vec = pp_vec + return out + + +# =========================================================================== +# Fix A: get_retrograde_concentrations crash +# =========================================================================== +class TestRetrogradeNoCrash(unittest.TestCase): + """Fix A: get_retrograde_concentrations must not raise TypeError.""" + + @patch('phasetools.models.garnet_growth.MAGEMinGarnetCalculator.__init__', return_value=None) + def test_retrograde_no_crash(self, mock_init): + """Verify the fixed call to _gt_single_point_from_jl works.""" + from phasetools.models.garnet_growth import GarnetGenerator + + gen = GarnetGenerator.__new__(GarnetGenerator) + # Set up minimal state + gen.Pi = np.array([10.0, 11.0, 12.0, 13.0]) + gen.Ti = np.array([800.0, 810.0, 820.0, 830.0]) + gen.ti = np.array([0.0, 1.0, 2.0, 3.0]) + gen.gt_vol_frac = np.array([0.01, 0.02, 0.05, 0.05]) + gen.Mgi = np.array([0.3, 0.3, 0.3, 0.3]) + gen.Mni = np.array([0.05, 0.05, 0.05, 0.05]) + gen.Fei = np.array([0.4, 0.4, 0.4, 0.4]) + gen.Cai = np.array([0.25, 0.25, 0.25, 0.25]) + # X_along_path: last_growth_index=2 has different bulk than [-1] + gen.X_along_path = np.array([ + [50.0, 50.0], + [48.0, 52.0], + [45.0, 55.0], + [40.0, 60.0], + ], dtype=float) + gen.last_growth_index = 2 + gen.Xoxides = ['SiO2', 'Al2O3'] + gen.sys_in = 'mol' + gen.rm_list = None + + # Mock _gt_single_point_from_jl to return controlled output + mock_return = (0.05, 0.05, 0.05, 0.3, 0.05, 0.4, 0.25, MagicMock()) + gen._gt_single_point_from_jl = MagicMock(return_value=mock_return) + + # Should not raise TypeError + result = gen.get_retrograde_concentrations() + self.assertEqual(result.shape[0], 7) # t, T, P, Mn, Mg, Fe, Ca + self.assertTrue(gen._gt_single_point_from_jl.called) + + +# =========================================================================== +# Fix B: wt fraction uses wt basis +# =========================================================================== +class TestWtFractionUsesWtBasis(unittest.TestCase): + """Fix B: when sys_in='wt', fractionate_phase should use gt_wt, not gt_frac.""" + + @patch('phasetools.calculators.garnet.MAGEMin_C') + @patch('phasetools.calculators.garnet.jlconvert') + @patch('phasetools.calculators.phase_search.PhaseFunctions.__init__', return_value=None) + @patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__', return_value=None) + def test_wt_fraction_uses_wt_basis(self, mock_grid_init, mock_pf_init, mock_jlconvert, mock_magemin): + from phasetools.calculators.garnet import MAGEMinGarnetCalculator + + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + calc.db = 'ig' + calc.dataset = 636 + calc.verbose = False + calc.sys_in = 'wt' + calc.X = np.array([50.0, 50.0]) + calc.Xoxides = MagicMock() + calc._Xoxides_py = ['SiO2', 'Al2O3'] + calc.rm_list = None + calc.data = MagicMock() + + # gt_frac (mol) = 0.10, gt_wt = 0.05 — they differ + # Step 0 returns mol=0.10, wt=0.05; step 1 returns mol=0.12, wt=0.07 + mock_return_0 = (0.10, 0.05, 0.08, 0.3, 0.05, 0.4, 0.25, MagicMock()) + mock_return_1 = (0.12, 0.07, 0.10, 0.3, 0.05, 0.4, 0.25, MagicMock()) + calc._gt_single_point_from_jl = MagicMock(side_effect=[mock_return_0, mock_return_1]) + + # Mock PhaseFunctions + mock_pf = MagicMock() + mock_pf.fractionate_phase = MagicMock(return_value=np.array([50.0, 50.0])) + + # Patch PhaseFunctions at its source module (local import in gt_along_path) + with patch('phasetools.calculators.phase_search.PhaseFunctions', return_value=mock_pf): + calc._copy_state_to = MagicMock() + # jlconvert should pass through the array so np.array(X) works + mock_jlconvert.side_effect = lambda t, v: np.array(v, dtype=float) + + P = np.array([10.0, 11.0]) + T = np.array([800.0, 810.0]) + calc.gt_along_path(P, T, fractionate=True, normalise_start=True) + + # Step 0: normalise_start=True, so no fractionation at i=0 + # Step 1: i>0, frac_amount = gt_wt[1] - gt_wt[0] = 0.07 - 0.05 = 0.02 + calls = mock_pf.fractionate_phase.call_args_list + self.assertEqual(len(calls), 1, f"Expected 1 fractionation call, got {len(calls)}") + _, kwargs = calls[0] + self.assertAlmostEqual(kwargs['frac_amount'], 0.02, places=10, + msg="frac_amount should be based on wt fraction (0.07-0.05), not mol (0.12-0.10)") + + +# =========================================================================== +# Fix D: X_along_path normalised to 1 +# =========================================================================== +class TestXAlongPathNormalised(unittest.TestCase): + """Fix D: each row of X_along_path must sum to 1.0.""" + + @patch('phasetools.calculators.garnet.MAGEMin_C') + @patch('phasetools.calculators.garnet.jlconvert') + @patch('phasetools.calculators.phase_search.PhaseFunctions.__init__', return_value=None) + @patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__', return_value=None) + def test_x_along_path_normalised_to_one(self, mock_grid_init, mock_pf_init, mock_jlconvert, mock_magemin): + from phasetools.calculators.garnet import MAGEMinGarnetCalculator + + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + calc.db = 'ig' + calc.dataset = 636 + calc.verbose = False + calc.sys_in = 'mol' + calc.X = np.array([50.0, 50.0]) # Julia vector — jlconvert will wrap + calc.Xoxides = MagicMock() + calc._Xoxides_py = ['SiO2', 'Al2O3'] + calc.rm_list = None + calc.data = MagicMock() + + mock_return = (0.05, 0.05, 0.05, 0.3, 0.05, 0.4, 0.25, MagicMock()) + calc._gt_single_point_from_jl = MagicMock(return_value=mock_return) + mock_jlconvert.return_value = calc.X + + P = np.array([10.0, 11.0]) + T = np.array([800.0, 810.0]) + + _, _, _, _, _, _, _, X_along_path = calc.gt_along_path(P, T, fractionate=False) + + for i in range(len(P)): + self.assertAlmostEqual(np.sum(X_along_path[i]), 1.0, places=10, + msg=f"Row {i} of X_along_path does not sum to 1.0") + + +# =========================================================================== +# Fix C: pure-phase IndexError +# =========================================================================== +class TestFractionatePurePhase(unittest.TestCase): + """Fix C: fractionate_phase must handle pure phases (PP_vec) without IndexError.""" + + def test_fractionate_pure_phase(self): + from phasetools.calculators.phase_search import PhaseFunctions + + pf = PhaseFunctions.__new__(PhaseFunctions) + + # ph=['q', 'liq'], n_SS=1 => 'q' is a pure phase at index 0 in PP_vec + out = _make_mock_out( + phases=['q', 'liq'], + n_SS=1, + bulk_mol=[60.0, 40.0], + bulk_wt=[62.0, 38.0], + ph_frac_mol=[0.15, 0.85], + ph_frac_wt=[0.16, 0.84], + comps_mol=[[100.0, 0.0], [50.0, 50.0]], + comps_wt=[[100.0, 0.0], [48.0, 52.0]], + ) + + # 'q' is at index 0 in out.ph, n_SS=1, so it's a pure phase (PP_vec[0]) + result = pf.fractionate_phase('q', out, 'mol', frac_amount=0.1) + self.assertIsNotNone(result) + self.assertTrue(np.all(np.isfinite(result))) + + def test_fractionate_solution_phase_still_works(self): + """Verify solution-phase path (SS_vec) still works after the fix.""" + from phasetools.calculators.phase_search import PhaseFunctions + + pf = PhaseFunctions.__new__(PhaseFunctions) + + out = _make_mock_out( + phases=['liq', 'g'], + n_SS=2, + bulk_mol=[60.0, 40.0], + bulk_wt=[62.0, 38.0], + ph_frac_mol=[0.85, 0.15], + ph_frac_wt=[0.84, 0.16], + comps_mol=[[50.0, 50.0], [40.0, 60.0]], + comps_wt=[[48.0, 52.0], [38.0, 62.0]], + ) + + result = pf.fractionate_phase('g', out, 'mol', frac_amount=0.1) + self.assertIsNotNone(result) + self.assertTrue(np.all(np.isfinite(result))) + + +# =========================================================================== +# Zero guard: frac_amount=0 is a no-op +# =========================================================================== +class TestFractionateZeroIsNoop(unittest.TestCase): + """frac_amount=0 must return the bulk unchanged (not normalised to sum=1).""" + + def test_fractionate_zero_is_noop(self): + from phasetools.calculators.phase_search import PhaseFunctions + + pf = PhaseFunctions.__new__(PhaseFunctions) + + out = _make_mock_out( + phases=['g', 'liq'], + n_SS=2, + bulk_mol=[6000.0, 4000.0], # sum=10000, not 1 + bulk_wt=[6200.0, 3800.0], + ph_frac_mol=[0.15, 0.85], + ph_frac_wt=[0.16, 0.84], + comps_mol=[[40.0, 60.0], [50.0, 50.0]], + comps_wt=[[38.0, 62.0], [48.0, 52.0]], + ) + + result = pf.fractionate_phase('g', out, 'mol', frac_amount=0.0) + np.testing.assert_array_equal(result, np.array(out.bulk, dtype=float)) + + +# =========================================================================== +# Fix H1: self.X restored after run_fractional_stages +# =========================================================================== +class TestMagmaOceanXRestored(unittest.TestCase): + """Fix H1: run_fractional_stages must restore self.X after execution.""" + + @patch('phasetools.models.magma_ocean.MAGEMinBase.__init__', return_value=None) + def test_magma_ocean_x_restored(self, mock_base_init): + from phasetools.models.magma_ocean import MagmaOcean + import phasetools.models.magma_ocean as mo_module + + mo = MagmaOcean.__new__(MagmaOcean) + mo._Xoxides_py = ['SiO2', 'Al2O3'] + mo.sys_in = 'mol' + mo.data = MagicMock() + mo.Xoxides = MagicMock() + mo.rm_list = None + mo.X = np.array([50.0, 50.0]) + mo.radius_body = 1737.1 + mo.radius_core = 330.0 + mo.g = 1.62 + mo.rho_avg = 3350.0 + + saved_X = mo.X.copy() + + # Mock find_temperature_at_vol_frac to return a fixed T + mo.find_temperature_at_vol_frac = MagicMock(return_value=1200.0) + + # Build mock MAGEMin output + mock_out = MagicMock() + mock_out.ph = ['ol', 'liq'] + mock_out.n_SS = 2 + mock_out.ph_frac_vol = [0.3, 0.7] + + mock_ol = MagicMock() + mock_ol.rho = 3300.0 + mock_ol.Comp = np.array([30.0, 10.0]) + mock_ol.Comp_wt = np.array([28.0, 12.0]) + + mock_liq = MagicMock() + mock_liq.rho = 2800.0 + mock_liq.Comp = np.array([45.0, 55.0]) + mock_liq.Comp_wt = np.array([43.0, 57.0]) + + mock_out.SS_vec = [mock_ol, mock_liq] + mock_out.PP_vec = [] + + mock_magemin_c = MagicMock() + mock_magemin_c.single_point_minimization = MagicMock(return_value=mock_out) + + # Patch MAGEMin_C and jlconvert at the module level + with patch.object(mo_module, 'MAGEMin_C', mock_magemin_c), \ + patch.object(mo_module, 'jlconvert', side_effect=lambda t, v: np.array(v, dtype=float)): + mo.get_phase_chemistry_at_index = MagicMock(return_value=np.array([45.0, 55.0])) + mo.get_volume_between_radii = MagicMock(return_value=1e12) + mo.pressure_to_radius = MagicMock(return_value=1400.0) + mo.radius_to_pressure = MagicMock(return_value=5.0) + + starting_melt = np.array([45.0, 55.0]) + mo.run_fractional_stages( + starting_melt=starting_melt, + p_start=5.0, + p_end=0.001, + vol_step=0.05, + starting_vol_frac=0.5, + n_stages=2, + ) + + np.testing.assert_array_equal( + mo.X, saved_X, + err_msg="self.X was not restored after run_fractional_stages" + ) + + +# =========================================================================== +# Fix H4: starting_melt length validation +# =========================================================================== +class TestStartingMeltValidation(unittest.TestCase): + """Fix H4: run_fractional_stages must reject mismatched starting_melt length.""" + + @patch('phasetools.models.magma_ocean.MAGEMinBase.__init__', return_value=None) + def test_starting_melt_length_mismatch(self, mock_base_init): + from phasetools.models.magma_ocean import MagmaOcean + + mo = MagmaOcean.__new__(MagmaOcean) + mo._Xoxides_py = ['SiO2', 'Al2O3', 'MgO'] + mo.sys_in = 'mol' + mo.data = MagicMock() + mo.X = np.array([33.0, 33.0, 34.0]) + mo.radius_body = 1737.1 + mo.radius_core = 330.0 + mo.g = 1.62 + mo.rho_avg = 3350.0 + + # starting_melt has 2 elements, but _Xoxides_py has 3 + with self.assertRaises(ValueError) as ctx: + mo.run_fractional_stages( + starting_melt=np.array([50.0, 50.0]), + p_start=5.0, + p_end=0.001, + ) + self.assertIn("does not match", str(ctx.exception)) + + +if __name__ == '__main__': + unittest.main() From 2abd7dbc48681eab10dabb7adb28396ceedc8463 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:44:34 +0800 Subject: [PATCH 10/10] test: remove MagmaOcean and garnet_growth-dependent tests from garnet-fe-basis sub-PR --- tests/test_fe_basis.py | 11 --- tests/test_fractionation.py | 147 ------------------------------------ 2 files changed, 158 deletions(-) diff --git a/tests/test_fe_basis.py b/tests/test_fe_basis.py index 1565e5c..125a2a5 100644 --- a/tests/test_fe_basis.py +++ b/tests/test_fe_basis.py @@ -85,16 +85,5 @@ def test_default_is_feot(self): self.assertEqual(calc.fe_basis, 'feot') -class TestGarnetGeneratorFeBasis(unittest.TestCase): - """GarnetGenerator forwards fe_basis to the garnet calculator.""" - - @patch('phasetools.models.garnet_growth.MAGEMinGarnetCalculator.__init__') - def test_forwards_fe_basis(self, mock_init): - from phasetools.models.garnet_growth import GarnetGenerator - GarnetGenerator(db='mpe', fe_basis='Fe2+') - _, kwargs = mock_init.call_args - self.assertEqual(kwargs.get('fe_basis'), 'Fe2+') - - if __name__ == '__main__': unittest.main() diff --git a/tests/test_fractionation.py b/tests/test_fractionation.py index a67bc8d..297f9ff 100644 --- a/tests/test_fractionation.py +++ b/tests/test_fractionation.py @@ -45,49 +45,6 @@ def _make_mock_out(phases, n_SS, bulk_mol, bulk_wt, ph_frac_mol, ph_frac_wt, com return out -# =========================================================================== -# Fix A: get_retrograde_concentrations crash -# =========================================================================== -class TestRetrogradeNoCrash(unittest.TestCase): - """Fix A: get_retrograde_concentrations must not raise TypeError.""" - - @patch('phasetools.models.garnet_growth.MAGEMinGarnetCalculator.__init__', return_value=None) - def test_retrograde_no_crash(self, mock_init): - """Verify the fixed call to _gt_single_point_from_jl works.""" - from phasetools.models.garnet_growth import GarnetGenerator - - gen = GarnetGenerator.__new__(GarnetGenerator) - # Set up minimal state - gen.Pi = np.array([10.0, 11.0, 12.0, 13.0]) - gen.Ti = np.array([800.0, 810.0, 820.0, 830.0]) - gen.ti = np.array([0.0, 1.0, 2.0, 3.0]) - gen.gt_vol_frac = np.array([0.01, 0.02, 0.05, 0.05]) - gen.Mgi = np.array([0.3, 0.3, 0.3, 0.3]) - gen.Mni = np.array([0.05, 0.05, 0.05, 0.05]) - gen.Fei = np.array([0.4, 0.4, 0.4, 0.4]) - gen.Cai = np.array([0.25, 0.25, 0.25, 0.25]) - # X_along_path: last_growth_index=2 has different bulk than [-1] - gen.X_along_path = np.array([ - [50.0, 50.0], - [48.0, 52.0], - [45.0, 55.0], - [40.0, 60.0], - ], dtype=float) - gen.last_growth_index = 2 - gen.Xoxides = ['SiO2', 'Al2O3'] - gen.sys_in = 'mol' - gen.rm_list = None - - # Mock _gt_single_point_from_jl to return controlled output - mock_return = (0.05, 0.05, 0.05, 0.3, 0.05, 0.4, 0.25, MagicMock()) - gen._gt_single_point_from_jl = MagicMock(return_value=mock_return) - - # Should not raise TypeError - result = gen.get_retrograde_concentrations() - self.assertEqual(result.shape[0], 7) # t, T, P, Mn, Mg, Fe, Ca - self.assertTrue(gen._gt_single_point_from_jl.called) - - # =========================================================================== # Fix B: wt fraction uses wt basis # =========================================================================== @@ -255,109 +212,5 @@ def test_fractionate_zero_is_noop(self): np.testing.assert_array_equal(result, np.array(out.bulk, dtype=float)) -# =========================================================================== -# Fix H1: self.X restored after run_fractional_stages -# =========================================================================== -class TestMagmaOceanXRestored(unittest.TestCase): - """Fix H1: run_fractional_stages must restore self.X after execution.""" - - @patch('phasetools.models.magma_ocean.MAGEMinBase.__init__', return_value=None) - def test_magma_ocean_x_restored(self, mock_base_init): - from phasetools.models.magma_ocean import MagmaOcean - import phasetools.models.magma_ocean as mo_module - - mo = MagmaOcean.__new__(MagmaOcean) - mo._Xoxides_py = ['SiO2', 'Al2O3'] - mo.sys_in = 'mol' - mo.data = MagicMock() - mo.Xoxides = MagicMock() - mo.rm_list = None - mo.X = np.array([50.0, 50.0]) - mo.radius_body = 1737.1 - mo.radius_core = 330.0 - mo.g = 1.62 - mo.rho_avg = 3350.0 - - saved_X = mo.X.copy() - - # Mock find_temperature_at_vol_frac to return a fixed T - mo.find_temperature_at_vol_frac = MagicMock(return_value=1200.0) - - # Build mock MAGEMin output - mock_out = MagicMock() - mock_out.ph = ['ol', 'liq'] - mock_out.n_SS = 2 - mock_out.ph_frac_vol = [0.3, 0.7] - - mock_ol = MagicMock() - mock_ol.rho = 3300.0 - mock_ol.Comp = np.array([30.0, 10.0]) - mock_ol.Comp_wt = np.array([28.0, 12.0]) - - mock_liq = MagicMock() - mock_liq.rho = 2800.0 - mock_liq.Comp = np.array([45.0, 55.0]) - mock_liq.Comp_wt = np.array([43.0, 57.0]) - - mock_out.SS_vec = [mock_ol, mock_liq] - mock_out.PP_vec = [] - - mock_magemin_c = MagicMock() - mock_magemin_c.single_point_minimization = MagicMock(return_value=mock_out) - - # Patch MAGEMin_C and jlconvert at the module level - with patch.object(mo_module, 'MAGEMin_C', mock_magemin_c), \ - patch.object(mo_module, 'jlconvert', side_effect=lambda t, v: np.array(v, dtype=float)): - mo.get_phase_chemistry_at_index = MagicMock(return_value=np.array([45.0, 55.0])) - mo.get_volume_between_radii = MagicMock(return_value=1e12) - mo.pressure_to_radius = MagicMock(return_value=1400.0) - mo.radius_to_pressure = MagicMock(return_value=5.0) - - starting_melt = np.array([45.0, 55.0]) - mo.run_fractional_stages( - starting_melt=starting_melt, - p_start=5.0, - p_end=0.001, - vol_step=0.05, - starting_vol_frac=0.5, - n_stages=2, - ) - - np.testing.assert_array_equal( - mo.X, saved_X, - err_msg="self.X was not restored after run_fractional_stages" - ) - - -# =========================================================================== -# Fix H4: starting_melt length validation -# =========================================================================== -class TestStartingMeltValidation(unittest.TestCase): - """Fix H4: run_fractional_stages must reject mismatched starting_melt length.""" - - @patch('phasetools.models.magma_ocean.MAGEMinBase.__init__', return_value=None) - def test_starting_melt_length_mismatch(self, mock_base_init): - from phasetools.models.magma_ocean import MagmaOcean - - mo = MagmaOcean.__new__(MagmaOcean) - mo._Xoxides_py = ['SiO2', 'Al2O3', 'MgO'] - mo.sys_in = 'mol' - mo.data = MagicMock() - mo.X = np.array([33.0, 33.0, 34.0]) - mo.radius_body = 1737.1 - mo.radius_core = 330.0 - mo.g = 1.62 - mo.rho_avg = 3350.0 - - # starting_melt has 2 elements, but _Xoxides_py has 3 - with self.assertRaises(ValueError) as ctx: - mo.run_fractional_stages( - starting_melt=np.array([50.0, 50.0]), - p_start=5.0, - p_end=0.001, - ) - self.assertIn("does not match", str(ctx.exception)) - - if __name__ == '__main__': unittest.main()