# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
__all__ = [
'Gauss2DPDF',
'MultiGauss2D',
'gaussian_sum_moments',
]
__doctest_requires__ = {('gaussian_sum_moments'): ['uncertainties']}
[docs]class Gauss2DPDF(object):
"""2D symmetric Gaussian PDF.
Reference: http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Bivariate_case
Parameters
----------
sigma : float
Gaussian width.
"""
def __init__(self, sigma=1):
self.sigma = np.asarray(sigma, np.float64)
@property
def _sigma2(self):
"""Sigma squared (float)"""
return self.sigma * self.sigma
@property
def amplitude(self):
"""PDF amplitude at the center (float)"""
return self.__call(0, 0)
[docs] def __call__(self, x, y=0):
"""dp / (dx dy) at position (x, y)
Parameters
----------
x : `~numpy.ndarray`
x coordinate
y : `~numpy.ndarray`, optional
y coordinate
Returns
-------
dpdxdy : `~numpy.ndarray`
dp / (dx dy)
"""
x = np.asarray(x, dtype=np.float64)
y = np.asarray(y, dtype=np.float64)
theta2 = x * x + y * y
amplitude = 1 / (2 * np.pi * self._sigma2)
exponent = -0.5 * theta2 / self._sigma2
return amplitude * np.exp(exponent)
[docs] def dpdtheta2(self, theta2):
"""dp / dtheta2 at position theta2 = theta ^ 2
Parameters
----------
theta2 : `~numpy.ndarray`
Offset squared
Returns
-------
dpdtheta2 : `~numpy.ndarray`
dp / dtheta2
"""
theta2 = np.asarray(theta2, dtype=np.float64)
amplitude = 1 / (2 * self._sigma2)
exponent = -0.5 * theta2 / self._sigma2
return amplitude * np.exp(exponent)
[docs] def containment_fraction(self, theta):
"""Containment fraction.
Parameters
----------
theta : `~numpy.ndarray`
Offset
Returns
-------
containment_fraction : `~numpy.ndarray`
Containment fraction
"""
theta = np.asarray(theta, dtype=np.float64)
return 1 - np.exp(-0.5 * theta ** 2 / self._sigma2)
[docs] def containment_radius(self, containment_fraction):
"""Containment angle for a given containment fraction.
Parameters
----------
containment_fraction : `~numpy.ndarray`
Containment fraction
Returns
-------
containment_radius : `~numpy.ndarray`
Containment radius
"""
containment_fraction = np.asarray(containment_fraction, dtype=np.float64)
return self.sigma * np.sqrt(-2 * np.log(1 - containment_fraction))
[docs] def gauss_convolve(self, sigma):
"""Convolve with another Gaussian 2D PDF.
Parameters
----------
sigma : `~numpy.ndarray` or float
Gaussian width of the new Gaussian 2D PDF to covolve with.
Returns
-------
gauss_convolve : `~gammapy.image.models.Gauss2DPDF`
Convolution of both Gaussians.
"""
sigma = np.asarray(sigma, dtype=np.float64)
new_sigma = np.sqrt(self._sigma2 + sigma ** 2)
return Gauss2DPDF(new_sigma)
[docs]class MultiGauss2D(object):
"""Sum of multiple 2D Gaussians.
Parameters
----------
sigmas : `~numpy.ndarray`
widths of the Gaussians to add
norms : `~numpy.ndarray`, optional
normalizations of the Gaussians to add
Notes
-----
* This sum is no longer a PDF, it is not normalized to 1.
* The "norm" of each component represents the 2D integral,
not the amplitude at the origin.
"""
def __init__(self, sigmas, norms=None):
# If no norms are given, you have a PDF.
sigmas = np.asarray(sigmas, dtype=np.float64)
self.components = [Gauss2DPDF(sigma) for sigma in sigmas]
if norms is None:
self.norms = np.ones(len(self.components))
else:
self.norms = np.asarray(norms, dtype=np.float64)
[docs] def __call__(self, x, y=0):
"""dp / (dx dy) at position (x, y)
Parameters
----------
x : `~numpy.ndarray`
x coordinate
y : `~numpy.ndarray`, optional
y coordinate
Returns
-------
total : `~numpy.ndarray`
dp / (dx dy)
"""
x = np.asarray(x, dtype=np.float64)
y = np.asarray(y, dtype=np.float64)
total = np.zeros_like(x)
for norm, component in zip(self.norms, self.components):
total += norm * component(x, y)
return total
@property
def n_components(self):
"""Number of components (int)"""
return len(self.components)
@property
def sigmas(self):
"""Array of Gaussian widths (`~numpy.ndarray`)"""
return np.array([_.sigma for _ in self.components])
@property
def integral(self):
"""Integral as sum of norms (`~numpy.ndarray`)"""
return np.nansum(self.norms)
@property
def amplitude(self):
"""Amplitude at the center (float)"""
return self.__call__(0, 0)
@property
def max_sigma(self):
"""Largest Gaussian width (float)"""
return self.sigmas.max()
@property
def eff_sigma(self):
r"""Effective Gaussian width for single-Gauss approximation (float)
Notes
-----
The effective Gaussian width is given by:
.. math:: \sigma_\mathrm{eff} = \sqrt{\sum_i N_i \sigma_i^2}
where ``N`` is normalization and ``sigma`` is width.
"""
sigma2s = np.array([component._sigma2 for component
in self.components])
return np.sqrt(np.sum(self.norms * sigma2s))
[docs] def dpdtheta2(self, theta2):
"""dp / dtheta2 at position theta2 = theta ^ 2
Parameters
----------
theta2 : `~numpy.ndarray`
Offset squared
Returns
-------
dpdtheta2 : `~numpy.ndarray`
dp / dtheta2
"""
# Actually this is only a PDF if sum(norms) == 1
theta2 = np.asarray(theta2, dtype=np.float64)
total = np.zeros_like(theta2)
for norm, component in zip(self.norms, self.components):
total += norm * component.dpdtheta2(theta2)
return total
[docs] def normalize(self):
"""Normalize function.
Returns
-------
norm_multigauss : `~gammapy.image.models.MultiGauss2D`
normalized function
"""
sum = self.integral
if sum!=0:
self.norms /= sum
return self
[docs] def containment_fraction(self, theta):
"""Containment fraction.
Parameters
----------
theta : `~numpy.ndarray`
Offset
Returns
-------
containment_fraction : `~numpy.ndarray`
Containment fraction
"""
theta = np.asarray(theta, dtype=np.float64)
total = np.zeros_like(theta)
for norm, component in zip(self.norms, self.components):
total += norm * component.containment_fraction(theta)
return total
[docs] def containment_radius(self, containment_fraction):
"""Containment angle for a given containment fraction.
Parameters
----------
containment_fraction : `~numpy.ndarray`
Containment fraction
Returns
-------
containment_radius : `~numpy.ndarray`
Containment radius
"""
# I had big problems with fsolve running into negative thetas.
# So instead I'll find a theta_max myself so that theta
# is in the interval [0, theta_max] and then use good ol brentq
if not containment_fraction < self.integral:
raise ValueError('containment_fraction = {} not possible for integral = {}'
''.format(containment_fraction, self.integral))
from scipy.optimize import brentq
def f(theta):
# positive if theta too large
return self.containment_fraction(theta) - containment_fraction
# TODO: if it is an array we have to loop by hand!
# containment = np.asarray(containment, dtype=np.float64)
# Inital guess for theta
theta_max = self.eff_sigma
# Expand until we really find a theta_max
while f(theta_max) < 0:
theta_max *= 2
return brentq(f, a=0, b=theta_max)
[docs] def match_sigma(self, containment_fraction):
"""Compute equivalent Gauss width.
Find the sigma of a single-Gaussian distribution that
approximates this one, such that theta matches for a given
containment.
Parameters
----------
containment_fraction : `~numpy.ndarray`
Containment fraction
Returns
-------
sigma : `~numpy.ndarray`
Equivalent containment radius
"""
theta1 = self.containment_radius(containment_fraction)
theta2 = Gauss2DPDF(sigma=1).containment_radius(containment_fraction)
return theta1 / theta2
[docs] def gauss_convolve(self, sigma, norm=1):
"""Convolve with another Gauss.
Compute new norms and sigmas of all the components such that
the new distribution represents the convolved old distribution
by a Gaussian of width sigma and then multiplied by norm.
This MultiGauss2D is unchanged, a new one is created and returned.
This is useful if you need to e.g. compute theta for one PSF
and many sigmas.
Parameters
----------
sigma : `~numpy.ndarray` or float
Gaussian width of the new Gaussian 2D PDF to covolve with.
norm : `~numpy.ndarray` or float
Normalization of the new Gaussian 2D PDF to covolve with.
Returns
-------
new_multi_gauss_2d : `~gammapy.image.models.MultiGauss2D`
Convolution as new MultiGauss2D
"""
sigma = np.asarray(sigma, dtype=np.float64)
norm = np.asarray(norm, dtype=np.float64)
sigmas, norms = [], []
for ii in range(self.n_components):
sigmas.append(self.components[ii].gauss_convolve(sigma).sigma)
norms.append(self.norms[ii] * norm)
return MultiGauss2D(sigmas, norms)
[docs]def gaussian_sum_moments(F, sigma, x, y, cov_matrix, shift=0.5):
"""Compute image moments with uncertainties for sum of Gaussians.
The moments are computed analytically, the formulae are documented below.
Calls ``uncertainties.correlated_values`` to propagate the errors.
Parameters
----------
F : array
Integral norms of the Gaussian components.
sigmas : array
Widths of the Gaussian components.
x : array
x positions of the Gaussian components.
y : array
y positions of the Gaussian components.
cov_matrix : array
Covariance matrix of the parameters. The columns have to follow the order:
[sigma_1, x_1, y_1, F_1, sigma_2, x_2, y_2, F_2, ..., sigma_N, x_N, y_N, F_N]
shift : float (default = 0.5)
Depending on where the image values are given, the grid has to be
shifted. If the values are given at the center of the pixel
shift = 0.5.
Returns
-------
nominal_values : list
List of image moment nominal values:
[F_sum, x_sum, y_sum, x_sigma, y_sigma, sqrt(x_sigma * y_sigma)]
All values are given in pixel coordinates.
std_devs : list
List of image moment standard deviations.
Examples
--------
A simple example for an image consisting of three Gaussians
with zero covariance matrix:
>>> import numpy as np
>>> from gammapy.image.models.gauss import gaussian_sum_moments
>>> cov_matrix = np.zeros((12, 12))
>>> F = [100, 200, 300]
>>> sigma = [15, 10, 5]
>>> x = [100, 120, 70]
>>> y = [100, 90, 120]
>>> nominal_values, std_devs = gaussian_sum_moments(F, sigma, x, y, cov_matrix)
Notes
-----
The 0th moment (total flux) is given by:
.. math::
F_{\\Sigma} = \\int_{-\\infty}^{\\infty}f_{\\Sigma}(x, y)dx dy =
\\sum_i^N F_i
The 1st moments (position) are given by:
.. math::
x_{\\Sigma} = \\frac{1}{F_{\\Sigma}} \\int_{-\\infty}^{\\infty}x
f_{\\Sigma}(x, y)dx dy = \\frac{1}{F_{\\Sigma}}\\sum_i^N x_iF_i
y_{\\Sigma} = \\frac{1}{F_{\\Sigma}} \\int_{-\\infty}^{\\infty}y
f_{\\Sigma}(x, y)dx dy = \\frac{1}{F_{\\Sigma}}\\sum_i^N y_iF_i
The 2nd moments (extension) are given by:
.. math::
\\sigma_{\\Sigma_x}^2 = \\frac{1}{F_{\\Sigma}} \\sum_i^N F_i
\\cdot (\\sigma_i^2 + x_i^2) - x_{\\Sigma}^2
\\sigma_{\\Sigma_y}^2 = \\frac{1}{F_{\\Sigma}} \\sum_i^N F_i
\\cdot (\\sigma_i^2 + y_i^2) - y_{\\Sigma}^2
"""
import uncertainties
# Check input arrays
if len(F) != len(sigma) or len(F) != len(x) or len(F) != len(x):
raise Exception("Input arrays have to have the same size")
# Order parameter values
values = []
for i in range(len(F)):
values += [sigma[i], x[i], y[i], F[i]]
# Set up parameters with uncertainties
parameter_list = uncertainties.correlated_values(values, cov_matrix)
# Reorder parameters by splitting into 4-tuples
parameters = list(zip(*[iter(parameter_list)] * 4))
def zero_moment(parameters):
"""0th moment of the sum of Gaussian components."""
F_sum = 0
for component in parameters:
F_sum += component[3]
return F_sum
def first_moment(parameters, F_sum, shift):
"""1st moment of the sum of Gaussian components."""
x_sum, y_sum = 0, 0
for component in parameters:
x_sum += (component[1] + shift) * component[3]
y_sum += (component[2] + shift) * component[3]
return x_sum / F_sum, y_sum / F_sum
def second_moment(parameters, F_sum, x_sum, y_sum, shift):
"""2nd moment of the sum of Gaussian components."""
var_x_sum, var_y_sum = 0, 0
for p in parameters:
var_x_sum += ((p[1] + shift) ** 2 + p[0] ** 2) * p[3]
var_y_sum += ((p[2] + shift) ** 2 + p[0] ** 2) * p[3]
return var_x_sum / F_sum - x_sum ** 2, var_y_sum / F_sum - y_sum ** 2
# Compute moments
F_sum = zero_moment(parameters)
x_sum, y_sum = first_moment(parameters, F_sum, shift)
var_x_sum, var_y_sum = second_moment(parameters, F_sum, x_sum, y_sum, shift)
# Return values and stddevs separately
values = [F_sum, x_sum, y_sum, var_x_sum ** 0.5, var_y_sum ** 0.5, (var_x_sum * var_y_sum) ** 0.25]
nominal_values = [_.nominal_value for _ in values]
std_devs = [float(_.std_dev) for _ in values]
return nominal_values, std_devs