{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modelling XShooter data with xtool ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from xtool.data import XShooterData, Order\n",
    "from xtool.model import OrderModel, GenericBackground, MoffatTrace, VirtualPixelWavelength\n",
    "\n",
    "from scipy import sparse\n",
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Reading XShooter data ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xd = XShooterData('xtool_ds/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "current_order = xd[17]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generating a virtual pixel table for \"Wavelength\"-pixels ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "virt_pix = VirtualPixelWavelength.from_order(current_order)\n",
    "pixel_table = virt_pix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#### Initializing the two Models ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "background_mdl = GenericBackground(pixel_table, virt_pix.wavelength_pixels)\n",
    "trace_mdl = MoffatTrace(pixel_table, virt_pix.wavelength_pixels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "order_model = OrderModel([background_mdl, trace_mdl])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Show fittable parameters ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<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])>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "order_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Change fittable parameters ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Parameter('trace_pos', value=0.0, bounds=(-6, 6))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "order_model.trace_pos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "order_model.trace_pos = 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Parameter('trace_pos', value=10.0, bounds=(-6, 6))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "order_model.trace_pos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generating a model ####\n",
    "\n",
    "1. We generate a design matrix (https://en.wikipedia.org/wiki/Design_matrix)\n",
    "2. We solve the design matrix\n",
    "\n",
    "The evaluate does both of these steps at the same time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Generating the design matrix often depicted as capital A\n",
    "\n",
    "A, model_widths = order_model.generate_design_matrix(current_order, trace_pos=-5, sigma=1.5)\n",
    "\n",
    "# adding the uncertainties to the design matrix\n",
    "A.data /= current_order.uncertainty.compressed()[A.row]\n",
    "\n",
    "# making a vector of the result pixels often depicted as lower-case b\n",
    "\n",
    "b = current_order.data.compressed() / current_order.uncertainty.compressed()\n",
    "result = sparse.linalg.lsmr(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ -1.84607069e+03,   5.93872551e+05,   3.70733542e+06, ...,\n",
       "          1.52388567e+00,   1.24504355e+00,   6.61715836e-01]),\n",
       " 2,\n",
       " 8639,\n",
       " 4613.470857972648,\n",
       " 0.0024695645829768912,\n",
       " 0.5353327067182627,\n",
       " 13429.46405837754,\n",
       " 2820921905.777216)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy.sparse import coo_matrix\n",
    "\n",
    "#from http://stackoverflow.com/questions/22961541/python-matplotlib-plot-sparse-matrix-pattern\n",
    "\n",
    "def plot_coo_matrix(m):\n",
    "    if not isinstance(m, coo_matrix):\n",
    "        m = coo_matrix(m)\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111, axisbg='white')\n",
    "    ax.plot(m.col, m.row, 's', color='black', ms=1)\n",
    "    ax.set_xlim(0, m.shape[1])\n",
    "    ax.set_ylim(0, m.shape[0])\n",
    "    ax.set_aspect('auto')\n",
    "    for spine in ax.spines.values():\n",
    "        spine.set_visible(False)\n",
    "    ax.invert_yaxis()\n",
    "    ax.set_aspect('auto')\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10b2b6950>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADxdJREFUeJzt3U+I1eUXx/FnxilmEkFGsxlNjRaS1W9ENAskF2ZQgtFQ\nJCNI9AcSFCxsEdTGhS5y40IXBVmLWpSrNtNGaTMQ6GomAl3EKAw4hXHDzb2gzm1R/sav3/O9c/98\nv885z/O8X8szOXOCw+HzPN/7p6/ZbDoAgB/92g0AQEpYugDgEUsXADxi6QKARyxdAPCIpQsAHg0s\n8XNeTwYz+vr6crWdO3e6qakphW6ArAfns9ls5gfWkXQRiKtXr4r1r776ynMnQJ4QCBYK/9sl3hxB\n0oUJUsp1zjne3ANtBbPZaDabQ9IPSLowj4ULq44cOVL0o/8V/YCkC9NYuLBMms9169a5ubk5eXAd\nSReG/f3332L9hx9+8NwJkFcUCObm5lr/O5IurCLlwqo2ZpOki7CwcGHVd999J9ZnZ2fb+vckXZjz\n8ssvuwsXLuTqLFxY0GYgIOkiHNLC/eCDDxQ6AbLKOIGRdGEK1wqw6qGHHnJ37tzJ1Qtmk6QL+1i4\nsExauIcPH+7497B0YcKXX34p1uv1uudOgLyiQHDmzJnOfxfXC7BAGuoVK1a4W7duKXQDLOryBMb1\nAuwqGmoWLrQdP35crPdy5UXShapHHnlEvELgHhcWSIFg06ZNhZ96d/8/LfoBSReqpIXbzcMJoGxF\nJ7A2Fm7r30vShRZerQCrSphNki5sYeEiNJ9//nkpv4ekC+9Onz7tPvroo1ydhQsLSgoEhUmXpQvv\nSLmwqsTZ5HoBNrBwYdVnn30m1sueTZIuvNmyZYubmZnJ1Vm4sEAKBI8++qj7888/u/p1RT8g6cIb\naeGeOHFCoRMgq+gE1uXCbf23SLrwgWsFWPXiiy+6qampXL3H2STpQg8LF5ZJC/fbb7+t7O+RdFE5\naemycGFBhYGApAsd0lAPDAwodAJkaZ3AWLqoTNFQ375923MnQFaj0RDrN27cqPxvs3RRiZ9++kms\nc60AC4aGhnK1/v5+NzIyUvnf5k4XlZBS7pNPPul+//13hW6ARZ6uFbjThT9FQ83Chbbz58+LdZ8n\nMJIuSjUyMuL++OOPXJ1rBVggBYLt27e7y5cvl/6nCn/A0kWZpKG+dOmSe+655xS6ARZ5frUCSxfV\n400QsGrz5s3uypUruXqFs8nSRbUGBgbc3bt3c3UWLiyQAsGNGzeqfLUCD9JQnfn5eRYuzCo6gfl4\neZiEpIueSUPd19fnFhYWFLoBFileeZF0UY2ioWbhwirtExhLF13jXWewzOrnfrB00bW9e/fmaseO\nHVPoBMiy/Lkf3OmiK7w8DFadPXvWHTlyJFf3PJvc6aI84+PjYp2FCwukhfvee+8pdCIj6aJjUsqd\nnJx0r776qkI3wCJDJzDeHIFyGBpqIGP37t3u559/ztWVZpOli96xcGGZNJ/T09NubGxMoRuWLkrA\nd53BKoOBgAdp6I001IODgwqdAFkGF25LJF0sKbShRlqMnsBIuujOyZMnxbqBoQbEhbtq1SqFTtpH\n0kVL0lC/9dZb7vvvv1foBlhk/ATGgzR0zvhQI2E//vije/3113N1Q7PJ9QI68+abb4p1Q0ONhEkL\nd2JiQqGTzpF0IZJSbq1WcytXrlToBlgUyAmMpIv2FQ01CxfaXnnlFbFubOG2RNJFRiApAomS5nN2\ndtY98cQT/ptpjQdpaI/R1zwCoQUCrhewNGmotb68D7hff7+8qowu3JZIunDOBZcikJgAT2AkXRQ7\ndeqUWDc+1EiEtHDXr1+v0Ek5SLoQh/rtt99233zzjf9mgPsEfAIj6UJWNNQsXGj7+uuvxXoAC7cl\nlm7CDh48KNZDH2rE4d13383Vjh49qtBJubheSFiADyeQiICvFe7hegFZRUMNaNu/f79YD2jhtkTS\nTVAEKQIRk+azXq+H9k0lvCMNi7hWgFURBQKuF/Avaaj37dun0AmQFdHCbYmkm5BUhhrhaTQabmho\nKFcPeDZJuqmL9TWPiIO0cHft2qXQSfVIuokI6CPxkJhIT2Ak3ZQVDTULF9o+/vhjsR74wm2JpBu5\nZcuWuYWFhVw95qFGOCL+WiheMpYqXh4GqyK9VriH64UUSUP9+OOPK3QCZK1Zs0asR7JwWyLpRiry\nFIGARfjyMAlJNyXXrl0T65ENNQIlLdynnnpKoRMdJN0ISSl3enrajY2NKXQDLEroBMaDtFQkNNQI\nzNWrV8VEG+lscr2Qgm3bton1SIcagZEWbq1WU+hEF0k3Irw8DFYleAIj6cZOGuqRkRGFToCsjRs3\nivWIF25LJN0IJJgiEJBET2Ak3VjNz8+L9QSGGgGQFu6qVasUOrGDpRu40dHRXO3cuXMKnQBZRSew\nmzdveu7EFpZuwIqG+p133vHcCZB19uxZsc4JjDvdYI2OjopXCww1LJACwYULF9xLL72k0I0K3hwR\nm0QfTiAAPNh1zvEgLS7SUK9du1ahEyBreHhYrCe2cFsi6QaGFAHLOIH9H0k3BjMzM2I90aGGMdLC\nnZiYUOjENpJuQKShrtfrbnBwUKEbYBEnsBySbuiKhpqFC22cwDpD0g3A5s2b3ZUrV3J1hhoWcI8r\n4iVjIWOoYRXXCoW4XgiVNNR79uxR6ATIYuF2h6RrGEMNyziBtUTSDc2pU6fEOkMNC6SFOzk5qdBJ\neEi6RpEiYBUnsLaQdEMiDbX0tdWAb4cOHRLrLNz2kXSNIUXAMk5gbSPphoyhhgXSwt2wYYNCJ2Fj\n6RoiDfX09LRCJ0BW0Qns+vXrnjsJH0vXiKKhHhsb89wJkMXbfMvF0jVg27ZtYp2hhgVbtmzJ1S5d\nuqTQSRx4kGYADydgFQ92u8aDNKukod60aZNCJ0BW0UMyFm5vSLqKSBGwjBNYT0i61vBwApZJC/fT\nTz9V6CQ+JF0lpAhYxQmsFCRdS4qGGtD24YcfinUWbnlIup4NDw+7Wq2WqzPUsIATWGlIulawcGGV\ntHCffvpphU7ixtL1SBrqc+fOKXQCZBVdef3222+eO4kf1wue8HAClnGtUDquFzTt379frDPUsICF\n6xdJ1wOGGlZxAqsMSVeLNNQTExMKnQBZRd9GwsKtFkm3QqQIWMYJrFIkXd9Onjwp1hlqWCAtXD6u\n0Q+SbkVIEbCKE5gXJF2fpKHesWOHQidA1u7du8U6C9cfkm7JSBGwjBOYNyRdHxqNhlhnqGEBH9do\nA0m3RKQIWMUJzDuSbtWkoV62bJlCJ0DWJ598ItZZuDpIuiUgRcAyTmAqSLq+MdSwQFq4b7zxhkIn\nuIek2yNSBKziBKaKpFsFvnYHVs3Pz4t1Fq4+lm6XnnnmGbHOUMOC0dHRXI3ZtIHrhS5xrQCrpNlc\nsWKFu3XrlkI3yeJ6oUzSUJ84cUKhEyCr6MqLhWsHSbdDPJyAZZzAzCDplmHr1q1inaGGBSzcMJB0\nO8BQwypOYOaQdHslDfWZM2cUOgGyVq5cKdZZuDaRdNtAioBVjUZD/K4zZlMdSbdbfFgILJMWbq1W\nU+gE7SLpLoF7XFjFCcw0km43+NBnWDU4OCjWWbj2kXQLkCJgGScw80i6nTh+/LhYZ6hhgbRw6/W6\nQifoBklXQIqAVZzAgkHSbZc01Hv37lXoBMh67LHHxDoLNywk3fuQImAZJ7CgkHSXMj4+LtYZaljA\nPW48SLr/IUXAKk5gQSLptiIN9RdffKHQCZA1MjIi1lm44Uo+6ZIiYBknsGCRdCWnT58W6ww1LGDh\nxinppMtQwypOYMEj6T5IGurDhw8rdAJksXDjlmTSZahhGSewKJB07zl06JBYZ6hhAQs3fsklXYYa\nVnECiwpJ1zl5qKenpxU6AbJYuOlIJuky1LCK7zmLUtpJ97XXXhPrDDUsYOGmJYmkyz0urJJm84UX\nXnC//PKLQjcoUWHSjX7pFn06U9F3TAG+cOUVtTSvF4qGmoULbY1GQ6yzcOMXbdJ9+OGH3e3bt3N1\nhhoWcOUVvfSSLgsXVkkLd9euXQqdQEOUSZcUAau4x01GOklXGurly5crdAJk8RZ0OBdZ0iVFwDJO\nYEmJP+lOTU2JdYYaFkgLt1arKXQCbdEkXVIErOIElqS4k6401EePHlXoBMhav369WGfhpiv4pEuK\ngGWcwJIVZ9LdunWrWGeoYQELF5Kgky5DDauk2dywYYO7fv26QjdQEN8H3rBwYRVXXnCxXS9IQ93f\nH+T/CiIzMzMj1lm4uCe4pEuKgGWcwPCfOJLutWvXxDpDDQukhTs5OanQCSwLKumSImAVJzA8IPyk\nKw317OysQidA1urVq8U6CxeSIJIuKQKWcQKDINykOzo6KtYZaljAwkWnzC/d+fn5XI2hhgVF3+QL\ntGL6eoEUAau48sISwrtekIZ63759Cp0AWTt37hTrLFy0w2TSJUXAMk5gaEM4SZfvkYJlLFz0ylzS\nZahhlTSbzz77rPv1118VuoFxYXzKGAsXVnHlhQ7Zv16QhrronT6AT++//75YZ+GiGyaSLikClnEC\nQxfsJt3nn39erDPUsICFi7KpJ12GGlZJs7lu3To3Nzen0A0CY/NBGgsXVnHlhR7Zu16Qhnr58uUK\nnQBZa9euFessXJRBJemSImAZJzCUwE7SZeHCMhYuqqb+6gXnGGrYIC3cHTt2KHSCmHm9XiBFwCpO\nYCiZ/vWCNNT1et3XnwcKsXDhk5elWzTUg4ODPv48UOjixYtinYWLqlR+vUCKgGVceaEiOtcL58+f\nF+sMNSxg4UJDpUmXoYZVnMBQMf9Jl4ULq1i40FTJ0pWGevv27VX8KaAjLFxoK/16gaGGZZzA4Imf\n6wUWLixj4cKC0pbu+Pi4WGeoYQELF1aUdr3AUMMqTmBQUO31AgsXVrFwYU3PS1ca6mPHjvX6a4Ge\nsXBhUU/XCww1rLp48aLbs2dPrs5swpPyvyONhQvLuPKCsnLvdFm4sIyFC8s6XrqXL18W6ww1LODb\nH2Bdx9cLpAhYxQkMhpRzvcDChVUsXISi7aUrDfXw8HCpzQDdYOEiJG0t3aKh/uuvv0ptBujUgQMH\nxDoLF1YteadLioBlXHnBqO7udFm4sIyFixB1/JIxhhoW8PZzhKrl9UJfX1/mhyxcWEHKhXGF1wsD\nbf6CWrPZ5KUKsOSu+/ektuCcm2s2mxuV+wHastSDNABAiSr7NmAAQB5LFwA8YukCgEcsXQDwiKUL\nAB79A/DFIPPmWuHnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11bcdab50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plot_coo_matrix(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}