API Reference

Complete API documentation for the GlobalSplineFit package.

Model Classes

The main interface to GlobalSplineFit functionality.

globalsplinefit.GSFEnergy

Bases: GSFBase

GSF model expecting total energy per nucleus in GeV.

This class implements the GSF model for cosmic ray calculations using total energy per nucleus as input. The energy should be the total kinetic + rest mass energy of the nucleus.

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 expecting total energy per nucleus in GeV.

    This class implements the GSF model for cosmic ray calculations using
    total energy per nucleus as input. The energy should be the total
    kinetic + rest mass energy of the nucleus.

    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: str | int | list[int],  # 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(np.atleast_1d(energy).astype(float))

    def flux(
        self,
        energy: ArrayLike,
        target: str | int | list[int],
        *,
        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)
        energy = self._transform_energy(energy, target)

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

    def jacobian(
        self,
        energy: ArrayLike,
        target: str | int | list[int],
        *,
        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)
        energy = self._transform_energy(energy, target)

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

    def covariance(
        self,
        target1: str | int | list[int],
        target2: str | int | list[int],
        energy: ArrayLike,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Calculate covariance matrix of flux."""
        zlist1, leader1 = self._resolve_z(target1)
        zlist2, leader2 = self._resolve_z(target2)
        # NOTE: pass the RAW input through — jacobian() applies
        # _transform_energy itself. Transforming here too double-applies the
        # kinetic->total rest-mass shift (GSFKineticEnergy) and squares the
        # energy_scale factor.

        # Calculate total jacobian for each group (sum over all elements)
        jac1 = self.jacobian(
            energy,
            target1,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        jac2 = jac1
        if zlist2 != zlist1:
            jac2 = self.jacobian(
                energy,
                target2,
                time_interval=time_interval,
                rigidity_cutoff=rigidity_cutoff,
            )

        # Use group leaders for covariance matrix lookup
        # 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))

        # Check if covariance matrix entry exists
        if cov_key in self.cov:
            return self._propagate_cov(jac1, jac2, self.cov[cov_key])
        else:
            n_energies = len(np.atleast_1d(energy))
            return np.zeros((n_energies, n_energies))

flux(energy: ArrayLike, target: str | int | list[int], *, 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	str | int | list[int]	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: str | int | list[int],
    *,
    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)
    energy = self._transform_energy(energy, target)

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

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

Calculate covariance matrix of flux.

Source code in src/globalsplinefit/model.py

def covariance(
    self,
    target1: str | int | list[int],
    target2: str | int | list[int],
    energy: ArrayLike,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Calculate covariance matrix of flux."""
    zlist1, leader1 = self._resolve_z(target1)
    zlist2, leader2 = self._resolve_z(target2)
    # NOTE: pass the RAW input through — jacobian() applies
    # _transform_energy itself. Transforming here too double-applies the
    # kinetic->total rest-mass shift (GSFKineticEnergy) and squares the
    # energy_scale factor.

    # Calculate total jacobian for each group (sum over all elements)
    jac1 = self.jacobian(
        energy,
        target1,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    jac2 = jac1
    if zlist2 != zlist1:
        jac2 = self.jacobian(
            energy,
            target2,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )

    # Use group leaders for covariance matrix lookup
    # 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))

    # Check if covariance matrix entry exists
    if cov_key in self.cov:
        return self._propagate_cov(jac1, jac2, self.cov[cov_key])
    else:
        n_energies = len(np.atleast_1d(energy))
        return np.zeros((n_energies, n_energies))

jacobian(energy: ArrayLike, target: str | int | list[int], *, 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: str | int | list[int],
    *,
    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)
    energy = self._transform_energy(energy, target)

    jac = 0.0
    for sid in self._target_sids(zlist):
        mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
        jac += (
            self._element_flux_jacobian(sid, 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: str | int | list[int]
    ) -> np.ndarray:
        """Convert kinetic energy per nucleus to total energy per nucleus."""
        kinetic_energy = self._scale_energy(np.atleast_1d(kinetic_energy).astype(float))
        zlist, leader = self._resolve_z(target)
        if len(zlist) == 1 and isinstance(zlist[0], tuple):
            # Explicit single-species target ("D", (1, 2.014)): its own mass.
            sid = zlist[0]
        else:
            # leader is a charge; resolve to its species id (z, a) so a
            # multi-species charge (p + D at Z=1) — which has no bare-int
            # z_to_a alias — still works.
            sid = self._leader_by_charge.get(leader, leader)
        rest_mass = self.z_to_a[sid] * NUCLEON_MASS_GEV
        return kinetic_energy + rest_mass

globalsplinefit.GSFRigidity

Bases: GSFBase

GSF model expecting rigidity in GV.

This class implements the GSF model for cosmic ray calculations using magnetic rigidity as input parameter. Rigidity is defined as R = pc/Z where p is momentum, c is speed of light, and Z is charge.

The model supports both Local Interstellar Spectrum (LIS) calculations and solar modulation effects using the force-field approximation to transform between Earth and interstellar rigidity spectra.

Note: the global energy_scale parameter is intentionally a no-op on this class — the rigidity/LIS model is the fit's low-energy anchor, and a rigidity scaling is not an energy scaling.

Examples:

>>> from globalsplinefit import GSFRigidity
>>> import numpy as np
>>> model = GSFRigidity()
>>> rigidity = np.logspace(0, 3, 100)  # 1 GV to 1 TV
>>> proton_flux_lis = model.flux(rigidity, "p")  # 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 expecting rigidity in GV.

    This class implements the GSF model for cosmic ray calculations using
    magnetic rigidity as input parameter. Rigidity is defined as
    R = pc/Z where p is momentum, c is speed of light, and Z is charge.

    The model supports both Local Interstellar Spectrum (LIS) calculations
    and solar modulation effects using the force-field approximation to
    transform between Earth and interstellar rigidity spectra.

    Note: the global ``energy_scale`` parameter is intentionally a no-op on
    this class — the rigidity/LIS model is the fit's low-energy anchor, and
    a rigidity scaling is not an energy scaling.

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

    # ---------- new helpers ----------
    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 (abs keeps unphysical R<0 non-negative)
            Lambda = (E_is / E) * (np.abs(rigidity) / (R_is + _ZERO_GUARD)) ** 3
        return R_is, Lambda

    # ---------- public API ----------
    def flux(
        self,
        rigidity: ArrayLike,
        target: str | int | list[int],
        *,
        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.

        Notes
        -----
            Solar modulation is now supported for rigidity-based calculations.
            The transformation uses the force-field approximation to convert
            between Earth and interstellar rigidity spectra.
        """
        time_interval = self._resolve_time_interval(time_interval)
        rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
        zlist, _ = self._resolve_z(target)
        rigidity = np.atleast_1d(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):  # loop version (safe, readable)
            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: str | int | list[int],
        *,
        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 = np.atleast_1d(rigidity)

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

                def _tilt(R, slope=slope, sid=sid):
                    # extrapolation tilt above the species' top knot
                    # (clamped below xmax and saturated at R_sat;
                    # 1.0 for slope 0)
                    if slope == 0.0:
                        return 1.0
                    xmax = self.kx[sid][-1]
                    dx_sat = max(SUBLEADING_SAT_LNR - xmax, 0.0)
                    with np.errstate(divide="ignore"):
                        return np.exp(slope * np.clip(
                            np.log(R) - xmax, 0.0, dx_sat))

                if phi == 0.0:
                    contrib = ratio * self._rigidity_flux_jacobian(
                        leading, rigidity)
                    if slope != 0.0:
                        contrib = contrib * _tilt(rigidity)[:, None]
                else:
                    R_is, fac = self._rigidity_phi_transform(sid, rigidity, phi)
                    contrib = (
                        ratio
                        * self._rigidity_flux_jacobian(leading, R_is)
                        * (fac * _tilt(R_is))[:, 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

    def covariance(
        self,
        target1: str | int | list[int],
        target2: str | int | list[int],
        rigidity: ArrayLike,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Calculate covariance matrix of flux."""
        zlist1, leader1 = self._resolve_z(target1)
        zlist2, leader2 = self._resolve_z(target2)

        jac1 = self.jacobian(
            rigidity,
            target1,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        jac2 = (
            jac1
            if zlist2 == zlist1
            else self.jacobian(
                rigidity,
                target2,
                time_interval=time_interval,
                rigidity_cutoff=rigidity_cutoff,
            )
        )

        # Use group leaders for covariance lookup
        # 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))

        # Check if covariance matrix entry exists
        if cov_key in self.cov:
            return self._propagate_cov(jac1, jac2, self.cov[cov_key])
        else:
            # If no covariance matrix entry exists, return zero covariance
            # This happens for elements not in the main groups (H, He, O, Fe)
            rigidity = np.atleast_1d(rigidity)
            n_rigidities = len(rigidity)
            return np.zeros((n_rigidities, n_rigidities))

flux(rigidity: ArrayLike, target: str | int | list[int], *, 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	str | int | list[int]	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.

Notes

Solar modulation is now supported for rigidity-based calculations.
The transformation uses the force-field approximation to convert
between Earth and interstellar rigidity spectra.

Source code in src/globalsplinefit/model.py

def flux(
    self,
    rigidity: ArrayLike,
    target: str | int | list[int],
    *,
    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.

    Notes
    -----
        Solar modulation is now supported for rigidity-based calculations.
        The transformation uses the force-field approximation to convert
        between Earth and interstellar rigidity spectra.
    """
    time_interval = self._resolve_time_interval(time_interval)
    rigidity_cutoff = self._resolve_rigidity_cutoff(rigidity_cutoff)
    zlist, _ = self._resolve_z(target)
    rigidity = np.atleast_1d(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):  # loop version (safe, readable)
        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

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

Calculate covariance matrix of flux.

Source code in src/globalsplinefit/model.py

def covariance(
    self,
    target1: str | int | list[int],
    target2: str | int | list[int],
    rigidity: ArrayLike,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Calculate covariance matrix of flux."""
    zlist1, leader1 = self._resolve_z(target1)
    zlist2, leader2 = self._resolve_z(target2)

    jac1 = self.jacobian(
        rigidity,
        target1,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    jac2 = (
        jac1
        if zlist2 == zlist1
        else self.jacobian(
            rigidity,
            target2,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
    )

    # Use group leaders for covariance lookup
    # 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))

    # Check if covariance matrix entry exists
    if cov_key in self.cov:
        return self._propagate_cov(jac1, jac2, self.cov[cov_key])
    else:
        # If no covariance matrix entry exists, return zero covariance
        # This happens for elements not in the main groups (H, He, O, Fe)
        rigidity = np.atleast_1d(rigidity)
        n_rigidities = len(rigidity)
        return np.zeros((n_rigidities, n_rigidities))

jacobian(rigidity: ArrayLike, target: str | int | list[int], *, 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: str | int | list[int],
    *,
    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 = np.atleast_1d(rigidity)

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

            def _tilt(R, slope=slope, sid=sid):
                # extrapolation tilt above the species' top knot
                # (clamped below xmax and saturated at R_sat;
                # 1.0 for slope 0)
                if slope == 0.0:
                    return 1.0
                xmax = self.kx[sid][-1]
                dx_sat = max(SUBLEADING_SAT_LNR - xmax, 0.0)
                with np.errstate(divide="ignore"):
                    return np.exp(slope * np.clip(
                        np.log(R) - xmax, 0.0, dx_sat))

            if phi == 0.0:
                contrib = ratio * self._rigidity_flux_jacobian(
                    leading, rigidity)
                if slope != 0.0:
                    contrib = contrib * _tilt(rigidity)[:, None]
            else:
                R_is, fac = self._rigidity_phi_transform(sid, rigidity, phi)
                contrib = (
                    ratio
                    * self._rigidity_flux_jacobian(leading, R_is)
                    * (fac * _tilt(R_is))[:, 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).

This class implements the GSF model for calculating nucleon (proton + neutron) flux from cosmic ray nuclei. Input energy is specified per nucleon, and the output separates proton and neutron contributions.

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).

    This class implements the GSF model for calculating nucleon (proton + neutron)
    flux from cosmic ray nuclei. Input energy is specified per nucleon, and the
    output separates proton and neutron contributions.

    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: str | int | list[int],  # noqa: ARG002
    ) -> np.ndarray:
        """Transform input energy per nucleon. Subclasses override for kinetic energy."""
        return self._scale_energy(np.atleast_1d(energy_per_nucleon).astype(float))

    def p_and_n_flux(
        self,
        energy_per_nucleon: ArrayLike,
        target: str | int | list[int],
        *,
        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):
            zi, ai = sid[0], self.z_to_a[sid]
            energy = energy_per_nucleon * ai
            mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
            fl = self._element_flux(sid, energy, time_interval)
            flux[0] += fl * ai * zi * mask  # protons
            flux[1] += fl * ai * (ai - zi) * mask  # neutrons
        return flux

    def flux(
        self,
        energy_per_nucleon: ArrayLike,
        target: str | int | list[int],
        *,
        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: str | int | list[int],
        *,
        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):
            zi, ai = sid[0], self.z_to_a[sid]
            energy = energy_per_nucleon * ai
            mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
            j = self._element_flux_jacobian(sid, energy, time_interval)
            jac_p += j * (ai * zi * mask)[:, np.newaxis]
            jac_n += j * (ai * (ai - zi) * mask)[:, np.newaxis]
        return np.asarray(jac_p), np.asarray(jac_n)

    def jacobian(
        self,
        energy_per_nucleon: ArrayLike,
        target: str | int | list[int],
        *,
        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: str | int | list[int],
        target2: str | int | list[int],
        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,
            )
        )

        # Use group leaders for covariance lookup
        # 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))

        # Check if covariance matrix entry exists
        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:
            # If no covariance matrix entry exists, return zero covariance
            # This happens for elements not in the main groups (H, He, O, Fe)
            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 covariance(
        self,
        target1: str | int | list[int],
        target2: str | int | list[int],
        energy_per_nucleon: ArrayLike,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """Covariance matrix of total nucleon flux (protons + neutrons).

        Uses the total Jacobian ``J = J_p + J_n``, which correctly includes all
        cross-terms ``Cov(P,P) + Cov(P,N) + Cov(N,P) + Cov(N,N)``.

        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
        -------
            Covariance matrix of shape ``(N, N)``.
        """
        zlist1, leader1 = self._resolve_z(target1)
        zlist2, leader2 = self._resolve_z(target2)

        # Use total Jacobian J = J_p + J_n directly.
        # Cov(P+N, P+N) = J_total @ C @ J_total.T which correctly includes
        # all four terms: Cov(P,P) + Cov(P,N) + Cov(N,P) + Cov(N,N).
        jac1 = self.jacobian(
            energy_per_nucleon,
            target1,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        jac2 = (
            jac1
            if zlist2 == zlist1
            else self.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, jac2, self.cov[cov_key])
        else:
            energy_per_nucleon = np.atleast_1d(energy_per_nucleon)
            n_energies = len(energy_per_nucleon)
            return np.zeros((n_energies, n_energies))

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

        Sums the nucleon flux from all active cosmic ray groups
        (protons, helium, oxygen group, iron group).

        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 total nucleon flux in particles/(m²·s·sr·GeV).
        """
        energy_or_rigidity = np.atleast_1d(energy_or_rigidity)
        # Initialize with proper shape for total nucleon flux (N,)
        total_flux = np.zeros(len(energy_or_rigidity), dtype=float)

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

        return total_flux

    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 total_error(
        self,
        energy_or_rigidity: ArrayLike,
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """1-sigma uncertainty of total nucleon flux from all element groups.

        Accounts for correlations between groups through the full covariance
        matrix.

        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 total nucleon flux uncertainties (1-sigma), shape ``(N,)``.
        """
        energy_or_rigidity = np.atleast_1d(energy_or_rigidity)
        n_points = len(energy_or_rigidity)

        # Initialize covariance matrix for total nucleon flux
        total_cov = np.zeros((n_points, n_points), dtype=float)

        # Sum covariances across all group pairs
        for l1 in self.active_groups:
            for l2 in self.active_groups:
                cov = self.covariance(
                    l1,
                    l2,
                    energy_or_rigidity,
                    time_interval=time_interval,
                    rigidity_cutoff=rigidity_cutoff,
                )
                total_cov += cov

        # Return uncertainties as (N,) array
        return np.sqrt(np.diag(total_cov))

    def p_and_n_error(
        self,
        energy_per_nucleon: ArrayLike,
        target: str | int | list[int],
        *,
        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.
        """
        cov_pp, cov_nn = self.p_and_n_covariance(
            target,
            target,
            energy_per_nucleon,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        return np.array([np.sqrt(np.diag(cov_pp)), np.sqrt(np.diag(cov_nn))])

    def error(
        self,
        energy_or_rigidity: ArrayLike,
        target: str | int | list[int],
        *,
        time_interval: tuple[int, int] | str | None = None,
        rigidity_cutoff: float | None = None,
    ) -> np.ndarray:
        """1-sigma uncertainty of total nucleon flux (protons + neutrons).

        Use `p_and_n_error` to get separate proton and neutron uncertainties.

        Parameters
        ----------
        energy_or_rigidity
            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 total nucleon flux uncertainties (1-sigma), shape ``(N,)``.
        """
        # Get covariance for total nucleon flux
        total_cov = self.covariance(
            target,
            target,
            energy_or_rigidity,
            time_interval=time_interval,
            rigidity_cutoff=rigidity_cutoff,
        )
        return np.sqrt(np.diag(total_cov))

flux(energy_per_nucleon: ArrayLike, target: str | int | list[int], *, 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	str | int | list[int]	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: str | int | list[int],
    *,
    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]

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

1-sigma uncertainty of total nucleon flux (protons + neutrons).

Use p_and_n_error to get separate proton and neutron uncertainties.

Parameters:

Name	Type	Description	Default
energy_or_rigidity	ArrayLike	Energy per nucleon in GeV.	required
target	str | int | list[int]	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 total nucleon flux uncertainties (1-sigma), shape ``(N,)``.

Source code in src/globalsplinefit/model.py

def error(
    self,
    energy_or_rigidity: ArrayLike,
    target: str | int | list[int],
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """1-sigma uncertainty of total nucleon flux (protons + neutrons).

    Use `p_and_n_error` to get separate proton and neutron uncertainties.

    Parameters
    ----------
    energy_or_rigidity
        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 total nucleon flux uncertainties (1-sigma), shape ``(N,)``.
    """
    # Get covariance for total nucleon flux
    total_cov = self.covariance(
        target,
        target,
        energy_or_rigidity,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    return np.sqrt(np.diag(total_cov))

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

Covariance matrix of total nucleon flux (protons + neutrons).

Uses the total Jacobian J = J_p + J_n, which correctly includes all cross-terms Cov(P,P) + Cov(P,N) + Cov(N,P) + Cov(N,N).

Parameters:

Name	Type	Description	Default
target1	str | int | list[int]	Target specifications for the two targets being correlated.	required
target2	str | int | list[int]	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
Covariance matrix of shape ``(N, N)``.

Source code in src/globalsplinefit/model.py

def covariance(
    self,
    target1: str | int | list[int],
    target2: str | int | list[int],
    energy_per_nucleon: ArrayLike,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """Covariance matrix of total nucleon flux (protons + neutrons).

    Uses the total Jacobian ``J = J_p + J_n``, which correctly includes all
    cross-terms ``Cov(P,P) + Cov(P,N) + Cov(N,P) + Cov(N,N)``.

    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
    -------
        Covariance matrix of shape ``(N, N)``.
    """
    zlist1, leader1 = self._resolve_z(target1)
    zlist2, leader2 = self._resolve_z(target2)

    # Use total Jacobian J = J_p + J_n directly.
    # Cov(P+N, P+N) = J_total @ C @ J_total.T which correctly includes
    # all four terms: Cov(P,P) + Cov(P,N) + Cov(N,P) + Cov(N,N).
    jac1 = self.jacobian(
        energy_per_nucleon,
        target1,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    jac2 = (
        jac1
        if zlist2 == zlist1
        else self.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, jac2, self.cov[cov_key])
    else:
        energy_per_nucleon = np.atleast_1d(energy_per_nucleon)
        n_energies = len(energy_per_nucleon)
        return np.zeros((n_energies, n_energies))

jacobian(energy_per_nucleon: ArrayLike, target: str | int | list[int], *, 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	str | int | list[int]	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: str | int | list[int],
    *,
    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

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

Total nucleon flux summed over all element groups.

Sums the nucleon flux from all active cosmic ray groups (protons, helium, oxygen group, iron group).

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 total nucleon flux in particles/(m²·s·sr·GeV).

Source code in src/globalsplinefit/model.py

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

    Sums the nucleon flux from all active cosmic ray groups
    (protons, helium, oxygen group, iron group).

    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 total nucleon flux in particles/(m²·s·sr·GeV).
    """
    energy_or_rigidity = np.atleast_1d(energy_or_rigidity)
    # Initialize with proper shape for total nucleon flux (N,)
    total_flux = np.zeros(len(energy_or_rigidity), dtype=float)

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

    return total_flux

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

1-sigma uncertainty of total nucleon flux from all element groups.

Accounts for correlations between groups through the full covariance matrix.

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 total nucleon flux uncertainties (1-sigma), shape ``(N,)``.

Source code in src/globalsplinefit/model.py

def total_error(
    self,
    energy_or_rigidity: ArrayLike,
    *,
    time_interval: tuple[int, int] | str | None = None,
    rigidity_cutoff: float | None = None,
) -> np.ndarray:
    """1-sigma uncertainty of total nucleon flux from all element groups.

    Accounts for correlations between groups through the full covariance
    matrix.

    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 total nucleon flux uncertainties (1-sigma), shape ``(N,)``.
    """
    energy_or_rigidity = np.atleast_1d(energy_or_rigidity)
    n_points = len(energy_or_rigidity)

    # Initialize covariance matrix for total nucleon flux
    total_cov = np.zeros((n_points, n_points), dtype=float)

    # Sum covariances across all group pairs
    for l1 in self.active_groups:
        for l2 in self.active_groups:
            cov = self.covariance(
                l1,
                l2,
                energy_or_rigidity,
                time_interval=time_interval,
                rigidity_cutoff=rigidity_cutoff,
            )
            total_cov += cov

    # Return uncertainties as (N,) array
    return np.sqrt(np.diag(total_cov))

p_and_n_flux(energy_per_nucleon: ArrayLike, target: str | int | list[int], *, 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	str | int | list[int]	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: str | int | list[int],
    *,
    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):
        zi, ai = sid[0], self.z_to_a[sid]
        energy = energy_per_nucleon * ai
        mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
        fl = self._element_flux(sid, energy, time_interval)
        flux[0] += fl * ai * zi * mask  # protons
        flux[1] += fl * ai * (ai - zi) * mask  # neutrons
    return flux

p_and_n_jacobian(energy_per_nucleon: ArrayLike, target: str | int | list[int], *, 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	str | int | list[int]	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: str | int | list[int],
    *,
    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):
        zi, ai = sid[0], self.z_to_a[sid]
        energy = energy_per_nucleon * ai
        mask = self._rigidity_cutoff_mask(sid, energy, rigidity_cutoff)
        j = self._element_flux_jacobian(sid, energy, time_interval)
        jac_p += j * (ai * zi * mask)[:, np.newaxis]
        jac_n += j * (ai * (ai - zi) * mask)[:, np.newaxis]
    return np.asarray(jac_p), np.asarray(jac_n)

p_and_n_covariance(target1: str | int | list[int], target2: str | int | list[int], 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	str | int | list[int]	Target specifications for the two targets being correlated.	required
target2	str | int | list[int]	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: str | int | list[int],
    target2: str | int | list[int],
    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,
        )
    )

    # Use group leaders for covariance lookup
    # 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))

    # Check if covariance matrix entry exists
    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:
        # If no covariance matrix entry exists, return zero covariance
        # This happens for elements not in the main groups (H, He, O, Fe)
        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: str | int | list[int], *, 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	str | int | list[int]	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: str | int | list[int],
    *,
    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.
    """
    cov_pp, cov_nn = self.p_and_n_covariance(
        target,
        target,
        energy_per_nucleon,
        time_interval=time_interval,
        rigidity_cutoff=rigidity_cutoff,
    )
    return np.array([np.sqrt(np.diag(cov_pp)), np.sqrt(np.diag(cov_nn))])

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: str | int | list[int],  # noqa: ARG002
    ) -> np.ndarray:
        """Convert kinetic energy per nucleon to total energy per nucleon."""
        return (
            self._scale_energy(np.atleast_1d(kinetic_energy_per_nucleon).astype(float))
            + NUCLEON_MASS_GEV
        )

Data Management

globalsplinefit.data_management.Parameters

Parameter loading and validation class for GSF model data.

This class handles loading and validation of all GSF model parameters including knots, spline parameters, covariance matrices, nuclear data, and solar modulation data.

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	Model version to use. "2026" (the default) is the promoted fit: the mixture covering -- an equal-weight combination of the Auger FD-2026 SIBYLL-2.3e and EPOS-LHC-R interpretations -- with the Ghelfi-Maurin-Derome modulation potential. "2026-USO" is the one sanctioned alternative: the same mixture fit with the Usoskin 2017 potential (about 65 MV lower, giving a 10-14% lower interstellar spectrum below 2 GV). "2026-UHE-S23e" and "2026-EPOS-LHCR" are the single-interpretation variants (SIBYLL-2.3e and EPOS-LHC-R, respectively). "2025", "2019" and "2017" are superseded historical releases, kept only so older work can be reproduced. See :data:MODEL_VERSIONS. If specified, overrides data_path and uses the corresponding package data directory.	None
solar_cycle_average_bins	int or None	Number of bins used to period-average the solar modulation. The monthly phi values of the requested interval are histogrammed into this many weighted representatives (each bin's mean phi, weighted by its month count) and the modulated flux is averaged over them -- the fitter's treatment. Default 6, worst-case 0.50% from the full monthly average over E >= 1 GeV (12 bins gives 0.10%). None averages over every month explicitly; 1 evaluates once at the mean phi (the former "approximate" mode, off by 5.2% at 1 GeV because flux is nonlinear in phi). See :class:~globalsplinefit.model.GSFBase for the full table.	6

Source code in src/globalsplinefit/data_management.py

class Parameters:
    """Parameter loading and validation class for GSF model data.

    This class handles loading and validation of all GSF model parameters
    including knots, spline parameters, covariance matrices, nuclear data,
    and solar modulation data.

    Parameters
    ----------
    data_path : str or Path, optional
        Path to directory containing GSF data files. If None, uses
        :data:`DEFAULT_VERSION`.
    version : str, optional
        Model version to use. "2026" (the default) is the promoted fit: the
        mixture covering -- an equal-weight combination of the Auger FD-2026
        SIBYLL-2.3e and EPOS-LHC-R interpretations -- with the
        Ghelfi-Maurin-Derome modulation potential. "2026-USO" is the one
        sanctioned alternative: the same mixture fit with the Usoskin 2017
        potential (about 65 MV lower, giving a 10-14% lower interstellar
        spectrum below 2 GV). "2026-UHE-S23e" and "2026-EPOS-LHCR" are the
        single-interpretation variants (SIBYLL-2.3e and EPOS-LHC-R,
        respectively). "2025", "2019" and "2017" are superseded
        historical releases, kept only so older work can be reproduced. See
        :data:`MODEL_VERSIONS`. If specified, overrides data_path and uses the
        corresponding package data directory.
    solar_cycle_average_bins : int or None, optional
        Number of bins used to period-average the solar modulation. The monthly
        phi values of the requested interval are histogrammed into this many
        weighted representatives (each bin's mean phi, weighted by its month
        count) and the modulated flux is averaged over them -- the fitter's
        treatment. Default 6, worst-case 0.50% from the full monthly average
        over E >= 1 GeV (12 bins gives 0.10%). ``None`` averages over every
        month explicitly; ``1`` evaluates once at the mean phi (the former
        "approximate" mode, off by 5.2% at 1 GeV because flux is nonlinear in
        phi). See :class:`~globalsplinefit.model.GSFBase` for the full table.
    """

    def __init__(
        self,
        data_path: str | Path | None = None,
        version: str | None = None,
        solar_cycle_average_bins: int | None = 6,
    ):
        """Initialize GSF parameters."""
        if solar_cycle_average_bins is not None and (
            not isinstance(solar_cycle_average_bins, (int, np.integer))
            or isinstance(solar_cycle_average_bins, bool)
            or solar_cycle_average_bins < 1
        ):
            raise ValueError(
                "solar_cycle_average_bins must be a positive int or None, got "
                f"{solar_cycle_average_bins!r}"
            )
        self.solar_cycle_average_bins = (
            None if solar_cycle_average_bins is None else int(solar_cycle_average_bins)
        )
        self.data_path = self._setup_data_path(data_path, version)
        self._load_all_data()

    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, or None when loading an arbitrary data_path.
        #: Unlike the constructor argument this is filled in for the default, so
        #: a model built with no arguments still reports what it loaded.
        self.version = None
        if version is not None:
            version_str = str(version).strip()
            if not version_str:
                raise ValueError(
                    f"Version '' not found. Available versions: {get_available_versions()}"
                )
            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
        if data_path is None:
            data_dir = Path(__file__).parent / "data" / DEFAULT_VERSION
            if not data_dir.exists():
                raise OSError(f"Default data directory not found: {data_dir}")
            self.version = DEFAULT_VERSION
            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:
        """What this parameter set actually is, read from its ``fit_result.json``.

        Returns the recorded fit provenance -- most importantly ``covering`` (the
        air-shower interpretation, e.g. the SIBYLL/EPOS mixture) and
        ``solar_modulation_source`` (GMD or USO), the two axes that distinguish
        the fits. For a registered version the :data:`MODEL_VERSIONS` entry is
        merged in under ``registry``.

        Returns an empty dict for legacy sets, which predate the metadata file.
        """
        info: dict = {}
        meta_file = self.data_path / "fit_result.json"
        if meta_file.exists():
            try:
                loaded = json.loads(meta_file.read_text())
            except (OSError, ValueError):
                loaded = {}
            meta = loaded.get("metadata", loaded)
            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._calculate_flux_ratios()

    @staticmethod
    def _ncols(path) -> int:
        """Whitespace-token count of the first data (non-#) line (``:`` -> space).
        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):
        """The unique species id ``(Z, A)`` at charge ``z`` (v1: one species per
        charge). Used to normalise v1 files, which key by charge only."""
        return self.z_to_sids[z][0]

    def _load_nuclei_data(self):
        """Load nuclear data from nuclei.dat. Species are identified by ``(Z, A)``
        so ISOTOPES (e.g. deuteron D at Z=1 alongside p) coexist. The leader of a
        charge group is the lightest-A species at the leader charge (p for the
        proton group; the single species otherwise)."""
        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]

        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
        #   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 l: the lightest-A species with z == l
        self.leader_sid = {}                                 # charge -> leader sid
        for z, a, l in rows:
            if z == l and (l not in self.leader_sid or a < self.leader_sid[l][1]):
                self.leader_sid[l] = (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, l in rows:
            if z not in self.z_group.get(l, []):
                self.z_group.setdefault(l, []).append(z)
            self.z_ungroup[z] = l

    def _load_knots(self):
        """Load knot data. 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.log([10 ** float(v) 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))

        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
            self.pars[sid][i - self.offset[sid]] = val

    def _load_covariance(self):
        """Load covariance. v2: ``i j Z1 A1 Z2 A2 val``; v1: ``i j Z1 Z2 val``.
        Keyed by species-id pairs ``((Z1,A1), (Z2,A2))``."""
        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, 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)][i - self.offset[s1], j - self.offset[s2]] = val
            self.cov[(s2, s1)][j - self.offset[s2], i - self.offset[s1]] = val

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

        Precedence: a VERSION-LOCAL ``<data_path>/solar_modulation.dat`` if the
        parameter set ships its own potential, else the shared table at the
        package data root. This matters because the LIS is demodulated with a
        specific phi(t): the fitted LIS must be re-modulated by the SAME
        potential to recover a flux at Earth. Sets whose LIS was demodulated
        with a non-default potential (e.g. 2026, Ghelfi-Maurin-Derome) ship
        their table alongside the parameters; legacy sets (2017/2019/2025) and
        the Usoskin variant (2026-USO) fall back to the shared Usoskin table.

        Robust to both the shared file (UTF-16-BOM, Usoskin) and version-local
        UTF-8 files: encoding is detected from the byte-order mark, header lines
        are '#'-commented, and the columns are Year followed by the 12 monthly
        values (a trailing Annual column, present in the Usoskin file, is
        ignored). Rows with any NaN month (e.g. Jan 1951 in the Usoskin table)
        are dropped.
        """
        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 the OPTIONAL subleading.dat extrapolation table: per sub-leading
        species ``(Z, A, norm, slope)`` with ratio(R > Rmax) = norm *
        (R/Rmax)**slope above the species' top knot Rmax. Absent in pre-slope
        bundles (2017/2019/2025) -> empty table, and the constant-ratio
        extrapolation is recomputed from the splines exactly as before."""
        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)):
            self._stored_sub[(int(r["z"]), float(r["a"]))] = (
                float(r["norm"]), float(r["slope"]))

    def _calculate_flux_ratios(self):
        """Calculate flux ratios (and extrapolation slopes) for subleading
        species (keyed by species id). When subleading.dat supplied stored
        (norm, slope) values, those are used verbatim — they are the values the
        fit itself used; otherwise the ratio is recomputed from the splines at
        the species' top knot and the slope defaults to 0 (the historical
        constant-ratio extrapolation)."""
        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):
        """Expose the per-species dicts under a bare integer charge as well, for
        every charge with a SINGLE species (all charges in a v1 model). Lets
        charge-indexed access (existing callers/tests) coexist with (Z, A) keys.
        ``species`` (not ``z_to_a``) is the species iterator, so aliasing z_to_a
        is safe. Multi-species charges (e.g. Z=1 with p+D) get no int alias."""
        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.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 tuple (start, end) in YYYYMM format for Solar Cycle 24
            (December 2008 to December 2019)
        """
        return SOLAR_CYCLE_24_START[0], SOLAR_CYCLE_24_END[0]

    def get_solar_cycle_24_phi_average(self) -> float:
        """Get the average solar modulation potential for Solar Cycle 24.

        This method implements the approximate averaging by calculating the
        mean of all monthly phi values during Solar Cycle 24.

        Returns
        -------
        float
            Average solar modulation potential in GV for Solar Cycle 24
        """
        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 tuple (start, end) in YYYYMM format for Solar Cycle 24 (December 2008 to December 2019)

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 tuple (start, end) in YYYYMM format for Solar Cycle 24
        (December 2008 to December 2019)
    """
    return SOLAR_CYCLE_24_START[0], SOLAR_CYCLE_24_END[0]

get_solar_cycle_24_phi_average() -> float

Get the average solar modulation potential for Solar Cycle 24.

This method implements the approximate averaging by calculating the mean of all monthly phi values during Solar Cycle 24.

Returns:

Type	Description
float	Average solar modulation potential in GV for Solar Cycle 24

Source code in src/globalsplinefit/data_management.py

def get_solar_cycle_24_phi_average(self) -> float:
    """Get the average solar modulation potential for Solar Cycle 24.

    This method implements the approximate averaging by calculating the
    mean of all monthly phi values during Solar Cycle 24.

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

PCA

globalsplinefit.pca.HybridPCA

Hybrid PCA dimensionality reduction for any GSF model.

Provides the same interface as the underlying model (flux, error, covariance, jacobian) but uses the low-rank + per-energy-block approximation for uncertainty propagation. Flux values are delegated unchanged to the underlying model; only the uncertainty methods use the reduction.

Parameters:

Name	Type	Description	Default
model	GSFBase	Any GSF model instance (GSFEnergy, GSFRigidity, GSFEnergyPerNucleon, etc.).	required
n_components	int	Number of retained components. Default 8. Fixed-energy observables (group bands, all-particle flux, total nucleon flux) are exact on the reference grid for ANY value; n_components controls the fidelity of cross-energy correlations and of smooth samples.	8
energy_grid	array - like	Reference energy grid for the decomposition. Default: 300 log-spaced points from 1 to 10^11 GeV.	None
gauge	(correlation, covariance)	Metric in which the low-rank factor is extracted. The default "correlation" diagonalizes the correlation matrix of the relative flux covariance (every grid row weighted equally); "covariance" diagonalizes the covariance itself (classic PCA, dominated by the rows with the largest relative variance).	"correlation"
**kwargs		Passed to model methods (e.g. time_interval, rigidity_cutoff).	{}

Attributes:

Name	Type	Description
L_param	(ndarray, shape(n_params, n_components))	Parameter-space low-rank factor.
B	(ndarray, shape(n_ref, n_sub, n_sub))	Per-energy cross-group residual blocks on the reference grid (n_sub = 4, or 8 for nucleon models: p/n per group), PSD-projected.
n_components	int	Number of retained components.
variance_explained	float	Fraction of the total variance in the chosen gauge captured by the low-rank part (in the correlation gauge: average fraction of the per-row variance).

Examples:

>>> from globalsplinefit import GSFEnergyPerNucleon
>>> from globalsplinefit.pca import HybridPCA
>>> gsf = GSFEnergyPerNucleon()          # promoted default (2026)
>>> pca = HybridPCA(gsf, n_components=8)
>>> E = np.logspace(1, 6, 50)
>>> flux = pca.flux(E, "p")             # same as gsf.flux(E, "p")
>>> sigma = pca.error(E, "p")           # absolute 1-sigma (like gsf.error)
>>> cov = pca.covariance("p", "p", E)   # absolute covariance matrix

Source code in src/globalsplinefit/pca.py

class HybridPCA:
    """Hybrid PCA dimensionality reduction for any GSF model.

    Provides the same interface as the underlying model (``flux``, ``error``,
    ``covariance``, ``jacobian``) but uses the low-rank + per-energy-block
    approximation for uncertainty propagation.  Flux values are delegated
    unchanged to the underlying model; only the uncertainty methods use the
    reduction.

    Parameters
    ----------
    model : GSFBase
        Any GSF model instance (GSFEnergy, GSFRigidity, GSFEnergyPerNucleon, etc.).
    n_components : int, optional
        Number of retained components. Default 8.  Fixed-energy observables
        (group bands, all-particle flux, total nucleon flux) are exact on the
        reference grid for ANY value; ``n_components`` controls the fidelity
        of cross-energy correlations and of smooth samples.
    energy_grid : array-like, optional
        Reference energy grid for the decomposition.
        Default: 300 log-spaced points from 1 to 10^11 GeV.
    gauge : {"correlation", "covariance"}, optional
        Metric in which the low-rank factor is extracted.  The default
        "correlation" diagonalizes the correlation matrix of the relative
        flux covariance (every grid row weighted equally); "covariance"
        diagonalizes the covariance itself (classic PCA, dominated by the
        rows with the largest relative variance).
    **kwargs
        Passed to model methods (e.g. ``time_interval``, ``rigidity_cutoff``).

    Attributes
    ----------
    L_param : ndarray, shape (n_params, n_components)
        Parameter-space low-rank factor.
    B : ndarray, shape (n_ref, n_sub, n_sub)
        Per-energy cross-group residual blocks on the reference grid
        (n_sub = 4, or 8 for nucleon models: p/n per group), PSD-projected.
    n_components : int
        Number of retained components.
    variance_explained : float
        Fraction of the total variance in the chosen gauge captured by the
        low-rank part (in the correlation gauge: average fraction of the
        per-row variance).

    Examples
    --------
    >>> from globalsplinefit import GSFEnergyPerNucleon
    >>> from globalsplinefit.pca import HybridPCA
    >>> gsf = GSFEnergyPerNucleon()          # promoted default (2026)
    >>> pca = HybridPCA(gsf, n_components=8)
    >>> E = np.logspace(1, 6, 50)
    >>> flux = pca.flux(E, "p")             # same as gsf.flux(E, "p")
    >>> sigma = pca.error(E, "p")           # absolute 1-sigma (like gsf.error)
    >>> cov = pca.covariance("p", "p", E)   # absolute covariance matrix
    """

    def __init__(
        self,
        model: GSFBase,
        n_components: int = 8,
        energy_grid: np.ndarray | list | float | None = None,
        gauge: str = "correlation",
        **kwargs,
    ):
        if gauge not in ("correlation", "covariance"):
            raise ValueError(
                f"gauge must be 'correlation' or 'covariance', got {gauge!r}"
            )
        if energy_grid is None:
            energy_grid = np.logspace(np.log10(1.0), 11, 300)
        energy_grid = np.atleast_1d(np.asarray(energy_grid, dtype=float))

        self.model = model
        self.n_components = n_components
        self.gauge = gauge
        self._energy_grid = energy_grid
        self._kwargs = kwargs
        self._is_nucleon_model = isinstance(model, _NUCLEON_MODELS)

        # Build the stacked system on the reference grid
        (jac_stack, cov_stack, par_stack, central_flux, group_leaders, npar_trimmed) = (
            _build_stacked_system(model, energy_grid, **kwargs)
        )

        self.par_stack = par_stack
        self.central_flux = central_flux
        self._group_leaders = group_leaders
        self._npar_trimmed = npar_trimmed

        # Mask out rows with zero central flux (e.g. neutron flux from H group)
        nonzero = central_flux != 0.0
        self._nonzero_mask = nonzero

        # Relative Jacobian: J_rel = J / f_central (only on non-zero rows)
        rel_jac = np.zeros_like(jac_stack)
        rel_jac[nonzero] = jac_stack[nonzero] / central_flux[nonzero, np.newaxis]

        # Relative flux covariance
        flux_rel_cov = rel_jac @ cov_stack @ rel_jac.T

        # Low-rank factor in the chosen gauge
        k = n_components
        if gauge == "correlation":
            s = np.sqrt(np.maximum(np.diag(flux_rel_cov), 0.0))
            snz = s > 0
            corr = np.zeros_like(flux_rel_cov)
            corr[np.ix_(snz, snz)] = flux_rel_cov[np.ix_(snz, snz)] / np.outer(
                s[snz], s[snz]
            )
            eigvals, eigvecs = np.linalg.eigh(corr)
            idx = np.argsort(eigvals)[::-1]
            eigvals = eigvals[idx]
            eigvecs = eigvecs[:, idx]
            # un-whiten the retained factor back to covariance scale
            L = (s[:, np.newaxis] * eigvecs[:, :k]) * np.sqrt(
                np.maximum(eigvals[:k], 0.0)
            )
        else:
            eigvals, eigvecs = np.linalg.eigh(flux_rel_cov)
            idx = np.argsort(eigvals)[::-1]
            eigvals = eigvals[idx]
            eigvecs = eigvecs[:, idx]
            L = eigvecs[:, :k] * np.sqrt(np.maximum(eigvals[:k], 0.0))

        self._eigvals = eigvals

        # Project to parameter space via least-squares (robust to rank-deficiency
        # when the energy grid has fewer points than parameters)
        self.L_param, _, _, _ = np.linalg.lstsq(rel_jac, L, rcond=None)

        # Store layout info
        self._ref_log_energy = np.log(energy_grid)
        self._n_ref = len(energy_grid)
        self._n_groups = len(_GROUPS)
        self._n_sub_per_group = 2 if self._is_nucleon_model else 1
        self._n_sub = self._n_groups * self._n_sub_per_group

        # Precompute row layout: rows_per_group for the reference grid
        self._rows_per_group_ref = (
            2 * self._n_ref if self._is_nucleon_model else self._n_ref
        )

        # Per-energy cross-group residual blocks: everything the low-rank part
        # misses BETWEEN sub-rows (group x p/n) at the SAME energy, kept
        # exactly (PSD-projected).  Cross-energy residuals are dropped.
        L_grid = rel_jac @ self.L_param
        n_ref = self._n_ref
        nsub = self._n_sub
        B = np.empty((n_ref, nsub, nsub))
        sub_offsets = self._sub_row_offsets(n_ref)
        for i in range(n_ref):
            rows = sub_offsets + i
            resid = flux_rel_cov[np.ix_(rows, rows)] - L_grid[rows] @ L_grid[rows].T
            B[i] = _nearest_psd(0.5 * (resid + resid.T))
        self.B = B

    def _sub_row_offsets(self, n_e: int) -> np.ndarray:
        """Row offset of each sub-row (group-major, p before n) for n_e energies."""
        rpg = self._n_sub_per_group * n_e
        return np.array(
            [
                ig * rpg + q * n_e
                for ig in range(self._n_groups)
                for q in range(self._n_sub_per_group)
            ]
        )

    # ------------------------------------------------------------------
    # Backward-compatible views of the block correction
    # ------------------------------------------------------------------

    @property
    def D(self) -> np.ndarray:
        """Diagonal of the block correction in stacked-row layout.

        Equivalent to the scalar diagonal correction of the original
        construction (relative variance not carried by the low-rank part).
        """
        n_ref = self._n_ref
        D = np.empty(self._n_sub * n_ref)
        offsets = self._sub_row_offsets(n_ref)
        for s, off in enumerate(offsets):
            D[off : off + n_ref] = self.B[:, s, s]
        return D

    @property
    def D_pn(self) -> np.ndarray | None:
        """p-n cross-term correction per group (nucleon models), from the block."""
        if not self._is_nucleon_model:
            return None
        n_ref = self._n_ref
        D_pn = np.empty(self._n_groups * n_ref)
        for ig in range(self._n_groups):
            D_pn[ig * n_ref : (ig + 1) * n_ref] = self.B[:, 2 * ig, 2 * ig + 1]
        return D_pn

    @property
    def variance_explained(self) -> float:
        """Fraction of gauge-metric variance captured by the low-rank part.

        In the correlation gauge this is the average fraction of the
        per-row variance resolved by the components; in the covariance
        gauge, the classic explained-variance fraction of the trace.
        """
        total = np.sum(np.maximum(self._eigvals, 0.0))
        captured = np.sum(np.maximum(self._eigvals[: self.n_components], 0.0))
        return float(captured / total) if total > 0 else 1.0

    # ------------------------------------------------------------------
    # Internal helpers
    # ------------------------------------------------------------------

    def _reduced_jacobian_full(
        self, energy: np.ndarray, **kwargs
    ) -> tuple[np.ndarray, np.ndarray]:
        """Compute stacked reduced relative Jacobian M and central flux."""
        kw = {**self._kwargs, **kwargs}
        jac_stack, central_flux = _build_jacobian_at_energy(self.model, energy, **kw)
        nonzero = central_flux != 0.0
        rel_jac = np.zeros_like(jac_stack)
        rel_jac[nonzero] = jac_stack[nonzero] / central_flux[nonzero, np.newaxis]
        M = rel_jac @ self.L_param
        return M, central_flux

    def _interpolate_B(self, energy: np.ndarray) -> np.ndarray:
        """Interpolate the residual blocks to arbitrary energies.

        Linear interpolation in log-energy; a convex combination of PSD
        blocks is PSD.  Outside the reference grid the edge block is held
        constant.
        """
        log_e = np.log(energy)
        idx = np.clip(
            np.searchsorted(self._ref_log_energy, log_e) - 1, 0, self._n_ref - 2
        )
        lo = self._ref_log_energy[idx]
        hi = self._ref_log_energy[idx + 1]
        t = np.clip((log_e - lo) / (hi - lo), 0.0, 1.0)
        return (1.0 - t)[:, np.newaxis, np.newaxis] * self.B[idx] + t[
            :, np.newaxis, np.newaxis
        ] * self.B[idx + 1]

    def _group_row_slice(self, ig: int, n_e: int) -> slice:
        """Return the row slice for group index ig in the stacked layout."""
        rpg = 2 * n_e if self._is_nucleon_model else n_e
        return slice(ig * rpg, (ig + 1) * rpg)

    # ------------------------------------------------------------------
    # Public interface — mirrors the original GSF model API
    # ------------------------------------------------------------------

    def flux(
        self,
        energy: ArrayLike,
        target: str | int | list[int],
        **kwargs,
    ) -> np.ndarray:
        """Calculate flux (delegated to the underlying model).

        Identical to ``model.flux(energy, target, ...)``.
        """
        kw = {**self._kwargs, **kwargs}
        return self.model.flux(energy, target, **kw)

    def jacobian(
        self,
        energy: ArrayLike,
        target: str | int | list[int],
        **kwargs,
    ) -> np.ndarray:
        """Calculate flux Jacobian (delegated to the underlying model).

        Identical to ``model.jacobian(energy, target, ...)``.
        """
        kw = {**self._kwargs, **kwargs}
        return self.model.jacobian(energy, target, **kw)

    def _group_pair_absolute_covariance(
        self,
        ig1: int,
        ig2: int,
        energy: np.ndarray,
        **kwargs,
    ) -> np.ndarray:
        """Compute full absolute covariance between two group indices.

        Returns the full abs_cov matrix including all nucleon blocks if
        applicable.  The low-rank part carries all cross-energy structure;
        the interpolated residual block is added on the same-energy
        diagonal for every sub-row pair of the two groups.
        """
        n_e = len(energy)
        nspg = self._n_sub_per_group

        M_full, cf = self._reduced_jacobian_full(energy, **kwargs)
        Bx = self._interpolate_B(energy)

        s1 = self._group_row_slice(ig1, n_e)
        s2 = self._group_row_slice(ig2, n_e)

        M1 = M_full[s1]
        M2 = M_full[s2]
        cf1 = cf[s1]
        cf2 = cf[s2]

        # Relative covariance for this group pair: low rank ...
        rel_cov = M1 @ M2.T
        # ... plus the same-energy residual block entries
        diag_idx = np.arange(n_e)
        for a in range(nspg):
            for b in range(nspg):
                rel_cov[a * n_e + diag_idx, b * n_e + diag_idx] += Bx[
                    :, ig1 * nspg + a, ig2 * nspg + b
                ]

        # Convert to absolute: Cov_abs[i,j] = rel_cov[i,j] * f1[i] * f2[j]
        return rel_cov * cf1[:, np.newaxis] * cf2[np.newaxis, :]

    def covariance(
        self,
        target1: str | int | list[int],
        target2: str | int | list[int],
        energy: ArrayLike,
        **kwargs,
    ) -> np.ndarray:
        """Approximate absolute flux covariance via low rank + energy blocks.

        Same signature and return units as the original model's ``covariance``.

        Parameters
        ----------
        target1, target2 : str, int, or list[int]
            Target specification (group name, element Z, or list of Z values).
        energy : array-like
            Energy values.
        **kwargs
            Override kwargs for model methods.

        Returns
        -------
        cov : ndarray, shape (N, N)
            Absolute flux covariance matrix.
        """
        energy = np.atleast_1d(np.asarray(energy, dtype=float))
        n_e = len(energy)
        ig1 = _group_index(target1, self.model)
        ig2 = _group_index(target2, self.model)

        abs_cov = self._group_pair_absolute_covariance(ig1, ig2, energy, **kwargs)

        if self._is_nucleon_model:
            # For nucleon models, the stacked layout is [p-rows, n-rows].
            # Sum the 4 blocks: Cov(p+n, p+n) = Cov(p,p) + Cov(p,n) + Cov(n,p) + Cov(n,n)
            cov_pp = abs_cov[:n_e, :n_e]
            cov_pn = abs_cov[:n_e, n_e:]
            cov_np = abs_cov[n_e:, :n_e]
            cov_nn = abs_cov[n_e:, n_e:]
            return cov_pp + cov_pn + cov_np + cov_nn
        return abs_cov

    def error(
        self,
        energy: ArrayLike,
        target: str | int | list[int],
        **kwargs,
    ) -> np.ndarray:
        """Approximate absolute 1-sigma flux uncertainty via the reduction.

        Same signature and return units as the original model's ``error``.

        Parameters
        ----------
        energy : array-like
            Energy values.
        target : str, int, or list[int]
            Target specification.
        **kwargs
            Override kwargs for model methods.

        Returns
        -------
        sigma : ndarray, shape (N,)
            Absolute 1-sigma flux uncertainties.
        """
        cov = self.covariance(target, target, energy, **kwargs)
        return np.sqrt(np.diag(cov))

    def total_flux(
        self,
        energy: ArrayLike,
        **kwargs,
    ) -> np.ndarray:
        """Total flux summed over all groups (delegated to model)."""
        kw = {**self._kwargs, **kwargs}
        return self.model.total_flux(energy, **kw)

    def total_error(
        self,
        energy: ArrayLike,
        **kwargs,
    ) -> np.ndarray:
        """Total flux uncertainty summed over all group pairs.

        Same as the original model's ``total_error`` but using the reduced
        representation.  Exact on the reference grid (the residual blocks
        carry all same-energy cross-group covariance).
        """
        energy = np.atleast_1d(np.asarray(energy, dtype=float))
        n_e = len(energy)
        total_cov = np.zeros((n_e, n_e))
        for g1 in _GROUPS:
            for g2 in _GROUPS:
                total_cov += self.covariance(g1, g2, energy, **kwargs)
        return np.sqrt(np.diag(total_cov))

    # ------------------------------------------------------------------
    # Nucleon-model-specific methods
    # ------------------------------------------------------------------

    def p_and_n_flux(
        self,
        energy: ArrayLike,
        target: str | int | list[int],
        **kwargs,
    ) -> np.ndarray:
        """Separate proton and neutron flux (delegated to model).

        Only available for nucleon-model backends.
        """
        kw = {**self._kwargs, **kwargs}
        return self.model.p_and_n_flux(energy, target, **kw)

    def p_and_n_covariance(
        self,
        target1: str | int | list[int],
        target2: str | int | list[int],
        energy: ArrayLike,
        **kwargs,
    ) -> tuple[np.ndarray, np.ndarray]:
        """Separate proton and neutron covariance matrices via the reduction.

        Returns
        -------
        cov_pp, cov_nn : ndarray, shape (N, N) each
            Proton-proton and neutron-neutron absolute covariance matrices.
        """
        energy = np.atleast_1d(np.asarray(energy, dtype=float))
        n_e = len(energy)
        ig1 = _group_index(target1, self.model)
        ig2 = _group_index(target2, self.model)

        abs_cov = self._group_pair_absolute_covariance(ig1, ig2, energy, **kwargs)

        # Extract p-p and n-n blocks
        cov_pp = abs_cov[:n_e, :n_e]
        cov_nn = abs_cov[n_e:, n_e:]
        return cov_pp, cov_nn

    def p_and_n_error(
        self,
        energy: ArrayLike,
        target: str | int | list[int],
        **kwargs,
    ) -> np.ndarray:
        """Separate proton and neutron uncertainties via the reduction.

        Returns
        -------
        errors : ndarray, shape (2, N)
            [0, :] proton flux uncertainties, [1, :] neutron flux uncertainties.
        """
        cov_pp, cov_nn = self.p_and_n_covariance(target, target, energy, **kwargs)
        return np.array([np.sqrt(np.diag(cov_pp)), np.sqrt(np.diag(cov_nn))])

    # ------------------------------------------------------------------
    # PCA-specific extras
    # ------------------------------------------------------------------

    def reduced_jacobian(
        self,
        energy: ArrayLike,
        target: str | int | list[int] | None = None,
        **kwargs,
    ) -> np.ndarray:
        """Reduced relative Jacobian M(E) = J_rel(E) @ L_param.

        Exact at any energy — no interpolation needed.

        Parameters
        ----------
        energy : array-like
            Energy values.
        target : str, int, list[int], or None
            If given, return only rows for this group.
            If None, return the full stacked matrix over all groups.
        **kwargs
            Override kwargs for model methods.

        Returns
        -------
        M : ndarray
            Reduced relative Jacobian. Shape depends on target:
            - target=None: (n_groups * rows_per_group, n_components)
            - target given: (rows_per_group, n_components) for that group
        """
        energy = np.atleast_1d(np.asarray(energy, dtype=float))
        M, _ = self._reduced_jacobian_full(energy, **kwargs)
        if target is not None:
            ig = _group_index(target, self.model)
            s = self._group_row_slice(ig, len(energy))
            return M[s]
        return M

    def sample(
        self,
        n_samples: int,
        energy: ArrayLike,
        rng: np.random.Generator | None = None,
        residual_noise: bool = True,
        **kwargs,
    ) -> np.ndarray:
        """Draw random flux realizations from the reduced model.

        By default, samples are drawn from the same low-rank + per-energy
        block covariance used by :meth:`covariance` and :meth:`error`:

            f_i = f_central,i * (1 + M_i @ phi + eps_i)

        Here ``phi ~ N(0, I_k)`` is shared by all energies and
        ``eps_i ~ N(0, B_i)`` is independent between requested energies.
        The residual term makes the sample covariance match the exact
        cross-group covariance of the approximation at every energy.

        Set ``residual_noise=False`` to draw only the correlated low-rank
        modes.  Those draws are smooth, but they do not sample the covariance
        returned by :meth:`covariance` and generally underestimate the
        pointwise variance.

        Parameters
        ----------
        n_samples : int
            Number of realizations to draw.
        energy : array-like
            Energy values.
        rng : numpy.random.Generator, optional
            Random number generator. Default: ``np.random.default_rng()``.
        residual_noise : bool, optional
            If True, include per-energy noise from the residual blocks B so
            that samples match the same-energy covariance returned by the
            uncertainty methods. Default True. Set to False for smooth
            low-rank-only draws.
        **kwargs
            Override kwargs for model methods.

        Returns
        -------
        samples : ndarray, shape (n_samples, rows)
            Absolute flux realizations (stacked over all groups).
        """
        if rng is None:
            rng = np.random.default_rng()

        energy = np.atleast_1d(np.asarray(energy, dtype=float))
        n_e = len(energy)

        M, central_flux = self._reduced_jacobian_full(energy, **kwargs)

        n_rows = M.shape[0]
        phi = rng.standard_normal((n_samples, self.n_components))
        rel_variation = phi @ M.T

        if residual_noise:
            Bx = self._interpolate_B(energy)  # (n_e, nsub, nsub)
            # matrix square root per energy (PSD by construction)
            w, u = np.linalg.eigh(Bx)
            A = u * np.sqrt(np.maximum(w, 0.0))[:, np.newaxis, :]
            z = rng.standard_normal((n_samples, n_e, self._n_sub))
            # eps[s, i, a] = sum_b A[i, a, b] * z[s, i, b]
            eps_sub = np.einsum("iab,sib->sia", A, z)
            eps = np.zeros((n_samples, n_rows))
            offsets = self._sub_row_offsets(n_e)
            for s_idx, off in enumerate(offsets):
                eps[:, off : off + n_e] = eps_sub[:, :, s_idx]
            rel_variation += eps

        return central_flux[np.newaxis, :] * (1.0 + rel_variation)

D: np.ndarray property

Diagonal of the block correction in stacked-row layout.

Equivalent to the scalar diagonal correction of the original construction (relative variance not carried by the low-rank part).

D_pn: np.ndarray | None property

p-n cross-term correction per group (nucleon models), from the block.

variance_explained: float property

Fraction of gauge-metric variance captured by the low-rank part.

In the correlation gauge this is the average fraction of the per-row variance resolved by the components; in the covariance gauge, the classic explained-variance fraction of the trace.

flux(energy: ArrayLike, target: str | int | list[int], **kwargs) -> np.ndarray

Calculate flux (delegated to the underlying model).

Identical to model.flux(energy, target, ...).

Source code in src/globalsplinefit/pca.py

def flux(
    self,
    energy: ArrayLike,
    target: str | int | list[int],
    **kwargs,
) -> np.ndarray:
    """Calculate flux (delegated to the underlying model).

    Identical to ``model.flux(energy, target, ...)``.
    """
    kw = {**self._kwargs, **kwargs}
    return self.model.flux(energy, target, **kw)

jacobian(energy: ArrayLike, target: str | int | list[int], **kwargs) -> np.ndarray

Calculate flux Jacobian (delegated to the underlying model).

Identical to model.jacobian(energy, target, ...).

Source code in src/globalsplinefit/pca.py

def jacobian(
    self,
    energy: ArrayLike,
    target: str | int | list[int],
    **kwargs,
) -> np.ndarray:
    """Calculate flux Jacobian (delegated to the underlying model).

    Identical to ``model.jacobian(energy, target, ...)``.
    """
    kw = {**self._kwargs, **kwargs}
    return self.model.jacobian(energy, target, **kw)

covariance(target1: str | int | list[int], target2: str | int | list[int], energy: ArrayLike, **kwargs) -> np.ndarray

Approximate absolute flux covariance via low rank + energy blocks.

Same signature and return units as the original model's covariance.

Parameters:

Name	Type	Description	Default
target1	str, int, or list[int]	Target specification (group name, element Z, or list of Z values).	required
target2	str, int, or list[int]	Target specification (group name, element Z, or list of Z values).	required
energy	array - like	Energy values.	required
**kwargs		Override kwargs for model methods.	{}

Returns:

Name	Type	Description
cov	(ndarray, shape(N, N))	Absolute flux covariance matrix.

Source code in src/globalsplinefit/pca.py

def covariance(
    self,
    target1: str | int | list[int],
    target2: str | int | list[int],
    energy: ArrayLike,
    **kwargs,
) -> np.ndarray:
    """Approximate absolute flux covariance via low rank + energy blocks.

    Same signature and return units as the original model's ``covariance``.

    Parameters
    ----------
    target1, target2 : str, int, or list[int]
        Target specification (group name, element Z, or list of Z values).
    energy : array-like
        Energy values.
    **kwargs
        Override kwargs for model methods.

    Returns
    -------
    cov : ndarray, shape (N, N)
        Absolute flux covariance matrix.
    """
    energy = np.atleast_1d(np.asarray(energy, dtype=float))
    n_e = len(energy)
    ig1 = _group_index(target1, self.model)
    ig2 = _group_index(target2, self.model)

    abs_cov = self._group_pair_absolute_covariance(ig1, ig2, energy, **kwargs)

    if self._is_nucleon_model:
        # For nucleon models, the stacked layout is [p-rows, n-rows].
        # Sum the 4 blocks: Cov(p+n, p+n) = Cov(p,p) + Cov(p,n) + Cov(n,p) + Cov(n,n)
        cov_pp = abs_cov[:n_e, :n_e]
        cov_pn = abs_cov[:n_e, n_e:]
        cov_np = abs_cov[n_e:, :n_e]
        cov_nn = abs_cov[n_e:, n_e:]
        return cov_pp + cov_pn + cov_np + cov_nn
    return abs_cov

error(energy: ArrayLike, target: str | int | list[int], **kwargs) -> np.ndarray

Approximate absolute 1-sigma flux uncertainty via the reduction.

Same signature and return units as the original model's error.

Parameters:

Name	Type	Description	Default
energy	array - like	Energy values.	required
target	str, int, or list[int]	Target specification.	required
**kwargs		Override kwargs for model methods.	{}

Returns:

Name	Type	Description
sigma	(ndarray, shape(N))	Absolute 1-sigma flux uncertainties.

Source code in src/globalsplinefit/pca.py

def error(
    self,
    energy: ArrayLike,
    target: str | int | list[int],
    **kwargs,
) -> np.ndarray:
    """Approximate absolute 1-sigma flux uncertainty via the reduction.

    Same signature and return units as the original model's ``error``.

    Parameters
    ----------
    energy : array-like
        Energy values.
    target : str, int, or list[int]
        Target specification.
    **kwargs
        Override kwargs for model methods.

    Returns
    -------
    sigma : ndarray, shape (N,)
        Absolute 1-sigma flux uncertainties.
    """
    cov = self.covariance(target, target, energy, **kwargs)
    return np.sqrt(np.diag(cov))

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

Total flux summed over all groups (delegated to model).

Source code in src/globalsplinefit/pca.py

def total_flux(
    self,
    energy: ArrayLike,
    **kwargs,
) -> np.ndarray:
    """Total flux summed over all groups (delegated to model)."""
    kw = {**self._kwargs, **kwargs}
    return self.model.total_flux(energy, **kw)

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

Total flux uncertainty summed over all group pairs.

Same as the original model's total_error but using the reduced representation. Exact on the reference grid (the residual blocks carry all same-energy cross-group covariance).

Source code in src/globalsplinefit/pca.py

def total_error(
    self,
    energy: ArrayLike,
    **kwargs,
) -> np.ndarray:
    """Total flux uncertainty summed over all group pairs.

    Same as the original model's ``total_error`` but using the reduced
    representation.  Exact on the reference grid (the residual blocks
    carry all same-energy cross-group covariance).
    """
    energy = np.atleast_1d(np.asarray(energy, dtype=float))
    n_e = len(energy)
    total_cov = np.zeros((n_e, n_e))
    for g1 in _GROUPS:
        for g2 in _GROUPS:
            total_cov += self.covariance(g1, g2, energy, **kwargs)
    return np.sqrt(np.diag(total_cov))

p_and_n_flux(energy: ArrayLike, target: str | int | list[int], **kwargs) -> np.ndarray

Separate proton and neutron flux (delegated to model).

Only available for nucleon-model backends.

Source code in src/globalsplinefit/pca.py

def p_and_n_flux(
    self,
    energy: ArrayLike,
    target: str | int | list[int],
    **kwargs,
) -> np.ndarray:
    """Separate proton and neutron flux (delegated to model).

    Only available for nucleon-model backends.
    """
    kw = {**self._kwargs, **kwargs}
    return self.model.p_and_n_flux(energy, target, **kw)

p_and_n_covariance(target1: str | int | list[int], target2: str | int | list[int], energy: ArrayLike, **kwargs) -> tuple[np.ndarray, np.ndarray]

Separate proton and neutron covariance matrices via the reduction.

Returns:

Type	Description
cov_pp, cov_nn : ndarray, shape (N, N) each	Proton-proton and neutron-neutron absolute covariance matrices.

Source code in src/globalsplinefit/pca.py

def p_and_n_covariance(
    self,
    target1: str | int | list[int],
    target2: str | int | list[int],
    energy: ArrayLike,
    **kwargs,
) -> tuple[np.ndarray, np.ndarray]:
    """Separate proton and neutron covariance matrices via the reduction.

    Returns
    -------
    cov_pp, cov_nn : ndarray, shape (N, N) each
        Proton-proton and neutron-neutron absolute covariance matrices.
    """
    energy = np.atleast_1d(np.asarray(energy, dtype=float))
    n_e = len(energy)
    ig1 = _group_index(target1, self.model)
    ig2 = _group_index(target2, self.model)

    abs_cov = self._group_pair_absolute_covariance(ig1, ig2, energy, **kwargs)

    # Extract p-p and n-n blocks
    cov_pp = abs_cov[:n_e, :n_e]
    cov_nn = abs_cov[n_e:, n_e:]
    return cov_pp, cov_nn

p_and_n_error(energy: ArrayLike, target: str | int | list[int], **kwargs) -> np.ndarray

Separate proton and neutron uncertainties via the reduction.

Returns:

Name	Type	Description
errors	(ndarray, shape(2, N))	[0, :] proton flux uncertainties, [1, :] neutron flux uncertainties.

Source code in src/globalsplinefit/pca.py

def p_and_n_error(
    self,
    energy: ArrayLike,
    target: str | int | list[int],
    **kwargs,
) -> np.ndarray:
    """Separate proton and neutron uncertainties via the reduction.

    Returns
    -------
    errors : ndarray, shape (2, N)
        [0, :] proton flux uncertainties, [1, :] neutron flux uncertainties.
    """
    cov_pp, cov_nn = self.p_and_n_covariance(target, target, energy, **kwargs)
    return np.array([np.sqrt(np.diag(cov_pp)), np.sqrt(np.diag(cov_nn))])

reduced_jacobian(energy: ArrayLike, target: str | int | list[int] | None = None, **kwargs) -> np.ndarray

Reduced relative Jacobian M(E) = J_rel(E) @ L_param.

Exact at any energy — no interpolation needed.

Parameters:

Name	Type	Description	Default
energy	array - like	Energy values.	required
target	str, int, list[int], or None	If given, return only rows for this group. If None, return the full stacked matrix over all groups.	None
**kwargs		Override kwargs for model methods.	{}

Returns:

Name	Type	Description
M	ndarray	Reduced relative Jacobian. Shape depends on target: - target=None: (n_groups * rows_per_group, n_components) - target given: (rows_per_group, n_components) for that group

Source code in src/globalsplinefit/pca.py

def reduced_jacobian(
    self,
    energy: ArrayLike,
    target: str | int | list[int] | None = None,
    **kwargs,
) -> np.ndarray:
    """Reduced relative Jacobian M(E) = J_rel(E) @ L_param.

    Exact at any energy — no interpolation needed.

    Parameters
    ----------
    energy : array-like
        Energy values.
    target : str, int, list[int], or None
        If given, return only rows for this group.
        If None, return the full stacked matrix over all groups.
    **kwargs
        Override kwargs for model methods.

    Returns
    -------
    M : ndarray
        Reduced relative Jacobian. Shape depends on target:
        - target=None: (n_groups * rows_per_group, n_components)
        - target given: (rows_per_group, n_components) for that group
    """
    energy = np.atleast_1d(np.asarray(energy, dtype=float))
    M, _ = self._reduced_jacobian_full(energy, **kwargs)
    if target is not None:
        ig = _group_index(target, self.model)
        s = self._group_row_slice(ig, len(energy))
        return M[s]
    return M

sample(n_samples: int, energy: ArrayLike, rng: np.random.Generator | None = None, residual_noise: bool = True, **kwargs) -> np.ndarray

Draw random flux realizations from the reduced model.

By default, samples are drawn from the same low-rank + per-energy block covariance used by :meth:covariance and :meth:error:

f_i = f_central,i * (1 + M_i @ phi + eps_i)

Here phi ~ N(0, I_k) is shared by all energies and eps_i ~ N(0, B_i) is independent between requested energies. The residual term makes the sample covariance match the exact cross-group covariance of the approximation at every energy.

Set residual_noise=False to draw only the correlated low-rank modes. Those draws are smooth, but they do not sample the covariance returned by :meth:covariance and generally underestimate the pointwise variance.

Parameters:

Name	Type	Description	Default
n_samples	int	Number of realizations to draw.	required
energy	array - like	Energy values.	required
rng	Generator	Random number generator. Default: np.random.default_rng().	None
residual_noise	bool	If True, include per-energy noise from the residual blocks B so that samples match the same-energy covariance returned by the uncertainty methods. Default True. Set to False for smooth low-rank-only draws.	True
**kwargs		Override kwargs for model methods.	{}

Returns:

Name	Type	Description
samples	(ndarray, shape(n_samples, rows))	Absolute flux realizations (stacked over all groups).

Source code in src/globalsplinefit/pca.py

def sample(
    self,
    n_samples: int,
    energy: ArrayLike,
    rng: np.random.Generator | None = None,
    residual_noise: bool = True,
    **kwargs,
) -> np.ndarray:
    """Draw random flux realizations from the reduced model.

    By default, samples are drawn from the same low-rank + per-energy
    block covariance used by :meth:`covariance` and :meth:`error`:

        f_i = f_central,i * (1 + M_i @ phi + eps_i)

    Here ``phi ~ N(0, I_k)`` is shared by all energies and
    ``eps_i ~ N(0, B_i)`` is independent between requested energies.
    The residual term makes the sample covariance match the exact
    cross-group covariance of the approximation at every energy.

    Set ``residual_noise=False`` to draw only the correlated low-rank
    modes.  Those draws are smooth, but they do not sample the covariance
    returned by :meth:`covariance` and generally underestimate the
    pointwise variance.

    Parameters
    ----------
    n_samples : int
        Number of realizations to draw.
    energy : array-like
        Energy values.
    rng : numpy.random.Generator, optional
        Random number generator. Default: ``np.random.default_rng()``.
    residual_noise : bool, optional
        If True, include per-energy noise from the residual blocks B so
        that samples match the same-energy covariance returned by the
        uncertainty methods. Default True. Set to False for smooth
        low-rank-only draws.
    **kwargs
        Override kwargs for model methods.

    Returns
    -------
    samples : ndarray, shape (n_samples, rows)
        Absolute flux realizations (stacked over all groups).
    """
    if rng is None:
        rng = np.random.default_rng()

    energy = np.atleast_1d(np.asarray(energy, dtype=float))
    n_e = len(energy)

    M, central_flux = self._reduced_jacobian_full(energy, **kwargs)

    n_rows = M.shape[0]
    phi = rng.standard_normal((n_samples, self.n_components))
    rel_variation = phi @ M.T

    if residual_noise:
        Bx = self._interpolate_B(energy)  # (n_e, nsub, nsub)
        # matrix square root per energy (PSD by construction)
        w, u = np.linalg.eigh(Bx)
        A = u * np.sqrt(np.maximum(w, 0.0))[:, np.newaxis, :]
        z = rng.standard_normal((n_samples, n_e, self._n_sub))
        # eps[s, i, a] = sum_b A[i, a, b] * z[s, i, b]
        eps_sub = np.einsum("iab,sib->sia", A, z)
        eps = np.zeros((n_samples, n_rows))
        offsets = self._sub_row_offsets(n_e)
        for s_idx, off in enumerate(offsets):
            eps[:, off : off + n_e] = eps_sub[:, :, s_idx]
        rel_variation += eps

    return central_flux[np.newaxis, :] * (1.0 + rel_variation)

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() (the promoted 2026 set, Solar Cycle 24 average).	None
n_pivots	int	Number of log-spaced pivot energies per species. If neither this nor pivot_energies is given (and energy_range is left at its default), the published grid for the model version is used (:data:RECOMMENDED_PIVOTS; 12 pivots for "2026", worst-case coverage factor 1.24) — the reproducible, citable default. Model versions with no shipped entry get an equivalent grid derived at construction time by :func:optimize_pivots (deterministic, ~10 s, with a :exc:UserWarning showing the literal to pin for repeated runs) rather than a silent drop to a naive log-spaced grid. Passing n_pivots explicitly requests a log-spaced grid instead.	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).	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 ["p_p", "p_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()                    # 20 parameters (2 x 10)
>>> 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)             # -> 1.0 if uncorrelated

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()`` (the
        promoted 2026 set, Solar Cycle 24 average).
    n_pivots : int, optional
        Number of log-spaced pivot energies per species.  If neither this
        nor ``pivot_energies`` is given (and ``energy_range`` is left at
        its default), the **published grid** for the model version is used
        (:data:`RECOMMENDED_PIVOTS`; 12 pivots for "2026", worst-case
        coverage factor 1.24) — the reproducible, citable default.  Model
        versions with no shipped entry get an equivalent grid derived at
        construction time by :func:`optimize_pivots` (deterministic, ~10 s,
        with a :exc:`UserWarning` showing the literal to pin for repeated
        runs) rather than a silent drop to a naive log-spaced grid.
        Passing ``n_pivots`` explicitly requests a log-spaced grid instead.
    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).
    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 ``["p_p", "p_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()                    # 20 parameters (2 x 10)
    >>> 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)             # -> 1.0 if uncorrelated
    """

    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:
                if model.version in RECOMMENDED_PIVOTS:
                    # the published, citable grid for this model version
                    pivot_energies = np.array(RECOMMENDED_PIVOTS[model.version])
                else:
                    # no shipped table for this version — derive an equivalent
                    # grid now and tell the caller how to pin it
                    pivot_energies = _derive_pivots(model, per_group, basis, kwargs)
            else:
                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.atleast_1d(np.asarray(pivot_energies, dtype=float))
        if len(pivot_energies) < 2 or np.any(np.diff(pivot_energies) <= 0):
            raise ValueError("pivot_energies must be >= 2 strictly 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 = kwargs

        n_piv = len(pivot_energies)
        jac_rel, cov_par, central, self.species = _relative_species_system(
            model, pivot_energies, per_group, **kwargs
        )
        self.cov = jac_rel @ cov_par @ jac_rel.T
        self.cov = 0.5 * (self.cov + self.cov.T)
        self._central_pivots = central.reshape(len(self.species), n_piv)
        self._chol = None  # lazy Cholesky of cov for penalty()

        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
        ]

    # ------------------------------------------------------------------
    # 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
        return self.cov / np.outer(s, s)

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

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

        Cardinal in the pivot values (``H`` is the identity at the pivots)
        and a partition of unity, for either ``basis_type``:

        - ``"spline"`` (default): local cubic (Catmull-Rom) spline in
          log-energy — smooth (C1) deformations, each component confined to
          its two neighboring intervals.
        - ``"hat"``: piecewise-linear hat functions in log-energy — also
          local, non-negative, but kinked.

        Outside the pivot range the edge pivot's deformation is held
        constant.
        """
        log_e = np.log(np.atleast_1d(np.asarray(energy, dtype=float)))
        return _interp_basis(self._log_pivots, log_e, self.basis_type)

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

    def _central_flux(self, energy: np.ndarray, **kwargs) -> np.ndarray:
        """Central flux per species, shape (S, n_E)."""
        kw = {**self._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 = np.atleast_1d(np.asarray(energy, dtype=float))
        central = self._central_flux(energy, **kwargs)
        if theta is None:
            return central
        theta = np.asarray(theta, dtype=float).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 = np.atleast_1d(np.asarray(energy, dtype=float))
        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; hat-interpolated in between.

        Returns
        -------
        sigma : ndarray, shape (S, n_E)
        """
        energy = np.atleast_1d(np.asarray(energy, dtype=float))
        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()
        w, u = np.linalg.eigh(self.cov)
        a = u * np.sqrt(np.maximum(w, 0.0))
        return rng.standard_normal((n_samples, self.n_params)) @ a.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.
        """
        if self._chol is None:
            self._chol = np.linalg.cholesky(self.cov)
        z = np.linalg.solve(self._chol, np.asarray(theta, dtype=float))
        return float(z @ z)

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

    def to_dict(self) -> dict:
        """JSON-serializable description of the reduction.

        Contains everything a downstream fitter needs besides the central
        flux itself: pivot energies, species labels, the central flux at
        the pivots, and the component covariance.
        """
        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).

Cardinal in the pivot values (H is the identity at the pivots) and a partition of unity, for either basis_type:

• "spline" (default): local cubic (Catmull-Rom) spline in log-energy — smooth (C1) deformations, each component confined to its two neighboring intervals.
• "hat": piecewise-linear hat functions in log-energy — also local, non-negative, but kinked.

Outside the pivot range the edge pivot's deformation is held constant.

Source code in src/globalsplinefit/reduced.py

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

    Cardinal in the pivot values (``H`` is the identity at the pivots)
    and a partition of unity, for either ``basis_type``:

    - ``"spline"`` (default): local cubic (Catmull-Rom) spline in
      log-energy — smooth (C1) deformations, each component confined to
      its two neighboring intervals.
    - ``"hat"``: piecewise-linear hat functions in log-energy — also
      local, non-negative, but kinked.

    Outside the pivot range the edge pivot's deformation is held
    constant.
    """
    log_e = np.log(np.atleast_1d(np.asarray(energy, dtype=float)))
    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 = np.atleast_1d(np.asarray(energy, dtype=float))
    central = self._central_flux(energy, **kwargs)
    if theta is None:
        return central
    theta = np.asarray(theta, dtype=float).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 = np.atleast_1d(np.asarray(energy, dtype=float))
    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; hat-interpolated in between.

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; hat-interpolated in between.

    Returns
    -------
    sigma : ndarray, shape (S, n_E)
    """
    energy = np.atleast_1d(np.asarray(energy, dtype=float))
    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()
    w, u = np.linalg.eigh(self.cov)
    a = u * np.sqrt(np.maximum(w, 0.0))
    return rng.standard_normal((n_samples, self.n_params)) @ a.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.

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.
    """
    if self._chol is None:
        self._chol = np.linalg.cholesky(self.cov)
    z = np.linalg.solve(self._chol, np.asarray(theta, dtype=float))
    return float(z @ z)

to_dict() -> dict

JSON-serializable description of the reduction.

Contains everything a downstream fitter needs besides the central flux itself: pivot energies, species labels, the central flux at the pivots, and the component covariance.

Source code in src/globalsplinefit/reduced.py

def to_dict(self) -> dict:
    """JSON-serializable description of the reduction.

    Contains everything a downstream fitter needs besides the central
    flux itself: pivot energies, species labels, the central flux at
    the pivots, and the component covariance.
    """
    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).

Runtime is dominated by the exchange loop (roughly half a minute for the defaults; scales with n_pivots * n_grid * n_restarts and is a few times slower with per_group=True). 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).

    Runtime is dominated by the exchange loop (roughly half a minute for
    the defaults; scales with ``n_pivots * n_grid * n_restarts`` and is a
    few times slower with ``per_group=True``).  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)

    # anti-aliasing guards: candidates on every second grid point (so the
    # objective always samples between pivots) and a minimum separation
    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

Dictionary mapping group names to atomic numbers.

Maps string group names ("p", "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