Example Dysmalpy 1D fitting, using fitting wrapper
In this example, we use dysmalpy to measure the kinematics of galaxy GS4_43501 at \(z=1.613\) in 1D, using a fitting wrapper which simplifies the implementation of the fitting algorithm to the user. In this specific case, the fittign method is \(\texttt{MPFIT}\), as specified at the bottom of the fitting_1D_mpfit.params file. The notebook shows how to find the best fit models for the one-dimensional velocity and velocity dispersion profiles (\(v(r)\) and \(\sigma(r)\)). These fits allow us to measure quantities such as the total mass (disk+bulge), the effective radius \(r_\mathrm{eff}\), dark matter fraction \(f_\mathrm{DM}\) and velocity dispersion \(\sigma_0\).
The fitting includes the following components:
Disk + Bulge
NFW halo
Constant velocity dispersion
The structure of the notebook is the following:
Setup steps (and load the params file)
Run
Dysmalpyfitting: 1D wrapper, with fit method= MPFITExamine results
1) Setup steps
First import modules
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from dysmalpy import fitting
from dysmalpy.fitting_wrappers import dysmalpy_fit_single, utils_io
import os
import numpy as np
INFO:numexpr.utils:Note: NumExpr detected 10 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
INFO:numexpr.utils:NumExpr defaulting to 8 threads.
Setup notebook
# Setup plotting
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
mpl.rcParams['figure.dpi']= 300
mpl.rc("savefig", dpi=300)
from IPython.core.display import Image
Set data, output paths
Note this will override the
datadirandoutdirspecified in the param file.(This is useful for the example here. When running from command line, it’s recommended to properly set the directories in the param file.)
# Data directory (datadir = /YOUR/DATA/PATH/)
filepath = os.path.abspath(fitting.__file__)
datadir = os.sep.join(os.path.dirname(filepath).split(os.sep)[:-1]+["tests", "test_data", ""])
# Load the parameters file from the examples directory
param_path = os.sep.join(os.path.dirname(filepath).split(os.sep)[:-1]+["examples", "examples_param_files", ""])
param_filename = param_path+'fitting_1D_mpfit.params'
# Where to save output files (output = /YOUR/OUTPUTS/PATH)
outdir = '/Users/sedona/data/dysmalpy_test_examples/JUPYTER_OUTPUT_1D_FITTING_WRAPPER/'
outdir = '/Users/jespejo/Dropbox/Postdoc/Data/dysmalpy_test_examples/JUPYTER_OUTPUT_1D_FITTING_WRAPPER/'
outdir_mpfit = outdir + 'MPFIT/'
Settings in parameter file:
Note there are many commented out options / parameters. These given an more complete overview of the settings & parameters that can be specified with the fitting wrapper parameter files.
with open(param_filename, 'r') as f:
print(f.read())
# Example parameters file for fitting a single object with 1D data
# Note: DO NOT CHANGE THE NAMES IN THE 1ST COLUMN AND KEEP THE COMMAS!!
# See README for a description of each parameter and its available options.
# ******************************* OBJECT INFO **********************************
galID, GS4_43501 # Name of your object
z, 1.613 # Redshift
# ****************************** DATA INFO *************************************
datadir, None # Optional: Full path to data directory.
fdata, GS4_43501.obs_prof.txt # Full path to your data. Alternatively, just the filename if 'datadir' is set.
data_inst_corr, True # Is the dispersion corrected for
# instrumental broadening?
slit_width, 0.55 # arcsecs
slit_pa, 142. # Degrees from N towards blue
symmetrize_data, False # Symmetrize data before fitting?
######## Data extraction apertures: ########
## Rectangular apertures:
profile1d_type, rect_ap_cube # 1D aperture extraction shape
pix_perp, 5 # Aperture width (e.g., slit width), in number of pixels
# So arcsec width = pix_perp * pixscale (under 'Instrument Setup')
pix_parallel, 3 # Aperture length (e.g., how many pixels along the slit direction), in npix
# ## Circular apertures:
# profile1d_type, circ_ap_cube # 1D aperture extraction shape
# aperture_radius, 0.275 # Circular aperture radius, in ARCSEC. Have used half slit width in past
# -- Eg, aperture diam = slit width
############################################
# ***************************** OUTPUT *****************************************
outdir, GS4_43501_1D_out/ # Full path for output directory
# ***************************** OBSERVATION SETUP ******************************
# Instrument Setup
# ------------------
pixscale, 0.125 # Pixel scale in arcsec/pixel
fov_npix, 37 # Number of pixels on a side of model cube
spec_type, velocity # DON'T CHANGE!
spec_start, -1000. # Starting value for spectral axis // generally don't change
spec_step, 10. # Step size for spectral axis in km/s // generally don't change
nspec, 201 # Number of spectral steps // generally don't change
# LSF Setup
# ---------
use_lsf, True # True/False if using an LSF
sig_inst_res, 51.0 # Instrumental dispersion in km/s
# PSF Setup
# ---------
psf_type, Gaussian # Gaussian, Moffat, or DoubleGaussian
psf_fwhm, 0.55 # PSF FWHM in arcsecs
psf_beta, -99. # Beta parameter for a Moffat PSF
# ## ELLIPTICAL PSF:
# psf_type, Gaussian # Gaussian, Moffat, or DoubleGaussian
# psf_fwhm_major, 0.55 # PSF major axis FWHM in arcsecs
# psf_fwhm_minor, 0.25 # PSF minor axis FWHM in arcsecs
# psf_PA, 0. # PA of PSF major axis, in deg E of N. (0=N, 90=E)
# psf_beta, -99. # Beta parameter for a Moffat PSF
# # DoubleGaussian: settings instead of psf_fwhm
# psf_type, DoubleGaussian
# psf_fwhm1, 0.16 # FWHM of PSF component 1, in arcsecs. SINFONI AO: 0.16
# psf_fwhm2, 0.48 # FWHM of PSF component 1, in arcsecs. SINFONI AO: 0.48
# psf_scale1, 0.368 # Flux scaling (*not* peak height) of component 1. SINFONI AO: 0.368
# psf_scale2, 0.632 # Flux scaling (*not* peak height) of component 2. SINFONI AO: 0.632
# **************************** SETUP MODEL *************************************
# Model Settings
# -------------
# List of components to use: SEPARATE WITH SPACES
## MUST always keep: geometry zheight_gaus
## RECOMMENDED: always keep: disk+bulge const_disp_prof
components_list, disk+bulge const_disp_prof geometry zheight_gaus halo
# possible options:
# disk+bulge, sersic, blackhole
# const_disp_prof, geometry, zheight_gaus, halo,
# radial_flow, uniform_planar_radial_flow, uniform_bar_flow, uniform_wedge_flow,
# unresolved_outflow, biconical_outflow,
# CAUTION: azimuthal_planar_radial_flow, variable_bar_flow, spiral_flow
# List of components that emit light. SEPARATE WITH SPACES
## Current options: disk+bulge / bulge / disk [corresponding to the mass disk+bulge component],
## also: light_sersic, light_gaussian_ring
light_components_list, disk
# NOTE: if a separate light profile (eg light_sersic) is used,
# this MUST be changed to e.g., 'light_components_list, light_sersic'
adiabatic_contract, False # Apply adiabatic contraction?
pressure_support, True # Apply assymmetric drift correction?
noord_flat, True # Apply Noordermeer flattenning?
oversample, 1 # Spatial oversample factor
oversize, 1 # Oversize factor
moment_calc, False # If False, observed profiles fit with GAUSSIANS
zcalc_truncate, True # Truncate in zgal direction when calculating or not
n_wholepix_z_min, 3 # Minimum number of whole pixels in zgal dir, if zcalc_truncate=True
# ********************************************************************************
# DISK + BULGE
# ------------
# Initial Values
total_mass, 11.0 # Total mass of disk and bulge log(Msun)
bt, 0.3 # Bulge-to-Total Ratio
r_eff_disk, 5.0 # Effective radius of disk in kpc
n_disk, 1.0 # Sersic index for disk
invq_disk, 5.0 # disk scale length to zheight ratio for disk
n_bulge, 4.0 # Sersic index for bulge
invq_bulge, 1.0 # disk scale length to zheight ratio for bulge
r_eff_bulge, 1.0 # Effective radius of bulge in kpc
# Fixed? True if its a fixed parameter, False otherwise
total_mass_fixed, False
r_eff_disk_fixed, False
bt_fixed, True
n_disk_fixed, True
r_eff_bulge_fixed, True
n_bulge_fixed, True
# Parameter bounds. Lower and upper bounds
total_mass_bounds, 10.0 13.0
bt_bounds, 0.0 1.0
r_eff_disk_bounds, 0.1 30.0
n_disk_bounds, 1.0 8.0
r_eff_bulge_bounds, 1.0 5.0
n_bulge_bounds, 1.0 8.0
# # ********************************************************************************
# # BLACK HOLE
# # ------------
#
# # Initial Values
# BH_mass, 11. # log(Msun)
#
# # Fixed? True if its a fixed parameter, False otherwise
# BH_mass_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# BH_mass_bounds, 6. 18.
# # ********************************************************************************
# # Separate light profile: (Truncated) Sersic profile
# # ------------
# # Initial values
# L_tot_sersic, 1. # arbitrary units
# lr_eff, 4. # kpc
# lsersic_n, 1. # Sersic index of light profile
# lsersic_rinner, 0. # [kpc] Inner truncation radius of sersic profile. 0 = no truncation
# lsersic_router, inf # [kpc] Outer truncation radius of sersic profile. inf = no truncation
#
# # Fixed? True if its a fixed parameter, False otherwise
# L_tot_sersic_fixed, True
# lr_eff_fixed, False
# lsersic_n_fixed, True
# lsersic_rinner_fixed, True
# lsersic_router_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# L_tot_sersic_bounds, 0. 2. # arbitrary units
# lr_eff_bounds, 0.5 15. # kpc
# lsersic_n_bounds, 0.5 8.
# lsersic_rinner_bounds, 0. 5. # kpc
# lsersic_router_bounds, 4. 20. # kpc
# # ********************************************************************************
# # Separate light profile: Gaussian ring
# # ------------
# # Initial values
# L_tot_gaus_ring, 1. # arbitrary units
# R_peak_gaus_ring, 6. # kpc
# FWHM_gaus_ring, 1. # kpc
#
# # Fixed? True if its a fixed parameter, False otherwise
# L_tot_gaus_ring_fixed, True
# R_peak_gaus_ring_fixed, True
# FWHM_gaus_ring_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# L_tot_gaus_ring_bounds, 0. 2. # arbitrary units
# R_peak_gaus_ring_bounds, 0. 15. # kpc
# FWHM_gaus_ring_bounds, 0.1 10. # kpc
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# DARK MATTER HALO
# ----------------
# Halo type: options: NFW / twopowerhalo / burkert / einasto / dekelzhao
halo_profile_type, NFW
# ** NOTE **: Uncomment the section below corresponding to the selected halo type.
# ********************************************************************************
# NFW halo
# Initial Values
mvirial, 11.5 # Halo virial mass in log(Msun)
halo_conc, 5.0 # Halo concentration parameter
fdm, 0.5 # Dark matter fraction at r_eff_disk
# Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
mvirial_fixed, False
halo_conc_fixed, True
fdm_fixed, False
# Parameter bounds. Lower and upper bounds
mvirial_bounds, 10.0 13.0
halo_conc_bounds, 1.0 20.0
fdm_bounds, 0.0 1.0
# Tie the parameters?
fdm_tied, True # for NFW, fdm_tied=True determines fDM from Mvirial (+baryons)
mvirial_tied, False # for NFW, mvirial_tied=True determines Mvirial from fDM (+baryons)
# ********************************************************************************
# # ********************************************************************************
# # Two-power halo
#
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# alpha, 1. # TPH: inner slope. NFW has alpha=1
# beta, 3. # TPH: outer slope. NFW has beta=3
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# alpha_fixed, False
# beta_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# alpha_bounds, 0.0 3.0
# beta_bounds, 1.0 4.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# alpha_tied, False # for TPH, alpha_tied=True determines alpha from free fDM + other parameters.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Burkert halo
#
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# rB, 10. # Burkert: Halo core radius, in kpc
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# rB_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# rB_bounds, 1.0 20.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# rB_tied, False # for Burkert, rB_tied=True determines rB from free fDM + other parameters.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Einasto halo
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# alphaEinasto, 1. # Einasto: Halo profile index
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# alphaEinasto_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# alphaEinasto_bounds, 0.0 2.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# alphaEinasto_tied, False # for Einasto, alphaEinasto_tied=True determines alphaEinasto from free fDM + other params.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Dekel-Zhao halo
# # Initial Values
# mvirial, 12.0 # Halo virial mass in log(Msun)
# s1, 1.5 # Inner logarithmic slope (at resolution r1=0.01*Rvir)
# c2, 25.0 # Concentration parameter (defined relative to c, a)
# fdm, 0.5 # Dark matter fraction at r_eff_disk
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# s1_fixed, False
# c2_fixed, False
# fdm_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0 # log(Msun)
# s1_bounds, 0.0 2.0
# c2_bounds, 0.0 40.0
# fdm_bounds, 0.0 1.0
#
# # Tie the parameters?
# mvirial_tied, True # mvirial_tied=True determines Mvirial from fDM, s1, c2.
# s1_tied, True # Tie the s1 to M*/Mvir using best-fit Freundlich+20 (Eqs 45, 47, 48, Table 1)
# c2_tied, True # Tie the c2 to M*/Mvir using best-fit Freundlich+20 (Eqs 47, 49, Table 1)
# fdm_tied, False # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
#
# ### OTHER SETTINGS:
# lmstar, 10.5 # Used to infer s1, c2 if s1_tied or c2_tied = True
#
# # ********************************************************************************
# ********************************************************************************
# INTRINSIC DISPERSION PROFILE
# ------------------
# Initial Values
sigma0, 39.0 # Constant intrinsic dispersion value
# Fixed? True if its a fixed parameter, False otherwise
sigma0_fixed, False
# Parameter bounds. Lower and upper bounds
sigma0_bounds, 5.0 300.0
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# HIGHER ORDER COMPONENTS: INFLOW, OUTFLOW
# ----------------------------------------
# # ********************************************************************************
# # UNIFORM SPHERICAL RADIAL FLOW -- in rhat direction in spherical coordinates
# # radial_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# # ********************************************************************************
# # UNIFORM PLANAR RADIAL FLOW -- in Rhat direction in cylindrical coordinates
# # (eg, radial in galaxy midplane)
# # uniform_planar_radial_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# # ********************************************************************************
# # UNIFORM BAR FLOW -- in xhat direction along bar in cartesian coordinates,
# # with bar at an angle relative to galaxy major axis (blue)
# # uniform_bar_flow
# # -------------------
#
# vbar, -90. # Bar flow [km/s]. Positive: Outflow. Negative: Inflow.
# phi, 90. # Azimuthal angle of bar [degrees], counter-clockwise from blue major axis.
# # Default is 90 (eg, along galaxy minor axis)
# bar_width, 2 # Width of the bar perpendicular to bar direction.
# # Bar velocity only is nonzero between -bar_width/2 < ygal < bar_width/2.
# # ********************************************************************************
# # UNIFORM WEDGE FLOW -- in planar radial flow in cylindrical coordinates, restricted to pos, neg wedges
# # uniform_wedge_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# theta, 60. # Opening angle of wedge [deg]. (the full angular span)
# phi, 90. # Angle offset relative to the galaxy angle, so the wedge center is at phi.
# # Default: 90 deg, so centered along minor axis
# # ********************************************************************************
# # UNRESOLVED OUTFLOW -- at galaxy center (ie, AGN unresolved outflow)
# # unresolved_outflow
# # -------------------
#
# vcenter, 0. # Central velocity of the Gaussian in km/s
# fwhm, 1000. # FWHM of the Gaussian in km/s
# amplitude, 1.e12 # Amplitude of the Gaussian, for flux in ~M/L=1 luminosity units
# # with the dimming applied ... roughly ....
# # ********************************************************************************
# # BICONICAL OUTFLOW
# # biconical_outflow
# # -------------------
#
# n, 0.5 # Power law index
# vmax, 500. # Maximum velocity of the outflow in km/s
# rturn, 5. # Turn-over radius in kpc of the velocty profile
# thetain, 30. # Half inner opening angle in degrees. Measured from the bicone axis
# dtheta, 20. # Difference between inner and outer opening angle in degrees
# rend, 10. # Maximum radius of the outflow in kpc
# norm_flux, 8. # Log flux amplitude of the outflow at r = 0.
# # Need to check dimming/flux conventions
# tau_flux, 1. # Exponential decay rate of the flux
# biconical_profile_type, both # Type of velocity profile:
# # 'both', 'increase', 'decrease', 'constant'
# biconical_outflow_dispersion, 80. # Dispersion (stddev of gaussian) of biconical outflow, km/s
# # ********************************************************************************
# # VARIABLE BAR FLOW -- in xhat direction along bar in cartesian coordinates,
# # with bar at an angle relative to galaxy major axis (blue)
# # CAUTION!!!
# # variable_bar_flow
# # -------------------
#
# vbar_func_bar_flow, -90.*np.exp(-R/5.) # Bar flow FUNCTION [km/s]. Positive: Outflow. Negative: Inflow.
# phi, 90. # Azimuthal angle of bar [degrees], counter-clockwise from blue major axis.
# # Default is 90 (eg, along galaxy minor axis)
# bar_width, 2 # Width of the bar perpendicular to bar direction.
# # Bar velocity only is nonzero between -bar_width/2 < ygal < bar_width/2.
# # ********************************************************************************
# # AZIMUTHAL PLANAR RADIAL FLOW -- in Rhat direction in cylindrical coordinates
# # (eg, radial in galaxy midplane), with an added azimuthal term
# # CAUTION!!!
# # azimuthal_planar_radial_flow
# # -------------------
#
# vr_func_azimuthal_planar_flow, -90.*np.exp(-R/5.) # Radial flow [km/s].
# # Positive: Outflow. Negative: Inflow.
# m, 2 # Number of modes in the azimuthal pattern. m=0 gives a purely radial profile.
# phi0, 0. # Angle offset relative to the galaxy angle [deg],
# # so the azimuthal variation goes as cos(m(phi_gal - phi0))
# # ********************************************************************************
# # SPIRAL DENSIY WAVE FLOW -- as in Davies et al. 2009, ApJ, 702, 114
# # Here assuming CONSTANT velocity -- try to match real Vrot...
# # CAUTION!!! NO SPACES IN FUNCTION DEFINITONS!
# # spiral_flow
# # -------------------
#
# Vrot_func_spiral_flow, 150.+0.*R # Unperturbed rotation velocity of the galaxy
# dVrot_dR_func_spiral_flow, 0.*R # Derivative of Vrot(R) -- ideally evaluated analytically, otherwise very slow.
# rho0_func_spiral_flow, 1.e11*np.exp(-R/5.) # Unperturbed midplane density profile of the galaxy
# f_func_spiral_flow, (np.sqrt(m**2-2.)*Vrot(R)/cs)*np.log(R) # Function describing the spiral shape, m*phi = f(R)
# # with k = df/dR
# k_func_spiral_flow, (np.sqrt(m**2-2.)*Vrot(R)/cs)/R # Function for radial wavenumber
#
# m, 2 # Number of photometric/density spiral arms.
# cs, 10. # Sound speed of medium, in km/s.
# epsilon, 1. # Density contrast of perturbation (unitless).
# Om_p, 0. # Angular speed of the driving force, Omega_p
# phi0, 0. # Angle offset of the arm winding, in degrees. Default: 0.
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# ZHEIGHT PROFILE
# ---------------
# Initial Values
sigmaz, 0.9 # Gaussian width of the galaxy in z, in kpc
# Fixed? True if its a fixed parameter, False otherwise
sigmaz_fixed, False
# Parameter bounds. Lower and upper bounds
sigmaz_bounds, 0.1 1.0
# Tie the zheight to the effective radius of the disk?
# If set to True, make sure sigmaz_fixed is False
zheight_tied, True
# GEOMETRY
# --------
# Initial Values
inc, 62. # Inclination of galaxy, 0=face-on, 90=edge-on
# Fixed? True if its a fixed parameter, False otherwise
inc_fixed, True
# Parameter bounds. Lower and upper bounds
inc_bounds, 0.0 90.0
# **************************** Fitting Settings ********************************
fit_method, mpfit # mcmc, nested, or mpfit
do_plotting, True # Produce all output plots?
fitdispersion, True # Simultaneously fit the velocity and dispersion?
fitflux, False # Also fit for the flux?
# MPFIT Settings
#----------------
maxiter, 200 # Maximum number of iterations before mpfit quits
Add some settings for this notebook example:
plot_type = 'png'
2) Run Dysmalpy fitting: 1D wrapper, with fit method= MPFIT
dysmalpy_fit_single.dysmalpy_fit_single(param_filename=param_filename,
datadir=datadir, outdir=outdir_mpfit,
plot_type=plot_type, overwrite=True)
INFO:DysmalPy:*************************************
INFO:DysmalPy: Fitting: GS4_43501 using MPFIT
INFO:DysmalPy: obs: OBS
INFO:DysmalPy: velocity file: /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501.obs_prof.txt
INFO:DysmalPy: nSubpixels: 1
INFO:DysmalPy: mvirial_tied: False
INFO:DysmalPy:
MPFIT Fitting:
Start: 2024-01-03 10:10:44.647732
INFO:DysmalPy:Iter 1 CHI-SQUARE = 82.00732713 DOF = 32
disk+bulge:total_mass = 11
disk+bulge:r_eff_disk = 5
halo:mvirial = 11.5
dispprof_LINE:sigma0 = 39
INFO:DysmalPy:Iter 2 CHI-SQUARE = 71.78611963 DOF = 32
disk+bulge:total_mass = 10.98404781
disk+bulge:r_eff_disk = 5.431735616
halo:mvirial = 11.69606647
dispprof_LINE:sigma0 = 40.2897216
INFO:DysmalPy:Iter 3 CHI-SQUARE = 68.0789783 DOF = 32
disk+bulge:total_mass = 10.93780574
disk+bulge:r_eff_disk = 5.645223846
halo:mvirial = 12.13248856
dispprof_LINE:sigma0 = 43.27809695
INFO:DysmalPy:Iter 4 CHI-SQUARE = 64.55272712 DOF = 32
disk+bulge:total_mass = 10.87549154
disk+bulge:r_eff_disk = 5.025363644
halo:mvirial = 12.21393822
dispprof_LINE:sigma0 = 43.71883553
INFO:DysmalPy:Iter 5 CHI-SQUARE = 60.90507325 DOF = 32
disk+bulge:total_mass = 10.76106807
disk+bulge:r_eff_disk = 3.544651118
halo:mvirial = 12.39674668
dispprof_LINE:sigma0 = 42.85289749
INFO:DysmalPy:Iter 6 CHI-SQUARE = 58.86424672 DOF = 32
disk+bulge:total_mass = 10.66780383
disk+bulge:r_eff_disk = 2.712423548
halo:mvirial = 12.54275972
dispprof_LINE:sigma0 = 39.47363979
INFO:DysmalPy:Iter 7 CHI-SQUARE = 58.59151209 DOF = 32
disk+bulge:total_mass = 10.67734926
disk+bulge:r_eff_disk = 2.762632506
halo:mvirial = 12.49244494
dispprof_LINE:sigma0 = 38.10677514
INFO:DysmalPy:Iter 8 CHI-SQUARE = 58.59027979 DOF = 32
disk+bulge:total_mass = 10.67811434
disk+bulge:r_eff_disk = 2.772714672
halo:mvirial = 12.49263336
dispprof_LINE:sigma0 = 38.18435198
INFO:DysmalPy:Iter 9 CHI-SQUARE = 58.59022498 DOF = 32
disk+bulge:total_mass = 10.67826899
disk+bulge:r_eff_disk = 2.775022027
halo:mvirial = 12.4926511
dispprof_LINE:sigma0 = 38.2000539
INFO:DysmalPy:Iter 10 CHI-SQUARE = 58.59022442 DOF = 32
disk+bulge:total_mass = 10.67828168
disk+bulge:r_eff_disk = 2.775243281
halo:mvirial = 12.49266183
dispprof_LINE:sigma0 = 38.20157538
INFO:DysmalPy:Iter 11 CHI-SQUARE = 58.59022439 DOF = 32
disk+bulge:total_mass = 10.67828485
disk+bulge:r_eff_disk = 2.775292025
halo:mvirial = 12.49266282
dispprof_LINE:sigma0 = 38.201916
INFO:DysmalPy:
End: 2024-01-03 10:10:47.219127
******************
Time= 2.57 (sec), 0:2.57 (m:s)
MPFIT Status = 1
MPFIT Error/Warning Message = None
******************
------------------------------------------------------------------
Dysmalpy MPFIT fitting complete for: GS4_43501
output folder: /Users/jespejo/Dropbox/Postdoc/Data/dysmalpy_test_examples/JUPYTER_OUTPUT_1D_FITTING_WRAPPER/MPFIT/
------------------------------------------------------------------
3) Examine results
Result plots
Read in parameter file
params = utils_io.read_fitting_params(fname=param_filename)
# Override data + output paths:
params['datadir'] = datadir
params['outdir'] = outdir_mpfit
# Add the plot type:
params['plot_type'] = plot_type
f_galmodel = params['outdir'] + '{}_model.pickle'.format(params['galID'])
f_results = params['outdir'] + '{}_{}_results.pickle'.format(params['galID'],
params['fit_method'])
Best-fit plot
# Look at best-fit saved plot:
filepath = outdir_mpfit+"{}_{}_bestfit_{}.{}".format(params['galID'],
params['fit_method'],
params['obs_1_name'],
params['plot_type'])
Image(filename=filepath, width=600)
Directly generating result plots
Reload the galaxy, results files:
gal, results = fitting.reload_all_fitting(filename_galmodel=f_galmodel,
filename_results=f_results,
fit_method=params['fit_method'])
Plot the best-fit results:
results.plot_results(gal)
Results reports
We now look at the results reports, which include the best-fit values and uncertainties (as well as other fitting settings and output).
# Print report
print(results.results_report(gal=gal))
###############################
Fitting for GS4_43501
Date: 2024-01-03 10:10:48.930775
obs: OBS
Datafiles:
vel : /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501.obs_prof.txt
apertures: RectApertures
fit_velocity: True
fit_dispersion: True
fit_flux: False
moment: False
partial_weight: True
zcalc_truncate: True
n_wholepix_z_min: 3
oversample: 1
oversize: 1
Fitting method: MPFIT
fit status: 1
pressure_support: True
pressure_support_type: 1
###############################
Fitting results
-----------
disk+bulge
total_mass 10.6783 +/- 0.0400
r_eff_disk 2.7753 +/- 0.3250
n_disk 1.0000 [FIXED]
r_eff_bulge 1.0000 [FIXED]
n_bulge 4.0000 [FIXED]
bt 0.3000 [FIXED]
mass_to_light 1.0000 [FIXED]
noord_flat True
-----------
halo
mvirial 12.4927 +/- 0.1405
fdm 0.2684 [TIED]
conc 5.0000 [FIXED]
-----------
dispprof_LINE
sigma0 38.2020 +/- 4.4594
-----------
zheightgaus
sigmaz 0.4714 [TIED]
-----------
geom_1
inc 62.0000 [FIXED]
pa 142.0000 [FIXED]
xshift 0.0000 [FIXED]
yshift 0.0000 [FIXED]
vel_shift 0.0000 [FIXED]
-----------
Adiabatic contraction: False
-----------
Red. chisq: 1.8309
To directly save the results report to a file, we can use the following:
# Save report to file:
f_report = params['outdir'] + '{}_fit_report.txt'.format(params['galID'])
results.results_report(gal=gal, filename=f_report)
Also note the fitting wrappers automatically save two versions of the report files:
fbase = '{}_{}_bestfit_results'.format(params['galID'], params['fit_method'])
f_report_pretty = params['outdir'] + fbase + '_report.info'
f_report_machine = params['outdir'] + fbase + '.dat'
The “pretty” version, automatically saved as *_best_fit_results_report.info, is formatted to be human-readable, and includes more information on the fit settings at the beginning (for reference).
with open(f_report_pretty, 'r') as f:
lines = [line.rstrip() for line in f]
for line in lines: print(line)
###############################
Fitting for GS4_43501
Date: 2024-01-03 10:10:47.831461
obs: OBS
Datafiles:
vel : /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501.obs_prof.txt
apertures: RectApertures
fit_velocity: True
fit_dispersion: True
fit_flux: False
moment: False
partial_weight: True
zcalc_truncate: True
n_wholepix_z_min: 3
oversample: 1
oversize: 1
Fitting method: MPFIT
fit status: 1
pressure_support: True
pressure_support_type: 1
###############################
Fitting results
-----------
disk+bulge
total_mass 10.6783 +/- 0.0400
r_eff_disk 2.7753 +/- 0.3250
n_disk 1.0000 [FIXED]
r_eff_bulge 1.0000 [FIXED]
n_bulge 4.0000 [FIXED]
bt 0.3000 [FIXED]
mass_to_light 1.0000 [FIXED]
noord_flat True
-----------
halo
mvirial 12.4927 +/- 0.1405
fdm 0.2684 [TIED]
conc 5.0000 [FIXED]
-----------
dispprof_LINE
sigma0 38.2020 +/- 4.4594
-----------
zheightgaus
sigmaz 0.4714 [TIED]
-----------
geom_1
inc 62.0000 [FIXED]
pa 142.0000 [FIXED]
xshift 0.0000 [FIXED]
yshift 0.0000 [FIXED]
vel_shift 0.0000 [FIXED]
-----------
Adiabatic contraction: False
-----------
Red. chisq: 1.8309
The “machine” version, automatically saved as *_best_fit_results.dat, is formatted as a machine-readable space-separated ascii file. It includes key parameter fit information, as well as the best-fit reduced chisq.
with open(f_report_machine, 'r') as f:
lines = [line.rstrip() for line in f]
for line in lines: print(line)
# component param_name fixed best_value l68_err u68_err
disk+bulge total_mass False 10.6783 0.0400 0.0400
disk+bulge r_eff_disk False 2.7753 0.3250 0.3250
disk+bulge n_disk True 1.0000 -99.0000 -99.0000
disk+bulge r_eff_bulge True 1.0000 -99.0000 -99.0000
disk+bulge n_bulge True 4.0000 -99.0000 -99.0000
disk+bulge bt True 0.3000 -99.0000 -99.0000
disk+bulge mass_to_light True 1.0000 -99.0000 -99.0000
halo mvirial False 12.4927 0.1405 0.1405
halo fdm TIED 0.2684 -99.0000 -99.0000
halo conc True 5.0000 -99.0000 -99.0000
dispprof_LINE sigma0 False 38.2020 4.4594 4.4594
zheightgaus sigmaz TIED 0.4714 -99.0000 -99.0000
geom_1 inc True 62.0000 -99.0000 -99.0000
geom_1 pa True 142.0000 -99.0000 -99.0000
geom_1 xshift True 0.0000 -99.0000 -99.0000
geom_1 yshift True 0.0000 -99.0000 -99.0000
geom_1 vel_shift True 0.0000 -99.0000 -99.0000
mvirial ----- ----- 12.4927 -99.0000 -99.0000
fit_status ----- ----- 1 -99.0000 -99.0000
adiab_contr ----- ----- False -99.0000 -99.0000
redchisq ----- ----- 1.8309 -99.0000 -99.0000
noord_flat ----- ----- True -99.0000 -99.0000
pressure_support ----- ----- True -99.0000 -99.0000
pressure_support_type ----- ----- 1 -99.0000 -99.0000
obs:OBS:apertures ----- ----- RectApertures -99.0000 -99.0000
obs:OBS:moment ----- ----- False -99.0000 -99.0000
obs:OBS:partial_weight ----- ----- True -99.0000 -99.0000