From 611308d58d69f2c6feaec08422f1a32c91029757 Mon Sep 17 00:00:00 2001 From: Patti Carroll Date: Tue, 11 Aug 2015 18:40:00 -0700 Subject: [PATCH 1/3] Add convolve_model function. --- astropy/convolution/convolve.py | 45 +++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/astropy/convolution/convolve.py b/astropy/convolution/convolve.py index 4d346e52acae..426c7908e089 100644 --- a/astropy/convolution/convolve.py +++ b/astropy/convolution/convolve.py @@ -563,3 +563,48 @@ 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 a + `astropy.modeling.CompoundModel` object. + + 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.model.CompoundModel` + Returns a CompoundModel instance that convolves``model`` and ``kernel`` + when evaluated. + + Examples + -------- + + """ + + # 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 From 01c2eafc2ef77764848735e547a4973e5ecaf2db Mon Sep 17 00:00:00 2001 From: Patti Carroll Date: Wed, 12 Aug 2015 19:11:37 -0700 Subject: [PATCH 2/3] Add required import statements for convolve_model function. --- astropy/convolution/convolve.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/astropy/convolution/convolve.py b/astropy/convolution/convolve.py index 426c7908e089..b7b364571421 100644 --- a/astropy/convolution/convolve.py +++ b/astropy/convolution/convolve.py @@ -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 From 218c3e8b8c40acdb5ab8709172c8c8ea09af19ac Mon Sep 17 00:00:00 2001 From: Patti Carroll Date: Tue, 18 Aug 2015 12:46:52 -0700 Subject: [PATCH 3/3] Add docs for model convolution. --- astropy/convolution/__init__.py | 2 +- astropy/convolution/convolve.py | 10 +- docs/modeling/convolve-models.rst | 150 ++++++++++++++++++++++++++++++ docs/modeling/index.rst | 6 +- 4 files changed, 159 insertions(+), 9 deletions(-) create mode 100644 docs/modeling/convolve-models.rst diff --git a/astropy/convolution/__init__.py b/astropy/convolution/__init__.py index a188e1d5f741..3fd4ff565c59 100644 --- a/astropy/convolution/__init__.py +++ b/astropy/convolution/__init__.py @@ -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 diff --git a/astropy/convolution/convolve.py b/astropy/convolution/convolve.py index b7b364571421..ffd23d7b8c7a 100644 --- a/astropy/convolution/convolve.py +++ b/astropy/convolution/convolve.py @@ -569,8 +569,8 @@ def convolve_fft(array, kernel, boundary='fill', fill_value=0, crop=True, def convolve_model(model, kernel): """ - Convolve a model with another model. Returns a - `astropy.modeling.CompoundModel` object. + Convolve a model with another model. Returns an + ``astropy.modeling.CompoundModel`` instance. Parameters ---------- @@ -585,13 +585,9 @@ def convolve_model(model, kernel): Returns ------- - model : `astropy.model.CompoundModel` + model : astropy.modeling.CompoundModel Returns a CompoundModel instance that convolves``model`` and ``kernel`` when evaluated. - - Examples - -------- - """ # Check model and kernel are Model instances diff --git a/docs/modeling/convolve-models.rst b/docs/modeling/convolve-models.rst new file mode 100644 index 000000000000..ea1c296b6f87 --- /dev/null +++ b/docs/modeling/convolve-models.rst @@ -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') diff --git a/docs/modeling/index.rst b/docs/modeling/index.rst index f8c050e7d1cc..b794d011627d 100644 --- a/docs/modeling/index.rst +++ b/docs/modeling/index.rst @@ -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 `. @@ -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` ======================== @@ -326,6 +329,7 @@ Using `astropy.modeling` parameters fitting compound-models + convolve-models new bounding-boxes algorithms