fiberassign/FiberAssign.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fiber Assignment\n",
    "\n",
    "This tutorial will teach you:\n",
    "* The required inputs to fiber assignment.\n",
    "* How to run `fiberassign` on those.\n",
    "* How to interpret the fiber assignment output files, including understanding coverage.\n",
    "\n",
    "This tutorial covers the high level `fiberassign` script.  See https://fiberassign.readthedocs.io/en/latest/ for installation instructions and descriptions of the lower level `fba_run`, `fba_merge`, `fba_run_qa` etc. scripts.\n",
    "\n",
    "In this notebook the inputs to fiberassign are random targets. The notebooks `FiberAssignDECaLS.ipynb` and `FiberAssignMocks.ipnynb` show how \n",
    "to use inputs from observations and mock files, respectively.\n",
    "\n",
    "See https://desi.lbl.gov/trac/wiki/Computing/JupyterAtNERSC for instructions on configuring jupyter kernels with pre-installed DESI software at NERSC.  This tutorial was last tested with the 22.2 kernel on April 29 2022\n",
    "\n",
    "Stephen Bailey  \n",
    "December 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic python imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys, subprocess\n",
    "from collections import Counter\n",
    "from datetime import datetime\n",
    "\n",
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "from astropy.io import fits\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import desimodel.io\n",
    "import desimodel.focalplane\n",
    "import desimodel.footprint\n",
    "from desitarget.targetmask import desi_mask, obsconditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an output directory in $SCRATCH/desi/test/fiberassign"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can set this to a local path if you have the software installed and a copy of the data.\n",
    "# workdir = os.path.join(os.environ['HOME'], 'scratch', 'desi', 'tutorials', 'fiberassign_obs')\n",
    "workdir = os.path.join(os.environ['SCRATCH'], 'desi', 'test', 'fiberassign')\n",
    "os.makedirs(workdir, exist_ok=True)\n",
    "os.chdir(workdir)\n",
    "outdir = os.path.join(workdir, 'output')\n",
    "os.makedirs(outdir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fiber Assignment inputs\n",
    "\n",
    "Fiber assignment requires the following input files:\n",
    "\n",
    "**TODO**: document what each column is\n",
    "\n",
    "**TODO**: some of these are no longer required in fiberassign 1.0.0; identify and remove those (likely BGS_TARGET, MWS_TARGET, and BRICKNAME).\n",
    "    \n",
    "* \"Merged Target List\" (MTL): FITS file with columns:\n",
    "  * TARGETID\n",
    "  * RA\n",
    "  * DEC\n",
    "  * DESI_TARGET\n",
    "  * BGS_TARGET\n",
    "  * MWS_TARGET\n",
    "  * PRIORITY\n",
    "  * SUBPRIORITY\n",
    "  * BRICKNAME\n",
    "  * OBSCONDITIONS\n",
    "  * NUMOBS_MORE\n",
    "* sky locations; FITS file with columns:\n",
    "  * TARGETID\n",
    "  * RA\n",
    "  * DEC\n",
    "  * DESI_TARGET\n",
    "  * BGS_TARGET\n",
    "  * MWS_TARGET\n",
    "  * PRIORITY\n",
    "  * SUBPRIORITY\n",
    "  * BRICKNAME\n",
    "  * OBSCONDITIONS\n",
    "\n",
    "Standard stars can either be in the same file as the science targets, or they can\n",
    "be in a separate file (with the same required columns).  `fiberassign` uses `DESI_TARGET`\n",
    "mask bits `STD_FAINT`, `STD_BRIGHT`, and `STD_WD` to determine which targets may be\n",
    "used as standard stars.\n",
    "\n",
    "Optional inputs:\n",
    "\n",
    "* footprint tiles (default `$DESIMODEL/data/footprint/desi-tiles.fits`) with columns\n",
    "  * TILEID\n",
    "  * RA\n",
    "  * DEC\n",
    "  * PASS\n",
    "  * OBSCONDITIONS\n",
    "  * IN_DESI\n",
    "  * PROGRAM\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fiber Assignment from scratch\n",
    "\n",
    "First we'll run fiberassignment on a set of random targets and tiles to understand the required inputs.  Then we'll proceed with running it on real data for real DESI tile locations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom tiling\n",
    "\n",
    "We'll start by generting a set of offset tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ntiles = 5\n",
    "tiles = Table()\n",
    "tile_radius = desimodel.focalplane.get_tile_radius_deg()\n",
    "tiles['TILEID'] = np.arange(ntiles, dtype='i4')\n",
    "tiles['RA'] = 2 + np.arange(ntiles)*0.5*tile_radius\n",
    "tiles['DEC'] = np.zeros(ntiles)\n",
    "tiles['PASS'] = np.zeros(ntiles, dtype='i2')\n",
    "tiles['OBSCONDITIONS'] = np.ones(ntiles, dtype='i4') * obsconditions.DARK\n",
    "tiles['IN_DESI'] = np.ones(ntiles, dtype='i2')\n",
    "tiles['PROGRAM'] = np.full(ntiles, 'DARK', dtype='S6')\n",
    "tiles.write('tiles.fits', format='fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sky_circle(x, y, radius, alpha=0.2):\n",
    "    theta = np.linspace(0, 2*np.pi, 50)\n",
    "    xx = x + np.cos(np.radians(y))*radius*np.cos(theta)\n",
    "    yy = y + radius*np.sin(theta)\n",
    "    plt.fill(xx, yy, alpha=alpha)\n",
    "\n",
    "def plot_tile(ra, dec, alpha=0.2):\n",
    "    tile_radius = desimodel.focalplane.get_tile_radius_deg()\n",
    "    plot_sky_circle(ra, dec, tile_radius)\n",
    "\n",
    "def plot_tiles(tiles):\n",
    "    for i in range(len(tiles)):\n",
    "        plot_tile(tiles['RA'][i], tiles['DEC'][i])\n",
    "    plt.xlabel('RA [degrees]')\n",
    "    plt.ylabel('dec [degrees]')\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plot_tiles(tiles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare that tiling coverage to the fiber coverage.  We'll see that the coverage of\n",
    "actual fibers (positioners) is less than that of the tile, and has a more complex spatial structure.  By default, desimodel returns the realistic state of the focalplane at the time you run this notebook.  That is great for working in realtime on real data.  For this exercise, we will be using a focalplane from a time prior to commissioning so that we will have all positioners working and with their nominal properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rundatestr = \"2022-01-01T00:00:00+00:00\"\n",
    "rundate = datetime.strptime(rundatestr, \"%Y-%m-%dT%H:%M:%S+00:00\")\n",
    "print(rundate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_positioners(tilera, tiledec, alpha=0.1):\n",
    "    fp, exclude, state, tmstr = desimodel.io.load_focalplane(time=rundate)\n",
    "    # Select just science positioners\n",
    "    rows = np.where(fp['DEVICE_TYPE'] == 'POS')[0]\n",
    "    fp = fp[rows]\n",
    "    ra, dec = desimodel.focalplane.xy2radec(tilera, tiledec, fp['OFFSET_X'], fp['OFFSET_Y'])\n",
    "    plt.plot(ra, dec, '.', alpha=alpha)\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "for i in range(len(tiles)):\n",
    "    plot_positioners(tiles['RA'][i], tiles['DEC'][i])\n",
    "\n",
    "plt.xlabel('RA [degrees]')\n",
    "plt.ylabel('dec [degrees]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random targets\n",
    "\n",
    "Now we'll generate randomly distributed targets.  To facilitate counting, we'll trim those\n",
    "to just targets that are potentially covered by a tile (but might still not be covered by\n",
    "any fibers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_random_targets(density, tiles):\n",
    "    #- Get basic bounds; don't worry about RA wraparound for this example\n",
    "    tile_radius = desimodel.focalplane.get_tile_radius_deg()\n",
    "    ramin = np.min(tiles['RA'] - tile_radius*np.cos(np.radians(tiles['DEC'])))\n",
    "    ramax = np.max(tiles['RA'] + tile_radius*np.cos(np.radians(tiles['DEC'])))\n",
    "    decmin = np.min(tiles['DEC']) - tile_radius\n",
    "    decmax = np.max(tiles['DEC']) + tile_radius\n",
    "    \n",
    "    area = (ramax-ramin) * np.degrees((np.sin(np.radians(decmax)) - np.sin(np.radians(decmin))))\n",
    "    n = int(area*density)\n",
    "\n",
    "    #- Iterate if needed to get unique TARGETIDs\n",
    "    while True:\n",
    "        targetids = np.random.randint(0, 2**62-1, n)\n",
    "        if len(set(targetids)) == n:\n",
    "            break\n",
    "\n",
    "    #- Create targets table\n",
    "    targets = Table()\n",
    "    targets['TARGETID'] = targetids\n",
    "    targets['RA'] = np.random.uniform(ramin, ramax, n)\n",
    "    phimin = np.radians(90-decmin)\n",
    "    phimax = np.radians(90-decmax)\n",
    "    targets['DEC'] = 90-np.degrees(np.arccos(np.random.uniform(np.cos(phimin), np.cos(phimax), n)))\n",
    "    targets['DESI_TARGET'] = np.zeros(n, dtype='i8')\n",
    "    targets['BGS_TARGET'] = np.zeros(n, dtype='i8')\n",
    "    targets['MWS_TARGET'] = np.zeros(n, dtype='i8')\n",
    "    targets['SUBPRIORITY'] = np.random.uniform(0, 1, n)\n",
    "    targets['BRICKNAME'] = np.full(n, '000p0000')    #- required !?!\n",
    "    targets['BRICKID'] = np.full(n, 0)    #- required !?!\n",
    "    targets['BRICK_OBJID'] = np.arange(n)\n",
    "    \n",
    "    #- dummy values for fluxes\n",
    "    for filt in ['G', 'R', 'Z']:\n",
    "        targets['FIBERFLUX_'+filt] = np.zeros(n, dtype='f4')\n",
    "        targets['FIBERFLUX_IVAR_'+filt] = np.ones(n, dtype='f4')\n",
    "    \n",
    "    #- Trim to targets that are covered by a tile\n",
    "    ii = desimodel.footprint.is_point_in_desi(tiles, targets['RA'], targets['DEC'])\n",
    "    targets = targets[ii]\n",
    "    \n",
    "    return targets\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- ELG-like sample: high density, only one requested observation per target\n",
    "density = 2400\n",
    "targets = generate_random_targets(density, tiles)\n",
    "n = len(targets)\n",
    "targets['PRIORITY'] = 1000\n",
    "targets['SUBPRIORITY'] = np.random.uniform(0, 1, n)\n",
    "targets['DESI_TARGET'] = desi_mask.ELG\n",
    "targets['OBSCONDITIONS'] = np.ones(n, dtype='i4') * obsconditions.DARK\n",
    "targets['NUMOBS_MORE'] = np.ones(n, dtype='i8')\n",
    "targets.meta['EXTNAME'] = 'MTL'\n",
    "targets.write('mtl.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Sky targets at 4x the density of fibers\n",
    "fiber_density = 5000 / 7.5\n",
    "sky_density = 4*fiber_density\n",
    "sky = generate_random_targets(sky_density, tiles)\n",
    "nsky = len(sky)\n",
    "sky['DESI_TARGET'] = desi_mask.SKY\n",
    "sky['OBSCONDITIONS'] = np.ones(nsky, dtype='i4') * obsconditions.mask('DARK|GRAY|BRIGHT')\n",
    "sky.meta['EXTNAME'] = 'SKY'\n",
    "sky.write('sky.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_density = 50\n",
    "stdstars = generate_random_targets(std_density, tiles)\n",
    "nstd = len(stdstars)\n",
    "stdstars['SUBPRIORITY'] = np.random.uniform(0, 1, nstd)\n",
    "stdstars['DESI_TARGET'] = desi_mask.STD_FAINT\n",
    "stdstars['OBSCONDITIONS'] = np.ones(nstd, dtype='i4') * obsconditions.DARK\n",
    "stdstars.meta['EXTNAME'] = 'STD'\n",
    "stdstars.write('stdstars.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(targets['RA'], targets['DEC'], 'k,', alpha=0.3, label='targets')\n",
    "plt.plot(stdstars['RA'], stdstars['DEC'], 'rx', alpha=0.5, label='stdstars')\n",
    "plot_tiles(tiles)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run fiberassign on those targets/standards/sky\n",
    "\n",
    "Fiberassign calls desimodel to get focalplane properties.  By default this uses the current time when you run this notebook.  Instead, we will use the same historical time above in order to get a nominal focalplane prior to commissioning will all positioners working."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd = 'fiberassign --overwrite --mtl mtl.fits --stdstar stdstars.fits --sky sky.fits'\n",
    "cmd += ' --rundate {}'.format(rundatestr)\n",
    "cmd += ' --footprint ./tiles.fits'\n",
    "cmd += ' --outdir ./output/'\n",
    "cmd = cmd.format(outdir=outdir)\n",
    "\n",
    "print('RUNNING: '+cmd)\n",
    "try:\n",
    "    results = subprocess.check_output(cmd.split(), stderr=subprocess.STDOUT)\n",
    "    print(results.decode())\n",
    "    print('--- SUCCESS ---')\n",
    "except subprocess.CalledProcessError as ex:\n",
    "    print('--- ERROR {} ---'.format(ex.returncode))\n",
    "    print(ex.output.decode())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fiberassign output\n",
    "\n",
    "`fiberassign` output files contain the following HDUs:\n",
    "* `FIBERASSIGN`: fiber assignments for the 5000 positioners with fibers connected to a spectrograph\n",
    "* `POTENTIAL_ASSIGNMENTS`: a table of what targets were covered by what fibers, whether or not they were assigned\n",
    "* `SKY_MONITOR`: assignments of the 20 fibers that go to the sky monitor camera (not the spectrographs)\n",
    "* `TARGETS`: (deprecated) a copy of the input targets table with all columns\n",
    "\n",
    "The latest version of fiberassign also has HDUS `FASSIGN`, `FTARGETS`, and `FAVAIL`, which we won't cover here (yet).\n",
    "\n",
    "Note: previously the primary fiberassign output files were called `tile-*.fits`, which caused confusion\n",
    "with the \"tile file\" defining the DESI footprint, so these have been renamed `fiberassign-*.fits`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fafile = '{}/fiberassign-{:06d}.fits'.format(outdir, tiles['TILEID'][0])\n",
    "fx = fits.open(fafile)\n",
    "print(fx.info())\n",
    "print('\\nFIBERASSIGN columns:')\n",
    "print(fx['FIBERASSIGN'].data.dtype.names)\n",
    "print('\\nPOTENTIAL_ASSIGNMENTS columns:')\n",
    "print(fx['POTENTIAL_ASSIGNMENTS'].data.dtype.names)\n",
    "fx.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Which targets were assigned?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assignments = list()\n",
    "for tileid in tiles['TILEID']:\n",
    "    fafile = '{}/fiberassign-{:06d}.fits'.format(outdir, tileid)\n",
    "    assignments.append(Table.read(fafile, 'FIBERASSIGN'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assigned_targetids = np.concatenate([tmp['TARGETID'] for tmp in assignments])\n",
    "isAssigned = np.isin(targets['TARGETID'], assigned_targetids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(targets['RA'][isAssigned], targets['DEC'][isAssigned], 'k,')\n",
    "plt.title('Targets assigned to fibers')\n",
    "plot_tiles(tiles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## And which targets weren't assigned?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(targets['RA'][~isAssigned], targets['DEC'][~isAssigned], 'k,')\n",
    "plt.title('Targets not assigned to fibers')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the number that were assigned to the simple estimate of counting fibers.\n",
    "For this density of targets, nearly all fibers are assigned to a target.\n",
    "i.e. if a target wasn't assigned, it was because either it wasn't covered by a fiber at all,\n",
    "or because that fiber was given to a different target instead.  Good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_assigned = np.count_nonzero(isAssigned)\n",
    "max_possible_assigned = int(len(tiles) * 5000 * 0.9)  #- 10% of fibers are reserved for std and sky calibrators\n",
    "print('Targets assigned        ', num_assigned)\n",
    "print('Max possible assignments', max_possible_assigned)\n",
    "print('Ratio                   ', num_assigned/max_possible_assigned)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assigning a target more than once - part I\n",
    "\n",
    "Now let's make a sample like Survey Validation, with a high target density but requesting 4 assignments per target when possible (`NUMOBS_MORE = 4`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Science Verification sample: high density, with multiple requested observations per target\n",
    "density = 2400\n",
    "targets = generate_random_targets(density, tiles)\n",
    "n = len(targets)\n",
    "targets['PRIORITY'] = 1000\n",
    "targets['SUBPRIORITY'] = np.random.uniform(0, 1, n)\n",
    "targets['DESI_TARGET'] = desi_mask.ELG\n",
    "targets['OBSCONDITIONS'] = np.ones(n, dtype='i4') * obsconditions.DARK\n",
    "targets['NUMOBS_MORE'] = np.ones(n, dtype='i8') * 4  #- THIS IS WHAT IS DIFFERENT THAN BEFORE\n",
    "targets.meta['EXTNAME'] = 'MTL'\n",
    "targets.write('mtl.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Rerun fiberassign\n",
    "#- Rerun fiberassign\n",
    "print('RUNNING: '+cmd)\n",
    "try:\n",
    "    results = subprocess.check_output(cmd.split(), stderr=subprocess.STDOUT)\n",
    "    print(results.decode())\n",
    "    print('--- SUCCESS ---')\n",
    "except subprocess.CalledProcessError as ex:\n",
    "    print('--- ERROR {} ---'.format(ex.returncode))\n",
    "    print(ex.output.decode())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Histogram how many times each target was assigned"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_assigned_covered(tiles):\n",
    "    assigned_ids = list()\n",
    "    covered_ids = list()\n",
    "    for tileid in tiles['TILEID']:\n",
    "        fafile = '{}/fiberassign-{:06d}.fits'.format(outdir, tileid)\n",
    "        assignments = Table.read(fafile, 'FIBERASSIGN')\n",
    "        assigned_ids.extend(assignments['TARGETID'])\n",
    "\n",
    "        #- Note: a single target could be covered by more than one fiber on the same\n",
    "        #- tile, but it still can only be assigned once, so use np.unique() per tile\n",
    "        #- to not double count those cases\n",
    "        potential = Table.read(fafile, 'POTENTIAL_ASSIGNMENTS')\n",
    "        covered_ids.extend(np.unique(potential['TARGETID']))\n",
    "\n",
    "    #- Trim to just science targets (i.e. not SKY, not non-science STDSTAR)\n",
    "    assigned_ids = np.array(assigned_ids)\n",
    "    ii = np.in1d(assigned_ids, targets['TARGETID'])\n",
    "    assigned_ids = assigned_ids[ii]\n",
    "\n",
    "    covered_ids = np.array(covered_ids)\n",
    "    ii = np.in1d(covered_ids, targets['TARGETID'])\n",
    "    covered_ids = covered_ids[ii]\n",
    "\n",
    "    return assigned_ids, covered_ids\n",
    "\n",
    "#- Do a bit of counting magic to include targets that were never assigned or covered\n",
    "def count_coverage(ids):\n",
    "    c = Counter(ids)\n",
    "    c.update(targets['TARGETID'])\n",
    "    count = np.array(list(c.values())) - 1\n",
    "    return count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assigned_ids, covered_ids = get_assigned_covered(tiles)\n",
    "assigned_count = count_coverage(assigned_ids)\n",
    "covered_count = count_coverage(covered_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(211)\n",
    "plt.hist(assigned_count, 11, (-0.5, 10.5), rwidth=0.8)\n",
    "plt.xticks(np.arange(0,11))\n",
    "plt.xlabel('number of times assigned')\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.hist(covered_count, 11, (-0.5, 10.5), rwidth=0.8)\n",
    "plt.xticks(np.arange(0,11))\n",
    "plt.xlabel('number of times covered')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why are the number of assignments so much less than the coverage?\n",
    "The current fiberassign algorithm prioritizes observing more targets before\n",
    "getting more observations of a previously observed target.  This is what we\n",
    "want for Lyman-alpha QSOs during the main survey, but it isn't what we want\n",
    "for Survey Validation (or LRGs).  GitHub [fiberassign ticket #140](https://github.com/desihub/fiberassign/issues/140) is tracking this.\n",
    "\n",
    "**TODO** Add an option to fiberassign to give the flexibility to finish targets\n",
    "before starting new ones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assigning a target more than once - part II\n",
    "\n",
    "Let's repeat the `NUMOBS_MORE>1` exercise with a low density of science targets so that\n",
    "they don't compete with each other so much."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- QSO-like sample: low density, multiple observations per target\n",
    "density = 50  #- \n",
    "targets = generate_random_targets(density, tiles)\n",
    "n = len(targets)\n",
    "targets['PRIORITY'] = 1000\n",
    "targets['SUBPRIORITY'] = np.random.uniform(0, 1, n)\n",
    "targets['DESI_TARGET'] = desi_mask.ELG\n",
    "targets['OBSCONDITIONS'] = np.ones(n, dtype='i4') * obsconditions.DARK\n",
    "targets['NUMOBS_MORE'] = np.ones(n, dtype='i8') * 4  #- THIS IS WHAT IS DIFFERENT THAN BEFORE\n",
    "targets.meta['EXTNAME'] = 'MTL'\n",
    "targets.write('mtl.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Rerun fiberassign\n",
    "print('RUNNING: '+cmd)\n",
    "try:\n",
    "    results = subprocess.check_output(cmd.split(), stderr=subprocess.STDOUT)\n",
    "    print(results.decode())\n",
    "    print('--- SUCCESS ---')\n",
    "except subprocess.CalledProcessError as ex:\n",
    "    print('--- ERROR {} ---'.format(ex.returncode))\n",
    "    print(ex.output.decode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assigned_ids, covered_ids = get_assigned_covered(tiles)\n",
    "assigned_count = count_coverage(assigned_ids)\n",
    "covered_count = count_coverage(covered_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(211)\n",
    "plt.hist(assigned_count, 11, (-0.5, 10.5), rwidth=0.8)\n",
    "plt.xticks(np.arange(0,11))\n",
    "plt.xlabel('number of times assigned')\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.hist(covered_count, 11, (-0.5, 10.5), rwidth=0.8)\n",
    "plt.xticks(np.arange(0,11))\n",
    "plt.xlabel('number of times covered')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks much better.  Let's confirm that the targets with low assignments are the ones with low coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n2d = np.zeros((10, 10), dtype=int)\n",
    "\n",
    "count_assigned = Counter(assigned_ids)\n",
    "count_covered = Counter(covered_ids)\n",
    "\n",
    "for tid in targets['TARGETID']:\n",
    "    n = count_assigned[tid]\n",
    "    m = count_covered[tid]\n",
    "    n2d[n,m] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(n2d)\n",
    "plt.xlabel('Number of times Covered')\n",
    "plt.ylabel('Number of times Assigned')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, nearly all cases of low assignments are due to low coverage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "1. Run an intermediate density case where the number of targets times the number of requested observations is approximately equal to the total number of fibers.\n",
    "2. Track down a case where a fiber received fewer than its maximum number of possible assignments based upon coverage.  Determine if the \"missing\" case was legitimate (e.g. if the fiber was assigned to a higher SUBPRIORITY target instead).\n",
    "3. Generate an input target list with two different kinds of targets at the same density but different `PRIORITY`.  Run `fiberassign` on those and study the outputs.  Do the higher priority targets get more assigned fibers?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO\n",
    "* Add examples for running on real data, including using `desitarget.mtl.make_mtl` to convert targets -> mtl.\n",
    "* When running fiberassign, capture stderr/stdout to a log and print an informative error message\n",
    "  if it crashes.\n",
    "* Add tutorials about lower-level `fba_run` and `fba_merge` commands."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "DESI 22.2",
   "language": "python",
   "name": "desi-22.2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}