Skip to content

Add convolve_model function. #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
patti wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Add convolve_model function. #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
611308d
Add convolve_model function.
patti Aug 12, 2015
01c2eaf
Add required import statements for convolve_model function.
patti Aug 13, 2015
218c3e8
Add docs for model convolution.
patti Aug 18, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion astropy/convolution/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@

try:
# Not guaranteed available at setup time
from .convolve import convolve, convolve_fft
from .convolve import convolve, convolve_fft, convolve_model
except ImportError:
if not _ASTROPY_SETUP_:
raise
46 changes: 45 additions & 1 deletion astropy/convolution/convolve.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,10 @@
from .core import Kernel, Kernel1D, Kernel2D, MAX_NORMALIZATION
from ..utils.exceptions import AstropyUserWarning
from ..utils.console import human_file_size

from ..modeling.core import Model, _make_arithmetic_operator, \
BINARY_OPERATORS, _CompoundModelMeta
from scipy.signal import fftconvolve
from functools import partial

# Disabling all doctests in this module until a better way of handling warnings
# in doctests can be determined
Expand Down Expand Up @@ -563,3 +566,44 @@ def convolve_fft(array, kernel, boundary='fill', fill_value=0, crop=True,
return result
else:
return rifft.real

def convolve_model(model, kernel):
"""
Convolve a model with another model. Returns an
``astropy.modeling.CompoundModel`` instance.

Parameters
----------
model : `astropy.modeling.Model`
Source model to be convolved.
kernel : `astropy.modeling.Model`
Kernal model. Assumed to be normalized.

See Also
--------
convolve_fft, convolve

Returns
-------
model : astropy.modeling.CompoundModel
Returns a CompoundModel instance that convolves``model`` and ``kernel``
when evaluated.
"""

# Check model and kernel are Model instances
if not (isinstance(model, Model) & isinstance(kernel, Model)):
raise TypeError("Input model and kernal should be "
"`astropy.modeling.Model` instances.")

# Check that the number of dimensions is compatible
if model.n_inputs != kernel.n_inputs:
raise ValueError("Model and kernel must have same number of "
"dimensions")

convolve_func = partial(fftconvolve, mode='same')

BINARY_OPERATORS['convolve'] = _make_arithmetic_operator(convolve_func)

conv = _CompoundModelMeta._from_operator('convolve', model, kernel)

return conv
150 changes: 150 additions & 0 deletions docs/modeling/convolve-models.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,150 @@
.. _convolve-models:

Model Convolution
=================

.. versionadded:: 1.1

Models can be convolved using the :func:`~astropy.convolution.convolve_model` function. The result is a CompoundModel that performs the convolution on evaluation. For example a source model may be convolved with a instrument response function or point spread function. The convolved model can then be fit to data.

It is the user's responsibility to make sure the kernel is both centered and normalized as no errors or warnings will be raised. Usage and fitting are illustrated in the following examples

An Absorption Line Convolved with a Gaussian Instrument Response
----------------------------------------------------------------
The absorption line and Instrument Response Function (IRF) are each modeled by `astropy.modeling.functional_models.Gaussian1D`. In this example, the IRF parameters are fixed and the convolved compound model is fit to simulated data to recover the source model parameters.

.. plot::
:include-source:

import numpy as np
import matplotlib.pyplot as plot
from astropy.modeling import models, fitting
from astropy.convolution.convolve import convolve_model
from scipy.signal import fftconvolve

np.random.seed(100)

# Absorption Line Model
mean = 0
amp = -2
stddev = .01
src = models.Gaussian1D(amp, mean, stddev, name='Source')

# Normalized Gaussian PSF Model
irf_width = .3
amplitude = 1 / np.sqrt(2 * np.pi * irf_width ** 2)
irf = models.Gaussian1D(mean=mean, stddev=irf_width,
amplitude=amplitude, name='IRF')

# Fix the IRF parameters
for param in irf.param_names: irf.fixed[param] = True

# Create convolved CompoundModel object
cmod = convolve_model(src, irf).rename('Convolved')

# Simulate data image
# Note that the inputs are odd-shaped and centered on the kernel mean
x = np.linspace(-10, 10, 1001)
expected = cmod(x)
noise = np.random.normal(0, 1, (expected.shape))
cim = expected + noise

# Offset model parameters to test fitting
cmod.mean_0 = .1
cmod.stddev_0 = .02
cmod.amplitude_0 = -2.5

# Fit
fitter = fitting.LevMarLSQFitter()
fit_mod = fitter(cmod, x, cim)
fit = fit_mod(x)

# Plot
plt.figure(figsize=(16,12))
for i,(t,im) in enumerate([('Source', src(x)),
('IRF', irf(x)), ('Data', cim),
('Fit', fit), ('Residual', cim-fit),
('Actual - Fit', expected-fit)]):
plt.subplot(3, 2, i+1)
plt.plot(x, im)
if i>3:
plt.ylim(-.4, .4)
plt.title(t)
plt.xlabel('x')
plt.ylabel('Amplitude')

An Elliptical Galaxy convolved with a Gaussian PSF
--------------------------------------------------
The galaxy model is given by the de Vancouleur's profile using `astropy.modeling.functional_models.Sersic2D` with sersic index, ``n=4``. The Point Spread Function (PSF) is a `astropy.modeling.functional_models.Gaussian2D` model instance with fixed parameters. The convolved compound model is fit to simulated data to recover the source model parameters.

.. plot::
:include-source:

import numpy as np
import matplotlib.pyplot as plot
from astropy.modeling import models, fitting
from astropy.convolution.convolve import convolve_model
from scipy.signal import fftconvolve

np.random.seed(100)

# DeVanCoulier's Galaxy Source Model
x0, y0 = 50, 50
r_eff = 5000
n = 4
b0 = .1
ellip = .5
theta = np.pi / 3
src = models.Sersic2D(b0, r_eff, n, x0, y0, ellip, theta, name='Source')

# Normalized Gaussian PSF Model
psf_width = 5
amplitude = 1 / (2 * np.pi * psf_width ** 2)
psf = models.Gaussian2D(x_mean = x0, y_mean=y0,
x_stddev=psf_width, y_stddev=psf_width,
amplitude=amplitude,
name='PSF')
for param in psf.param_names: psf.fixed[param] = True

# Create convolved CompoundModel object
cmod = convolve_model(src, psf).rename('Convolved')

# Simulate data image
# Note that the inputs are odd-shaped and centered on the kernel mean
y,x = np.indices((101, 101))
expected = cmod(x, y)
noise = np.random.normal(0, .1, (expected.shape))
cim = expected + noise

# Offset model parameters to test fitting
cmod.amplitude_0 = .15
cmod.r_eff_0 = 5500
cmod.n_0 = 3.6
cmod.ellip_0 = .64
cmod.theta_0 = np.pi/4
cmod.x_0 = 45
cmod.y_0=55

# Fit
fitter = fitting.LevMarLSQFitter()
fit_mod = fitter(cmod, x, y, cim)
fit = fit_mod(x, y)

# Plot
plt.figure(figsize=(16,8))
for i,(t,im) in enumerate([('Source',src(x, y)),
('PSF', psf(x, y)), ('Data', cim),
('Fit', fit), ('Residual', cim-fit),
('Actual - Fit', expected-fit)]):
plt.subplot(2, 3, i+1)
if i in [0, 2, 3]:
vmin, vmax = 0, 60
elif i==1:
vmin, vmax = 0, psf.amplitude.value
else:
vmin, vmax = -5, 5
plt.imshow(im, origin='lower', vmin=vmin, vmax=vmax)
plt.colorbar()
plt.title(t)
plt.xlabel('x')
plt.ylabel('y')
6 changes: 5 additions & 1 deletion docs/modeling/index.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -277,7 +277,7 @@ The resulting object ``g1_plus_2`` is itself a new model. Evaluating, say,

This model can be further combined with other models in new expressions. It is
also possible to define entire new model *classes* using arithmetic expressions
of other model classes. This allows general compound models to be created
of other model classes. This allows general compound models to be created
without specifying any parameter values up front. This more advanced usage is
explained in more detail in the :ref:`compound model documentation
<compound-model-classes>`.
Expand Down Expand Up @@ -315,6 +315,9 @@ are some complexities involved in correctly matching up the inputs and outputs
of all models used to build a compound model. You can learn more details in
the :doc:`compound-models` documentation.

Likewise, models can be convolved to form a new compound model using the
:func:`~astropy.convolution.convolve_model` function. See
:doc:`convolve-models` for examples.

Using `astropy.modeling`
========================
Expand All @@ -326,6 +329,7 @@ Using `astropy.modeling`
parameters
fitting
compound-models
convolve-models
new
bounding-boxes
algorithms
Expand Down