Profiling galaxies¶
Table of contents¶
Finding the effective Sersic profile given a Bulge+disc profile
Mathematics of profile functions
## Import the necessary modules
import numpy as np
import scipy, scipy.optimize
import time
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
Sersic profile¶
The radial light profile is given by
\begin{align} I_n(r) &= I_0 \exp{\left[-\left(\frac{r}{r_{s}}\right)^{1/n}\right]} \\ &= I_{1/2}\exp{\left(-b_n\left[(\frac{r}{r_n})^{1/n}-1\right]\right)}, \end{align}where $r_s$ is the scalelength, related to $r_n$, the half-light radius as $r_s = \frac{r_n}{{b_n}^n}$
$I_0$ is the central peak intensity, and $I_{1/2} = I_0 e^{-b_n}$
$n$ is the Sersic index, $b_n$ is an $n$-dependent function.
The total flux, where integrated upto a radius $R$ is given by
$F( In the first parametrisation,
$F( In the second parametrisation,
$F( The total flux $F = \lim_{R\rightarrow \infty} F( In terms of the total flux, we can write the Sersic profile as $I_n(r) = \frac{F}{2\pi n r_s^2 \Gamma(2n)} \exp{\left[-\left(\frac{r}{r_{s}}\right)^{1/n}\right]}$ The amplitude of the intensity profile expressed in terms of total flux is more robust than expressing it as a peak intensity or intensity along an isophote. By the definition of the half-light radius,
$F( The function $b_n$ is defined by the relation: $\Gamma_\text{inc}(2n,b_n) = \frac{1}{2}\Gamma(2n)$. There are several excellent approximations to $b_n$. It is largely linear in the range of values of $n$ we are interested in. Capaciolli (1989) provided one such approximation as $b_n \approx 1.992n - 0.3271$ for $0.5 < n < 10$.
Ciotti & Bertin (1999) provided an asympotic form for $b_n$, which is what we use here. The approximations are compared against the actual values below.Finding $b_n$¶
References¶
## Fast approximation for b_n
def bn_capiciolli(n):
""" Linear approximation for b_n"""
return 1.9992*n - 0.3271
def bn_as4(n):
""" Asymptotic form for b_n"""
return 2*n - 1./3 + (4./405)*(n**(-1)) + (46./25515)*(n**(-2)) + (131./1148175)*(n**(-3)) - (2194697./30690717750)*(n**(-4))
## Approximating b_n as a function of n
n = np.arange(0.5,6.2,0.1)
bn = [ ]
for nn in n:
bn_func = lambda x: (scipy.special.gammainc(2.*nn,x)-0.5)**2
bn_star = scipy.optimize.minimize_scalar(bn_func,tol=1e-30).x
bn.append( bn_star )
fig, ax_array = plt.subplots(2,sharex=True)
ax = ax_array[0]
ax.plot(n,bn,'k.')
ax.plot(n,bn_capiciolli(n),'r--',label='Linear') ## Capiciolli (1989) approximation
ax.plot(n,bn_as4(n),'g-',label='Asymptotic')
ax.set_ylabel(r'$b_n$')
ax = ax_array[1]
ax.plot(n,np.abs(bn-bn_capiciolli(n)),'r--')
ax.plot(n,np.abs(bn-bn_as4(n)),'g-')
ax.set_yscale('log')
ax.set_ylabel('| Residual |')
ax.set_xlabel(r'Sersic $n$')
_lgnd = ax_array[0].legend(loc='best')
## Define Sersic profiles
def sersic(r,n,h):
## h is the scalelength
return np.exp(-(1.*r/h)**(1./n))/(2*np.pi*n*h*h*scipy.special.gamma(2.*n))
def sersic_hlr(r,n,hlr, exact=False):
## hlr is the half-light radius
## Exact calculation
if exact:
bn_func = lambda x: (scipy.special.gammainc(2.*n,x)-0.5)**2
bn = scipy.optimize.minimize(bn_func,x0=4.675).x[0]
else:
## Excellent approximation
bn = bn_as4(n)
return np.exp(-bn*((r/hlr)**(1./n)))*(bn**(2*n))/(2*np.pi*n*hlr*hlr*scipy.special.gamma(2.*n))
How do Sersic profiles look?¶
## Visualising radial profiles
fig, ax_array = plt.subplots(1,2)
fig.set_size_inches(12,6)
X = np.linspace(0,10,1001)
h = 1 ## scale radius
for n in xrange(1,5):
## Flux normalised
ax = ax_array[1]
Y = sersic(X,n,h)
ax.plot(X,Y,label=r'$n={0}$'.format(n))
## Peak-intensity normalised
ax = ax_array[0]
ax.plot(X,Y*(2*np.pi*n*h*h*scipy.special.gamma(2.*n)))
ax_array[0].axvline(x=h,color='k',ls=':')
ax_array[0].set_xlabel(r'$x = r/r_s$'); ax_array[1].set_xlabel(r'$x = r/r_s$')
ax_array[0].set_ylabel(r'$I_n(r)$'); ax_array[1].set_ylabel(r'$I_n(r)$')
ax_array[0].set_xscale('log')
ax_array[1].set_yscale('log')
_lgnd = ax_array[1].legend(loc='best')
How concentrated is the flux?¶
## Visualise F(<R)
fig, ax = plt.subplots()
X = np.logspace(-2,6,19) ## R/r_s
for n in xrange(1,5):
Yn = [scipy.special.gammainc(2.*n,x**(1./n)) for x in X]
ax.plot(X,Yn,label=r'$n={0}$'.format(n))
ax.set_xscale('log')
ax.set_ylabel(r'$F\, (<x)$',rotation=90)
ax.set_xlabel(r'$x=r/r_s$')
_lgnd = ax.legend(loc='best').draggable()
The flux in the bulge $(n=4)$ component is spread out to very very large radii. For the disk component, $F(
Bulge+Disc profile¶
The bulge profile is given a sum of a Devaucolers profile ($n=4$) and an Exponential disc ($n=1$). Thus, for the bulge+disc profile with bulge fraction $f$,
\begin{equation} I_\text{bd}(r) \propto \frac{f}{8\pi r_b^2 \Gamma(8)} \exp{\left[-\left(\frac{r}{r_{b}}\right)^{1/4}\right]} + \frac{1-f}{2\pi r_d^2 \Gamma(2)} \exp{\left[-\left(\frac{r}{r_{d}}\right)\right]}, \end{equation}where the constant of proportionality is the total flux of the profile.
The scale radius $r_d$ of an exponential disc is related to its half-light radius $r_1$ as $r_1 = r_d \times b_1$, where $b_1 = 1.678$. Similarly, the scale radius $r_b$ of a Devaucolers bulge is related to its half-light radius $r_4$ as $r_4 = r_b \times b_4^4$, where $b_4 = 7.669$.
In lensfit, the half-light radius of the bulge component is set equal to the exponential scalelength of the disk component, which we set to 1 in this exercise. So, the corresponding half-light radius of the bulge component is 1.678 and the scale radius of the bulge component is $1 / b_4^4 \approx 0.00029$. For numerical stability reasons, we shall deal only with half-light radii which has smaller dynamic range compared to the scaleradii.
def bulge(r,h):
return sersic(r,4.,h)
def bulge_hlr(r,hlr):
return sersic_hlr(r,4.,hlr)
def bulge_lf(r,hlr):
#b4 = 7.66924835630383
#h = hlr/((b4)**4)
#return np.exp(-b4*((r/hlr)**(0.25)))/(8*np.pi*h*h*scipy.special.gamma(8))
return sersic_hlr(r,4.,hlr)
def disk(r,h):
return sersic(r,1.,h)
def bulgedisk(r,f,k):
b = bulge(r,k)
d = disk(r,1.)
return f*b+(1.-f)*d
def bulgedisk_lf(r,f):
b = bulge_lf(r,1.)
d = disk(r,1.)
return f*b + (1.-f)*d
Finding the effective Sersic profile given a Bulge+disc profile¶
Scale the flux so that $F=1$ and scale the coordinates such that $r_d = 1$. Let $r_b = \kappa r_d = \kappa$. One "intuitive" way to define an effective Sersic profile is to a least squares method:
\begin{align} n^*, r_s^* &= \text{argmin}_{n,r_s} \int_0^\infty \mathrm{d}r\, r\, (I_n(r) - I_\text{bd}(r))^2. \;\; [incorrect] \end{align}Note the Jacobian factor $r$ (a Lebesgue measure here), since the profile is actually 2 dimensional.
For a fixed $\kappa$, $n^*$ and $r_s^*$ are a function of $f$
This is quite incorrect! Such an objective function is valid when $I_n(r)$ is a model that we fit to a 'noisy data' $I_{bd}$ that we assume comes from the model with an additive Gaussian noise. However, we are evaluating the closeness of two different models here, and that needs to be defined rigorously in the space of models.
Mathematics of profile functions¶
The collection of generalised intensity profiles (without the requirement of positive amplitude), along with a field of real numbers, can be thought of a vector space of functions. The $\ell 1$ norm of the vectors is then defined as \begin{equation} || v || := \int_0^\infty \mathrm d\, r\, r |I(r)| = F, \end{equation} which has the interpretation of the total flux when $I(r) \ge 0 \, \forall\, r$. The $\ell 1$ ball $B_1(F)$ is then the collection of all generalised profiles which have norm $F$. The $\ell 1$ ball is convex due to linearity (unlike the $\ell 2$ ball). It suffices to restrict our attention to $B_1(1)$ since other vectors can be obtained by scaling.
Sersic profiles form a continuous subset of vectors, lying on a curve in the space. Without loss of generality, we will consider Sersic profiles with normalised flux. In particular, $I_1$ and $I_4$ are two points on $B_1(1)$. Bulge+disc profiles lie on the line joining $I_1$ and $I_4$, parametrised by the bulge fraction $f$. While $B_1(1)$ is convex, the set of Sersic profiles is not and for $0 < f < 1$, $I_\text{bd}$ is outside this subset.
The task at hand to find the 'closest' vector from the subset of Sersic profiles, given an interior point on the line joining $I_1$ and $I_4$. To define what 'closest' means, we first need to introduce a metric in this space of functions. The induced metric in the space as radial intensity profile functions, based on the norm defined as above is then \begin{equation} d(u,v) = \int_0^\infty \mathrm d\, r\, r |u-v|. \end{equation}
Based on the above metric,
\begin{align} n^*, r_s^* &= \text{argmin}_{n,r_s} d(I_n,I_\text{bd}) \equiv || I_n - I_\text{bd} || \\ &= \text{argmin}_{n,r_s} \int_0^\infty \mathrm{d}r\, r\, \left|(I_n(r) - I_\text{bd}(r)) \right|. \end{align}## Set parameters
flux = 1.
kappa = 0.0003 # 1.4 ## dummy
disk_scalelength = 1.
bulge_scalength = kappa
def integrand(r,n,h,f):
ignd = r*(bulgedisk(r,f,kappa)-sersic(r,n,h))**2
return ignd
def lf_integrand(r,n,h,f):
ignd = r*(bulgedisk_lf(r,f) - sersic(r,n,h))**2
return ignd
def hlr_integrand(r,n,hlr,f):
#ignd = r*(bulgedisk_lf(r,f) - sersic_hlr(r,n,hlr))**2
ignd = r*np.abs(bulgedisk_lf(r,f) - sersic_hlr(r,n,hlr))
return ignd
def logintegrand(r,n,h,f):
ignd = np.exp(r/h)*(bulgedisk(r,f,kappa)-sersic(r,n,h))**2
return ignd
Solver options¶
Type of solver. Should be one of
- 'Nelder-Mead' :ref:`(see here) <optimize.minimize-neldermead>`
- 'Powell' :ref:`(see here) <optimize.minimize-powell>`
- 'CG' :ref:`(see here) <optimize.minimize-cg>`
- 'BFGS' :ref:`(see here) <optimize.minimize-bfgs>`
- 'Newton-CG' :ref:`(see here) <optimize.minimize-newtoncg>`
- 'L-BFGS-B' :ref:`(see here) <optimize.minimize-lbfgsb>`
- 'TNC' :ref:`(see here) <optimize.minimize-tnc>`
- 'COBYLA' :ref:`(see here) <optimize.minimize-cobyla>`
- 'SLSQP' :ref:`(see here) <optimize.minimize-slsqp>`
- 'dogleg' :ref:`(see here) <optimize.minimize-dogleg>`
- 'trust-ncg' :ref:`(see here) <optimize.minimize-trustncg>`
- custom - a callable object (added in version 0.14.0),tol = 1e-9
methods = ['Nelder-Mead','Powell','CG','BFGS','Newton-CG','L-BFGS-B','TNC','COBYLA','SLSQP','dogleg','trust-ncg']
method = methods[0] ## BFGS by default
Optimisation¶
btt_range = np.linspace(0.,1.,11)
nstars = []
hstars = []
hlrstars = []
btts = []
solns = []
t1 = time.time()
for f in btt_range:
objective = lambda nh: scipy.integrate.quad(hlr_integrand, 0, np.inf, args=(nh[0],nh[1],f))[0]
#soln = scipy.optimize.minimize(objective, x0=np.array([1.+3.*f, 1.+(kappa-1.)*f]),bounds=[(1.,1.),(4.,kappa)],
# method=method,tol=tol) # for scale radius
#soln = scipy.optimize.minimize(objective, x0=np.array([1.+3.*f, 1.+(kappa-1.)*f]),bounds=[(1.,1.),(4.,kappa)],
# method=method,tol=tol)
#nstar, hstar = soln.x
soln = scipy.optimize.minimize(objective, x0=np.array([1.+3.*f, 1.5]), method=method,tol=tol) # for half-light radius
nstar, hlrstar = soln.x
solns.append(soln)
if soln.success:
nstars.append(nstar)
#hstars.append(hstar)
hlrstars.append(hlrstar)
btts.append(f)
#print f, nstar, hlrstar
else:
print f, soln.message
t2 = time.time()
print "Optimisation done in ", np.round(t2-t1,2), " seconds."
btts_array = np.array(btts)
Effective Sersic parameters as a function of BTT¶
# Plot
fig, ax = plt.subplots(1,2)
fig.set_size_inches(15,6)
ax[0].plot(btts,nstars,'r.-')
ax[0].plot(btts,1+3*btts_array,'r--')
ax[0].axhspan(ymin=1.0,ymax=4.0,color='k',alpha=0.2)
ax[0].set_ylabel(r'$n_{eff}$')
ax[0].set_xlabel(r'B/T')
ax[1].plot(btts,hlrstars,'b.-')
ax[1].plot(btts_array,1.678+(1-1.678)*btts_array,'b--')
ax[1].axhspan(ymin=1.0,ymax=1.678,color='k',ls='--',alpha=0.2)
ax[1].set_ylabel('half-light radius')
ax[1].set_xlabel('B/T')
### Saving the output
output_array = np.vstack([btts, nstars, hlrstars]).T
np.savetxt('/disks/shear14/arunkannawadi/imsim/effective_sersic_l1.txt',output_array,header="btt\t n\t hlr", delimiter='\t',fmt=['%.2f', '%.4f','%.4f'])
Visualise the fits¶
## Plot the profiles
ncols = 2
nrows = int(np.ceil(len(btt_range)/2.))
fig, ax_array = plt.subplots(nrows,ncols)
fig.set_size_inches(12,6*nrows)
x = np.arange(0.,2.,0.001)
for f_index in xrange(len(btt_range)):
ax = ax_array[f_index/ncols,f_index%ncols]
f = btt_range[f_index]
## Draw the Sersic profile
soln = solns[f_index]
if soln.success:
nstar = nstars[f_index]
hlrstar = hlrstars[f_index]
f0 = [sersic_hlr(r,nstar,hlrstar) for r in x]
ax.plot(x,f0,'b-',lw=2,label='Effective Sersic')
## Draw the bulge+disk profiles
f1b = [f*bulge_lf(r,1.) for r in x]
f1d = [(1.-f)*disk(r,1.) for r in x]
f1 = [bulgedisk_lf(r,f) for r in x]
ax.plot(x,f1b,'m--',label='Bulge component')
ax.plot(x,f1d,'c-.',label='Disk component')
ax.plot(x,f1,'r-',label='B+D')
ax.set_title('Bulge fraction: '+str(f))
ax.set_yscale('log')
ax.set_xscale('log')
_lgnd = ax.legend(loc='best').draggable()
## Visualise the goodness of fits
fig, ax = plt.subplots()
ax.axhline(0.0,color='g',ls=':')
for f_index in xrange(len(btt_range)):
f = btt_range[f_index]
soln = solns[f_index]
if soln.success:
dist_from_bulge = scipy.integrate.quad(hlr_integrand, 0, np.inf, args=(4.,1.,f))[0]
dist_from_disc = scipy.integrate.quad(hlr_integrand, 0, np.inf, args=(1.,1.678,f))[0]
min_dist = min(dist_from_bulge, dist_from_disc)
ax.plot(f, dist_from_bulge, 'b.')
ax.plot(f, dist_from_disc, 'r.')
ax.plot(f,soln.fun,'k.')
_ttl = fig.suptitle(r'$\chi^2$')
#ax.set_yscale('log')
ax.set_xlabel('B/T')
## Tests
## Test fluxes
def calculate_flux(profile, n,h, R=np.inf):
def rprof(r,n,h):
ignd = 2*np.pi*r*sersic(r,n,h)
return ignd
flux = scipy.integrate.quad(rprof, 0, R, args=(n,h))
return flux
calculate_flux(sersic, 0.5 ,4.)
ii = 2
rr = 1
solns[ii], scipy.integrate.quad(integrand, a=0,b=np.inf, args=(nstars[ii],hstars[ii],btts[ii]))
print integrand(1,nstars[ii],hstars[ii],btts[ii])
print btts[ii], bulgedisk(rr,btts[ii],kappa), sersic(rr,nstars[ii],hstars[ii])
np.array([scipy.special.gammainc(2.*n,x**(1./n)) for x in X])#/scipy.special.gamma(2*n)
np.exp(-1.0*((0.0)**(1./4.)))
b4_func = lambda x: (scipy.special.gammainc(2.*4,x) - 0.5)**2
soln = scipy.optimize.minimize(b4_func,x0=7.675)
b4 = soln.x[0]
rb = 1./(b4**4)
print b4, rb
print b4**8/(8*np.pi*scipy.special.gamma(8)), 1./(2.*np.pi*scipy.special.gamma(2.))
print sersic_hlr(0,4,1.), sersic(0,1,1.)
print bulge_lf(0,1.), disk(0,1)
print bulgedisk_lf(100.,1.), bulgedisk_lf(100.,0.)
fig, ax = plt.subplots()
ax.plot(X,bulgedisk_lf(X,1.),'r-')
ax.plot(X,bulgedisk_lf(X,0.),'b--')
ax.set_xscale('log')
print solns[0]
print " --- "
print scipy.integrate.quad(tmp_integrand, 0, np.inf, args=(1.,7.669,0.))[0]
print scipy.integrate.quad(tmp_integrand, 0, np.inf, args=(1.,1.678,0.))[0]
def tmp_integrand(r,n,hlr,f):
ignd = r*(sersic(r,1.,1.) - sersic_hlr(r,n,hlr))**2
return ignd
ff = 0.1
print scipy.integrate.quad(hlr_integrand, 0, np.inf, args=(4,1.,ff))[0] + scipy.integrate.quad(hlr_integrand, 0, np.inf, args=(1,1.678,ff))[0]