# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-
"""
Implements a set of line profile parameterizations.
----
.. include license and copyright
.. include:: ../include/copy.rst
----
.. include common links, assuming primary doc root is up one directory
.. include:: ../include/links.rst
"""
import numpy
from scipy import special, fftpack
from astropy.modeling import FittableModel, Parameter
#-----------------------------------------------------------------------
[docs]
class GaussianLineProfile(FittableModel):
r"""
Define a Gaussian line profile as parameterized by its zeroth
moment, mean, and standard deviation:
.. math::
\mathcal{N}(x|f,\mu,\sigma) = \frac{f}{\sqrt{2\pi}\sigma}
\exp\left(\frac{-(x-\mu)^2}{2\sigma^2}\right)
The base class is `astropy.modeling.FittableModel`_, which
facilitates its use in combining multiple components and other
models in the `astropy.modeling`_ suite.
"""
# Inputs the x positions
inputs = ('x',)
# Returns the renormalized Gaussian PDF at each x position
outputs = ('y',)
# Parameters are the peak value of the profile, the profile mean,
# and the profile standard deviation (sigma)
zmom=Parameter(default=1.0)
mean=Parameter(default=0.0)
sigma=Parameter(default=1.0)
[docs]
@staticmethod
def evaluate(x, zmom, mean, sigma):
# Discrete samples
return zmom * numpy.exp(-0.5*numpy.square((x - mean)/sigma)) / numpy.sqrt(2*numpy.pi)/sigma
# # Integrated over a pixel
# diff = (x-mean)/pixelscale
# denom = numpy.sqrt(2.0)*sigma/pixelscale
# return zmom * (erf((diff+0.5)/denom) - erf((diff-0.5)/denom))/2.
[docs]
@staticmethod
def fit_deriv(x, zmom, mean, sigma):
y = GaussianLineProfile.evaluate(x,zmom,mean,sigma)
d_zmom = y/zmom
d_mean = y*(x-mean)/numpy.square(sigma)
d_sigma = y*(numpy.square(x-mean)/numpy.power(sigma,3) - 1.0/sigma)
return [d_zmom, d_mean, d_sigma]
[docs]
@staticmethod
def flux(zmom, mean, sigma):
return zmom
[docs]
@staticmethod
def flux_err(zmom, mean, sigma, zmome, meane, sigmae):
return zmome
[docs]
@staticmethod
def moment(order, zmom, mean, sigma):
if order == 0:
return zmom
if order == 1:
return mean
if order == 2:
return sigma
[docs]
@staticmethod
def moment_err(order, zmom, mean, sigma, zmome, meane, sigmae):
if order == 0:
return zmome
if order == 1:
return meane
if order == 2:
return sigmae
[docs]
@staticmethod
def integral(zmom, mean, sigma):
return integrate.quad( GaussianLineProfile.evaluate, -sigma, sigma, args=(zmom,mean,sigma))
[docs]
@staticmethod
def scale_flux(zmom, mean, sigma, fac):
return [zmom*fac, mean, sigma]
[docs]
@staticmethod
def fix_flux():
return [ True, False, False ]
[docs]
@staticmethod
def mean_indx():
return 1
[docs]
@staticmethod
def shift_mean(zmom, mean, sigma, shift):
return [zmom, mean+shift, sigma ]
[docs]
@staticmethod
def fix_mean():
return [ False, True, False ]
[docs]
@staticmethod
def scale_stddev(zmom, mean, sigma, fac):
return [zmom, mean, sigma*fac ]
[docs]
@staticmethod
def fix_stddev():
return [ False, False, True ]
[docs]
class GaussianLSF:
r"""
Define a Gaussian line profile, *sampled* over the width of the
sampling step, parameterized by its integral (:math:`F`), center
(:math:`\mu`), and standard deviation (:math:`\sigma`). I.e:
.. math::
\mathcal{N}(x|f,\mu,\sigma) = \frac{f}{\sqrt{2\pi}\sigma}
\exp\left(\frac{-\Delta^2}{2\sigma^2}\right)
where :math:`\Delta = x-\mu`. The coordinate vector :math:`x` does
not need to be uniformly sampled.
Args:
p (array-like): (**Optional**) Input parameters ordered as the
total integral of the profile, the profile center, and the
profile standard deviation. Assumed to be (1.0, 0.0, 1.0)
by default.
Attributes:
p (numpy.ndarray): Most recently used parameters
Raises:
ValueError: Raised if the provided parameter vector is not 3
elements long.
"""
def __init__(self, p=None):
self.set_par(p)
def __call__(self, x, p):
"""Calculate the profile.
Args:
x (array-like): Independent variable.
p (array-like): LSF parameters.
"""
self.set_par(p)
return self.sample(x)
[docs]
@staticmethod
def npar():
return 3
[docs]
def set_par(self, p):
"""
Set the internal parameters to the provided set.
Args:
p (array-like): LSF parameters.
Raises:
ValueError: Raised if the provided parameter vector is not 3
elements long.
"""
if p is None:
self.p = numpy.array([1.0, 0.0, 1.0])
return
if len(p) != GaussianLSF.npar():
raise ValueError('Must provide {0} parameters.'.format(GaussianLSF.npar()))
self.p = numpy.asarray(p)
[docs]
def sample(self, x):
"""
Sample the profile.
Args:
x (array-like): Independent variable.
"""
return self.p[0] * numpy.exp(-numpy.square((x-self.p[1])/self.p[2])/2.) \
/ numpy.sqrt(2.0*numpy.pi) / self.p[2]
[docs]
def parameters_from_moments(self, mom0, mom1, mom2):
"""
Provided the 0th, 1st, and 2nd moments, produce a set of
parameters for the profile.
"""
return numpy.array([mom0, mom1, mom2])
[docs]
class IntegratedGaussianLSF(GaussianLSF):
r"""
Define a Gaussian line profile, *integrated* over the width of the
sampling step, parameterized by its integral (:math:`F`), center
(:math:`\mu`), and standard deviation (:math:`\sigma`). I.e:
.. math::
\mathcal{N}(x|F,\mu,\sigma) = \frac{F}{2} \left[
{\rm erf}\left(\frac{\Delta+\delta_x/2}{\sqrt{2}\sigma}\right) -
{\rm erf}\left(\frac{\Delta-\delta_x/2}{\sqrt{2}\sigma}\right)\right]
where :math:`{\rm erf}(x)` is the error function, :math:`\Delta =
x-\mu`, and :math:`\delta_x` is the sampling step. The sampling
*must* be uniform in :math:`x`.
Args:
p (array-like): (**Optional**) Input parameters ordered as the
total integral of the profile, the profile center, and the
profile standard deviation. Assumed to be (1.0, 0.0, 1.0)
by default.
dx (float): (**Optional**) Sampling width. Default is 1.
Attributes:
p (numpy.ndarray): Most recently used parameters
dx (float): Assumed sampling.
Raises:
ValueError: Raised if the provided parameter vector is not 3
elements long.
"""
def __init__(self, p=None, dx=None):
self.set_par(p)
self.dx = 1.0 if dx is None else dx
[docs]
def sample(self, x):
"""
Sample the profile.
.. warning::
Does **not** check if the provided :math:`x` values are
sampled at :attr:`dx`.
Args:
x (array-like): Independent variable.
"""
n = numpy.sqrt(2.)*self.p[2]
d = numpy.asarray(x)-self.p[1]
return self.p[0] * (special.erf((d+self.dx/2.)/n) - special.erf((d-self.dx/2.)/n))/2.
[docs]
class FFTGaussianLSF(GaussianLSF):
r"""
Define a Gaussian line profile by first constructing the analytic
FFT of the profile and then returning the inverse real FFT. See
ppxf_util.emline by M. Cappellari. The sampling *must* be uniform
in :math:`x`.
Args:
p (array-like): (**Optional**) Input parameters ordered as the
total integral of the profile, the profile center, and the
profile standard deviation. Assumed to be (1.0, 0.0, 1.0)
by default.
dx (float): (**Optional**) Sampling width. Default is 1.
pixel (bool): (**Optional**) Flag to produce profile integrated
over the sampling width.
Attributes:
p (numpy.ndarray): Most recently used parameters
dx (float): Assumed sampling.
pixel (bool): Flag to produce profile integrated over the
sampling width.
Raises:
ValueError: Raised if the provided parameter vector is not 3
elements long.
"""
def __init__(self, p=None, dx=None, pixel=True):
self.set_par(p)
self.dx = 1.0 if dx is None else dx
self.pixel = pixel
[docs]
def sample(self, x):
"""
Sample the profile.
.. warning::
Does **not** check if the provided :math:`x` values are
sampled at :attr:`dx`.
Args:
x (array-like): Independent variable.
"""
# Require that the center of the line be within the range of x
if self.p[1] < x[0] or self.p[1] > x[-1]:
return numpy.zeros_like(x)
xsig = self.p[2]/self.dx
x0 = (self.p[1]-x[0])/self.dx
npad = 2**int(numpy.ceil(numpy.log2(x.size)))
# npad = fftpack.next_fast_len(x.size)
w = numpy.linspace(0,numpy.pi,npad//2+1)
rfft = self.p[0]*numpy.exp(-0.5*numpy.square(w*xsig) - 1j*w*x0)
if self.pixel:
rfft *= numpy.sinc(w/(2*numpy.pi))
lsf = numpy.fft.irfft(rfft, n=npad)[:x.size]
return lsf if self.pixel else lsf/self.dx
#-----------------------------------------------------------------------
[docs]
class NCompLineProfile:
"""
Construct a single line profile from many components with the same
profile parameterization.
"""
def __init__(self, ncomp, par=None, err=None, profile=GaussianLineProfile):
# Check the type of the profile to use
if not issubclass(profile, FittableModel):
raise TypeError('Profile must be a subclass of FittableModel.')
self.profile_class = profile
# Ensure that the profile can be defined
if not ncomp > 0:
raise ValueError('Number of input components must be 1 or higher!')
# Calculate and save the number of parameters
self.ncomp = ncomp
self.npar_per_component = len(self.profile_class.param_names)
self.npar = self.ncomp*self.npar_per_component
# Check the parameter vector, if provided
if par is not None and len(par.shape) > 2:
raise ValueError('Input parameter array should only be one or two dimensional.')
_par = None if par is None else par.ravel()
if _par is not None and len(_par) != self.npar:
raise ValueError('Incorrect number of parameters provided')
self.par = None if _par is None else _par.reshape(self.ncomp,-1).astype(float)
# Check the parameter vector, if provided
if err is not None and len(err.shape) > 2:
raise ValueError('Input parameter error array should only be one or two dimensional.')
_err = None if err is None else err.ravel()
if _err is not None and len(_err) != self.npar:
raise ValueError('Incorrect number of parameter errors provided')
self.err = None if _err is None else _err.reshape(self.ncomp,-1).astype(float)
# Build the profile function
self.profile = self.profile_class() \
if self.par is None else self.profile_class(*self.par[0])
for i in range(1,ncomp):
self.profile += self.profile_class() \
if self.par is None else self.profile_class(*self.par[i])
def __call__(self, x, par=None):
return self.sample(x,par)
[docs]
def _quick_sample(self, x):
"""
Sample without providing/checking new input parameters.
"""
return self.profile.evaluate(x, *self.par.ravel())
[docs]
def _mom1_integrand(self, x):
return x*self._quick_sample(x)
[docs]
def _mom2_integrand(self, x):
return x*x*self._quick_sample(x)
[docs]
def assign_par(self, par):
_par = par.ravel() if len(par.shape) == 2 else par
if len(_par) != self.npar:
raise ValueError('Incorrect number of parameters provided')
self.par = _par.copy().reshape(self.ncomp,-1)
[docs]
def assign_err(self, err):
_err = err.ravel() if len(err.shape) == 2 else err
if len(_err) != self.npar:
raise ValueError('Incorrect number of parameters provided')
self.err = _err.copy().reshape(self.ncomp,-1)
[docs]
def sample(self, x, par=None):
if par is not None:
self.assign_par(par)
return self._quick_sample(x)
[docs]
def flux(self, par=None):
if par is not None:
self.assign_par(par)
return numpy.sum( [ self.profile_class.flux(*p) \
for p in numpy.take(self.par,numpy.arange(self.ncomp),axis=0) ] )
[docs]
def flux_err(self, par=None, err=None):
if par is not None:
self.assign_par(par)
if err is not None:
self.assign_err(err)
return numpy.sum([ self.profile_class.flux_err(*p, *e) \
for p,e in zip(numpy.take(self.par,numpy.arange(self.ncomp),axis=0),
numpy.take(self.err,numpy.arange(self.ncomp),axis=0)) ] )
[docs]
def moment(self, order=0, par=None):
"""
.. todo::
impose some reasonable limits for the integrals
"""
if par is not None:
self.assign_par(par)
if self.ncomp == 1:
return self.profile_class.moment(order, *self.par[0,:])
flux = self.flux(par=par)
if order == 0:
return flux
mom1, _ = integrate.quad(self._mom1_integrand, -numpy.inf, numpy.inf)
if order == 1:
return mom1/flux
mom2, _ = integrate.quad(self._mom2_integrand, -numpy.inf, numpy.inf)
return numpy.sqrt(mom2/flux-numpy.square(mom1/flux))
[docs]
def moment_err(self, order=0, par=None, err=None):
"""
.. todo::
impose some reasonable limits for the integrals
"""
if self.err is None and err is None:
return 0.0
if par is not None:
self.assign_par(par)
if err is not None:
self.assign_err(err)
if self.ncomp == 1:
return self.profile_class.moment_err(order, *self.par[0,:], *self.err[0,:])
raise NotImplementedError('Unable to calculate error for multiple components.')
[docs]
def scale_flux(self, fac):
for i in range(self.ncomp):
self.par[i,:] = self.profile_class.scale_flux(*self.par[i,:], fac)
[docs]
def set_flux(self, flx):
self.scale_flux( flx/self.flux() )
[docs]
def fix_flux(self):
return numpy.array(self.profile_class.fix_flux()*self.ncomp)
[docs]
def mean_indx(self):
return [ self.profile_class.mean_indx()+i*self.ncomp for i in range(self.ncomp) ]
[docs]
def shift_mean(self, shift):
for i in range(self.ncomp):
self.par[i,:] = self.profile_class.shift_mean(*self.par[i,:], shift)
[docs]
def set_mean(self, mean):
self.shift_mean( mean-self.moment(order=1) )
[docs]
def fix_mean(self):
return numpy.array(self.profile_class.fix_mean()*self.ncomp)
[docs]
def scale_stddev(self, fac):
if self.ncomp != 1:
raise NotImplementedError('Cannot scale the standard deviation of multi-component profiles.')
self.par[0,:] = self.profile_class.scale_stddev(*self.par[0,:], fac)
[docs]
def set_stddev(self, stddev):
self.scale_stddev( stddev/self.moment(order=2) )
[docs]
def fix_stddev(self):
if self.ncomp != 1:
raise NotImplementedError('Cannot fix the standard deviation of multi-component profiles.')
return numpy.array(self.profile_class.fix_stddev()*self.ncomp)