Circular imports

arclines/holy/fitting.py

@@ -6,9 +6,12 @@
 import numpy as np
 import pdb
 
+from arclines.utils import calc_fit_rms, func_val, robust_polyfit
+from arclines.holy.qa import arc_fit_qa
+
 
 def iterative_fitting(spec, tcent, ifit, IDs, llist, disp, plot_fil=None,
-                      verbose=False, load_pypit=False, aparm=None):
+                      verbose=False, aparm=None):
 
     if aparm is None:
         aparm = dict(llist='',
@@ -20,15 +23,6 @@ def iterative_fitting(spec, tcent, ifit, IDs, llist, disp, plot_fil=None,
                     n_final=4,           # Order of polynomial for final fit
                     nsig_rej=2.,         # Number of sigma for rejection
                     nsig_rej_final=3.0)  # Number of sigma for rejection (final fit)
-    # PYPIT
-    if load_pypit:
-        from pypit import pyputils
-        msgs = pyputils.get_dummy_logger()
-    from pypit import arutils
-    from pypit import arqa
-    if load_pypit:
-        arutils.dummy_settings()
-
     npix = spec.size
 
     # Setup for fitting
@@ -44,9 +38,9 @@ def iterative_fitting(spec, tcent, ifit, IDs, llist, disp, plot_fil=None,
     while (n_order <= aparm['n_final']) and (flg_quit is False):
         # Fit with rejection
         xfit, yfit = tcent[ifit], all_ids[ifit]
-        mask, fit = arutils.robust_polyfit(xfit, yfit, n_order, function=aparm['func'], sigma=aparm['nsig_rej'], minv=fmin, maxv=fmax)
+        mask, fit = robust_polyfit(xfit, yfit, n_order, function=aparm['func'], sigma=aparm['nsig_rej'], minv=fmin, maxv=fmax)
 
-        rms_ang = arutils.calc_fit_rms(xfit[mask==0], yfit[mask==0],
+        rms_ang = calc_fit_rms(xfit[mask==0], yfit[mask==0],
                                        fit, aparm['func'], minv=fmin, maxv=fmax)
         rms_pix = rms_ang/disp
         if verbose:
@@ -55,7 +49,7 @@ def iterative_fitting(spec, tcent, ifit, IDs, llist, disp, plot_fil=None,
         # Reject but keep originals (until final fit)
         ifit = list(ifit[mask == 0]) + sv_ifit
         # Find new points (should we allow removal of the originals?)
-        twave = arutils.func_val(fit, tcent, aparm['func'], minv=fmin, maxv=fmax)
+        twave = func_val(fit, tcent, aparm['func'], minv=fmin, maxv=fmax)
         for ss,iwave in enumerate(twave):
             mn = np.min(np.abs(iwave-llist['wave']))
             if mn/aparm['disp'] < aparm['match_toler']:
@@ -78,14 +72,14 @@ def iterative_fitting(spec, tcent, ifit, IDs, llist, disp, plot_fil=None,
     # Final fit (originals can now be rejected)
     fmin, fmax = 0., 1.
     xfit, yfit = tcent[ifit]/(npix-1), all_ids[ifit]
-    mask, fit = arutils.robust_polyfit(xfit, yfit, n_order, function=aparm['func'], sigma=aparm['nsig_rej_final'], minv=fmin, maxv=fmax)#, debug=True)
+    mask, fit = robust_polyfit(xfit, yfit, n_order, function=aparm['func'], sigma=aparm['nsig_rej_final'], minv=fmin, maxv=fmax)#, debug=True)
     irej = np.where(mask==1)[0]
     if len(irej) > 0:
         xrej = xfit[irej]
         yrej = yfit[irej]
         if verbose:
             for kk,imask in enumerate(irej):
-                wave = arutils.func_val(fit, xrej[kk], aparm['func'], minv=fmin, maxv=fmax)
+                wave = func_val(fit, xrej[kk], aparm['func'], minv=fmin, maxv=fmax)
                 print('Rejecting arc line {:g}; {:g}'.format(yfit[imask], wave))
     else:
         xrej = []
@@ -94,7 +88,7 @@ def iterative_fitting(spec, tcent, ifit, IDs, llist, disp, plot_fil=None,
     yfit = yfit[mask==0]
     ions = all_idsion[ifit][mask==0]
     # Final RMS
-    rms_ang = arutils.calc_fit_rms(xfit, yfit, fit, aparm['func'],
+    rms_ang = calc_fit_rms(xfit, yfit, fit, aparm['func'],
                                    minv=fmin, maxv=fmax)
     rms_pix = rms_ang/disp
     #
@@ -125,8 +119,9 @@ def iterative_fitting(spec, tcent, ifit, IDs, llist, disp, plot_fil=None,
         shift=0., tcent=tcent, rms=rms_pix)
     # QA
     if plot_fil is not None:
-        arqa.arc_fit_qa(None, final_fit, outfile=plot_fil)
+        arc_fit_qa(None, final_fit, plot_fil)
     # Return
     return final_fit
 
 
+

arclines/holy/grail.py

@@ -10,7 +10,6 @@
 from arclines.holy import patterns as arch_patt
 from arclines.holy import fitting as arch_fit
 from arclines.holy import utils as arch_utils
-from arclines.pypit_utils import find_peaks
 
 
 def basic(spec, lines, wv_cen, disp, siglev=20., min_ampl=300.,

arclines/holy/qa.py

@@ -0,0 +1,113 @@
+""" QA for arclines.holy
+"""
+from __future__ import (print_function, absolute_import, division, unicode_literals)
+
+import numpy as np
+
+from matplotlib import pyplot as plt
+from matplotlib import gridspec
+
+from arclines.utils import func_val
+
+def arc_fit_qa(slf, fit, outfile, ids_only=False, title=None):
+    """
+    QA for Arc spectrum
+
+    Parameters
+    ----------
+    fit : Wavelength fit
+    arc_spec : ndarray
+      Arc spectrum
+    outfile : str, optional
+      Name of output file
+    """
+
+    plt.rcdefaults()
+    plt.rcParams['font.family']= 'times new roman'
+
+    arc_spec = fit['spec']
+
+    # Begin
+    if not ids_only:
+        plt.figure(figsize=(8, 4.0))
+        plt.clf()
+        gs = gridspec.GridSpec(2, 2)
+        idfont = 'xx-small'
+    else:
+        plt.figure(figsize=(11, 8.5))
+        plt.clf()
+        gs = gridspec.GridSpec(1, 1)
+        idfont = 'small'
+
+    # Simple spectrum plot
+    ax_spec = plt.subplot(gs[:,0])
+    ax_spec.plot(np.arange(len(arc_spec)), arc_spec)
+    ymin, ymax = 0., np.max(arc_spec)
+    ysep = ymax*0.03
+    for kk, x in enumerate(fit['xfit']*fit['xnorm']):
+        yline = np.max(arc_spec[int(x)-2:int(x)+2])
+        # Tick mark
+        ax_spec.plot([x,x], [yline+ysep*0.25, yline+ysep], 'g-')
+        # label
+        ax_spec.text(x, yline+ysep*1.3,
+                     '{:s} {:g}'.format(fit['ions'][kk], fit['yfit'][kk]), ha='center', va='bottom',
+                     size=idfont, rotation=90., color='green')
+    ax_spec.set_xlim(0., len(arc_spec))
+    ax_spec.set_ylim(ymin, ymax*1.2)
+    ax_spec.set_xlabel('Pixel')
+    ax_spec.set_ylabel('Flux')
+    if title is not None:
+        ax_spec.text(0.04, 0.93, title, transform=ax_spec.transAxes,
+                     size='x-large', ha='left')#, bbox={'facecolor':'white'})
+    if ids_only:
+        plt.tight_layout(pad=0.2, h_pad=0.0, w_pad=0.0)
+        plt.savefig(outfile, dpi=800)
+        plt.close()
+        return
+
+    # Arc Fit
+    ax_fit = plt.subplot(gs[0, 1])
+    # Points
+    ax_fit.scatter(fit['xfit']*fit['xnorm'], fit['yfit'], marker='x')
+    if len(fit['xrej']) > 0:
+        ax_fit.scatter(fit['xrej']*fit['xnorm'], fit['yrej'], marker='o',
+                       edgecolor='gray', facecolor='none')
+    # Solution
+    xval = np.arange(len(arc_spec))
+    wave = func_val(fit['fitc'], xval/fit['xnorm'], 'legendre',
+                            minv=fit['fmin'], maxv=fit['fmax'])
+    ax_fit.plot(xval, wave, 'r-')
+    xmin, xmax = 0., len(arc_spec)
+    ax_fit.set_xlim(xmin, xmax)
+    ymin,ymax = np.min(wave)*.95,  np.max(wave)*1.05
+    ax_fit.set_ylim(np.min(wave)*.95,  np.max(wave)*1.05)
+    ax_fit.set_ylabel('Wavelength')
+    ax_fit.get_xaxis().set_ticks([]) # Suppress labeling
+    # Stats
+    wave_fit = func_val(fit['fitc'], fit['xfit'], 'legendre',
+                                minv=fit['fmin'], maxv=fit['fmax'])
+    rms = np.sqrt(np.sum((fit['yfit']-wave_fit)**2)/len(fit['xfit'])) # Ang
+    dwv_pix = np.median(np.abs(wave-np.roll(wave,1)))
+    ax_fit.text(0.1*len(arc_spec), 0.90*ymin+(ymax-ymin),
+                r'$\Delta\lambda$={:.3f}$\AA$ (per pix)'.format(dwv_pix), size='small')
+    ax_fit.text(0.1*len(arc_spec), 0.80*ymin+(ymax-ymin),
+                'RMS={:.3f} (pixels)'.format(rms/dwv_pix), size='small')
+    # Arc Residuals
+    ax_res = plt.subplot(gs[1,1])
+    res = fit['yfit']-wave_fit
+    ax_res.scatter(fit['xfit']*fit['xnorm'], res/dwv_pix, marker='x')
+    ax_res.plot([xmin,xmax], [0.,0], 'k--')
+    ax_res.set_xlim(xmin, xmax)
+    ax_res.set_xlabel('Pixel')
+    ax_res.set_ylabel('Residuals (Pix)')
+
+    # Finish
+    plt.tight_layout(pad=0.2, h_pad=0.0, w_pad=0.0)
+    plt.savefig(outfile, dpi=800)
+    plt.close()
+
+    plt.rcdefaults()
+
+    return
+
+

arclines/io.py

@@ -89,7 +89,7 @@ def load_line_list(line_file, add_path=False, use_ion=False, NIST=False):
                     except ValueError:
                         reli.append(0.)
         line_list.remove_column('Rel.')
-        line_list['Rel.'] = reli
+        line_list['RelInt'] = reli
         #
         gdrows = line_list['Observed'] > 0.  # Eliminate dummy lines
         line_list = line_list[gdrows]

arclines/load_source.py

@@ -69,7 +69,6 @@ def load_pypit(version, src_file, ions, plot=False, **kwargs):
       Additional lines
 
     """
-    from pypit import arutils
     # Load
     if version != 1:
         raise IOError("Unimplemented version!")
@@ -93,7 +92,7 @@ def load_pypit(version, src_file, ions, plot=False, **kwargs):
     mn_ID, mx_ID = min(ID_lines['wave']), max(ID_lines['wave'])
     # Unknown
     # Use PYPIT to decode
-    wave = arutils.func_val(pypit_fit['fitc'],
+    wave = arcl_utils.func_val(pypit_fit['fitc'],
                             np.array(pypit_fit['tcent'])/(npix-1),
                             pypit_fit['function'],
                             minv=pypit_fit['fmin'], maxv=pypit_fit['fmax'])

arclines/pypit_utils.py

@@ -6,6 +6,8 @@
 import numpy as np
 import pdb
 
+from arclines import utils as al_utils
+
 
 def find_peaks(censpec):
     """
@@ -21,7 +23,6 @@ def find_peaks(censpec):
     tampl, tcent, twid, w, detns
 
     """
-    from pypit import ararc
     fitp = 7  #slf._argflag['arc']['calibrate']['nfitpix']
     if len(censpec.shape) == 3:
         detns = censpec[:, 0].flatten()
@@ -34,8 +35,38 @@ def find_peaks(censpec):
                     (detns > np.roll(detns, 1)) & (detns >= np.roll(detns, -1)) &
                     (np.roll(detns, 1) > np.roll(detns, 2)) & (np.roll(detns, -1) > np.roll(detns, -2)) &
                     (np.roll(detns, 2) > np.roll(detns, 3)) & (np.roll(detns, -2) > np.roll(detns, -3)))[0]
-    tampl, tcent, twid = ararc.fit_arcspec(xrng, detns, pixt, fitp)
+    tampl, tcent, twid = fit_arcspec(xrng, detns, pixt, fitp)
     w = np.where((np.isnan(twid) == False) & (twid > 0.0) & (twid < 10.0/2.35) & (tcent > 0.0) & (tcent < xrng[-1]))
     # Return
     return tampl, tcent, twid, w, detns
 
+
+def fit_arcspec(xarray, yarray, pixt, fitp):
+
+    # Setup the arrays with fit parameters
+    sz_p = pixt.size
+    sz_a = yarray.size
+    ampl, cent, widt = -1.0*np.ones(sz_p, dtype=np.float),\
+                       -1.0*np.ones(sz_p, dtype=np.float),\
+                       -1.0*np.ones(sz_p, dtype=np.float)
+
+    for p in range(sz_p):
+        pmin = pixt[p]-(fitp-1)//2
+        pmax = pixt[p]-(fitp-1)//2 + fitp
+        if pmin < 0:
+            pmin = 0
+        if pmax > sz_a:
+            pmax = sz_a
+        if pmin == pmax:
+            continue
+        if pixt[p]-pmin <= 1 or pmax-pixt[p] <= 1:
+            continue  # Probably won't be a good solution
+        # Fit the gaussian
+        try:
+            popt = al_utils.func_fit(xarray[pmin:pmax], yarray[pmin:pmax], "gaussian", 3)
+            ampl[p] = popt[0]
+            cent[p] = popt[1]
+            widt[p] = popt[2]
+        except RuntimeError:
+            pass
+    return ampl, cent, widt

arclines/setup_package.py

@@ -2,4 +2,4 @@
 
 
 def get_package_data():
-    return {'pypit.tests': ['files/*']}
+    return {'arclines.tests': ['files/*']}

arclines/tests/test_io.py

@@ -45,6 +45,6 @@ def test_load_nist():
     # Stack
     line_lists = vstack(lists, join_type='exact')
     # Test
-    for key in ['wave','Aki','Rel.','Ion','NIST']:
+    for key in ['wave','Aki','RelInt','Ion','NIST']:
         assert key in line_lists.keys()