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 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
|
.. _modeling-units:
********************************
Support for units and quantities
********************************
.. note:: The functionality presented here was recently added. If you run into
any issues, please don't hesitate to open an issue in the `issue
tracker <https://github.com/astropy/astropy/issues>`_.
The `astropy.modeling` package includes partial support for the use of units and
quantities in model parameters, models, and during fitting. At this time, only
some of the built-in models (such as
:class:`~astropy.modeling.functional_models.Gaussian1D`) support units, but this
will be extended in future to all models where this is appropriate.
Setting parameters to quantities
================================
Models can take :class:`~astropy.units.Quantity` objects as parameters::
>>> from astropy import units as u
>>> from astropy.modeling.models import Gaussian1D
>>> g1 = Gaussian1D(mean=3 * u.m, stddev=2 * u.cm, amplitude=3 * u.Jy)
Accessing the parameter then returns a Parameter object that contains the value
and the unit::
>>> g1.mean
Parameter('mean', value=3.0, unit=m)
It is then possible to access the individual properties of the parameter::
>>> g1.mean.name
'mean'
>>> g1.mean.value
3.0
>>> g1.mean.unit
Unit("m")
If a parameter has been initialized as a Quantity, it should always be set to a
quantity, but the units don't have to be compatible with the initial ones::
>>> g1.mean = 3 * u.s
>>> g1 # doctest: +FLOAT_CMP
<Gaussian1D(amplitude=3. Jy, mean=3. s, stddev=2. cm)>
To change the value of a parameter and not the unit, simply set the value
property::
>>> g1.mean.value = 2
>>> g1 # doctest: +FLOAT_CMP
<Gaussian1D(amplitude=3. Jy, mean=2. s, stddev=2. cm)>
Setting a parameter which was originally set to a quantity to a scalar doesn't
work because it's ambiguous whether the user means to change just the value and
preserve the unit, or get rid of the unit::
>>> g1.mean = 2 # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
UnitsError : The 'mean' parameter should be given as a Quantity because it
was originally initialized as a Quantity
On the other hand, if a parameter previously defined without units is given a
Quantity with a unit, this works because it is unambiguous::
>>> g2 = Gaussian1D(mean=3)
>>> g2.mean = 3 * u.m
In other words, once units are attached to a parameter, they can't be removed
due to ambiguous meaning.
Evaluating models with quantities
=================================
Quantities can be passed to model during evaluation::
>>> g3 = Gaussian1D(mean=3 * u.m, stddev=5 * u.cm)
>>> g3(2.9 * u.m) # doctest: +FLOAT_CMP
<Quantity 0.1353352832366122>
>>> g3(2.9 * u.s) # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
UnitsError : Units of input 'x', s (time), could not be converted to
required input units of m (length)
In this case, since the mean and standard deviation have units, the value passed
during evaluation also needs units::
>>> g3(3) # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
UnitsError : Units of input 'x', (dimensionless), could not be converted to
required input units of m (length)
Equivalencies
-------------
Equivalencies require special care - a Gaussian defined in frequency space is
not a Gaussian in wavelength space for example. For this reason, we don't allow
equivalencies to be attached to the parameters themselves. Instead, we take the
approach of converting the input data to the parameter space, and any
equivalencies should be applied at evaluation time to the data (not the
parameters).
Let's consider a model that is Gaussian in wavelength space::
>>> g4 = Gaussian1D(mean=3 * u.micron, stddev=1 * u.micron, amplitude=3 * u.Jy)
By default, passing a frequency will not work:
>>> g4(1e2 * u.THz) # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
UnitsError : Units of input 'x', THz (frequency), could not be converted to
required input units of micron (length)
But you can pass a dictionary of equivalencies to the equivalencies argument
(this needs to be a dictionary since some models can contain multiple inputs)::
>>> g4(110 * u.THz, equivalencies={'x': u.spectral()}) # doctest: +FLOAT_CMP
<Quantity 2.888986819525229 Jy>
The key of the dictionary should be the name of the inputs according to::
>>> g4.inputs
('x',)
It is also possible to set default equivalencies for the input parameters using
the input_units_equivalencies property::
>>> g4.input_units_equivalencies = {'x': u.spectral()}
>>> g4(110 * u.THz) # doctest: +FLOAT_CMP
<Quantity 2.888986819525229 Jy>
Fitting models with units to data
=================================
Fitting models with units to data with units should be seamless provided that
the model supports fitting with units. To demonstrate this, we start off by
generating synthetic data:
.. plot::
:context: reset
:include-source:
import numpy as np
from astropy import units as u
import matplotlib.pyplot as plt
x = np.linspace(1, 5, 30) * u.micron
y = np.exp(-0.5 * (x - 2.5 * u.micron)**2 / (200 * u.nm)**2) * u.mJy
plt.plot(x, y, 'ko')
plt.xlabel('Wavelength (microns)')
plt.ylabel('Flux density (mJy)')
and we then define the initial guess for the fitting and we carry out the fit as
we would without any units:
.. plot::
:context:
:include-source:
from astropy.modeling import models, fitting
g5 = models.Gaussian1D(mean=3 * u.micron, stddev=1 * u.micron, amplitude=1 * u.Jy)
fitter = fitting.LevMarLSQFitter()
g5_fit = fitter(g5, x, y)
plt.plot(x, y, 'ko')
plt.plot(x, g5_fit(x), 'r-')
plt.xlabel('Wavelength (microns)')
plt.ylabel('Flux density (mJy)')
Fitting with equivalencies
--------------------------
Let's now consider the case where the data is not equivalent to those of the
parameters, but they are convertible via equivalencies. In this case, the
equivalencies can either be passed via a dictionary as shown higher up for the
evaluation examples:
.. plot::
:context:
:include-source:
g6 = models.Gaussian1D(mean=110 * u.THz, stddev=10 * u.THz, amplitude=1 * u.Jy)
g6_fit = fitter(g6, x, y, equivalencies={'x': u.spectral()})
plt.plot(x, g6_fit(x, equivalencies={'x': u.spectral()}), 'b-')
plt.xlabel('Wavelength (microns)')
plt.ylabel('Flux density (mJy)')
In this case, the fit (in blue) is slightly worse, because a Gaussian in
frequency space (blue) is not a Gaussian in wavelength space (red). As mentioned
previously, you can also set input_units_equivalencies on the model itself to
avoid having to pass extra arguments to the fitter::
g6.input_units_equivalencies = {'x': u.spectral()}
g6_fit = fitter(g6, x, y)
.. _units-mapping:
Support for units in otherwise unitless models
==============================================
Some models, like polynomials, do not work intrinsically with units. Instead,
the :meth:`~astropy.modeling.core.Model.coerce_units` method provides a way to add input and return units to
unitless models by enclosing the unitless model with two instances of :class:`~astropy.modeling.mappings.UnitsMapping`.
Internally the inputs are stripped of the units before passed
to the model and units are attached to the result if ``return_units`` is specified.
The method returns a new composite model::
>>> from astropy.modeling import models
>>> from astropy import units as u
>>> model = models.Polynomial1D(1, c0=1, c1=2)
>>> new_model = model.coerce_units(input_units={'x': u.Hz}, return_units={'y': u.s},
... input_units_equivalencies={'x':u.spectral()})
>>> new_model(10 * u.Hz)
<Quantity 21. s>
|