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

from mltier1 import (get_center, Field, parallel_process, describe)

%load_ext autoreload

%autoreload

from IPython.display import clear_output

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) 

%pylab inline

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)])

"""
Explanation: Analysis of the final catalogue of matched sources
We will analyse the changes in the classification using the new sigma and the new catalogue without the galaxies that went to LGZ
Configuration
Load libraries and setup
End of explanation
"""


save_intermediate = True
plot_intermediate = True

idp = "idata/final_analysis_pdf_v1.0"

if not os.path.isdir(idp):
    os.makedirs(idp)

"""
Explanation: General configuration
End of explanation
"""


pwli = Table.read("lofar_pw_pdf.fits")

pwli.colnames

lofar_all = Table.read("data/LOFAR_HBA_T1_DR1_merge_ID_optical_v1.1b.fits")

"""
Explanation: Load data
End of explanation
"""


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'])

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]

bin_list, centers, Q_0_colour, n_m, q_m = pickle.load(open("lofar_params.pckl", "rb"))

"""
Explanation: Join data tables
End of explanation
"""


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

pwl["colour"][(pwl["colour"] == 1e+20)] = np.nan

describe(pwl["colour"])

"""
Explanation: Explore and repair data
End of explanation
"""


pwl["AllWISE_input"][pwl["AllWISE_input"] == ""] = "N/A"

"""
Explanation: Change the AllWISE_input that are '' to 'N/A'. That comes from a previous error with the fill value.
End of explanation
"""


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: Plots for the paper
q(m) over n(m) depending on the colour (not used)
End of explanation
"""


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")

"""
Explanation: Explanation of the KDE used in the computing of the q(m) and n(m):
* https://jakevdp.github.io/blog/2013/12/01/kernel-density-estimation/
Figure 3 of the paper
Help from:
* https://stackoverflow.com/questions/19877666/add-legends-to-linecollection-plot
* https://matplotlib.org/users/mathtext.html#mathtext-tutorial
End of explanation
"""


q0 = np.sum(Q_0_colour)

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])

pwl["lrt"] = pwl["lr"]
pwl["lrt"][np.isnan(pwl["lr"])] = 0

n_test = 100
threshold_mean = np.percentile(pwl["lrt"], 100*(1 - q0))

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

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

thresholds_fine[np.argmax(average_t_fine)]

thresholds_fine[np.argmin(np.abs(completeness_t_fine-reliability_t_fine))]

np.sum(pwl["lrt"] >= 0.358)

np.sum(pwl["lrt"] >= 0.639)

np.sum(pwl["lrt"] >= 0.639)

len(pwl)

len(pwl[pwl['ID_flag'] == 1])

n0 = np.sum((pwl["lrt"] >= 0.639) & (pwl['ID_flag'] == 1))
print(n0, n0/len(pwl[pwl['ID_flag'] == 1]))

n0 = np.sum((pwl["lrt"] >= 0.81268) & (pwl['ID_flag'] == 1))
print(n0, n0/len(pwl[pwl['ID_flag'] == 1]))

n0 = np.sum((pwl["ML_LR"] >= 0.81268) & (pwl['ID_flag'] == 1))
print(n0, n0/len(pwl[pwl['ID_flag'] == 1]))

np.sum(pwl["lrt"] >= 0.639)/len(pwl)

threshold_sel = 0.639

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)

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)

"""
Explanation: Figure 4 of the paper: completeness and reliability
End of explanation
"""


cond_mlr = (pwl['ID_flag'] == 1)

len(pwl)

np.sum(np.isnan(pwl["category"]))

np.sum(np.isnan(pwl["category"][cond_mlr]))

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])

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]
         ))

total_n_c

"""
Explanation: N_LOFAR
End of explanation
"""


field = Field(170.0, 190.0, 46.8, 55.9)

combined = Table.read("pw.fits")

combined = field.filter_catalogue(combined, 
                               colnames=("ra", "dec"))

combined["colour"] = combined["i"] - combined["W1mag"]

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]

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

# 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])})

combined["category"] = np.nan
for i in range(len(colour_bin_def)):
    combined["category"][colour_bin_def[i]["condition"]] = i

numbers_combined_bins = np.array([np.sum(a["condition"]) for a in colour_bin_def])

numbers_combined_bins

pwlf = field.filter_catalogue(pwl[cond_mlr], 
                               colnames=("ra", "dec"))

pwlf_all = field.filter_catalogue(pwl, 
                               colnames=("ra", "dec"))

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]))

np.sum(np.array(lofar_c_all)/np.array(pw_c))

np.sum(np.array(lofar_c_all))/np.sum(np.array(pw_c))

len(pwl)/len(combined)

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)

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)

"""
Explanation: Fractions
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.
End of explanation
"""


threshold_sel = 0.639

cond_mlr = (pwl['ID_flag'] == 1) & (pwl['Maj'] < 30.)

pwlaux = pwl[cond_mlr].filled()

"""
Explanation: Create the additional columns for the types of matches
End of explanation
"""


pwlaux_match = pwlaux[~np.isnan(pwlaux['ML_LR'])]

len(pwlaux_match)

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)

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)

cond_nomatch = (
    np.isnan(pwlaux_match['lr']) |
    (pwlaux_match['lr'] < threshold_sel)
)
m_nm = np.sum(cond_nomatch)
print(m_nm)

m_nm + m_dm + m_m

217070+512+1021

"""
Explanation: Matched sources
End of explanation
"""


pwlaux_nomatch = pwlaux[np.isnan(pwlaux['ML_LR'])]

len(pwlaux_nomatch)

cond2_match = (
    ~np.isnan(pwlaux_nomatch['lr']) &
    (pwlaux_nomatch['lr'] >= threshold_sel)   
)
m2_m = np.sum(cond2_match)
print(m2_m)

cond2_nomatch = (
    np.isnan(pwlaux_nomatch['lr']) |
    (pwlaux_nomatch['lr'] < threshold_sel)
)
m2_nm = np.sum(cond2_nomatch)
print(m2_nm)

m2_nm + m2_m

m2_nm + m2_m + m_nm + m_dm + m_m

"""
Explanation: Non-matched sources
End of explanation
"""


pwl['match_code'] = 0

pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & ~np.isnan(pwl['ML_LR'])][cond_diffmatch]["Source_Name"])] = 2

pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & ~np.isnan(pwl['ML_LR'])][cond_match]["Source_Name"])]=1

pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & ~np.isnan(pwl['ML_LR'])][cond_nomatch]["Source_Name"])] = 3

pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & np.isnan(pwl['ML_LR'])][cond2_match]["Source_Name"])] = 4

pwl['match_code'][np.isin(pwl["Source_Name"], pwl[cond_mlr & np.isnan(pwl['ML_LR'])][cond2_nomatch]["Source_Name"])] = 5

for i in range(6):
    print(i, np.sum(pwl['match_code'] == i))

"""
Explanation: Diagnostic columns
End of explanation
"""


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

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"]]

pwl[pwl["objID"] == 164861629587942860]

"""
Explanation: Study the 3 repeated sources
3 sources that are in group 1 and 2
End of explanation
"""


most_common(pwl["AllWISE_input"].filled(), n=10)

"""
Explanation: Analyse changes in the matches
End of explanation
"""


np.sum(~np.isnan(pwl["colour"]) & (pwl["match_code"] == 3))

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()

"""
Explanation: Save data for tests
End of explanation
"""


describe(pwlaux['ML_LR'])

describe(pwlaux['lr'])

len(pwlaux)

np.sum(
    (pwlaux["AllWISE"] != pwlaux["AllWISE_input"]) |
    (pwlaux["objID"] != pwlaux["objID_input"])  
      )

np.sum(
    (pwlaux["AllWISE"] != pwlaux["AllWISE_input"]) &
    (pwlaux["objID"] != pwlaux["objID_input"])  
      )

np.sum(
    (pwl[cond_mlr].filled()["AllWISE"] != pwl[cond_mlr]["AllWISE_input"]) &
    (pwl[cond_mlr].filled()["objID"] == pwl[cond_mlr]["objID_input"])  
      )

np.sum(
    (pwl[cond_mlr].filled()["AllWISE"] == pwl[cond_mlr]["AllWISE_input"]) &
    (pwl[cond_mlr].filled()["objID"] != pwl[cond_mlr]["objID_input"])  
      )

for i in ["AllWISE", "AllWISE_input", "objID", "objID_input"]:
    print(i)
    most_common(pwl[cond_mlr][i].filled(), n=3)

Q_0_colour

np.sum(Q_0_colour)

"""
Explanation: Additional description of the data
End of explanation
"""


