Skip to content

Commit

Permalink
correct style
Browse files Browse the repository at this point in the history
  • Loading branch information
@bastiencarreres
bastiencarreres committed Feb 22, 2024
1 parent 3beb6c4 commit 63b56b9
Show file tree
Hide file tree
Showing 3 changed files with 64 additions and 51 deletions.
64 changes: 36 additions & 28 deletions sncosmo/models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2034,7 +2034,7 @@ def propagate(self, wave, flux, phase=None):

class G10(PropagationEffect):
"""Guy (2010) SNe Ia non-coherent scattering.
Implementation is done following arxiv:1209.2482."""

_param_names = ['L0', 'F0', 'F1', 'dL']
Expand All @@ -2048,7 +2048,6 @@ def __init__(self, SALTsource):
self._maxwave = SALTsource.maxwave()

self._seed = np.random.SeedSequence()


def compute_sigma_nodes(self):
"""Computes the sigma nodes."""
Expand All @@ -2057,16 +2056,21 @@ def compute_sigma_nodes(self):
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

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

Expand All @@ -2086,38 +2090,40 @@ def __init__(self):
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])
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
# vUBVRI correlation matrix extract from SNANA 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]
]
[
[+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._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
Expand All @@ -2128,16 +2134,18 @@ def propagate(self, wave, 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)
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)

magscat[~inf_mask & ~sup_mask] = sine_interp(
wave[~inf_mask & ~sup_mask],
self._lam_nodes, siglam_values)

return flux * 10**(-0.4 * magscat)
37 changes: 20 additions & 17 deletions sncosmo/tests/test_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -377,27 +377,29 @@ def test_G10():

G10 = sncosmo.models.G10(SALT2Source)
ModelWithG10 = sncosmo.Model(source=SALT2Source,
effects=[G10],
effect_frames=['rest'],
effect_names=['G10'])
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.]))
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]))

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))
assert_allclose(ModelWithG10.flux(0, lam_nodes),
ModelRef.flux(0, lam_nodes) *
10**(-0.4 * siglam_values * random))


def test_C11():
Expand All @@ -407,9 +409,9 @@ def test_C11():

C11 = sncosmo.models.C11()
ModelWithC11 = sncosmo.Model(source=SALT2Source,
effects=[C11],
effect_frames=['rest'],
effect_names=['C11_'])
effects=[C11],
effect_frames=['rest'],
effect_names=['C11_'])

for CvU in [-1, 0, 1]:
ModelWithC11.set(C11_CvU=CvU, C11_Sf=1)
Expand All @@ -422,7 +424,8 @@ def test_C11():
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)
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))
ModelRef.flux(0, C11._lam_nodes) * 10**(-0.4 * random))
14 changes: 8 additions & 6 deletions sncosmo/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -549,15 +549,17 @@ def sine_interp(x_new, fun_x, fun_y):
raise ValueError('x_new is out of range of fun_x')

Check warning on line 549 in sncosmo/utils.py

View check run for this annotation

Codecov / codecov/patch

sncosmo/utils.py#L549 

Added line #L549 was not covered by tests

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

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

0 comments on commit 63b56b9

Please sign in to comment.