Xtool¶

This is a package designed to model spectra from the XShooter spectrograph.

A example dataset is available from here

Installation¶

For installation of xtool we fully rely on the anaconda installation environment.

The anaconda environment already brings a number of packages - for a barebone installation please refer to miniconda.

Currently, xtool is designed to work with Python 2.7. We have provided an environment that contains all the necessary dependencies for xtool:

curl -O https://raw.githubusercontent.com/eso/xtool/master/xtool_env.yml
conda env create -n <yourname_for_the_environment> --file xtool_env.yml
source activate <yourname_for_the_environment>

Most of the analysis is performed in notebooks (one can think of them as cookbooks):

curl -O https://raw.githubusercontent.com/eso/xtool/master/docs/notebooks/cookbook_reading.ipynb
jupyter notebook cookbook_reading.ipynb

Generating Data for XTool¶

(provided by Sabine Möhler)

download science data with associated raw calibrations

(UVB XSHOO.2013-03-04T07:06:20.957, VIS XSHOO.2013-03-04T07:06:26.098)
process with reflex as STARE
look for SCI_SLIT_DIVFF, SCI_SLIT_WAVE_MAP, SCI_SLIT_SLIT_MAP in

<reflex_data_dir>/reflex_tmp_products/xshooter/xsh_scired_slit_stare_1/<date-time>
the slit wavelength calibration raw data are delivered together with

the other raw data, but their processing is not supported by reflex. I used the calSelector version on my disk to get the master calibrations used by QC to process these data and processed them manually. This way I got the “ON” frame in the correct orientation.

Reading XShooter data with xtool¶

In [1]:

from xtool.data import XShooterData

In [2]:

xd = XShooterData('xtool_ds') # with the absolute or relative path given
xd

Out[2]:

<XShooterData XSHOO.2013-03-04T07:06:29.289.fits NIR>

You can access the individual orders given in the dataset by just indexing the dataset¶

In [3]:

xd[11]

Out[3]:

<Order 11 XSHOO.2013-03-04T07:06:29.289.fits NIR>

In [4]:

xd[5]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-4-dee1a2b27f4d> in <module>()
----> 1 xd[5]

/Users/wkerzend/projects/eso/xtool/xtool/data.pyc in __getitem__(self, item)
     52         if item not in self.order_table.index:
     53             raise IndexError('Dataset only contains orders {0} - {1}'.format(
---> 54                 self.order_table.index.min(), self.order_table.index.max()))
     55         else:
     56             return Order.from_xshooter_data(self, item)

IndexError: Dataset only contains orders 11 - 26

Each order contains the cutout of the order in masked arrays that cutout intraorder regions.¶

In [5]:

order = xd[11]

Each one of the order objects also contains a WCS that can transform between pixelspace and woorld coordinate system (angstrom, slit position)

In [6]:

order.wcs(51, 52)

Out[6]:

(2474.587158203125, 0.5863650441169739)

Modelling XShooter data with xtool¶

In [1]:

from xtool.data import XShooterData, Order
from xtool.model import OrderModel, GenericBackground, MoffatTrace, VirtualPixelWavelength

from scipy import sparse
from scipy import optimize

Reading XShooter data¶

In [2]:

xd = XShooterData('xtool_ds/')

In [3]:

current_order = xd[17]

Generating a virtual pixel table for “Wavelength”-pixels¶

In [4]:

virt_pix = VirtualPixelWavelength.from_order(current_order)
pixel_table = virt_pix()

Initializing the two Models¶

In [7]:

background_mdl = GenericBackground(pixel_table, virt_pix.wavelength_pixels)
trace_mdl = MoffatTrace(pixel_table, virt_pix.wavelength_pixels)

In [8]:

order_model = OrderModel([background_mdl, trace_mdl])

Show fittable parameters¶

In [9]:

order_model

Out[9]:

<OrderModel(background_level=[ 0.  0.  0. ...,  0.  0.  0.], amplitude=[ nan  nan  nan ...,  nan  nan  nan], trace_pos=0.0, sigma=1.0, beta=1.5 [f])>

Change fittable parameters¶

In [10]:

order_model.trace_pos

Out[10]:

Parameter('trace_pos', value=0.0, bounds=(-6, 6))

In [11]:

order_model.trace_pos = 10.

In [12]:

order_model.trace_pos

Out[12]:

Parameter('trace_pos', value=10.0, bounds=(-6, 6))

Generating a model¶

We generate a design matrix (https://en.wikipedia.org/wiki/Design_matrix)
We solve the design matrix

The evaluate does both of these steps at the same time

In [15]:

# Generating the design matrix often depicted as capital A

A, model_widths = order_model.generate_design_matrix(current_order, trace_pos=-5, sigma=1.5)

# adding the uncertainties to the design matrix
A.data /= current_order.uncertainty.compressed()[A.row]

# making a vector of the result pixels often depicted as lower-case b

b = current_order.data.compressed() / current_order.uncertainty.compressed()
result = sparse.linalg.lsmr(A, b)

In [16]:

result

Out[16]:

(array([ -1.84607069e+03,   5.93872551e+05,   3.70733542e+06, ...,
          1.52388567e+00,   1.24504355e+00,   6.61715836e-01]),
 2,
 8639,
 4613.470857972648,
 0.0024695645829768912,
 0.5353327067182627,
 13429.46405837754,
 2820921905.777216)

In [17]:

import matplotlib.pyplot as plt
from scipy.sparse import coo_matrix

#from http://stackoverflow.com/questions/22961541/python-matplotlib-plot-sparse-matrix-pattern

def plot_coo_matrix(m):
    if not isinstance(m, coo_matrix):
        m = coo_matrix(m)
    fig = plt.figure()
    ax = fig.add_subplot(111, axisbg='white')
    ax.plot(m.col, m.row, 's', color='black', ms=1)
    ax.set_xlim(0, m.shape[1])
    ax.set_ylim(0, m.shape[0])
    ax.set_aspect('auto')
    for spine in ax.spines.values():
        spine.set_visible(False)
    ax.invert_yaxis()
    ax.set_aspect('auto')
    ax.set_xticks([])
    ax.set_yticks([])
    return ax

In [19]:

%matplotlib inline

plot_coo_matrix(A)

Out[19]:

<matplotlib.axes._subplots.AxesSubplot at 0x10b2b6950>

_images/notebooks_cookbook_modeling_20_1.png

Fitting XShooter data with xtool¶

In [1]:

from xtool.data import XShooterData, Order
from xtool.model import OrderModel, VirtualPixelWavelength, GenericBackground, MoffatTrace, SlopedMoffatTrace

from collections import OrderedDict
from scipy import sparse
from scipy import optimize
import numpy as np

Reading XShooter data¶

In [2]:

xd = XShooterData('xtool_ds/')

In [3]:

current_order = xd[17]
current_order.enable_flags_as_mask() # important to update the data to include the bad pixels from the pipeline

Generating a virtual pixel table for “Wavelength”-pixels¶

In [4]:

virt_pix = VirtualPixelWavelength.from_order(current_order, wavelength_sampling=0.03)
pixel_table = virt_pix()

Initializing the two Models¶

In [5]:

background_mdl = GenericBackground(pixel_table, virt_pix.wavelength_pixels)
trace_mdl = SlopedMoffatTrace(pixel_table, virt_pix.wavelength_pixels)

In [6]:

order_model = OrderModel([background_mdl, trace_mdl])

Show fittable parameters¶

In [7]:

order_model

Out[7]:

<OrderModel(background_level=[    0.     0.  2023. ...,  2008.  2008.  2008.], amplitude=[ nan  nan  nan ...,  nan  nan  nan], trace_pos=0.0, trace_slope=0.0, sigma=1.0, beta=1.5 [f])>

Initializing the fitter¶

In [8]:

from xtool.model.fitters import Fitter

In [9]:

fitter = Fitter(order_model, current_order)

Differential Evolution fit¶

Differential evolution is a very thorough but slow processs. It does not use starting values but moves within the bounds of the parameters.

In [31]:

order_model.beta.fixed = False
dresult = fitter.fit_differential_evolution(disp=True)

differential_evolution step 1: f(x)= 4.22785e+06
differential_evolution step 2: f(x)= 4.22785e+06
differential_evolution step 3: f(x)= 2.95083e+06
differential_evolution step 4: f(x)= 2.95083e+06
differential_evolution step 5: f(x)= 2.26008e+06
differential_evolution step 6: f(x)= 2.26008e+06
differential_evolution step 7: f(x)= 2.26008e+06
differential_evolution step 8: f(x)= 2.26008e+06
differential_evolution step 9: f(x)= 2.06608e+06
differential_evolution step 10: f(x)= 2.06231e+06
differential_evolution step 11: f(x)= 2.00095e+06
differential_evolution step 12: f(x)= 1.96445e+06
differential_evolution step 13: f(x)= 1.96445e+06
differential_evolution step 14: f(x)= 1.96445e+06
differential_evolution step 15: f(x)= 1.89723e+06
differential_evolution step 16: f(x)= 1.89723e+06
differential_evolution step 17: f(x)= 1.89723e+06
differential_evolution step 18: f(x)= 1.89723e+06
differential_evolution step 19: f(x)= 1.89654e+06
differential_evolution step 20: f(x)= 1.89654e+06
differential_evolution step 21: f(x)= 1.89654e+06
Fit finished in 1475.61700892 s

In [32]:

dresult

Out[32]:

fun: 1896544.4149950482
 message: 'Optimization terminated successfully.'
    nfev: 1425
     nit: 21
 success: True
       x: array([-2.44188896, -0.26832235,  0.19794569,  2.0046534 ])

Nelder Mead / Simplex¶

This one is much faster algorithm that will use starting values, but is unbounded

In [10]:

order_model.trace_pos = -2
order_model.sigma = 0.1
order_model.trace_slope = 0.
order_model.beta.fixed = False
order_model.beta = 2.0
result = fitter.fit_scipy_minimize('Nelder-Mead')

/Users/wkerzend/anaconda3/envs/xtool/lib/python2.7/site-packages/scipy/optimize/_minimize.py:394: RuntimeWarning: Method Nelder-Mead cannot handle constraints nor bounds.
  RuntimeWarning)

Fit finished in 851.559898853 s

/Users/wkerzend/anaconda3/envs/xtool/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:243: SparseEfficiencyWarning: splu requires CSC matrix format
  warn('splu requires CSC matrix format', SparseEfficiencyWarning)
/Users/wkerzend/anaconda3/envs/xtool/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:161: SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format
  'is in the CSC matrix format', SparseEfficiencyWarning)

In [15]:

spec = order_model.model_list[1].to_spectrum()

In [19]:

%pylab inline
plot(spec.wavelength, spec.flux)

Populating the interactive namespace from numpy and matplotlib

Out[19]:

[<matplotlib.lines.Line2D at 0x1231a7a50>]

_images/notebooks_cookbook_fitting_21_2.png

Plotting in DS9¶

In [11]:

import pyds9

In [12]:

d = pyds9.DS9()

In [14]:

order_model.trace_pos = result.x[0]
order_model.trace_slope = result.x[1]
order_model.sigma = result.x[2]
order_model.beta = result.x[3]
model = order_model.evaluate_to_frame(current_order, trace_pos=order_model.trace_pos,
                                      trace_slope=order_model.trace_slope, sigma=order_model.sigma,
                                      beta=order_model.beta.value)
d.set('frame 1')
d.set_np2arr(current_order.data.filled())

d.set('frame 2')
d.set_np2arr(model.filled())

d.set('frame 3')
d.set_np2arr((current_order.data.filled() - model.filled()) / current_order.uncertainty.filled())

Out[14]:

Navigation

Xtool¶

Installation¶

Generating Data for XTool¶

Reading XShooter data with xtool¶

You can access the individual orders given in the dataset by just indexing the dataset¶

Each order contains the cutout of the order in masked arrays that cutout intraorder regions.¶

Modelling XShooter data with xtool¶

Reading XShooter data¶

Generating a virtual pixel table for “Wavelength”-pixels¶

Initializing the two Models¶

Show fittable parameters¶

Change fittable parameters¶

Generating a model¶

Fitting XShooter data with xtool¶

Reading XShooter data¶

Generating a virtual pixel table for “Wavelength”-pixels¶

Initializing the two Models¶

Show fittable parameters¶

Initializing the fitter¶

Differential Evolution fit¶

Nelder Mead / Simplex¶

Plotting in DS9¶