astropy:docs

Source code for astropy.modeling.functional_models

# Licensed under a 3-clause BSD style license - see LICENSE.rst

"""Mathematical models."""

from __future__ import division

import collections

from textwrap import dedent

import numpy as np

from .core import (ParametricModel, Parametric1DModel, Parametric2DModel,
                   Model, format_input, ModelDefinitionError)
from .parameters import Parameter, InputParameterError
from ..utils import find_current_module


__all__ = sorted([
    'AiryDisk2D', 'Beta1D', 'Beta2D', 'Box1D',
    'Box2D', 'Const1D', 'Const2D', 'Disk2D',
    'Gaussian1D', 'Gaussian2D', 'Linear1D', 'Lorentz1D',
    'MexicanHat1D', 'MexicanHat2D', 'Scale', 'Shift',
    'Sine1D', 'Trapezoid1D', 'TrapezoidDisk2D', 'Ring2D',
    'custom_model_1d'
])


[docs]class Gaussian1D(Parametric1DModel): """ One dimensional Gaussian model. Parameters ---------- amplitude : float Amplitude of the gaussian mean : float Mean of the gaussian stddev : float Standard deviation of the gaussian Notes ----- Model formula: .. math:: f(x) = A e^{- \\frac{\\left(x - x_{0}\\right)^{2}}{2 \\sigma^{2}}} Examples -------- >>> from astropy.modeling import models >>> def tie_center(model): ... mean = 50 * model.stddev ... return mean >>> tied_parameters = {'mean': tie_center} Specify that 'mean' is a tied parameter in one of two ways: >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3, ... tied=tied_parameters) or >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3) >>> g1.mean.tied False >>> g1.mean.tied = tie_center >>> g1.mean.tied <function tie_center at 0x...> Fixed parameters: >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3, ... fixed={'stddev': True}) >>> g1.stddev.fixed True or >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3) >>> g1.stddev.fixed False >>> g1.stddev.fixed = True >>> g1.stddev.fixed True See Also -------- Gaussian2D, Box1D, Beta1D, Lorentz1D """ amplitude = Parameter('amplitude') mean = Parameter('mean') stddev = Parameter('stddev') def __init__(self, amplitude, mean, stddev, **constraints): try: param_dim = len(amplitude) except TypeError: param_dim = 1 super(Gaussian1D, self).__init__(param_dim=param_dim, amplitude=amplitude, mean=mean, stddev=stddev, **constraints) @staticmethod
[docs] def eval(x, amplitude, mean, stddev): """ Model function Gauss1D """ return amplitude * np.exp(- 0.5 * (x - mean) ** 2 / stddev ** 2)
@staticmethod
[docs] def deriv(x, amplitude, mean, stddev): """ Model function derivatives Gauss1D """ d_amplitude = np.exp(-0.5 / stddev ** 2 * (x - mean) ** 2) d_mean = amplitude * d_amplitude * (x - mean) / stddev ** 2 d_stddev = amplitude * d_amplitude * (x - mean) ** 2 / stddev ** 3 return [d_amplitude, d_mean, d_stddev]
[docs]class Gaussian2D(Parametric2DModel): """ Two dimensional Gaussian model. Parameters ---------- amplitude : float Amplitude of the gaussian x_mean : float Mean of the gaussian in x y_mean : float Mean of the gaussian in y x_stddev : float Standard deviation of the gaussian in x Either x_fwhm or x_stddev must be specified y_stddev : float Standard deviation of the gaussian in y Either y_fwhm or y_stddev must be specified theta : float Rotation angle in radians. Note: increases clockwise. cov_matrix : ndarray A 2x2 covariance matrix. If specified, overrides stddev, fwhm, and theta specification. Notes ----- Model formula: .. math:: f(x, y) = A e^{-a\\left(x - x_{0}\\right)^{2} -b\\left(x - x_{0}\\right) \\left(y - y_{0}\\right) -c\\left(y - y_{0}\\right)^{2}} Using the following definitions: .. math:: a = \\left(- \\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}} - \\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}}\\right) b = \\left(\\frac{\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{y}^{2}} - \\frac{\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{x}^{2}}\\right) c = \\left(\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}} + \\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}}\\right) See Also -------- Gaussian1D, Box2D, Beta2D """ amplitude = Parameter('amplitude') x_mean = Parameter('x_mean') y_mean = Parameter('y_mean') x_stddev = Parameter('x_stddev') y_stddev = Parameter('y_stddev') theta = Parameter('theta') def __init__(self, amplitude, x_mean, y_mean, x_stddev=None, y_stddev=None, theta=0.0, cov_matrix=None, **constraints): if y_stddev is None and cov_matrix is None: raise InputParameterError( "Either y_stddev must be specified, or a " "covariance matrix.") elif x_stddev is None and cov_matrix is None: raise InputParameterError( "Either x_stddev must be specified, or a " "covariance matrix.") elif cov_matrix is not None and (x_stddev is not None or y_stddev is not None): raise InputParameterError( "Cannot specify both cov_matrix and x/y_stddev") # Compute principle coordinate system transformation elif cov_matrix is not None: cov_matrix = np.array(cov_matrix) assert cov_matrix.shape == (2, 2), "Covariance matrix must be 2x2" eig_vals, eig_vecs = np.linalg.eig(cov_matrix) x_stddev, y_stddev = np.sqrt(eig_vals) y_vec = eig_vecs[:, 0] theta = np.arctan2(y_vec[1], y_vec[0]) super(Gaussian2D, self).__init__( amplitude=amplitude, x_mean=x_mean, y_mean=y_mean, x_stddev=x_stddev, y_stddev=y_stddev, theta=theta, **constraints) @staticmethod
[docs] def eval(x, y, amplitude, x_mean, y_mean, x_stddev, y_stddev, theta): """Two dimensional Gaussian function""" a = 0.5 * ((np.cos(theta) / x_stddev) ** 2 + (np.sin(theta) / y_stddev) ** 2) b = 0.5 * (np.cos(theta) * np.sin(theta) * (1. / x_stddev ** 2 - 1. / y_stddev ** 2)) c = 0.5 * ((np.sin(theta) / x_stddev) ** 2 + (np.cos(theta) / y_stddev) ** 2) return amplitude * np.exp(-(a * (x - x_mean) ** 2 + b * (x - x_mean) * (y - y_mean) + c * (y - y_mean) ** 2))
@staticmethod
[docs] def deriv(x, y, amplitude, x_mean, y_mean, x_stddev, y_stddev, theta): """Two dimensional Gaussian function derivative""" # Helper quantities # Derivatives are not checked yet a = 0.5 * ((np.cos(theta) / x_stddev) ** 2 + (np.sin(theta) / y_stddev) ** 2) b = 0.5 * (np.cos(theta) * np.sin(theta) * (1. / x_stddev ** 2 - 1. / y_stddev ** 2)) c = 0.5 * ((np.sin(theta) / x_stddev) ** 2 + (np.cos(theta) / y_stddev) ** 2) da_dtheta = np.sin(theta) * np.cos(theta) * (1/y_stddev**2 - 1/x_stddev**2) da_dx_stddev = -np.cos(theta)**2/x_stddev**3 da_dy_stddev = -np.sin(theta)**2/y_stddev**3 db_dtheta = 0.5*np.cos(2*theta) * (1/x_stddev**2 - 1/y_stddev**2) db_dx_stddev = -np.cos(theta)*np.sin(theta)/x_stddev**3 db_dy_stddev = np.cos(theta)*np.sin(theta)/y_stddev**3 dc_dtheta = np.cos(theta)*np.sin(theta)*(1/x_stddev**2 - 1/y_stddev**2) dc_dx_stddev = -np.sin(theta)**2/x_stddev**3 dc_dy_stddev = -np.cos(theta)**2/y_stddev**3 d_A = np.exp(- a * (x - x_mean) ** 2 - b * (x - x_mean) * (y - y_mean) - c * (y - y_mean) ** 2) d_theta = amplitude * (-(x - x_mean) ** 2 * da_dtheta - (x - x_mean) * (y - y_mean) * db_dtheta - (y - y_mean) ** 2 * dc_dtheta) * d_A d_x_stddev = amplitude * (-(x - x_mean) ** 2 *da_dx_stddev - (x - x_mean) * (y - y_mean) * db_dx_stddev - (y - y_mean) ** 2 * dc_dx_stddev) * d_A d_y_stddev = amplitude * (-(x - x_mean) ** 2 *da_dy_stddev - (x - x_mean) * (y - y_mean) * db_dy_stddev - (y - y_mean) ** 2 * dc_dy_stddev) * d_A d_y_mean = amplitude * (+(x - x_mean) * b + 2 * (y - y_mean) * c) * d_A d_x_mean = amplitude * ((y - y_mean) * b + 2 * (x - x_mean) * a) * d_A return [d_A, d_x_mean, d_y_mean, d_x_stddev, d_y_stddev, d_theta]
[docs]class Shift(Model): """ Shift a coordinate. Parameters ---------- offsets : float or a list of floats offsets to be applied to a coordinate if a list - each value in the list is an offset to be applied to a column in the input coordinate array """ offsets = Parameter('offsets') def __init__(self, offsets, param_dim=1): if not isinstance(offsets, collections.Sequence): param_dim = 1 else: param_dim = len(offsets) self._offsets = offsets super(Shift, self).__init__(param_dim=param_dim)
[docs] def inverse(self): if self.param_dim == 1: return Shift(offsets=(-1) * self._offsets) else: return Shift(offsets=[off * (-1) for off in self._offsets])
@format_input def __call__(self, x): """ Transforms data using this model. Parameters ---------- x : array like or a number input """ return self._offsets + x
[docs]class Scale(Model): """ Multiply a model by a factor. Parameters ---------- factors : float or a list of floats scale for a coordinate """ factors = Parameter('factors') def __init__(self, factors, param_dim=1): if not isinstance(factors, collections.Sequence): param_dim = 1 else: param_dim = len(factors) self._factors = factors super(Scale, self).__init__(param_dim=param_dim)
[docs] def inverse(self): if self.param_dim == 1: return Scale(factors=1. / self._factors) else: return Scale(factors=[1 / factor for factor in self._factors])
@format_input def __call__(self, x): """ Transforms data using this model. Parameters ---------- x : array like or a number input """ return self._factors * x
[docs]class Sine1D(Parametric1DModel): """ One dimensional Sine model. Parameters ---------- amplitude : float Oscillation amplitude frequency : float Oscillation frequency See Also -------- Const1D, Linear1D Notes ----- Model formula: .. math:: f(x) = A \\sin(2 \\pi f x) """ amplitude = Parameter('amplitude') frequency = Parameter('frequency') def __init__(self, amplitude, frequency, **constraints): super(Sine1D, self).__init__(amplitude=amplitude, frequency=frequency, **constraints) @staticmethod
[docs] def eval(x, amplitude, frequency): """One dimensional Sine model function""" return amplitude * np.sin(2 * np.pi * frequency * x)
@staticmethod
[docs] def deriv(x, amplitude, frequency): """One dimensional Sine model derivative""" d_amplitude = np.sin(2 * np.pi * frequency * x) d_frequency = (2 * np.pi * x * amplitude * np.cos(2 * np.pi * frequency * x)) return [d_amplitude, d_frequency]
[docs]class Linear1D(Parametric1DModel): """ One dimensional Line model. Parameters ---------- slope : float Slope of the straight line intercept : float Intercept of the straight line See Also -------- Const1D Notes ----- Model formula: .. math:: f(x) = a x + b """ slope = Parameter('slope') intercept = Parameter('intercept') linear = True def __init__(self, slope, intercept, **constraints): super(Linear1D, self).__init__(slope=slope, intercept=intercept, **constraints) @staticmethod
[docs] def eval(x, slope, intercept): """One dimensional Line model function""" return slope * x + intercept
@staticmethod
[docs] def deriv(x, slope, intercept): """One dimensional Line model derivative""" d_slope = x d_intercept = np.ones_like(x) return [d_slope, d_intercept]
[docs]class Lorentz1D(Parametric1DModel): """ One dimensional Lorentzian model. Parameters ---------- amplitude : float Peak value x_0 : float Position of the peak fwhm : float Full width at half maximum See Also -------- Gaussian1D, Box1D, MexicanHat1D Notes ----- Model formula: .. math:: f(x) = \\frac{A \\gamma^{2}}{\\gamma^{2} + \\left(x - x_{0}\\right)^{2}} """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') fwhm = Parameter('fwhm') def __init__(self, amplitude, x_0, fwhm, **constraints): super(Lorentz1D, self).__init__(amplitude=amplitude, x_0=x_0, fwhm=fwhm, **constraints) @staticmethod
[docs] def eval(x, amplitude, x_0, fwhm): """One dimensional Lorentzian model function""" return (amplitude * ((fwhm / 2.) ** 2) / ((x - x_0) ** 2 + (fwhm / 2.) ** 2))
@staticmethod
[docs] def deriv(x, amplitude, x_0, fwhm): """One dimensional Lorentzian model derivative""" d_amplitude = fwhm ** 2 / (fwhm ** 2 + (x - x_0) ** 2) d_x_0 = (amplitude * d_amplitude * (2 * x - 2 * x_0) / (fwhm ** 2 + (x - x_0) ** 2)) d_fwhm = 2 * amplitude * d_amplitude / fwhm * (1 - d_amplitude) return [d_amplitude, d_x_0, d_fwhm]
[docs]class Const1D(Parametric1DModel): """ One dimensional Constant model. Parameters ---------- amplitude : float Value of the constant function See Also -------- Const2D Notes ----- Model formula: .. math:: f(x) = A """ amplitude = Parameter('amplitude') def __init__(self, amplitude, **constraints): super(Const1D, self).__init__(amplitude=amplitude, **constraints) @staticmethod
[docs] def eval(x, amplitude): """One dimensional Constant model function""" return amplitude * np.ones_like(x)
@staticmethod
[docs] def deriv(x, amplitude): """One dimensional Constant model derivative""" d_amplitude = np.ones_like(x) return [d_amplitude]
[docs]class Const2D(Parametric2DModel): """ Two dimensional Constant model. Parameters ---------- amplitude : float Value of the constant function See Also -------- Const1D Notes ----- Model formula: .. math:: f(x, y) = A """ amplitude = Parameter('amplitude') def __init__(self, amplitude, **constraints): super(Const2D, self).__init__(amplitude=amplitude, **constraints) @staticmethod
[docs] def eval(x, y, amplitude): """Two dimensional Constant model function""" return amplitude * np.ones_like(x)
[docs]class Disk2D(Parametric2DModel): """ Two dimensional radial symmetric Disk model. Parameters ---------- amplitude : float Value of the disk function x_0 : float x position center of the disk y_0 : float y position center of the disk R_0 : float Radius of the disk See Also -------- Box2D, TrapezoidDisk2D Notes ----- Model formula: .. math:: f(r) = \\left \\{ \\begin{array}{ll} A & : r \\leq R_0 \\\\ 0 & : r > R_0 \\end{array} \\right. """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') y_0 = Parameter('y_0') R_0 = Parameter('R_0') def __init__(self, amplitude, x_0, y_0, R_0, **constraints): super(Disk2D, self).__init__(amplitude=amplitude, x_0=x_0, y_0=y_0, R_0=R_0, **constraints) @staticmethod
[docs] def eval(x, y, amplitude, x_0, y_0, R_0): """Two dimensional Disk model function""" rr = (x - x_0) ** 2 + (y - y_0) ** 2 return np.select([rr <= R_0 ** 2], [amplitude])
[docs]class Ring2D(Parametric2DModel): """ Two dimensional radial symmetric Ring model. Parameters ---------- amplitude : float Value of the disk function x_0 : float x position center of the disk y_0 : float y position center of the disk r_in : float Inner radius of the ring width : float Width of the ring. r_out : float Outer Radius of the ring. Can be specified instead of width. See Also -------- Disk2D, TrapezoidDisk2D Notes ----- Model formula: .. math:: f(r) = \\left \\{ \\begin{array}{ll} A & : r_{in} \\leq r \\leq r_{out} \\\\ 0 & : \\textnormal{else} \\end{array} \\right. Where :math:`r_{out} = r_{in} + r_{width}`. """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') y_0 = Parameter('y_0') r_in = Parameter('r_in') width = Parameter('width') def __init__(self, amplitude, x_0, y_0, r_in, width=None, r_out=None, **constraints): if r_out is not None: width = r_out - r_in if r_out is None and width is None: raise ModelDefinitionError("Either specify width or r_out.") super(Ring2D, self).__init__(amplitude=amplitude, x_0=x_0, y_0=y_0, r_in=r_in, width=width, **constraints) @staticmethod
[docs] def eval(x, y, amplitude, x_0, y_0, r_in, width): """Two dimensional Ring model function.""" rr = (x - x_0) ** 2 + (y - y_0) ** 2 r_range = np.logical_and(rr >= r_in ** 2, rr <= (r_in + width) ** 2) return np.select([r_range], [amplitude])
class Delta1D(Parametric1DModel): """One dimensional Dirac delta function.""" def __init__(self): raise ModelDefinitionError("Not implemented") class Delta2D(Parametric2DModel): """Two dimensional Dirac delta function.""" def __init__(self): raise ModelDefinitionError("Not implemented")
[docs]class Box1D(Parametric1DModel): """ One dimensional Box model. Parameters ---------- amplitude : float Amplitude A x_0 : float Position of the center of the box function width : float Width of the box See Also -------- Box2D, TrapezoidDisk2D Notes ----- Model formula: .. math:: f(x) = \\left \\{ \\begin{array}{ll} A & : x_0 - w/2 \\geq x \\geq x_0 + w/2 \\\\ 0 & : \\textnormal{else} \\end{array} \\right. """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') width = Parameter('width') def __init__(self, amplitude, x_0, width, **constraints): super(Box1D, self).__init__(amplitude=amplitude, x_0=x_0, width=width, **constraints) @staticmethod
[docs] def eval(x, amplitude, x_0, width): """One dimensional Box model function""" return np.select([np.logical_and(x >= x_0 - width / 2., x <= x_0 + width / 2.)], [amplitude], 0)
@classmethod
[docs] def deriv(cls, x, amplitude, x_0, width): """One dimensional Box model derivative""" d_amplitude = cls.eval(x, 1, x_0, width) d_x_0 = np.zeros_like(x) d_width = np.zeros_like(x) return [d_amplitude, d_x_0, d_width]
[docs]class Box2D(Parametric2DModel): """ Two dimensional Box model. Parameters ---------- amplitude : float Amplitude A x_0 : float x position of the center of the box function x_width : float Width in x direction of the box y_0 : float y position of the center of the box function y_width : float Width in y direction of the box See Also -------- Box1D, Gaussian2D, Beta2D Notes ----- Model formula: .. math:: f(x, y) = \\left \\{ \\begin{array}{ll} A & : x_0 - w_x/2 \\geq x \\geq x_0 + w_x/2 \\\\ A & : y_0 - w_y/2 \\geq y \\geq y_0 + w_y/2 \\\\ 0 & : \\textnormal{else} \\end{array} \\right. """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') y_0 = Parameter('y_0') x_width = Parameter('x_width') y_width = Parameter('y_width') def __init__(self, amplitude, x_0, y_0, x_width, y_width, **constraints): super(Box2D, self).__init__(amplitude=amplitude, x_0=x_0, y_0=y_0, x_width=x_width, y_width=y_width, **constraints) @staticmethod
[docs] def eval(x, y, amplitude, x_0, y_0, x_width, y_width): """Two dimensional Box model function""" x_range = np.logical_and(x >= x_0 - x_width / 2., x <= x_0 + x_width / 2.) y_range = np.logical_and(y >= y_0 - y_width / 2., y <= y_0 + y_width / 2.) return np.select([np.logical_and(x_range, y_range)], [amplitude], 0)
[docs]class Trapezoid1D(Parametric1DModel): """ One dimensional Trapezoid model. Parameters ---------- amplitude : float Amplitude of the trapezoid x_0 : float Center position of the trapezoid width : float Width of the constant part of the trapezoid. slope : float Slope of the tails of the trapezoid See Also -------- Box1D, Gaussian1D, Beta1D """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') width = Parameter('width') slope = Parameter('slope') def __init__(self, amplitude, x_0, width, slope, **constraints): super(Trapezoid1D, self).__init__(amplitude=amplitude, x_0=x_0, width=width, slope=slope, **constraints) @staticmethod
[docs] def eval(x, amplitude, x_0, width, slope): """One dimensional Trapezoid model function""" # Compute the four points where the trapezoid changes slope # x1 <= x2 <= x3 <= x4 x2 = x_0 - width / 2. x3 = x_0 + width / 2. x1 = x2 - amplitude / slope x4 = x3 + amplitude / slope # Compute model values in pieces between the change points range_a = np.logical_and(x >= x1, x < x2) range_b = np.logical_and(x >= x2, x < x3) range_c = np.logical_and(x >= x3, x < x4) val_a = slope * (x - x1) val_b = amplitude val_c = slope * (x4 - x) return np.select([range_a, range_b, range_c], [val_a, val_b, val_c])
[docs]class TrapezoidDisk2D(Parametric2DModel): """ Two dimensional circular Trapezoid model. Parameters ---------- amplitude : float Amplitude of the trapezoid x_0 : float x position of the center of the trapezoid y_0 : float y position of the center of the trapezoid R_0 : float Radius of the constant part of the trapezoid. slope : float Slope of the tails of the trapezoid in x direction. See Also -------- Disk2D, Box2D """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') y_0 = Parameter('y_0') R_0 = Parameter('R_0') slope = Parameter('slope') def __init__(self, amplitude, x_0, y_0, R_0, slope, **constraints): super(TrapezoidDisk2D, self).__init__(amplitude=amplitude, x_0=x_0, y_0=y_0, R_0=R_0, slope=slope, **constraints) @staticmethod
[docs] def eval(x, y, amplitude, x_0, y_0, R_0, slope): """Two dimensional Trapezoid Disk model function""" r = np.sqrt((x - x_0) ** 2 + (y - y_0) ** 2) range_1 = r <= R_0 range_2 = np.logical_and(r > R_0, r <= R_0 + amplitude / slope) val_1 = amplitude val_2 = amplitude + slope * (R_0 - r) return np.select([range_1, range_2], [val_1, val_2])
[docs]class MexicanHat1D(Parametric1DModel): """ One dimensional Mexican Hat model. Parameters ---------- amplitude : float Amplitude x_0 : float Position of the peak sigma : float Width of the Mexican hat See Also -------- MexicanHat2D, Box1D, Gaussian1D, Trapezoid1D Notes ----- Model formula: .. math:: f(x) = {A \\left(1 - \\frac{\\left(x - x_{0}\\right)^{2}}{\\sigma^{2}}\\right) e^{- \\frac{\\left(x - x_{0}\\right)^{2}}{2 \\sigma^{2}}}} """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') sigma = Parameter('sigma') def __init__(self, amplitude, x_0, sigma, **constraints): super(MexicanHat1D, self).__init__(amplitude=amplitude, x_0=x_0, sigma=sigma, **constraints) @staticmethod
[docs] def eval(x, amplitude, x_0, sigma): """One dimensional Mexican Hat model function""" xx_ww = (x - x_0) ** 2 / (2 * sigma ** 2) return amplitude * (1 - 2 * xx_ww) * np.exp(-xx_ww)
[docs]class MexicanHat2D(Parametric2DModel): """ Two dimensional symmetric Mexican Hat model. Parameters ---------- amplitude : float Amplitude x_0 : float x position of the peak y_0 : float y position of the peak sigma : float Width of the Mexican hat See Also -------- MexicanHat1D, Gaussian2D Notes ----- Model formula: .. math:: f(x, y) = A \\left(1 - \\frac{\\left(x - x_{0}\\right)^{2} + \\left(y - y_{0}\\right)^{2}}{\\sigma^{2}}\\right) e^{\\frac{- \\left(x - x_{0}\\right)^{2} - \\left(y - y_{0}\\right)^{2}}{2 \\sigma^{2}}} """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') y_0 = Parameter('y_0') sigma = Parameter('sigma') def __init__(self, amplitude, x_0, y_0, sigma, **constraints): super(MexicanHat2D, self).__init__(amplitude=amplitude, x_0=x_0, y_0=y_0, sigma=sigma, **constraints) @staticmethod
[docs] def eval(x, y, amplitude, x_0, y_0, sigma): """Two dimensional Mexican Hat model function""" rr_ww = ((x - x_0) ** 2 + (y - y_0) ** 2) / (2 * sigma ** 2) return amplitude * (1 - rr_ww) * np.exp(- rr_ww)
[docs]class AiryDisk2D(Parametric2DModel): """ Two dimensional Airy disk model. Parameters ---------- amplitude : float Amplitude of the Airy function. x_0 : float x position of the maximum of the Airy function. y_0 : float y position of the maximum of the Airy function. width : float Width of the Airy function. See Also -------- Box2D, TrapezoidDisk2D, Gaussian2D Notes ----- Model formula: .. math:: f(r) = A \\frac{J_1(2 \\pi r)}{\\pi r} Where J1 is the first order Bessel function of first kind. """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') y_0 = Parameter('y_0') width = Parameter('width') _j1 = None def __init__(self, amplitude, x_0, y_0, width, **constraints): if self._j1 is None: try: from scipy.special import j1 self.__class__._j1 = j1 # add a ValueError here for python3 + scipy < 0.12 except (ValueError, ImportError): raise ImportError("AiryDisk2D model requires scipy > 0.11.") super(AiryDisk2D, self).__init__(amplitude=amplitude, x_0=x_0, y_0=y_0, width=width, **constraints) def __deepcopy__(self, memo): new_model = self.__class__(self.amplitude.value, self.x_0.value, self.y_0.value, self.width.value) return new_model def __copy__(self): new_model = self.__class__(self.amplitude.value, self.x_0.value, self.y_0.value, self.width.value) return new_model @classmethod
[docs] def eval(cls, x, y, amplitude, x_0, y_0, width): """Two dimensional Airy model function""" r = np.sqrt((x - x_0) ** 2 + (y - y_0) ** 2) / width # Since r can be zero, we have to take care to treat that case # separately so as not to raise a Numpy warning z = np.ones(r.shape) z[r > 0] = (amplitude * (cls._j1(2 * np.pi * r[r > 0]) / (np.pi * r[r > 0])) ** 2) return z
[docs]class Beta1D(Parametric1DModel): """ One dimensional Beta model. Parameters ---------- amplitude : float Amplitude of the model. x_0 : float x position of the maximum of the Beta model. gamma : float Core width of the Beta model. alpha : float Power index of the beta model. See Also -------- Gaussian1D, Box1D Notes ----- Model formula: .. math:: f(x) = A \\left(1 + \\frac{\\left(x - x_{0}\\right)^{2}}{\\gamma^{2}}\\right)^{- \\alpha} """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') gamma = Parameter('gamma') alpha = Parameter('alpha') def __init__(self, amplitude, x_0, gamma, alpha, **constraints): super(Beta1D, self).__init__(amplitude=amplitude, x_0=x_0, gamma=gamma, alpha=alpha, **constraints) @staticmethod
[docs] def eval(x, amplitude, x_0, gamma, alpha): """One dimensional Beta model function""" return amplitude * (1 + ((x - x_0) / gamma) ** 2) ** (-alpha)
@staticmethod
[docs] def deriv(x, amplitude, x_0, gamma, alpha): """One dimensional Beta model derivative""" d_A = (1 + (x - x_0) ** 2 / gamma ** 2) ** (-alpha) d_x_0 = (-amplitude * alpha * d_A * (-2 * x + 2 * x_0) / (gamma ** 2 * d_A ** alpha)) d_gamma = (2 * amplitude * alpha * d_A * (x - x_0) ** 2 / (gamma ** 3 * d_A ** alpha)) d_alpha = -amplitude * d_A * np.log(1 + (x - x_0) ** 2 / gamma ** 2) return [d_A, d_x_0, d_gamma, d_alpha]
[docs]class Beta2D(Parametric2DModel): """ Two dimensional Beta model. Parameters ---------- amplitude : float Amplitude of the model. x_0 : float x position of the maximum of the Beta model. y_0 : float y position of the maximum of the Beta model. gamma : float Core width of the Beta model. alpha : float Power index of the beta model. See Also -------- Gaussian2D, Box2D Notes ----- Model formula: .. math:: f(x, y) = A \\left(1 + \\frac{\\left(x - x_{0}\\right)^{2} + \\left(y - y_{0}\\right)^{2}}{\\gamma^{2}}\\right)^{- \\alpha} """ amplitude = Parameter('amplitude') x_0 = Parameter('x_0') y_0 = Parameter('y_0') gamma = Parameter('gamma') alpha = Parameter('alpha') def __init__(self, amplitude, x_0, y_0, gamma, alpha, **constraints): super(Beta2D, self).__init__(amplitude=amplitude, x_0=x_0, y_0=y_0, gamma=gamma, alpha=alpha, **constraints) @staticmethod
[docs] def eval(x, y, amplitude, x_0, y_0, gamma, alpha): """Two dimensional Beta model function""" rr_gg = ((x - x_0) ** 2 + (y - y_0) ** 2) / gamma ** 2 return amplitude * (1 + rr_gg) ** (-alpha)
@staticmethod
[docs] def deriv(x, y, amplitude, x_0, y_0, gamma, alpha): """Two dimensional Beta model derivative""" rr_gg = ((x - x_0) ** 2 + (y - y_0) ** 2) / gamma ** 2 d_A = (1 + rr_gg) ** (-alpha) d_x_0 = (-amplitude * alpha * d_A * (-2 * x + 2 * x_0) / (gamma ** 2 * (1 + rr_gg))) d_y_0 = (-amplitude * alpha * d_A * (-2 * y + 2 * y_0) / (gamma ** 2 * (1 + rr_gg))) d_alpha = -amplitude * d_A * np.log(1 + rr_gg) d_gamma = 2 * amplitude * alpha * d_A * (rr_gg / (gamma * (1 + rr_gg))) return [d_A, d_x_0, d_y_0, d_gamma, d_alpha]
[docs]def custom_model_1d(func, func_deriv=None): """ Create a one dimensional model from a user defined function. The parameters of the model will be inferred from the arguments of the function. ..note :: All model parameters have to be defined as keyword arguments with default values in the model function. If you want to use parameter sets in the model, the parameters should be treated as lists or arrays. Parameters ---------- func : function Function which defines the model. It should take one positional argument (the independent variable in the model), and any number of keyword arguments (the parameters). It must return the value of the model (typically as an array, but can also be a scalar for scalar inputs). This corresponds to the `ParametricModel.eval` method. func_deriv : function, optional Function which defines the Jacobian derivative of the model. I.e., the derivive with respect to the *parameters* of the model. It should have the same argument signature as `func`, but should return a sequence where each element of the sequence is the derivative with respect to the correseponding argument. This corresponds to the `ParametricModel.deriv` method. Examples -------- Define a sinusoidal model function as a custom 1D model: >>> from astropy.modeling.models import custom_model_1d >>> import numpy as np >>> def sine_model(x, amplitude=1., frequency=1.): ... return amplitude * np.sin(2 * np.pi * frequency * x) >>> def sine_deriv(x, amplitude=1., frequency=1.): ... return 2 * np.pi * amplitude * np.cos(2 * np.pi * frequency * x) >>> SineModel = custom_model_1d(sine_model, func_deriv=sine_deriv) Create an instance of the custom model and evaluate it: >>> model = SineModel() >>> model(0.25) 1.0 This model instance can now be used like a usual astropy model. """ if not callable(func): raise ModelDefinitionError("Not callable. Must be function") if func_deriv is not None and not callable(func_deriv): raise ModelDefinitionError("func_deriv not callable. Must be function") model_name = func.__name__ param_values = func.func_defaults # Check if all parameters are keyword arguments nparams = len(param_values) if func_deriv is not None and len(func_deriv.func_defaults) != nparams: raise ModelDefinitionError("derivative function should accept" " same number of parameters as func.") if func.func_code.co_argcount == nparams + 1: param_names = func.func_code.co_varnames[1:nparams + 1] else: raise ModelDefinitionError( "All parameters must be keyword arguments") params = dict((name, Parameter(name, default=default)) for name, default in zip(param_names, param_values)) arg_signature_1 = ', '.join('{0}=None'.format(name) for name in param_names) arg_signature_2 = ', '.join('{0}={0}'.format(name) for name in param_names) mod = find_current_module(2) if mod: filename = mod.__file__ modname = mod.__name__ else: filename = '<string>' modname = '__main__' members = {'eval': staticmethod(func)} eval_globals = {} init_code_string = dedent(""" def __init__(self, {0}, **constraints): super(self.__class__, self).__init__({1}, **constraints) """).format(arg_signature_1, arg_signature_2) eval(compile(init_code_string, filename, 'single'), eval_globals) if func_deriv is not None: members['deriv'] = staticmethod(func_deriv) members['__init__'] = eval_globals['__init__'] members.update(params) cls = type(model_name, (Parametric1DModel,), members) cls.__module__ = modname return cls

Page Contents