Merge branch 'main' into dev_beam

.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

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
 

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"""
 

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

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",