Merge pull request #3 from wkerzendorf/bolton_gaussian

Bolton gaussian

xtool/data.py

@@ -1,6 +1,8 @@
 import os
 from glob import glob
 
+import matplotlib as mpl
+import matplotlib.cm
 import numpy as np
 from astropy import units as u
 from astropy.io import fits
@@ -65,7 +67,12 @@ def xbin(self):
         -------
             : int
         """
-        return self.science_header.get("HIERARCH ESO DET WIN1 BINX", 2)
+
+        if "HIERARCH ESO DET WIN1 BINX" not in self.science_header:
+            print "WARNING: XBin not known defaulting to 1"
+            return 1
+        else:
+            return self.science_header.get("HIERARCH ESO DET WIN1 BINX")
 
     @property
     def ybin(self):
@@ -74,8 +81,12 @@ def ybin(self):
         -------
             : int
         """
-
-        return self.science_header.get("HIERARCH ESO DET WIN1 BINY", 2)
+
+        if "HIERARCH ESO DET WIN1 BINY" not in self.science_header:
+            print "WARNING: YBin not known defaulting to 1"
+            return 1
+        else:
+            return self.science_header.get("HIERARCH ESO DET WIN1 BINY")
 
     def read_data(self, xshooter_dir, data_fname, ext=0):
         """
@@ -190,7 +201,12 @@ def read_spectral_format_table(self):
         spectral_format = Table.read(
             os.path.join(xtool_data_path, '{0}_{1}.fits'.format(
                 self.spectral_format_fname, self.instrument_arm))).to_pandas()
-        return spectral_format.set_index['ORDER']
+        spectral_format['ORDER'] = spectral_format['ORDER'].astype(int)
+
+        spectral_format = spectral_format.set_index('ORDER')
+        spectral_format.index = spectral_format.index.rename('ABSORDER')
+        return spectral_format
+
 
     def read_order_table(self):
         """
@@ -209,7 +225,6 @@ def __repr__(self):
         return "<XShooterData {0} {1}>".format(self.science_header['ARCFILE'],
                                                self.instrument_arm)
 
-
     def get_order_coefficients(self, order_id):
         """
         Getting the polynomial coefficients for a specific order
@@ -233,9 +248,16 @@ def get_order_coefficients(self, order_id):
         """
 
         order_data = self.order_table.ix[order_id]
+        spectral_format_table = self.spectral_format_table.ix[order_id]
         deg = np.int(order_data['DEGY'])
-        slice_y = slice(np.int(order_data['STARTY']) - 1,
-                        np.int(order_data['ENDY']) - 1)
+        if self.instrument_arm == 'NIR' and order_id > 15:
+            slice_y = slice(np.int(spectral_format_table['YMIN']) - 1 + 30,
+                            np.int(spectral_format_table['YMAX']) - 1 - 30)
+        elif 'YMIN' in spectral_format_table:
+            slice_y = slice(np.int(spectral_format_table.get('YMIN') - 1),
+                            np.int(spectral_format_table.get('YMAX') - 1))
+        else:
+            slice_y = slice(None, None)
 
         edg_up_coef = [order_data['EDGUPCOEF{0}'.format(i)]
                        for i in xrange(deg + 1)]
@@ -265,7 +287,7 @@ def calculate_edges(self, order_id):
         """
 
         slice_y, edg_up_coef, edg_lo_coef = self.get_order_coefficients(order_id)
-        yrange = (np.arange(self.science_data.shape[0]))
+        yrange = (np.arange(self.science_data.shape[0]) * self.ybin)
 
 
         lower_edges = np.polynomial.polynomial.polyval(
@@ -336,9 +358,12 @@ def cutout_and_mask_order(self, order_id, data_array, extra_sample=5):
         slice_y, lower_edge, upper_edge = self.calculate_edges(order_id)
         y, x = np.mgrid[
                :cutout_data_array.shape[0], :cutout_data_array.shape[1]]
+
         mask = ((x > (lower_edge - edge_slice.start)[None].T) &
-                (x < (upper_edge - edge_slice.start)[None].T) &
-                (y > slice_y.start) & (y < slice_y.stop))
+                (x < (upper_edge - edge_slice.start)[None].T))
+
+        if slice_y.start is not None and slice_y.stop is not None:
+            mask = mask & (y > slice_y.start) & (y < slice_y.stop)
 
         return cutout_data_array, mask
 
@@ -355,11 +380,41 @@ def get_mask(self, data_array):
         return mask
 
 
+    def display_data(self, figure, vmin=0, vmax=500,
+                     cmap=matplotlib.cm.gray):
+        """
+        Displau the XShooter data in
 
+        Parameters
+        ----------
+        figure : matplotlib.Figure
+        vmin : minimum cut [default=0]
+        vmax : maximum cut [default=500]
+        cmap : matplotlib colormaps [default 'gray']
 
+        """
+        ax = figure.add_subplot(111)
+        ax.imshow(self.science_data, vmin=vmin, vmax=vmax, cmap=cmap)
+        for order_id in self.order_table.index:
+            slice_y, low_edge, up_edge = self.calculate_edges(order_id)
+#            slice_y2 = slice(
+#                self.spectral_format_table.loc[order_id, 'YMIN'] - 1 + 50,
+#                self.spectral_format_table.loc[order_id, 'YMAX'] - 1 - 50)
+
+            ax.plot(low_edge, np.arange(len(low_edge)), color='blue', lw=2)
+            ax.plot(up_edge, np.arange(len(up_edge)), color='red', lw=2)
+            if slice_y.start is not None and slice_y.stop is not None:
+                ax.plot([low_edge[slice_y.start], up_edge[slice_y.start]],
+                     [slice_y.start, slice_y.start], color='purple', lw=2)
+                ax.plot([low_edge[slice_y.stop], up_edge[slice_y.stop]],
+                     [slice_y.stop, slice_y.stop], color='purple', lw=2)
+            ax.text(0.5 * (low_edge[len(low_edge)/2] +
+                           up_edge[len(up_edge)/2]), len(low_edge)/2, str(order_id),
+                 verticalalignment='center', horizontalalignment='center',
+                 bbox=dict(facecolor='white'))
 
-class Order(object):
 
+class Order(object):
     @classmethod
     def from_xshooter_data(cls, xshooter_data, order_id):
         """
@@ -446,20 +501,35 @@ def _update_mask(self, mask):
         ----------
         mask : numpy.ndarray
         """
-        self.data.mask = mask
-        self.uncertainty.mask = mask
-        self.flags.mask = mask
-        self.wcs._update_mask(mask)
+        updated_mask = self.data.mask | mask
+
+        self.data.mask = updated_mask
+        self.uncertainty.mask = updated_mask
+        self.flags.mask = updated_mask
+        self.wcs._update_mask(updated_mask)
 
 
 
     def enable_flags_as_mask(self):
         flags_mask = (self.flags.filled() == 0) & self.order_mask
         self._update_mask(~flags_mask)
 
+    def enable_instrument_model_masking(self):
+        instrument_model_mask = ((self.wcs.pix_to_wave_ma.data > 0.0) &
+                                 self.order_mask)
+        self._update_mask(~instrument_model_mask)
 
+    def standard_masking(self):
+        """
+        Flags everything with QUAL flag != 0 and everything less
+        than 0.
+
+        """
 
+        self.enable_flags_as_mask()
+        self.enable_instrument_model_masking()
 
+        self._update_mask(~((self.data.filled() > 0) & self.order_mask))
 
 
     def __repr__(self):

xtool/instrument/base.py

@@ -0,0 +1 @@
+### This module lets

xtool/model/__init__.py

@@ -1,3 +1,4 @@
 from xtool.model.base import OrderModel, VirtualPixelWavelength
-from xtool.model.order_models import (GenericBackground, MoffatTrace,
-                                      SlopedMoffatTrace)
+from xtool.model.virtual_pixel_order_models import (
+    GenericBackground, MoffatTrace, SlopedMoffatTrace, PolynomialBackground,
+    PolynomialMoffatTrace)

xtool/model/base.py

@@ -5,6 +5,7 @@
 from astropy import modeling
 from scipy import sparse
 from scipy.spatial import cKDTree
+from scipy.optimize import lsq_linear
 
 from xtool.wcs import PolynomialOrderWCS
 
@@ -23,7 +24,8 @@ class VirtualPixelWavelength(modeling.Model):
                                     'NIR' : 0.03}
 
     @classmethod
-    def from_order(cls, order, poly_order=(2, 3), wavelength_sampling=None):
+    def from_order(cls, order, poly_order=(2, 3), wavelength_sampling=None,
+                   sub_sampling=5):
         """
         Instantiate a Virtualpixel table from the order raw Lookup Table WCS
         and then fitting a Polynomial WCS with given orde
@@ -36,7 +38,7 @@ def from_order(cls, order, poly_order=(2, 3), wavelength_sampling=None):
             float for wavelength spacing. If `None` will use for different arms
             * UVB - 0.04 nm
             * VIS - 0.04 nm
-            * NIR - 0.1 nm
+            * NIR - 0.03 nm
 
         Returns
         -------
@@ -55,7 +57,8 @@ def from_order(cls, order, poly_order=(2, 3), wavelength_sampling=None):
             wavelength_sampling)
 
         return VirtualPixelWavelength(polynomial_wcs, order.wcs,
-                                      wavelength_bins)
+                                      wavelength_bins,
+                                      sub_sampling=sub_sampling)
 
     def __init__(self, polynomial_wcs, lut_wcs, raw_wavelength_pixels,
                  sub_sampling=5):
@@ -179,10 +182,10 @@ def evaluate(self, wave_transform_coef):
                 pixel_table.wavelength_pixel_id))
 
         pixel_table['wavelength'] = self.wavelength_pixels[
-            pixel_table.wavelength_pixel_id]
+            pixel_table.wavelength_pixel_id].astype(np.float64)
 
         pixel_table['slit_pos'] = self.lut_wcs.pix_to_slit[
-            pixel_table.pixel_id]
+            pixel_table.pixel_id].astype(np.float64)
 
         pixel_table['normed_wavelength'] = (
                 (pixel_table.wavelength - pixel_table.wavelength.min()) /
@@ -194,21 +197,27 @@ def evaluate(self, wave_transform_coef):
 
 class OrderModel(object):
 
-    def __init__(self, model_list):
+    def __init__(self, model_list, virtual_pixel=None):
         self.parameter_dict = self._generate_parameter_dict(model_list)
         self.model_list = model_list
+        self.virtual_pixel = virtual_pixel
 
 
     def __getitem__(self, item):
         return self.model_list[item]
 
     @property
     def fittable_parameter_names(self):
-        return self._get_variable_normal_parameters()
+
+        fittable_parameter_names = self._get_variable_normal_parameters()
+        if self.virtual_pixel is not None:
+            fittable_parameter_names += ['wave_transform_coef']
+        return fittable_parameter_names
+
 
     @property
     def fittable_parameter_dict(self):
-        return OrderedDict([(param_name, getattr(self, param_name).value)
+        return OrderedDict([(param_name, getattr(self, param_name))
                             for param_name in self.fittable_parameter_names])
 
     def _get_variable_normal_parameters(self):
@@ -293,6 +302,29 @@ def result_dict(self):
         return result_dict
 
     def generate_design_matrix(self, order, **kwargs):
+        """
+        Generate the Design matrix for solving the Linear Least Squares problem
+        Parameters
+        ----------
+        order : xtool.data.Order objects
+        kwargs : fittable kwargs
+
+        Returns
+        -------
+        """
+
+        if 'wave_transform_coef' in kwargs:
+
+            wave_transform_coef = (
+                kwargs['wave_transform_coef'].reshape(
+                    self.virtual_pixel.wave_transform_coef.shape))
+            self.virtual_pixel.wave_transform_coef = wave_transform_coef
+            pixel_table = self.virtual_pixel()
+            for model in self.model_list:
+                if hasattr(model, 'update_pixel_table'):
+                    model.update_pixel_table(pixel_table)
+
+
         design_matrices = []
         uncertainties = order.uncertainty.compressed()
         for model in self.model_list:
@@ -334,20 +366,26 @@ def solve_design_matrix(self, dmatrix, order, solver='lsmr',
         b = order.data.compressed() / order.uncertainty.compressed()
         if solver == 'lsmr':
             result = sparse.linalg.lsmr(dmatrix.tobsr(), b, **solver_dict)
+        elif solver == 'lsq':
+            lsq_dict = dict(bounds=(0, np.inf), lsmr_tol='auto', verbose=1)
+            lsq_dict.update(solver_dict)
+            result = lsq_linear(dmatrix, b,  **lsq_dict)
+            result = [result.x, result]
+
         else:
             raise NotImplementedError('Solver {0} is not implemented'.format(
                 solver))
         return result
 
-    def set_matrix_parameters(self, b, model_widths):
+    def set_matrix_parameters(self, result, model_widths):
         matrix_model_columns = np.cumsum(model_widths)
         for i, model in enumerate(self.model_list):
             current_matrix_column = (
                 matrix_model_columns[i] - matrix_model_columns[0])
             for (matrix_parameter,
                  matrix_slice) in model.matrix_parameter_slices.items():
                 setattr(model, matrix_parameter,
-                        b[matrix_slice.start + current_matrix_column:
+                        result[matrix_slice.start + current_matrix_column:
                         matrix_slice.stop + current_matrix_column])
 
     def set_matrix_uncertainties(self, order):
@@ -372,7 +410,6 @@ def evaluate(self, order, solver='lsmr', solver_dict={}, **kwargs):
         dmatrix, model_widths = self.generate_design_matrix(order, **kwargs)
         result = self.solve_design_matrix(dmatrix, order, solver=solver,
                                           solver_dict=solver_dict)
-
         self.set_matrix_parameters(result[0], model_widths)
         return (dmatrix * result[0]) * order.uncertainty.compressed()
 
@@ -384,10 +421,18 @@ def evaluate_to_frame(self, order, solver='lsmr', solver_dict={}, **kwargs):
         return result_frame
 
     def evaluate_residuals(self, order, solver='lsmr', solver_dict={},
-                           **kwargs):
+                           ylim=None, **kwargs):
         result = self.evaluate(order, solver, solver_dict, **kwargs)
         b = order.data.compressed()
-        return ((result - b) / order.uncertainty.compressed())
+
+        if ylim is not None:
+            row_filter = (order.wcs.y > ylim[0]) & (order.wcs.y < ylim[1])
+            result = result[row_filter]
+            b = b[row_filter]
+            uncertainty = order.uncertainty.compressed()[row_filter]
+        else:
+            uncertainty = order.uncertainty.compressed()
+        return ((result - b) / uncertainty)
 
 
     def evaluate_chi2(self, order, solver='lsmr', solver_dict={}, **kwargs):