import os, glob
import fitsio
import numpy as np
import matplotlib.pyplot as plt
from astropy.table import vstack, Table, hstack

import seaborn as sns
sns.set(context='talk', style='ticks', font_scale=1.0)
%matplotlib inline

lslgaver = b'L6'
lslgafile = '/global/cfs/cdirs/desi/users/ioannis/LSLGA-data-DR9fg/LSLGA-model-v6.0.fits'
#lslgafile = '/global/project/projectdirs/cosmo/staging/largegalaxies/v6.0/LSLGA-model-v6.0.fits'
#lslgafile = os.getenv('LARGEGALAXIES_CAT')
lslga = Table(fitsio.read(lslgafile))
print('Read {} galaxies from {}'.format(len(lslga), lslgafile))
#lslga

"""
Explanation: Check the geometry of LSLGA galaxies derived from ellipse-fitting against the original geometry from Hyperleda.
End of explanation
"""


plslga = lslga[lslga['REF_CAT'] == lslgaver]

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(6, 14))
ax1.hexbin(np.log10(plslga['D25']), np.log10(plslga['D25_ORIG']), mincnt=1,
          cmap=plt.cm.Blues_r)
#ax1.scatter(np.log10(lslga['D25']), np.log10(lslga['D25_ORIG']), s=15)
ax1.set_xlim(-0.7, 1)
ax1.set_ylim(-0.7, 1)
ax1.plot([-1, 1], [-1, 1], color='k', ls='--', lw=3, alpha=0.7)
ax1.set_xlabel(r'$\log_{10}\ D_{25}$ [ellipse-fit, arcmin]')
ax1.set_ylabel(r'$\log_{10}\ D_{25}$ [Hyperleda, arcmin]')

ax2.hexbin(plslga['PA'], plslga['PA_ORIG'], mincnt=1,
           cmap=plt.cm.Blues_r)
ax2.set_xlabel('PA [ellipse-fit, degrees]')
ax2.set_ylabel('PA [Hyperleda, degrees]')
ax2.plot([0, 180], [0, 180], color='k', ls='--', lw=3, alpha=0.7)

ax3.hexbin(plslga['BA'], plslga['BA_ORIG'], mincnt=1,
           cmap=plt.cm.Blues_r)
ax3.set_xlabel('b/a [ellipse-fit]')
ax3.set_ylabel('b/a [Hyperleda]')
ax3.plot([0, 1], [0, 1], color='k', ls='--', lw=3, alpha=0.7)

fig.subplots_adjust(hspace=0.3)

lslga[(lslga['LSLGA_ID'] > 0) * (lslga['D25'] == 0)]

ff = Table.read('/global/cfs/cdirs/desi/users/ioannis/LSLGA-data-DR9fg/LSLGA-reject-v6.0.fits')
ff

rr = Table.read('/global/cfs/cdirs/desi/users/ioannis/LSLGA-data-DR9fg/LSLGA-inspect-v6.0.fits')
rr

"""
Explanation: Geometry comparison of LSLGA galaxies
End of explanation
"""


plslga.colnames

def qa_radii(ax, xx, yy, xlim, ylim):
    ax.hexbin(np.log10(xx), np.log10(yy), mincnt=1, cmap=plt.cm.Blues_r)
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    #ax.plot([-1, 1], [-1, 1], color='k', ls='--', lw=3, alpha=0.7)
    #ax.set_xlabel(r'$\log_{10}\ D_{25}$ [ellipse-fit, arcmin]')
    #ax.set_ylabel(r'$\log_{10}\ D_{25}$ [Hyperleda, arcmin]')    

fig, ax = plt.subplots(2, 2, figsize=(12, 8), sharey=True, sharex=True)
for ax1, typ in zip(ax.flat, ['REX', 'DEV', 'EXP', 'SER']):
    ww = np.where(plslga['TYPE'] == typ)[0]
    xx = plslga['RADIUS_SB26'][ww] / 60 
    yy = plslga['RADIUS_SB26'][ww] / plslga['RADIUS_SB25'][ww]
    print(typ, np.median(yy))
    ax1.set_xlim(-1, 0.7)
    ax1.set_ylim(-0.1, 0.3)
    ax1.axhline(y=0, ls='-', lw=2, color='gray')
    ax1.axhline(y=np.log10(1.25), ls='--', lw=2, color='red')
    if typ == 'SER':
        ax1.hexbin(np.log10(xx), np.log10(yy), mincnt=1, cmap=plt.cm.Blues_r)
    else:
        ax1.scatter(np.log10(xx), np.log10(yy), s=5, color='cornflowerblue')
    ax1.text(0.9, 0.9, '{} (N={})'.format(typ, len(ww)), 
             transform=ax1.transAxes, ha='right', va='center',
             fontsize=14)
ax[0, 0].set_ylabel(r'$\log_{10}\ R(26) / R(25)$')
ax[1, 0].set_ylabel(r'$\log_{10}\ R(26) / R(25)$')
ax[1, 0].set_xlabel(r'$\log_{10}\ R(26)$ (arcmin)')
ax[1, 1].set_xlabel(r'$\log_{10}\ R(26)$ (arcmin)')
fig.subplots_adjust(hspace=0.1, wspace=0.1)
print('Median over all types: {:.4f}'.format(np.median(plslga['RADIUS_SB26'] / plslga['RADIUS_SB25'])))

fig, ax = plt.subplots(2, 2, figsize=(12, 8), sharey=True, sharex=True)
for ax1, typ in zip(ax.flat, ['REX', 'DEV', 'EXP', 'SER']):
    ww = np.where(plslga['TYPE'] == typ)[0]
    xx = plslga['RADIUS_SB26'][ww] / 60 
    yy = plslga['RADIUS_SB26'][ww] / plslga['RADIUS_SB24'][ww]
    print(typ, np.median(yy))
    ax1.set_xlim(-1, 0.7)
    ax1.set_ylim(-0.1, 0.6)
    ax1.axhline(y=0, ls='-', lw=2, color='gray')
    ax1.axhline(y=np.log10(1.25), ls='--', lw=2, color='red')
    if typ == 'SER':
        ax1.hexbin(np.log10(xx), np.log10(yy), mincnt=1, cmap=plt.cm.Blues_r)
    else:
        ax1.scatter(np.log10(xx), np.log10(yy), s=5, color='cornflowerblue')
    ax1.text(0.9, 0.9, '{} (N={})'.format(typ, len(ww)), 
             transform=ax1.transAxes, ha='right', va='center',
             fontsize=14)
ax[0, 0].set_ylabel(r'$\log_{10}\ R(26) / R(24)$')
ax[1, 0].set_ylabel(r'$\log_{10}\ R(26) / R(24)$')
ax[1, 0].set_xlabel(r'$\log_{10}\ R(26)$ (arcmin)')
ax[1, 1].set_xlabel(r'$\log_{10}\ R(26)$ (arcmin)')
fig.subplots_adjust(hspace=0.1, wspace=0.1)
print('Median over all types: {:.4f}'.format(np.median(plslga['RADIUS_SB26'] / plslga['RADIUS_SB24'])))

"""
Explanation: Compare the various sizes
End of explanation
"""


