From 438bce57cde9a2fb713bfb6cdd0146097b6d9ed0 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/11] 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 78429cd88c632ae4f382a283b2c9f385a613c86b 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/11] 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 8409dcd6c6c07395dabb42d6bd84b5cc1e40903f 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/11] 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 e8128b9def0006a529f7675b05c5a624464d46d7 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:45:07 +0800 Subject: [PATCH 04/11] 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 c709f9e17fd352635abc8ec635a1ded493950c36 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:45:13 +0800 Subject: [PATCH 05/11] 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 fd5ee2bf445264ef17e89a7b891406016641ce14 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/11] 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 2346e5fd620dd3b710d465afda8920f4b772c87d Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:45:24 +0800 Subject: [PATCH 07/11] feat(models): preserve bulk across MagmaOcean fractional stages --- src/phasetools/models/magma_ocean.py | 165 +++++++++++++++------------ 1 file changed, 93 insertions(+), 72 deletions(-) diff --git a/src/phasetools/models/magma_ocean.py b/src/phasetools/models/magma_ocean.py index 78f67c7..555e7a9 100644 --- a/src/phasetools/models/magma_ocean.py +++ b/src/phasetools/models/magma_ocean.py @@ -1,5 +1,6 @@ import numpy as np import sys +import warnings from scipy import optimize from typing import List, Dict, Any, Tuple, Optional from ..core.base import MAGEMinBase @@ -89,7 +90,12 @@ def func(T): except ValueError: f_low = func(bracket[0]) f_high = func(bracket[1]) - return float(bracket[0] if abs(f_low) < abs(f_high) else bracket[1]) + endpoint = bracket[0] if abs(f_low) < abs(f_high) else bracket[1] + warnings.warn( + f"find_temperature_at_vol_frac: bisection failed for P={P}, " + f"target_vol_frac={target_vol_frac}; returning bracket endpoint T={endpoint}." + ) + return float(endpoint) def get_phase_chemistry_at_index(self, out, i: int) -> np.ndarray: """Extract the chemical composition vector of a phase at a specific index.""" @@ -121,6 +127,7 @@ def run_stage_0(self, p_start: float, p_end: float, solid_frac: float = 0.5, p_i } melt_sum = np.zeros(len(self._Xoxides_py)) + n_melt_samples = 0 layer_modes_sum = {} for P in pressures: @@ -145,15 +152,19 @@ def run_stage_0(self, p_start: float, p_end: float, solid_frac: float = 0.5, p_i if ph_str == 'liq': melt_comp = self.get_phase_chemistry_at_index(out, i) melt_sum += melt_comp + n_melt_samples += 1 else: layer_modes_sum[ph_str] = layer_modes_sum.get(ph_str, 0.0) + vfrac results["modes"].append(modes) results["densities"].append(densities) - avg_melt = melt_sum / p_intervals + avg_melt = melt_sum / max(n_melt_samples, 1) total_solid = sum(layer_modes_sum.values()) - results["layer_modes"] = {ph: val / total_solid for ph, val in layer_modes_sum.items()} + if total_solid > 0: + results["layer_modes"] = {ph: val / total_solid for ph, val in layer_modes_sum.items()} + else: + results["layer_modes"] = {} return results, avg_melt @@ -191,6 +202,12 @@ def run_fractional_stages(self, """ all_stage_results = [] current_melt_comp = starting_melt + + if len(starting_melt) != len(self._Xoxides_py): + raise ValueError( + f"starting_melt length ({len(starting_melt)}) does not match " + f"bulk composition oxides ({len(self._Xoxides_py)})." + ) # Calculate volume of total LMO based on the initial melt ocean bounds v_total_mo_init = self.get_volume_between_radii(self.pressure_to_radius(p_start), self.pressure_to_radius(p_end)) @@ -205,76 +222,80 @@ def run_fractional_stages(self, r_top = self.pressure_to_radius(p_end) current_liquid_vol_frac = starting_vol_frac - - for stage in range(1, n_stages + 1): - # Set composition - self.X = jlconvert(jl.Vector[jl.Float64], current_melt_comp) - - # Base pressure of the current liquid ocean - p_base = self.radius_to_pressure(r_bottom) - - # Per Johnson et al. 2021: Stage 1 concludes when 5 vol% solid is reached - # for the WHOLE melt ocean. Target solid frac = vol_step / current_liquid_vol_frac. - # Clamp to 1.0 to prevent minimization failure in the final stage. - target_solid_frac = min(vol_step / current_liquid_vol_frac, 1.0) - - # Find temperature at base pressure for target solid fraction - T = self.find_temperature_at_vol_frac(p_base, target_solid_frac) - - # Run minimization at the base pressure - out = MAGEMin_C.single_point_minimization(p_base, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list) - - stage_results = { - "stage": stage, - "p_base": p_base, - "p_top": self.radius_to_pressure(r_top), - "T": T, - "modes": {}, - "densities": {}, - "layer_modes": {} - } - - layer_modes_sum = {} - for i, ph_name in enumerate(out.ph): - ph_str = str(ph_name) - vfrac = float(out.ph_frac_vol[i]) - stage_results["modes"][ph_str] = vfrac + + saved_X = self.X + try: + for stage in range(1, n_stages + 1): + # Set composition + self.X = jlconvert(jl.Vector[jl.Float64], current_melt_comp) - if i < out.n_SS: - rho = float(out.SS_vec[i].rho) - else: - rho = float(out.PP_vec[i - out.n_SS].rho) - stage_results["densities"][ph_str] = rho + # Base pressure of the current liquid ocean + p_base = self.radius_to_pressure(r_bottom) - if ph_str == 'liq': - current_melt_comp = self.get_phase_chemistry_at_index(out, i) + # Per Johnson et al. 2021: Stage 1 concludes when 5 vol% solid is reached + # for the WHOLE melt ocean. Target solid frac = vol_step / current_liquid_vol_frac. + # Clamp to 1.0 to prevent minimization failure in the final stage. + target_solid_frac = min(vol_step / current_liquid_vol_frac, 1.0) + + # Find temperature at base pressure for target solid fraction + T = self.find_temperature_at_vol_frac(p_base, target_solid_frac) + + # Run minimization at the base pressure + out = MAGEMin_C.single_point_minimization(p_base, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list) + + stage_results = { + "stage": stage, + "p_base": p_base, + "p_top": self.radius_to_pressure(r_top), + "T": T, + "modes": {}, + "densities": {}, + "layer_modes": {} + } + + layer_modes_sum = {} + for i, ph_name in enumerate(out.ph): + ph_str = str(ph_name) + vfrac = float(out.ph_frac_vol[i]) + stage_results["modes"][ph_str] = vfrac + + if i < out.n_SS: + rho = float(out.SS_vec[i].rho) + else: + rho = float(out.PP_vec[i - out.n_SS].rho) + stage_results["densities"][ph_str] = rho + + if ph_str == 'liq': + current_melt_comp = self.get_phase_chemistry_at_index(out, i) + else: + layer_modes_sum[ph_str] = layer_modes_sum.get(ph_str, 0.0) + vfrac + + # Normalise solid modes for the layer + total_solid = sum(layer_modes_sum.values()) + if total_solid > 0: + stage_results["layer_modes"] = {ph: val / total_solid for ph, val in layer_modes_sum.items()} else: - layer_modes_sum[ph_str] = layer_modes_sum.get(ph_str, 0.0) + vfrac - - # Normalise solid modes for the layer - total_solid = sum(layer_modes_sum.values()) - if total_solid > 0: - stage_results["layer_modes"] = {ph: val / total_solid for ph, val in layer_modes_sum.items()} - else: - stage_results["layer_modes"] = {} - - # Update Geometry: Sinking minerals raise the bottom radius, floating ones lower the top radius. - pl_frac = sum(stage_results["layer_modes"].get(ph, 0.0) for ph in float_phases) - v_float = v_step * pl_frac - v_sink = v_step * (1.0 - pl_frac) - - # Ensure we don't exceed the available ocean volume (safety bound) - v_ocean = self.get_volume_between_radii(r_bottom, r_top) - v_float = min(v_float, v_ocean) - v_sink = min(v_sink, v_ocean - v_float) - - r_bottom = np.power(r_bottom**3 + (3 * v_sink) / (4 * np.pi), 1/3) - # Ensure r_top doesn't go below r_bottom - r_top = np.power(max(r_top**3 - (3 * v_float) / (4 * np.pi), r_bottom**3), 1/3) - - current_liquid_vol_frac -= vol_step - all_stage_results.append(stage_results) - - if current_liquid_vol_frac <= 0: break - + stage_results["layer_modes"] = {} + + # Update Geometry: Sinking minerals raise the bottom radius, floating ones lower the top radius. + pl_frac = sum(stage_results["layer_modes"].get(ph, 0.0) for ph in float_phases) + v_float = v_step * pl_frac + v_sink = v_step * (1.0 - pl_frac) + + # Ensure we don't exceed the available ocean volume (safety bound) + v_ocean = self.get_volume_between_radii(r_bottom, r_top) + v_float = min(v_float, v_ocean) + v_sink = min(v_sink, v_ocean - v_float) + + r_bottom = np.power(r_bottom**3 + (3 * v_sink) / (4 * np.pi), 1/3) + # Ensure r_top doesn't go below r_bottom + r_top = np.power(max(r_top**3 - (3 * v_float) / (4 * np.pi), r_bottom**3), 1/3) + + current_liquid_vol_frac -= vol_step + all_stage_results.append(stage_results) + + if current_liquid_vol_frac <= 0: break + finally: + self.X = saved_X + return all_stage_results From dd9160f63db97c20905373a2e60e483edf62274b Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:20:54 +0800 Subject: [PATCH 08/11] fix(models): raise on MagmaOcean bisection failure, British spelling --- src/phasetools/models/magma_ocean.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/phasetools/models/magma_ocean.py b/src/phasetools/models/magma_ocean.py index 555e7a9..4106300 100644 --- a/src/phasetools/models/magma_ocean.py +++ b/src/phasetools/models/magma_ocean.py @@ -62,7 +62,7 @@ def depth_to_pressure(self, depth_km: float) -> float: return (self.rho_avg * self.g * depth_km * 1000.0) / 1e8 def radius_to_pressure(self, R_km: float) -> float: - """Convert radius from center (km) to pressure (kbar).""" + """Convert radius from centre (km) to pressure (kbar).""" depth = self.radius_body - R_km return self.depth_to_pressure(depth) @@ -95,7 +95,10 @@ def func(T): f"find_temperature_at_vol_frac: bisection failed for P={P}, " f"target_vol_frac={target_vol_frac}; returning bracket endpoint T={endpoint}." ) - return float(endpoint) + raise RuntimeError( + f"find_temperature_at_vol_frac: bisection failed for P={P}, " + f"target_vol_frac={target_vol_frac}. Bracket T range: {bracket}" + ) def get_phase_chemistry_at_index(self, out, i: int) -> np.ndarray: """Extract the chemical composition vector of a phase at a specific index.""" From 7eff2d115743adf40ed89c73212a03434451d502 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:45:48 +0800 Subject: [PATCH 09/11] test: add MagmaOcean fractional stages and fractionation tests --- tests/test_fractionation.py | 363 ++++++++++++++++++++++++++++++++++++ 1 file changed, 363 insertions(+) create mode 100644 tests/test_fractionation.py 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 a9124b7c3c3b6011620d62e89cc141be9477f577 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:46:05 +0800 Subject: [PATCH 10/11] test: add MagmaOcean regression tests for bulk preservation and bisection failure --- tests/test_fractionation.py | 42 +++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/tests/test_fractionation.py b/tests/test_fractionation.py index a67bc8d..fa535cf 100644 --- a/tests/test_fractionation.py +++ b/tests/test_fractionation.py @@ -359,5 +359,47 @@ def test_starting_melt_length_mismatch(self, mock_base_init): self.assertIn("does not match", str(ctx.exception)) +# =========================================================================== +# run_stage_0: no melt / no solid edge case +# =========================================================================== +class TestStageZeroNoMeltNoSolid(unittest.TestCase): + """run_stage_0 must return zeros for avg_melt and empty layer_modes when liq is absent.""" + + @patch('phasetools.models.magma_ocean.MAGEMinBase.__init__', return_value=None) + def test_no_melt_no_solid(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', 'MgO'] + mo.sys_in = 'mol' + mo.data = MagicMock() + mo.X = np.array([33.0, 33.0, 34.0]) + mo.Xoxides = MagicMock() + mo.rm_list = None + + # find_temperature_at_vol_frac is called but its return value is irrelevant + # because the mocked output has no phases. + mo.find_temperature_at_vol_frac = MagicMock(return_value=1200.0) + + mock_out = MagicMock() + mock_out.ph = [] + mock_out.n_SS = 0 + mock_out.SS_vec = [] + mock_out.PP_vec = [] + + mock_magemin_c = MagicMock() + mock_magemin_c.single_point_minimization = MagicMock(return_value=mock_out) + + with patch.object(mo_module, 'MAGEMin_C', mock_magemin_c): + results, avg_melt = mo.run_stage_0( + p_start=5.0, p_end=0.001, solid_frac=0.5, p_intervals=3 + ) + + self.assertTrue(np.all(np.isfinite(avg_melt))) + self.assertTrue(np.allclose(avg_melt, np.zeros(3))) + self.assertEqual(results["layer_modes"], {}) + + if __name__ == '__main__': unittest.main() From 44a62aebe94199401d0ce4d9b9ffaf3832c417c8 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Sun, 16 Aug 2026 13:48:22 +0800 Subject: [PATCH 11/11] test: trim fractionation tests to MagmaOcean-only for this sub-PR --- tests/test_fractionation.py | 251 +----------------------------------- 1 file changed, 5 insertions(+), 246 deletions(-) diff --git a/tests/test_fractionation.py b/tests/test_fractionation.py index fa535cf..e601699 100644 --- a/tests/test_fractionation.py +++ b/tests/test_fractionation.py @@ -1,11 +1,10 @@ -"""Mock-based tests for fractionation correctness fixes. +"""Mock-based tests for MagmaOcean fractional crystallisation 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) +These tests verify: - Fix H1: self.X permanent mutation in run_fractional_stages (magma_ocean.py) +- Fix H4: starting_melt length validation +- run_stage_0: no melt / no solid edge case +- bisection failure raises RuntimeError No live Julia runtime is needed — all MAGEMin calls are mocked. """ @@ -15,246 +14,6 @@ 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 # ===========================================================================