#!/usr/bin/env python

# Initial setup...
import numpy as np
import pandas as pd
from astropy.io import fits
import fitsio
from scipy import interpolate
import glob
import math
import os
import shutil
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

#--------------------------------------------------------------------------
# Main code.

def main():

    import argparse

    """Create command line arguments"""
    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument('--planFile', help='name of the input plan file', default='Zeropoints_for_QA_plan.par')
    parser.add_argument('--verbose', help='verbosity level of output to screen (0,1,2,...)', default=0, type=int)
    args = parser.parse_args()

    if args.verbose > 0: print args

    status = zeropoints_for_qa(args)

    return status

#--------------------------------------------------------------------------

def zeropoints_for_qa(args):

    import paramFile

    # Read contents of planFile into python dictionary...
    planDict = paramFile.readParamFile(args.planFile, args.verbose)

    # Check that all the required keywords were found...
    requiredKeywordList = ['matchFileListFile', 
                           'qaDir', 
                           'blancoOpticsFile', 
                           'procEpochsFile',
                           'combEpochsFile']

    flag = 0
    for requiredKeyword in requiredKeywordList:
        if requiredKeyword not in planDict:
            print """Required keyword '%s' not in planFile %s""" % (requiredKeyword, args.planFile)
            flag = flag + 1
    if flag > 0:
        print 'Returning with error code 1 now...'
        return 1

    # Read in the matchFileListFile...
    matchFileListFile = planDict['matchFileListFile']
    if os.path.isfile(matchFileListFile)==False:
        print """matchFileListFile %s does not exist...""" % (matchFileListFile)
        print 'Returning with error code 1 now...'
        return 1
    # np.genfromtxt yielded some issues downstream when 
    #  matchFileListFile only contained a single entry, 
    #  so here we go back to the basics!...
    inputFileList = []
    with open(matchFileListFile,'r') as fin:
        for lin in fin:
            # Skip commented out lines...
            if lin[0] == '#': continue
            lins=lin.strip().split()
            inputFileList.append(lins[0])
    inputFileArray = np.array(inputFileList)

    # Create qaDir if it does not already exist...
    # See 
    # http://stackoverflow.com/questions/273192/how-to-check-if-a-directory-exists-and-create-it-if-necessary/14364249#14364249 
    # for avoiding the race condition beteween checking  a directory exists 
    # and creating the directory... 
    qaDir = planDict['qaDir']
    try: 
        os.makedirs(qaDir)
    except OSError:
        if not os.path.isdir(qaDir):
            raise

    # Read in the Blanco optics cleaning history file...
    blancoOpticsFile = planDict['blancoOpticsFile']
    if blancoOpticsFile.lower() == 'default':
        # Grab path and name of blanco_optics_cleaning_events.csv file 
        #  in the DECam_PGCM data directory...
        #  Is there a better way to do this?
        #  See also:
        #  http://stackoverflow.com/questions/779495/python-access-data-in-package-subdirectory
        # Absolute path for the directory containing this module:
        moduleAbsPathName = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
        blancoOpticsFile = os.path.join(moduleAbsPathName, "data", "blanco_optics_cleaning_events.csv")
    if os.path.isfile(blancoOpticsFile)==False:
        print """blancoOpticsFile %s does not exist...""" % (blancoOpticsFile)
        print 'Returning with error code 1 now...'
        return 1
    blancoOpticsDF = pd.read_csv(blancoOpticsFile, comment='#')
    if args.verbose > 1:
        print blancoOpticsDF


    # Read in the data processing epochs history file...
    procEpochsFile = planDict['procEpochsFile']
    if procEpochsFile.lower() == 'default':
        # Grab path and name of processing_epochs.csv file 
        #  in the DECam_PGCM data directory...
        #  Is there a better way to do this?
        #  See also:
        #  http://stackoverflow.com/questions/779495/python-access-data-in-package-subdirectory
        # Absolute path for the directory containing this module:
        moduleAbsPathName = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
        procEpochsFile = os.path.join(moduleAbsPathName, "data", "processing_epochs.csv")
    if os.path.isfile(procEpochsFile)==False:
        print """procEpochsFile %s does not exist...""" % (procEpochsFile)
        print 'Returning with error code 1 now...'
        return 1
    procEpochsDF = pd.read_csv(procEpochsFile, comment='#')
    if args.verbose > 1:
        print procEpochsDF

    # Read in the combined epochs history file...
    combEpochsFile = planDict['combEpochsFile']
    if combEpochsFile.lower() == 'default':
        # Grab path and name of combined_epochs.csv file 
        #  in the DECam_PGCM data directory...
        #  Is there a better way to do this?
        #  See also:
        #  http://stackoverflow.com/questions/779495/python-access-data-in-package-subdirectory
        # Absolute path for the directory containing this module:
        moduleAbsPathName = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
        combEpochsFile = os.path.join(moduleAbsPathName, "data", "combined_epochs.csv")
    if os.path.isfile(combEpochsFile)==False:
        print """combEpochsFile %s does not exist...""" % (combEpochsFile)
        print 'Returning with error code 1 now...'
        return 1
    combEpochsDF = pd.read_csv(combEpochsFile, comment='#')
    if args.verbose > 1:
        print combEpochsDF


    # Loop through the input files in inpuFileArray
    for inputFile in inputFileArray:

        print inputFile

        baseName = os.path.basename(inputFile)
        baseNameNoExt = os.path.splitext(baseName)[0]

        df = pd.read_csv(inputFile, usecols=['EXPNUM_2',
                                             'BAND_2',
                                             'MJD_OBS_2',
                                             'AIRMASS_2',
                                             'FLUX_PSF_2', 
                                             'EXPTIME_2', 
                                             'MAG_PSF_MEDIAN_1'])

        df.rename(columns={'EXPNUM_2':'expnum', 
                           'BAND_2':'band', 
                           'MJD_OBS_2':'mjd', 
                           'AIRMASS_2':'airmass', 
                           'FLUX_PSF_2':'flux_psf',  
                           'EXPTIME_2':'exptime', 
                           'MAG_PSF_MEDIAN_1':'mag_std'
                           }, inplace=True)
    
        # Check to ensure that there is only one type of filter b
and in this data file...
        bandsList = df.band.unique() 
        if bandsList.size > 1:
            print 
            print """inputFile %s appears to contain multiple bands (%s)...  Skipping""" % \
                (inputFile, bandsList)
            print 
            continue
        else:
            band = bandsList[0]
            print """inputFile %s contains data for band (%s)...""" % \
                (inputFile, band)

        # Add dmag [= mag_std - mag_inst] 
column
        df['dmag'] = df['mag_std'] + 2.5*np.log10(df['flux_psf']/df['exptime'])

        # Remove any entries with weird or extreme airmasses...
        mask = (df.airmass >=1.0) & (df.airmass < 3.0)
        df = df[mask]

        # Remove exposures with noticeably large variation in ccd-to-ccd ZPs...

        # Let's plot a 2D histogram of log(Nobs), binned by dmag and airmass...
        x=df.mjd
        y=df.dmag
        fig, axs = plt.subplots(ncols=1)
        ax=axs
        hb = ax.hexbin(x, y, gridsize=100, bins='log', cmap='inferno')
        ax.set_title(baseName)
        ax.set_xlabel("MJD")
        ax.set_ylabel("dmag")
        cb = fig.colorbar(hb, ax=ax)
        cb.set_label('log10(N)')
        plt.g