Commit 47df8d7c authored by bepinat's avatar bepinat

Update with minor corrections

parent 542d229f
import numpy as np
from astropy.io import fits
from scipy.interpolate import interp2d
import Tools.tools as tools
class Cube:
"""
Stocks a fits file and determines all parameters like size etc ...
:param string or ndarray filename: path+name of the file, can be a numpy array.
"""
def __init__(self, filename, mask=None):
self.data_rebin = None
if type(filename) is str:
self.header = fits.getheader(filename)
self.data = fits.getdata(filename)
self.length = self.header['NAXIS1']
self.height = self.header['NAXIS2']
self.sizez = self.header['NAXIS3']
self.size = np.array(np.shape(self.data))
self.len = self.length*self.height
self.oversample = 1
if mask is None:
self.mask = np.logical_not(np.isnan(self.data))
else:
self.mask = mask
self.nan_to_num()
# Add for create model without Image's class
if type(filename) is np.ndarray:
self.header = None
self.data = np.nan_to_num(filename)
self.size = np.array(np.shape(self.data))
self.height = self.size[2]
self.length = self.size[1]
self.sizez = self.size[0]
self.len = self.length*self.height
self.oversample = 1
if mask is None:
self.mask = np.logical_not(np.isnan(self.data))
else:
self.mask = mask
self.nan_to_num()
def nan_to_num(self):
"""
Converts 'nan' into 0
"""
self.data = np.nan_to_num(self.data)
def get_size(self):
return self.size
def get_lenght(self):
return self.length
def get_height(self):
return self.height
def set_oversamp(self, oversamp):
self.oversample = oversamp
def get_oversamp(self):
return self.oversample
......@@ -16,8 +16,8 @@ class Image:
if type(filename) is str:
self.header = fits.getheader(filename)
self.data = fits.getdata(filename)
self.height = self.header['NAXIS1']
self.length = self.header['NAXIS2']
self.length = self.header['NAXIS1']
self.height = self.header['NAXIS2']
self.size = np.array(np.shape(self.data))
self.len = self.length*self.height
self.oversample = 1
......
import math
import numpy as np
import Tools.tools as tools
from scipy import constants as ct
from Class.PSF import PSF
from Class.Images import Image
......@@ -73,15 +74,6 @@ class Model2D:
self.model_RCparams = None
def set_constraints(self, config):
"""
Sets the value of the model parameters
:param float xc: abscissa of the center in pixel
"""
return
def set_parameters(self, xc, yc, pa, inc, vs, *args):
"""
Sets the value of the model parameters
......@@ -228,6 +220,197 @@ class Model2D:
return np.sum(chi2)
#class Model3D: # XXX do the same thing: initiate, then produce velocity map, intensity/flux map and velocity dispersion map at high resolution and then produce a cube with the same sampling as data and adjusted with respect to a reference line and wavelengths in the cube. Then convolve with PSF and rebin. Everything needed is in the program with made with Debora
class Model3D: # XXX do the same thing: initiate, then produce velocity map, intensity/flux map and velocity dispersion map at high resolution and then produce a cube with the same sampling as data and adjusted with respect to a reference line and wavelengths in the cube. Then convolve with PSF and rebin. Everything needed is in the program with made with Debora
"""
Models in 3D of the datacube around a line of galaxies
Also compute the velocity field and the dispersion of a model using observed cube in input.
This model can be used with any rotational curve models (in 'velocity_model.py').
"""
def __init__(self, cube, vcube, flux, psf, vel_model, psfz=0, sig0=0, slope=0., header=None, z=0., lbda0=6562.78, cont=0.):
"""
:param Image cube: represents the cube to fit
:param Image vcube: the variance cube
:param Image flux: flux distribution at the same or in high resolution than the velocity
:param PSF psf: the class of the psf
:param func vel_model: velocity function used to create the model
:param float psfz: velocity dispersion corresponding to the spectral resolution
:param float sig0: velocity dispersion of the model
:param float slope: slope of the velocity dispersion
:param hdufits header header: header associated to the datacube
:param float z: systemic redshift of the source
:param float lbda0: restframe wavelength of the line used (could be a list of lines and then need to add line ratios, etc.)
:param float cont: continuum level (assumed constant for now, but could be polynomial)
"""
# Parameters which must be initialized
self.vel_model = vel_model.vel_model
self.model_sig0 = sig0
self.model_slope = slope
self.psfz = psfz
# XXX a reprendre en 3D
self.cube_mod = np.zeros(cube.shape)
#self.vel_map = np.zeros(cube.shape[1:])
self.vel_map_hd = np.zeros(np.array(flux.get_size()))
self.disp_map_hd = np.zeros(np.array(flux.get_size()))
#self.cube_hd = np.zeros((cube.shape[0], np.array(flux.get_size()))) # Check def works (dimension declaration)
#
self.cube = cube
self.vcube = vcube
self.header = header
self.flux = flux
self.psf = psf
self.oversamp = self.flux.get_oversamp()
self.indices = tools.lowres_indices_highres_frame(flux.get_size(), self.oversamp)
self.z = z
self.lbda0 = lbda0
self.lbda1 = lbda0 * (1 + z)
self.cont = cont # XXX devrait être une image en fait!
# We initiate the wavelength range of the datacube
self.lbda = (np.arange(header['NAXIS3']) - header['CRPIX3']) * header['CDELT3'] + header['CRVAL3'] # XXX Check units!
self.vind = (self.lbda - self.lbda1) / self.lbda1 * ct.c * 1e-3 # km/s ; XXX Ne fonctionne que s'il n'y a qu'une seule raie
# We initiate the geometrical model parameters
self.model_pa = None
self.model_inc = None
self.model_xc = None
self.model_yc = None
self.model_vs = None
# We initiate the associated coordinate maps in the galaxy frame
self.model_radius = None
self.model_theta = None
# We write an ordered list containing the names of the model parameters
self.model_parname = ['xc', 'yc', 'pa', 'inc', 'vs']
self.model_modname = ['geometrical', 'geometrical', 'geometrical', 'geometrical', 'geometrical']
self.model_compname = ['geometrical', 'geometrical', 'geometrical', 'geometrical', 'geometrical']
# We initiate a dictionary containing the names of rotation curve parameters. The number of arguments adjusts automatically
#ii = 0
for name in vel_model.model_parname:
self.model_parname.append(name)
for name in vel_model.model_modname:
self.model_modname.append(name)
for name in vel_model.model_compname:
self.model_compname.append(name)
#ii += 1
#print(self.model_parname)
#for name in self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount]:
#self.model_parname.append(name)
#self.model_RCparams_names = self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount] # can be useful to check the parameters are provided with the correct order
#self.model_RCparams_values = [None for par in self.model_RCparams_names]
#self.model_RCparams = dict((par, None) for par in self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount])
self.model_RCparams = None
def set_parameters(self, xc, yc, pa, inc, vs, *args):
"""
Sets the value of the model parameters
:param float xc: abscissa of the center in pixel
:param float yc: ordinate of the center in pixel
:param float pa: position angle of the major axis in degree
:param float inc: inclination of the disk in degree
:param float vs: systemic velocity in km/s
:param list of float args: list that contains other model arguments
"""
self.model_pa = pa
self.model_inc = inc
self.model_xc = xc
self.model_yc = yc
self.model_vs = vs
#self.model_radius, self.model_theta = tools.sky_coord_to_galactic(self.model_xc, self.model_yc, self.model_pa, self.model_inc, im_size=np.shape(self.vel_map_hd), oversamp=self.oversamp)
self.model_radius, self.model_theta = tools.sky_coord_to_galactic(self.indices[1], self.indices[0], self.model_xc, self.model_yc, self.model_pa, self.model_inc)
self.model_RCparams = args
def get_parameter(self):
"""
Gets the actual parameters of the model (in low resolution scale)
:return ndarray:
"""
return [self.model_xc, self.model_yc, self.model_pa, self.model_inc, self.model_vs, *self.model_RCparams]
def disk_velocity(self):
"""
Computes the velocity field
:return ndarray:
"""
vr = self.vel_model(self.model_radius, *self.model_RCparams)
# Calculation of the velocity field
v = vr * math.sin(math.radians(self.model_inc)) * self.model_theta + self.model_vs
return v
def linear_velocity_dispersion(self):
"""
Returns the velocity dispersion map needed for the fit process
:return ndarray: velocity dispersion map
"""
# Calculation of the velocity dispersion
sig = self.model_sig0 + self.model_slope * np.abs(self.model_radius)
sig[np.where(sig <= 0)] = 0
return sig
def cube_model(self):
"""
Computes a high resolution cube
: return cube:
"""
self.vel_map_hd = self.disk_velocity()
self.disp_map_hd = self.linear_velocity_dispersion() # besoin du centre du modèle, donc ici! En principe pour la mmodélisation du cube, il faudrait que sigma et slope soient des paramètres libre du modèle
# Il faudrait en principe donner la LSF en input (de la même manière que la PSF). Dans un premier temps, on peut faire l'hypothèse que la LSF est Gaussienne, c'est donc inclus dans sigma
self.cube_hd = self.flux * np.exp(-np.subtract(self.vind.reshape(len(vind), 1, 1), self.vel_map_hd.reshape(1, self.flux.shape[0], self.flux.shape[1])) ** 2 / (2 * self.disp_map_hd ** 2)) + self.cont
self.cube_mod = tools.rebin_data(self.psf.convolution(self.cube_hd), self.oversamp)
def least_square(self, p, fjac=None):
"""
Function minimized by mpfit.
Returns (data-model)^2/err^2 XXX Why no square???
:param ndarray p: array of parameters
:return ndarray:
"""
self.set_parameters(*p)
self.cube_model()
return [0, np.reshape((self.cube - self.cube_mod)/np.sqrt(self.vcube), -1)] # sqrt, to have uncertainty, i.e. stdev
def log_likelihood(self, cube, ndim, nparams):
"""
Log likelihood function which maximized by multinest
Returns -sum[(data-model)^2/(2*err^2)]
:param ndarray cube: data whith n_params dimension
:param int ndim: number of dimension if different of the number of paramerters
:param int nparams: number of parameters
:return float:
"""
#self.set_parameters(cube[0], cube[1], cube[2], cube[3], cube[4], cube[5], cube[6])
self.set_parameters(cube[0], cube[1], cube[2], cube[3], cube[4], *cube[5:ndim])
self.cube_model()
chi2 = -(self.cube_mod - self.cube)**2 / (2*self.vcube) # no **2 because we use variance
return np.sum(chi2)
#class Model1D:
import math
import numpy as np
import Tools.tools as tools
from scipy import constants as ct
from Class.PSF import PSF
from Class.Images import Image
class Model2D:
"""
Models in 2D of the velocity field of galaxies
Compute the velocity field and the dispersion of a model using observed data in input.
This model can be used with any rotational curve models (in 'velocity_model.py').
"""
def __init__(self, vel, errvel, flux, psf, vel_model, psfz=0, sig0=0, slope=0.):
"""
:param Image vel: represent the velocity field to fit
:param Image errvel: the error map of the velocity (vel)
:param Image flux: flux distribution at the same or in high resolution than the velocity
:param PSF psf: the class of the psf
:param func vel_model: velocity function used to create the model
:param float psfz: velocity dispersion corresponding to the spectral resolution
:param float sig0: velocity dispersion of the model
:param float slope: slope of the velocity dispersion
"""
# Parameters which must be initialized
#self.vel_model = vel_model
self.vel_model = vel_model.vel_model
self.model_sig0 = sig0
self.model_slope = slope
self.psfz = psfz
self.vel_map = np.zeros(vel.get_size())
self.vel_disp = np.zeros(vel.get_size())
self.vel_map_hd = np.zeros(np.array(flux.get_size()))
self.vel = vel
self.errvel = errvel
self.flux = flux
self.psf = psf
self.oversamp = self.flux.get_oversamp()
self.indices = tools.lowres_indices_highres_frame(flux.get_size(), self.oversamp)
# We initiate the geometrical model parameters
self.model_pa = None
self.model_inc = None
self.model_xc = None
self.model_yc = None
self.model_vs = None
# We initiate the associated coordinate maps in the galaxy frame
self.model_radius = None
self.model_theta = None
# We write an ordered list containing the names of the model parameters
self.model_parname = ['xc', 'yc', 'pa', 'inc', 'vs']
self.model_modname = ['geometrical', 'geometrical', 'geometrical', 'geometrical', 'geometrical']
self.model_compname = ['geometrical', 'geometrical', 'geometrical', 'geometrical', 'geometrical']
# We initiate a dictionary containing the names of rotation curve parameters. The number of arguments adjusts automatically
#ii = 0
for name in vel_model.model_parname:
self.model_parname.append(name)
for name in vel_model.model_modname:
self.model_modname.append(name)
for name in vel_model.model_compname:
self.model_compname.append(name)
#ii += 1
#print(self.model_parname)
#for name in self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount]:
#self.model_parname.append(name)
#self.model_RCparams_names = self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount] # can be useful to check the parameters are provided with the correct order
#self.model_RCparams_values = [None for par in self.model_RCparams_names]
#self.model_RCparams = dict((par, None) for par in self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount])
self.model_RCparams = None
def set_parameters(self, xc, yc, pa, inc, vs, *args):
"""
Sets the value of the model parameters
:param float xc: abscissa of the center in pixel
:param float yc: ordinate of the center in pixel
:param float pa: position angle of the major axis in degree
:param float inc: inclination of the disk in degree
:param float vs: systemic velocity in km/s
:param list of float args: list that contains other model arguments
"""
self.model_pa = pa
self.model_inc = inc
self.model_xc = xc
self.model_yc = yc
self.model_vs = vs
#self.model_radius, self.model_theta = tools.sky_coord_to_galactic(self.model_xc, self.model_yc, self.model_pa, self.model_inc, im_size=np.shape(self.vel_map_hd), oversamp=self.oversamp)
self.model_radius, self.model_theta = tools.sky_coord_to_galactic(self.indices[1], self.indices[0], self.model_xc, self.model_yc, self.model_pa, self.model_inc)
self.model_RCparams = args
def get_parameter(self):
"""
Gets the actual parameters of the model (in low resolution scale)
:return ndarray:
"""
return [self.model_xc, self.model_yc, self.model_pa, self.model_inc, self.model_vs, *self.model_RCparams]
def disk_velocity(self):
"""
Computes the velocity field
:return ndarray:
"""
vr = self.vel_model(self.model_radius, *self.model_RCparams)
# Calculation of the velocity field
v = vr * math.sin(math.radians(self.model_inc)) * self.model_theta + self.model_vs
return v
def velocity_map(self):
"""
Computes the velocity map of the model
"""
self.vel_map_hd = self.disk_velocity()
vel_times_flux = tools.rebin_data(self.psf.convolution(self.flux.data * self.vel_map_hd), self.oversamp)
self.vel_map[self.vel.mask] = vel_times_flux[self.vel.mask] / self.flux.data_rebin[self.vel.mask]
def linear_velocity_dispersion(self):
"""
Returns the velocity dispersion map needed for the fit process
:return ndarray: velocity dispersion map
"""
# Calculation of the velocity dispersion
sig = self.model_sig0 + self.model_slope * np.abs(self.model_radius)
sig[np.where(sig <= 0)] = 0
return sig
def vel_disp_map(self):
"""
Computes the velocity dispersion map from the velocity field (with bestfit parameters) and the flux distribution
:return ndarray: velocity dispersion map
"""
term2 = np.zeros(self.vel.size)
term3 = np.zeros(self.vel.size)
term2[self.vel.mask] = tools.rebin_data(self.psf.convolution(self.vel_map_hd ** 2 * self.flux.data), self.flux.oversample)[self.vel.mask] \
/ self.flux.data_rebin[self.vel.mask]
term3[self.vel.mask] = (tools.rebin_data(self.psf.convolution(self.vel_map_hd * self.flux.data), self.flux.oversample)[self.vel.mask] /
self.flux.data_rebin[self.vel.mask]) ** 2
self.vel_disp = np.sqrt(term2 - term3 + self.psfz**2)
def vel_disp_map_model(self):
"""
Computes the velocity dispersion map from the velocity field (with bestfit parameters) and the flux distribution
:return ndarray: velocity dispersion map
"""
term1 = np.zeros(self.vel.size)
term2 = np.zeros(self.vel.size)
term3 = np.zeros(self.vel.size)
sig = self.linear_velocity_dispersion()
term1[self.vel.mask] = tools.rebin_data(self.psf.convolution(sig ** 2 * self.flux.data), self.flux.oversample)[self.vel.mask] \
/ self.flux.data_rebin[self.vel.mask]
term2[self.vel.mask] = tools.rebin_data(self.psf.convolution(self.vel_map_hd ** 2 * self.flux.data), self.flux.oversample)[self.vel.mask] \
/ self.flux.data_rebin[self.vel.mask]
term3[self.vel.mask] = (tools.rebin_data(self.psf.convolution(self.vel_map_hd * self.flux.data), self.flux.oversample)[self.vel.mask] /
self.flux.data_rebin[self.vel.mask]) ** 2
self.vel_disp = np.sqrt(term1 + term2 - term3 + self.psfz**2)
def least_square(self, p, fjac=None):
"""
Function minimized by mpfit.
Returns (data-model)^2/err^2
:param ndarray p: array of parameters
:return ndarray:
"""
self.set_parameters(*p)
self.velocity_map()
return [0, np.reshape((self.vel.data[self.vel.mask]-self.vel_map[self.vel.mask])/self.errvel.data[self.vel.mask], -1)]
def log_likelihood(self, cube, ndim, nparams):
"""
Log likelihood function which maximized by multinest
Returns -sum[(data-model)^2/(2*err^2)]
:param ndarray cube: data whith n_params dimension
:param int ndim: number of dimension if different of the number of paramerters
:param int nparams: number of parameters
:return float:
"""
#self.set_parameters(cube[0], cube[1], cube[2], cube[3], cube[4], cube[5], cube[6])
self.set_parameters(cube[0], cube[1], cube[2], cube[3], cube[4], *cube[5:ndim])
self.velocity_map()
chi2 = -(self.vel_map[self.vel.mask] - self.vel.data[self.vel.mask])**2 / (2*self.errvel.data[self.vel.mask]**2)
return np.sum(chi2)
class Model3D: # XXX do the same thing: initiate, then produce velocity map, intensity/flux map and velocity dispersion map at high resolution and then produce a cube with the same sampling as data and adjusted with respect to a reference line and wavelengths in the cube. Then convolve with PSF and rebin. Everything needed is in the program with made with Debora
"""
Models in 3D of the datacube around a line of galaxies
Also compute the velocity field and the dispersion of a model using observed cube in input.
This model can be used with any rotational curve models (in 'velocity_model.py').
"""
def __init__(self, cube, vcube, flux, psf, vel_model, psfz=0, sig0=0, slope=0., z=0., lbda0=6562.78, cont=0.):
"""
:param Image cube: represents the cube to fit
:param Image vcube: the variance cube
:param Image flux: flux distribution at the same or in high resolution than the velocity
:param PSF psf: the class of the psf
:param func vel_model: velocity function used to create the model
:param float psfz: velocity dispersion corresponding to the spectral resolution
:param float sig0: velocity dispersion of the model
:param float slope: slope of the velocity dispersion
:param hdufits header header: header associated to the datacube
:param float z: systemic redshift of the source
:param float lbda0: restframe wavelength of the line used (could be a list of lines and then need to add line ratios, etc.)
:param float cont: continuum level (assumed constant for now, but could be polynomial)
"""
# Parameters which must be initialized
self.vel_model = vel_model.vel_model
self.model_sig0 = sig0
self.model_slope = slope
self.psfz = psfz
self.cube_mod = np.zeros(cube.get_size())
#self.vel_map = np.zeros(cube.shape[1:])
self.vel_mod_hd = np.zeros(np.array(flux.get_size()))
self.disp_mod_hd = np.zeros(np.array(flux.get_size()))
self.disp_map_hd = np.zeros(np.array(flux.get_size()))
#self.cube_hd = np.zeros((cube.shape[0], np.array(flux.get_size()))) # Check def works (dimension declaration)
#
self.cube = cube
self.vcube = vcube
self.flux = flux
self.psf = psf
self.oversamp = self.flux.get_oversamp()
self.indices = tools.lowres_indices_highres_frame(flux.get_size(), self.oversamp)
self.z = z
self.lbda0 = lbda0
self.lbda1 = lbda0 * (1 + z)
self.cont = cont # XXX devrait être une image en fait!
# We initiate the wavelength range of the datacube
header = cube.header
self.lbda = (np.arange(header['NAXIS3']) + 1 - header['CRPIX3']) * header['CDELT3'] + header['CRVAL3'] # XXX Check units!
self.vind = (self.lbda - self.lbda1) / self.lbda1 * ct.c * 1e-3 # km/s ; XXX Ne fonctionne que s'il n'y a qu'une seule raie
# We initiate the geometrical model parameters
self.model_pa = None
self.model_inc = None
self.model_xc = None
self.model_yc = None
self.model_vs = None
# We initiate the associated coordinate maps in the galaxy frame
self.model_radius = None
self.model_theta = None
# We write an ordered list containing the names of the model parameters
self.model_parname = ['xc', 'yc', 'pa', 'inc', 'vs']
self.model_modname = ['geometrical', 'geometrical', 'geometrical', 'geometrical', 'geometrical']
self.model_compname = ['geometrical', 'geometrical', 'geometrical', 'geometrical', 'geometrical']
# We initiate a dictionary containing the names of rotation curve parameters. The number of arguments adjusts automatically
#ii = 0
for name in vel_model.model_parname:
self.model_parname.append(name)
for name in vel_model.model_modname:
self.model_modname.append(name)
for name in vel_model.model_compname:
self.model_compname.append(name)
#ii += 1
#print(self.model_parname)
#for name in self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount]:
#self.model_parname.append(name)
#self.model_RCparams_names = self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount] # can be useful to check the parameters are provided with the correct order
#self.model_RCparams_values = [None for par in self.model_RCparams_names]
#self.model_RCparams = dict((par, None) for par in self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount])
self.model_RCparams = None
def set_parameters(self, xc, yc, pa, inc, vs, *args):
"""
Sets the value of the model parameters
:param float xc: abscissa of the center in pixel
:param float yc: ordinate of the center in pixel
:param float pa: position angle of the major axis in degree
:param float inc: inclination of the disk in degree
:param float vs: systemic velocity in km/s
:param list of float args: list that contains other model arguments
"""
self.model_pa = pa
self.model_inc = inc
self.model_xc = xc
self.model_yc = yc
self.model_vs = vs
#self.model_radius, self.model_theta = tools.sky_coord_to_galactic(self.model_xc, self.model_yc, self.model_pa, self.model_inc, im_size=np.shape(self.vel_mod_hd), oversamp=self.oversamp)
self.model_radius, self.model_theta = tools.sky_coord_to_galactic(self.indices[1], self.indices[0], self.model_xc, self.model_yc, self.model_pa, self.model_inc)
self.model_RCparams = args
def get_parameter(self):
"""
Gets the actual parameters of the model (in low resolution scale)
:return ndarray:
"""