Commit 76fd2e22 authored by bepinat's avatar bepinat

new project

parent 0bc485a4
......@@ -6,10 +6,10 @@ import Tools.tools as tools
class Image:
"""
Stock fits file and determine all parameters like size etc ...
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.
"""
: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
......@@ -44,7 +44,7 @@ class Image:
def nan_to_num(self):
"""
convert 'nan' into 0
Converts 'nan' into 0
"""
self.data = np.nan_to_num(self.data)
......@@ -65,8 +65,7 @@ class Image:
def interpolation(self):
"""
Perform an interpolation of the image at higher resolution and stock the new image
Performs an interpolation of the image at higher resolution and stock the new image
"""
x = np.linspace(0, self.length, self.length) # + 0.5*self.oversample
y = np.linspace(0, self.high, self.high) # + 0.5*self.oversample
......@@ -82,7 +81,7 @@ class Image:
def conv_inter_flux(self, psf):
"""
Do a convolution of the interpolated image with a psf and stock the result
Does a convolution of the interpolated image with a psf and stock the result
:param PSF psf: psf object
"""
......
This diff is collapsed.
......@@ -3,8 +3,13 @@
####
gal name: test
#In this section write the name of fits' files from the directory which contain this file
#The program can search file in subdirectories
# maybe add a key nobj to identify the number of galaxies to be fitted
# then for each galaxy, need a flux_HD.fits file
# for each galaxy add a key ncomp to identify the number of components to be used
# for each component, only rotation curve parameters change as well as model: model should be in each component rather than in 'config fit'
#In this section write the name of fits' files from the directory which contains this file
#The program can search files in subdirectories
files:
flux: flux_HD.fits
vel: vel.fits
......@@ -13,59 +18,64 @@ files:
psf:
#This section is to set parameters needed by the model used
init fit:
#Inside 'parname' write the parameters which used by the model in the same order than you declared your model
parname: ['xc', 'yc', 'pa', 'inc', 'vs', 'vm', 'rt']
#For each parameters, set the initial 'value', 'limits' and if is 'fixed' (1) or not (0)
#The order of parameters is no important
xc:
desc: center abscissa in pixel
value: 14.8
limits: [10, 20]
fixed: 1
yc:
desc: center ordinate in pixel
value: 14.9
limits: [10, 20]
fixed: 1
pa:
desc: positon angle in degree
value: 0
limits: [0, 180]
fixed: 0
inc:
desc: inclination in degree
value: 45
limits: [10, 80]
fixed: 1
vs:
desc: systemic velocity in km/s
value: 0
limits: [-10, 10]
fixed: 0
vm:
desc: model velocity in km/s
value: 209
limits: [0, 500]
fixed: 0
rt:
desc: transition radius
value: 4
limits: [1, 12]
fixed: 0
#Inside 'parname' write the parameters which are used by the model in the same order than you declared your model
parname: ['xc', 'yc', 'pa', 'inc', 'vs', 'vm', 'rt']
#For each parameter, set the initial 'value', 'limits' and if is 'fixed' (1) or not (0)
#The order of parameters is no important
xc:
desc: center abscissa in pixel
value: 14.8
limits: [10, 20]
fixed: 1
yc:
desc: center ordinate in pixel
value: 14.9
limits: [10, 20]
fixed: 1
pa:
desc: positon angle in degree
value: 0
limits: [0, 360]
fixed: 0
inc:
desc: inclination in degree
value: 45
limits: [10, 80]
fixed: 1
vs:
desc: systemic velocity in km/s
value: 0
limits: [-10, 10]
fixed: 0
vm:
desc: model velocity in km/s
value: 209
limits: [0, 400]
fixed: 0
rt:
desc: transition radius
value: 4
limits: [1, 12]
fixed: 0
#This section is for parameters of the model which are not fitted
conf model:
sig: 40
slope: 0
psfx: 3.5
psfz: 49
smooth: 0
#Selection the 'model' wanted, the 'method' and others things about the running of the program
sig: 40
slope: 0
psfx: 3.5
psfz: 49
smooth: 0
#Selection of the desired 'model', 'method' and others things about running the program
config fit:
model: flat
method: mpfit
verbose: False
#parameters about PyMultiNest running
PyMultiNest:
nbp: 0
plt stats: False
plt xy: True
model: flat
method: mpfit
verbose: False
#parameters about PyMultiNest running
PyMultiNest:
nbp: 0
plt stats: False
plt xy: True
#parameters about mpfit running
mpfit:
ftol: 1e-10
gtol: 1e-10
xtol: 1e-10
......@@ -2,30 +2,37 @@ import numpy as np
from scipy.special import i0, i1, k0, k1
"""
Librairy which contains elementary rotation curve models
"""
def exponential_velocity(r, rt, vm):
"""
Velocity function for an exponential disk
Rotation curve for an exponential disk
:param ndarray r: 2D array which contain the radius
:param int rt: radius at which the maximum velocity is reached
:param ndarray r: 1D or 2D array which contains the radius
:param float rt: radius at which the maximum velocity is reached
:param float vm: Maximum velocity of the model
"""
rd = rt / 2.15 # disk scale length
rd = rt / 2.15 # disk scale length, maximum is reached at 2.15 rd
vr = np.zeros(np.shape(r))
q = np.where(r != 0) # To prevent any problem in the center
vr[q] = r[q] / rd * vm / 0.88 * np.sqrt(i0(0.5 * r[q] / rd) * k0(0.5 * r[q] / rd) - i1(0.5 * r[q] / rd) * k1(0.5 * r[q] / rd))
vr[q] = r[q] / rd * vm / 0.8798243 * np.sqrt(i0(0.5 * r[q] / rd) * k0(0.5 * r[q] / rd) - i1(0.5 * r[q] / rd) * k1(0.5 * r[q] / rd))
return vr
def flat_velocity(r, rt, vm):
"""
Velocity function for flat disk
Flat rotation curve:
v = r/rt * vm when r<=rt
v = vm when r>=rt
:param ndarray r: 2D array which contain the radius
:param int rt: radius at which the maximum velocity is reached
:param ndarray r: 1D or 2D array which contains the radius
:param float rt: radius at which the maximum velocity is reached
:param float vm: maximum velocity of the model
"""
......@@ -39,14 +46,162 @@ def flat_velocity(r, rt, vm):
def arctan_velocity(r, rt, vm):
"""
:param ndarray r:
:param int rt:
:param float vm:
:return:
Arctangent rotation curve
:param ndarray r: 1D or 2D array which contains the radius
:param float rt: scalelength
:param float vm: maximum velocity of the model
"""
return 2*vm/np.pi*np.arctan(2*r/rt)
# Must be at th end of the file
list_model = {'exp': exponential_velocity, 'flat': flat_velocity, 'arctan': arctan_velocity}
def kent_velocity(r, rt, vm):
"""
Rotation curve for a Kent or pseudo isothermal sphere density profile (e.g. Korsaga et al. 2018)
:param ndarray r: 1D or 2D array which contains the radius
:param float rt: scalelength
:param float vm: maximum velocity of the model
"""
vr = np.zeros(np.shape(r))
q = np.where(r != 0) # To prevent any problem in the center
vr[q] = vm * np.sqrt(1 - rt/r[q] * np.arctan(r[q]/rt))
return vr
def hubblemodified_velocity(r, rt, vm):
"""
Rotation curve for a Hubble modified density profile (see Binney & Tremaine)
:param ndarray r: 1D or 2D array which contains the radius
:param float rt: scalelength
:param float vm: maximum velocity of the model
"""
rd = rt / 2.92 # scalelength, maximum is reached at 2.15 rd
vr = np.zeros(np.shape(r))
q = np.where(r != 0) # To prevent any problem in the center
vr[q] = vm / 0.53851 * np.sqrt(1 / (r[q]/rd) * (np.log(r[q]/rd + np.sqrt(1 + (r[q]/rd)**2)) - r[q]/rd * (1 + (r[q]/rd)**2)**(-0.5)))
return vr
def nfw_velocity(r, rt, vm):
"""
Rotation curve for a Navaro Frenk & White (NFW) density profile (see Binney & Tremaine)
:param ndarray r: 1D or 2D array which contains the radius
:param float rt: scalelength
:param float vm: maximum velocity of the model
"""
rd = rt / 2.163 # scalelength, maximum is reached at 2.163 rd
vr = np.zeros(np.shape(r))
q = np.where(r != 0) # To prevent any problem in the center
vr[q] = vm / 0.46499096 * np.sqrt((np.log(1 + r[q]/rd) - (r[q]/rd) / (1 + r[q]/rd)) / (r[q]/rd))
return vr
#def iso_velocity(r, rt, vm):
#"""
#Rotation curve for a Isothermal sphere density profile (see Binney & Tremaine)
#XXX See function in mass model
#:param ndarray r: 1D or 2D array which contains the radius
#:param float rt: scalelength
#:param float vm: maximum velocity of the model
#"""
#def hernquist_velocity():
#"""
#Rotation curve for a Hernquist density profile (see Binney & Tremaine)
#:param ndarray r: 1D or 2D array which contains the radius
#:param float rt: scalelength
#:param float vm: maximum velocity of the model
#"""
#def jaffe_velocity():
#"""
#Rotation curve for a Jaffe density profile (see Binney & Tremaine)
#:param ndarray r: 1D or 2D array which contains the radius
#:param float rt: scalelength
#:param float vm: maximum velocity of the model
#"""
#def einasto_velocity():
#"""
#Rotation curve for a Einasto density profile
#:param ndarray r: 1D or 2D array which contains the radius
#:param float rt: scalelength
#:param float vm: maximum velocity of the model
#:param float
#"""
#def adiabatic_contraction():
#"""
#Should be applicable to any model and have to include baryon distribution
#"""
#def sersic():
def courteau_velocity(r, rt, vm, a):
"""
Courteau rotation curve
:param ndarray r: 1D or 2D array which contains the radius
:param float rt: scalelength
:param float vm: typical velocity
:param float a: Courteau index
"""
vr = np.zeros(np.shape(r))
q = np.where(r != 0) # To prevent any problem in the center
vr[q] = vm / (1 + (1 / (r[q]/rt))**a) ** (1 / a)
return vr
def zhao_velocity(r, rt, vt, a, b, g):
"""
Zhao rotation curve (see Epinat et al. 2008)
:param ndarray r: 1D or 2D array which contains the radius
:param float rt: scalelength
:param float vt: typical velocity
:param float a: first Zhao index
:param float b: second Zhao index
:param float g: third Zhao index
"""
if a == 0:
exp = 1
else:
exp = b/a
return vt * (r/rt)**g / (1 + (r/rt)**a)**exp
#def nfw_exponential_velocity(r, rt1, vm1, rt2, vm2):
#return np.sqrt(nfw_velocity(r, rt1, vm1) ** 2 + exponential_velocity(r, rt1, vm1) ** 2)
# Must be at the end of the file
list_model = {'exp': exponential_velocity, 'flat': flat_velocity, 'arctan': arctan_velocity, 'hubmod': hubblemodified_velocity, 'nfw': nfw_velocity, 'kent': kent_velocity,
#'iso': iso_velocity, 'hernquist': hernquist_velocity, 'jaffe': jaffe_velocity, 'einasto': einasto_velocity,
'courteau': courteau_velocity, 'zhao': zhao_velocity}
# GalKin: Galaxy Kinematics
# MocKinG: Mocking/Modeling Kinematics of Galaxies
## What is it?
**GalKin** fits a velocity model on an observed velocity field of galaxy.
**MocKinG** fits a velocity model on an observed velocity field of galaxy.
## Model
......@@ -18,22 +18,21 @@ But one can add more rotational curves in "velocity_model.py" or create another
## Installation
The program needs PyMultiNest to be installed, please refer to [this link](https://johannesbuchner.github.io/PyMultiNest/) for
its installation guide.
The program needs PyMultiNest to be installed, please refer to [this link](https://johannesbuchner.github.io/PyMultiNest/) for its installation guide.
The program also needs **astropy** and **yaml** libraries.
## How to launch the program?
**Galkin** can be launch from prompt with
**MocKinG** can be launched from prompt with
```
galkin.py path filename
mocking.py path filename
```
or in a python script/console by importing first **galkin** and then calling the main function
or in a python script/console by importing first **MocKinG** and then calling the main function
```
import galkin.py
galkin.galkin(path, filename, rank=0)
import mocking.py
mocking.mocking(path, filename, rank=0)
```
where path can be absolute or relative, filename is the name of the **YAML** configuration file. The parameter rank is the id of the thread.
......@@ -44,7 +43,7 @@ For example, if you run the program with mpi and 4 core, rank may be equal from
For more information and how to install it, follow [this link](http://pythonhosted.org/mpi4py/).
To execute the program with **mpi4py**:
```
mpiexec -n (nbcore) galkin.py path filename
mpiexec -n (nbcore) mocking.py path filename
```
## Output
......@@ -54,9 +53,13 @@ method_model_fixedparams. A summary of the parameters of the model is written in
## Example
There is a Example directory in which you have fits files and a **YAML** file as example to try the program.
There is an Example directory in which you have fits files and a **YAML** file as example to try the program.
## Create model
The subprocess to create map was used to validate the program (cross test) with existing programs.
This subprocess needs to be updated from the last changes of the architecture and the addition of the data language **YAML**.
### Credits
Written by Jeremy Dumoulin under supervision of Benoît Epinat
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname('../Models/'), '..')))
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname('../Class/'), '..')))
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname('../Tools/'), '..')))
import numpy as np
from numpy.fft import fftshift, fft2, ifft2
import math
......@@ -15,6 +12,9 @@ from Class.PSF import PSF
from Class.Model2D import Model2D
import argparse
from astropy.io import ascii
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname('../Models/'), '..')))
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname('../Class/'), '..')))
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname('../Tools/'), '..')))
def main(parser):
......
......@@ -6,16 +6,19 @@ from Class.Model2D import Model2D
import Tools.cap_mpfit as mpfit
import logging
logger = logging.getLogger('__galkin__')
logger = logging.getLogger('__mocking__')
def use_mpfit(model, params, quiet=False):
def use_mpfit(model, params, quiet=False, gtol=1e-10, ftol=1e-10, xtol=1e-10):
"""
Function call MPFIt to perform the fit of the model
Function call mpfit to perform the fit of the model
:param Model2D model: model initialised
:param dict params: dictionary which contain all parameters with the limits and if are fixed or not
:param Bool quiet: print or not verbose from the fit method
:param float gtol: orthogonality desired between the function vector and the columns of the jacobian.
:param float ftol: relative error desired in the sum of squares.
:param float xtol: relative error desired in the approximate solution.
"""
if quiet is True:
......@@ -24,6 +27,7 @@ def use_mpfit(model, params, quiet=False):
verbose = 1
# create the array of dictionary 'parinfo' needed by mpfit
# xxx eventuellement utiliser la liste définie dans velocity_models, ça permettrait de faire les vérifications + éviter d'avoir la liste de noms
p0 = [params[key] for key in params['parname']]
parinfo = [{'value': 0., 'fixed': 0, 'limited': [1, 1], 'limits': [0., 0.], 'parname': 0., 'step': 0.} for i in range(len(p0))]
......@@ -40,7 +44,7 @@ def use_mpfit(model, params, quiet=False):
t1 = time.time()
# ### Call mpfit
model_fit = mpfit.mpfit(model.least_square, parinfo=parinfo, autoderivative=1, gtol=1e-10, ftol=1e-10, xtol=1e-10, quiet=verbose)
model_fit = mpfit.mpfit(model.least_square, parinfo=parinfo, autoderivative=1, gtol=gtol, ftol=ftol, xtol=xtol, quiet=verbose)
# ###
t2 = time.time()
......
......@@ -10,7 +10,7 @@ import matplotlib.pyplot as plt
import logging
logger = logging.getLogger('__galkin__')
logger = logging.getLogger('__mocking__')
def use_pymultinest(model, params, quiet=True, nbp=0, pltstats=False, rank=0, path=None, whd=None):
......
......@@ -6,12 +6,12 @@ import sys
import yaml
import logging
logger = logging.getLogger('__galkin__')
logger = logging.getLogger('__mocking__')
def sky_coord_to_galactic(xc, yc, pa, incl, im_size=None):
"""
Convert position from Sky coordinates to Galactic coordinates
Converts position from Sky coordinates to Galactic coordinates
:param float xc: position of the center in arcsec
:param float yc: position of the center in arcsec
......@@ -37,7 +37,7 @@ def sky_coord_to_galactic(xc, yc, pa, incl, im_size=None):
def rebin_data(data, factor):
"""
Rebin an image.
Rebins an image.
example: For rebin an image from 240x240 to 60x60 pixels, factor=5
......@@ -81,7 +81,7 @@ def write_fits(data, filename, config, results, mask=None):
def search_file(path, filename):
"""
Search a file in a directory and all its subdirectories. Return the path of the file
Searches a file in a directory and all its subdirectories. Return the path of the file
:param str path: path where to search
:param str filename: name of the file to search (with its extension)
......@@ -110,8 +110,8 @@ def search_file(path, filename):
def make_dir(path, config):
"""
Create the directory where results will be written
The name of the directory depend of fixed parameters
Creates the directory where results will be written
The name of the directory depends on fixed parameters
:param path: path where fits files are
:param dict config: YAML config dictionary
......@@ -141,8 +141,7 @@ def make_dir(path, config):
def write_yaml(path, params, galname, whd):
"""
Setup the stream and write the YAML file which contain the results
Sets up the stream and writes the YAML file which contains the results
:param str path: path where to write the YAML file
:param dict params: dictionary of the bestfit parameters
......
#!/usr/bin/env python3
#!/usr/bin/env python
import argparse
import os
import sys
......@@ -26,35 +26,35 @@ except ImportError:
# For debug use logging.DEBUG instead of logging.INFO
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger('__galkin__')
logger = logging.getLogger('__mocking__')
def galkin(path=None, filename=None, rank=0):
def mocking(path=None, filename=None, rank=0):
"""
(name not found) fit model to velocity field of galaxies.
For more information see the help or refer to the git repository:
MockinG (Mocking/Modeling Kinematics of Galaxies) fits a model to a velocity field of galaxies.
For more information see the help or refer to the git repository:
https://gitlab.lam.fr/bepinat/mocking
https://github.com/Meirdrarel/batman
developed on python 3.6
developed on python 3.6
This program can be launched from console:
mocking.py (path) (filename)
This program can be lunch from console:
or in a python shell/console by importing this file and call:
galkin.py (path) (filename)
mocking.mocking(path, filename, rank)
or in a python shell/console by importing this file and call:
To run with mpi type in the prompt:
galkin(path, filename, rank)
mpiexec -n (nbcore) mocking.py (path) (filename)
For run with mpi type in the prompt:
Outputs are written in the directory where the YAML file is,
the directory name is: (method)_(model)_(paramfix)
mpiexec -n (nbcore) galkin.py (path) (filename)
Output are written in a directory where the YAML file is,
the directory name is: (method)_(model)_(paramfix)
:param str path: path where the config file is
:param str filename: name of the config file
:param str path: path where the YAML configuration file is and where outputs will be written
:param str filename: name of the YAML file which contains all parameters
:param int rank: id of the thread when the program is run with MPI4PY
"""
......@@ -80,7 +80,7 @@ def galkin(path=None, filename=None, rank=0):
disp = Image(tools.search_file(path, files['disp']))
###################################################################################################
# Test the size of the flux image with de velocity field and perform an interpolation if needed
# Test the size of the flux image with the velocity field and perform an interpolation if needed
length_rap = flux.length / vel.get_lenght()
high_rap = flux.high / vel.get_high()
if rank == 0:
......@@ -96,8 +96,8 @@ def galkin(path=None, filename=None, rank=0):
else:
flux.set_oversamp(int(length_rap))
whd = '_whd'
####################################################################################################
# Check psf file existence and import it or create a gaussian
img_psf = None
if files['psf'] is not None:
......@@ -113,18 +113,26 @@ def galkin(path=None, filename=None, rank=0):
# Initialise the model's class
model = Model2D(vel, errvel, flux, psf, vm.list_model[config['config fit']['model']], confmod['sig'], confmod['slope'])
# Retrieve model parameters
# XXX ajout pour trouver automatiquement le nombre de paramètres du modèle