1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
"""
===========================
Creating a new Source class
===========================
Extending sncosmo with a custom type of Source.
A ``Source`` is something that specifies a spectral timeseries as
a function of an arbitrary number of parameters. For example, the SALT2
model has three parameters (``x0``, ``x1`` and ``c``) that determine a
unique spectrum as a function of phase. The ``SALT2Source`` class implements
the behavior of the model: how the spectral time series depends on those
parameters.
If you have a spectral timeseries model that follows the behavior of one of
the existing classes, such as ``TimeSeriesSource``, great! There's no need to
write a custom class. However, suppose you want to implement a model that
has some new parameterization. In this case, you need a new class that
implements the behavior.
In this example, we implement a new type of source model. Our model is a linear
combination of two spectral time series, with a parameter ``w`` that
determines the relative weight of the models.
"""
import numpy as np
from scipy.interpolate import RectBivariateSpline
import sncosmo
class ComboSource(sncosmo.Source):
_param_names = ['amplitude', 'w']
param_names_latex = ['A', 'w'] # used in plotting display
def __init__(self, phase, wave, flux1, flux2, name=None, version=None):
self.name = name
self.version = version
self._phase = phase
self._wave = wave
# ensure that fluxes are on the same scale
flux2 = flux1.max() / flux2.max() * flux2
self._model_flux1 = RectBivariateSpline(phase, wave, flux1, kx=3, ky=3)
self._model_flux2 = RectBivariateSpline(phase, wave, flux2, kx=3, ky=3)
self._parameters = np.array([1., 0.5]) # initial parameters
def _flux(self, phase, wave):
amplitude, w = self._parameters
return amplitude * ((1.0 - w) * self._model_flux1(phase, wave) +
w * self._model_flux2(phase, wave))
########################################################################
# ... and that's all that we need to define!: A couple class attributes
# (``_param_names`` and ``param_names_latex``, an ``__init__`` method,
# and a ``_flux`` method. The ``_flux`` method is guaranteed to be passed
# numpy arrays for phase and wavelength.
#
# We can now initialize an instance of this source from two spectral time
# series:
#Just as an example, we'll use some undocumented functionality in
# sncosmo to download the Nugent Ia and 2p templates. Don't rely on this
# the `DATADIR` object, or these paths in your code though, as these are
# subject to change between version of sncosmo!
from sncosmo.builtins import DATADIR
phase1, wave1, flux1 = sncosmo.read_griddata_ascii(
DATADIR.abspath('models/nugent/sn1a_flux.v1.2.dat'))
phase2, wave2, flux2 = sncosmo.read_griddata_ascii(
DATADIR.abspath('models/nugent/sn2p_flux.v1.2.dat'))
# In our __init__ method we defined above, the two fluxes need to be on
# the same grid, so interpolate the second onto the first:
flux2_interp = RectBivariateSpline(phase2, wave2, flux2)(phase1, wave1)
source = ComboSource(phase1, wave1, flux1, flux2_interp, name='sn1a+sn2p')
##########################################################################
# We can get a summary of the Source we created:
print(source)
##########################################################################
# Get a spectrum at phase 10 for different parameters:
from matplotlib import pyplot as plt
wave = np.linspace(2000.0, 10000.0, 500)
for w in (0.0, 0.2, 0.4, 0.6, 0.8, 1.0):
source.set(w=w)
plt.plot(wave, source.flux(10., wave), label='w={:3.1f}'.format(w))
plt.legend()
plt.show()
##########################################################################
# The w=0 spectrum is that of the Ia model, the w=1 spectrum is that of
# the IIp model, while intermediate spectra are weighted combinations.
#
# We can even fit the model to some data!
model = sncosmo.Model(source=source)
data = sncosmo.load_example_data()
result, fitted_model = sncosmo.fit_lc(data, model,
['z', 't0', 'amplitude', 'w'],
bounds={'z': (0.2, 1.0),
'w': (0.0, 1.0)})
sncosmo.plot_lc(data, model=fitted_model, errors=result.errors)
##########################################################################
# The fact that the fitted value of w is closer to 0 than 1 indicates that
# the light curve looks more like the Ia template than the IIp template.
# This is generally what we expected since the example data here was
# generated from a Ia template (although not the Nugent template!).
|
