Merge branch 'dev' of https://github.com/bastiencarreres/snsim into dev

.github/workflows/python-package.yml

@@ -5,9 +5,9 @@ name: Python package
 
 on:
   push:
-    branches: [ "dev_2" ]
+    branches: [ "main" ]
   pull_request:
-    branches: [ "dev_2" ]
+    branches: [ "main" ]
 
 jobs:
   build:

Examples/SN_simulation.ipynb

@@ -948,7 +948,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 17,
    "id": "2d78191a",
    "metadata": {},
    "outputs": [
@@ -1290,7 +1290,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.11.9"
   }
  },
  "nbformat": 4,

README.md

@@ -4,7 +4,7 @@
 <img src="docs/_static/snsimlogo.svg" width=350>
 
 # Code for simulation of transient surveys
-[![Documentation Status](https://readthedocs.org/projects/snsim/badge/?version=dev)](https://snsim.readthedocs.io/en/main/?badge=dev) ![Tests](https://github.com/bastiencarreres/snsim/actions/workflows/python-package.yml/badge.svg?branch=dev_2)
+[![Documentation Status](https://readthedocs.org/projects/snsim/badge/?version=dev)](https://snsim.readthedocs.io/en/main/?badge=dev) ![Tests](https://github.com/bastiencarreres/snsim/actions/workflows/python-package.yml/badge.svg?branch=dev)
 ## Installation
 In the setup.py directory use:
 ```

docs/configfile.rst

@@ -26,8 +26,8 @@ The input file is a .yml with the following structure:
       end_day: MJD NUMBER or 'YYYY-MM-DD'  # Opt, default given by survey file
       duration: SURVEY DURATION (DAYS)  # Opt, default given by survey file
       field_map: FIELD MAP FILE  # Opt, default is rectangle field
-      fake_skynoise: [VALUE, 'add' or 'replace'] # Opt, default is no fake_skynoise
       sub_field: 'sub_field_key' # Used to divided observation in CCD quadrant for example
+      cdd_noise : sig_ccd_noise # Optional, default is 0 ADU
    sim_par:
       z_range: [ZMIN, ZMAX] REDSHIFT RANGE
       randseed: RANDSEED TO REPRODUCE SIMULATION  # Optional
@@ -114,7 +114,7 @@ This section contains informations about the survey configuration :
    first day of the SQL database.
 -  **end_day** same as **start_day** but for the end of the survey.
    *type* : float or str. *Optional* : default is the last day of the
-   SQL database.
+   observations.
 -  **duration** : instead of setting an **end_day** you can specify a
    duration in **days**. *type* : float. *Optional* : the **duration**
    is ignored if an **end_day** is configured.
@@ -135,15 +135,7 @@ This section contains informations about the survey configuration :
    *Optional*
 -  **add_data** is a list of database key that you want to retrieve in
    lightcurves tables. *type* : list(str). *Optional*
--  **fake_skynoise** allow to add or replace the skynoise term. The fake
-   skynoise is multiply by the **PSF** if there is one given. This is a
-   list : [VALUE, ‘add’ or ‘replace’] the VALUE is the skynoise value in
-   ADU, if you use ‘add’ the fake_skynoise is added to skynoise from the
-   SQL database, else, if you use ‘replace’ the skynoise from SQL
-   database is just ignored. Note that if you set **fake_skynoise** with
-   ‘replace’ option and **sig_psf** = 0, the skynoise is exactly the
-   **fake_skynoise** value. *type* : list(float, str). *Optional*
-   default is no **fake_skynoise**
+-  **ccd_noise** is the noise from instrument in ADU / pixels. *Optional* : default is 0
 
 sim_par
 -------

docs/obsfile.rst

@@ -1,7 +1,7 @@
 Observation database file
 =========================
 
-It’s a **SQL DataBase** (.db) or a **Comma-separated values** (.csv)
+It’s a **Comma-separated values** (.csv) or **parquet** (.parquet)
 file which contain observations informations. It’s used to find obs
 epoch and their noise.
 
@@ -43,11 +43,11 @@ setting the two additional columns:
 +-----------------------------------+-----------------------------------+
 
 In addition you can take into account the variation of the PSF as the
-**F**\ ull **W**\ idth at **H**\ alf **M**\ aximum
-:math:`FWHM = 2 \sqrt{2 \log(2)} \sigma_{skynoise}`
+**F**\ ull **W**\ idth at **H**\ alf **M**\ aximum 
+:math:`FWHM = 2 \sqrt{2 \log(2)} \sigma_\mathrm{psf}`
 
 +-------------------------------------------+
-| FWHMeff                                   |
+| fwhm_psf                                  |
 +===========================================+
 | The Full Width at Half Maximum of the PSF |
 +-------------------------------------------+
@@ -56,7 +56,7 @@ And you can set a different gain for each observation by giving the
 **gain** column :
 
 +------------------------+
-| **gain**               |
+| gain                   |
 +========================+
 | The CCD gain in e-/ADU |
 +------------------------+
@@ -79,19 +79,28 @@ a 4 x 4 grid, you have to put something like that in your .dat file :
 If a sub field is not observed you should set the ID value to -1.
 
 In addition, you can add space between subfield by adding a header
-(begin line with %) that defines some “space-symbols”: \```pseudocode %
-#:ra:0.13 % @:dec:0.13
+(begin line with %) that defines some “space-symbols”: 
 
-ID01:ID02:#:ID03:ID04 ID05:ID06:#:ID07:ID08 @ ID09:ID10:#:ID11:ID12
-ID13:ID14:#:ID15:ID16 \``\` In the previous example the symbol # is used
+.. code:: pseudocode 
+
+   % #:ra:0.13 
+   % @:dec:0.13
+
+   ID01:ID02:#:ID03:ID04 
+   ID05:ID06:#:ID07:ID08 
+   @ 
+   ID09:ID10:#:ID11:ID12
+   ID13:ID14:#:ID15:ID16 
+
+In the previous example the symbol # is used
 has a ra space of 0.13 degrees and the @ is used has a dec space of 0.13
 degrees.
 
 You can show the sub filed map by :
 
 .. code:: python
 
-   sim.survey.fields.show_map()
+   sim.survey.show_map()
 
 .. figure:: _static/show_map.png
    :alt: show_map

docs/simmath.rst

@@ -53,24 +53,17 @@ The flux error is computed as :
 .. math::
 
 
-   \sigma^2_F = \frac{F}{G} + \sigma_{skynoise}^2 + \left(\frac{\ln(10)}{2.5}F\sigma_{zp}\right)^2
+   \sigma^2_F = \frac{F}{G} + \sigma_{skynoise}^2 + \left(\sigma_\mathrm{CCD}^2\right) \times 4\pi\sigma_{PSF}^2+ \left(\frac{\ln(10)}{2.5}F\sigma_{zp}\right)^2
 
 The first term is the Poisson noise with **G** the gain in $ e^- $ /
 ADU.
 
-The second term is the noise from sky flux. If there is a PSF given this
-term take into account the PSF by applying :
-
-.. math::
-
-
-   \sigma_{skynoise}^2  *= 4\pi\sigma_{PSF}^2
+The second term is the noise from sky flux.
 
 If you use limiting magnitude at 5σ, sky-noise is computed as :
 
 .. math::
 
-
    \sigma_{skynoise} = \frac{1}{5}10^{0.4\left(ZP - m_{5\sigma}\right)}
 
 The last term is the propagation of the zero point incertitude.

snsim/__init__.py

@@ -13,7 +13,7 @@
 
 __snsim_dir_path__ = os.path.dirname(__file__)
 
-__version__ = "0.5.0+dev"
+__version__ = "0.5.0"
 
 from .simu import Simulator
 from .sample import SimSample

snsim/astrobj.py

@@ -108,31 +108,23 @@ def _init_model(self, effects):
             SN simulation model.
 
         """
-
-        if "model_version" not in self._sim_par:
-            version = None
-        else:
-            version = self._sim_par["model_version"]
-
-        snc_source = snc.get_source(
-            name=self._sim_par["model_name"], version=version
-        )
+        snc_source = self.source
 
         if "model_version" not in self._sim_par:
             self._sim_par["model_version"] = snc_source.version
 
         if effects is not None:
-            eff = [eff["source"] for eff in effects]
+            eff_sources = [eff["source"] for eff in effects]
             eff_names = [eff["name"] for eff in effects]
             eff_frames = [eff["frame"] for eff in effects]
         else:
-            eff = None
+            eff_sources = None
             eff_names = None
             eff_frames = None
 
         model = snc.Model(
             source=snc_source,
-            effects=eff,
+            effects=eff_sources,
             effect_names=eff_names,
             effect_frames=eff_frames,
         )
@@ -266,17 +258,34 @@ def gen_flux(self, obs, mod_fcov=False, seed=None):
             if k not in sim_lc.columns:
                 sim_lc[k] = obs[k].values
 
+        snc_par = {k: v for k, v in zip(self.sim_model.param_names, self.sim_model.parameters) if k!= 'z'}
         sim_lc.attrs = {
             "mu": self.mu,
             "zobs": self.zobs,
             "zCMB": self.zCMB,
-            **self._sim_par,
+            "effects": self.sim_model.effect_names,
+            **snc_par,
+            **self._sim_par
         }
 
         sim_lc.reset_index(inplace=True, drop=True)
         sim_lc.index.set_names("epochs", inplace=True)
         return sim_lc
 
+    def mag_restframeband_to_amp(self, mag, band, magsys, amp_param_name='x0'):
+        source = self.source
+        m_current = source.peakmag(band, magsys)
+        return 10.**(0.4 * (m_current - mag)) * source.get(amp_param_name)
+
+    @property
+    def source(self):
+        if "model_version" not in self._sim_par:
+            version = None
+        else:
+            version = self._sim_par["model_version"]
+
+        return snc.get_source(name=self._sim_par["model_name"], version=version)
+
     @property
     def zpec(self):
         """Get peculiar velocity redshift."""
@@ -383,12 +392,11 @@ def _set_model_par(self, model):
 
         self._sim_par["mb"] = mb
 
-        # Set x1 and c
-        model.set(x1=self._sim_par["x1"], c=self._sim_par["c"])
-
         # Compute the x0 parameter
-        model.set_source_peakmag(self._sim_par["mb"], "bessellb", "ab")
-        self._sim_par["x0"] = model.get("x0")
+        self._sim_par["x0"] = self.mag_restframeband_to_amp(self._sim_par["mb"], 'bessellb', 'ab')
+
+        # Set x1 and c
+        model.set(x0=self._sim_par["x0"], x1=self._sim_par["x1"], c=self._sim_par["c"])
         return model
 
     @staticmethod

snsim/constants.py

@@ -5,6 +5,8 @@
 from astropy import constants as cst
 import numpy as np
 import shapely.geometry as shp_geo
+from . import dust_utils as dst_ut
+from . import scatter as sct
 
 path_location = Path(__file__).absolute().parent
 init_location = path_location / "__init__.py"
@@ -34,4 +36,4 @@
 h_article = {"jla": 0.70, "li11": 0.73, "sullivan06": 0.70}
 
 # value of fitted parameter of SNIa-Host_galaxy from Sullivan et al 2006 https://iopscience.iop.org/article/10.1086/506137/pdf
-sullivan_para = {"mass": 5.3 * 1.0e-14, "SFR": 3.9 * 1.0e-4}
+sullivan_para = {"mass": 5.3 * 1.0e-14, "SFR": 3.9 * 1.0e-4}

snsim/generators.py

@@ -783,33 +783,21 @@ def _add_print(self):
 
     def _add_effects(self):
         effects = []
+        args = []
         # Add scattering model if needed
         if "sct_model" in self._params:
             if self._params["sct_model"] == "G10":
                 if len(self.sim_sources["model_name"]) > 1 or self.sim_sources[
                     "model_name"
                 ][0] not in ["salt2", "salt3"]:
                     raise ValueError("G10 cannot be used")
-                effects.append(
-                    {
-                        "source": sct.G10(
-                            snc.get_source(
+                args = [
+                    snc.get_source(
                                 name=self.sim_sources["model_name"][0],
                                 version=self.sim_sources["model_version"][0],
                             )
-                        ),
-                        "frame": "rest",
-                        "name": "G10_",
-                    }
-                )
-            elif self._params["sct_model"] in ["C11_0", "C11_1", "C11_2"]:
-                effects.append({"source": sct.C11(),
-                                "frame": "rest",
-                                "name": "C11_"})
-            elif self._params["sct_model"].lower() == "bs20":
-                effects.append({"source": snc.CCM89Dust(),
-                                "frame": "rest", 
-                                "name": "BS20_"})
+                        ]
+        effects.append(sct.init_sn_sct_model(self._params["sct_model"], *args))
         return effects
 
     def _update_header(self):

snsim/salt_utils.py

@@ -6,6 +6,11 @@
 from scipy.interpolate import RectBivariateSpline as spline2d
 
 
+def x0_to_mb(x0, source_name, source_version):
+    source = snc.get_source(source_name, version=source_version)
+    mb = -2.5 * np.log10(x0 / source.get('x0')) + source.peakmag('bessellb', 'ab')
+    return mb
+
 def n21_x1_model(z, seed=None):
     """X1 distribution redshift dependant model from  Nicolas et al. 2021.
 
@@ -263,3 +268,42 @@ def compute_salt_fit_error(fit_model, cov, band, time_th, zp, magsys="ab"):
     J = np.asarray([[d1, d2, d3] for d1, d2, d3 in zip(dfdx0, dfdx1, dfdc)])
     err_th = np.sqrt(np.einsum("ki,ki->k", J, np.einsum("ij,kj->ki", cov, J)))
     return err_th
+
+
+def salt_fitter(lc, fit_model, fit_par, **kwargs):
+    """Fit a given lightcurve with sncosmo salt model.
+
+    Parameters
+    ----------
+    lc : astropy.Table
+        The SN lightcurve.
+    fit_model : sncosmo.Model
+        Model used to fit the ligthcurve.
+    fit_par : list(str)
+        The parameters to fit.
+
+    Returns
+    -------
+    sncosmo.utils.Result (numpy.nan if no result)
+        sncosmo dict of fit results.
+
+    """
+    try:
+        res = snc.fit_lc(data=lc, model=fit_model, vparam_names=fit_par, **kwargs)
+        if res[0]["covariance"] is None:
+            res[0]["covariance"] = np.empty(
+                (len(res[0]["vparam_names"]), len(res[0]["vparam_names"]))
+            )
+            res[0]["covariance"][:] = np.nan
+
+        res[0]["param_names"] = np.append(res[0]["param_names"], "mb")
+
+        res[0]["parameters"] = np.append(
+            res[0]["parameters"], x0_to_mb(res[0]['parameters'][2], res[1].source.name, res[1].source.name)
+        )
+
+        res_dic = {k: v for k, v in zip(res[0]["param_names"], res[0]["parameters"])}
+        res = np.append(res, res_dic)
+    except (RuntimeError, snc.fitting.DataQualityError):
+        res = ["NaN", "NaN", "NaN"]
+    return res