diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml new file mode 100644 index 0000000..f9828b0 --- /dev/null +++ b/.github/workflows/publish.yml @@ -0,0 +1,75 @@ +name: Publish Python 🐍 distribution 📦 to PyPI and TestPyPI + +on: + release: + types: [published] + +jobs: + build: + name: Build distribution 📦 + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.x" + + - name: Install pypa/build + run: >- + python3 -m + pip install + build + --user + - name: Build a binary wheel and a source tarball + run: python3 -m build + - name: Store the distribution packages + uses: actions/upload-artifact@v3 + with: + name: python-package-distributions + path: dist/ + + testpypi-publish: + name: Upload release to Test-PyPI + needs: + - build + runs-on: ubuntu-latest + + environment: + name: testpypi + url: https://test.pypi.org/p/fgspectra + + permissions: + id-token: write # IMPORTANT: this permission is mandatory for trusted publishing + steps: + - name: Download all the dists + uses: actions/download-artifact@v3 + with: + name: python-package-distributions + path: dist/ + - name: Publish package distributions to TestPyPI + uses: pypa/gh-action-pypi-publish@release/v1 + with: + repository-url: https://test.pypi.org/legacy/ + + pypi-publish: + name: Upload release to PyPI + needs: + - build + runs-on: ubuntu-latest + + environment: + name: publish_pypi + url: https://pypi.org/p/fgspectra + + permissions: + id-token: write # IMPORTANT: this permission is mandatory for trusted publishing + steps: + - name: Download all the dists + uses: actions/download-artifact@v3 + with: + name: python-package-distributions + path: dist/ + - name: Publish package distributions to PyPI + uses: pypa/gh-action-pypi-publish@release/v1 diff --git a/fgspectra/cross.py b/fgspectra/cross.py index 1c97bd3..777c78a 100644 --- a/fgspectra/cross.py +++ b/fgspectra/cross.py @@ -82,7 +82,7 @@ class FactorizedCrossSpectrum(Model): sed : callable :math:`f(\nu)`. It returns an array with shape ``(..., freq)``. It can be :class:`fgspectra.frequency.SED`. - cl_args : callable + cl : callable :math:`C_\ell`. It returns an array with shape ``(..., ell)``. It can be :class:`fgspectra.power.PowerSpectrum` @@ -242,6 +242,121 @@ def eval(self, sedT_kwargs={}, sedE_kwargs={}, cl_kwargs={}): return np.einsum("l...i,l...j,...l->...ijl", fT_nu, fE_nu, cl) +class DecorrelatedFactorizedCrossSpectrum(Model): + r""" Decorrelated factorized cross spectrum + + Cross-spectrum of **one** component for which the scaling in frequency + and in multipoles are factorizable. The scaling in frequency is not exact + and some decorrelation is specified. + + .. math:: xC_{\ell}^{(ij)} = f_{decor}(\nu_i, \nu_j) f(\nu_j) f(\nu_i) C_{\ell} + + This model implements a decorrelation by rigidely rescalling a decorelation ratio + with respect to a reference cross-term. For example, for 3 frequencies, with the + reference taken at $\nu_1\times\nu_3$: + + .. math:: f_{decor} = I(3) + \Delta_{decor} \begin{pmatrix} + 0 & r_{12} & 1 \\ + r_{12} & 0 & r_{23} \\ + 1 & r_{23} & 0 + \end{pmatrix} + + Parameters + ---------- + sed : callable + :math:`f(\nu)`. It returns an array with shape ``(..., freq)``. + It can be :class:`fgspectra.frequency.SED`. + cl : callable + :math:`C_\ell`. It returns an array with shape ``(..., ell)``. + It can be :class:`fgspectra.power.PowerSpectrum` + + Note + ---- + The two (optional) sets of extra dimensions ``...`` must be + broadcast-compatible. + """ + + def __init__(self, sed, cl, **kwargs): + self._sed = sed + self._cl = cl + self.set_defaults(**kwargs) + + def set_defaults(self, **kwargs): + if "sed_kwargs" in kwargs: + self._sed.set_defaults(**kwargs["sed_kwargs"]) + if "cl_kwargs" in kwargs: + self._cl.set_defaults(**kwargs["cl_kwargs"]) + kw_decor = { + key: kwargs[key] + for key in kwargs.keys() + if key not in ["sed_kwargs", "cl_kwargs"] + } + super().set_defaults(**kw_decor) + + @property + def defaults(self): + kw_decor = { + key: super().defaults[key] + for key in super().defaults.keys() + if key not in ["kwargs"] + } + return { + **kw_decor, + "sed_kwargs": self._sed.defaults, + "cl_kwargs": self._cl.defaults, + } + + def _get_repr(self): + decor_repr = super()._get_repr() + key = list(decor_repr.keys())[0] + decor_repr[key] = { + k: decor_repr[key][k] + for k in decor_repr[key].keys() + if k not in ["sed_kwargs", "cl_kwargs"] + } + decor_repr["Decorrelation"] = decor_repr.pop(key) + + sed_repr = self._sed._get_repr() + key = list(sed_repr.keys())[0] + sed_repr[key + " (SED)"] = sed_repr.pop(key) + + cl_repr = self._cl._get_repr() + key = list(cl_repr.keys())[0] + cl_repr[key + " (Cl)"] = cl_repr.pop(key) + + return {type(self).__name__: [decor_repr, sed_repr, cl_repr]} + + def eval(self, decor=None, f_decor=None, sed_kwargs={}, cl_kwargs={}): + """Compute the model at frequency and ell combinations. + + Parameters + ---------- + sed_args : list + Arguments for which the `sed` is evaluated. + cl_args : list + Arguments for which the `cl` is evaluated. + decor : float + Decorelation facot, by which f_decor is rescaled + f_decor : ndarray + Array of the decorelation scaling in frequency. + Shape is ``(freq, freq)``. + Returns + ------- + cross : ndarray + Cross-spectrum. The shape is ``(..., freq, freq, ell)``. + """ + + decorrelation = np.eye(f_decor.shape[0]) + decor * f_decor + + f_nu = self._sed(**sed_kwargs) + cl = self._cl(**cl_kwargs) + if f_nu.shape[0] != cl.shape[-1] or (f_nu.shape[0] == 1 and cl.shape[-1] == 1): + f_nu = f_nu[np.newaxis] + + res = np.einsum("ij,l...i,l...j,...l->...ijl", decorrelation, f_nu, f_nu, cl) + + return res + class CorrelatedFactorizedCrossSpectrum(FactorizedCrossSpectrum): r"""Factorized cross-spectrum of correlated components diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py index cf2a853..42eabe0 100644 --- a/fgspectra/frequency.py +++ b/fgspectra/frequency.py @@ -20,6 +20,25 @@ H_OVER_KT_CMB = constants.h * 1e9 / constants.k / T_CMB +def _flux2cmb(nu): + """Converts flux to thermodynamics units""" + x = H_OVER_KT_CMB * nu + g2_min1 = ( + 2.0 + * constants.k**3 + * T_CMB**2 + * x**4 + * np.exp(x) + / (constants.h * constants.c * np.expm1(x)) ** 2 + ) + return 1.0 / g2_min1 + + +def _rj2cmb(nu): + x = H_OVER_KT_CMB * nu + return (np.expm1(x) / x) ** 2 / np.exp(x) + + def _bandpass_integration(): """Bandpass integrated version of the caller @@ -100,11 +119,6 @@ def _bandpass_integration(): return res -def _rj2cmb(nu): - x = H_OVER_KT_CMB * nu - return (np.expm1(x) / x) ** 2 / np.exp(x) - - class PowerLaw(Model): r"""Power Law @@ -313,6 +327,37 @@ def eval(self, nu=None, amp=1.0): return amp * np.ones_like(np.array(nu)) +class FreeSED(Model): + """Frequency-dependent component for which every entries of the SED are specifified.""" + + def eval(self, nu=None, sed=None): + """Evaluation of the SED + + Parameters + ---------- + nu: float or array + It just determines the shape of the output. + sed: float or array + Values of the SED. Must be the same shape as nu. + There is no normalisation, entries are used as is. + + Returns + ------- + sed: ndarray + If `nu` is an array, the shape is ``amp.shape + (freq)``. + If `nu` is scalar, the shape is ``amp.shape + (1)``. + Note that the last dimension is guaranteed to be the frequency. + """ + if isinstance(nu, (list, np.ndarray)): + try: + assert len(nu) == len(sed[..., :]) + except: + print("SED and nu must have the same shape.") + if isinstance(nu, list): + return _bandpass_integration() + return np.asarray(sed) + + class Join(Model): """Join several SED models together""" diff --git a/fgspectra/power.py b/fgspectra/power.py index 69d1d3e..ce0303e 100644 --- a/fgspectra/power.py +++ b/fgspectra/power.py @@ -82,7 +82,12 @@ def __init__(self, filenames, **kwargs): def eval(self, ell=None, ell_0=None, amp=1.0): """Compute the power spectrum with the given ell and parameters.""" - return amp * self._cl[..., ell] / self._cl[..., ell_0, np.newaxis] + #### TEMPORARY #### + if ell_0 is None: + print("NO NORMALISATION") + return amp * self._cl[..., ell] + else: + return amp * self._cl[..., ell] / self._cl[..., ell_0, np.newaxis] class PowerLawExtendedTemplate(PowerSpectrumFromFile): diff --git a/pyproject.toml b/pyproject.toml index 5b7fe82..917c78d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,11 +4,12 @@ build-backend = "setuptools.build_meta" [project] name = "fgspectra" -version = "1.1.0" +version = "1.2.1" description = "Foreground SED and power spectrum library for building cross frequency power spectrum models" readme = "README.md" authors = [ { name="Zack Li", email="xzackli@gmail.com" }, + { name="Benjamin Beringue", email="beringueb@gmail.com"} ] classifiers = [ "License :: OSI Approved :: MIT License",