<jupyter_start><jupyter_text>ML match for LOFAR and the combined PanSTARRS WISE catalogue: Source catalogue  Configuration Load libraries and setup<jupyter_code>import numpy as np
from astropy.table import Table, hstack
from astropy import units as u
from astropy.coordinates import SkyCoord, search_around_sky
from IPython.display import clear_output
from scipy.optimize import curve_fit
import pickle
import os
from pymoc import MOC
from mltier import (get_center, get_n_m, estimate_q_m, Field, SingleMLEstimator, MultiMLEstimator,
                     parallel_process, get_sigma, get_q_m, get_threshold, q0_min_level, q0_min_numbers)
%load_ext autoreload
%autoreload
%pylab inline<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>General configuration<jupyter_code>save_intermediate = True
plot_intermediate = True
idp = "idata/main"
if not os.path.isdir(idp):
    os.makedirs(idp)<jupyter_output><empty_output><jupyter_text>Area limits<jupyter_code>ra_down0 = 170.0
ra_up0 = 190.0
dec_down0 = 45.5
dec_up0 = 56.5
# Busy week Edinburgh 2017
ra_down = 172.09
ra_up = 187.5833
dec_down = 46.106
dec_up = 56.1611
# Busy week Hatfield 2017
ra_down = 170.
ra_up = 190.
dec_down = 46.8
dec_up = 55.9
# Lockman Swire deep region south
ra_down = 150
ra_up = 170
dec_down = 50
dec_up = 70

#load the moc that describes the area where the maximum likelihood matching will occur
#MOC_loc = '/lustre/scratch/astro/im281/data/Lockman_swire_deep/dmu13-lockman-hole-1_MOC.fits'
MOC_loc = '../dmu2/dmu2_field_coverages/Lockman-SWIRE_MOC.fits'
MHz150_MOC = MOC(filename=MOC_loc)
MHz150_MOC.area* (180/np.pi)**2
field = Field(ra_down, ra_up, dec_down, dec_up, MHz150_MOC)
field_no_moc = Field(ra_down, ra_up, dec_down, dec_up)<jupyter_output>field is using a MOC
290483669.7022219
field is not using MOC
2578860555.96<jupyter_text>Load data<jupyter_code>lofar_150mhz = Table.read('data/deep-lockman-v0.1.cat.fits')
len(lofar_150mhz)
np.array(lofar_150mhz.colnames)
#irac_data = Table.read('/users/ianmccheyne/Documents/phd/data/master_catalogue_lockman-swire_iracall_20170710.fits')
irac_data = Table.read('../dmu1/dmu1_ml_Lockman-SWIRE/data/master_catalogue_lockman-swire_20170710.fits')
len(irac_data)
np.array(irac_data.colnames)
len(irac_data['f_irac1'][irac_data['f_irac1'] == np.nan])
a = irac_data['f_irac1'][0]
len(irac_data['f_irac1'][irac_data['f_irac1'] == a])
b = irac_data['f_irac1'][1]
irac0 = []
iracnon0 = []
for n in range(len(irac_data['f_irac1'])):
    total = 0
    if np.isnan(irac_data['f_irac1'][n]) == True:
        if np.isnan(irac_data['f_irac2'][n]) == False:
            total = total + irac_data['f_irac2'][n]
        if np.isnan(irac_data['f_irac3'][n]) == False:
            total = total + irac_data['f_irac3'][n]
        if np.isnan(irac_data['f_irac4'][n]) == False:
            total = total + irac_data['f_irac4'][n]
        irac0.append(total)
    
    else:
        total = total + irac_data['f_irac1'][n]
        if np.isnan(irac_data['f_irac2'][n]) == False:
            total = total + irac_data['f_irac2'][n]
        if np.isnan(irac_data['f_irac3'][n]) == False:
            total = total + irac_data['f_irac3'][n]
        if np.isnan(irac_data['f_irac4'][n]) == False:
            total = total + irac_data['f_irac4'][n]
        iracnon0.append(total)
        
for n in range(len(irac0)):
    if irac0[n] <=0:
        irac0[n] = 0.01
for n in range(len(iracnon0)):
    if iracnon0[n] <=0:
        iracnon0[n] = 0.01
plt.hist(np.log(irac0),bins=100)
plt.hist(np.log(iracnon0),bins = 100,alpha= 0.3)
plt.xlim(-3,10)
plt.ylim(0,100000)
plt.xlabel('log(combined irac flux in all four bands)')<jupyter_output><empty_output><jupyter_text>Filter catalogues<jupyter_code>lofar = field.filter_catalogue(lofar_150mhz, colnames=('RA','DEC'))
lofar = lofar[lofar['Maj']<30]
len(lofar)
# currently filters in a rectangular area around the sources and does not 
#use a moc as the filter in a moc function wasn't working for some reason
combined = field.filter_catalogue(irac_data, colnames=('ra','dec'))
len(combined)
plt.scatter(irac_data['ra'],irac_data['dec'],s=0.1)
plt.scatter(lofar_150mhz['RA'],lofar_150mhz['DEC'],s=0.1,alpha= 0.4)
#plt.xlim(150,170)
#plt.ylim(50,70)
plt.scatter(combined['ra'],combined['dec'],s=0.1)
plt.scatter(lofar['RA'],lofar['DEC'],s=0.1)
#plt.xlim(159.5,165)
#plt.ylim(55,63)<jupyter_output><empty_output><jupyter_text>Additional data  Sky coordinates<jupyter_code>coords_combined = SkyCoord(combined['ra'], 
                           combined['dec'], 
                           unit=(u.deg, u.deg), 
                           frame='icrs')
coords_lofar = SkyCoord(lofar['RA'], 
                       lofar['DEC'], 
                       unit=(u.deg, u.deg), 
                       frame='icrs')<jupyter_output><empty_output><jupyter_text>Class of sources in the combined catalogueThe sources are grouped depending on the available photometric data.  Description  Sky coverage<jupyter_code>plt.rcParams["figure.figsize"] = (15,5)
plot(lofar["RA"],
     lofar["DEC"],
     ls="", marker=",");
plt.rcParams["figure.figsize"] = (8,5)
plot(lofar["RA"],
     lofar["DEC"],
     ls="", marker=",");
len(lofar)
combined["f_irac1"][np.isnan(combined["f_irac1"])] = 0
sum(combined['f_irac1'])<jupyter_output><empty_output><jupyter_text>Summary of galaxy types in the combined catalogue<jupyter_code>plt.rcParams["figure.figsize"] = (15,5)
subplot(1,3,1)
hist(combined["f_irac1"], bins=np.arange(0,1000,1))
xlabel("i")
#subplot(1,3,2)
#hist(combined["W1mag"][combined_matched], bins=50)
#xlabel("W1")
#subplot(1,3,3)
#hist((combined["i"] - combined["W1mag"])[combined_matched], bins=50)
#xlabel("(i - W1)");<jupyter_output><empty_output><jupyter_text>Maximum Likelihood 1st  i-band preparation<jupyter_code>catalogue_i = combined
bin_list_i = np.linspace(0., 200., 401) # Bins of 0.5
center_i = get_center(bin_list_i)
# get the mid point of each bin i.e if there is a bin from 1 to 4 then the centre will be at 2.5
np.percentile(catalogue_i['f_irac1'],99)
#although the field was filtered in a rectangle around the sources the area is calcukated from the moc so is accurate
n_m_i = get_n_m(catalogue_i["f_irac1"], bin_list_i, field.area)
# get a cumulative distribution of the magnitudes divided by the field area
plt.rcParams["figure.figsize"] = (5,5)
plot(center_i, n_m_i);
q_m_i = estimate_q_m(catalogue_i["f_irac1"], bin_list_i, n_m_i, coords_lofar, coords_combined, radius=5)
plt.rcParams["figure.figsize"] = (5,5)
plot(center_i, q_m_i);<jupyter_output><empty_output><jupyter_text>$Q_0$ and likelihood estimators<jupyter_code># Old
Q0_i = 0.545
Q0_w1 = 0.77
# Better Q0
Q0_i = 0.877
Q0_w1 = 0.699
q0 = Q0_i
likelihood_ratio_i = SingleMLEstimator(Q0_i, n_m_i, q_m_i, center_i)<jupyter_output><empty_output><jupyter_text>i-band match<jupyter_code>radius = 10
idx_lofar, idx_i, d2d, d3d = search_around_sky(
    coords_lofar, coords_combined, radius*u.arcsec)
idx_lofar_unique = np.unique(idx_lofar)
idx_i[idx_lofar == 3]
idx_lofar_unique
lofar["lr_i"] = np.nan                   # Likelihood ratio
lofar["lr_dist_i"] = np.nan              # Distance to the selected source
lofar["lr_index_i"] = np.nan             # Index of the PanSTARRS source in combined
total_sources = len(idx_lofar_unique)
combined_aux_index = np.arange(len(combined))
def ml(i):

#    print(i)
#    print(idx_lofar==i)
    
    idx_0 = idx_i[idx_lofar == i]
#    print(idx_0)
    d2d_0 = d2d[idx_lofar == i]
#    print(d2d_0)
#    i_mag = catalogue_i["flux_aper_1"][idx_0]
#    print(i_mag)
#    i_mag = catalogue_i["flux_aper_2"][idx_0]
#    print(i_mag)
#    i_mag = catalogue_i["flux_aper_3"][idx_0]
#    print(i_mag)
#    i_mag = catalogue_i["flux_aper_4"][idx_0]
#    print(i_mag)
#    i_mag = catalogue_i["flux_aper_5"][idx_0]
#    print(i_mag)
    i_mag = catalogue_i["f_irac1"][idx_0]
    
    
    lofar_ra = lofar[i]["RA"]
    lofar_dec = lofar[i]["DEC"]
    lofar_pa = lofar[i]["PA"]
    lofar_maj_err = lofar[i]["E_Maj"]
    lofar_min_err = lofar[i]["E_Min"]
    c_ra = catalogue_i["ra"][idx_0]
    c_dec = catalogue_i["dec"][idx_0]
#    c_ra_err = catalogue_i["raErr"][idx_0]
    c_ra_err = 1/3600
#    c_dec_err = catalogue_i["decErr"][idx_0]
    c_dec_err = 1/3600
    
    sigma = get_sigma(lofar_maj_err, lofar_min_err, lofar_pa, 
                      lofar_ra, lofar_dec, 
                      c_ra, c_dec, c_ra_err, c_dec_err)
    
    lr_0 = likelihood_ratio_i(i_mag, d2d_0.arcsec, sigma)
#    print (lr_0)
    
    if len(lr_0) == 0:
        result = [nan,nan,nan,nan]
        return result
    chosen_index = np.argmax(lr_0)
    result = [combined_aux_index[idx_0[chosen_index]], # Index
              (d2d_0.arcsec)[chosen_index],                        # distance
              lr_0[chosen_index],(i_mag[chosen_index])]                                  # LR

    return result
res = []
nomatches = 0
for i in range(len(lofar)):
    result = ml(i)
    if result == [nan,nan,nan,nan]:
        nomatches = nomatches + 1
    res.append(result)

print(nomatches)
print(len(res))
print(len(lofar["lr_index_i"][idx_lofar_unique]))
(lofar["lr_index_i"][idx_lofar_unique], 
 lofar["lr_dist_i"][idx_lofar_unique], 
 lofar["lr_i"][idx_lofar_unique],
 lofar['lr_i_flux']) = list(map(list, zip(*res)))<jupyter_output><empty_output><jupyter_text>Threshold and selection for i-band<jupyter_code>lofar["lr_i"][np.isnan(lofar["lr_i"])] = 0
#threshold is taken for sources whose likelihood is greater than 0
threshold_i = np.percentile(lofar[lofar['lr_i'] !=0]["lr_i"], 100*(1 - Q0_i))
#threshold_i = np.percentile(["lr_i"], 100*(1 - Q0_i))
Q0_i
threshold_i
plt.rcParams["figure.figsize"] = (15,6)
subplot(1,2,1)
hist(lofar[lofar["lr_i"] != 0]["lr_i"], bins=200)
vlines([threshold_i], 0, 1000)
ylim([0,200])
subplot(1,2,2)
hist(np.log10(lofar[lofar["lr_i"] != 0]["lr_i"]+1), bins=200)
vlines(np.log10(threshold_i+1), 0, 1000)
ticks, _ = xticks()
xticks(ticks, ["{:.1f}".format(10**t-1) for t in ticks])
ylim([0,200]);
lofar["lrt"] = lofar["lr_i"]
lofar["lrt"][np.isnan(lofar["lr_i"])] = 0
def completeness(lr, threshold, q0):
    n = len(lr)
    lrt = lr[lr < threshold]
    return 1. - np.sum((q0 * lrt)/(q0 * lrt + (1 - q0)))/float(n)

def reliability(lr, threshold, q0):
    n = len(lr)
    lrt = lr[lr > threshold]
    return 1. - np.sum((1. - q0)/(q0 * lrt + (1 - q0)))/float(n)

completeness_v = np.vectorize(completeness, excluded=[0])
reliability_v = np.vectorize(reliability, excluded=[0])
n_test = 100
threshold_mean = np.percentile(lofar["lrt"], 100*(1 - q0))
thresholds = np.arange(0., 10., 0.01)
thresholds_fine = np.arange(0.1, 1., 0.001)
completeness_t = completeness_v(lofar["lrt"], thresholds, q0)
reliability_t = reliability_v(lofar["lrt"], thresholds, q0)
average_t = (completeness_t + reliability_t)/2
completeness_t_fine = completeness_v(lofar["lrt"], thresholds_fine, q0)
reliability_t_fine = reliability_v(lofar["lrt"], thresholds_fine, q0)
average_t_fine = (completeness_t_fine + reliability_t_fine)/2
threshold_sel = thresholds_fine[np.argmax(average_t_fine)]
plt.rcParams["figure.figsize"] = (15,6)
subplot(1,2,1)
plot(thresholds, completeness_t, "r-")
plot(thresholds, reliability_t, "g-")
plot(thresholds, average_t, "k-")
vlines(threshold_sel, 0.9, 1., "k", linestyles="dashed")
vlines(threshold_mean, 0.9, 1., "y", linestyles="dashed")
ylim([0.9, 1.])
xlabel("Threshold")
ylabel("Completeness/Reliability")
subplot(1,2,2)
plot(thresholds_fine, completeness_t_fine, "r-")
plot(thresholds_fine, reliability_t_fine, "g-")
plot(thresholds_fine, average_t_fine, "k-")
vlines(threshold_sel, 0.9, 1., "k", linestyles="dashed")
#vlines(threshold_mean, 0.9, 1., "y", linestyles="dashed")
ylim([0.97, 1.])
xlabel("Threshold")
ylabel("Completeness/Reliability")
plt.rcParams["figure.figsize"] = (15,6)
plot(completeness_t,reliability_t)
print(threshold_sel)
plt.rcParams["figure.figsize"] = (15,6)
subplot(1,2,1)
hist(lofar[lofar["lrt"] != 0]["lrt"], bins=200)
vlines([threshold_sel], 0, 500)
ylim([0,300])
subplot(1,2,2)
hist(np.log10(lofar[lofar["lrt"] != 0]["lrt"]+1), bins=200)
vlines(np.log10(threshold_sel+1), 0, 500)
ticks, _ = xticks()
xticks(ticks, ["{:.1f}".format(10**t-1) for t in ticks])
ylim([0,300]);<jupyter_output><empty_output><jupyter_text>create a table containing the lofar sources whose likelihood matches were above the new threshold and create a table of the data for the counterparts and then merge then using hstack<jupyter_code>lofar_match = lofar[~np.isnan(lofar['lr_dist_i'])]
lofar_match = lofar_match[lofar_match['lr_i']>threshold_sel]
len(lofar_match)
lockman_match = Table(rows=combined[int(lofar_match['lr_index_i'][0])],names=combined.colnames)
for index in lofar_match['lr_index_i'][1:]:
    lockman_match.add_row(combined[int(index)])

new_table = hstack([lofar_match,lockman_match],join_type='outer')
Table.write(new_table,'data/merged_lofar_masterlist_lockman-swire.fits')
photoz = Table.read('../dmu24/dmu24_Lockman-SWIRE/data/master_catalogue_Lockman-SWIRE_20170710_photoz_20170802_r_and_irac1_optimised.fits')
photoz[:8]
plt.scatter(photoz['z1_area'],photoz['z2_area'],s=0.2)
plt.plot([0,1],[1,0])
plt.plot([0,1],[0,1])
plt.xlim(0,1)
plt.ylim(0,0.5)
photoz.colnames
c_index = lofar['lr_index_i']
help_id = []
for n in c_index:
    if np.isnan(n)==True:
        continue
    help_id.append(catalogue_i['help_id'][int(n)])
len(help_id)
photoz_filter = field.filter_catalogue(photoz, colnames=('RA','DEC'))
photoz_id = photoz['help_id']
len(photoz_id)
sel = np.in1d(photoz_id,help_id)
len(sel)
photoz_match = photoz[sel]<jupyter_output><empty_output><jupyter_text>number of matches with photoz data<jupyter_code>photoz_match_id = photoz_match['help_id']
len(photoz_match)
xid = Table.read('/Users/ianmccheyne/Documents/phd/data/dmu26_XID+SPIRE_CDFS-SWIRE_cat_20170919.fits')
xid.colnames
xid_match = field.filter_catalogue(xid, colnames=('RA','Dec'))
xid_match_id = xid_match['HELP_ID']
print(len(xid_match_id))
print(xid_match_id[1])
print(photoz_match_id[1])
sel = np.in1d(xid_id,help_id)
sum(sel)
c_index = lofar['lr_index_i']
irac_flux1 = []
irac_flux2 = []
irac_flux3 = []
irac_flux4 = []
lofar_flux = []
for n in range(len(lofar['lrt'])):
    irac_flux1.append(catalogue_i["f_irac1"][int(c_index[n])])
    irac_flux2.append(catalogue_i["f_irac2"][int(c_index[n])])
    irac_flux3.append(catalogue_i["f_irac3"][int(c_index[n])])
    irac_flux4.append(catalogue_i["f_irac4"][int(c_index[n])])
    lofar_flux.append(lofar['Peak_flux'][n])
irac_flux1 = np.array(irac_flux1)
irac_flux2 = np.array(irac_flux2)
irac_flux3 = np.array(irac_flux3)
irac_flux4 = np.array(irac_flux4)
lofar_flux = []
print(len(irac_flux1[np.isnan(irac_flux1)]))
print(len(irac_flux2[np.isnan(irac_flux2)]))
print(len(irac_flux3[np.isnan(irac_flux3)]))
print(len(irac_flux4[np.isnan(irac_flux4)]))
test1 = lofar['lrt'] < threshold_i
test2 = lofar['lrt'] < threshold_sel
lofar_flux = np.array(lofar_flux)
irac_flux = np.array(lofar['lr_i_flux'])
thing2 = irac_flux2[test2]
print(len(irac_flux1[test2]))
print(len(thing2[np.isnan(thing2)]))
thing3 = irac_flux3[test2]
print(len(thing3[np.isnan(thing3)]))
thing4 = irac_flux4[test2]
print(len(thing4[np.isnan(thing4)]))
and23 = 0
and24 = 0
and34 = 0
and234 = 0
for n in range(len(thing2)):
    if np.isnan(thing2[n]) == True and np.isnan(thing3[n]) == True:
        and23 = and23 + 1
    if np.isnan(thing2[n]) == True and np.isnan(thing4[n]) == True:
        and24 = and24 + 1
    if np.isnan(thing3[n]) == True and np.isnan(thing4[n]) == True:
        and34 = and34 + 1
    if np.isnan(thing2[n]) == True and np.isnan(thing3[n]) == True and np.isnan(thing4[n]):
        and234 = and234 + 1
print(and23)
print(and24)
print(and34)
print(and234)
and23 = 0
and24 = 0
and34 = 0
and234 = 0
for n in range(len(irac_flux1)):
    if np.isnan(irac_flux2[n]) == True and np.isnan(irac_flux3[n]) == True:
        and23 = and23 + 1
    if np.isnan(irac_flux2[n]) == True and np.isnan(irac_flux4[n]) == True:
        and24 = and24 + 1
    if np.isnan(irac_flux3[n]) == True and np.isnan(irac_flux4[n]) == True:
        and34 = and34 + 1
    if np.isnan(irac_flux2[n]) == True and np.isnan(irac_flux3[n]) == True and np.isnan(irac_flux4[n]):
        and234 = and234 + 1
print(and23)
print(and24)
print(and34)
print(and234)
lofar_below_original = lofar_flux[(test1)]
lofar_below_new = lofar_flux[lofar['lrt'] > threshold_sel]
irac_below_original = irac_flux[(test1)]
irac_below_new = irac_flux[lofar['lrt'] > threshold_sel]
plt.scatter(np.log(lofar_flux),np.log(irac_flux),s=0.4)
plt.scatter(np.log(lofar_below_original),np.log(irac_below_original),s =0.2)
#plt.scatter(np.log(lofar_below_new),np.log(irac_below_new),s=0.4)
#xlim(-0.001,0.1)
#ylim(-15,500)
plt.scatter(np.log(lofar_flux),np.log(irac_flux),s=1)
#plt.scatter(np.log(lofar_below_original),np.log(irac_below_original),s =0.2)
plt.scatter(np.log(lofar_below_new),np.log(irac_below_new),s=0.4)
ratio = irac_flux/lofar_flux
for n in range(len(ratio)):
    if ratio[n] <=0:
        ratio[n] = 0.00001

plt.hist(np.log(ratio),bins=100)
plt.hist(np.log(ratio[lofar['lrt'] < threshold_i]),bins=100)
plt.hist(np.log(ratio[lofar['lrt'] < threshold_sel]),bins=100)
plt.xlim(0,20)
plt.xlabel('log( irac flux / lofar flux)')
print(2)
plt.hist(np.log(ratio[lofar['lrt'] < threshold_sel]),bins=100)
plt.plot(np.arange(0.0,20.0,0.01),gauss(np.arange(0.0,20.0,0.01),16.307,0.224,250))
ratio_hist = np.histogram(np.log(ratio),bins=100)
ratio_hist_old = np.histogram(np.log(ratio[lofar['lrt'] < threshold_i]),bins=100)
ratio_hist_new = np.histogram(np.log(ratio[lofar['lrt'] < threshold_sel]),bins=100)
ratio_hist[0][0] = 0
ratio_hist_old[0][0] = 0
ratio_hist_new[0][0] = 0
def gauss(x,mean,sig):
    amp = 1 / (sig*np.sqrt(2*3.1419))
    exp = -0.5 * (x-mean)**2/(sig**2)
    return(amp * np.exp(exp))
popt1, pcov1 = curve_fit(gauss,ratio_hist[1][:99],ratio_hist[0][:99])
popt2, pcov2 = curve_fit(gauss,ratio_hist_old[1][:99],ratio_hist_old[0][:99])
popt3, pcov3 = curve_fit(gauss,ratio_hist_new[1][:99],ratio_hist_new[1][:99])
print(popt1)
print(popt2)
print(popt3)
max(lofar['lr_i'])
percent = np.linspace(0,60.531,201)
hist = np.histogram(lofar[lofar['lr_i'] !=0]['lr_i'], bins=201)
#print (hist)
cum = np.cumsum(hist[0])
print (cum)
plt.plot(percent,cum/cum[-1])
plt.xlabel('likelihood')
plt.ylabel('cumulative fraction')
lofar["lr_index_sel_i"] = lofar["lr_index_i"]
lofar["lr_index_sel_i"][lofar["lr_i"] < threshold_i] = np.nan<jupyter_output><empty_output><jupyter_text>W1-band match<jupyter_code>idx_lofar, idx_i, d2d, d3d = search_around_sky(
    coords_lofar, coords_combined[combined_w1], radius*u.arcsec)
idx_lofar_unique = np.unique(idx_lofar)
lofar["lr_w1"] = np.nan                   # Likelihood ratio
lofar["lr_dist_w1"] = np.nan              # Distance to the selected source
lofar["lr_index_w1"] = np.nan             # Index of the PanSTARRS source in combined
def ml_w1(i):
    idx_0 = idx_i[idx_lofar == i]
    d2d_0 = d2d[idx_lofar == i]
    w1_mag = catalogue_w1["W1mag"][idx_0]
    
    lofar_ra = lofar[i]["RA"]
    lofar_dec = lofar[i]["DEC"]
    lofar_pa = lofar[i]["PA"]
    lofar_maj_err = lofar[i]["E_Maj"]
    lofar_min_err = lofar[i]["E_Min"]
    c_ra = catalogue_w1["ra"][idx_0]
    c_dec = catalogue_w1["dec"][idx_0]
    c_ra_err = catalogue_w1["raErr"][idx_0]
    c_dec_err = catalogue_w1["decErr"][idx_0]
    
    sigma = get_sigma(lofar_maj_err, lofar_min_err, lofar_pa, 
                      lofar_ra, lofar_dec, 
                      c_ra, c_dec, c_ra_err, c_dec_err)
    
    lr_0 = likelihood_ratio_w1(w1_mag, d2d_0.arcsec, sigma)
    chosen_index = np.argmax(lr_0)
    result = [combined_aux_index[combined_w1][idx_0[chosen_index]], # Index
              (d2d_0.arcsec)[chosen_index],                        # distance
              lr_0[chosen_index]]                                  # LR
    return result
res = parallel_process(idx_lofar_unique, ml_w1)
(lofar["lr_index_w1"][idx_lofar_unique], 
 lofar["lr_dist_w1"][idx_lofar_unique], 
 lofar["lr_w1"][idx_lofar_unique]) = list(map(list, zip(*res)))<jupyter_output><empty_output><jupyter_text>Threshold and selection for W1 band<jupyter_code>lofar["lr_w1"][np.isnan(lofar["lr_w1"])] = 0
threshold_w1 = np.percentile(lofar["lr_w1"], 100*(1 - Q0_w1))
threshold_w1
plt.rcParams["figure.figsize"] = (15,6)
subplot(1,2,1)
hist(lofar[lofar["lr_w1"] != 0]["lr_w1"], bins=200)
vlines([threshold_i], 0, 1000)
ylim([0,1000])
subplot(1,2,2)
hist(np.log10(lofar[lofar["lr_w1"] != 0]["lr_w1"]+1), bins=200)
vlines(np.log10(threshold_i+1), 0, 1000)
ticks, _ = xticks()
xticks(ticks, ["{:.1f}".format(10**t-1) for t in ticks])
ylim([0,1000]);
lofar["lr_index_sel_w1"] = lofar["lr_index_w1"]
lofar["lr_index_sel_w1"][lofar["lr_w1"] < threshold_w1] = np.nan<jupyter_output><empty_output><jupyter_text>Final selection of the matchWe combine the ML matching done in i-band and W1-band. All the galaxies were the LR is above the selection ratio for the respective band are finally selected.<jupyter_code># lr_i_and_w1 = (lofar["lr_i"] != 0) & (lofar["lr_w1"] != 0)
# lr_only_i = (lofar["lr_i"] != 0) & (lofar["lr_w1"] == 0)
# lr_only_w1 = (lofar["lr_i"] == 0) & (lofar["lr_w1"] != 0)
# lr_no_match = (lofar["lr_i"] == 0) & (lofar["lr_w1"] == 0)
lr_i_and_w1 = ~np.isnan(lofar["lr_index_sel_i"]) & ~np.isnan(lofar["lr_index_sel_w1"])
lr_only_i = ~np.isnan(lofar["lr_index_sel_i"]) & np.isnan(lofar["lr_index_sel_w1"])
lr_only_w1 = np.isnan(lofar["lr_index_sel_i"]) & ~np.isnan(lofar["lr_index_sel_w1"])
lr_no_match = np.isnan(lofar["lr_index_sel_i"]) & np.isnan(lofar["lr_index_sel_w1"])
print(np.sum(lr_i_and_w1))
print(np.sum(lr_only_i))
print(np.sum(lr_only_w1))
print(np.sum(lr_no_match))
lofar["lr_index_1"] = np.nan
lofar["lr_dist_1"] = np.nan
lofar["lr_1"] = np.nan
lofar["lr_type_1"] = 0
# Only i matches
lofar["lr_1"][lr_only_i] = lofar["lr_i"][lr_only_i]
lofar["lr_index_1"][lr_only_i] = lofar["lr_index_i"][lr_only_i]
lofar["lr_dist_1"][lr_only_i] = lofar["lr_dist_i"][lr_only_i]
lofar["lr_type_1"][lr_only_i] = 1

# Only w1 matches
lofar["lr_1"][lr_only_w1] = lofar["lr_w1"][lr_only_w1]
lofar["lr_index_1"][lr_only_w1] = lofar["lr_index_w1"][lr_only_w1]
lofar["lr_dist_1"][lr_only_w1] = lofar["lr_dist_w1"][lr_only_w1]
lofar["lr_type_1"][lr_only_w1] = 2

# Both matches
lofar["lr_1"][lr_i_and_w1] = np.max([lofar["lr_i"][lr_i_and_w1], lofar["lr_w1"][lr_i_and_w1]], axis=0)
lofar["lr_type_1"][lr_i_and_w1] = np.argmax([lofar["lr_i"][lr_i_and_w1], lofar["lr_w1"][lr_i_and_w1]], axis=0) + 1

c1 = (lofar["lr_type_1"] == 1)
c2 = (lofar["lr_type_1"] == 2)
lofar["lr_index_1"][lr_i_and_w1 & c1] = lofar["lr_index_i"][lr_i_and_w1 & c1]
lofar["lr_index_1"][lr_i_and_w1 & c2] = lofar["lr_index_w1"][lr_i_and_w1 & c2]
lofar["lr_dist_1"][lr_i_and_w1 & c1] = lofar["lr_dist_i"][lr_i_and_w1 & c1]
lofar["lr_dist_1"][lr_i_and_w1 & c2] = lofar["lr_dist_w1"][lr_i_and_w1 & c2]<jupyter_output><empty_output><jupyter_text>Summary of the number of sources matches of each type<jupyter_code>print("match    sel-i: ", np.sum(lofar["lr_type_1"][lr_i_and_w1] == 1))
print("match   sel-W1: ", np.sum(lofar["lr_type_1"][lr_i_and_w1] == 2))
print("match     both: ", np.sum(lofar["lr_type_1"][lr_i_and_w1] == 1) + 
                          np.sum(lofar["lr_type_1"][lr_i_and_w1] == 2))
print("match   i-only: ", np.sum(lofar["lr_type_1"] == 1) - np.sum(lofar["lr_type_1"][lr_i_and_w1] == 1))
print("match  W1-only: ", np.sum(lofar["lr_type_1"] == 2) - np.sum(lofar["lr_type_1"][lr_i_and_w1] == 2))
print("match      all: ", np.sum(lofar["lr_type_1"] == 1) + 
                          np.sum(lofar["lr_type_1"] == 2))
print("         Total: ", len(lofar))<jupyter_output><empty_output><jupyter_text>The number of sources for which the match in i-band and W1-band are above the threshold but gives a different match to the combined catalogue.<jupyter_code>print(np.sum(lofar["lr_index_i"][lr_i_and_w1] != lofar["lr_index_w1"][lr_i_and_w1]))<jupyter_output><empty_output><jupyter_text>Duplicated sourcesThis is the nymber of sources of the combined catalogue that are combined to multiple LOFAR sources. In the case of the catalogue of Gaussians the number can be very high.<jupyter_code>values, counts = np.unique(lofar[lofar["lr_type_1"] != 0]["lr_index_1"], return_counts=True)
len(values[counts > 1])
n_dup, n_sour = np.unique(counts[counts > 1], return_counts=True)
plt.rcParams["figure.figsize"] = (6,6)
semilogy(n_dup, n_sour, marker="x")
xlabel("Number of multiple matches")
ylabel("Number of sources in the category")<jupyter_output><empty_output><jupyter_text>Save intermediate data<jupyter_code>if save_intermediate:
    pickle.dump([bin_list_i, center_i, Q0_i, n_m_i, q_m_i], 
                open("{}/lofar_params_1i.pckl".format(idp), 'wb'))
    pickle.dump([bin_list_w1, center_w1, Q0_w1, n_m_w1, q_m_w1], 
                open("{}/lofar_params_1w1.pckl".format(idp), 'wb'))
    lofar.write("{}/lofar_m1.fits".format(idp), format="fits")<jupyter_output><empty_output><jupyter_text>Second iteration using colourFrom now on we will take into account the effect of the colour. The sample was distributed in several categories according to the colour of the source and this is considered here. Rusable parameters for all the iterationsThese parameters are derived from the underlying population and will not change.First we compute the number of galaxies in each bin for the combined catalogue<jupyter_code>bin_list = [bin_list_w1 if i == 0 else bin_list_i for i in range(len(colour_bin_def))]
centers = [center_w1 if i == 0 else center_i for i in range(len(colour_bin_def))]
numbers_combined_bins = np.array([np.sum(a["condition"]) for a in colour_bin_def])
numbers_combined_bins<jupyter_output><empty_output><jupyter_text>Get the colour category and magnitudes for the matched LOFAR sources<jupyter_code>n_m = []

# W1 only sources
n_m.append(get_n_m(combined["W1mag"][combined["category"] == 0], bin_list_w1, field.area))

# Rest of the sources
for i in range(1, len(colour_bin_def)):
    n_m.append(get_n_m(combined["i"][combined["category"] == i], bin_list_i, field.area))
plt.rcParams["figure.figsize"] = (15,15)
for i, n_m_k in enumerate(n_m):
    subplot(5,5,i+1)
    plot(centers[i], n_m_k)<jupyter_output><empty_output><jupyter_text>Parameters of the matched sampleThe parameters derived from the matched LOFAR galaxies: $q_0$, q(m) and the number of sources per category.The columns "category", "W1mag" and "i" will contain the properties of the matched galaxies and will be updated in each iteration to save space.<jupyter_code>lofar["category"] = np.nan
lofar["W1mag"] = np.nan
lofar["i"] = np.nan
c = ~np.isnan(lofar["lr_index_1"])
indices = lofar["lr_index_1"][c].astype(int)
lofar["category"][c] = combined[indices]["category"]
lofar["W1mag"][c] = combined[indices]["W1mag"]
lofar["i"][c] = combined[indices]["i"]<jupyter_output><empty_output><jupyter_text>The next parameter represent the number of matched LOFAR sources in each colour category.<jupyter_code>numbers_lofar_combined_bins = np.array([np.sum(lofar["category"] == c) 
                                        for c in range(len(numbers_combined_bins))])
numbers_lofar_combined_bins<jupyter_output><empty_output><jupyter_text>The $Q_0$ for each category are obtained by dividing the number of sources in the category by the total number of sources in the sample.<jupyter_code>Q_0_colour = numbers_lofar_combined_bins/len(lofar) ### Q_0
q0_total = np.sum(Q_0_colour)
q0_total<jupyter_output><empty_output><jupyter_text>The q(m) is not estimated with the method of Fleuren et al. but with the most updated distributions and numbers for the matches.<jupyter_code>q_m = []
radius = 15. 

# W1 only sources
q_m.append(get_q_m(lofar["W1mag"][lofar["category"] == 0], 
                   bin_list_w1, 
                   numbers_lofar_combined_bins[0], 
                   n_m[0], 
                   field.area, 
                   radius=radius))

# Rest of the sources
for i in range(1, len(numbers_lofar_combined_bins)):
    q_m.append(get_q_m(lofar["i"][lofar["category"] == i], 
                   bin_list_i, 
                   numbers_lofar_combined_bins[i], 
                   n_m[i], 
                   field.area, 
                   radius=radius))
plt.rcParams["figure.figsize"] = (15,15)
for i, q_m_k in enumerate(q_m):
    subplot(5,5,i+1)
    plot(centers[i], q_m_k)<jupyter_output><empty_output><jupyter_text>Save intermediate parameters<jupyter_code>if save_intermediate:
    pickle.dump([bin_list, centers, Q_0_colour, n_m, q_m], 
                open("{}/lofar_params_2.pckl".format(idp), 'wb'))<jupyter_output><empty_output><jupyter_text>Prepare for ML<jupyter_code>selection = ~np.isnan(combined["category"]) # Avoid the two dreaded sources with no actual data
catalogue = combined[selection]
radius = 15
def apply_ml(i, likelihood_ratio_function):
    idx_0 = idx_i[idx_lofar == i]
    d2d_0 = d2d[idx_lofar == i]
    
    category = catalogue["category"][idx_0].astype(int)
    mag = catalogue["i"][idx_0]
    mag[category == 0] = catalogue["W1mag"][idx_0][category == 0]
    
    lofar_ra = lofar[i]["RA"]
    lofar_dec = lofar[i]["DEC"]
    lofar_pa = lofar[i]["PA"]
    lofar_maj_err = lofar[i]["E_Maj"]
    lofar_min_err = lofar[i]["E_Min"]
    c_ra = catalogue["ra"][idx_0]
    c_dec = catalogue["dec"][idx_0]
    c_ra_err = catalogue["raErr"][idx_0]
    c_dec_err = catalogue["decErr"][idx_0]
    
    sigma = get_sigma(lofar_maj_err, lofar_min_err, lofar_pa, 
                      lofar_ra, lofar_dec, 
                      c_ra, c_dec, c_ra_err, c_dec_err)

    lr_0 = likelihood_ratio_function(mag, d2d_0.arcsec, sigma, category)
    
    chosen_index = np.argmax(lr_0)
    result = [combined_aux_index[selection][idx_0[chosen_index]], # Index
              (d2d_0.arcsec)[chosen_index],                        # distance
              lr_0[chosen_index]]                                  # LR
    return result<jupyter_output><empty_output><jupyter_text>Run the cross-matchThis will not need to be repeated after<jupyter_code>idx_lofar, idx_i, d2d, d3d = search_around_sky(
    coords_lofar, coords_combined[selection], radius*u.arcsec)
idx_lofar_unique = np.unique(idx_lofar)<jupyter_output><empty_output><jupyter_text>Run the ML matching<jupyter_code>likelihood_ratio = MultiMLEstimator(Q_0_colour, n_m, q_m, centers)
def ml(i):
    return apply_ml(i, likelihood_ratio)
res = parallel_process(idx_lofar_unique, ml, n_jobs=4)
lofar["lr_index_2"] = np.nan
lofar["lr_dist_2"] = np.nan
lofar["lr_2"] = np.nan
(lofar["lr_index_2"][idx_lofar_unique], 
 lofar["lr_dist_2"][idx_lofar_unique], 
 lofar["lr_2"][idx_lofar_unique]) = list(map(list, zip(*res)))<jupyter_output><empty_output><jupyter_text>Get the new threshold for the ML matching. FIX THIS<jupyter_code>lofar["lr_2"][np.isnan(lofar["lr_2"])] = 0
threshold = np.percentile(lofar["lr_2"], 100*(1 - q0_total))
#manual_q0 = 0.65
#threshold = np.percentile(lofar["lr_2"], 100*(1 - manual_q0))
threshold
plt.rcParams["figure.figsize"] = (15,6)
subplot(1,2,1)
hist(lofar[lofar["lr_2"] != 0]["lr_2"], bins=200)
vlines([threshold_i], 0, 1000)
ylim([0,1000])
subplot(1,2,2)
hist(np.log10(lofar[lofar["lr_2"] != 0]["lr_2"]+1), bins=200)
vlines(np.log10(threshold_i+1), 0, 1000)
ticks, _ = xticks()
xticks(ticks, ["{:.1f}".format(10**t-1) for t in ticks])
ylim([0,1000]);
lofar["lr_index_sel_2"] = lofar["lr_index_2"]
lofar["lr_index_sel_2"][lofar["lr_2"] < threshold] = np.nan
n_changes = np.sum((lofar["lr_index_sel_2"] != lofar["lr_index_1"]) & 
                   ~np.isnan(lofar["lr_index_sel_2"]) &
                   ~np.isnan(lofar["lr_index_1"]))
n_changes<jupyter_output><empty_output><jupyter_text>Enter the results<jupyter_code># Clear aux columns
lofar["category"] = np.nan
lofar["W1mag"] = np.nan
lofar["i"] = np.nan

c = ~np.isnan(lofar["lr_index_sel_2"])
indices = lofar["lr_index_sel_2"][c].astype(int)
lofar["category"][c] = combined[indices]["category"]
lofar["W1mag"][c] = combined[indices]["W1mag"]
lofar["i"][c] = combined[indices]["i"]<jupyter_output><empty_output><jupyter_text>Save intermediate data<jupyter_code>if save_intermediate:
    lofar.write("{}/lofar_m2.fits".format(idp), format="fits")<jupyter_output><empty_output><jupyter_text>Iterate until convergence<jupyter_code>radius = 15. 
from matplotlib import pyplot as plt
for j in range(10):
    iteration = j+3 
    print("Iteration {}".format(iteration))
    print("=============")
    ## Get new parameters
    # Number of matched sources per bin
    numbers_lofar_combined_bins = np.array([np.sum(lofar["category"] == c) 
                                            for c in range(len(numbers_combined_bins))])
    print("numbers_lofar_combined_bins")
    print(numbers_lofar_combined_bins)
    # q_0
    Q_0_colour_est = numbers_lofar_combined_bins/len(lofar) ### Q_0
    Q_0_colour = q0_min_numbers(Q_0_colour_est, numbers_combined_bins)
    print("Q_0_colour")
    print(Q_0_colour)
    q0_total = np.sum(Q_0_colour)
    print("Q_0_total: ", q0_total)
    # q_m
    q_m = []
    # W1 only sources
    q_m.append(get_q_m(lofar["W1mag"][lofar["category"] == 0], 
                   bin_list_w1, 
                   numbers_lofar_combined_bins[0], 
                   n_m[0], 
                   field.area, 
                   radius=radius))
    # Rest of the sources
    for i in range(1, len(numbers_lofar_combined_bins)):
        q_m.append(get_q_m(lofar["i"][lofar["category"] == i], 
                       bin_list_i, 
                       numbers_lofar_combined_bins[i], 
                       n_m[i], 
                       field.area, 
                       radius=radius))
    # Save new parameters
    if save_intermediate:
        pickle.dump([bin_list, centers, Q_0_colour, n_m, q_m], 
                    open("{}/lofar_params_{}.pckl".format(idp, iteration), 'wb'))
    if plot_intermediate:
        fig = plt.figure(figsize=(15,15))
        for i, q_m_k in enumerate(q_m):
            plt.subplot(5,5,i+1)
            plt.plot(centers[i], q_m_k)
        plt.savefig('{}/q0_{}.png'.format(idp, iteration))
        del fig
    ## Define new likelihood_ratio
    likelihood_ratio = MultiMLEstimator(Q_0_colour, n_m, q_m, centers)
    def ml(i):
        return apply_ml(i, likelihood_ratio)
    ## Run the ML
    res = parallel_process(idx_lofar_unique, ml, n_jobs=4)
    lofar["lr_index_{}".format(iteration)] = np.nan
    lofar["lr_dist_{}".format(iteration)] = np.nan
    lofar["lr_{}".format(iteration)] = np.nan
    (lofar["lr_index_{}".format(iteration)][idx_lofar_unique], 
     lofar["lr_dist_{}".format(iteration)][idx_lofar_unique], 
     lofar["lr_{}".format(iteration)][idx_lofar_unique]) = list(map(list, zip(*res)))
    lofar["lr_{}".format(iteration)][np.isnan(lofar["lr_{}".format(iteration)])] = 0
    ## Get and apply the threshold
    threshold = np.percentile(lofar["lr_{}".format(iteration)], 100*(1 - q0_total))
    #threshold = get_threshold(lofar[lofar["lr_{}".format(iteration)] != 0]["lr_{}".format(iteration)])
    print("Threshold: ", threshold)
    if plot_intermediate:
        fig = plt.figure(figsize=(15,6))
        plt.subplot(1,2,1)
        plt.hist(lofar[lofar["lr_{}".format(iteration)] != 0]["lr_{}".format(iteration)], bins=200)
        plt.vlines([threshold], 0, 1000)
        plt.ylim([0,1000])
        plt.subplot(1,2,2)
        plt.hist(np.log10(lofar[lofar["lr_{}".format(iteration)] != 0]["lr_{}".format(iteration)]+1), bins=200)
        plt.vlines(np.log10(threshold+1), 0, 1000)
        ticks, _ = plt.xticks()
        plt.xticks(ticks, ["{:.1f}".format(10**t-1) for t in ticks])
        plt.ylim([0,1000])
        plt.savefig('{}/lr_distribution_{}.png'.format(idp, iteration))
        del fig
    ## Apply the threshold
    lofar["lr_index_sel_{}".format(iteration)] = lofar["lr_index_{}".format(iteration)]
    lofar["lr_index_sel_{}".format(iteration)][lofar["lr_{}".format(iteration)] < threshold] = np.nan
    ## Enter changes into the catalogue
    # Clear aux columns
    lofar["category"] = np.nan
    lofar["W1mag"] = np.nan
    lofar["i"] = np.nan
    # Update data
    c = ~np.isnan(lofar["lr_index_sel_{}".format(iteration)])
    indices = lofar["lr_index_sel_{}".format(iteration)][c].astype(int)
    lofar["category"][c] = combined[indices]["category"]
    lofar["W1mag"][c] = combined[indices]["W1mag"]
    lofar["i"][c] = combined[indices]["i"]
    # Save the data
    if save_intermediate:
        lofar.write("{}/lofar_m{}.fits".format(idp, iteration), format="fits")
    ## Compute number of changes
    n_changes = np.sum((
            lofar["lr_index_sel_{}".format(iteration)] != lofar["lr_index_sel_{}".format(iteration-1)]) & 
            ~np.isnan(lofar["lr_index_sel_{}".format(iteration)]) &
            ~np.isnan(lofar["lr_index_sel_{}".format(iteration-1)]))
    print("N changes: ", n_changes)
    ## Check changes
    if n_changes == 0:
        break
    else:
        print("******** continue **********")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16140911429422436, 'fraction_numerical': 0.0387750128618234, 'mean_word_length': 4.807831420034597, 'pattern_counts': {'":': 0, '<': 155, '<?xml version=': 0, '>': 145, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '35', 'dir': 'jupyter-structured-clean-dedup', 'id': '298108', 'license': 'MIT', 'n_tokens_mistral': 16161, 'n_tokens_neox': 15427, 'n_words': 2399, 'path': 'dmu8/Match_LOFAR_combined_final.ipynb', 'repo_name': 'djbsmith/dmu_products', 'size': '36931'}