adding G10 and C11 intrinsic scattering models

[bastiencarreres]

Hi, I propose this implementation for the G10 scattering model, it's a cleaner implementation of one I use in snsim. It follows the implementation used in SNANA as described in arxiv:1209.2482.

I create the G10 class in models.py that inherit from PropagationEffect. To be initialized the class need a SALTsource:

import sncosmo as snc

SALTSource = snc.models.get_source(name='salt2')
G10 = snc.models.G10(Source)
SALTwithG10 = snc.Model('salt2')
SALTwithG10.add_effect(G10, name='G10', frame='rest')

The compute_sigma_nodes class function compute a sigma value from source.minwave() to source.maxwave() each 800 A, then for each sigma an intrinsic scattering is drawn with gaussian pdf. The result is two arrays : the wavelength at values lam_nodes and the scattering value for each wavelength siglam_values.

In propagate class function the wavelengths are interpolated by a sinus function on (lam_nodes , siglam_values).

Possible problem: The Gaussian drawing is done each time propagate is called, resulting in different results for different calls of SALTwithG10.bandflux.

[MickaelRigault]

I would like to comment, that I like that very much.
We should also have others like C11 etc. and an example on how to build your own.
Thanks @bastiencarreres

[bastiencarreres]

Hi @MickaelRigault, I also have the C11 but need some extra-test before push I think

[bastiencarreres]

I modfied the PR, I've added the C11 class in model.py.

import sncosmo as snc

C11 = snc.models.C11()
ModelC11 = snc.Model('salt3')
ModelC11.add_effect(C11, name='C11',frame='rest')

I followed the SNANA implementation also described in arxiv:1209.2482.
I added a sine_interp(x_new, fun_x, fun_y) function in utils.py. This function is used by G10 and C11.

The implementation of C11 works as follows:

• There is 2 parameters inside self._parameters : the first one is C_vU is the correlation between the v and U bands and is usually -1, 0 or 1(maybe need to add an error on this); the second parameter S_f is a scaling factor for the covariance fixed to 1.3 in arxiv:1209.2482.
• The self._lam_nodes are defined as the central wavelengths of vUBVRI bands.
• The correlation matrix from Chotard et al. 2011 between vUBVRI bands self._corr_matrix is converted into a covariance self._cov_matrix matrix by multiplying by sigma values self._siglam_values.
• The correlation matrix self._cov_matrix is then rescale by the scaling factor S_f
• When propagate(self, wave, flux) is called siglam_values are drawn on a multivariate normal centered on 0 and with the defined covariance matrix, then the wave are interpolated with sine_interp on (self.lam_nodes, siglam_values), wavelengths outside the vUBVRI range are interpolated as the siglam_values[0] or siglam_values[-1] constant.

As for G10 the Gaussian drawing is done each time propagate is called, giving different results at each call.

[MickaelRigault]

@benjaminrose

Any update on this ? Could that be accepted or you think it should not be made like that ?

[benjaminrose]

Sorry. I started a new job and am slowly getting settled in. I'll work on this soon.

[benjaminrose]

After a quick review, SALTwithG10.bandflux should be repeatable and not vary by call. I would like to see an example in the docs. Either in the model & effects section or in the simulation section (or both?). Also I want to run the test suite. Can you add a few tests similar to the ones found here

sncosmo/sncosmo/tests/test_models.py

Lines 342 to 369 in 9037060

	def test_effect_frame_free():
	"""Test Model with PropagationEffect with a 'free' frame."""
	model = sncosmo.Model(source=flatsource(),
	effects=[StepEffect(800., 20000.)],
	effect_frames=['free'],
	effect_names=['screen'])
	model.set(z=1.0, screenz=0.5, screenstepwave=10000.)
	# maxwave is set by the effect at z=0.5
	assert model.maxwave() == (1. + 0.5) * 20000.
	wave = np.array([12000., 13000., 14000., 14999., 15000., 16000.])
	assert_allclose(model.flux(0., wave), [0., 0., 0., 0., 0.5, 0.5])
	def test_effect_phase_dependent():
	"""Test Model with PropagationEffect with phase dependence"""
	model = sncosmo.Model(source=flatsource())
	flux = model.flux(50., 5000.).flatten()
	model_micro = sncosmo.Model(source=flatsource(),
	effects=[AchromaticMicrolensing([40., 60.],
	[10., 10.])],
	effect_frames=['rest'],
	effect_names=['micro'])
	assert_approx_equal(model_micro.flux(50., 5000.).flatten(), flux*10.)

.

[benjaminrose]

Would moving this to __init__ make this SALTwithG10.bandflux the same on repeated calls?

[bastiencarreres]

Yes it should works since the random part is in self.compute_sigma_nodes()

[bastiencarreres]

Moving this part in __init__ leads to a non-functional Model.set() since the updates on self._parameters will not be taken into account. I set a self._seed parameter in the __init__ to fix the randomness when init the model.

[benjaminrose]

I think these first three statements could be moved to __init__. Am I missing something?

[bastiencarreres]

With respect to @benjaminrose comments, I moved all random variables in __init__.

[bastiencarreres]

I'va added tests for G10 and C11 in test_models.py as two test functions test_G10 and test_C11. I've also done some debugging relative to the previous changes that prevented Model.set() from working.

[benjaminrose]

#373

[benjaminrose]

Sorry this is taking a long time. This is a fairly major change in functionality and I want to use it as a template for future additions. I think we are really close to something extremely useful.

[benjaminrose]

There are two issues found in the tests. 1) the tests don't work on python 3.7 and we should discuss what to do to fix that. and 2) there are some code style issues.

[bastiencarreres]

I've corrected the code style. For python 3.7 I think the problem comes from np.random.SeedSequence(). A way to avoid that is to use the standard library function secrets.randbits as shown here.

[benjaminrose]

Python 3.7 is also end of life, https://devguide.python.org/versions/, maybe we should just bump our minimum dependancies?

[bastiencarreres]

I've re-corrected style. Yes Python 3.7 is end of life and not supported by numpy since 1.22 so changing minimum dependencies is a good solution !

[benjaminrose]

Let me drop 3.7 support, then we should be able to merge this in.