API Reference

Model Classes

globalsplinefit.GSFEnergy

Bases: GSFBase

GSF model evaluated at total energy per nucleus in GeV.

Examples:

>>> from globalsplinefit import GSFEnergy
>>> import numpy as np
>>> model = GSFEnergy()
>>> energy = np.logspace(0, 3, 100)  # 1 GeV to 1 TeV
>>> proton_flux = model.flux(energy, "p")
>>> he_flux = model.flux(energy, "He")
>>> total_flux = model.total_flux(energy)

Source code in src/globalsplinefit/model.py

class GSFEnergy(GSFBase):
    """GSF model evaluated at total energy per nucleus in GeV.

    Examples
    --------
        >>> from globalsplinefit import GSFEnergy
        >>> import numpy as np
        >>> model = GSFEnergy()
        >>> energy = np.logspace(0, 3, 100)  # 1 GeV to 1 TeV
        >>> proton_flux = model.flux(energy, "p")
        >>> he_flux = model.flux(energy, "He")
        >>> total_flux = model.total_flux(energy)
    """

    def _transform_energy(
        self,
        energy: ArrayLike,
        target: Target,  # noqa: ARG002
    ) -> np.ndarray:
        """Transform input energy to total energy per nucleus.

        Subclasses override this to convert from kinetic energy, etc.
        """
        return self._scale_energy(self._as_1d_values(energy, "energy"))

    def _total_energy_for_sid(self, energy: np.ndarray, sid) -> np.ndarray:  # noqa: ARG002
        """Convert prepared input values to total energy for one species."""
        return energy

    def flux(
        self,
        energy: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Calculate flux for target (group or elements).

        Parameters
        ----------
        energy
            Total energy per nucleus in GeV. Can be scalar, list, or numpy array.
        target
            Target specification:
            - String: Group name ("p", "proton", "H", "He", "O*", "Fe*", etc.)
              or species name ("D" for deuterium, in model versions that
              carry it as its own species)
            - Integer: Single element atomic number (e.g., 1 for H, 2 for He)
            - Species id ``(Z, A)``: a single isotope, e.g. ``(1, 2.014)``
            - List: Multiple elements from same group (e.g., [6,7,8] for CNO)
        time_interval
            Time period specification:
            - None: Uses the default_time_interval set during initialization
            - "LIS": Local Interstellar Spectrum (no modulation)
            - tuple: (start, end) in YYYYMM format, end month EXCLUSIVE, e.g. (200901, 201001) for calendar year 2009
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. Nuclei with rigidity
            below this value are excluded. None uses the default.

        Returns
        -------
            Array of differential flux values in units of particles/(m²·s·sr·GeV).
            Shape matches the input energy array.

        Examples
        --------
            >>> flux_p = model.flux([1, 10, 100], "p")  # proton flux at 1, 10, 100 GeV
            >>> flux_he = model.flux(energy_array, "He")  # helium flux
            >>> flux_cno = model.flux(energy_array, [6, 7, 8])  # combined CNO flux
            >>> flux_lis = model.flux(energy_array, "p", time_interval="LIS")  # LIS flux
        """
        rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
        zlist, _group_leader = self._resolve_z(target)
        input_energy = self._transform_energy(energy, target)

        flux = np.zeros_like(input_energy, dtype=float)
        for sid in self._target_sids(zlist):  # charges expand to species (p, D, …)
            total_energy = self._total_energy_for_sid(input_energy, sid)
            mask = self._rigidity_cutoff_mask(sid, total_energy, rigidity_cutoff)
            flux += self._element_flux(sid, total_energy, time_interval) * mask
        return flux

    def jacobian(
        self,
        energy: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Calculate Jacobian of flux for uncertainty propagation."""
        rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
        zlist, _group_leader = self._resolve_z(target)
        input_energy = self._transform_energy(energy, target)

        jac = 0.0
        for sid in self._target_sids(zlist):
            total_energy = self._total_energy_for_sid(input_energy, sid)
            mask = self._rigidity_cutoff_mask(sid, total_energy, rigidity_cutoff)
            jac += (
                self._element_flux_jacobian(sid, total_energy, time_interval)
                * mask[:, np.newaxis]
            )
        return np.asarray(jac)

flux(energy: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Calculate flux for target (group or elements).

Parameters:

Name	Type	Description	Default
energy	ArrayLike	Total energy per nucleus in GeV. Can be scalar, list, or numpy array.	required
target	Target	Target specification: - String: Group name ("p", "proton", "H", "He", "O", "Fe", etc.) or species name ("D" for deuterium, in model versions that carry it as its own species) - Integer: Single element atomic number (e.g., 1 for H, 2 for He) - Species id (Z, A): a single isotope, e.g. (1, 2.014) - List: Multiple elements from same group (e.g., [6,7,8] for CNO)	required
time_interval	tuple[int, int] | str | None	Time period specification: - None: Uses the default_time_interval set during initialization - "LIS": Local Interstellar Spectrum (no modulation) - tuple: (start, end) in YYYYMM format, end month EXCLUSIVE, e.g. (200901, 201001) for calendar year 2009	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. Nuclei with rigidity below this value are excluded. None uses the default.	None

Returns:

Type	Description
Array of differential flux values in units of particles/(m²·s·sr·GeV).	Shape matches the input energy array.

Examples:

>>> flux_p = model.flux([1, 10, 100], "p")  # proton flux at 1, 10, 100 GeV
>>> flux_he = model.flux(energy_array, "He")  # helium flux
>>> flux_cno = model.flux(energy_array, [6, 7, 8])  # combined CNO flux
>>> flux_lis = model.flux(energy_array, "p", time_interval="LIS")  # LIS flux

Source code in src/globalsplinefit/model.py

def flux(
    self,
    energy: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Calculate flux for target (group or elements).

    Parameters
    ----------
    energy
        Total energy per nucleus in GeV. Can be scalar, list, or numpy array.
    target
        Target specification:
        - String: Group name ("p", "proton", "H", "He", "O*", "Fe*", etc.)
          or species name ("D" for deuterium, in model versions that
          carry it as its own species)
        - Integer: Single element atomic number (e.g., 1 for H, 2 for He)
        - Species id ``(Z, A)``: a single isotope, e.g. ``(1, 2.014)``
        - List: Multiple elements from same group (e.g., [6,7,8] for CNO)
    time_interval
        Time period specification:
        - None: Uses the default_time_interval set during initialization
        - "LIS": Local Interstellar Spectrum (no modulation)
        - tuple: (start, end) in YYYYMM format, end month EXCLUSIVE, e.g. (200901, 201001) for calendar year 2009
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. Nuclei with rigidity
        below this value are excluded. None uses the default.

    Returns
    -------
        Array of differential flux values in units of particles/(m²·s·sr·GeV).
        Shape matches the input energy array.

    Examples
    --------
        >>> flux_p = model.flux([1, 10, 100], "p")  # proton flux at 1, 10, 100 GeV
        >>> flux_he = model.flux(energy_array, "He")  # helium flux
        >>> flux_cno = model.flux(energy_array, [6, 7, 8])  # combined CNO flux
        >>> flux_lis = model.flux(energy_array, "p", time_interval="LIS")  # LIS flux
    """
    rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
    zlist, _group_leader = self._resolve_z(target)
    input_energy = self._transform_energy(energy, target)

    flux = np.zeros_like(input_energy, dtype=float)
    for sid in self._target_sids(zlist):  # charges expand to species (p, D, …)
        total_energy = self._total_energy_for_sid(input_energy, sid)
        mask = self._rigidity_cutoff_mask(sid, total_energy, rigidity_cutoff)
        flux += self._element_flux(sid, total_energy, time_interval) * mask
    return flux

jacobian(energy: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Calculate Jacobian of flux for uncertainty propagation.

Source code in src/globalsplinefit/model.py

def jacobian(
    self,
    energy: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Calculate Jacobian of flux for uncertainty propagation."""
    rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
    zlist, _group_leader = self._resolve_z(target)
    input_energy = self._transform_energy(energy, target)

    jac = 0.0
    for sid in self._target_sids(zlist):
        total_energy = self._total_energy_for_sid(input_energy, sid)
        mask = self._rigidity_cutoff_mask(sid, total_energy, rigidity_cutoff)
        jac += (
            self._element_flux_jacobian(sid, total_energy, time_interval)
            * mask[:, np.newaxis]
        )
    return np.asarray(jac)

globalsplinefit.GSFKineticEnergy

Bases: GSFEnergy

GSF model expecting kinetic energy per nucleus in GeV.

Converts kinetic energy to total energy (kinetic + rest mass) internally before using the GSFEnergy calculation methods.

Examples:

>>> from globalsplinefit import GSFKineticEnergy
>>> import numpy as np
>>> model = GSFKineticEnergy()
>>> kinetic_energy = np.logspace(0, 3, 100)  # 1 GeV to 1 TeV kinetic energy
>>> proton_flux = model.flux(kinetic_energy, "p")
>>> he_flux = model.flux(kinetic_energy, "He")
>>> total_flux = model.total_flux(kinetic_energy)

Source code in src/globalsplinefit/model.py

class GSFKineticEnergy(GSFEnergy):
    """GSF model expecting kinetic energy per nucleus in GeV.

    Converts kinetic energy to total energy (kinetic + rest mass) internally
    before using the GSFEnergy calculation methods.

    Examples
    --------
        >>> from globalsplinefit import GSFKineticEnergy
        >>> import numpy as np
        >>> model = GSFKineticEnergy()
        >>> kinetic_energy = np.logspace(0, 3, 100)  # 1 GeV to 1 TeV kinetic energy
        >>> proton_flux = model.flux(kinetic_energy, "p")
        >>> he_flux = model.flux(kinetic_energy, "He")
        >>> total_flux = model.total_flux(kinetic_energy)
    """

    def _transform_energy(
        self, kinetic_energy: ArrayLike, _target: Target
    ) -> np.ndarray:
        """Prepare kinetic energy; species rest masses are added during summation."""
        return self._scale_energy(self._as_1d_values(kinetic_energy, "kinetic energy"))

    def _total_energy_for_sid(self, energy: np.ndarray, sid) -> np.ndarray:
        """Convert kinetic to total energy using the evaluated species' mass."""
        return energy + self.z_to_a[sid] * NUCLEON_MASS_GEV

globalsplinefit.GSFRigidity

Bases: GSFBase

GSF model evaluated at magnetic rigidity in GV.

Solar modulation uses the force-field approximation. energy_scale applies to the energy-based models only.

Examples:

>>> from globalsplinefit import GSFRigidity
>>> import numpy as np
>>> model = GSFRigidity()
>>> rigidity = np.logspace(0, 3, 100)  # 1 GV to 1 TV
>>> proton_flux = model.flux(rigidity, "p")  # Solar Cycle 24 average
>>> proton_flux_lis = model.flux(rigidity, "p", time_interval="LIS")
>>> proton_flux_2009 = model.flux(rigidity, "p", time_interval=(200901, 201001))
>>> total_flux = model.total_flux(rigidity)

Source code in src/globalsplinefit/model.py

class GSFRigidity(GSFBase):
    """GSF model evaluated at magnetic rigidity in GV.

    Solar modulation uses the force-field approximation. ``energy_scale``
    applies to the energy-based models only.

    Examples
    --------
        >>> from globalsplinefit import GSFRigidity
        >>> import numpy as np
        >>> model = GSFRigidity()
        >>> rigidity = np.logspace(0, 3, 100)  # 1 GV to 1 TV
        >>> proton_flux = model.flux(rigidity, "p")  # Solar Cycle 24 average
        >>> proton_flux_lis = model.flux(rigidity, "p", time_interval="LIS")
        >>> proton_flux_2009 = model.flux(rigidity, "p", time_interval=(200901, 201001))
        >>> total_flux = model.total_flux(rigidity)
    """

    def _rigidity_phi_transform(
        self, sid, rigidity: np.ndarray, phi: float
    ) -> tuple[np.ndarray, np.ndarray]:
        """Convert Earth rigidity to interstellar rigidity under force-field phi.

        Force field: total-energy loss Z*phi (Gleeson-Axford). Returns R_IS and the
        dN/dR prefactor Lambda = (E_IS/E) * (R/R_IS)**3.
        """
        sid = self._as_sid(sid)
        nucleon_mass = NUCLEON_MASS_GEV
        z = sid[0]
        a = self.z_to_a[sid]
        m = a * nucleon_mass

        # Earth energy from input R; interstellar energy after Z*phi loss
        E = np.sqrt((z * rigidity) ** 2 + m**2)
        E_is = E + z * phi

        # Interstellar rigidity
        with np.errstate(divide="ignore", invalid="ignore"):
            p2_is = E_is**2 - m**2
            p2_is[p2_is < 0] = 0.0
            R_is = np.sqrt(p2_is) / z

            # dN/dR force-field prefactor. Public inputs are non-negative.
            Lambda = (E_is / E) * (rigidity / (R_is + _ZERO_GUARD)) ** 3
        return R_is, Lambda

    def flux(
        self,
        rigidity: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Calculate flux for target (group or elements).

        Parameters
        ----------
        rigidity
            Magnetic rigidity in GV. Can be scalar, list, or numpy array.
        target
            Target specification (same as GSFEnergy.flux).
        time_interval
            Time period specification:
            - None: Uses the default_time_interval set during initialization
            - "LIS": Local Interstellar Spectrum (no modulation)
            - tuple: (start, end) in YYYYMM format, end month EXCLUSIVE, e.g. (200901, 201001) for calendar year 2009
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. Flux at rigidities
            below this value is set to zero. None uses the default.

        Returns
        -------
            Array of differential flux values in units of particles/(m²·s·sr·GV).
            Shape matches the input rigidity array.
        """
        time_interval = self._resolve_time_interval(time_interval)
        rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
        zlist, _ = self._resolve_z(target)
        rigidity = self._as_1d_values(rigidity, "rigidity")

        # φ list and period-average weights (length 1 with 0.0 for LIS)
        phis, phi_weights = self._phi_list(time_interval)

        flux = np.zeros_like(rigidity, dtype=float)
        for phi, w in zip(phis, phi_weights, strict=True):
            for sid in self._target_sids(zlist):  # charges expand to species (p, D, …)
                if phi == 0.0:
                    flux += w * self._rigidity_flux_lis(sid, rigidity)
                else:
                    R_is, fac = self._rigidity_phi_transform(sid, rigidity, phi)
                    flux += w * self._rigidity_flux_lis(sid, R_is) * fac

        # Apply rigidity cutoff (in rigidity space, cutoff is Z-independent)
        if rigidity_cutoff is not None:
            if self.cutoff_width > 0:
                mask = _sigmoid((rigidity - rigidity_cutoff) / self.cutoff_width)
                flux *= mask
            else:
                flux[rigidity < rigidity_cutoff] = 0.0

        return flux

    def jacobian(
        self,
        rigidity: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Jacobian including solar modulation."""
        time_interval = self._resolve_time_interval(time_interval)
        rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
        zlist, _ = self._resolve_z(target)
        rigidity = self._as_1d_values(rigidity, "rigidity")

        phis, phi_weights = self._phi_list(time_interval)
        jac = 0.0
        for phi, w in zip(phis, phi_weights, strict=True):
            for sid in self._target_sids(zlist):
                leading, _ratio = self.flux_ratio[sid]

                if phi == 0.0:
                    scale = self._subleading_scale(sid, leading, rigidity)
                    contrib = self._rigidity_flux_jacobian(leading, rigidity)
                    if sid != leading:
                        contrib = contrib * np.asarray(scale)[:, None]
                else:
                    R_is, fac = self._rigidity_phi_transform(sid, rigidity, phi)
                    scale = self._subleading_scale(sid, leading, R_is)
                    contrib = (
                        self._rigidity_flux_jacobian(leading, R_is)
                        * (fac * scale)[:, None]
                    )
                jac += w * contrib
        jac = np.asarray(jac)

        # Apply rigidity cutoff (in rigidity space, cutoff is Z-independent)
        if rigidity_cutoff is not None:
            if self.cutoff_width > 0:
                mask = _sigmoid((rigidity - rigidity_cutoff) / self.cutoff_width)
                jac *= mask[:, None]
            else:
                jac[rigidity < rigidity_cutoff, :] = 0.0

        return jac

flux(rigidity: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Calculate flux for target (group or elements).

Parameters:

Name	Type	Description	Default
rigidity	ArrayLike	Magnetic rigidity in GV. Can be scalar, list, or numpy array.	required
target	Target	Target specification (same as GSFEnergy.flux).	required
time_interval	tuple[int, int] | str | None	Time period specification: - None: Uses the default_time_interval set during initialization - "LIS": Local Interstellar Spectrum (no modulation) - tuple: (start, end) in YYYYMM format, end month EXCLUSIVE, e.g. (200901, 201001) for calendar year 2009	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. Flux at rigidities below this value is set to zero. None uses the default.	None

Returns:

Type	Description
Array of differential flux values in units of particles/(m²·s·sr·GV).	Shape matches the input rigidity array.

Source code in src/globalsplinefit/model.py

def flux(
    self,
    rigidity: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Calculate flux for target (group or elements).

    Parameters
    ----------
    rigidity
        Magnetic rigidity in GV. Can be scalar, list, or numpy array.
    target
        Target specification (same as GSFEnergy.flux).
    time_interval
        Time period specification:
        - None: Uses the default_time_interval set during initialization
        - "LIS": Local Interstellar Spectrum (no modulation)
        - tuple: (start, end) in YYYYMM format, end month EXCLUSIVE, e.g. (200901, 201001) for calendar year 2009
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. Flux at rigidities
        below this value is set to zero. None uses the default.

    Returns
    -------
        Array of differential flux values in units of particles/(m²·s·sr·GV).
        Shape matches the input rigidity array.
    """
    time_interval = self._resolve_time_interval(time_interval)
    rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
    zlist, _ = self._resolve_z(target)
    rigidity = self._as_1d_values(rigidity, "rigidity")

    # φ list and period-average weights (length 1 with 0.0 for LIS)
    phis, phi_weights = self._phi_list(time_interval)

    flux = np.zeros_like(rigidity, dtype=float)
    for phi, w in zip(phis, phi_weights, strict=True):
        for sid in self._target_sids(zlist):  # charges expand to species (p, D, …)
            if phi == 0.0:
                flux += w * self._rigidity_flux_lis(sid, rigidity)
            else:
                R_is, fac = self._rigidity_phi_transform(sid, rigidity, phi)
                flux += w * self._rigidity_flux_lis(sid, R_is) * fac

    # Apply rigidity cutoff (in rigidity space, cutoff is Z-independent)
    if rigidity_cutoff is not None:
        if self.cutoff_width > 0:
            mask = _sigmoid((rigidity - rigidity_cutoff) / self.cutoff_width)
            flux *= mask
        else:
            flux[rigidity < rigidity_cutoff] = 0.0

    return flux

jacobian(rigidity: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Jacobian including solar modulation.

Source code in src/globalsplinefit/model.py

def jacobian(
    self,
    rigidity: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Jacobian including solar modulation."""
    time_interval = self._resolve_time_interval(time_interval)
    rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
    zlist, _ = self._resolve_z(target)
    rigidity = self._as_1d_values(rigidity, "rigidity")

    phis, phi_weights = self._phi_list(time_interval)
    jac = 0.0
    for phi, w in zip(phis, phi_weights, strict=True):
        for sid in self._target_sids(zlist):
            leading, _ratio = self.flux_ratio[sid]

            if phi == 0.0:
                scale = self._subleading_scale(sid, leading, rigidity)
                contrib = self._rigidity_flux_jacobian(leading, rigidity)
                if sid != leading:
                    contrib = contrib * np.asarray(scale)[:, None]
            else:
                R_is, fac = self._rigidity_phi_transform(sid, rigidity, phi)
                scale = self._subleading_scale(sid, leading, R_is)
                contrib = (
                    self._rigidity_flux_jacobian(leading, R_is)
                    * (fac * scale)[:, None]
                )
            jac += w * contrib
    jac = np.asarray(jac)

    # Apply rigidity cutoff (in rigidity space, cutoff is Z-independent)
    if rigidity_cutoff is not None:
        if self.cutoff_width > 0:
            mask = _sigmoid((rigidity - rigidity_cutoff) / self.cutoff_width)
            jac *= mask[:, None]
        else:
            jac[rigidity < rigidity_cutoff, :] = 0.0

    return jac

globalsplinefit.GSFEnergyPerNucleon

Bases: GSFBase

GSF model for nucleon flux calculations (energy per nucleon).

The nucleon flux is calculated by: - Proton flux = sum over nuclei: flux(nucleus) × A × Z - Neutron flux = sum over nuclei: flux(nucleus) × A × (A-Z)

Where A is atomic mass number and Z is atomic number.

Examples:

>>> from globalsplinefit import GSFEnergyPerNucleon
>>> import numpy as np
>>> model = GSFEnergyPerNucleon()
>>> energy_per_nucleon = np.logspace(0, 2, 50)  # GeV/nucleon
>>> total_nucleons = model.flux(energy_per_nucleon, "He")  # shape (N,)
>>> p_and_n = model.p_and_n_flux(energy_per_nucleon, "He")  # shape (2, N)
>>> proton_nucleons = p_and_n[0]
>>> neutron_nucleons = p_and_n[1]

Source code in src/globalsplinefit/model.py

class GSFEnergyPerNucleon(GSFBase):
    """GSF model for nucleon flux calculations (energy per nucleon).

    The nucleon flux is calculated by:
    - Proton flux = sum over nuclei: flux(nucleus) × A × Z
    - Neutron flux = sum over nuclei: flux(nucleus) × A × (A-Z)

    Where A is atomic mass number and Z is atomic number.

    Examples
    --------
        >>> from globalsplinefit import GSFEnergyPerNucleon
        >>> import numpy as np
        >>> model = GSFEnergyPerNucleon()
        >>> energy_per_nucleon = np.logspace(0, 2, 50)  # GeV/nucleon
        >>> total_nucleons = model.flux(energy_per_nucleon, "He")  # shape (N,)
        >>> p_and_n = model.p_and_n_flux(energy_per_nucleon, "He")  # shape (2, N)
        >>> proton_nucleons = p_and_n[0]
        >>> neutron_nucleons = p_and_n[1]
    """

    def _transform_energy_per_nucleon(
        self,
        energy_per_nucleon: ArrayLike,
        target: Target,  # noqa: ARG002
    ) -> np.ndarray:
        """Transform input energy per nucleon. Subclasses override for kinetic energy."""
        return self._scale_energy(
            self._as_1d_values(energy_per_nucleon, "energy per nucleon")
        )

    def p_and_n_flux(
        self,
        energy_per_nucleon: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Calculate separate proton and neutron flux from target (group or elements).

        Parameters
        ----------
        energy_per_nucleon
            Energy per nucleon in GeV. Can be scalar, list, or array.
        target
            Target specification (same as GSFEnergy.flux).
        time_interval
            Time period specification:
            - None: Uses the default_time_interval set during initialization
            - "LIS": Local Interstellar Spectrum (no modulation)
            - tuple: (start, end) in YYYYMM format
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. Nuclei with rigidity
            below this value are excluded. None uses the default.

        Returns
        -------
            Array of shape (2, N) where N is the number of energy points:
            - [0, :]: Proton nucleon flux in particles/(m²·s·sr·GeV)
            - [1, :]: Neutron nucleon flux in particles/(m²·s·sr·GeV)

        Examples
        --------
            >>> nucleon_flux = model.p_and_n_flux([1, 10, 100], "He")
            >>> proton_flux = nucleon_flux[0]  # 2 protons per He nucleus
            >>> neutron_flux = nucleon_flux[1]  # 2 neutrons per He nucleus
        """
        time_interval = self._resolve_time_interval(time_interval)
        rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
        zlist, _group_leader = self._resolve_z(target)
        energy_per_nucleon = self._transform_energy_per_nucleon(
            energy_per_nucleon, target
        )

        flux = np.zeros((2, len(energy_per_nucleon)))
        for sid in self._target_sids(zlist):
            charge = sid[0]
            mass_scale = self.z_to_a[sid]
            nucleons = self.mass_number[sid]
            energy = energy_per_nucleon * mass_scale
            mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
            fl = self._element_flux(sid, energy, time_interval)
            flux[0] += fl * mass_scale * charge * mask
            flux[1] += fl * mass_scale * (nucleons - charge) * mask
        return flux

    def flux(
        self,
        energy_per_nucleon: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Calculate total nucleon flux (protons + neutrons).

        Equivalent to the sum of the two components returned by `p_and_n_flux`.

        Parameters
        ----------
        energy_per_nucleon
            Energy per nucleon in GeV.
        target
            Target specification (group name or list of Z values).
        time_interval
            Time period for solar modulation, or "LIS" for unmodulated.
            ``None`` uses the default set during initialization.
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

        Returns
        -------
            Array of total nucleon flux in particles/(m²·s·sr·GeV).
        """
        p_and_n = self.p_and_n_flux(
            energy_per_nucleon,
            target,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        return p_and_n[0] + p_and_n[1]

    def p_and_n_jacobian(
        self,
        energy_per_nucleon: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> tuple[np.ndarray, np.ndarray]:
        """Jacobian of separate proton and neutron flux.

        Parameters
        ----------
        energy_per_nucleon
            Energy per nucleon in GeV.
        target
            Target specification (group name or list of Z values).
        time_interval
            Time period for solar modulation, or "LIS" for unmodulated.
            ``None`` uses the default set during initialization.
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

        Returns
        -------
            Tuple ``(jac_p, jac_n)`` of proton and neutron flux Jacobians,
            each of shape ``(N, P)`` where ``P`` is the number of parameters.
        """
        time_interval = self._resolve_time_interval(time_interval)
        rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
        zlist, _group_leader = self._resolve_z(target)
        energy_per_nucleon = self._transform_energy_per_nucleon(
            energy_per_nucleon, target
        )

        jac_p = 0.0
        jac_n = 0.0
        for sid in self._target_sids(zlist):
            charge = sid[0]
            mass_scale = self.z_to_a[sid]
            nucleons = self.mass_number[sid]
            energy = energy_per_nucleon * mass_scale
            mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
            j = self._element_flux_jacobian(sid, energy, time_interval)
            jac_p += j * (mass_scale * charge * mask)[:, np.newaxis]
            jac_n += j * (mass_scale * (nucleons - charge) * mask)[:, np.newaxis]
        return np.asarray(jac_p), np.asarray(jac_n)

    def jacobian(
        self,
        energy_per_nucleon: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Jacobian of total nucleon flux (protons + neutrons).

        Parameters
        ----------
        energy_per_nucleon
            Energy per nucleon in GeV.
        target
            Target specification (group name or list of Z values).
        time_interval
            Time period for solar modulation, or "LIS" for unmodulated.
            ``None`` uses the default set during initialization.
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

        Returns
        -------
            Jacobian of total nucleon flux, shape ``(N, P)``.
        """
        jac_p, jac_n = self.p_and_n_jacobian(
            energy_per_nucleon,
            target,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        return jac_p + jac_n

    def p_and_n_covariance(
        self,
        target1: Target,
        target2: Target,
        energy_per_nucleon: ArrayLike,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> tuple[np.ndarray, np.ndarray]:
        """Separate covariance matrices for proton and neutron flux.

        Parameters
        ----------
        target1, target2
            Target specifications for the two targets being correlated.
        energy_per_nucleon
            Energy per nucleon in GeV.
        time_interval
            Time period for solar modulation. ``None`` uses the default.
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

        Returns
        -------
            Tuple ``(cov_pp, cov_nn)`` of the proton-proton and
            neutron-neutron covariance matrices, each shape ``(N, N)``.
        """
        zlist1, leader1 = self._resolve_z(target1)
        zlist2, leader2 = self._resolve_z(target2)

        # Calculate jacobians using the helper method
        jac1_p, jac1_n = self.p_and_n_jacobian(
            energy_per_nucleon,
            target1,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        jac2_p, jac2_n = (
            (jac1_p, jac1_n)
            if zlist2 == zlist1
            else self.p_and_n_jacobian(
                energy_per_nucleon,
                target2,
                time_interval=time_interval,
                rigidity_cutoff=rigidity_cutoff,
            )
        )

        # Key the covariance by the leader SPECIES id (Z, A). _resolve_z returns a
        # charge; a charge with a single species has a cov alias under the bare int,
        # but a charge carrying >1 species (p + D at Z=1) does not — so resolve to
        # the leader sid, which is always a real cov key.
        cov_key = (
            self._leader_by_charge.get(leader1, leader1),
            self._leader_by_charge.get(leader2, leader2),
        )

        if cov_key in self.cov:
            return (
                self._propagate_cov(jac1_p, jac2_p, self.cov[cov_key]),
                self._propagate_cov(jac1_n, jac2_n, self.cov[cov_key]),
            )
        else:
            # Elements not in the main groups (H, He, O, Fe) carry no covariance
            energy_per_nucleon = np.atleast_1d(energy_per_nucleon)
            n_energies = len(energy_per_nucleon)
            zero_cov = np.zeros((n_energies, n_energies))
            return zero_cov, zero_cov

    def p_and_n_total_flux(
        self,
        energy_or_rigidity: ArrayLike,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Separate proton and neutron flux summed over all element groups.

        Parameters
        ----------
        energy_or_rigidity
            Energy per nucleon in GeV.
        time_interval
            Time period for solar modulation. ``None`` uses the default.
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

        Returns
        -------
            Array of shape ``(2, N)``: row 0 is total proton flux,
            row 1 is total neutron flux, both in particles/(m²·s·sr·GeV).
        """
        energy_or_rigidity = np.atleast_1d(energy_or_rigidity)
        # Initialize with proper shape for nucleon flux (2, N)
        total_flux = np.zeros((2, len(energy_or_rigidity)), dtype=float)

        for group in self.active_groups:
            group_flux = self.p_and_n_flux(
                energy_or_rigidity,
                group,
                time_interval=time_interval,
                rigidity_cutoff=rigidity_cutoff,
            )
            total_flux += group_flux

        return total_flux

    def p_and_n_error(
        self,
        energy_per_nucleon: ArrayLike,
        target: Target,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Separate proton and neutron flux uncertainties (1-sigma).

        Parameters
        ----------
        energy_per_nucleon
            Energy per nucleon in GeV.
        target
            Target specification (group name or list of Z values).
        time_interval
            Time period for solar modulation. ``None`` uses the default.
        rigidity_cutoff
            Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

        Returns
        -------
            Array of shape ``(2, N)``: row 0 is proton uncertainties,
            row 1 is neutron uncertainties.
        """
        jac_p, jac_n = self.p_and_n_jacobian(
            energy_per_nucleon,
            target,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        block = self._covariance_block(target, target)
        if block is None:
            return np.zeros((2, jac_p.shape[0]))
        var_p = np.einsum("ni,ij,nj->n", jac_p, block, jac_p)
        var_n = np.einsum("ni,ij,nj->n", jac_n, block, jac_n)
        return np.sqrt(np.clip(np.vstack([var_p, var_n]), 0.0, None))

flux(energy_per_nucleon: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Calculate total nucleon flux (protons + neutrons).

Equivalent to the sum of the two components returned by p_and_n_flux.

Parameters:

Name	Type	Description	Default
energy_per_nucleon	ArrayLike	Energy per nucleon in GeV.	required
target	Target	Target specification (group name or list of Z values).	required
time_interval	tuple[int, int] | str | None	Time period for solar modulation, or "LIS" for unmodulated. None uses the default set during initialization.	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. None uses the default.	None

Returns:

Type	Description
Array of total nucleon flux in particles/(m²·s·sr·GeV).

Source code in src/globalsplinefit/model.py

def flux(
    self,
    energy_per_nucleon: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Calculate total nucleon flux (protons + neutrons).

    Equivalent to the sum of the two components returned by `p_and_n_flux`.

    Parameters
    ----------
    energy_per_nucleon
        Energy per nucleon in GeV.
    target
        Target specification (group name or list of Z values).
    time_interval
        Time period for solar modulation, or "LIS" for unmodulated.
        ``None`` uses the default set during initialization.
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

    Returns
    -------
        Array of total nucleon flux in particles/(m²·s·sr·GeV).
    """
    p_and_n = self.p_and_n_flux(
        energy_per_nucleon,
        target,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    return p_and_n[0] + p_and_n[1]

jacobian(energy_per_nucleon: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Jacobian of total nucleon flux (protons + neutrons).

Parameters:

Name	Type	Description	Default
energy_per_nucleon	ArrayLike	Energy per nucleon in GeV.	required
target	Target	Target specification (group name or list of Z values).	required
time_interval	tuple[int, int] | str | None	Time period for solar modulation, or "LIS" for unmodulated. None uses the default set during initialization.	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. None uses the default.	None

Returns:

Type	Description
Jacobian of total nucleon flux, shape ``(N, P)``.

Source code in src/globalsplinefit/model.py

def jacobian(
    self,
    energy_per_nucleon: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Jacobian of total nucleon flux (protons + neutrons).

    Parameters
    ----------
    energy_per_nucleon
        Energy per nucleon in GeV.
    target
        Target specification (group name or list of Z values).
    time_interval
        Time period for solar modulation, or "LIS" for unmodulated.
        ``None`` uses the default set during initialization.
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

    Returns
    -------
        Jacobian of total nucleon flux, shape ``(N, P)``.
    """
    jac_p, jac_n = self.p_and_n_jacobian(
        energy_per_nucleon,
        target,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    return jac_p + jac_n

p_and_n_flux(energy_per_nucleon: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Calculate separate proton and neutron flux from target (group or elements).

Parameters:

Name	Type	Description	Default
energy_per_nucleon	ArrayLike	Energy per nucleon in GeV. Can be scalar, list, or array.	required
target	Target	Target specification (same as GSFEnergy.flux).	required
time_interval	tuple[int, int] | str | None	Time period specification: - None: Uses the default_time_interval set during initialization - "LIS": Local Interstellar Spectrum (no modulation) - tuple: (start, end) in YYYYMM format	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. Nuclei with rigidity below this value are excluded. None uses the default.	None

Returns:

Type	Description
Array of shape (2, N) where N is the number of energy points:	[0, :]: Proton nucleon flux in particles/(m²·s·sr·GeV) [1, :]: Neutron nucleon flux in particles/(m²·s·sr·GeV)

Examples:

>>> nucleon_flux = model.p_and_n_flux([1, 10, 100], "He")
>>> proton_flux = nucleon_flux[0]  # 2 protons per He nucleus
>>> neutron_flux = nucleon_flux[1]  # 2 neutrons per He nucleus

Source code in src/globalsplinefit/model.py

def p_and_n_flux(
    self,
    energy_per_nucleon: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Calculate separate proton and neutron flux from target (group or elements).

    Parameters
    ----------
    energy_per_nucleon
        Energy per nucleon in GeV. Can be scalar, list, or array.
    target
        Target specification (same as GSFEnergy.flux).
    time_interval
        Time period specification:
        - None: Uses the default_time_interval set during initialization
        - "LIS": Local Interstellar Spectrum (no modulation)
        - tuple: (start, end) in YYYYMM format
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. Nuclei with rigidity
        below this value are excluded. None uses the default.

    Returns
    -------
        Array of shape (2, N) where N is the number of energy points:
        - [0, :]: Proton nucleon flux in particles/(m²·s·sr·GeV)
        - [1, :]: Neutron nucleon flux in particles/(m²·s·sr·GeV)

    Examples
    --------
        >>> nucleon_flux = model.p_and_n_flux([1, 10, 100], "He")
        >>> proton_flux = nucleon_flux[0]  # 2 protons per He nucleus
        >>> neutron_flux = nucleon_flux[1]  # 2 neutrons per He nucleus
    """
    time_interval = self._resolve_time_interval(time_interval)
    rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
    zlist, _group_leader = self._resolve_z(target)
    energy_per_nucleon = self._transform_energy_per_nucleon(
        energy_per_nucleon, target
    )

    flux = np.zeros((2, len(energy_per_nucleon)))
    for sid in self._target_sids(zlist):
        charge = sid[0]
        mass_scale = self.z_to_a[sid]
        nucleons = self.mass_number[sid]
        energy = energy_per_nucleon * mass_scale
        mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
        fl = self._element_flux(sid, energy, time_interval)
        flux[0] += fl * mass_scale * charge * mask
        flux[1] += fl * mass_scale * (nucleons - charge) * mask
    return flux

p_and_n_jacobian(energy_per_nucleon: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> tuple[np.ndarray, np.ndarray]

Jacobian of separate proton and neutron flux.

Parameters:

Name	Type	Description	Default
energy_per_nucleon	ArrayLike	Energy per nucleon in GeV.	required
target	Target	Target specification (group name or list of Z values).	required
time_interval	tuple[int, int] | str | None	Time period for solar modulation, or "LIS" for unmodulated. None uses the default set during initialization.	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. None uses the default.	None

Returns:

Type	Description
Tuple ``(jac_p, jac_n)`` of proton and neutron flux Jacobians,	each of shape (N, P) where P is the number of parameters.

Source code in src/globalsplinefit/model.py

def p_and_n_jacobian(
    self,
    energy_per_nucleon: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Jacobian of separate proton and neutron flux.

    Parameters
    ----------
    energy_per_nucleon
        Energy per nucleon in GeV.
    target
        Target specification (group name or list of Z values).
    time_interval
        Time period for solar modulation, or "LIS" for unmodulated.
        ``None`` uses the default set during initialization.
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

    Returns
    -------
        Tuple ``(jac_p, jac_n)`` of proton and neutron flux Jacobians,
        each of shape ``(N, P)`` where ``P`` is the number of parameters.
    """
    time_interval = self._resolve_time_interval(time_interval)
    rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
    zlist, _group_leader = self._resolve_z(target)
    energy_per_nucleon = self._transform_energy_per_nucleon(
        energy_per_nucleon, target
    )

    jac_p = 0.0
    jac_n = 0.0
    for sid in self._target_sids(zlist):
        charge = sid[0]
        mass_scale = self.z_to_a[sid]
        nucleons = self.mass_number[sid]
        energy = energy_per_nucleon * mass_scale
        mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
        j = self._element_flux_jacobian(sid, energy, time_interval)
        jac_p += j * (mass_scale * charge * mask)[:, np.newaxis]
        jac_n += j * (mass_scale * (nucleons - charge) * mask)[:, np.newaxis]
    return np.asarray(jac_p), np.asarray(jac_n)

p_and_n_covariance(target1: Target, target2: Target, energy_per_nucleon: ArrayLike, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> tuple[np.ndarray, np.ndarray]

Separate covariance matrices for proton and neutron flux.

Parameters:

Name	Type	Description	Default
target1	Target	Target specifications for the two targets being correlated.	required
target2	Target	Target specifications for the two targets being correlated.	required
energy_per_nucleon	ArrayLike	Energy per nucleon in GeV.	required
time_interval	tuple[int, int] | str | None	Time period for solar modulation. None uses the default.	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. None uses the default.	None

Returns:

Type	Description
Tuple ``(cov_pp, cov_nn)`` of the proton-proton and	neutron-neutron covariance matrices, each shape (N, N).

Source code in src/globalsplinefit/model.py

def p_and_n_covariance(
    self,
    target1: Target,
    target2: Target,
    energy_per_nucleon: ArrayLike,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Separate covariance matrices for proton and neutron flux.

    Parameters
    ----------
    target1, target2
        Target specifications for the two targets being correlated.
    energy_per_nucleon
        Energy per nucleon in GeV.
    time_interval
        Time period for solar modulation. ``None`` uses the default.
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

    Returns
    -------
        Tuple ``(cov_pp, cov_nn)`` of the proton-proton and
        neutron-neutron covariance matrices, each shape ``(N, N)``.
    """
    zlist1, leader1 = self._resolve_z(target1)
    zlist2, leader2 = self._resolve_z(target2)

    # Calculate jacobians using the helper method
    jac1_p, jac1_n = self.p_and_n_jacobian(
        energy_per_nucleon,
        target1,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    jac2_p, jac2_n = (
        (jac1_p, jac1_n)
        if zlist2 == zlist1
        else self.p_and_n_jacobian(
            energy_per_nucleon,
            target2,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
    )

    # Key the covariance by the leader SPECIES id (Z, A). _resolve_z returns a
    # charge; a charge with a single species has a cov alias under the bare int,
    # but a charge carrying >1 species (p + D at Z=1) does not — so resolve to
    # the leader sid, which is always a real cov key.
    cov_key = (
        self._leader_by_charge.get(leader1, leader1),
        self._leader_by_charge.get(leader2, leader2),
    )

    if cov_key in self.cov:
        return (
            self._propagate_cov(jac1_p, jac2_p, self.cov[cov_key]),
            self._propagate_cov(jac1_n, jac2_n, self.cov[cov_key]),
        )
    else:
        # Elements not in the main groups (H, He, O, Fe) carry no covariance
        energy_per_nucleon = np.atleast_1d(energy_per_nucleon)
        n_energies = len(energy_per_nucleon)
        zero_cov = np.zeros((n_energies, n_energies))
        return zero_cov, zero_cov

p_and_n_error(energy_per_nucleon: ArrayLike, target: Target, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Separate proton and neutron flux uncertainties (1-sigma).

Parameters:

Name	Type	Description	Default
energy_per_nucleon	ArrayLike	Energy per nucleon in GeV.	required
target	Target	Target specification (group name or list of Z values).	required
time_interval	tuple[int, int] | str | None	Time period for solar modulation. None uses the default.	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. None uses the default.	None

Returns:

Type	Description
Array of shape ``(2, N)``: row 0 is proton uncertainties,	row 1 is neutron uncertainties.

Source code in src/globalsplinefit/model.py

def p_and_n_error(
    self,
    energy_per_nucleon: ArrayLike,
    target: Target,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Separate proton and neutron flux uncertainties (1-sigma).

    Parameters
    ----------
    energy_per_nucleon
        Energy per nucleon in GeV.
    target
        Target specification (group name or list of Z values).
    time_interval
        Time period for solar modulation. ``None`` uses the default.
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

    Returns
    -------
        Array of shape ``(2, N)``: row 0 is proton uncertainties,
        row 1 is neutron uncertainties.
    """
    jac_p, jac_n = self.p_and_n_jacobian(
        energy_per_nucleon,
        target,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    block = self._covariance_block(target, target)
    if block is None:
        return np.zeros((2, jac_p.shape[0]))
    var_p = np.einsum("ni,ij,nj->n", jac_p, block, jac_p)
    var_n = np.einsum("ni,ij,nj->n", jac_n, block, jac_n)
    return np.sqrt(np.clip(np.vstack([var_p, var_n]), 0.0, None))

p_and_n_total_flux(energy_or_rigidity: ArrayLike, *, time_interval: tuple[int, int] | str | None = None, rigidity_cutoff: float | None = None) -> np.ndarray

Separate proton and neutron flux summed over all element groups.

Parameters:

Name	Type	Description	Default
energy_or_rigidity	ArrayLike	Energy per nucleon in GeV.	required
time_interval	tuple[int, int] | str | None	Time period for solar modulation. None uses the default.	None
rigidity_cutoff	float | None	Geomagnetic rigidity cutoff in GV. None uses the default.	None

Returns:

Type	Description
Array of shape ``(2, N)``: row 0 is total proton flux,	row 1 is total neutron flux, both in particles/(m²·s·sr·GeV).

Source code in src/globalsplinefit/model.py

def p_and_n_total_flux(
    self,
    energy_or_rigidity: ArrayLike,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Separate proton and neutron flux summed over all element groups.

    Parameters
    ----------
    energy_or_rigidity
        Energy per nucleon in GeV.
    time_interval
        Time period for solar modulation. ``None`` uses the default.
    rigidity_cutoff
        Geomagnetic rigidity cutoff in GV. ``None`` uses the default.

    Returns
    -------
        Array of shape ``(2, N)``: row 0 is total proton flux,
        row 1 is total neutron flux, both in particles/(m²·s·sr·GeV).
    """
    energy_or_rigidity = np.atleast_1d(energy_or_rigidity)
    # Initialize with proper shape for nucleon flux (2, N)
    total_flux = np.zeros((2, len(energy_or_rigidity)), dtype=float)

    for group in self.active_groups:
        group_flux = self.p_and_n_flux(
            energy_or_rigidity,
            group,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        total_flux += group_flux

    return total_flux

globalsplinefit.GSFKineticEnergyPerNucleon

Bases: GSFEnergyPerNucleon

GSF model for nucleon flux calculations using kinetic energy per nucleon.

Converts kinetic energy per nucleon to total energy per nucleon internally before using the GSFEnergyPerNucleon calculation methods.

Examples:

>>> from globalsplinefit import GSFKineticEnergyPerNucleon
>>> import numpy as np
>>> model = GSFKineticEnergyPerNucleon()
>>> kinetic_energy_per_nucleon = np.logspace(0, 2, 50)  # GeV/nucleon kinetic
>>> total_nucleons = model.flux(kinetic_energy_per_nucleon, "He")  # shape (N,)
>>> p_and_n = model.p_and_n_flux(kinetic_energy_per_nucleon, "He")

Source code in src/globalsplinefit/model.py

class GSFKineticEnergyPerNucleon(GSFEnergyPerNucleon):
    """GSF model for nucleon flux calculations using kinetic energy per nucleon.

    Converts kinetic energy per nucleon to total energy per nucleon internally
    before using the GSFEnergyPerNucleon calculation methods.

    Examples
    --------
        >>> from globalsplinefit import GSFKineticEnergyPerNucleon
        >>> import numpy as np
        >>> model = GSFKineticEnergyPerNucleon()
        >>> kinetic_energy_per_nucleon = np.logspace(0, 2, 50)  # GeV/nucleon kinetic
        >>> total_nucleons = model.flux(kinetic_energy_per_nucleon, "He")  # shape (N,)
        >>> p_and_n = model.p_and_n_flux(kinetic_energy_per_nucleon, "He")
    """

    def _transform_energy_per_nucleon(
        self,
        kinetic_energy_per_nucleon: ArrayLike,
        target: Target,  # noqa: ARG002
    ) -> np.ndarray:
        """Convert kinetic energy per nucleon to total energy per nucleon."""
        return (
            self._scale_energy(
                self._as_1d_values(
                    kinetic_energy_per_nucleon, "kinetic energy per nucleon"
                )
            )
            + NUCLEON_MASS_GEV
        )

Data Management

globalsplinefit.data_management.Parameters

Immutable bundle of loaded GSF model data.

Holds knots, spline parameters, covariance, nuclei and the solar-modulation table.

Parameters:

Name	Type	Description	Default
data_path	str or Path	Path to directory containing GSF data files. If None, uses :data:DEFAULT_VERSION.	None
version	str	Registered model version. Defaults to :data:DEFAULT_VERSION. See :data:MODEL_VERSIONS for provenance and available variants.	None

Source code in src/globalsplinefit/data_management.py

class Parameters:
    """Immutable bundle of loaded GSF model data.

    Holds knots, spline parameters, covariance, nuclei and the
    solar-modulation table.

    Parameters
    ----------
    data_path : str or Path, optional
        Path to directory containing GSF data files. If None, uses
        :data:`DEFAULT_VERSION`.
    version : str, optional
        Registered model version. Defaults to :data:`DEFAULT_VERSION`.
        See :data:`MODEL_VERSIONS` for provenance and available variants.
    """

    def __setattr__(self, name, value):
        if getattr(self, "_frozen", False):
            raise AttributeError("Parameters instances are immutable")
        object.__setattr__(self, name, value)

    def __init__(
        self,
        data_path: str | Path | None = None,
        version: str | None = None,
    ):
        self.data_path = self._setup_data_path(data_path, version)
        self._load_all_data()
        self._provenance = self._read_provenance()
        self._validate_loaded_data()
        self._freeze()

    @classmethod
    def for_model(
        cls,
        data_path: str | Path | None = None,
        version: str | None = None,
    ) -> "Parameters":
        """Return a shared immutable parameter bundle for model instances."""
        if data_path is not None:
            return cls(data_path=data_path, version=version)
        return cls._cached(resolve_version(version))

    @staticmethod
    @lru_cache(maxsize=32)
    def _cached(version: str) -> "Parameters":
        return Parameters(version=version)

    def _setup_data_path(
        self, data_path: str | Path | None, version: str | None
    ) -> Path:
        """Set up the data path, recording which version was resolved."""
        #: Resolved version name; None when loading an arbitrary data_path.
        self.version = None
        if version is not None and data_path is not None:
            raise ValueError("pass either version or data_path, not both")
        if version is not None or data_path is None:
            version_str = resolve_version(version)
            data_dir = Path(__file__).parent / "data" / version_str
            if not data_dir.exists():
                raise ValueError(
                    f"Version '{version_str}' not found. Available versions: {get_available_versions()}"
                )
            self.version = version_str
            return data_dir
        path = Path(data_path)
        if not path.exists():
            raise OSError(f"Data path does not exist: {path}")
        # An explicit path may still point at a packaged version directory.
        if (
            path.resolve().parent == (Path(__file__).parent / "data").resolve()
            and path.name in MODEL_VERSIONS
        ):
            self.version = path.name
        return path

    @property
    def provenance(self) -> dict:
        """Return a copy of the fit provenance metadata."""
        return deepcopy(self._provenance)

    def _read_provenance(self) -> dict:
        """Read and validate ``fit_result.json`` when present."""
        info: dict = {}
        meta_file = self.data_path / "fit_result.json"
        if meta_file.exists():
            try:
                loaded = json.loads(meta_file.read_text(encoding="utf-8"))
            except (OSError, ValueError) as exc:
                raise ValueError(f"invalid provenance file: {meta_file}") from exc
            if not isinstance(loaded, dict):
                raise ValueError(f"invalid provenance file: {meta_file}")
            meta = loaded.get("metadata", loaded)
            if not isinstance(meta, dict):
                raise ValueError(f"invalid provenance metadata: {meta_file}")
            for key in (
                "covering",
                "solar_modulation_source",
                "mixture",
                "mixture_components",
                "mixture_weights",
                "mixture_chi2_note",
                "slope_freeze",
                "norm_penalty",
                "written",
            ):
                if key in meta:
                    info[key] = meta[key]
        if self.version is not None:
            info["registry"] = dict(MODEL_VERSIONS[self.version])
            info["version"] = self.version
        return info

    def _load_all_data(self):
        """Load all required GSF data files."""
        self._load_nuclei_data()  # first: defines the species (Z, A) + groups
        self._load_knots()
        self._load_parameters()
        self._load_covariance()
        self._load_solar_modulation()
        self._load_subleading()
        self._load_reduced_pivots()
        self._calculate_flux_ratios()

    def _load_reduced_pivots(self):
        """Load the optional per-bundle ``reduced_pivots.dat`` grid.

        An absent file leaves ``reduced_pivots = None``.
        """
        self.reduced_pivots = None
        pivots_file = self.data_path / "reduced_pivots.dat"
        if not pivots_file.exists():
            return
        pivots = np.atleast_1d(np.loadtxt(pivots_file, dtype=float))
        if (
            pivots.ndim != 1
            or len(pivots) < 2
            or not np.all(np.isfinite(pivots))
            or np.any(pivots <= 0)
            or np.any(np.diff(pivots) <= 0)
        ):
            raise ValueError(
                f"invalid pivot table {pivots_file}: need at least two positive, "
                "finite, strictly increasing energies"
            )
        pivots.setflags(write=False)
        self.reduced_pivots = pivots

    def _validate_loaded_data(self) -> None:
        """Validate cross-file invariants before a model can use the bundle."""
        species = set(self.species)
        if not species:
            raise ValueError(f"no species found in {self.data_path / 'nuclei.dat'}")
        knot_species = {key for key in self.kx if isinstance(key, tuple)}
        parameter_species = {key for key in self.pars if isinstance(key, tuple)}
        if knot_species != species or parameter_species != species:
            missing_knots = sorted(species - knot_species)
            missing_parameters = sorted(species - parameter_species)
            raise ValueError(
                "incomplete model bundle: "
                f"missing knots for {missing_knots}, parameters for {missing_parameters}"
            )
        for sid in species:
            if sid[0] <= 0 or not np.isfinite(sid[1]) or sid[1] <= 0:
                raise ValueError(f"invalid species identifier {sid!r}")
            if self.z_ungroup[sid[0]] not in self.leader_sid:
                raise ValueError(f"species {sid!r} refers to an undefined group leader")
            knots = self.kx[sid][3:-3]
            if not np.all(np.isfinite(knots)) or np.any(np.diff(knots) <= 0):
                raise ValueError(f"knots for species {sid!r} are not finite/increasing")
            if len(self.pars[sid]) != self.npar[sid] + 4:
                raise ValueError(f"wrong parameter count for species {sid!r}")
            if not np.all(np.isfinite(self.pars[sid])):
                raise ValueError(f"non-finite parameters for species {sid!r}")
        for (sid1, sid2), block in self.cov.items():
            if not isinstance(sid1, tuple) or not isinstance(sid2, tuple):
                continue
            expected = (self.npar[sid1], self.npar[sid2])
            if block.shape != expected or not np.all(np.isfinite(block)):
                raise ValueError(
                    f"invalid covariance block {(sid1, sid2)!r}: "
                    f"expected {expected}, got {block.shape}"
                )
            reverse = self.cov.get((sid2, sid1))
            if reverse is None or not np.allclose(
                block, reverse.T, rtol=1e-12, atol=1e-12
            ):
                raise ValueError(f"asymmetric covariance blocks for {sid1!r}, {sid2!r}")

    def _freeze(self) -> None:
        """Make loaded numerical state safe to share between model instances."""
        array_maps = (self.kx, self.pars, self.cov, self.phi)
        for mapping in array_maps:
            for value in mapping.values():
                value.setflags(write=False)
        self.species = tuple(self.species)
        self._leaders = frozenset(self._leaders)
        self.z_to_sids = MappingProxyType(
            {charge: tuple(sids) for charge, sids in self.z_to_sids.items()}
        )
        self.z_group = MappingProxyType(
            {leader: tuple(charges) for leader, charges in self.z_group.items()}
        )
        for name in (
            "kx",
            "npar",
            "pars",
            "offset",
            "z_to_a",
            "mass_number",
            "z_ungroup",
            "leader_sid",
            "cov",
            "phi",
            "flux_ratio",
            "flux_slope",
        ):
            setattr(self, name, MappingProxyType(dict(getattr(self, name))))
        object.__setattr__(self, "_frozen", True)

    @staticmethod
    def _ncols(path) -> int:
        """Count whitespace-delimited fields in the first data line.

        Used to auto-detect the v1 vs v2 .dat column layout.
        """
        for line in Path(path).read_text().splitlines():
            if line.strip() and not line.startswith("#"):
                return len(line.replace(":", " ").split())
        return 0

    def _sid_for_z(self, z: int):
        """Return the unique species ID for a charge.

        Version 1 files contain one species per charge and use charge-only keys.
        """
        return self.z_to_sids[z][0]

    def _load_nuclei_data(self):
        """Load nuclear species and group membership.

        Species use ``(Z, A)`` IDs, allowing isotopes at the same charge. The
        lightest isotope at a leader charge is the group leader.
        """
        nuclei_file = self.data_path / "nuclei.dat"
        if not nuclei_file.exists():
            raise FileNotFoundError(f"Nuclei file not found: {nuclei_file}")

        data_array = np.atleast_1d(
            np.loadtxt(nuclei_file, dtype=[("z", int), ("a", float), ("l", int)])
        )
        rows = [(int(r["z"]), float(r["a"]), int(r["l"])) for r in data_array]
        species_rows = [(z, a) for z, a, _leader in rows]
        if len(species_rows) != len(set(species_rows)):
            raise ValueError(f"duplicate species in {nuclei_file}")
        if any(z <= 0 or not np.isfinite(a) or round(a) < z for z, a, _ in rows):
            raise ValueError(f"invalid charge or mass in {nuclei_file}")

        self.species = sorted({(z, a) for z, a, _ in rows})  # all species ids
        self.z_to_a = {(z, a): a for z, a, _ in rows}  # sid -> A (charge
        self.mass_number = {(z, a): int(round(a)) for z, a, _ in rows}
        #   aliases added by _add_charge_aliases; iterate self.species, not this)
        self.z_to_sids = {}  # charge -> [sids]
        for z, a, _ in rows:
            self.z_to_sids.setdefault(z, []).append((z, a))
        # Leader sid per group charge: lightest-A species at that charge.
        self.leader_sid = {}  # charge -> leader sid
        for z, a, leader in rows:
            if z == leader and (
                leader not in self.leader_sid or a < self.leader_sid[leader][1]
            ):
                self.leader_sid[leader] = (z, a)
        self._leaders = set(self.leader_sid.values())  # leader sids
        # Public group structure stays CHARGE-based (one entry per element charge)
        # so charge-indexed callers/tests are unchanged; isotopes are expanded
        # from a charge to its species ids (z_to_sids) inside the flux loops.
        self.z_group = {}  # int leader charge -> [member element charges]
        self.z_ungroup = {}  # element charge -> leader charge
        for z, _a, leader in rows:
            if z not in self.z_group.get(leader, []):
                self.z_group.setdefault(leader, []).append(z)
            self.z_ungroup[z] = leader

    def _load_knots(self):
        """Load spline knots.

        V2 lines are ``Z A: …``; v1 lines ``Z: …`` (keyed by
        the unique sid at that charge). Keys are species ids ``(Z, A)``.
        """
        knots_file = self.data_path / "knots.dat"
        if not knots_file.exists():
            raise FileNotFoundError(f"Knots file not found: {knots_file}")

        self.kx = {}
        self.npar = {}

        with open(knots_file) as f:
            for line in f:
                if line.startswith("#") or not line.strip():
                    continue
                head, k = line.split(":")
                toks = head.split()
                z = int(toks[0])
                sid = (z, float(toks[1])) if len(toks) >= 2 else self._sid_for_z(z)
                x = np.asarray([float(v) * np.log(10.0) for v in k.split()])
                self.npar[sid] = len(x) + 2
                # splev requires extended knot vector
                x = np.append((x[0], x[0], x[0]), x)
                x = np.append(x, (x[-1], x[-1], x[-1]))
                self.kx[sid] = x

    def _load_parameters(self):
        """Load spline parameters. v2: ``i Z A val``; v1: ``i Z val``."""
        params_file = self.data_path / "parameters.dat"
        if not params_file.exists():
            raise FileNotFoundError(f"Parameters file not found: {params_file}")

        self.pars = {}
        self.offset = {}
        v2 = self._ncols(params_file) >= 4
        dt = (
            [("i", int), ("z", int), ("a", float), ("val", float)]
            if v2
            else [("i", int), ("z", int), ("val", float)]
        )
        data_array = np.atleast_1d(np.loadtxt(params_file, dtype=dt))
        seen: dict[tuple[int, float], set[int]] = {}

        for r in data_array:
            i, z, val = int(r["i"]), int(r["z"]), float(r["val"])
            sid = (z, float(r["a"])) if v2 else self._sid_for_z(z)
            if sid not in self.pars:
                # 4 extra zeros at the end are needed by splev
                self.pars[sid] = np.zeros(self.npar[sid] + 4)
                self.offset[sid] = i
                seen[sid] = set()
            local_index = i - self.offset[sid]
            if local_index in seen[sid] or not 0 <= local_index < self.npar[sid]:
                raise ValueError(
                    f"invalid or duplicate parameter index {i} for {sid!r}"
                )
            seen[sid].add(local_index)
            self.pars[sid][local_index] = val

        for sid in self.species:
            if seen.get(sid) != set(range(self.npar[sid])):
                raise ValueError(f"incomplete parameter vector for species {sid!r}")

    def _load_covariance(self):
        """Load the parameter covariance.

        V2 rows contain ``i j Z1 A1 Z2 A2 value``; v1 rows contain
        ``i j Z1 Z2 value``.

        Keyed by species-id pairs ``((Z1,A1), (Z2,A2))``, group leaders only.
        The fit freezes the sub-leading splines, so uncertainty propagation
        resolves every target to its group leader (:meth:`_covariance_block`)
        and a sub-leading species inherits the leader's uncertainty. The 2026
        sets still carry per-species blocks for all 28 charges; they are read
        and dropped here rather than kept as an attractive nuisance.
        """
        cov_file = self.data_path / "covariance.dat"
        if not cov_file.exists():
            raise FileNotFoundError(f"Covariance file not found: {cov_file}")

        self.cov = {}
        v2 = self._ncols(cov_file) >= 7
        dt = (
            [
                ("i", int),
                ("j", int),
                ("z1", int),
                ("a1", float),
                ("z2", int),
                ("a2", float),
                ("val", float),
            ]
            if v2
            else [("i", int), ("j", int), ("z1", int), ("z2", int), ("val", float)]
        )
        data_array = np.atleast_1d(np.loadtxt(cov_file, dtype=dt))

        for r in data_array:
            i, j, val = int(r["i"]), int(r["j"]), float(r["val"])
            if v2:
                s1 = (int(r["z1"]), float(r["a1"]))
                s2 = (int(r["z2"]), float(r["a2"]))
            else:
                s1, s2 = self._sid_for_z(int(r["z1"])), self._sid_for_z(int(r["z2"]))
            if s1 not in self.npar or s2 not in self.npar:
                raise ValueError(
                    f"covariance references unknown species {s1!r}, {s2!r}"
                )
            if s1 not in self._leaders or s2 not in self._leaders:
                continue
            local_i = i - self.offset[s1]
            local_j = j - self.offset[s2]
            if not 0 <= local_i < self.npar[s1] or not 0 <= local_j < self.npar[s2]:
                raise ValueError(
                    f"covariance index {(i, j)!r} is outside the parameter blocks "
                    f"for {s1!r}, {s2!r}"
                )
            if (s1, s2) not in self.cov:
                self.cov[(s1, s2)] = np.zeros((self.npar[s1], self.npar[s2]))
            if (s2, s1) not in self.cov:
                self.cov[(s2, s1)] = np.zeros((self.npar[s2], self.npar[s1]))
            self.cov[(s1, s2)][local_i, local_j] = val
            self.cov[(s2, s1)][local_j, local_i] = val

    def _load_solar_modulation(self):
        """Load the monthly solar-modulation potential table (phi, MV).

        Precedence: a VERSION-LOCAL ``<data_path>/solar_modulation.dat``,
        else the shared Usoskin table at the package data root. The LIS is
        demodulated with a specific phi(t), so the fitted LIS must be
        re-modulated by the SAME potential to recover a flux at Earth: the
        GMD sets (2026 mixture and single-interpretation variants) store
        their Ghelfi-Maurin-Derome table alongside the parameters, while
        2017/2019/2025 and 2026-USO use the shared Usoskin table.
        """
        local = self.data_path / "solar_modulation.dat"
        shared = Path(__file__).parent / "data" / "solar_modulation.dat"
        solar_file = local if local.exists() else shared
        if not solar_file.exists():
            raise FileNotFoundError(f"Solar modulation file not found: {solar_file}")

        # Detect encoding from the byte-order mark (shared file is UTF-16-LE
        # with BOM; version-local files are plain UTF-8).
        with open(solar_file, "rb") as fh:
            encoding = "utf-16" if fh.read(2) == b"\xff\xfe" else "utf-8"

        # '#'-commented header of arbitrary length; data rows are numeric.
        data_array = np.loadtxt(solar_file, comments="#", encoding=encoding)

        self.phi = {}
        for row in data_array:
            months = row[1:13]  # Year, Jan..Dec[, Annual] -> keep the 12 months
            if np.isnan(months).any():  # drop incomplete years (e.g. Jan 1951)
                continue
            self.phi[int(row[0])] = months * 1e-3  # MV -> GV

    def _load_subleading(self):
        """Load optional subleading-species extrapolation parameters."""
        self._stored_sub = {}
        sub_file = self.data_path / "subleading.dat"
        if not sub_file.exists():
            return
        dt = [("z", int), ("a", float), ("norm", float), ("slope", float)]
        for r in np.atleast_1d(np.loadtxt(sub_file, dtype=dt)):
            sid = (int(r["z"]), float(r["a"]))
            values = (float(r["norm"]), float(r["slope"]))
            if sid not in self.species or sid in self._stored_sub:
                raise ValueError(f"invalid or duplicate species {sid!r} in {sub_file}")
            if not np.all(np.isfinite(values)) or values[0] < 0:
                raise ValueError(f"invalid extrapolation values for {sid!r}")
            self._stored_sub[sid] = values

    def _calculate_flux_ratios(self):
        """Calculate subleading-species flux ratios and slopes.

        Stored values take precedence. Otherwise the ratio is evaluated at the
        species' top knot and the slope is zero.
        """
        self.flux_ratio = {}
        self.flux_slope = {}

        for sid in self.species:
            leader_sid = self.leader_sid[self.z_ungroup[sid[0]]]
            xmax = self.kx[sid][-1]
            if sid != leader_sid:
                if sid in self._stored_sub:
                    ratio, slope = self._stored_sub[sid]
                else:
                    # Subleading species - ratio to its group leader
                    ratio = splev(
                        xmax, (self.kx[sid], self.pars[sid], SPLINE_DEGREE)
                    ) / splev(
                        xmax,
                        (self.kx[leader_sid], self.pars[leader_sid], SPLINE_DEGREE),
                    )
                    slope = 0.0
            else:
                ratio, slope = 1.0, 0.0
            self.flux_ratio[sid] = (leader_sid, ratio)
            self.flux_slope[sid] = slope
        self._add_charge_aliases()

    def _add_charge_aliases(self):
        """Add integer aliases for charges with exactly one species.

        Multi-isotope charges require explicit ``(Z, A)`` keys.
        """
        for z, sids in self.z_to_sids.items():
            if len(sids) != 1:
                continue
            s = sids[0]
            self.kx[z] = self.kx[s]
            self.npar[z] = self.npar[s]
            self.pars[z] = self.pars[s]
            self.offset[z] = self.offset[s]
            self.z_to_a[z] = self.z_to_a[s]
            self.mass_number[z] = self.mass_number[s]
            self.flux_ratio[z] = self.flux_ratio[s]
            self.flux_slope[z] = self.flux_slope[s]
        for s1, s2 in list(self.cov):
            z1, z2 = s1[0], s2[0]
            if len(self.z_to_sids[z1]) == 1 and len(self.z_to_sids[z2]) == 1:
                self.cov[(z1, z2)] = self.cov[(s1, s2)]

    def get_solar_cycle_24_interval(self) -> tuple[int, int]:
        """Get the time interval for Solar Cycle 24.

        Returns
        -------
        tuple[int, int]
            Time interval (start, end) in YYYYMM format.
        """
        return SOLAR_CYCLE_24_START[0], SOLAR_CYCLE_24_END[0]

    def get_solar_cycle_24_phi_average(self) -> float:
        """Mean of the monthly modulation potentials over Solar Cycle 24, in GV.

        Returns
        -------
        float
            Average solar modulation potential in GV.
        """
        start, end = self.get_solar_cycle_24_interval()
        phis = _collect_phi_values(self.phi, start, end)
        return float(np.mean(phis))

get_solar_cycle_24_interval() -> tuple[int, int]

Get the time interval for Solar Cycle 24.

Returns:

Type	Description
tuple[int, int]	Time interval (start, end) in YYYYMM format.

Source code in src/globalsplinefit/data_management.py

def get_solar_cycle_24_interval(self) -> tuple[int, int]:
    """Get the time interval for Solar Cycle 24.

    Returns
    -------
    tuple[int, int]
        Time interval (start, end) in YYYYMM format.
    """
    return SOLAR_CYCLE_24_START[0], SOLAR_CYCLE_24_END[0]

get_solar_cycle_24_phi_average() -> float

Mean of the monthly modulation potentials over Solar Cycle 24, in GV.

Returns:

Type	Description
float	Average solar modulation potential in GV.

Source code in src/globalsplinefit/data_management.py

def get_solar_cycle_24_phi_average(self) -> float:
    """Mean of the monthly modulation potentials over Solar Cycle 24, in GV.

    Returns
    -------
    float
        Average solar modulation potential in GV.
    """
    start, end = self.get_solar_cycle_24_interval()
    phis = _collect_phi_values(self.phi, start, end)
    return float(np.mean(phis))

Reduced representation

globalsplinefit.reduced.ReducedGSF

Pivot-component reduction of a GSF nucleon-flux model.

Parameters:

Name	Type	Description	Default
model	GSFEnergyPerNucleon or GSFKineticEnergyPerNucleon	Nucleon model to reduce. Default: GSFEnergyPerNucleon() using the 2026 set and Solar Cycle 24 average.	None
n_pivots	int	Number of log-spaced pivot energies per species. By default the model version's published pivot grid (model.params.reduced_pivots) is used.	None
energy_range	tuple of float	(E_min, E_max) of the pivot grid in GeV per nucleon. Default (1.0, 1e9) — beyond ~1e9 the heavy-group fluxes underflow and relative deviations lose meaning.	(1.0, 1000000000.0)
pivot_energies	array - like	Explicit pivot energies (overrides n_pivots/energy_range). Must be strictly increasing.	None
per_group	bool	If False (default), the species are the total proton and total neutron flux — sufficient for atmospheric-cascade applications, and much lower-dimensional. If True, every mass group contributes its own p and n species (8 species; for composition-sensitive users). The eight species come from four leader amplitude blocks, so their component covariance is rank-deficient (68 of 96 on the default grid) and :meth:penalty leaves that nullspace unconstrained. Use it to inspect composition correlations, not as a fit prior.	False
basis	(spline, hat)	Interpolation between pivots (in log-energy). "spline" (default) is a local cubic (Catmull-Rom) spline: smooth (C1) deformations, each component confined to its two neighboring intervals, at the cost of small side lobes there (the cardinal functions dip to ~-0.12). "hat" is piecewise-linear: strictly local and non-negative, but the deformations are kinked at the pivots. The component covariance is identical for both — only the behavior between pivots differs, and the coverage of the full model's variance is comparable.	"spline"
**kwargs		Passed to model methods (e.g. time_interval, rigidity_cutoff).	{}

Attributes:

Name	Type	Description
pivot_energies	(ndarray, shape(N))	The pivot grid.
species	list of str	Species row labels: ["p", "n"], or ["H_p", "H_n", "He_p", ...] (group then p/n) with per_group=True.
n_params	int	len(species) * N — the length of theta.
labels	list of str	One quotable name per component ("p_9TeV", "n_30PeV", ...), in theta order (species-major: all pivots of species 0, then species 1, ...).
cov	(ndarray, shape(n_params, n_params))	Exact covariance of theta (relative-flux units).
sigma	(ndarray, shape(n_params))	sqrt(diag(cov)) — the 1-sigma prior width of each component.
correlation	(ndarray, shape(n_params, n_params))	The correlation matrix of theta.

Examples:

>>> from globalsplinefit.reduced import ReducedGSF
>>> red = ReducedGSF()                    # 24 parameters (2 x 12)
>>> E = np.logspace(1, 6, 50)
>>> f_central = red.flux(E)               # (2, 50): [p, n], theta = 0
>>> theta = np.zeros(red.n_params)
>>> theta[3] = red.sigma[3]               # +1 sigma on one component
>>> f_varied = red.flux(E, theta)
>>> chi2 = red.penalty(theta)             # Gaussian prior contribution

Source code in src/globalsplinefit/reduced.py

class ReducedGSF:
    """Pivot-component reduction of a GSF nucleon-flux model.

    Parameters
    ----------
    model : GSFEnergyPerNucleon or GSFKineticEnergyPerNucleon, optional
        Nucleon model to reduce. Default: ``GSFEnergyPerNucleon()`` using the
        2026 set and Solar Cycle 24 average.
    n_pivots : int, optional
        Number of log-spaced pivot energies per species.  By default the
        model version's published pivot grid
        (``model.params.reduced_pivots``) is used.
    energy_range : tuple of float, optional
        ``(E_min, E_max)`` of the pivot grid in GeV per nucleon.
        Default ``(1.0, 1e9)`` — beyond ~1e9 the heavy-group fluxes
        underflow and relative deviations lose meaning.
    pivot_energies : array-like, optional
        Explicit pivot energies (overrides ``n_pivots``/``energy_range``).
        Must be strictly increasing.
    per_group : bool, optional
        If False (default), the species are the total proton and total
        neutron flux — sufficient for atmospheric-cascade applications, and
        much lower-dimensional.  If True, every mass group contributes its
        own p and n species (8 species; for composition-sensitive users).
        The eight species come from four leader amplitude blocks, so their
        component covariance is rank-deficient (68 of 96 on the default
        grid) and :meth:`penalty` leaves that nullspace unconstrained.  Use
        it to inspect composition correlations, not as a fit prior.
    basis : {"spline", "hat"}, optional
        Interpolation between pivots (in log-energy).  ``"spline"``
        (default) is a local cubic (Catmull-Rom) spline: smooth (C1)
        deformations, each component confined to its two neighboring
        intervals, at the cost of small side lobes there (the cardinal
        functions dip to ~-0.12).  ``"hat"`` is piecewise-linear: strictly
        local and non-negative, but the deformations are kinked at the
        pivots.  The component covariance is identical for both — only the
        behavior *between* pivots differs, and the coverage of the full
        model's variance is comparable.
    **kwargs
        Passed to model methods (e.g. ``time_interval``,
        ``rigidity_cutoff``).

    Attributes
    ----------
    pivot_energies : ndarray, shape (N,)
        The pivot grid.
    species : list of str
        Species row labels: ``["p", "n"]``, or ``["H_p", "H_n", "He_p", ...]``
        (group then p/n) with ``per_group=True``.
    n_params : int
        ``len(species) * N`` — the length of ``theta``.
    labels : list of str
        One quotable name per component (``"p_9TeV"``, ``"n_30PeV"``, ...),
        in ``theta`` order (species-major: all pivots of species 0, then
        species 1, ...).
    cov : ndarray, shape (n_params, n_params)
        Exact covariance of ``theta`` (relative-flux units).
    sigma : ndarray, shape (n_params,)
        ``sqrt(diag(cov))`` — the 1-sigma prior width of each component.
    correlation : ndarray, shape (n_params, n_params)
        The correlation matrix of ``theta``.

    Examples
    --------
    >>> from globalsplinefit.reduced import ReducedGSF
    >>> red = ReducedGSF()                    # 24 parameters (2 x 12)
    >>> E = np.logspace(1, 6, 50)
    >>> f_central = red.flux(E)               # (2, 50): [p, n], theta = 0
    >>> theta = np.zeros(red.n_params)
    >>> theta[3] = red.sigma[3]               # +1 sigma on one component
    >>> f_varied = red.flux(E, theta)
    >>> chi2 = red.penalty(theta)             # Gaussian prior contribution
    """

    def __init__(
        self,
        model: GSFEnergyPerNucleon | None = None,
        n_pivots: int | None = None,
        energy_range: tuple[float, float] = (1.0, 1e9),
        pivot_energies: ArrayLike | None = None,
        per_group: bool = False,
        basis: str = "spline",
        **kwargs,
    ):
        if model is None:
            model = GSFEnergyPerNucleon()
        if not isinstance(model, _NUCLEON_MODELS):
            raise TypeError(
                "ReducedGSF requires a nucleon model "
                "(GSFEnergyPerNucleon or GSFKineticEnergyPerNucleon), "
                f"got {type(model).__name__}"
            )
        if basis not in ("spline", "hat"):
            raise ValueError(f"basis must be 'spline' or 'hat', got {basis!r}")

        if pivot_energies is None:
            if n_pivots is None and energy_range == _DEFAULT_ENERGY_RANGE:
                pivot_energies = model.params.reduced_pivots
                if pivot_energies is None:
                    raise ValueError(
                        "this model bundle has no reduced_pivots.dat. Derive "
                        "a grid once with optimize_pivots(model, n_pivots=12) "
                        "and pass it via pivot_energies= (or store it as "
                        "reduced_pivots.dat in the bundle directory), or "
                        "request a log-spaced grid with n_pivots=."
                    )
            else:
                if n_pivots is not None and (
                    isinstance(n_pivots, bool)
                    or not isinstance(n_pivots, (int, np.integer))
                    or n_pivots < 2
                ):
                    raise ValueError("n_pivots must be an integer >= 2")
                if (
                    len(energy_range) != 2
                    or not np.all(np.isfinite(energy_range))
                    or energy_range[0] <= 0
                    or energy_range[0] >= energy_range[1]
                ):
                    raise ValueError(
                        "energy_range must be two positive increasing values"
                    )
                pivot_energies = np.logspace(
                    np.log10(energy_range[0]),
                    np.log10(energy_range[1]),
                    n_pivots if n_pivots is not None else 10,
                )
        pivot_energies = np.array(pivot_energies, dtype=float, copy=True)
        if (
            pivot_energies.ndim != 1
            or len(pivot_energies) < 2
            or not np.all(np.isfinite(pivot_energies))
            or np.any(pivot_energies <= 0)
            or np.any(np.diff(pivot_energies) <= 0)
        ):
            raise ValueError(
                "pivot_energies must be at least two positive, finite, increasing values"
            )

        self.model = model
        self.per_group = per_group
        self.pivot_energies = pivot_energies
        self.basis_type = basis
        self._log_pivots = np.log(pivot_energies)
        self._kwargs = dict(kwargs)

        n_piv = len(pivot_energies)
        jac_rel, cov_par, central, self.species = _relative_species_system(
            model, pivot_energies, per_group, **kwargs
        )
        covariance = jac_rel @ cov_par @ jac_rel.T
        covariance = 0.5 * (covariance + covariance.T)
        eigenvalues, eigenvectors = np.linalg.eigh(covariance)
        tolerance = (
            np.finfo(float).eps
            * max(covariance.shape)
            * max(float(np.max(np.abs(eigenvalues))), 1.0)
        )
        if np.min(eigenvalues) < -100 * tolerance:
            raise ValueError("reduced covariance is not positive semidefinite")
        self.cov = (eigenvectors * np.maximum(eigenvalues, 0.0)) @ eigenvectors.T
        self._central_pivots = central.reshape(len(self.species), n_piv)
        self._precision = None
        self._sample_factor = eigenvectors * np.sqrt(np.maximum(eigenvalues, 0.0))

        self.n_params = len(self.species) * n_piv
        self.labels = [
            f"{s}_{_format_energy(e)}" for s in self.species for e in pivot_energies
        ]
        for array in (
            self.pivot_energies,
            self._log_pivots,
            self.cov,
            self._central_pivots,
            self._sample_factor,
        ):
            array.setflags(write=False)

    # ------------------------------------------------------------------
    # Derived views
    # ------------------------------------------------------------------

    @property
    def sigma(self) -> np.ndarray:
        """1-sigma prior width of each component (relative flux units)."""
        return np.sqrt(np.diag(self.cov))

    @property
    def correlation(self) -> np.ndarray:
        """Correlation matrix of the components."""
        s = self.sigma
        denominator = np.outer(s, s)
        return np.divide(
            self.cov,
            denominator,
            out=np.zeros_like(self.cov),
            where=denominator > 0,
        )

    # ------------------------------------------------------------------
    # Basis
    # ------------------------------------------------------------------

    def basis(self, energy: ArrayLike) -> np.ndarray:
        """Interpolation basis H, shape (n_E, N).

        See the ``basis`` constructor parameter.
        """
        energy = self.model._as_1d_values(energy, "energy", positive=True)
        log_e = np.log(energy)
        return _interp_basis(self._log_pivots, log_e, self.basis_type)

    # ------------------------------------------------------------------
    # Model evaluation
    # ------------------------------------------------------------------

    def _effective_kwargs(self, overrides: dict) -> dict:
        for key, value in overrides.items():
            if key not in self._kwargs or self._kwargs[key] != value:
                raise ValueError(
                    f"cannot override {key!r} on an existing ReducedGSF; "
                    "construct a new reduction for different physical settings"
                )
        return self._kwargs

    def _central_flux(self, energy: np.ndarray, **kwargs) -> np.ndarray:
        """Central flux per species, shape (S, n_E)."""
        kw = self._effective_kwargs(kwargs)
        if self.per_group:
            return np.vstack(
                [self.model.p_and_n_flux(energy, g, **kw) for g in _GROUPS]
            )
        return np.asarray(self.model.p_and_n_total_flux(energy, **kw))

    def flux(
        self,
        energy: ArrayLike,
        theta: ArrayLike | None = None,
        **kwargs,
    ) -> np.ndarray:
        """Deformed flux for parameter vector ``theta``.

        Parameters
        ----------
        energy : array-like
            Energy per nucleon in GeV.
        theta : array-like, shape (n_params,), optional
            Component values (relative deviations at the pivots).
            ``None`` (default) returns the central flux.
        **kwargs
            Override kwargs for model methods.

        Returns
        -------
        flux : ndarray, shape (S, n_E)
            One row per entry of :attr:`species`.
        """
        energy = self.model._as_1d_values(energy, "energy", positive=True)
        central = self._central_flux(energy, **kwargs)
        if theta is None:
            return central
        theta = np.asarray(theta, dtype=float)
        if (
            theta.ndim != 1
            or theta.size != self.n_params
            or not np.all(np.isfinite(theta))
        ):
            raise ValueError(f"theta must be a finite vector of length {self.n_params}")
        theta = theta.reshape(len(self.species), len(self.pivot_energies))
        H = self.basis(energy)
        return central * (1.0 + theta @ H.T)

    def flux_jacobian(self, energy: ArrayLike, **kwargs) -> np.ndarray:
        """Compute the derivative of :meth:`flux` w.r.t. ``theta``.

        Returns
        -------
        jac : ndarray, shape (S, n_E, n_params)
            ``jac[s, i, j] = d flux[s, i] / d theta[j]``.  Species ``s``
            only responds to its own block of components.
        """
        energy = self.model._as_1d_values(energy, "energy", positive=True)
        central = self._central_flux(energy, **kwargs)
        H = self.basis(energy)
        n_s = len(self.species)
        n_piv = len(self.pivot_energies)
        jac = np.zeros((n_s, len(energy), self.n_params))
        for s in range(n_s):
            jac[s, :, s * n_piv : (s + 1) * n_piv] = central[s][:, np.newaxis] * H
        return jac

    def error(self, energy: ArrayLike, **kwargs) -> np.ndarray:
        """Absolute 1-sigma flux uncertainty of the reduced model.

        Exact at the pivot energies and interpolated with the selected basis.

        Returns
        -------
        sigma : ndarray, shape (S, n_E)
        """
        energy = self.model._as_1d_values(energy, "energy", positive=True)
        central = self._central_flux(energy, **kwargs)
        H = self.basis(energy)
        n_piv = len(self.pivot_energies)
        out = np.empty_like(central)
        for s in range(len(self.species)):
            block = self.cov[s * n_piv : (s + 1) * n_piv, s * n_piv : (s + 1) * n_piv]
            out[s] = np.sqrt(np.einsum("ij,jk,ik->i", H, block, H))
        return central * out

    # ------------------------------------------------------------------
    # Statistics
    # ------------------------------------------------------------------

    def sample(
        self,
        n_samples: int,
        rng: np.random.Generator | None = None,
    ) -> np.ndarray:
        """Draw component vectors ``theta ~ N(0, cov)``.

        Returns
        -------
        theta : ndarray, shape (n_samples, n_params)
        """
        if rng is None:
            rng = np.random.default_rng()
        if isinstance(n_samples, bool) or not isinstance(n_samples, (int, np.integer)):
            raise ValueError("n_samples must be a positive integer")
        if n_samples < 1:
            raise ValueError("n_samples must be a positive integer")
        return rng.standard_normal((n_samples, self.n_params)) @ self._sample_factor.T

    def penalty(self, theta: ArrayLike) -> float:
        """Gaussian penalty ``theta^T cov^-1 theta`` for a fit.

        Add this to the fit's chi-square to constrain the components to the
        GSF uncertainty.  The default (p, n) covariance is full rank.  With
        ``per_group=True`` it is not, and the pseudo-inverse charges nothing
        along the nullspace, so a fit is free to move there.
        """
        theta = np.asarray(theta, dtype=float)
        if (
            theta.ndim != 1
            or theta.size != self.n_params
            or not np.all(np.isfinite(theta))
        ):
            raise ValueError(f"theta must be a finite vector of length {self.n_params}")
        if self._precision is None:
            self._precision = np.linalg.pinv(self.cov, hermitian=True)
        return float(theta @ self._precision @ theta)

    # ------------------------------------------------------------------
    # Export
    # ------------------------------------------------------------------

    def to_dict(self) -> dict:
        """JSON-serializable description of the reduction."""
        return {
            "description": (
                "GSF reduced flux representation: theta are relative flux "
                "deviations at the pivot energies, interpolated in log(E) "
                f"with a {self.basis_type!r} basis; "
                "penalty = theta^T cov^-1 theta"
            ),
            "basis": self.basis_type,
            "model_version": self.model.version,
            "species": list(self.species),
            "pivot_energies_GeV": self.pivot_energies.tolist(),
            "labels": list(self.labels),
            "central_flux_at_pivots": self._central_pivots.tolist(),
            "cov": self.cov.tolist(),
        }

sigma: np.ndarray property

1-sigma prior width of each component (relative flux units).

correlation: np.ndarray property

Correlation matrix of the components.

basis(energy: ArrayLike) -> np.ndarray

Interpolation basis H, shape (n_E, N).

See the basis constructor parameter.

Source code in src/globalsplinefit/reduced.py

def basis(self, energy: ArrayLike) -> np.ndarray:
    """Interpolation basis H, shape (n_E, N).

    See the ``basis`` constructor parameter.
    """
    energy = self.model._as_1d_values(energy, "energy", positive=True)
    log_e = np.log(energy)
    return _interp_basis(self._log_pivots, log_e, self.basis_type)

flux(energy: ArrayLike, theta: ArrayLike | None = None, **kwargs) -> np.ndarray

Deformed flux for parameter vector theta.

Parameters:

Name	Type	Description	Default
energy	array - like	Energy per nucleon in GeV.	required
theta	(array - like, shape(n_params))	Component values (relative deviations at the pivots). None (default) returns the central flux.	None
**kwargs		Override kwargs for model methods.	{}

Returns:

Name	Type	Description
flux	(ndarray, shape(S, n_E))	One row per entry of :attr:species.

Source code in src/globalsplinefit/reduced.py

def flux(
    self,
    energy: ArrayLike,
    theta: ArrayLike | None = None,
    **kwargs,
) -> np.ndarray:
    """Deformed flux for parameter vector ``theta``.

    Parameters
    ----------
    energy : array-like
        Energy per nucleon in GeV.
    theta : array-like, shape (n_params,), optional
        Component values (relative deviations at the pivots).
        ``None`` (default) returns the central flux.
    **kwargs
        Override kwargs for model methods.

    Returns
    -------
    flux : ndarray, shape (S, n_E)
        One row per entry of :attr:`species`.
    """
    energy = self.model._as_1d_values(energy, "energy", positive=True)
    central = self._central_flux(energy, **kwargs)
    if theta is None:
        return central
    theta = np.asarray(theta, dtype=float)
    if (
        theta.ndim != 1
        or theta.size != self.n_params
        or not np.all(np.isfinite(theta))
    ):
        raise ValueError(f"theta must be a finite vector of length {self.n_params}")
    theta = theta.reshape(len(self.species), len(self.pivot_energies))
    H = self.basis(energy)
    return central * (1.0 + theta @ H.T)

flux_jacobian(energy: ArrayLike, **kwargs) -> np.ndarray

Compute the derivative of :meth:flux w.r.t. theta.

Returns:

Name	Type	Description
jac	(ndarray, shape(S, n_E, n_params))	jac[s, i, j] = d flux[s, i] / d theta[j]. Species s only responds to its own block of components.

Source code in src/globalsplinefit/reduced.py

def flux_jacobian(self, energy: ArrayLike, **kwargs) -> np.ndarray:
    """Compute the derivative of :meth:`flux` w.r.t. ``theta``.

    Returns
    -------
    jac : ndarray, shape (S, n_E, n_params)
        ``jac[s, i, j] = d flux[s, i] / d theta[j]``.  Species ``s``
        only responds to its own block of components.
    """
    energy = self.model._as_1d_values(energy, "energy", positive=True)
    central = self._central_flux(energy, **kwargs)
    H = self.basis(energy)
    n_s = len(self.species)
    n_piv = len(self.pivot_energies)
    jac = np.zeros((n_s, len(energy), self.n_params))
    for s in range(n_s):
        jac[s, :, s * n_piv : (s + 1) * n_piv] = central[s][:, np.newaxis] * H
    return jac

error(energy: ArrayLike, **kwargs) -> np.ndarray

Absolute 1-sigma flux uncertainty of the reduced model.

Exact at the pivot energies and interpolated with the selected basis.

Returns:

Name	Type	Description
sigma	(ndarray, shape(S, n_E))

Source code in src/globalsplinefit/reduced.py

def error(self, energy: ArrayLike, **kwargs) -> np.ndarray:
    """Absolute 1-sigma flux uncertainty of the reduced model.

    Exact at the pivot energies and interpolated with the selected basis.

    Returns
    -------
    sigma : ndarray, shape (S, n_E)
    """
    energy = self.model._as_1d_values(energy, "energy", positive=True)
    central = self._central_flux(energy, **kwargs)
    H = self.basis(energy)
    n_piv = len(self.pivot_energies)
    out = np.empty_like(central)
    for s in range(len(self.species)):
        block = self.cov[s * n_piv : (s + 1) * n_piv, s * n_piv : (s + 1) * n_piv]
        out[s] = np.sqrt(np.einsum("ij,jk,ik->i", H, block, H))
    return central * out

sample(n_samples: int, rng: np.random.Generator | None = None) -> np.ndarray

Draw component vectors theta ~ N(0, cov).

Returns:

Name	Type	Description
theta	(ndarray, shape(n_samples, n_params))

Source code in src/globalsplinefit/reduced.py

def sample(
    self,
    n_samples: int,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """Draw component vectors ``theta ~ N(0, cov)``.

    Returns
    -------
    theta : ndarray, shape (n_samples, n_params)
    """
    if rng is None:
        rng = np.random.default_rng()
    if isinstance(n_samples, bool) or not isinstance(n_samples, (int, np.integer)):
        raise ValueError("n_samples must be a positive integer")
    if n_samples < 1:
        raise ValueError("n_samples must be a positive integer")
    return rng.standard_normal((n_samples, self.n_params)) @ self._sample_factor.T

penalty(theta: ArrayLike) -> float

Gaussian penalty theta^T cov^-1 theta for a fit.

Add this to the fit's chi-square to constrain the components to the GSF uncertainty. The default (p, n) covariance is full rank. With per_group=True it is not, and the pseudo-inverse charges nothing along the nullspace, so a fit is free to move there.

Source code in src/globalsplinefit/reduced.py

def penalty(self, theta: ArrayLike) -> float:
    """Gaussian penalty ``theta^T cov^-1 theta`` for a fit.

    Add this to the fit's chi-square to constrain the components to the
    GSF uncertainty.  The default (p, n) covariance is full rank.  With
    ``per_group=True`` it is not, and the pseudo-inverse charges nothing
    along the nullspace, so a fit is free to move there.
    """
    theta = np.asarray(theta, dtype=float)
    if (
        theta.ndim != 1
        or theta.size != self.n_params
        or not np.all(np.isfinite(theta))
    ):
        raise ValueError(f"theta must be a finite vector of length {self.n_params}")
    if self._precision is None:
        self._precision = np.linalg.pinv(self.cov, hermitian=True)
    return float(theta @ self._precision @ theta)

to_dict() -> dict

JSON-serializable description of the reduction.

Source code in src/globalsplinefit/reduced.py

def to_dict(self) -> dict:
    """JSON-serializable description of the reduction."""
    return {
        "description": (
            "GSF reduced flux representation: theta are relative flux "
            "deviations at the pivot energies, interpolated in log(E) "
            f"with a {self.basis_type!r} basis; "
            "penalty = theta^T cov^-1 theta"
        ),
        "basis": self.basis_type,
        "model_version": self.model.version,
        "species": list(self.species),
        "pivot_energies_GeV": self.pivot_energies.tolist(),
        "labels": list(self.labels),
        "central_flux_at_pivots": self._central_pivots.tolist(),
        "cov": self.cov.tolist(),
    }

globalsplinefit.reduced.optimize_pivots(model: GSFEnergyPerNucleon | None = None, n_pivots: int = 12, energy_range: tuple[float, float] = (1.0, 1000000000.0), per_group: bool = False, basis: str = 'spline', n_grid: int = 300, n_restarts: int = 2, max_sweeps: int = 40, min_separation: float = 0.15, seed: int = 0, **kwargs) -> tuple[np.ndarray, float]

Optimize pivot placement for :class:ReducedGSF coverage.

Minimizes the worst-case mismatch between the reduced and the exact flux uncertainty over a dense log-energy grid,

max over (species, E) of |log(sigma_reduced / sigma_exact)|,

by coordinate exchange: candidate pivots are snapped to the dense grid, so the pivot covariance of every trial is a submatrix of one precomputed dense-grid covariance and each trial costs only linear algebra — no model re-evaluations. One pivot at a time is moved to its best available grid slot; sweeps repeat until no move improves the objective. The search restarts from the log-spaced grid and n_restarts random configurations, keeping the best result. The endpoint pivots stay fixed at energy_range.

Two guards keep the search honest: candidate pivots live on every second grid point, so the objective always samples between any two pivots (otherwise the exchange can hide spline ringing between its own evaluation points), and pivots must stay min_separation decades apart (near-duplicate pivots are statistically useless and make the cardinal spline ring violently).

Roughly half a minute with the defaults; scales as n_pivots * n_grid * n_restarts. The result is deterministic for a given seed.

Parameters:

Name	Type	Description	Default
model	GSFEnergyPerNucleon | None	As for :class:ReducedGSF.	None
energy_range	GSFEnergyPerNucleon | None	As for :class:ReducedGSF.	None
per_group	GSFEnergyPerNucleon | None	As for :class:ReducedGSF.	None
basis	GSFEnergyPerNucleon | None	As for :class:ReducedGSF.	None
**kwargs	GSFEnergyPerNucleon | None	As for :class:ReducedGSF.	None
n_pivots	int	Number of pivots to place. Default 12.	12
n_grid	int	Dense-grid resolution; pivots are quantized to this grid. Default 300 (about 0.03 decades over the default range).	300
n_restarts	int	Random restarts in addition to the log-spaced start. Default 2.	2
max_sweeps	int	Maximum exchange sweeps per start. Default 40.	40
min_separation	float	Minimum pivot separation in decades. Default 0.15.	0.15
seed	int	Seed for the restart configurations. Default 0.	0

Returns:

Name	Type	Description
pivot_energies	(ndarray, shape(n_pivots))	Optimized pivot grid — pass to ReducedGSF(pivot_energies=...) (with the same basis, per_group, and model kwargs).
max_ratio	float	Achieved worst-case coverage factor: sigma_reduced/sigma_exact lies within [1/max_ratio, max_ratio] on the dense grid.

Examples:

>>> pivots, worst = optimize_pivots(n_pivots=12)
>>> red = ReducedGSF(pivot_energies=pivots)

Source code in src/globalsplinefit/reduced.py

def optimize_pivots(
    model: GSFEnergyPerNucleon | None = None,
    n_pivots: int = 12,
    energy_range: tuple[float, float] = (1.0, 1e9),
    per_group: bool = False,
    basis: str = "spline",
    n_grid: int = 300,
    n_restarts: int = 2,
    max_sweeps: int = 40,
    min_separation: float = 0.15,
    seed: int = 0,
    **kwargs,
) -> tuple[np.ndarray, float]:
    """Optimize pivot placement for :class:`ReducedGSF` coverage.

    Minimizes the worst-case mismatch between the reduced and the exact
    flux uncertainty over a dense log-energy grid,

        max over (species, E) of |log(sigma_reduced / sigma_exact)|,

    by coordinate exchange: candidate pivots are snapped to the dense grid,
    so the pivot covariance of every trial is a submatrix of one
    precomputed dense-grid covariance and each trial costs only linear
    algebra — no model re-evaluations.  One pivot at a time is moved to its
    best available grid slot; sweeps repeat until no move improves the
    objective.  The search restarts from the log-spaced grid and
    ``n_restarts`` random configurations, keeping the best result.  The
    endpoint pivots stay fixed at ``energy_range``.

    Two guards keep the search honest: candidate pivots live on every
    *second* grid point, so the objective always samples between any two
    pivots (otherwise the exchange can hide spline ringing between its own
    evaluation points), and pivots must stay ``min_separation`` decades
    apart (near-duplicate pivots are statistically useless and make the
    cardinal spline ring violently).

    Roughly half a minute with the defaults; scales as
    ``n_pivots * n_grid * n_restarts``.  The result is deterministic for a
    given ``seed``.

    Parameters
    ----------
    model, energy_range, per_group, basis, **kwargs
        As for :class:`ReducedGSF`.
    n_pivots : int, optional
        Number of pivots to place.  Default 12.
    n_grid : int, optional
        Dense-grid resolution; pivots are quantized to this grid.
        Default 300 (about 0.03 decades over the default range).
    n_restarts : int, optional
        Random restarts in addition to the log-spaced start.  Default 2.
    max_sweeps : int, optional
        Maximum exchange sweeps per start.  Default 40.
    min_separation : float, optional
        Minimum pivot separation in decades.  Default 0.15.
    seed : int, optional
        Seed for the restart configurations.  Default 0.

    Returns
    -------
    pivot_energies : ndarray, shape (n_pivots,)
        Optimized pivot grid — pass to ``ReducedGSF(pivot_energies=...)``
        (with the same ``basis``, ``per_group``, and model kwargs).
    max_ratio : float
        Achieved worst-case coverage factor: ``sigma_reduced/sigma_exact``
        lies within ``[1/max_ratio, max_ratio]`` on the dense grid.

    Examples
    --------
    >>> pivots, worst = optimize_pivots(n_pivots=12)
    >>> red = ReducedGSF(pivot_energies=pivots)
    """
    if model is None:
        model = GSFEnergyPerNucleon()
    if not isinstance(model, _NUCLEON_MODELS):
        raise TypeError(
            "optimize_pivots requires a nucleon model "
            "(GSFEnergyPerNucleon or GSFKineticEnergyPerNucleon), "
            f"got {type(model).__name__}"
        )
    if basis not in ("spline", "hat"):
        raise ValueError(f"basis must be 'spline' or 'hat', got {basis!r}")
    if not 3 <= n_pivots < n_grid // 2:
        raise ValueError(f"n_pivots must be in [3, n_grid/2), got {n_pivots}")

    energies = np.logspace(np.log10(energy_range[0]), np.log10(energy_range[1]), n_grid)
    log_x = np.log(energies)
    jac_rel, cov_par, _, species = _relative_species_system(
        model, energies, per_group, **kwargs
    )
    S = jac_rel @ cov_par @ jac_rel.T
    sig2_exact = np.diag(S).reshape(len(species), n_grid)

    span_decades = np.log10(energy_range[1] / energy_range[0])
    min_sep = max(2, int(np.ceil(min_separation * (n_grid - 1) / span_decades)))
    candidates = np.arange(2, n_grid - 1, 2)
    if (n_pivots - 1) * min_sep >= n_grid - 1:
        raise ValueError(
            f"cannot place {n_pivots} pivots {min_separation} decades apart "
            f"within {span_decades:.2g} decades"
        )

    def objective(idx: np.ndarray) -> float:
        H = _interp_basis(log_x[idx], log_x, basis)
        worst = 0.0
        for s in range(len(species)):
            rows = s * n_grid + idx
            var = np.einsum("ij,jk,ik->i", H, S[np.ix_(rows, rows)], H)
            dev = np.abs(0.5 * np.log(var / sig2_exact[s]))
            worst = max(worst, float(dev.max()))
        return worst

    def exchange(idx: np.ndarray) -> tuple[np.ndarray, float]:
        best = objective(idx)
        for _ in range(max_sweeps):
            improved = False
            for j in range(1, n_pivots - 1):
                others = np.delete(idx, j)
                free = candidates[
                    np.min(np.abs(candidates[:, None] - others[None, :]), axis=1)
                    >= min_sep
                ]
                if len(free) == 0:
                    continue
                vals = [objective(np.sort(np.append(others, c))) for c in free]
                k = int(np.argmin(vals))
                if vals[k] < best - 1e-12:
                    idx = np.sort(np.append(others, free[k]))
                    best = vals[k]
                    improved = True
            if not improved:
                break
        return idx, best

    def random_start(rng: np.random.Generator) -> np.ndarray:
        # sequential rejection keeps the separation constraint satisfied
        idx = [0, n_grid - 1]
        pool = list(candidates)
        while len(idx) < n_pivots and pool:
            c = pool[rng.integers(len(pool))]
            if all(abs(c - i) >= min_sep for i in idx):
                idx.append(c)
            pool.remove(c)
        if len(idx) < n_pivots:
            raise ValueError("could not place pivots with the given separation")
        return np.sort(np.array(idx))

    rng = np.random.default_rng(seed)
    log_start = np.unique(np.round(np.linspace(0, n_grid - 1, n_pivots)).astype(int))
    starts = [log_start]
    starts += [random_start(rng) for _ in range(n_restarts)]

    results = [exchange(idx) for idx in starts]
    best_idx, best = min(results, key=lambda t: t[1])
    return energies[best_idx], float(np.exp(best))

Constants

GROUP_NAMES

Maps string group names ("H", "He", "O*", "Fe*", etc.) to their corresponding group leader atomic numbers.

from globalsplinefit.model import GSFBase
print(GSFBase.GROUP_NAMES)

SUBLEADING_SAT_LNR

Saturation rigidity for the sub-leading high-energy extrapolation, stored as \(\ln(R/\text{GV})\) — i.e. log(5e6), with \(R_{\text{sat}} = 5\) PV.

Above a sub-leading species' top knot the member-to-leader ratio is tilted by its fitted power-law slope and held constant beyond \(R_{\text{sat}}\): ratio(R) = norm * (min(R, R_sat)/Rmax)**slope. See Sub-leading Elements and High-Energy Extrapolation.

import numpy as np
from globalsplinefit.model import SUBLEADING_SAT_LNR
print(np.exp(SUBLEADING_SAT_LNR))  # 5e6 GV = 5 PV