#!/usr/bin/env python
# vim: set fileencoding=utf-8 ts=4 sts=4 sw=4 et tw=80 :
#
# Compare an image file and its associated uncertainty image.
#
# <NAME>
# Created:       2021-06-03
# Last modified: 2021-06-03
#--------------------------------------------------------------------------
#**************************************************************************
#--------------------------------------------------------------------------

## Logging setup:
import logging
#logging.basicConfig(level=logging.DEBUG)
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)
#logger.setLevel(logging.DEBUG)
logger.setLevel(logging.INFO)

## Current version:
__version__ = "0.0.1"

## Python version-agnostic module reloading:
try:
    reload                              # Python 2.7
except NameError:
    try:
        from importlib import reload    # Python 3.4+
    except ImportError:
        from imp import reload          # Python 3.0 - 3.3

## Modules:
#import argparse
#import shutil
import resource
import signal
#import glob
import gc
import os
import sys
import time
#import vaex
#import calendar
#import ephem
import numpy as np
#from numpy.lib.recfunctions import append_fields
#import datetime as dt
#from dateutil import parser as dtp
#import scipy.linalg as sla
#import scipy.signal as ssig
#import scipy.ndimage as ndi
#import scipy.optimize as opti
#import scipy.interpolate as stp
#import scipy.spatial.distance as ssd
import matplotlib.pyplot as plt
#import matplotlib.cm as cm
#import matplotlib.ticker as mt
#import matplotlib._pylab_helpers as hlp
#from matplotlib.colors import LogNorm
#import matplotlib.colors as mplcolors
#import matplotlib.collections as mcoll
#import matplotlib.gridspec as gridspec
#from functools import partial
#from collections import OrderedDict
#from collections.abc import Iterable
#import multiprocessing as mp
#np.set_printoptions(suppress=True, linewidth=160)
#import pandas as pd
#import statsmodels.api as sm
#import statsmodels.formula.api as smf
#from statsmodels.regression.quantile_regression import QuantReg
#import PIL.Image as pli
#import seaborn as sns
#import cmocean
import theil_sen as ts
#import window_filter as wf
#import itertools as itt
_have_np_vers = float('.'.join(np.__version__.split('.')[:2]))

##--------------------------------------------------------------------------##
## Disable buffering on stdout/stderr:
class Unbuffered(object):
   def __init__(self, stream):
       self.stream = stream
   def write(self, data):
       self.stream.write(data)
       self.stream.flush()
   def __getattr__(self, attr):
       return getattr(self.stream, attr)

sys.stdout = Unbuffered(sys.stdout)
sys.stderr = Unbuffered(sys.stderr)

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

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

## Home-brew robust statistics:
try:
    import robust_stats
    reload(robust_stats)
    rs = robust_stats
except ImportError:
    logger.error("module robust_stats not found!  Install and retry.")
    sys.stderr.write("\nError!  robust_stats module not found!\n"
           "Please install and try again ...\n\n")
    sys.exit(1)

## Home-brew KDE:
#try:
#    import my_kde
#    reload(my_kde)
#    mk = my_kde
#except ImportError:
#    logger.error("module my_kde not found!  Install and retry.")
#    sys.stderr.write("\nError!  my_kde module not found!\n"
#           "Please install and try again ...\n\n")
#    sys.exit(1)

## Fast FITS I/O:
#try:
#    import fitsio
#except ImportError:
#    logger.error("fitsio module not found!  Install and retry.")
#    sys.stderr.write("\nError: fitsio module not found!\n")
#    sys.exit(1)

## Various from astropy:
try:
#    import astropy.io.ascii as aia
    import astropy.io.fits as pf
#    import astropy.io.votable as av
#    import astropy.table as apt
#    import astropy.time as astt
#    import astropy.wcs as awcs
#    from astropy import constants as aconst
#    from astropy import coordinates as coord
#    from astropy import units as uu
except ImportError:
#    logger.error("astropy module not found!  Install and retry.")
    sys.stderr.write("\nError: astropy module not found!\n")
    sys.exit(1)

## Star extraction:
#try:
#    import easy_sep
#    reload(easy_sep)
#except ImportError:
#    logger.error("easy_sep module not found!  Install and retry.")
#    sys.stderr.write("Error: easy_sep module not found!\n\n")
#    sys.exit(1)
#pse = easy_sep.EasySEP()

##--------------------------------------------------------------------------##
## Colors for fancy terminal output:
NRED    = '\033[0;31m'   ;  BRED    = '\033[1;31m'
NGREEN  = '\033[0;32m'   ;  BGREEN  = '\033[1;32m'
NYELLOW = '\033[0;33m'   ;  BYELLOW = '\033[1;33m'
NBLUE   = '\033[0;34m'   ;  BBLUE   = '\033[1;34m'
NMAG    = '\033[0;35m'   ;  BMAG    = '\033[1;35m'
NCYAN   = '\033[0;36m'   ;  BCYAN   = '\033[1;36m'
NWHITE  = '\033[0;37m'   ;  BWHITE  = '\033[1;37m'
ENDC    = '\033[0m'

## Suppress colors in cron jobs:
if (os.getenv('FUNCDEF') == '--nocolors'):
    NRED    = ''   ;  BRED    = ''
    NGREEN  = ''   ;  BGREEN  = ''
    NYELLOW = ''   ;  BYELLOW = ''
    NBLUE   = ''   ;  BBLUE   = ''
    NMAG    = ''   ;  BMAG    = ''
    NCYAN   = ''   ;  BCYAN   = ''
    NWHITE  = ''   ;  BWHITE  = ''
    ENDC    = ''

## Fancy text:
degree_sign = u'\N{DEGREE SIGN}'

## Dividers:
halfdiv = '-' * 40
fulldiv = '-' * 80

##--------------------------------------------------------------------------##
## Save FITS image with clobber (astropy / pyfits):
#def qsave(iname, idata, header=None, padkeys=1000, **kwargs):
#    this_func = sys._getframe().f_code.co_name
#    parent_func = sys._getframe(1).f_code.co_name
#    sys.stderr.write("Writing to '%s' ... " % iname)
#    if header:
#        while (len(header) < padkeys):
#            header.append() # pad header
#    if os.path.isfile(iname):
#        os.remove(iname)
#    pf.writeto(iname, idata, header=header, **kwargs)
#    sys.stderr.write("done.\n")

##--------------------------------------------------------------------------##
## Save FITS image with clobber (fitsio):
#def qsave(iname, idata, header=None, **kwargs):
#    this_func = sys._getframe().f_code.co_name
#    parent_func = sys._getframe(1).f_code.co_name
#    sys.stderr.write("Writing to '%s' ... " % iname)
#    #if os.path.isfile(iname):
#    #    os.remove(iname)
#    fitsio.write(iname, idata, clobber=True, header=header, **kwargs)
#    sys.stderr.write("done.\n")

##--------------------------------------------------------------------------##
def ldmap(things):
    return dict(zip(things, range(len(things))))

def argnear(vec, val):
    return (np.abs(vec - val)).argmin()


##--------------------------------------------------------------------------##
## New-style string formatting (more at https://pyformat.info/):

#oldway = '%s %s' % ('one', 'two')
#newway = '{} {}'.format('one', 'two')

#oldway = '%d %d' % (1, 2)
#newway = '{} {}'.format(1, 2)

# With padding:
#oldway = '%10s' % ('test',)        # right-justified
#newway = '{:>10}'.format('test')   # right-justified
#oldway = '%-10s' % ('test',)       #  left-justified
#newway = '{:10}'.format('test')    #  left-justified

# Ordinally:
#newway = '{1} {0}'.format('one', 'two')     # prints "two one"

# Dictionarily:
#newway = '{lastname}, {firstname}'.format(firstname='Rob', lastname='Siverd')

# Centered (new-only):
#newctr = '{:^10}'.format('test')      # prints "   test   "

# Numbers:
#oldway = '%06.2f' % (3.141592653589793,)
#newway = '{:06.2f}'.format(3.141592653589793)

##--------------------------------------------------------------------------##
## Quick ASCII I/O:
#data_file = 'data.txt'
#gftkw = {'encoding':None} if (_have_np_vers >= 1.14) else {}
#gftkw.update({'names':True, 'autostrip':True})
#gftkw.update({'delimiter':'|', 'comments':'%0%0%0%0'})
#gftkw.update({'loose':True, 'invalid_raise':False})
#all_data = np.genfromtxt(data_file, dtype=None, **gftkw)
#all_data = aia.read(data_file)

#all_data = pd.read_csv(data_file)
#all_data = pd.read_table(data_file, delim_whitespace=True)
#all_data = pd.read_table(data_file, skipinitialspace=True)
#all_data = pd.read_table(data_file, sep='|')
#fields = all_data.dtype.names
#if not fields:
#    x = all_data[:, 0]
#    y = all_data[:, 1]
#else:
#    x = all_data[fields[0]]
#    y = all_data[fields[1]]

#vot_file = 'neato.xml'
#vot_data = av.parse_single_table(vot_file)
#vot_data = av.parse_single_table(vot_file).to_table()

##--------------------------------------------------------------------------##
## Quick FITS I/O:
ifile = 'SPITZER_I2_44772864_0004_0000_2_cbcd.fits'
ufile = 'SPITZER_I2_44772864_0004_0000_2_cbunc.fits'
idata, ihdrs = pf.getdata(ifile, header=True)
udata, uhdrs = pf.getdata(ufile, header=True)

gain = ihdrs['GAIN']
exptime = ihdrs['EXPTIME']
fluxconv = ihdrs['FLUXCONV']

ignore = np.isnan(idata) | np.isnan(udata)

isafe = idata[~ignore]
usafe = udata[~ignore]

ignore = (isafe <= 0.0)
iclean = isafe[~ignore]
uclean = usafe[~ignore]

ui_ratio = uclean / iclean

## Try to reproduce the idata:udata relationship ...
icounts = iclean / fluxconv * exptime * gain    # in electrons
ucounts = uclean / fluxconv * exptime * gain    # in electrons

#icounts -= np.median(icounts)

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

## Estimate icounts:ucounts relationship from bright pixels:
cutoff = 1e3
bright = (icounts >= cutoff)
ic_fit = icounts[bright]
uc_fit = ucounts[bright]
vc_fit = uc_fit**2

sys.stderr.write("Fitting variance(counts) for bright pixels ... ")
model = ts.linefit(ic_fit, vc_fit)
sys.stderr.write("done.\n")
#model = np.array([375., 1.05])

## A line for plotting:
pcounts = np.linspace(0.1, 3e4, 1000)
pcounts = np.logspace(-1.0, 4.5, 1000)
pvarian = model[0] + model[1] * pcounts

##--------------------------------------------------------------------------##
## Theil-Sen line-fitting (linear):
#model = ts.linefit(xvals, yvals)
#icept, slope = ts.linefit(xvals, yvals)

## Theil-Sen line-fitting (loglog):
#xvals, yvals = np.log10(original_xvals), np.log10(original_yvals)
#xvals, yvals = np.log10(df['x'].values), np.log10(df['y'].values)
#llmodel = ts.linefit(np.log10(xvals), np.log10(yvals))
#icept, slope = ts.linefit(xvals, yvals)
#fit_exponent = slope
#fit_multiplier = 10**icept
#bestfit_x = np.arange(5000)
#bestfit_y = fit_multiplier * bestfit_x**fit_exponent

## Log-log evaluator:
#def loglog_eval(xvals, model):
#    icept, slope = model
#    return 10**icept * xvals**slope
#def loglog_eval(xvals, icept, slope):
#    return 10**icept * xvals**slope

##--------------------------------------------------------------------------##
## Plot config:

# gridspec examples:
# https://matplotlib.org/users/gridspec.html

#gs1 = gridspec.GridSpec(4, 4)
#gs1.update(wspace=0.025, hspace=0.05)  # set axis spacing

#ax1 = plt.subplot2grid((3, 3), (0, 0), colspan=3) # top-left + center + right
#ax2 = plt.subplot2grid((3, 3), (1, 0), colspan=2) # mid-left + mid-center
#ax3 = plt.subplot2grid((3, 3), (1, 2), rowspan=2) # mid-right + bot-right
#ax4 = plt.subplot2grid((3, 3), (2, 0))            # bot-left
#ax5 = plt.subplot2grid((3, 3), (2, 1))            # bot-center


##--------------------------------------------------------------------------##
#plt.style.use('bmh')   # Bayesian Methods for Hackers style
fig_dims = (12, 10)
fig = plt.figure(1, figsize=fig_dims)
plt.gcf().clf()
#fig, axs = plt.subplots(2, 2, sharex=True, figsize=fig_dims, num=1)
# sharex='col' | sharex='row'
#fig.frameon = False # disable figure frame drawing
#fig.subplots_adjust(left=0.07, right=0.95)
#ax1 = plt.subplot(gs[0, 0])
ax1 = fig.add_subplot(111)
#ax1 = fig.add_axes([0, 0, 1, 1])
#ax1.patch.set_facecolor((0.8, 0.8, 0.8))
#ax1.grid(True)
#ax1.axis('off')
ax1.grid(True)
#ax1.scatter(iclean, uclean, lw=0, s=5)
ax1.scatter(icounts, ucounts**2, lw=0, s=5)

ax1.plot(pcounts, pvarian, c='r')
ax1.set_yscale('log')
ax1.set_xscale('log')

plot_name = 'gain_log.png'
fig.tight_layout() # adjust boundaries sensibly, matplotlib v1.1+
plt.draw()
fig.savefig(plot_name, bbox_inches='tight')

ax1.set_xscale('linear')
ax1.set_yscale('linear')
plot_name = 'gain_lin.png'
fig.tight_layout() # adjust boundaries sensibly, matplotlib v1.1+
plt.draw()
fig.savefig(plot_name, bbox_inches='tight')

## Disable axis offsets:
#ax1.xaxis.get_major_formatter().set_useOffset(False)
#ax1.yaxis.get_major_formatter().set_useOffset(False)

#ax1.plot(kde_pnts, kde_vals)

#blurb = "some text"
#ax1.text(0.5, 0.5, blurb, transform=ax1.transAxes)
#ax1.text(0.5, 0.5, blurb, transform=ax1.transAxes,
#      va='top', ha='left', bbox=dict(facecolor='white', pad=10.0))
#      fontdict={'family':'monospace'}) # fixed-width

#colors = cm.rainbow(np.linspace(0, 1, len(plot_list)))
#for camid, c in zip(plot_list, colors):
#    cam_data = subsets[camid]
#    xvalue = cam_data['CCDATEMP']
#    yvalue = cam_data['PIX_MED']
#    yvalue = cam_data['IMEAN']
#    ax1.scatter(xvalue, yvalue, color=c, lw=0, label=camid)

#mtickpos = [2,5,7]
#ndecades = 1.0   # for symlog, set width of linear portion in units of dex
#nonposx='mask' | nonposx='clip' | nonposy='mask' | nonposy='clip'
#ax1.set_xscale('log', basex=10, nonposx='mask', subsx=mtickpos)
#ax1.set_xscale('log', nonposx='clip', subsx=[3])
#ax1.set_yscale('symlog', basey=10, linthreshy=0.1, linscaley=ndecades)
#ax1.xaxis.set_major_formatter(formatter) # re-format x ticks
#ax1.set_ylim(ax1.get_ylim()[::-1])
#ax1.set_xlabel('whatever', labelpad=30)  # push X label down 

#ax1.set_xticks([1.0, 3.0, 10.0, 30.0, 100.0])
#ax1.set_xticks([1, 2, 3], ['Jan', 'Feb', 'Mar'])
#for label in ax1.get_xticklabels():
#    label.set_rotation(30)
#    label.set_fontsize(14) 

#ax1.xaxis.label.set_fontsize(18)
#ax1.yaxis.label.set_fontsize(18)

#ax1.set_xlim(nice_limits(xvec, pctiles=[1,99], pad=1.2))
#ax1.set_ylim(nice_limits(yvec, pctiles=[1,99], pad=1.2))

#spts = ax1.scatter(x, y, lw=0, s=5)
##cbar = fig.colorbar(spts, orientation='vertical')   # old way
#cbnorm = mplcolors.Normalize(*spts.get_clim())
#scm = plt.cm.ScalarMappable(norm=cbnorm, cmap=spts.cmap)
#scm.set_array([])
#cbar = fig.colorbar(scm, orientation='vertical')
#cbar = fig.colorbar(scm, ticks=cs.levels, orientation='vertical') # contours
#cbar.formatter.set_useOffset(False)
#cbar.update_ticks()

fig.tight_layout() # adjust boundaries sensibly, matplotlib v1.1+
plt.draw()
#fig.savefig(plot_name, bbox_inches='tight')




######################################################################
# CHANGELOG (compare_images.py):
#---------------------------------------------------------------------
#
#  2021-06-03:
#     -- Increased __version__ to 0.0.1.
#     -- First created compare_images.py.
#

#!/usr/bin/env python
# vim: set fileencoding=utf-8 ts=4 sts=4 sw=4 et tw=80 :
#
# Compare an image file and its associated uncertainty image.
#
# <NAME>
# Created:       2021-06-03
# Last modified: 2021-06-03
#--------------------------------------------------------------------------
#**************************************************************************
#--------------------------------------------------------------------------

## Logging setup:
import logging
#logging.basicConfig(level=logging.DEBUG)
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)
#logger.setLevel(logging.DEBUG)
logger.setLevel(logging.INFO)

## Current version:
__version__ = "0.0.1"

## Python version-agnostic module reloading:
try:
    reload                              # Python 2.7
except NameError:
    try:
        from importlib import reload    # Python 3.4+
    except ImportError:
        from imp import reload          # Python 3.0 - 3.3

## Modules:
#import argparse
#import shutil
import resource
import signal
#import glob
import gc
import os
import sys
import time
#import vaex
#import calendar
#import ephem
import numpy as np
#from numpy.lib.recfunctions import append_fields
#import datetime as dt
#from dateutil import parser as dtp
#import scipy.linalg as sla
#import scipy.signal as ssig
#import scipy.ndimage as ndi
#import scipy.optimize as opti
#import scipy.interpolate as stp
#import scipy.spatial.distance as ssd
import matplotlib.pyplot as plt
#import matplotlib.cm as cm
#import matplotlib.ticker as mt
#import matplotlib._pylab_helpers as hlp
#from matplotlib.colors import LogNorm
#import matplotlib.colors as mplcolors
#import matplotlib.collections as mcoll
#import matplotlib.gridspec as gridspec
#from functools import partial
#from collections import OrderedDict
#from collections.abc import Iterable
#import multiprocessing as mp
#np.set_printoptions(suppress=True, linewidth=160)
#import pandas as pd
#import statsmodels.api as sm
#import statsmodels.formula.api as smf
#from statsmodels.regression.quantile_regression import QuantReg
#import PIL.Image as pli
#import seaborn as sns
#import cmocean
import theil_sen as ts
#import window_filter as wf
#import itertools as itt
_have_np_vers = float('.'.join(np.__version__.split('.')[:2]))

##--------------------------------------------------------------------------##
## Disable buffering on stdout/stderr:
class Unbuffered(object):
   def __init__(self, stream):
       self.stream = stream
   def write(self, data):
       self.stream.write(data)
       self.stream.flush()
   def __getattr__(self, attr):
       return getattr(self.stream, attr)

sys.stdout = Unbuffered(sys.stdout)
sys.stderr = Unbuffered(sys.stderr)

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

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

## Home-brew robust statistics:
try:
    import robust_stats
    reload(robust_stats)
    rs = robust_stats
except ImportError:
    logger.error("module robust_stats not found!  Install and retry.")
    sys.stderr.write("\nError!  robust_stats module not found!\n"
           "Please install and try again ...\n\n")
    sys.exit(1)

## Home-brew KDE:
#try:
#    import my_kde
#    reload(my_kde)
#    mk = my_kde
#except ImportError:
#    logger.error("module my_kde not found!  Install and retry.")
#    sys.stderr.write("\nError!  my_kde module not found!\n"
#           "Please install and try again ...\n\n")
#    sys.exit(1)

## Fast FITS I/O:
#try:
#    import fitsio
#except ImportError:
#    logger.error("fitsio module not found!  Install and retry.")
#    sys.stderr.write("\nError: fitsio module not found!\n")
#    sys.exit(1)

## Various from astropy:
try:
#    import astropy.io.ascii as aia
    import astropy.io.fits as pf
#    import astropy.io.votable as av
#    import astropy.table as apt
#    import astropy.time as astt
#    import astropy.wcs as awcs
#    from astropy import constants as aconst
#    from astropy import coordinates as coord
#    from astropy import units as uu
except ImportError:
#    logger.error("astropy module not found!  Install and retry.")
    sys.stderr.write("\nError: astropy module not found!\n")
    sys.exit(1)

## Star extraction:
#try:
#    import easy_sep
#    reload(easy_sep)
#except ImportError:
#    logger.error("easy_sep module not found!  Install and retry.")
#    sys.stderr.write("Error: easy_sep module not found!\n\n")
#    sys.exit(1)
#pse = easy_sep.EasySEP()

##--------------------------------------------------------------------------##
## Colors for fancy terminal output:
NRED    = '\033[0;31m'   ;  BRED    = '\033[1;31m'
NGREEN  = '\033[0;32m'   ;  BGREEN  = '\033[1;32m'
NYELLOW = '\033[0;33m'   ;  BYELLOW = '\033[1;33m'
NBLUE   = '\033[0;34m'   ;  BBLUE   = '\033[1;34m'
NMAG    = '\033[0;35m'   ;  BMAG    = '\033[1;35m'
NCYAN   = '\033[0;36m'   ;  BCYAN   = '\033[1;36m'
NWHITE  = '\033[0;37m'   ;  BWHITE  = '\033[1;37m'
ENDC    = '\033[0m'

## Suppress colors in cron jobs:
if (os.getenv('FUNCDEF') == '--nocolors'):
    NRED    = ''   ;  BRED    = ''
    NGREEN  = ''   ;  BGREEN  = ''
    NYELLOW = ''   ;  BYELLOW = ''
    NBLUE   = ''   ;  BBLUE   = ''
    NMAG    = ''   ;  BMAG    = ''
    NCYAN   = ''   ;  BCYAN   = ''
    NWHITE  = ''   ;  BWHITE  = ''
    ENDC    = ''

## Fancy text:
degree_sign = u'\N{DEGREE SIGN}'

## Dividers:
halfdiv = '-' * 40
fulldiv = '-' * 80

##--------------------------------------------------------------------------##
## Save FITS image with clobber (astropy / pyfits):
#def qsave(iname, idata, header=None, padkeys=1000, **kwargs):
#    this_func = sys._getframe().f_code.co_name
#    parent_func = sys._getframe(1).f_code.co_name
#    sys.stderr.write("Writing to '%s' ... " % iname)
#    if header:
#        while (len(header) < padkeys):
#            header.append() # pad header
#    if os.path.isfile(iname):
#        os.remove(iname)
#    pf.writeto(iname, idata, header=header, **kwargs)
#    sys.stderr.write("done.\n")

##--------------------------------------------------------------------------##
## Save FITS image with clobber (fitsio):
#def qsave(iname, idata, header=None, **kwargs):
#    this_func = sys._getframe().f_code.co_name
#    parent_func = sys._getframe(1).f_code.co_name
#    sys.stderr.write("Writing to '%s' ... " % iname)
#    #if os.path.isfile(iname):
#    #    os.remove(iname)
#    fitsio.write(iname, idata, clobber=True, header=header, **kwargs)
#    sys.stderr.write("done.\n")

##--------------------------------------------------------------------------##
def ldmap(things):
    return dict(zip(things, range(len(things))))

def argnear(vec, val):
    return (np.abs(vec - val)).argmin()


##--------------------------------------------------------------------------##
## New-style string formatting (more at https://pyformat.info/):

#oldway = '%s %s' % ('one', 'two')
#newway = '{} {}'.format('one', 'two')

#oldway = '%d %d' % (1, 2)
#newway = '{} {}'.format(1, 2)

# With padding:
#oldway = '%10s' % ('test',)        # right-justified
#newway = '{:>10}'.format('test')   # right-justified
#oldway = '%-10s' % ('test',)       #  left-justified
#newway = '{:10}'.format('test')    #  left-justified

# Ordinally:
#newway = '{1} {0}'.format('one', 'two')     # prints "two one"

# Dictionarily:
#newway = '{lastname}, {firstname}'.format(firstname='Rob', lastname='Siverd')

# Centered (new-only):
#newctr = '{:^10}'.format('test')      # prints "   test   "

# Numbers:
#oldway = '%06.2f' % (3.141592653589793,)
#newway = '{:06.2f}'.format(3.141592653589793)

##--------------------------------------------------------------------------##
## Quick ASCII I/O:
#data_file = 'data.txt'
#gftkw = {'encoding':None} if (_have_np_vers >= 1.14) else {}
#gftkw.update({'names':True, 'autostrip':True})
#gftkw.update({'delimiter':'|', 'comments':'%0%0%0%0'})
#gftkw.update({'loose':True, 'invalid_raise':False})
#all_data = np.genfromtxt(data_file, dtype=None, **gftkw)
#all_data = aia.read(data_file)

#all_data = pd.read_csv(data_file)
#all_data = pd.read_table(data_file, delim_whitespace=True)
#all_data = pd.read_table(data_file, skipinitialspace=True)
#all_data = pd.read_table(data_file, sep='|')
#fields = all_data.dtype.names
#if not fields:
#    x = all_data[:, 0]
#    y = all_data[:, 1]
#else:
#    x = all_data[fields[0]]
#    y = all_data[fields[1]]

#vot_file = 'neato.xml'
#vot_data = av.parse_single_table(vot_file)
#vot_data = av.parse_single_table(vot_file).to_table()

##--------------------------------------------------------------------------##
## Quick FITS I/O:
ifile = 'SPITZER_I2_44772864_0004_0000_2_cbcd.fits'
ufile = 'SPITZER_I2_44772864_0004_0000_2_cbunc.fits'
idata, ihdrs = pf.getdata(ifile, header=True)
udata, uhdrs = pf.getdata(ufile, header=True)

gain = ihdrs['GAIN']
exptime = ihdrs['EXPTIME']
fluxconv = ihdrs['FLUXCONV']

ignore = np.isnan(idata) | np.isnan(udata)

isafe = idata[~ignore]
usafe = udata[~ignore]

ignore = (isafe <= 0.0)
iclean = isafe[~ignore]
uclean = usafe[~ignore]

ui_ratio = uclean / iclean

## Try to reproduce the idata:udata relationship ...
icounts = iclean / fluxconv * exptime * gain    # in electrons
ucounts = uclean / fluxconv * exptime * gain    # in electrons

#icounts -= np.median(icounts)

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

## Estimate icounts:ucounts relationship from bright pixels:
cutoff = 1e3
bright = (icounts >= cutoff)
ic_fit = icounts[bright]
uc_fit = ucounts[bright]
vc_fit = uc_fit**2

sys.stderr.write("Fitting variance(counts) for bright pixels ... ")
model = ts.linefit(ic_fit, vc_fit)
sys.stderr.write("done.\n")
#model = np.array([375., 1.05])

## A line for plotting:
pcounts = np.linspace(0.1, 3e4, 1000)
pcounts = np.logspace(-1.0, 4.5, 1000)
pvarian = model[0] + model[1] * pcounts

##--------------------------------------------------------------------------##
## Theil-Sen line-fitting (linear):
#model = ts.linefit(xvals, yvals)
#icept, slope = ts.linefit(xvals, yvals)

## Theil-Sen line-fitting (loglog):
#xvals, yvals = np.log10(original_xvals), np.log10(original_yvals)
#xvals, yvals = np.log10(df['x'].values), np.log10(df['y'].values)
#llmodel = ts.linefit(np.log10(xvals), np.log10(yvals))
#icept, slope = ts.linefit(xvals, yvals)
#fit_exponent = slope
#fit_multiplier = 10**icept
#bestfit_x = np.arange(5000)
#bestfit_y = fit_multiplier * bestfit_x**fit_exponent

## Log-log evaluator:
#def loglog_eval(xvals, model):
#    icept, slope = model
#    return 10**icept * xvals**slope
#def loglog_eval(xvals, icept, slope):
#    return 10**icept * xvals**slope

##--------------------------------------------------------------------------##
## Plot config:

# gridspec examples:
# https://matplotlib.org/users/gridspec.html

#gs1 = gridspec.GridSpec(4, 4)
#gs1.update(wspace=0.025, hspace=0.05)  # set axis spacing

#ax1 = plt.subplot2grid((3, 3), (0, 0), colspan=3) # top-left + center + right
#ax2 = plt.subplot2grid((3, 3), (1, 0), colspan=2) # mid-left + mid-center
#ax3 = plt.subplot2grid((3, 3), (1, 2), rowspan=2) # mid-right + bot-right
#ax4 = plt.subplot2grid((3, 3), (2, 0))            # bot-left
#ax5 = plt.subplot2grid((3, 3), (2, 1))            # bot-center


##--------------------------------------------------------------------------##
#plt.style.use('bmh')   # Bayesian Methods for Hackers style
fig_dims = (12, 10)
fig = plt.figure(1, figsize=fig_dims)
plt.gcf().clf()
#fig, axs = plt.subplots(2, 2, sharex=True, figsize=fig_dims, num=1)
# sharex='col' | sharex='row'
#fig.frameon = False # disable figure frame drawing
#fig.subplots_adjust(left=0.07, right=0.95)
#ax1 = plt.subplot(gs[0, 0])
ax1 = fig.add_subplot(111)
#ax1 = fig.add_axes([0, 0, 1, 1])
#ax1.patch.set_facecolor((0.8, 0.8, 0.8))
#ax1.grid(True)
#ax1.axis('off')
ax1.grid(True)
#ax1.scatter(iclean, uclean, lw=0, s=5)
ax1.scatter(icounts, ucounts**2, lw=0, s=5)

ax1.plot(pcounts, pvarian, c='r')
ax1.set_yscale('log')
ax1.set_xscale('log')

plot_name = 'gain_log.png'
fig.tight_layout() # adjust boundaries sensibly, matplotlib v1.1+
plt.draw()
fig.savefig(plot_name, bbox_inches='tight')

ax1.set_xscale('linear')
ax1.set_yscale('linear')
plot_name = 'gain_lin.png'
fig.tight_layout() # adjust boundaries sensibly, matplotlib v1.1+
plt.draw()
fig.savefig(plot_name, bbox_inches='tight')

## Disable axis offsets:
#ax1.xaxis.get_major_formatter().set_useOffset(False)
#ax1.yaxis.get_major_formatter().set_useOffset(False)

#ax1.plot(kde_pnts, kde_vals)

#blurb = "some text"
#ax1.text(0.5, 0.5, blurb, transform=ax1.transAxes)
#ax1.text(0.5, 0.5, blurb, transform=ax1.transAxes,
#      va='top', ha='left', bbox=dict(facecolor='white', pad=10.0))
#      fontdict={'family':'monospace'}) # fixed-width

#colors = cm.rainbow(np.linspace(0, 1, len(plot_list)))
#for camid, c in zip(plot_list, colors):
#    cam_data = subsets[camid]
#    xvalue = cam_data['CCDATEMP']
#    yvalue = cam_data['PIX_MED']
#    yvalue = cam_data['IMEAN']
#    ax1.scatter(xvalue, yvalue, color=c, lw=0, label=camid)

#mtickpos = [2,5,7]
#ndecades = 1.0   # for symlog, set width of linear portion in units of dex
#nonposx='mask' | nonposx='clip' | nonposy='mask' | nonposy='clip'
#ax1.set_xscale('log', basex=10, nonposx='mask', subsx=mtickpos)
#ax1.set_xscale('log', nonposx='clip', subsx=[3])
#ax1.set_yscale('symlog', basey=10, linthreshy=0.1, linscaley=ndecades)
#ax1.xaxis.set_major_formatter(formatter) # re-format x ticks
#ax1.set_ylim(ax1.get_ylim()[::-1])
#ax1.set_xlabel('whatever', labelpad=30)  # push X label down 

#ax1.set_xticks([1.0, 3.0, 10.0, 30.0, 100.0])
#ax1.set_xticks([1, 2, 3], ['Jan', 'Feb', 'Mar'])
#for label in ax1.get_xticklabels():
#    label.set_rotation(30)
#    label.set_fontsize(14) 

#ax1.xaxis.label.set_fontsize(18)
#ax1.yaxis.label.set_fontsize(18)

#ax1.set_xlim(nice_limits(xvec, pctiles=[1,99], pad=1.2))
#ax1.set_ylim(nice_limits(yvec, pctiles=[1,99], pad=1.2))

#spts = ax1.scatter(x, y, lw=0, s=5)
##cbar = fig.colorbar(spts, orientation='vertical')   # old way
#cbnorm = mplcolors.Normalize(*spts.get_clim())
#scm = plt.cm.ScalarMappable(norm=cbnorm, cmap=spts.cmap)
#scm.set_array([])
#cbar = fig.colorbar(scm, orientation='vertical')
#cbar = fig.colorbar(scm, ticks=cs.levels, orientation='vertical') # contours
#cbar.formatter.set_useOffset(False)
#cbar.update_ticks()

fig.tight_layout() # adjust boundaries sensibly, matplotlib v1.1+
plt.draw()
#fig.savefig(plot_name, bbox_inches='tight')




######################################################################
# CHANGELOG (compare_images.py):
#---------------------------------------------------------------------
#
#  2021-06-03:
#     -- Increased __version__ to 0.0.1.
#     -- First created compare_images.py.
#

#!/usr/bin/env python # vim: set fileencoding=utf-8 ts=4 sts=4 sw=4 et tw=80 : # # Compare an image file and its associated uncertainty image. # # <NAME> # Created:       2021-06-03 # Last modified: 2021-06-03 #-------------------------------------------------------------------------- #************************************************************************** #-------------------------------------------------------------------------- ## Logging setup: #logging.basicConfig(level=logging.DEBUG) #logger.setLevel(logging.DEBUG) ## Current version: ## Python version-agnostic module reloading: # Python 2.7 # Python 3.4+ # Python 3.0 - 3.3 ## Modules: #import argparse #import shutil #import glob #import vaex #import calendar #import ephem #from numpy.lib.recfunctions import append_fields #import datetime as dt #from dateutil import parser as dtp #import scipy.linalg as sla #import scipy.signal as ssig #import scipy.ndimage as ndi #import scipy.optimize as opti #import scipy.interpolate as stp #import scipy.spatial.distance as ssd #import matplotlib.cm as cm #import matplotlib.ticker as mt #import matplotlib._pylab_helpers as hlp #from matplotlib.colors import LogNorm #import matplotlib.colors as mplcolors #import matplotlib.collections as mcoll #import matplotlib.gridspec as gridspec #from functools import partial #from collections import OrderedDict #from collections.abc import Iterable #import multiprocessing as mp #np.set_printoptions(suppress=True, linewidth=160) #import pandas as pd #import statsmodels.api as sm #import statsmodels.formula.api as smf #from statsmodels.regression.quantile_regression import QuantReg #import PIL.Image as pli #import seaborn as sns #import cmocean #import window_filter as wf #import itertools as itt ##--------------------------------------------------------------------------## ## Disable buffering on stdout/stderr: ##--------------------------------------------------------------------------## ##--------------------------------------------------------------------------## ## Home-brew robust statistics: ## Home-brew KDE: #try: #    import my_kde #    reload(my_kde) #    mk = my_kde #except ImportError: #    logger.error("module my_kde not found!  Install and retry.") #    sys.stderr.write("\nError!  my_kde module not found!\n" #           "Please install and try again ...\n\n") #    sys.exit(1) ## Fast FITS I/O: #try: #    import fitsio #except ImportError: #    logger.error("fitsio module not found!  Install and retry.") #    sys.stderr.write("\nError: fitsio module not found!\n") #    sys.exit(1) ## Various from astropy: #    import astropy.io.ascii as aia #    import astropy.io.votable as av #    import astropy.table as apt #    import astropy.time as astt #    import astropy.wcs as awcs #    from astropy import constants as aconst #    from astropy import coordinates as coord #    from astropy import units as uu #    logger.error("astropy module not found!  Install and retry.") ## Star extraction: #try: #    import easy_sep #    reload(easy_sep) #except ImportError: #    logger.error("easy_sep module not found!  Install and retry.") #    sys.stderr.write("Error: easy_sep module not found!\n\n") #    sys.exit(1) #pse = easy_sep.EasySEP() ##--------------------------------------------------------------------------## ## Colors for fancy terminal output: ## Suppress colors in cron jobs: ## Fancy text: ## Dividers: ##--------------------------------------------------------------------------## ## Save FITS image with clobber (astropy / pyfits): #def qsave(iname, idata, header=None, padkeys=1000, **kwargs): #    this_func = sys._getframe().f_code.co_name #    parent_func = sys._getframe(1).f_code.co_name #    sys.stderr.write("Writing to '%s' ... " % iname) #    if header: #        while (len(header) < padkeys): #            header.append() # pad header #    if os.path.isfile(iname): #        os.remove(iname) #    pf.writeto(iname, idata, header=header, **kwargs) #    sys.stderr.write("done.\n") ##--------------------------------------------------------------------------## ## Save FITS image with clobber (fitsio): #def qsave(iname, idata, header=None, **kwargs): #    this_func = sys._getframe().f_code.co_name #    parent_func = sys._getframe(1).f_code.co_name #    sys.stderr.write("Writing to '%s' ... " % iname) #    #if os.path.isfile(iname): #    #    os.remove(iname) #    fitsio.write(iname, idata, clobber=True, header=header, **kwargs) #    sys.stderr.write("done.\n") ##--------------------------------------------------------------------------## ##--------------------------------------------------------------------------## ## New-style string formatting (more at https://pyformat.info/): #oldway = '%s %s' % ('one', 'two') #newway = '{} {}'.format('one', 'two') #oldway = '%d %d' % (1, 2) #newway = '{} {}'.format(1, 2) # With padding: #oldway = '%10s' % ('test',)        # right-justified #newway = '{:>10}'.format('test')   # right-justified #oldway = '%-10s' % ('test',)       #  left-justified #newway = '{:10}'.format('test')    #  left-justified # Ordinally: #newway = '{1} {0}'.format('one', 'two')     # prints "two one" # Dictionarily: #newway = '{lastname}, {firstname}'.format(firstname='Rob', lastname='Siverd') # Centered (new-only): #newctr = '{:^10}'.format('test')      # prints "   test   " # Numbers: #oldway = '%06.2f' % (3.141592653589793,) #newway = '{:06.2f}'.format(3.141592653589793) ##--------------------------------------------------------------------------## ## Quick ASCII I/O: #data_file = 'data.txt' #gftkw = {'encoding':None} if (_have_np_vers >= 1.14) else {} #gftkw.update({'names':True, 'autostrip':True}) #gftkw.update({'delimiter':'|', 'comments':'%0%0%0%0'}) #gftkw.update({'loose':True, 'invalid_raise':False}) #all_data = np.genfromtxt(data_file, dtype=None, **gftkw) #all_data = aia.read(data_file) #all_data = pd.read_csv(data_file) #all_data = pd.read_table(data_file, delim_whitespace=True) #all_data = pd.read_table(data_file, skipinitialspace=True) #all_data = pd.read_table(data_file, sep='|') #fields = all_data.dtype.names #if not fields: #    x = all_data[:, 0] #    y = all_data[:, 1] #else: #    x = all_data[fields[0]] #    y = all_data[fields[1]] #vot_file = 'neato.xml' #vot_data = av.parse_single_table(vot_file) #vot_data = av.parse_single_table(vot_file).to_table() ##--------------------------------------------------------------------------## ## Quick FITS I/O: ## Try to reproduce the idata:udata relationship ... # in electrons # in electrons #icounts -= np.median(icounts) ##--------------------------------------------------------------------------## ##--------------------------------------------------------------------------## ## Estimate icounts:ucounts relationship from bright pixels: #model = np.array([375., 1.05]) ## A line for plotting: ##--------------------------------------------------------------------------## ## Theil-Sen line-fitting (linear): #model = ts.linefit(xvals, yvals) #icept, slope = ts.linefit(xvals, yvals) ## Theil-Sen line-fitting (loglog): #xvals, yvals = np.log10(original_xvals), np.log10(original_yvals) #xvals, yvals = np.log10(df['x'].values), np.log10(df['y'].values) #llmodel = ts.linefit(np.log10(xvals), np.log10(yvals)) #icept, slope = ts.linefit(xvals, yvals) #fit_exponent = slope #fit_multiplier = 10**icept #bestfit_x = np.arange(5000) #bestfit_y = fit_multiplier * bestfit_x**fit_exponent ## Log-log evaluator: #def loglog_eval(xvals, model): #    icept, slope = model #    return 10**icept * xvals**slope #def loglog_eval(xvals, icept, slope): #    return 10**icept * xvals**slope ##--------------------------------------------------------------------------## ## Plot config: # gridspec examples: # https://matplotlib.org/users/gridspec.html #gs1 = gridspec.GridSpec(4, 4) #gs1.update(wspace=0.025, hspace=0.05)  # set axis spacing #ax1 = plt.subplot2grid((3, 3), (0, 0), colspan=3) # top-left + center + right #ax2 = plt.subplot2grid((3, 3), (1, 0), colspan=2) # mid-left + mid-center #ax3 = plt.subplot2grid((3, 3), (1, 2), rowspan=2) # mid-right + bot-right #ax4 = plt.subplot2grid((3, 3), (2, 0))            # bot-left #ax5 = plt.subplot2grid((3, 3), (2, 1))            # bot-center ##--------------------------------------------------------------------------## #plt.style.use('bmh')   # Bayesian Methods for Hackers style #fig, axs = plt.subplots(2, 2, sharex=True, figsize=fig_dims, num=1) # sharex='col' | sharex='row' #fig.frameon = False # disable figure frame drawing #fig.subplots_adjust(left=0.07, right=0.95) #ax1 = plt.subplot(gs[0, 0]) #ax1 = fig.add_axes([0, 0, 1, 1]) #ax1.patch.set_facecolor((0.8, 0.8, 0.8)) #ax1.grid(True) #ax1.axis('off') #ax1.scatter(iclean, uclean, lw=0, s=5) # adjust boundaries sensibly, matplotlib v1.1+ # adjust boundaries sensibly, matplotlib v1.1+ ## Disable axis offsets: #ax1.xaxis.get_major_formatter().set_useOffset(False) #ax1.yaxis.get_major_formatter().set_useOffset(False) #ax1.plot(kde_pnts, kde_vals) #blurb = "some text" #ax1.text(0.5, 0.5, blurb, transform=ax1.transAxes) #ax1.text(0.5, 0.5, blurb, transform=ax1.transAxes, #      va='top', ha='left', bbox=dict(facecolor='white', pad=10.0)) #      fontdict={'family':'monospace'}) # fixed-width #colors = cm.rainbow(np.linspace(0, 1, len(plot_list))) #for camid, c in zip(plot_list, colors): #    cam_data = subsets[camid] #    xvalue = cam_data['CCDATEMP'] #    yvalue = cam_data['PIX_MED'] #    yvalue = cam_data['IMEAN'] #    ax1.scatter(xvalue, yvalue, color=c, lw=0, label=camid) #mtickpos = [2,5,7] #ndecades = 1.0   # for symlog, set width of linear portion in units of dex #nonposx='mask' | nonposx='clip' | nonposy='mask' | nonposy='clip' #ax1.set_xscale('log', basex=10, nonposx='mask', subsx=mtickpos) #ax1.set_xscale('log', nonposx='clip', subsx=[3]) #ax1.set_yscale('symlog', basey=10, linthreshy=0.1, linscaley=ndecades) #ax1.xaxis.set_major_formatter(formatter) # re-format x ticks #ax1.set_ylim(ax1.get_ylim()[::-1]) #ax1.set_xlabel('whatever', labelpad=30)  # push X label down #ax1.set_xticks([1.0, 3.0, 10.0, 30.0, 100.0]) #ax1.set_xticks([1, 2, 3], ['Jan', 'Feb', 'Mar']) #for label in ax1.get_xticklabels(): #    label.set_rotation(30) #    label.set_fontsize(14) #ax1.xaxis.label.set_fontsize(18) #ax1.yaxis.label.set_fontsize(18) #ax1.set_xlim(nice_limits(xvec, pctiles=[1,99], pad=1.2)) #ax1.set_ylim(nice_limits(yvec, pctiles=[1,99], pad=1.2)) #spts = ax1.scatter(x, y, lw=0, s=5) ##cbar = fig.colorbar(spts, orientation='vertical')   # old way #cbnorm = mplcolors.Normalize(*spts.get_clim()) #scm = plt.cm.ScalarMappable(norm=cbnorm, cmap=spts.cmap) #scm.set_array([]) #cbar = fig.colorbar(scm, orientation='vertical') #cbar = fig.colorbar(scm, ticks=cs.levels, orientation='vertical') # contours #cbar.formatter.set_useOffset(False) #cbar.update_ticks() # adjust boundaries sensibly, matplotlib v1.1+ #fig.savefig(plot_name, bbox_inches='tight') ###################################################################### # CHANGELOG (compare_images.py): #--------------------------------------------------------------------- # #  2021-06-03: #     -- Increased __version__ to 0.0.1. #     -- First created compare_images.py. #