We will analyse the changes in the classification using the new sigma and the new catalogue without the galaxies that went to LGZ
In [1]:
import numpy as np
from astropy.table import Table, join
from astropy import units as u
from astropy.coordinates import SkyCoord, search_around_sky
from IPython.display import clear_output
import pickle
import os
In [3]:
from mltier1 import (get_center, Field, parallel_process, describe)
In [2]:
%load_ext autoreload
In [4]:
%autoreload
In [5]:
from IPython.display import clear_output
In [6]:
import matplotlib as mpl
# mpl.rc('figure', figsize=(6.64, 6.64*0.74), dpi=100)
# mpl.rc('figure.subplot', left=0.15, right=0.95, bottom=0.15, top=0.92)
# mpl.rc('lines', linewidth=1.75, markersize=8.0, markeredgewidth=0.75)
# mpl.rc('font', size=18.0, family="serif", serif="CM")
# mpl.rc('xtick', labelsize='small')
# mpl.rc('ytick', labelsize='small')
# mpl.rc('xtick.major', width=1.0, size=8)
# mpl.rc('ytick.major', width=1.0, size=8)
# mpl.rc('xtick.minor', width=1.0, size=4)
# mpl.rc('ytick.minor', width=1.0, size=4)
# mpl.rc('axes', linewidth=1.5)
# mpl.rc('legend', fontsize='small', numpoints=1, labelspacing=0.4, frameon=False)
# mpl.rc('text', usetex=True)
# mpl.rc('savefig', dpi=300)
In [7]:
%pylab inline
In [8]:
def most_common(a, n=2):
u, c = np.unique(a, return_counts=True)
order = np.argsort(c)
for i in range(n):
print(c[order][-(i+1)], u[order][-(i+1)])
In [9]:
save_intermediate = True
plot_intermediate = True
In [10]:
idp = "idata/final_analysis_pdf_v1.0"
In [11]:
if not os.path.isdir(idp):
os.makedirs(idp)
In [12]:
pwli = Table.read("lofar_pw_pdf.fits")
In [13]:
pwli.colnames
Out[13]:
In [14]:
lofar_all = Table.read("data/LOFAR_HBA_T1_DR1_merge_ID_optical_v1.1b.fits")
In [15]:
pwl = join(pwli, lofar_all[['Source_Name', 'AllWISE', 'objID', 'ML_LR',
'ID_flag', 'ID_name', 'ID_ra', 'ID_dec',
'LGZ_Size', 'LGZ_Width', 'LGZ_PA', 'LGZ_Assoc',
'LGZ_Assoc_Qual', 'LGZ_ID_Qual']],
join_type='left',
keys='Source_Name',
uniq_col_name='{col_name}{table_name}',
table_names=['', '_input'])
In [16]:
colour_limits = [0.0, 0.5, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.5, 4.0]
In [17]:
bin_list, centers, Q_0_colour, n_m, q_m = pickle.load(open("lofar_params.pckl", "rb"))
In [18]:
for col in pwl.colnames:
fv = pwl[col].fill_value
typ = pwl[col].dtype
print(col, fv, typ)
# Restore NaNs
if fv == 1e+20:
pwl[col][(pwl[col] == fv)] = np.nan
# if (isinstance(fv, np.float64) and (fv != 1e+20)):
# print(col, fv)
# pwl[col].fill_value = 1e+20
In [19]:
pwl["colour"][(pwl["colour"] == 1e+20)] = np.nan
In [20]:
describe(pwl["colour"])
Change the AllWISE_input that are '' to 'N/A'. That comes from a previous error with the fill value.
In [21]:
pwl["AllWISE_input"][pwl["AllWISE_input"] == ""] = "N/A"
In [22]:
plt.rcParams["figure.figsize"] = (12,10)
from matplotlib import cm
from matplotlib.collections import LineCollection
cm_subsection = linspace(0., 1., 16)
colors = [ cm.viridis(x) for x in cm_subsection ]
low = np.nonzero(centers[1] >= 15)[0][0]
high = np.nonzero(centers[1] >= 22.2)[0][0]
fig, a = plt.subplots()
for i, q_m_k in enumerate(q_m):
#plot(centers[i], q_m_old[i]/n_m_old[i])
a = subplot(4,4,i+1)
if i not in [-1]:
q_m_aux = q_m[i]/np.sum(q_m[i])
lwidths = (q_m_aux/np.max(q_m_aux)*10).astype(float)
#print(lwidths)
y_aux = q_m_k/n_m[i]
factor = np.max(y_aux[low:high])
y = y_aux
#print(y)
x = centers[i]
points = np.array([x, y]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
lc = LineCollection(segments, linewidths=lwidths, color=colors[i])
a.add_collection(lc)
#plot(centers[i], x/factor, color=colors[i-1])
xlim([12, 30])
if i == 0:
xlim([10, 23])
ylim([0, 1.2*factor])
subplots_adjust(left=0.125,
bottom=0.1,
right=0.9,
top=0.9,
wspace=0.4,
hspace=0.2)
Explanation of the KDE used in the computing of the q(m) and n(m):
In [29]:
save = True
save_pdf = True
plt.rcParams["figure.figsize"] = (6.64, 6.64*0.74)
plt.rcParams["figure.dpi"] = 100
plt.rcParams['lines.linewidth'] = 1.75
plt.rcParams['lines.markersize'] = 8.0
plt.rcParams['lines.markeredgewidth'] = 0.75
plt.rcParams['font.size'] = 15.0 ## not 18.0
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.serif'] = "CM"
plt.rcParams['xtick.labelsize'] = 'small'
plt.rcParams['ytick.labelsize'] = 'small'
plt.rcParams['xtick.major.width'] = 1.0
plt.rcParams['xtick.major.size'] = 8
plt.rcParams['ytick.major.width'] = 1.0
plt.rcParams['ytick.major.size'] = 8
plt.rcParams['xtick.minor.width'] = 1.0
plt.rcParams['xtick.minor.size'] = 4
plt.rcParams['ytick.minor.width'] = 1.0
plt.rcParams['ytick.minor.size'] = 4
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['legend.fontsize'] = 'small'
plt.rcParams['legend.numpoints'] = 1
plt.rcParams['legend.labelspacing'] = 0.4
plt.rcParams['legend.frameon'] = False
plt.rcParams['text.usetex'] = True
plt.rcParams['savefig.dpi'] = 300
from matplotlib import cm
from matplotlib.collections import LineCollection
from matplotlib.lines import Line2D
cm_subsection = linspace(0., 1., 16)
colors = [ cm.viridis(x) for x in cm_subsection ]
low = np.nonzero(centers[1] >= 15)[0][0]
high = np.nonzero(centers[1] >= 22.2)[0][0]
cm_subsection = linspace(0., 1., 14)
colors = [ cm.jet(x) for x in cm_subsection ]
fig, a = plt.subplots()
lcs = []
proxies = []
def make_proxy(zvalue, scalar_mappable, **kwargs):
color = scalar_mappable.cmap(scalar_mappable.norm(zvalue))
return Line2D([0, 1], [0, 1], color=color, **kwargs)
for i, q_m_k in enumerate(q_m):
#plot(centers[i], q_m_old[i]/n_m_old[i])
q_m_aux = q_m[i]/np.sum(q_m[i])
lwidths = (q_m_aux/np.max(q_m_aux)*10).astype(float)
if save_pdf and (i<6): # Solve problems with the line with in pdfs
lwidths[lwidths < 0.005] = 0
#print(lwidths)
y_aux = q_m_k/n_m[i]
factor = np.max(y_aux[low:high])
y = y_aux
#print(y)
x = centers[i]
points = np.array([x, y]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
#lc = LineCollection(segments, linewidths=lwidths, color=colors[i])
if i not in [0]:
if i == 1:
color = "k"
else:
color = colors[i-2]
lcs.append(LineCollection(segments, linewidths=lwidths, color=color))
proxies.append(Line2D([0, 1], [0, 1], color=color, lw=5))
a.add_collection(lcs[-1])
xlim([16, 26])
ylim([0, 700000])
xlabel("$i$ magnitude")
ylabel("$q(m,c)/n(m,c)$")
inset = plt.axes([0.285, 0.6, .2, .25])
q_m_aux = q_m[0]/np.sum(q_m[0])
lwidths = (q_m_aux/np.max(q_m_aux)*10).astype(float)
y_aux = q_m[0]/n_m[0]
factor = np.max(y_aux[low:high])
y = y_aux
x = centers[0]
points = np.array([x, y]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
lc = LineCollection(segments, linewidths=lwidths, color="k")
proxy = Line2D([0, 1], [0, 1], color="k", lw=5)
inset.add_collection(lc)
xlim([15, 23])
ylim([0, 10000])
xlabel("$W1$ magnitude")
ylabel("$q(m)/n(m)$")
inset.legend([proxy], ["only $W1$"], fontsize="xx-small")
#ylim([0, 1.2*factor])
#print(lcs)
# a.legend(proxies, ["$i-W1 < 0.0$",
# "$0.0 \leq i-W1 < 0.5$",
# "$0.5 \leq i-W1 < 1.0$",
# "$1.0 \leq i-W1 < 1.25$",
# "$1.25 \leq i-W1 < 1.5$",
# "$1.5 \leq i-W1 < 1.75$",
# "$1.75 \leq i-W1 < 2.0$",
# "$2.0 \leq i-W1 < 2.25$",
# "$2.25 \leq i-W1 < 2.5$",
# "$2.5 \leq i-W1 < 2.75$",
# "$2.75 \leq i-W1 < 3.0$",
# "$3.0 \leq i-W1 < 3.5$",
# "$3.5 \leq i-W1 < 4.0$",
# "$i-W1 \geq 4.0$"])
# a.legend(proxies, ["only $i$",
# "$i-W1 \in (-\infty, 0.0)$",
# "$i-W1 \in [0.0, 0.5)$",
# "$i-W1 \in [0.5, 1.0)$",
# "$i-W1 \in [1.0, 1.25)$",
# "$i-W1 \in [1.25, 1.5)$",
# "$i-W1 \in [1.5, 1.75)$",
# "$i-W1 \in [1.75, 2.0)$",
# "$i-W1 \in [2.0, 2.25)$",
# "$i-W1 \in [2.25, 2.5)$",
# "$i-W1 \in [2.5, 2.75)$",
# "$i-W1 \in [2.75, 3.0)$",
# "$i-W1 \in [3.0, 3.5)$",
# "$i-W1 \in [3.5, 4.0)$",
# "$i-W1 \in [4.0, \infty)$"],
# fontsize="xx-small")
a.legend(proxies, ["only $i$",
"$(-\infty, 0.0)$",
"$[0.0, 0.5)$",
"$[0.5, 1.0)$",
"$[1.0, 1.25)$",
"$[1.25, 1.5)$",
"$[1.5, 1.75)$",
"$[1.75, 2.0)$",
"$[2.0, 2.25)$",
"$[2.25, 2.5)$",
"$[2.5, 2.75)$",
"$[2.75, 3.0)$",
"$[3.0, 3.5)$",
"$[3.5, 4.0)$",
"$[4.0, \infty)$"],
fontsize="xx-small",
title="$i-W1$")
# subplots_adjust(left=0.15,
# bottom=0.1,
# right=0.9,
# top=0.9,
# wspace=0.4,
# hspace=0.2)
subplots_adjust(left=0.15,
bottom=0.15,
right=0.95,
top=0.92,
wspace=0.4,
hspace=0.2)
if save:
plt.savefig("idata/q_n_m.png")
plt.savefig("idata/q_n_m.svg")
plt.savefig("idata/q_n_m_high.png", dpi=800)
if save_pdf:
plt.savefig("idata/q_n_m.pdf")
In [30]:
q0 = np.sum(Q_0_colour)
In [31]:
def completeness(lr, threshold, q0):
n = len(lr)
lrt = lr[lr < threshold]
return 1. - np.sum((q0 * lrt)/(q0 * lrt + (1 - q0)))/float(n)/q0
def reliability(lr, threshold, q0):
n = len(lr)
lrt = lr[lr > threshold]
return 1. - np.sum((1. - q0)/(q0 * lrt + (1 - q0)))/float(n)/q0
completeness_v = np.vectorize(completeness, excluded=[0])
reliability_v = np.vectorize(reliability, excluded=[0])
In [32]:
pwl["lrt"] = pwl["lr"]
pwl["lrt"][np.isnan(pwl["lr"])] = 0
In [33]:
n_test = 100
threshold_mean = np.percentile(pwl["lrt"], 100*(1 - q0))
In [34]:
thresholds = np.arange(0., 10., 0.01)
thresholds_fine = np.arange(0.1, 2., 0.001)
completeness_t = completeness_v(pwl["lrt"], thresholds, q0)
reliability_t = reliability_v(pwl["lrt"], thresholds, q0)
average_t = (completeness_t + reliability_t)/2
In [35]:
completeness_t_fine = completeness_v(pwl["lrt"], thresholds_fine, q0)
reliability_t_fine = reliability_v(pwl["lrt"], thresholds_fine, q0)
average_t_fine = (completeness_t_fine + reliability_t_fine)/2
In [36]:
thresholds_fine[np.argmax(average_t_fine)]
Out[36]:
In [37]:
thresholds_fine[np.argmin(np.abs(completeness_t_fine-reliability_t_fine))]
Out[37]:
In [38]:
np.sum(pwl["lrt"] >= 0.358)
Out[38]:
In [39]:
np.sum(pwl["lrt"] >= 0.639)
Out[39]:
In [40]:
np.sum(pwl["lrt"] >= 0.639)
Out[40]:
In [41]:
len(pwl)
Out[41]:
In [42]:
len(pwl[pwl['ID_flag'] == 1])
Out[42]:
In [43]:
n0 = np.sum((pwl["lrt"] >= 0.639) & (pwl['ID_flag'] == 1))
print(n0, n0/len(pwl[pwl['ID_flag'] == 1]))
In [44]:
n0 = np.sum((pwl["lrt"] >= 0.81268) & (pwl['ID_flag'] == 1))
print(n0, n0/len(pwl[pwl['ID_flag'] == 1]))
In [45]:
n0 = np.sum((pwl["ML_LR"] >= 0.81268) & (pwl['ID_flag'] == 1))
print(n0, n0/len(pwl[pwl['ID_flag'] == 1]))
In [46]:
np.sum(pwl["lrt"] >= 0.639)/len(pwl)
Out[46]:
In [47]:
threshold_sel = 0.639
In [48]:
thresholds_fine = np.arange(0.0, 2., 0.001)
completeness_t_fine = completeness_v(pwl["lrt"], thresholds_fine, q0)
reliability_t_fine = reliability_v(pwl["lrt"], thresholds_fine, q0)
In [49]:
plt.rcParams["figure.figsize"] = (6.64, 6.64*0.74)
plt.rcParams["figure.dpi"] = 100
plt.rcParams['lines.linewidth'] = 1.75
plt.rcParams['lines.markersize'] = 8.0
plt.rcParams['lines.markeredgewidth'] = 0.75
plt.rcParams['font.size'] = 18.0
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.serif'] = "CM"
plt.rcParams['xtick.labelsize'] = 'small'
plt.rcParams['ytick.labelsize'] = 'small'
plt.rcParams['xtick.major.width'] = 1.0
plt.rcParams['xtick.major.size'] = 8
plt.rcParams['ytick.major.width'] = 1.0
plt.rcParams['ytick.major.size'] = 8
plt.rcParams['xtick.minor.width'] = 1.0
plt.rcParams['xtick.minor.size'] = 4
plt.rcParams['ytick.minor.width'] = 1.0
plt.rcParams['ytick.minor.size'] = 4
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['legend.fontsize'] = 'small'
plt.rcParams['legend.numpoints'] = 1
plt.rcParams['legend.labelspacing'] = 0.4
plt.rcParams['legend.frameon'] = False
plt.rcParams['text.usetex'] = True
plt.rcParams['savefig.dpi'] = 300
plot(thresholds_fine, completeness_t_fine, "-", label="Completeness")
plot(thresholds_fine, reliability_t_fine, "-", label="Reliability")
text(0.66, 0.971, "0.639")
#plot(thresholds_fine, average_t_fine, "-", label="average")
vlines(threshold_sel, 0.9, 1., "k", linestyles="dashed", label="Threshold\nselected")
#vlines(threshold_mean, 0.9, 1., "y", linestyles="dashed")
ylim([0.97, 1.])
xlim([0.0, 1.5])
legend(loc=4)
xlabel("Threshold")
ylabel("Completeness/Reliability")
# subplots_adjust(left=0.2,
# bottom=0.1,
# right=0.95,
# top=0.9,
# wspace=0.4,
# hspace=0.2)
subplots_adjust(left=0.15,
bottom=0.15,
right=0.95,
top=0.92,
wspace=0.4,
hspace=0.2)
save = True
if save:
plt.savefig("idata/completeness_reliability.png")
plt.savefig("idata/completeness_reliability.svg")
plt.savefig("idata/completeness_reliability.pdf")
plt.savefig("idata/completeness_reliability_high.png", dpi=800)
In [50]:
cond_mlr = (pwl['ID_flag'] == 1)
In [51]:
len(pwl)
Out[51]:
In [52]:
np.sum(np.isnan(pwl["category"]))
Out[52]:
In [53]:
np.sum(np.isnan(pwl["category"][cond_mlr]))
Out[53]:
In [54]:
n_c, n_c_mlr = [], []
for i in np.unique(pwl["category"][~np.isnan(pwl["category"])]):
n_c.append(np.sum((pwl["category"] == i)))
n_c_mlr.append(np.sum((pwl["category"][cond_mlr] == i)))
print(i, n_c[-1], n_c_mlr[-1])
In [55]:
total_n_c = np.sum(n_c)
total_n_c_mlr = np.sum(n_c_mlr)
print(total_n_c, total_n_c_mlr)
print(len(pwl)-np.sum(np.isnan(pwl["category"])), len(pwl)-np.sum(np.isnan(pwl["category"][cond_mlr])))
for i in range(16):
print("{:2d} {:6d} {:6d} {:6.3f} {:6.3f} {:.2%}".format(i,
n_c[i],
n_c_mlr[i],
n_c[i]/total_n_c*100,
n_c_mlr[i]/total_n_c_mlr*100,
(n_c[i]-n_c_mlr[i])/n_c[i]
))
In [56]:
total_n_c
Out[56]:
Load the full catalogue of galaxies and work out the fraction of matched LOFAR. It will be necessary to work in the restricted area with both the pw and the lofar matched catalogue. The fractions could be corrected if the relative fractions of sources of each type change between the full area and the restricted area for the LOFAR matches.
In [57]:
field = Field(170.0, 190.0, 46.8, 55.9)
In [58]:
combined = Table.read("pw.fits")
In [59]:
combined = field.filter_catalogue(combined,
colnames=("ra", "dec"))
In [60]:
combined["colour"] = combined["i"] - combined["W1mag"]
In [61]:
colour_limits = [0.0, 0.5, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.5, 4.0]
In [62]:
combined_panstarrs = (~np.isnan(combined["i"]) & np.isnan(combined["W1mag"])) # Sources with only i-band
combined_wise =(np.isnan(combined["i"]) & ~np.isnan(combined["W1mag"])) # Sources with only W1-band
In [63]:
# Start with the W1-only, i-only and "less than lower colour" bins
colour_bin_def = [{"name":"only W1", "condition": combined_wise},
{"name":"only i", "condition": combined_panstarrs},
{"name":"-inf to {}".format(colour_limits[0]),
"condition": (combined["colour"] < colour_limits[0])}]
# Get the colour bins
for i in range(len(colour_limits)-1):
name = "{} to {}".format(colour_limits[i], colour_limits[i+1])
condition = ((combined["colour"] >= colour_limits[i]) &
(combined["colour"] < colour_limits[i+1]))
colour_bin_def.append({"name":name, "condition":condition})
# Add the "more than higher colour" bin
colour_bin_def.append({"name":"{} to inf".format(colour_limits[-1]),
"condition": (combined["colour"] >= colour_limits[-1])})
In [64]:
combined["category"] = np.nan
for i in range(len(colour_bin_def)):
combined["category"][colour_bin_def[i]["condition"]] = i
In [65]:
numbers_combined_bins = np.array([np.sum(a["condition"]) for a in colour_bin_def])
In [66]:
numbers_combined_bins
Out[66]:
In [67]:
pwlf = field.filter_catalogue(pwl[cond_mlr],
colnames=("ra", "dec"))
In [68]:
pwlf_all = field.filter_catalogue(pwl,
colnames=("ra", "dec"))
In [69]:
lofar_c, pw_c, lofar_c_all = [], [], []
avg_colour = []
for i in range(16):
lofar_c.append(np.sum(pwlf["category"] == i))
lofar_c_all.append(np.sum(pwlf_all["category"] == i))
pw_c.append(np.sum(combined["category"] == i))
avg_colour.append(np.nanmedian(combined["colour"][combined["category"] == i]))
print("{:2d} {:6.2f} {:6d} {:7d} {:6.3f} +/- {:6.3f}".format(i, avg_colour[-1],
lofar_c_all[-1], pw_c[-1],
lofar_c_all[-1]/pw_c[-1], np.sqrt(lofar_c_all[-1])/pw_c[-1]))
#print("{:2d} {:6d} {:7d} {:5.3f}".format(i, lofar_c[-1], pw_c[-1], lofar_c[-1]/pw_c[-1]))
In [70]:
np.sum(np.array(lofar_c_all)/np.array(pw_c))
Out[70]:
In [71]:
np.sum(np.array(lofar_c_all))/np.sum(np.array(pw_c))
Out[71]:
In [72]:
len(pwl)/len(combined)
Out[72]:
In [73]:
lc = np.array(lofar_c_all)
pc = np.array(pw_c)
plt.rcParams["figure.figsize"] = (6.64, 6.64*0.74)
plt.rcParams["figure.dpi"] = 100
plt.rcParams['lines.linewidth'] = 1.75
plt.rcParams['lines.markersize'] = 8.0
plt.rcParams['lines.markeredgewidth'] = 0.75
plt.rcParams['font.size'] = 18.0
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.serif'] = "CM"
plt.rcParams['xtick.labelsize'] = 'small'
plt.rcParams['ytick.labelsize'] = 'small'
plt.rcParams['xtick.major.width'] = 1.0
plt.rcParams['xtick.major.size'] = 8
plt.rcParams['ytick.major.width'] = 1.0
plt.rcParams['ytick.major.size'] = 8
plt.rcParams['xtick.minor.width'] = 1.0
plt.rcParams['xtick.minor.size'] = 4
plt.rcParams['ytick.minor.width'] = 1.0
plt.rcParams['ytick.minor.size'] = 4
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['legend.fontsize'] = 'small'
plt.rcParams['legend.numpoints'] = 1
plt.rcParams['legend.labelspacing'] = 0.4
plt.rcParams['legend.frameon'] = False
plt.rcParams['text.usetex'] = True
plt.rcParams['savefig.dpi'] = 300
cm_subsection = linspace(0., 1., 14)
colors = [ cm.jet(x) for x in cm_subsection ]
# i-only
plot([0,0], [0,lc[1]/pc[1]*100], marker=",", ls="-", color="k")
scatter([0], lc[1]/pc[1]*100, s=lc[1]/10, c="k")
plot([0], lc[1]/pc[1]*100, marker=".", ls="", color="w", ms=2)
# w1-only
plot([15,15], [0,lc[0]/pc[0]*100], marker=",", ls="-", color="k")
scatter([15], lc[0]/pc[0]*100, s=lc[0]/10, c="k")
plot([15], lc[0]/pc[0]*100, marker=".", ls="", color="w", ms=2)
# colours
for i in range(14):
plot([i+1,i+1], [0,lc[i+2]/pc[i+2]*100], marker=",", ls="-", color=colors[i])
scatter([i+1], lc[i+2]/pc[i+2]*100, s=lc[i+2]/10, c=colors[i])
if i == 0:
plot([i+1], lc[i+2]/pc[i+2]*100, marker=".", ls="", color="w", ms=1)
else:
plot([i+1], lc[i+2]/pc[i+2]*100, marker=".", ls="", color="w", ms=2)
#text(0.66, 0.971, "0.639")
ylim([0., 14.])
xlim([-0.5, 15.5])
#legend(fontsize="small")
xlabel("Colour category")
ylabel("Fraction of galaxies detected by LoTSS\n(per cent)")
xt = ["$i$ only",
"$(-\infty, 0.0)$",
"$[0.0, 0.5)$",
"$[0.5, 1.0)$",
"$[1.0, 1.25)$",
"$[1.25, 1.5)$",
"$[1.5, 1.75)$",
"$[1.75, 2.0)$",
"$[2.0, 2.25)$",
"$[2.25, 2.5)$",
"$[2.5, 2.75)$",
"$[2.75, 3.0)$",
"$[3.0, 3.5)$",
"$[3.5, 4.0)$",
"$[4.0, \infty)$",
"$W1$ only"]
xticks(np.arange(16), xt, rotation=90)
subplots_adjust(left=0.2,
bottom=0.3,
right=0.95,
top=0.9,
wspace=0.4,
hspace=0.2)
save = True
if save:
plt.savefig("idata/fractiono.png")
plt.savefig("idata/fractiono.svg")
plt.savefig("idata/fractiono.pdf")
plt.savefig("idata/fractiono_high.png", dpi=800)
In [74]:
plt.rcParams["figure.figsize"] = (6.64, 6.64*0.74)
plt.rcParams["figure.dpi"] = 100
plt.rcParams['lines.linewidth'] = 1.75
plt.rcParams['lines.markersize'] = 8.0
plt.rcParams['lines.markeredgewidth'] = 0.75
plt.rcParams['font.size'] = 18.0
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.serif'] = "CM"
plt.rcParams['xtick.labelsize'] = 'small'
plt.rcParams['ytick.labelsize'] = 'small'
plt.rcParams['xtick.major.width'] = 1.0
plt.rcParams['xtick.major.size'] = 8
plt.rcParams['ytick.major.width'] = 1.0
plt.rcParams['ytick.major.size'] = 8
plt.rcParams['xtick.minor.width'] = 1.0
plt.rcParams['xtick.minor.size'] = 4
plt.rcParams['ytick.minor.width'] = 1.0
plt.rcParams['ytick.minor.size'] = 4
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['legend.fontsize'] = 'small'
plt.rcParams['legend.numpoints'] = 1
plt.rcParams['legend.labelspacing'] = 0.4
plt.rcParams['legend.frameon'] = False
plt.rcParams['text.usetex'] = True
plt.rcParams['savefig.dpi'] = 300
cm_subsection = linspace(0., 1., 14)
colors = [ cm.jet(x) for x in cm_subsection ]
# # i-only
# plot([0,0], [0,lc[1]/pc[1]*100], marker=",", ls="-", color="k")
# scatter([0], lc[1]/pc[1]*100, s=lc[1]/10, c="k")
# plot([0], lc[1]/pc[1]*100, marker=".", ls="", color="w", ms=2)
# # w1-only
# plot([15,15], [0,lc[0]/pc[0]*100], marker=",", ls="-", color="k")
# scatter([15], lc[0]/pc[0]*100, s=lc[0]/10, c="k")
# plot([15], lc[0]/pc[0]*100, marker=".", ls="", color="w", ms=2)
# colours
for i in range(14):
#plot([avg_colour[i+2],avg_colour[i+2]], [0,lc[i+2]/pc[i+2]*100], marker=",", ls="-", color=colors[i])
scatter(avg_colour[i+2], lc[i+2]/pc[i+2]*100, s=lc[i+2]/10, c=colors[i])
if i == 0:
plot(avg_colour[i+2], lc[i+2]/pc[i+2]*100, marker=".", ls="", color="w", ms=1)
else:
plot(avg_colour[i+2], lc[i+2]/pc[i+2]*100, marker=".", ls="", color="w", ms=2)
#text(0.66, 0.971, "0.639")
ylim([0., 14.])
#xlim([-0.5, 13.5])
#legend(fontsize="small")
xlabel("$i-W1$ (magnitude)")
ylabel("Fraction of galaxies detected\nby LoTSS (per cent)")
# xt = ["$(-\infty, 0.0)$",
# "$[0.0, 0.5)$",
# "$[0.5, 1.0)$",
# "$[1.0, 1.25)$",
# "$[1.25, 1.5)$",
# "$[1.5, 1.75)$",
# "$[1.75, 2.0)$",
# "$[2.0, 2.25)$",
# "$[2.25, 2.5)$",
# "$[2.5, 2.75)$",
# "$[2.75, 3.0)$",
# "$[3.0, 3.5)$",
# "$[3.5, 4.0)$",
# "$[4.0, \infty)$"]
# xticks(np.arange(14), xt, rotation=90)
# subplots_adjust(left=0.2,
# bottom=0.1,
# right=0.95,
# top=0.9,
# wspace=0.4,
# hspace=0.2)
subplots_adjust(left=0.15,
bottom=0.15,
right=0.95,
top=0.92,
wspace=0.4,
hspace=0.2)
save = True
if save:
plt.savefig("idata/fraction.png")
plt.savefig("idata/fraction.svg")
plt.savefig("idata/fraction.pdf")
plt.savefig("idata/fraction_high.png", dpi=800)
In [75]:
threshold_sel = 0.639
In [76]:
cond_mlr = (pwl['ID_flag'] == 1) & (pwl['Maj'] < 30.)
In [77]:
pwlaux = pwl[cond_mlr].filled()
Matched sources
In [78]:
pwlaux_match = pwlaux[~np.isnan(pwlaux['ML_LR'])]
In [79]:
len(pwlaux_match)
Out[79]:
In [80]:
cond_match = (
~np.isnan(pwlaux_match['lr']) &
(pwlaux_match['lr'] >= threshold_sel) &
(
(pwlaux_match["AllWISE_input"] != "N/A") |
~np.isnan(pwlaux_match['objID_input'])
) &
(
(
(pwlaux_match["AllWISE"] == pwlaux_match["AllWISE_input"]) &
(pwlaux_match["objID"] == pwlaux_match["objID_input"]) &
~np.isnan(pwlaux_match["objID"]) &
(pwlaux_match["AllWISE"] != "N/A")
) |
(
(pwlaux_match["AllWISE"] == pwlaux_match["AllWISE_input"]) &
np.isnan(pwlaux_match["objID"])
) |
(
(pwlaux_match["AllWISE"] == "N/A") &
(pwlaux_match["objID"] == pwlaux_match["objID_input"])
)
)
)
m_m = np.sum(cond_match)
print(m_m)
In [81]:
cond_diffmatch = (
~np.isnan(pwlaux_match['lr']) &
(pwlaux_match['lr'] >= threshold_sel) &
(
(pwlaux_match["AllWISE_input"] != "N/A") |
~np.isnan(pwlaux_match['objID_input'])
) &
(
(
(pwlaux_match["AllWISE"] != pwlaux_match["AllWISE_input"]) |
(pwlaux_match["objID"] != pwlaux_match["objID_input"])
)
)
)
m_dm = np.sum(cond_diffmatch)
print(m_dm)
In [82]:
cond_nomatch = (
np.isnan(pwlaux_match['lr']) |
(pwlaux_match['lr'] < threshold_sel)
)
m_nm = np.sum(cond_nomatch)
print(m_nm)
In [83]:
m_nm + m_dm + m_m
Out[83]:
In [84]:
217070+512+1021
Out[84]:
Non-matched sources
In [85]:
pwlaux_nomatch = pwlaux[np.isnan(pwlaux['ML_LR'])]
In [86]:
len(pwlaux_nomatch)
Out[86]:
In [87]:
cond2_match = (
~np.isnan(pwlaux_nomatch['lr']) &
(pwlaux_nomatch['lr'] >= threshold_sel)
)
m2_m = np.sum(cond2_match)
print(m2_m)
In [88]:
cond2_nomatch = (
np.isnan(pwlaux_nomatch['lr']) |
(pwlaux_nomatch['lr'] < threshold_sel)
)
m2_nm = np.sum(cond2_nomatch)
print(m2_nm)
In [89]:
m2_nm + m2_m
Out[89]:
In [90]:
m2_nm + m2_m + m_nm + m_dm + m_m
Out[90]:
Diagnostic columns
In [91]:
pwl['match_code'] = 0
In [92]:
pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & ~np.isnan(pwl['ML_LR'])][cond_diffmatch]["Source_Name"])] = 2
In [93]:
pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & ~np.isnan(pwl['ML_LR'])][cond_match]["Source_Name"])]=1
In [94]:
pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & ~np.isnan(pwl['ML_LR'])][cond_nomatch]["Source_Name"])] = 3
In [95]:
pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & np.isnan(pwl['ML_LR'])][cond2_match]["Source_Name"])] = 4
In [96]:
pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & np.isnan(pwl['ML_LR'])][cond2_nomatch]["Source_Name"])] = 5
In [97]:
for i in range(6):
print(i, np.sum(pwl['match_code'] == i))
3 sources that are in group 1 and 2
In [47]:
pwl['match_code2'] = 0
pwl['match_code2'][np.isin(pwl["Source_Name"], pwl[cond_mlr & ~np.isnan(pwl['ML_LR'])][cond_match]["Source_Name"])]=1
In [48]:
t = pwl[np.isin(pwl["Source_Name"], pwl[cond_mlr & ~np.isnan(pwl['ML_LR'])][cond_diffmatch]["Source_Name"])]
t[t['match_code2'] != 0][['Source_Name', "AllWISE", "AllWISE_input", "objID", "objID_input"]]
Out[48]:
In [49]:
pwl[pwl["objID"] == 164861629587942860]
Out[49]:
In [50]:
most_common(pwl["AllWISE_input"].filled(), n=10)
Save data for tests
In [51]:
np.sum(~np.isnan(pwl["colour"]) & (pwl["match_code"] == 3))
Out[51]:
In [52]:
subplot(1,2,1)
val, bins, _ = hist(pwl["colour"][~np.isnan(pwl["colour"]) & (pwl["match_code"] != 0)],
bins=50, normed=True, alpha=0.5, label="Total")
val, bins, _ = hist(pwl["colour"][~np.isnan(pwl["colour"]) & (pwl["match_code"] == 4)],
bins=bins, normed=True, alpha=0.5, label="New matches")
xlabel("Colour")
ylabel("N normed")
legend()
subplot(1,2,2)
val, bins, _ = hist(pwl["colour"][~np.isnan(pwl["colour"]) & (pwl["match_code"] != 0)],
bins=bins, normed=True, alpha=0.5, label="Total")
val, bins, _ = hist(pwl["colour"][~np.isnan(pwl["colour"]) & (pwl["match_code"] == 2)],
bins=bins, normed=True, alpha=0.5, label="Different match")
xlabel("Colour")
ylabel("N normed")
legend()
Out[52]:
In [69]:
describe(pwlaux['ML_LR'])
In [70]:
describe(pwlaux['lr'])
In [71]:
len(pwlaux)
Out[71]:
In [72]:
np.sum(
(pwlaux["AllWISE"] != pwlaux["AllWISE_input"]) |
(pwlaux["objID"] != pwlaux["objID_input"])
)
Out[72]:
In [73]:
np.sum(
(pwlaux["AllWISE"] != pwlaux["AllWISE_input"]) &
(pwlaux["objID"] != pwlaux["objID_input"])
)
Out[73]:
In [74]:
np.sum(
(pwl[cond_mlr].filled()["AllWISE"] != pwl[cond_mlr]["AllWISE_input"]) &
(pwl[cond_mlr].filled()["objID"] == pwl[cond_mlr]["objID_input"])
)
Out[74]:
In [75]:
np.sum(
(pwl[cond_mlr].filled()["AllWISE"] == pwl[cond_mlr]["AllWISE_input"]) &
(pwl[cond_mlr].filled()["objID"] != pwl[cond_mlr]["objID_input"])
)
Out[75]:
In [76]:
for i in ["AllWISE", "AllWISE_input", "objID", "objID_input"]:
print(i)
most_common(pwl[cond_mlr][i].filled(), n=3)
In [77]:
Q_0_colour
Out[77]:
In [78]:
np.sum(Q_0_colour)
Out[78]:
In [ ]: