adding G10 and C11 intrinsic scattering models

docs/models.rst

@@ -246,6 +246,47 @@ this value is perfectly known from the dust map. Therefore, when using
 a function such as `~sncosmo.fit_lc` to fit the parameters, be sure *not* to
 include ``'mwebv'`` in the list of parameters to vary.
 
+Adding color dependant scatter model
+====================================
+
+The intrinsic scattering of SNe Ia is color dependant it can be modelled for simulation purpose
+by G10 or C11 models. The implemention is based on Kessler et al. 2012.
+They both act as random variation in the spectra model of a `~sncosmo.SALT2Source` or `~sncosmo.SALT3Source`.
+
+The G10 model take a `SALTSource` as an argument and can be added to your `~sncosmo.Model` as:
+
+.. code:: python
+
+    >>> source = 'salt2'
+    >>> SALTSource = sncosmo.models.get_source(name=source)
+    >>> G10 = snc.models.G10(SALTsource=SALTSource)
+    >>> SALTwithG10 = sncosmo.Model(source='salt2',
+                                    effects=[G10],
+                                    effect_names=['G10'],
+                                    effect_frames=['rest'])
+
+The G10 model parameters are:   
+
+* ``L0``, ``F0`` and ``F1`` are used in the multiplicative factor introduced in Kessler et al. 2012. Their default values are ``L0=2157.3``, ``F0=0`` and ``F1=1.08e-4``.
+* ``dL`` the wavelength range between each scatter node. Following Kessler et al. 2012 it is set by default to 800A.
+
+The C11 model can be added to your `~sncosmo.Model` as:
+
+.. code:: python
+
+    >>> C11 = snc.models.C11()
+    >>> SALTwithC11 = sncosmo.Model(source='salt2',
+                                    effects=[C11],
+                                    effect_names=['C11'],
+                                    effect_frames=['rest'])
+
+The C11 model parameters are:
+
+* ``CvU`` the correlation coefficient between the v and U band that could be -1, 0 or 1.
+* ``Sf`` a scale factor fixed by default to ``S_f=1.3`` according to Kessler et al. 2012.
+
+
+
 Phase Dependant effects
 =======================
 
@@ -373,4 +414,4 @@ parameters:
     >>> source = sncosmo.SALT2Source(modeldir='/path/to/dir',
     ...                              m0file='mytemplate0file.dat')
 
-See `~sncosmo.SALT2Source` for more details.
+See  for more details.

sncosmo/models.py

@@ -26,7 +26,7 @@
 )
 from .magsystems import get_magsystem
 from .salt2utils import BicubicInterpolator, SALT2ColorLaw
-from .utils import integration_grid
+from .utils import integration_grid, sine_interp
 
 __all__ = ['get_source', 'Source', 'TimeSeriesSource', 'StretchSource',
            'SUGARSource', 'SALT2Source', 'SALT3Source', 'MLCS2k2Source',
@@ -2030,3 +2030,114 @@ def propagate(self, wave, flux, phase=None):
         """Propagate the flux."""
         ebv = self._parameters[0]
         return extinction.apply(self._f(wave, ebv * self._r_v), flux)
+
+
+class G10(PropagationEffect):
+    """Guy (2010) SNe Ia non-coherent scattering.
+
+    Implementation is done following arxiv:1209.2482."""
+
+    _param_names = ['L0', 'F0', 'F1', 'dL']
+    param_names_latex = [r'\lambda_0', 'F_0', 'F_1', 'd_L']
+
+    def __init__(self, SALTsource):
+        """Initialize G10 class."""
+        self._parameters = np.array([2157.3, 0.0, 1.08e-4, 800])
+        self._colordisp = SALTsource._colordisp
+        self._minwave = SALTsource.minwave()
+        self._maxwave = SALTsource.maxwave()
+
+        self._seed = np.random.SeedSequence()
+
+
+    def compute_sigma_nodes(self):
+        """Computes the sigma nodes."""
+        L0, F0, F1, dL = self._parameters
+        lam_nodes = np.arange(self._minwave, self._maxwave, dL)
+        if lam_nodes.max() < self._maxwave:
+            lam_nodes = np.append(lam_nodes, self._maxwave)
+        siglam_values = self._colordisp(lam_nodes) 
+
+        siglam_values[lam_nodes < L0] *= 1 + (lam_nodes[lam_nodes < L0] - L0) * F0
+        siglam_values[lam_nodes > L0] *= 1 + (lam_nodes[lam_nodes > L0] - L0) * F1
+
+        return lam_nodes, siglam_values
+
+    def propagate(self, wave, flux):
+        """Propagate the effect to the flux."""# Draw the scattering
+        lam_nodes, siglam_values = self.compute_sigma_nodes()
+        siglam_values *= np.random.default_rng(self._seed).normal(size=len(lam_nodes))
+        magscat = sine_interp(wave, lam_nodes, siglam_values)
+        return flux * 10**(-0.4 * magscat)
+
+
+class C11(PropagationEffect):
+    """C11 scattering effect for sncosmo.
+    Use COV matrix between the vUBVRI bands from N. Chottard thesis.
+        Implementation is done following arxiv:1209.2482."""
+
+    _param_names = ["CvU", 'Sf']
+    param_names_latex = ["\rho_\mathrm{vU}", 'S_f']
+    _minwave = 2000
+    _maxwave = 11000
+
+    def __init__(self):
+        """Initialise C11 class."""
+        self._parameters = np.array([0., 1.3])
+
+        # vUBVRI lambda eff
+        self._lam_nodes = np.array([2500.0, 3560.0, 4390.0, 5490.0, 6545.0, 8045.0])
+
+        # vUBVRI correlation matrix extract from SNANA, came from N.Chotard thesis
+        self._corr_matrix = np.array(
+            [
+                [+1.000000,  0.000000,  0.000000,  0.000000,  0.000000,  0.000000],
+                [ 0.000000, +1.000000, -0.118516, -0.768635, -0.908202, -0.219447],
+                [ 0.000000, -0.118516, +1.000000, +0.570333, -0.238470, -0.888611],
+                [ 0.000000, -0.768635, +0.570333, +1.000000, +0.530320, -0.399538],
+                [ 0.000000, -0.908202, -0.238470, +0.530320, +1.000000, +0.490134],
+                [ 0.000000, -0.219447, -0.888611, -0.399538, +0.490134, +1.000000]
+            ]
+            ) 
+
+        # vUBVRI sigma
+        self._variance = np.array([0.5900, 0.06001, 0.040034, 0.050014, 0.040017, 0.080007])
+
+        self._seed = np.random.SeedSequence()
+
+    def build_cov(self):
+        CvU, Sf = self._parameters
+
+        cov_matrix = self._corr_matrix.copy()
+
+        # Set up the vU correlation
+        cov_matrix[0, 1:] = CvU * self._corr_matrix[1, 1:]
+        cov_matrix[1:, 0] = CvU * self._corr_matrix[1:, 1]
+
+        # Convert corr to cov
+        cov_matrix *= np.outer(self._variance, 
+                               self._variance)
+
+        # Rescale covariance as in arXiv:1209.2482
+        cov_matrix *= Sf
+        return cov_matrix
+
+    def propagate(self, wave, flux):
+        """Propagate the effect to the flux."""
+
+        cov_matrix = self.build_cov()
+
+        # Draw the scattering
+        siglam_values = np.random.default_rng(self._seed).multivariate_normal(np.zeros(len(self._lam_nodes)), 
+                                                                              cov_matrix)
+
+        inf_mask = wave <= self._lam_nodes[0]
+        sup_mask = wave >= self._lam_nodes[-1]
+
+        magscat = np.zeros(len(wave))
+        magscat[inf_mask] = siglam_values[0]
+        magscat[sup_mask] = siglam_values[-1]
+        magscat[~inf_mask & ~sup_mask] = sine_interp(wave[~inf_mask & ~sup_mask], 
+                                                     self._lam_nodes, siglam_values)
+
+        return flux * 10**(-0.4 * magscat)

sncosmo/tests/test_models.py

@@ -367,3 +367,62 @@ def test_effect_phase_dependent():
                                 effect_frames=['rest'],
                                 effect_names=['micro'])
     assert_approx_equal(model_micro.flux(50., 5000.).flatten(), flux*10.)
+
+
+def test_G10():
+    """Test Model with G10 color dependant scatter"""
+
+    SALT2Source = sncosmo.models.get_source('salt2', version='2.4')
+    ModelRef = sncosmo.Model(SALT2Source)
+
+    G10 = sncosmo.models.G10(SALT2Source)
+    ModelWithG10 = sncosmo.Model(source=SALT2Source,
+                                effects=[G10],
+                                effect_frames=['rest'],
+                                effect_names=['G10'])
+
+    lam_nodes, siglam_values = G10.compute_sigma_nodes()
+
+    # Test how nodes are computed
+    assert_allclose(lam_nodes, np.array([2000., 2800., 3600., 4400., 
+                                            5200., 6000., 6800., 7600.,
+                                            8400., 9200.])) 
+
+    # Test how siglam values are computed
+    assert_allclose(siglam_values, np.array([1.308910000, 0.259717301, 0.078368072, 0.035382907, 
+                                            0.023921785, 0.024232781, 0.036799298 , 0.083808031, 
+                                            0.286347107, 1.468232113]))
+
+
+    # Compare with and without G10
+    random = np.random.default_rng(G10._seed).normal(size=len(lam_nodes))
+    assert_allclose(ModelWithG10.flux(0, lam_nodes), 
+                    ModelRef.flux(0, lam_nodes) * 10**(-0.4 * siglam_values * random))
+
+
+def test_C11():
+    """Test Model with C11 color dependant scatter"""
+    SALT2Source = sncosmo.models.get_source('salt2', version='2.4')
+    ModelRef = sncosmo.Model(SALT2Source)
+
+    C11 = sncosmo.models.C11()
+    ModelWithC11 = sncosmo.Model(source=SALT2Source,
+                                effects=[C11],
+                                effect_frames=['rest'],
+                                effect_names=['C11_'])
+
+    for CvU in [-1, 0, 1]:
+        ModelWithC11.set(C11_CvU=CvU, C11_Sf=1)
+        cov = ModelWithC11.effects[0].build_cov() 
+
+        # Test construvtion of cov matrix
+        corr = cov / np.outer(C11._variance, C11._variance)
+        assert_allclose(corr[0, 1:], CvU * C11._corr_matrix[1, 1:])
+        assert_allclose(corr[1:, 0], CvU * C11._corr_matrix[1, 1:])
+        assert_allclose(corr[1:, 1:], C11._corr_matrix[1:, 1:])
+
+        # Compare to model without C11 
+        random = np.random.default_rng(C11._seed).multivariate_normal(np.zeros(len(C11._lam_nodes)), 
+                                                                      cov)
+        assert_allclose(ModelWithC11.flux(0, C11._lam_nodes), 
+                        ModelRef.flux(0, C11._lam_nodes) * 10**(-0.4 * random))

sncosmo/utils.py

@@ -539,3 +539,25 @@ def warn_once(name, depver, rmver, extra=None):
             msg += " " + extra
         warnings.warn(msg, stacklevel=2)
         warned.append(name)
+
+
+def sine_interp(x_new, fun_x, fun_y):
+    """Sinus interpolation for intrinsic scattering models."""
+    if len(fun_x) != len(fun_y):
+        raise ValueError('x and y must have the same len')
+    if (x_new > fun_x[-1]).any() or (x_new < fun_x[0]).any():
+        raise ValueError('x_new is out of range of fun_x')
+
+    sup_bound = np.vstack([x_new >= x for x in fun_x])
+
+    idx_inf = np.sum(sup_bound, axis=0) - 1
+    idx_inf[idx_inf==len(fun_x) - 1] = -2
+
+    x_inf = fun_x[idx_inf]
+    x_sup = fun_x[idx_inf + 1]
+    fun_y_inf = fun_y[idx_inf]
+    fun_y_sup = fun_y[idx_inf + 1]
+
+    sin_interp = np.sin(np.pi * (x_new - 0.5 * (x_inf + x_sup)) / (x_sup - x_inf))
+    values = 0.5 * (fun_y_sup + fun_y_inf) + 0.5 * (fun_y_sup - fun_y_inf) * sin_interp
+    return values