{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses all the raw data from the masterlist, maps, PSF and relevant MOCs to create XID+ prior object and relevant tiling scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in MOCs\n",
    "The selection functions required are the main MOC associated with the masterlist. As the prior for XID+ is based on IRAC detected sources coming from two different surveys at different depths (SERVS and SWIRE) I will split the XID+ run into two different runs. Here we use the SERVS depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../data/Lockman-SWIRE/holes_Lockman-SWIRE_irac1_O16_MOC.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Masterlist\n",
    "Next step is to read in Masterlist and select only sources that are detected in mid-infrared and at least one other wavelength domain (i.e. optical or nir). This will remove most of the objects in the catalogue that are artefacts. We can do this by using the `flag_optnir_det` flag and selecting sources that have a binary value of $>= 5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "\n",
    "masterlist=fits.open('../data/Lockman-SWIRE/master_catalogue_lockman-swire_20171129.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=masterlist[1].data['flag_optnir_det']>=5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create uninformative (i.e. conservative) upper and lower limits based on IRAC fluxes\n",
    "As the default flux prior for XID+ is a uniform distribution, it makes sense to set reasonable upper and lower 24 micron flux limits based on the longest wavelength IRAC flux available. For a lower limit I take IRAC/500.0 and for upper limit I take IRACx500."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_lower=np.full(good.sum(),0.0)\n",
    "MIPS_upper=np.full(good.sum(),1E5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(250000,good.sum()):\n",
    "    if masterlist[1].data['f_irac_i4'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac_i4'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac_i4'][good][i]*500.0\n",
    "    elif masterlist[1].data['f_irac_i3'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac_i3'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac_i3'][good][i]*500.0\n",
    "    elif masterlist[1].data['f_irac_i2'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac_i2'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac_i2'][good][i]*500.0\n",
    "    elif masterlist[1].data['f_irac_i1'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac_i1'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac_i1'][good][i]*500.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "MIPS_Map=fits.open('../data/Lockman-SWIRE/wp4_lockman-swire_mips24_map_v1.0.fits.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_psf=fits.open('../../dmu17/dmu17_Lockman-SWIRE/dmu17_MIPS_Lockman-SWIRE_20171122.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "centre=np.long((MIPS_psf[1].header['NAXIS1']-1)/2)\n",
    "radius=20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x25b179d30>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASgAAAD7CAYAAADZ2gksAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+wXVWV5z/v5eXlB7xHwo8YKbubUsstWo62UKU9NBAB\nZ8QqTdNabYkBk0BTokNJ7BhoC6X80aVkQqA7PYQRCJEfUdog6jhFa1XTCK3D0BHF5seskdaedmBQ\n0oT8fCEk780f9weXe9e+d5977rnv3Pe+n6pXde56Z5+99znnrrv3WnuvNTQ1NYUQQpSR4elugBBC\nxJCCEkKUFikoIURpkYISQpQWKSghRGmRghJClJaR6W6AEGLmEUIYBm4A3gq8CFxsZk9lvY5GUEKI\nIvgjYL6Z/QFwJXBtNxeRghJCFMEfAn8LYGYPAad2c5GupnhZh287duyoL1d/85vfzOOPP87TTz/t\nnnvgwIEW2UsvvdQiGxoacsuPjo62yIaHW/Xw5OSkW96T79+/v0U2MTFRPz7//PPZtm0bhw4dcq/p\nyb0+HT582C3vnXvkyJEkGcT72o5169axfv366P+9+x97Jt7992Rz5851y4+MtL6mc+bMyVx+5cqV\nbN26FYD58+e3nHf00Ue3yJYsWeJe8+1vf7vTpqPcc4tk8eKF7NpV+c6ccMKY/wAyMDQ0lLy1ZGpq\nql1948Duhs9HQggjZua/5BG6HUF1PXxbsGBBl1WWl+OOO266m9BzXv3qV093E3rO8ccfP91N6Dkj\nI62KuiTsAcYaPg9nVU7QvYLqyfBNCFEuhoaGkv868CPgvQAhhHcC/9RVe7rZLBxCuBm428zurX7+\nV+C1MQ05MTExNRNHTkKUjNxTvJGRkWSFcPjw4Wh9DWagf1dt1yoz+1+Z25O1QJVMw7fHH3+8fnzq\nqaeyY8eOGWWDuuyyy9i0adOMskFdd911rFmzJvr/QbRBrV27lg0bNgAzxwZ1wgljPPfc3vpxXrzn\n0g1mNgl8LO91ulVQPwLeB/xNyvDtiSeeqB+feuqpPPHEE9Evo/dl8kZ5sZfRu673xfW+9OArE09p\nHjx48BWf9+zZE72mJ0+Vgd8n7z5lUbop1PrtvbSeMoq93KkKKvZOpCqj2P1rVHC1Pnl1pb47AA89\n9FCLzFNwsZmD16fXve7N7rn9JGHq1le6VVD3AO8OIfyY6vCtd00SQkwXM0JB9Wr4JoQoFzNCQQkh\nZiZSUEKI0uLZxqYTKSghRJ1ZOYLau3dvy+fY+qvUdVkxj02qdyfmMWr2zsVkzd6+iYmJ6DKD1Pqz\nLDPI4sVLpfnljPXHOxeyefG8X+rYr7cn9+5VrP7G8rUlI17fPFnsmbz44ostsoULF7bIPM8ewLx5\n81pkY2P/2iJbsuR33fJFMSsVlBBiMJCCEkKUFikoIURpkZFcCFFaZuUIat++fW0/d0NM06fuu4sZ\n2VON5M3lDxw4EDWoevLU/XUxeeqWoJg8ZS9drD+p5dvJUw3nWc5NMdLX9lB69z+LkTz1PfFkAGNj\nrfvm9uzZ0yKLbAUsjFmpoIQQg4EUlBCitEhBCSFKixSUEKK0zEovXrOh8ODBg8krxiHbqmXvXM+g\nnMX4mRK76dChQ5niOeUNOOfdv7z3NNae1F/VLEZyr0+x/qfGk0oJmFd7vl4QPE8Wa1MRq/sbgyBO\nFxpBCSFKixSUEKK0SEEJIUqLFJQQorTMGCN5COERKtldAH5lZopLLsSAMyNGUCGE+cCQmS1LOb/Z\nOzExMRH1OHk3yDs3S+yhLB6z1NhNzdsiDh06lGmrSur2lXbyZrrJcZhyvdTrZnm5s2yV8frfbTyq\n2rPM6xlNfSdjbfLSo3nZg/rNjFBQwFuBhSGEH1Sv8ZlqhmEhxABTNgXVbZa+A8AG4D9Sye5yZwhB\n9iwhBpzh4eHkv37QberzeVSyCU9UPz8MfMDMfu2d/+yzz04tXbo0V0OFEB3JPfw5+eSTkxXCk08+\nWfhwq9tRz2rgLcDHQwgnAuPA/4udfP3119ePv/KVr3DllVfOKBvUV7/6VS655JKBskF1Gsp//etf\n58Mf/nDy9VKu2encLCvRu7FBbdmyhdWrV0fP9eLZe6vLY+d6WYRjmYXHx8dbZK997WtbZO94x5lu\n+Rq9Tn0+U7x4twBbQwj/AEwBq80sGjzIMyjnNejm3VYR+9J3m+Dg8OHD0WvmjeeUSl77QXP52ude\nG9+zXjM1nlWKQbv2LIq41148qdgX3ku6oK0urXSbWfgQcH6P2yKEmGb6ZVtKRYZtIUSdGTGCEkLM\nTKSghBClZVZO8TyDchFG0hhZjOSpsZuay09OTkavWcRK7Lyk1FU7p2y/qjFi97nxudSO8z4TT+4Z\nxLPECIslWOgnM8WLJ4SYgZTtx0gKSghRZ1ZO8YQQg4FGUEKI0iIFJYQoLbPSSN7sBTty5EimbBdZ\nSN1CkmXfW4pscnKy51lVujnXo4itKh5l+/XNSpb3JHV/ZZY9n/v27WuRxfbX1fbf9RrZoIQQpaVs\nPzLlUpdCiGllaGgo+a9bQgjnhRC2pZyrEZQQok7RU7wQwl9SCXT5s6T2FNoaIcRA0YcR1I+BS1NP\n7ssIqtn4ODU1lclInjeYfZbYS3niNBU1f0+tv9d9ihl4IVvAudQ05f10BqS+J7E2pRrEvS0t4BvJ\nva0uDzzwgFv+5JN/35XnpVdevBDCRcCaJvEqM7srhLAs9Tqa4gkh6vTqR9bMbqES2DIXUlBCiDpa\nZiCEKC1lW2YgBSWEqNOPEZSZ3Q/cn3JukoIKIbwDuMbMloUQXg9spZIs4THgE2bW1uLtGcmzGLmz\nGMmLyBibSpbYQXnblGXVsndPvHOb66kF8U81aMcyoHgvvWeMjRloe/mlyZoIIubM8frvGcSzJE3w\nMgv/6le/cssXZSQv2wiq45MPIawDbgbmV0UbgavM7HQqebiWF9c8IUQ/mTNnTvJfP0j5afpn4I8b\nPp8C/LB6fC9wTq8bJYSYHvqxkjxTe1KGuiGEk4BvmNk7QwjPmNmJVflZVHLirWhX/plnnpk68cQT\ne9FeIUSc3Frjgx/8YLKNY/v27aXMLNw4KR8DXuhU4Itf/GL9ePPmzVx66aVRe8kg2qBuu+02Lrzw\nwky/Knn7VLQN6jvf+Q7Ll1dm74Nig+r07GrPKXZulsWjXlu9bMPz5s1zy8+fP79Fdswxx7TIXvOa\n17jl3/veyqSm15mFy2aD6kZB/TSEsKxqiT8X+PtOBbyMtbEbkXfVdJYvfiopL267L1Fqm2Krjj25\ntxLZy2wbOzfFSL5nzx4gPUFATEF5cu+L633BY+VTV6fHzkl9Jlne0yw/Gt4z9QzntWfQzGOP7QDg\nXe961yuO8zITFNSfATeFEEaBJ4HtvW2SEGK6GMiFmmb2L8A7q8f/GzizwDYJIaaJmTCCEkLMUAZy\nBCWEmB1oBCWEKC2zUkE1DxuHh4cL2+pShMcuxWPUrk+em9/z4sS8cLUtJ414sYOyePFiKbkb2b17\nd/R/Wbx4nkvdq3/BggVu+VSPX5bpSd7tR6mxo2KeWe/+ec/Pe/YAO3fudI/zMisVlBBiMJANSghR\nWjSCEkKUFikoIURpmZUKqtkgOGfOnNzGx7yxl4rI7Bvrk7fdwTOI7t+/3y3vxQnyjKfeVgno3kge\n22YB2YzkXrs8g3i7JA3NeM8ktlWm0a7Sqy9gatKFWJ9SEyx4zhCAXbt2ucd5mZUKSggxGEhBCSFK\ni7x4QojSohGUEKK0zEoF5RnJY8HoPeNrlizEHnmN7KmG91g7PYO4Z/jet2+fW94znnvG05iR3DO+\npsSD8tpYw5sKxIzknkHea1OWBAVe/bHpSWP5djsNUp9zTJ7lnfD6nxr3C175bNo9p6zMSgUlhBgM\npKCEEKVFRnIhRGnRCEoIUVoGUkE1ZRb+feB7wC+q/95sZncV1UAhRP8YOAVVzSx8AVBzJZ0CbDSz\na5MrafLujIyMRL0bRWx1yUJq+eZ2Tk5ORj0unsfN87zEvDGe3PPYZYkHlZKeK7bNArJ58by6PFnM\n/pG6rSYl7VWt3iJsLVm8eKnbYmLvVOPzj3lvu2HgFBQvZxa+vfr5FCCEEJZTGUVdbmZ7C2qfEKKP\nlE1BdfwZMbO7gUY1/jDwaTM7A/glcHVBbRNC9Jnh4eHkv37QTerzRWb2QlX+JmCTmZ3drvyzzz47\ntXTp0l60VwgRJ/fwZ926dckxs9evX1/K1OffDyFcZmYPA2cDP+lUYOPGjfXj9evXs27dutzxs7Ok\n+c5C6i9DY/1bt25l5cqVUVuAt0LcC2Wyd68/U/ZWkhdtg3r++ec59thj3etBNhvU6Ohokuzoo492\ny3vy8fHxFtnYmJ/6u1bX7bffzgUXXADkT6eeGsInVo8XGsaL3X7UUUe55RcvXgzAFVdcwTXXXFM/\nzkvZpnjdKKhLgU0hhJeAZ4FLOhVofhlHR0czxckp4qZluWZKmusjR45EA+R7xuZUGfjKqB9G8nbG\nV++LF/txSE1kEXMSeIrPU3CxpAuNRvbac8s7RcmSYMEjNU167J1qfNax594NA6mgmjILPwKcVmCb\nhBDTxEAqKCHE7EBbXYQQpUUjKCFEaZmVCsrLLFy2G9GOlBXC7VaSe0ZMTxYzSqeWjxlLUw2yzf1s\nl1jBW92dxUjsPX/vmpDuJIjdv0aPWe25pXrhsvQpb2bi1BhRzfLYOd1Qtu+lRlBCiDpSUEKI0lKk\nggohHAPcAYwDo8CnzOx/tCtTLpO9EGJamTNnTvJfF3wK+DszOxNYCfyXTgU0ghJC1Cl4incdUDMU\njgDxcBkNJwkhBNA7BRVCuAhY0yReZWb/GEJYSmWqd3mn60xbVpeUDBxlIdWLl+JxaSfLsv0nVRar\nK2WrS7tU5Fm8W5430Hv+Ma9h6v2L3f/GftSOU714MVL7n2UvX5bU6WXf6mJmtwC3NMtDCG8BvgGs\nNbMfdrqORlBCiDoFG8nfBHwT+JCZPZpSRgpKCFGn4BnMl4H5wF+GEAB2m9nydgWkoIQQdYrci9dJ\nGXlIQQkh6pTNBiwF1SWekTyvkTsWTymvkTw1aUGz4TZr8L9Y/amxo2L1pQYsjNXf2K/acd54Tp3q\n6URq/1PeiXbOjKwomoEQorRoBCWEKC1SUEKI0qIpnhCitAzUCCqEMBfYApwEzAO+BDwBbAWmgMeA\nT5hZW2tqc6ez3gRPqxeVWbhfZMlCm9egW4RBeLrJYlDOYyTPshK8iGeSkpk4byajRsr2/ek0nlsB\n/JuZnQ68B/hrYCNwVVU2BGRe2yCEKCdDQ0PJf/2gk4L6JvDZ6vEQcJhK6vPaHpp7gXOKaZoQot+U\nTUGlZhYeA74L3ARsMLMTq/KzgNVmtqJd+d/+9rdTS5Ys6UFzhRBtyK01brjhhuR56sc//vHpzywc\nQvgd4B7gBjPbFkJY3/DvMeCFTte48cYb68ef+9zn+MIXvuBm2wU/eWVqtmHIPx/3fhm8hXATExP1\n47vuuosPfehDbrZggF27drXIdu/e3SKLZRZurKuGF387SzSATvdpamqq7a+kZxeMBTHzsujOmzev\nReZl1gVYuHBhi8zLLFzLttvMokWLALj77rv5wAc+EK3L61OWmOJZIiR4daVmG4aXMw5v2LCBtWvX\n1o/z0mUgusJoO8ULIbwK+AFwhZltqYp/GkJYVj0+F3iwuOYJIfpJ2aZ4nUZQnwEWA58NIdRsUZ8E\n/iqEMAo8CWzvVEmWrC59m9tmqD9FlvWheedmiZGVt668GUyyUMQz9e5VP9OZy4vXH9oqKDP7JBWF\n1MyZxTRHCDGdaKGmEKK0DNQISggxu5CCEkKUllk5xWt2n86dO5eREb/qsmlwSDNoDw8PR120Xl+9\nc2Mvh3duqgy6D/Df7mXNsswg1aAdeydS+5ri+Mj6fhVlJM+Lt32nF8xKBSWEGAzKNkCQghJC1JGC\nEkKUFikoIURpKdtWl74oqLwryfNq9SxGxNS6vD7FjLzeHitPFivvyb39dbHyHt7+wub71O5lzWKk\n99rl9T9WfnR0NOmasf43Pqt2RuC8K8lnAhpBCSFKixSUEKK0SEEJIUqL1kEJIUqLRlBCiNIyK0dQ\nnscrS+yjIsjinUnd6uJ5psCPHpkqg2xp0j28e+15AZvxvGftrpnXixm7f57ca1usvUV48VLLl21E\n0olZqaCEEINB2RSqFJQQoo4UlBCitAyUgopkFv418D3gF9XTNpvZXQW2UQjRJwZKQfFyZuELQgjH\nAj8DvgBsNLNrUytp3sIwZ86cnictqJEaeD5WPvXcLFtdPOP3ggULWmRZUmllMch6RnZvW0nzNXtl\nJE81aMdSLHn3yrunMSO7ZyQvImlElvc09dx+G60HbS/eN3k5a0tjZuEQQlhOZRR1uZn5Cd2EEANF\n2UZQbdWzme0zs73VzMLbgauAh4FPm9kZwC+Bq4tvphCiH5QtL17H1OdNmYW3hBAWmdkL1f+9Cdhk\nZme3u8bOnTunjj/++F61WQjhk1tr3Hfffcnz3LPOOmt6U583ZBb+T2b2d1Xx90MIl5nZw8DZwE86\nVXLHHXfUjy+//HKuv/56Nx044KZE99KhHzp0yC2fmuY7r22g8Zpbt25l5cqV0Tbt37+/RealOY+l\nTvfKe+nQY/V7NijvPjX+WD3//PMce+yx7vWgGBuUl+IcYGxsrGsZvGyb2rZtG+effz7gt997T2IL\nYlPtVVnizGdJfV6TX3fddaxZs6Z+nJeyTfG6ySz8KeC6EMJLwLPAJZ0qaX4ZR0dHC8ui65FlhXCe\n+mMGRs+gm6o0oHvDfda6mu9J7Uvg1ZUlHpT3xfOUkWcMj53rfXH7GQ8qL3mdQXkSQWRt13TSbWbh\n04ppjhBiOilSQYUQjgK2URn0HAI+amZPtytTro03QohppWAj+Z8CP6k62O4A1nUqoJXkQog6RY6g\nzOz6EELNDvC7wAudykhBCSHq9EpBhRAuAtY0iVeZ2T+GEO4D3gK8u9N1pKCEEHV6paDM7Bbglsj/\nzgohvBH478Dr2l2nLwqq2Ys1b968TGmyixh25o0H1SybmpqKttPzLnleKM/bFqvfu38xL5a3hSYl\nxlTM7R+rP7Y1xvPieR67WH1HHXVUiyyLF6/x/tXer1TPaJYtVVk8q6np4FPK93I7TMFG8j8H/q+Z\n3Q7sA/wXvgGNoIQQdQpeZrAF+Fp1+jcHWNWpgBSUEKJOwUby3wDvyVJGCkoIUWegFmoKIWYXs1JB\nNRtEFyxYEI3d4xlfizCc9zoY/uTkZLRNqfuuPGOwV1fsmrGkC90aycfHx93rxeqPGclTYz912nfW\niGcQjxmL82xXKSqeU14jeVFbXZQ0QQhRWso2giqXuhRCiAY0ghJC1CnbCEoKSghRZ1YqqGZj6/j4\neCYjeaqRGNIy5rajW4Nq1nJZgtN7xnPP8JwliJ9nJG/uw6JFi4D8K9k9473X/lj5VIN4lt0BHlm+\nnKkG8dhzTnWcpNyT2DndMCsVlBBiMJAXTwhRWjSCEkKUFikoIURpGTgFVY2AdxMQgCngY8BBYGv1\n82PAJ8zMT38hhBBdkjKCeh+AmZ0WQlgG/AWV/FtXmdn9IYQbgeVUcue5HHPMMS2fY9syPI9E6vaX\nmLwfWWHabXVJrT/WJ8+7k+oFqrUtRdZM7bmlbsHI4rHKslUl9f71M3V56j3JkorLe36xZ9pYvl2K\n+qyUbQTV0WRvZt/m5dRSv0cljvApwA+rsnuBcwppnRCirwwPDyf/9YOOmYVrhBC+BpwHfBDYamYn\nVuVnAavNbEWs7P79+6diG2GFED0j9/DHzJKHoSGE6c0s3IiZfTSEcAXwP4HG8ARjdMjO8Mgjj9SP\nTz/9dB588EF+/vOfu+c+99xzLTIvC6+XbRjgxRdfbJGl7uaH7pI31jLW5p3iZVlo6E3RYn3qZop3\nzz33cN555wGDM8WL9akmv+2227jwwgvbnptKEVM8TxYzhdQiPHz+85/n6quvrh/nZeCmeCGEC6qx\nhAEOAJPAjqo9CuBc4MFimieE6CcF58XLTMoI6lvArSGEB4C5wOXAk8BNIYTR6vH2dheobZlo/Bz7\nZUg1FHqjIkj/ZYuNVlK2gHi0Oyc1wH6MvKnHuzUe1+J49TLAf6dzPVJHkFn6mTf1eOr9jxm5PXmW\nEVRj+VgdM4GOCsrM9gN/4vzrzN43RwgxnZRtq0u5WiOEEA1oJbkQok7ZjORSUEKIOrNSQS1delLL\n5wUL/GUGqStsY+5b79wsBtXUjLPdLEfolizG2zzlm2W1e5wa+yiv/SLLMotuDeK147x98t6/bleC\n18iSSKLxurFzuqFsCko2KCFEadEUTwhRR148IYRIRCMoIUSdstmgpk1BxTYPeytnsxgfUzLmQv4A\n+83nDQ8P5zbyxl6OIsKIFEHebM15z03ZC1mbwuRNBJF3Jbgnb87A3a58Y/tj53SDFJQQorSUTUHJ\nBiWEKC0aQQkh6pTNiycFJYSooymeEEIkMm0jqKOPPtqVpy73zxsRM29ExebyvfDiZSGvFzClfK2P\n3rA/dUtQForYvgOvbH/NU5ca/TLvVpWYh83z2C1cuLBFFovx1RgPrZfTsrKNoDTFE0LUkYISQpQW\nKSghRGkZOAUVySw8F/ge8IvqaZvN7K6iGimE6A8Dp6DwMwv/N2CjmV3bbcWxpAepxscsRvK8BnGP\n5muOjIxE60k1UufdflMUebfqFBHPKkvspkZDc83onZrgIGYk995Jz/DtyQAWL17cIvMcR7Hvyb59\n++rHvdzqUjZSkiZ8O4TwverHxszCIYSwnMoo6nIza01eJ4QYKPoxggohvJFKfs1XmZmf4LJKkn/S\nzA5XMwtvAu4EHgY+bWZnAL8Ers7XZCFEGSg6L14IYRy4FmjNsOu1J8v0IYSwlIrm+/dm9nRV9iZg\nk5mdHSt3+PCRqZERfz2HEKJn5B7+vPDCC8kKYdGiRZnqq6ZK/zrwZeA7wBs7jaBSjOQXAK8xsy/z\ncmbhb4UQLjOzh4GzgZ+0u8auXQfqxyecMMZzz+1lx45/cM/9zW9+0yLbs2dPi2z//v1u+YmJiRbZ\ngQMHWmSHDh1yy3up0w8fPtwia7Q33XrrraxatWrgbVCNv4p33HEHK1asaJGnlu8kz2KX6pUNavPm\nzVx66aWvkDUyiDaoiy++mJtvvrl+nJdeTfFCCBcBa5rE/wf4hpk9GkJIuk63mYV/DWwKIbwEPAtc\nktrwGjHDXl4juUfeL7j3Mjcrrblz50YVlCfPG6OqH0qrXYIBj5iCSFVGWTITZ8ms3CivvV+pq8az\nrAT3YpyNjY255U844YQWmbeSvNEY3kjjj27ZNvgCmNktwC2NshDCU8BFVeW1FPgBcEa76+TJLHxa\ncmuFEANBkUZyM3t97TiE8C/Af+hUpnyqVwghqmgluRCiTr8WaprZSSnnSUEJIeqUbSW5pnhCiNIy\nbSMoz1sXk3syz/WfhdgvhecJ8ly9zfXPnz8/6ln0PHadli40kneZQrfU7kXeZQael8k7N5ZBJTXb\nSsyL13huzSOc+p7F0op7cs8LF4t75nn3vKUHsexHjctsYud0g0ZQQgiRiGxQQog6ZRtBSUEJIeqU\nTUFpiieEKC3TNoKKGXQ942eWeDfedoVuDao1vH17zbJ2RnJP7tWTpXw/jOS1e5nXSO7d69SkBbHy\nqanHm+U147T3TqUm7IjV721/Oe6449zyb3vb21pkB9tum30lO3furB+Pj4+nF+yARlBCCJGIFJQQ\norTISC6EqFO2KZ4UlBCijhRUlTe84Q2u/PHHH2+ReatxYwHrPIO2F8QuJXZQjRTD+8KFC6Or2z25\nJyubkbxmSM4bcC41nlPMSJ531Xej8bq2WtszaHvlY+9J6r3eu9cP1f/oo4+2yLz3PPZONQZW9IIs\nzhQ0ghJC1CnbCEpGciFEadEISghRRyMoIYRIJGkEFUJYQiVzy7uBw8BWKmnQHwM+YWa9T90rhOg7\nZRtBpaSdmgv8V6DmCtsIXGVm94cQbgSWA/dkrTjmsUrdlhLz+Hhek9SsKjFStmosWLAgmiLI62uq\nDHxPTj9SvMe8YpA/K0uWFE/ethTPi+fFY4JXxktatGhRtLzXp9gzzftMvHty0NnrEnsGjeV7mdWl\nbAoqpWcbgBuBZ6qfTwF+WD2+FzingHYJIUT7zMIhhJVUknZ+KYRwP/Ax4D4zO7H6/7OA1Wa2ol0l\nyiwsRF/IPfyZyrCYbqgPw61OU7zVwFQI4RzgbcBtwJKG/48BL3SqxMss/NJL/gK2p556qkWWGoYX\n/Cmel1nYW7wZO9db/Nm4OO6KK67gmmuumVFTvA0bNrB27dro/wdxiveRj3yEO++8M1q+iClerE9e\nyN8s0RSef/55AN7//vfz3e9+t34802iroMysnvWzYQT1n0MIy8zsfuBc4O+LbKAQon+UzQbVdorX\nSIOCmgRuAkaBJ4E/NbP0PORCCJFIsoISQoh+o4WaQojSIgUlhCgtUlBCiNIiBSWEKC1SUEKI0tK3\ncCshhGHgBuCtwIvAxWbWuipzQAghvAO4xsyWhRBezwBvoK7ut9wCnATMA74EPMFg92kOleUwgUof\nPgYcZID7VGM2bd7v5wjqj4D5ZvYHwJXAtX2su6eEENYBNwO1Zb61DdSnU9lusHy62tYlK4B/q7b/\nPcBfM/h9eh+AmZ0GXAX8BYPfp3ab9we2T+3op4L6Q+BvAczsIeDUPtbda/4Z+OOGz4O+gfqbwGer\nx0NUfpUHuk9m9m3gkurH36OyJWug+1RlVm3e76eCGgd2N3w+EkIYyIieZnY30LhJa8jMaite9wLH\n9L9V3WNm+8xsbwhhDNhOZcQx0H0CMLPDIYSvAZuAOxnwPlU37z9nZt9vEA90nzrRTwW1h8rm4nrd\nZuanrBg8Guf8SRuoy0YI4Xeo7Ku83cy2MQP6BGBmHwXeQMUe1ZjKZRD7tBp4d3XbWdeb9weJfiqo\nHwHvBQghvBP4pz7WXTQ/DSEsqx6fCzw4jW3JTAjhVcAPgCvMbEtVPOh9uiCE8OfVjweoKNwdg9wn\nMzvDzM48kq+zAAAAi0lEQVQ0s2XAz4ALgXsHuU+d6OcU6x4q2v/HVOwcq/pYd9H8GXBTCKG2gXr7\nNLcnK58BFgOfDSHUbFGfBP5qgPv0LeDWEMIDwFzgcir9GOTn5DHo715btFlYCFFatFBTCFFapKCE\nEKVFCkoIUVqkoIQQpUUKSghRWqSghBClRQpKCFFapKCEEKXl/wPgKvtBO8twwQAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x259ca68d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab as plt\n",
    "plt.imshow(np.log10(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1]/np.max(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1])))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_MIPS=xidplus.prior(MIPS_Map[1].data,MIPS_Map[2].data,MIPS_Map[0].header,MIPS_Map[1].header,moc=Sel_func)\n",
    "prior_MIPS.prior_cat(masterlist[1].data['ra'][good],masterlist[1].data['dec'][good],'master_catalogue_lockman-swire_20170817.fits',flux_lower=MIPS_lower,\n",
    "                     flux_upper=MIPS_upper,ID=masterlist[1].data['help_id'][good])\n",
    "prior_MIPS.set_prf(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1]/1.0E6,np.arange(0,41/2.0,0.5),np.arange(0,41/2.0,0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate tiles\n",
    "As fitting the whole map would be too computationally expensive, I split based on HEALPix pixels. For MIPS, the optimum order is 11. So that I don't have to read the master prior based on the whole map into memory each time (which requires a lot more memory) I also create another layer of HEALPix pixels based at the lower order of 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 14127 tiles required for input catalogue and 73 large tiles\n"
     ]
    },
    {
     "ename": "OSError",
     "evalue": "[Errno 22] Invalid argument",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mOSError\u001b[0m                                   Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-27-4f593b9f51c1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0moutfile\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0moutput_folder\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'Master_prior.pkl'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutfile\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'wb'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m     \u001b[0mpickle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdump\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m'priors'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mprior_MIPS\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'tiles'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mtiles\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'order'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0morder\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'version'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mxidplus\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgit_version\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0moutfile\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0moutput_folder\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'Tiles.pkl'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutfile\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'wb'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mOSError\u001b[0m: [Errno 22] Invalid argument"
     ],
     "output_type": "error"
    }
   ],
   "source": [
    "import pickle\n",
    "#from moc, get healpix pixels at a given order\n",
    "from xidplus import moc_routines\n",
    "order=11\n",
    "tiles=moc_routines.get_HEALPix_pixels(order,prior_MIPS.sra,prior_MIPS.sdec,unique=True)\n",
    "order_large=7\n",
    "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior_MIPS.sra,prior_MIPS.sdec,unique=True)\n",
    "print('----- There are '+str(len(tiles))+' tiles required for input catalogue and '+str(len(tiles_large))+' large tiles')\n",
    "output_folder=''\n",
    "outfile=output_folder+'Master_prior.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'priors':[prior_MIPS],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)\n",
    "outfile=output_folder+'Tiles.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'tiles':tiles,'order':order,'tiles_large':tiles_large,'order_large':order_large,'version':xidplus.io.git_version()},f)\n",
    "raise SystemExit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('Master_prior.pkl', 'wb') as f:\n",
    "    pickle.dump({'priors':prior_MIPS,'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1e2d0dc (Sun Oct 8 14:24:36 2017 +0100) [with local modifications]'"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xidplus.io.git_version()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Large pickle problem\n",
    "It would appear that if the file trying to be saved is too big, then pickle will fail. To get round this, use the following fix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class MacOSFile(object):\n",
    "\n",
    "    def __init__(self, f):\n",
    "        self.f = f\n",
    "\n",
    "    def __getattr__(self, item):\n",
    "        return getattr(self.f, item)\n",
    "\n",
    "    def read(self, n):\n",
    "        # print(\"reading total_bytes=%s\" % n, flush=True)\n",
    "        if n >= (1 << 31):\n",
    "            buffer = bytearray(n)\n",
    "            idx = 0\n",
    "            while idx < n:\n",
    "                batch_size = min(n - idx, 1 << 31 - 1)\n",
    "                # print(\"reading bytes [%s,%s)...\" % (idx, idx + batch_size), end=\"\", flush=True)\n",
    "                buffer[idx:idx + batch_size] = self.f.read(batch_size)\n",
    "                # print(\"done.\", flush=True)\n",
    "                idx += batch_size\n",
    "            return buffer\n",
    "        return self.f.read(n)\n",
    "\n",
    "    def write(self, buffer):\n",
    "        n = len(buffer)\n",
    "        print(\"writing total_bytes=%s...\" % n, flush=True)\n",
    "        idx = 0\n",
    "        while idx < n:\n",
    "            batch_size = min(n - idx, 1 << 31 - 1)\n",
    "            print(\"writing bytes [%s, %s)... \" % (idx, idx + batch_size), end=\"\", flush=True)\n",
    "            self.f.write(buffer[idx:idx + batch_size])\n",
    "            print(\"done.\", flush=True)\n",
    "            idx += batch_size\n",
    "\n",
    "\n",
    "def pickle_dump(obj, file_path):\n",
    "    with open(file_path, \"wb\") as f:\n",
    "        return pickle.dump(obj, MacOSFile(f), protocol=pickle.HIGHEST_PROTOCOL)\n",
    "\n",
    "\n",
    "def pickle_load(file_path):\n",
    "    with open(file_path, \"rb\") as f:\n",
    "        return pickle.load(MacOSFile(f))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "writing total_bytes=2555266188...\n",
      "writing bytes [0, 1073741824)... done.\n",
      "writing bytes [1073741824, 2147483648)... done.\n",
      "writing bytes [2147483648, 2555266188)... done.\n"
     ]
    }
   ],
   "source": [
    "pickle_dump({'priors':[prior_MIPS],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},'Master_prior.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:new]",
   "language": "python",
   "name": "conda-env-new-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}