📣 An issue occured with the embedded container registry on October 25 2021, between 10:30 and 12:10 (UTC+2). Any persisting issues should be reported to CC-IN2P3 Support. 🐛

Commit fd22725d authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

No commit message

No commit message
parent 38f6a989
""" sz catalog evaluation tools
"""
## Copyright (C) 2009 APC CNRS Universite Paris Diderot <lejeune@apc.univ-paris7.fr>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see http://www.gnu.org/licenses/gpl.html
from pylab import *
import numpy
import re
import pyfits
from fitsfunc2 import *
#from kapteyn import wcs
import os as os
from filename import *
import glob
import math
def sextractor (atlasdir, catalog, thresh=5, assocdist=10, regexp="atlas*.fits"):
""" Run Sextractor on patches and merge catalogs.
Parameters
----------
atlasdir: location of the FITS files
catalog: name of the merged catalog
thresh: detection threshold (number of sigmas)
assocdist: distance of association in degrees
Return
------
The output catalog contains [glon glat flux map_index x_coord y_coord]
"""
## build template file
template = os.tmpnam()
os.system("echo X_WORLD > "+template+";echo Y_WORLD >> "+template+";echo FLUX_AUTO >> "+template+";echo Y_IMAGE >> "+template+";echo X_IMAGE >> "+template+"\n")
## build parameter file
parfile = os.tmpnam()
tmpargs = args2paramfile({"CATALOG_TYPE":"ASCII_HEAD","PARAMETERS_NAME":template, "BACK_SIZE":"256", "BACK_FILTERSIZE":"10", "BACKPHOTO_TYPE":"LOCAL", "BACKPHOTO_THICK":"20","DETECT_TYPE":"CCD", "DETECT_MINAREA":"5", "THRESH_TYPE":"RELATIVE","DETECT_THRESH":str(thresh),"ANALYSIS_THRESH":"10", "FILTER":"N", "FILTER_NAME":"default.conv", "CLEAN":"Y", "CLEAN_PARAM":"1.0", "DEBLEND_NTHRESH":"64", "DEBLEND_MINCONT":"0.0", "MASK_TYPE":"CORRECT", "VERBOSE_TYPE":"QUIET", "MEMORY_PIXSTACK":"1000000"})
os.system("sed 's/=//' "+tmpargs+" > "+parfile+"\n")
# apply sextractor
os.system("cd "+atlasdir+"/; for i in "+regexp+"; do sex -c "+parfile+" $i;mv test.cat "+atlasdir+"/$i.cat;done\n")
regtab = regexp.split("*")
# concatenate catalogs
cats = sort(glob.glob(atlasdir+"/*.cat"))
ns = 0
idx = -1
while ns==0:
idx = idx+1
try:
ncat = load(cats[idx])
ns = len(ncat)
except:
print "empty"
deg2rad = pi/180.0
assocdist = assocdist*deg2rad / 60
if isvector(ncat):
ncat = ncat.reshape(1,5)
fullcat = zeros((ns, 6))
fullcat[:,0] = glat2theta(ncat[:,1]) # glat ircs
fullcat[:,1] = glon2phi(ncat[:,0]) # glon ircs
fullcat[:,3] = int(cats[idx].split(regtab[0])[-1].split(regtab[1])[0]) # map index
fullcat[:,[2,4,5]] = ncat[:,[2,3,4]] # flux x and y coordinates
for i in range(idx+1, size(cats)):
try:
ncat = load(cats[i])
if isvector(ncat):
ncat = ncat.reshape(1,5)
ns = shape(ncat)[0]
nncat= zeros((ns, 6))
nncat[:,0] = glat2theta(ncat[:,1]) # glon ircs
nncat[:,1] = glon2phi(ncat[:,0]) # glat ircs
nncat[:,3] = int(cats[i].split(regtab[0])[-1].split(regtab[1])[0]) # map index
nncat[:,[2,4,5]] = ncat[:,[2,3,4]] # flux x and y coordinates
for s in range(ns):
vec1 = vstack((sin(fullcat[:,0])*cos(fullcat[:,1]), sin(fullcat[:,0])*sin(fullcat[:,1]), cos(fullcat[:,0])))
vec2 = vstack((sin(nncat[s,0])*cos(nncat[s,1]), sin(nncat[s,0])*sin(nncat[s,1]), cos(nncat[s,0])))
gamma = dot(vec1.transpose() , vec2.reshape(3,1))
for sj in range(len(gamma)):
gamma[sj] = math.acos(gamma[sj])
if len(find(gamma < assocdist))==0:
fullcat = vstack((fullcat,nncat[s,:]))
except:
print "empty"
## convert coordinates
tran = wcs.Transformation(wcs.icrs, wcs.galactic)
glon = phi2glon(fullcat[:,1])
glat = theta2glat(fullcat[:,0])
ndet = shape(fullcat)[0]
for j in range(ndet):
fullcat[j,0:2] = tran(array([glon[j], glat[j]]))
savetxt(catalog,fullcat)
os.system("rm "+atlasdir+"/atlas*.cat\n")
def make_histogram (x, y, xlog=False):
""" Build an histogram from two vectors.
Parameters
----------
x: array_like, shape (1, nx)
y: array_like, shape (1, n)
Returns
-------
hx: array-like, shape (1,2nx)
hy: array-like, shape (1,2nx)
jean-baptiste.melin@cea.fr
"""
if (xlog):
usex = log10(x)
else:
usex = x
nelts = size(x)
usex_mid = (usex[1:]+usex[0:-1])/2.
delta_usex_before = usex[1:]-usex_mid
delta_usex_after = usex_mid-usex[0:-1]
delta_usex_before = concatenate((array([delta_usex_after[0]]),delta_usex_before), axis=1)
delta_usex_after = concatenate((delta_usex_after,array([delta_usex_before[-1]])), axis=1)
delta_usex = delta_usex_before+delta_usex_after
dx = delta_usex
histo_factor = 10000
histo_usex = concatenate((array([ usex[0]-delta_usex_before[0]+delta_usex_before[0]/histo_factor]),array([usex[0]+delta_usex_after[0]-delta_usex_after[0]/histo_factor])))
for i in range(nelts-1):
histo_usex = concatenate((histo_usex, array([usex[i+1]-delta_usex_before[i+1]+delta_usex_before[i+1]/histo_factor]),array([usex[i+1]+delta_usex_after[i+1]-delta_usex_after[i+1]/histo_factor])))
if (xlog):
histo_x = exp(log(10)*histo_usex)
else:
histo_x = histo_usex
husex2usex = zeros((2*nelts), dtype=int)*-1
for i in range(nelts):
whi = intersect1d(find(histo_usex >= usex[i]-delta_usex_before[i]),
find(histo_usex < usex[i]+delta_usex_after[i]))
nwhi = size(whi)
if (nwhi >= 1):
husex2usex[whi] = i
if (xlog):
usey = log10(y)
else:
usey = y
dnydx = zeros(nelts)
for i in range(nelts):
whi = intersect1d(find(usey >= usex[i]-delta_usex_before[i]),
find(usey < usex[i]+delta_usex_after[i]))
nwhi = size(whi)
dnydx[i] = nwhi/delta_usex[i]
histo_dnydx = dnydx[husex2usex]
return [histo_x, histo_dnydx]
def associate_catalogs (out_table, in_table, das=5, cy_lim=2e-4, sdas=1, mdas=20):
""" Associate output clusters with input ones.
Parameters
----------
out_table: array_like, shape (nout, 2). Contains clusters coordinates (glon - glat) in degrees.
in_table : array_like, shape (nin, 4). Contains clusters description (glon - glat - flux - tv).
das: distance of association (arcmin)
cy_lim: integrated flux limit to use for the association
sdas: slope of the distance of association.
mdas: maximum value for the distance of association (arcmin)
Returns
-------
asso: array_like, shape (nout, 2). Contains associated cluster indice with
in_table in column 0 and its distance in column 1.
jean-baptiste.melin@cea.fr
"""
out_table = copy(out_table)
in_table = copy(in_table)
nin = shape(in_table)[0]
nout = shape(out_table)[0]
out = zeros((nout, 1)) ## boolean if 0 not associated if 1 already associated
asso = zeros((nin, 2), dtype=int)*-2
idx = find(in_table[:,2]>=cy_lim)
nidx = size(idx)
if (nidx>0):
asso[idx,:] = -1
## sort in descending fluxes
iso = flipud(argsort(in_table[:,2], axis=0))
in_table = in_table[iso, :]
## deg 2 radian
in_table [:,0:2] = in_table[:,0:2]*pi/180.0
out_table[:,0:2] = out_table[:,0:2]*pi/180.0
for i in range(nidx):
whkeep = find(out==0) ## searched within non associated
if (size(whkeep)==0):
break
else:
s = sin((out_table[whkeep, 1]-in_table[i, 1])/2.0) ##lat
l = sin((out_table[whkeep, 0]-in_table[i, 0])/2.0) ##lon
ci = cos(in_table [i, 1])
co = cos(out_table[whkeep, 1])
di = 2*arcsin(sqrt(s*s + co*ci*l*l))*180.0/pi
ma = max(das , sdas*in_table[i,3])/60.0
mi = min(mdas, sdas*in_table[i,3])/60.0
dilim = sdas*in_table[i,3]/60.0
#if (mdas>0):
# whd = intersect1d(find (di<=ma), find(di<=mi))
#else:
# whd = find( di<=ma)
if ((das/60.0)>= dilim):
dilim = das/60.0
if ((mdas/60.0)<= dilim):
dilim = mdas/60.0
whd = find(di<= dilim)
if (size(whd)>0):
whmax = argmin(di[whd])
asso[iso[i],0] = whkeep[whd[whmax]]
asso[iso[i],1] = di[whd[whmax]]
out[whkeep[whd[whmax]]] = 1
return asso
def plot_recovered (flux, idx, histx=None, figfile=None):
""" Plot simulated & recovered clusters.
Parameters
----------
flux: array-like, vector of input fluxes
idx: array-like, list of indices of recovered clusters
"""
if histx is None:
log10_cy = (arange(ceil((log10(max(flux))-log10(min(flux)))*10)+1)/10.0)+log10(min(flux))
capy = exp(log(10)*log10_cy)
else:
capy = histx
[histo_cy, histo_dndlogcy_input] = make_histogram (capy, flux, xlog=True)
[histo_cy, histo_dndlogcy_output] = make_histogram (capy, flux[idx], xlog=True)
figure()
loglog(histo_cy, histo_dndlogcy_input,color="r" )
loglog(histo_cy, histo_dndlogcy_output,color="b" )
axis([1e-4, 1, 3,5e4])
xlabel("$Y_{1r500}$ simulated [$arcmin^2$]", fontsize=18)
ylabel("$dN/dlogY_{1r500}$", fontsize=18)
title("Simulated & recovered clusters", fontsize=18)
legend(['Simulated (%4.0i clusters)' % size(flux), 'Recovered (%4.0i clusters)' % size(idx)])
if figfile is None:
show()
else:
savefig(figfile)
def plot_flux (fluxi, fluxo, figfile=None):
""" Plot simulated vs recovered flux.
Parameters
----------
fluxi: array-like, vector of input fluxes
fluxo: array-like, vector of output fluxes
"""
figure()
loglog(fluxi, fluxo,color="b" ,marker="p", linestyle=".")
loglog(arange(1e-5,1,1e-2), arange(1e-5,1,1e-2), color="k")
axis([1e-5, 1, 1e-5,1])
xlabel("$Y_{5r500}$ simulated [$arcmin^2$]", fontsize=18)
ylabel("$Y_{5r500}$ recovered [$arcmin^2$]", fontsize=18)
score = sqrt(sum((fluxo-fluxi)*(fluxo-fluxi)))
title("SZ flux : "+str(score), fontsize=18)
if figfile is None:
show()
else:
savefig(figfile)
def plot_radius (radiusi, radiuso, figfile=None):
""" Plot simulated vs recovered flux.
Parameters
----------
radiusi: array-like, vector of input radius
radiuso: array-like, vector of output radius
"""
figure()
loglog(radiusi, radiuso,color="b" ,marker="p", linestyle=".")
loglog(arange(1,1e3,1e2), arange(1,1e3,1e2), color="k")
axis([1, 1e3, 1,1e3])
xlabel("$\theta_{500}$ simulated [$arcmin$]", fontsize=18)
ylabel("$\theta_{500}$ recovered [$arcmin$]", fontsize=18)
score = sqrt(sum((radiuso-radiusi)*(radiuso-radiusi)))
title("SZ radius : "+ str(score), fontsize=18)
if figfile is None:
show()
else:
savefig(figfile)
def plot_completeness (flux, idx, fluxcut,histx=None, figfile=None ):
""" Plot completeness.
Parameters
----------
flux: array-like, vector of input fluxes
idx: array-like, list of indices of recovered clusters
fluxcut: limit on flux.
"""
if histx is None:
log10_cy = (arange(ceil((log10(max(flux))-log10(min(flux)))*10)+1)/10.0)+log10(min(flux))
capy = exp(log(10)*log10_cy)
else:
capy = histx
[histo_cy, histo_dndlogcy_input] = make_histogram (capy, flux, xlog=True)
[histo_cy, histo_dndlogcy_output] = make_histogram (capy, flux[idx], xlog=True)
whp = find(histo_dndlogcy_input>0)
completeness = zeros(shape(histo_dndlogcy_input))
completeness[whp] = histo_dndlogcy_output[whp]/histo_dndlogcy_input[whp]
whn1 = intersect1d(find(histo_dndlogcy_input == 0), find(histo_cy > fluxcut))
completeness[whn1] = 1
whn0 = intersect1d(find(histo_dndlogcy_input == 0), find(histo_cy <= fluxcut))
completeness[whn0] = 0
figure()
semilogx( histo_cy, completeness, color="k")
axis([1e-5,1,0,1])
xlabel("$Y_{1r500}$ simulated [$arcmin^2$]", fontsize=18)
ylabel("Recovered/Simulated", fontsize=18)
title("Completeness", fontsize=18)
if figfile is None:
show()
else:
savefig(figfile)
def plot_false (flux, idx,histx=None, figfile=None):
""" Plot detections and false detections
Parameters
----------
flux: array-like, vector of output fluxes
idx: array-like, list of indices of false detections
"""
if histx is None:
log10_cy = (arange(ceil((log10(max(flux))-log10(min(flux)))*10)+1)/10.0)+log10(min(flux))
capy = exp(log(10)*log10_cy)
else:
capy = histx
[histo_cy, histo_dndlogcy_output_all] = make_histogram (capy, flux, xlog=True)
[histo_cy, histo_dndlogcy_output_false] = make_histogram (capy, flux[idx], xlog=True)
figure()
loglog(histo_cy, histo_dndlogcy_output_all)
loglog(histo_cy, histo_dndlogcy_output_false)
axis([1e-5, 1, 1, 1e6])
xlabel("$Y_{1r500}$ recovered [$arcmin^2$]", fontsize=18)
ylabel("$dN/dlogY_{1r500}$", fontsize=18)
title("Detections & false detections", fontsize=18)
legend(['Detections (%4.0i clusters)' % shape(flux)[0], 'False (%4.0i clusters)' % size(idx)])
if figfile is None:
show()
else:
savefig(figfile)
def plot_contamination(flux, idx, histx=None,figfile=None):
""" Plot contamination
Parameters
----------
flux: array-like, vector of output fluxes
idx: array-like, list of indices of false detections
"""
if histx is None:
log10_cy = (arange(ceil((log10(max(flux))-log10(min(flux)))*10)+1)/10.0)+log10(min(flux))
capy = exp(log(10)*log10_cy)
else:
capy = histx
[histo_cy, histo_dndlogcy_output_all] = make_histogram (capy, flux, xlog=True)
[histo_cy, histo_dndlogcy_output_false] = make_histogram (capy, flux[idx], xlog=True)
whp = find(histo_dndlogcy_output_all >0)
contamination = zeros(shape(histo_dndlogcy_output_all))
contamination[whp] = histo_dndlogcy_output_false[whp]/histo_dndlogcy_output_all[whp]
figure()
semilogx(histo_cy, contamination)
axis([1e-5, 1, 0, 1])
xlabel("$Y_{1r500}$ recovered [$arcmin^2$]", fontsize=18)
ylabel("False/Detections", fontsize=18)
title("Contamination", fontsize=18)
legend(['Total %2.1f %%' % (size(idx)/float(size(flux))*100)])
if figfile is None:
show()
else:
savefig(figfile)
def sn_max (sn):
""" Return S/N ratio limit
Parameters
----------
sn: array-like
"""
min_counts = size(sn) / 2 ## TO DO : this is a parameter
n_steps = 1000
max_sn = (max(sn))
min_sn = (min(sn))
sn_step = (max_sn-min_sn)/n_steps
sn_tab = (arange(ceil((max_sn-min_sn)/sn_step))*sn_step)+min_sn
output_counts = zeros(size(sn_tab))
for i in range(n_steps):
idx = find(sn>= sn_tab[i])
output_counts[i] = size(idx)
use_max_sn = interp(min_counts, flipud(output_counts), flipud(sn_tab))
return use_max_sn
def plot_completeness_purity (purity ,completeness, figfile=None):
""" Plot completeness versus purity.
Parameters
----------
purity: array-like
completness: array-like
"""
figure()
plot(purity, completeness)
axis([0,1,0,0.15])
xlabel("Purity=1-Contamination", fontsize=18)
ylabel("Completeness", fontsize=18)
title ("Completeness vs. Purity", fontsize=18)
if figfile is None:
show()
else:
savefig(figfile)
def eval_catalogs (input_cat, output_cat, latcut=20, fluxcut=2e-4, purity_thres=0.9, dir_plot=None):
""" Evaluate catalog.
Parameters
----------
output_cat: file which contains output clusters description (pixel lon-deg lat-deg cy cyerr tv tverr sn)
input_cat : file which contains input clusters description (glon glat flux tv).
fluxcut: integrated flux limit to use for the association
latcut: latitude limit.
purity_thresh: purity limit.
"""
in_table = loadtxt (input_cat) #lon-deg lat-deg cy tv
out_table = loadtxt( output_cat) #pixel lon-deg lat-deg cy cyerr tv tverr sn
## select cluster above latitude cut
idx = find(abs(in_table[:,1])>=latcut)
in_table = in_table[idx, :]
idx = find(abs(out_table[:,2])>=latcut)
out_table = out_table[idx, :]
## select cluster above flux cut
idx = find(in_table[:,2]>=fluxcut)
in_table = in_table[idx]
## compute completeness purity
sn = out_table[:,7]
min_sn = min(sn)
max_sn = sn_max (sn)
n_steps = 20
sn_step = (max_sn-min_sn)/(n_steps)
sn_tab = (arange(ceil((max_sn-min_sn)/sn_step))*sn_step)+min_sn
comparison = zeros((n_steps, 3)) # sn purity completeness
comparison[:,0] = sn_tab
for i in range(n_steps):
ins = find(sn>sn_tab[i])
sub_table = out_table[ins,1:3]
asso = associate_catalogs (sub_table, in_table, cy_lim=fluxcut)
ia = find(asso[:,0]>=0)
comparison[i,1] = size(ia)/float(size(ins))
comparison[i,2] = size(ia)/float(shape(in_table)[0])
isort = argsort(comparison[:,1])
sn_thres = interp(purity_thres, comparison[isort,1], comparison[isort,0])
## select cluster above sn
idx = find(sn>sn_thres)
out_table = out_table[idx,:]
## do the association
asso = associate_catalogs(out_table[:,1:3], in_table, cy_lim=fluxcut)
wha = find(asso[:,0]>=0)
detected = zeros((shape(out_table)[0]), dtype='int')
detected[asso[wha,0]]=1
whna = find(detected == 0)
false_cat = out_table[whna,:]
## build histogram range
min_cy = 1e-5
max_cy = 1.0
step_log10_capy = 10.0
log10_cy = (arange(ceil((log10(max_cy)-log10(min_cy))*step_log10_capy)+1)/step_log10_capy)+log10(min_cy)
capy = exp(log(10)*log10_cy)
## plots
if dir_plot is None:
dir_plot = pwd
plot_recovered (in_table[:,2], wha, histx=capy, figfile=dir_plot+"/recovered.png")
plot_completeness (in_table[:,2], wha, fluxcut,histx=capy, figfile=dir_plot+"/completeness.png")
plot_false (out_table[:,3], whna, histx=capy,figfile=dir_plot+"/false.png" )
plot_contamination(out_table[:,3], whna, histx=capy,figfile=dir_plot+"/contamination.png" )
plot_completeness_purity (comparison[:,1],comparison[:,2],figfile=dir_plot+"/comp_purity.png")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment