Source code for hierarc.Likelihood.anisotropy_scaling
__author__ = 'sibirrer'
from scipy.interpolate import interp1d, interp2d
import numpy as np
[docs]class AnisotropyScalingSingleAperture(object):
"""
class to manage anisotropy scaling for single slit observation
"""
def __init__(self, ani_param_array, ani_scaling_array):
"""
:param ani_param_array: array of anisotropy parameter value
:param ani_scaling_array: array with the scaling of J() for single slit
"""
self._evalute_ani = False
if ani_param_array is not None and ani_scaling_array is not None:
if isinstance(ani_param_array, list):
self._dim_scaling = len(ani_param_array)
else:
self._dim_scaling = 1
if self._dim_scaling == 1:
self._f_ani = interp1d(ani_param_array, ani_scaling_array, kind='linear')
elif self._dim_scaling == 2:
self._f_ani = interp2d(ani_param_array[0], ani_param_array[1], ani_scaling_array.T)
else:
raise ValueError('anisotropy scaling with dimension %s not supported.' % self._dim_scaling)
self._evalute_ani = True
[docs] def ani_scaling(self, aniso_param_array):
"""
:param aniso_param_array: anisotropy parameter array
:return: scaling J(a_ani) for single slit
"""
if self._evalute_ani is not True or aniso_param_array is None:
return 1
if self._dim_scaling == 1:
return self._f_ani(aniso_param_array[0])
elif self._dim_scaling == 2:
return self._f_ani(aniso_param_array[0], aniso_param_array[1])[0]
[docs]class AnisotropyScalingIFU(object):
"""
class to manage anisotropy scalings for IFU data
"""
def __init__(self, anisotropy_model='NONE', ani_param_array=None, ani_scaling_array_list=None):
"""
:param anisotropy_model: string, either 'NONE', 'OM' or 'GOM'
:param ani_param_array: array of anisotropy parameter value (1d for 'OM' model, 2d for 'GOM' model)
:param ani_scaling_array_list: list of array with the scalings of J() for each IFU
"""
self._anisotropy_model = anisotropy_model
self._evalute_ani = False
if ani_param_array is not None and ani_scaling_array_list is not None and self._anisotropy_model != 'NONE':
self._evalute_ani = True
self._anisotropy_scaling_list = []
self._f_ani_list = []
for ani_scaling_array in ani_scaling_array_list:
self._anisotropy_scaling_list.append(AnisotropyScalingSingleAperture(ani_param_array=ani_param_array,
ani_scaling_array=ani_scaling_array))
if isinstance(ani_param_array, list):
self._dim_scaling = len(ani_param_array)
else:
self._dim_scaling = 1
if self._dim_scaling == 1 and anisotropy_model in ['OM', 'const']:
self._ani_param_min = np.min(ani_param_array)
self._ani_param_max = np.max(ani_param_array)
elif self._dim_scaling == 2 and anisotropy_model == 'GOM':
self._ani_param_min = [min(ani_param_array[0]), min(ani_param_array[1])]
self._ani_param_max = [max(ani_param_array[0]), max(ani_param_array[1])]
else:
raise ValueError('anisotropy scaling with dimension %s does not match anisotropy model %s'
% (self._dim_scaling, self._anisotropy_model))
[docs] def ani_scaling(self, aniso_param_array):
"""
:param aniso_param_array: anisotropy parameter array
:return: scaling J(a_ani) for the IFU's
"""
if self._evalute_ani is not True or aniso_param_array is None:
return [1]
scaling_list = []
for scaling_class in self._anisotropy_scaling_list:
scaling = scaling_class.ani_scaling(aniso_param_array)
scaling_list.append(scaling)
return np.array(scaling_list)
[docs] def draw_anisotropy(self, a_ani=None, a_ani_sigma=0, beta_inf=None, beta_inf_sigma=0):
"""
draw Gaussian distribution and re-sample if outside bounds
:param a_ani: mean of the distribution
:param a_ani_sigma: std of the distribution
:param beta_inf: anisotropy at infinity (relevant for GOM model)
:param beta_inf_sigma: std of beta_inf distribution
:return: random draw from the distribution
"""
if self._anisotropy_model in ['OM', 'const']:
if a_ani < self._ani_param_min or a_ani > self._ani_param_max:
raise ValueError('anisotropy parameter is out of bounds of the interpolated range!')
# we draw a linear gaussian for 'const' anisotropy and a scaled proportional one for 'OM
if self._anisotropy_model == 'OM':
a_ani_draw = np.random.normal(a_ani, a_ani_sigma*a_ani)
else:
a_ani_draw = np.random.normal(a_ani, a_ani_sigma)
if a_ani_draw < self._ani_param_min or a_ani_draw > self._ani_param_max:
return self.draw_anisotropy(a_ani, a_ani_sigma)
return np.array([a_ani_draw])
elif self._anisotropy_model in ['GOM']:
if a_ani < self._ani_param_min[0] or a_ani > self._ani_param_max[0] or beta_inf < self._ani_param_min[1] or beta_inf > self._ani_param_max[1]:
raise ValueError('anisotropy parameter is out of bounds of the interpolated range!')
a_ani_draw = np.random.normal(a_ani, a_ani_sigma*a_ani)
beta_inf_draw = np.random.normal(beta_inf, beta_inf_sigma)
if a_ani_draw < self._ani_param_min[0] or a_ani_draw > self._ani_param_max[0] or beta_inf_draw < self._ani_param_min[1] or beta_inf_draw > self._ani_param_max[1]:
return self.draw_anisotropy(a_ani, a_ani_sigma, beta_inf, beta_inf_sigma)
return np.array([a_ani_draw, beta_inf_draw])
return None