fiberassign/FiberAssignAlgorithms_Part1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assignment Algorithms : Part 1\n",
    "\n",
    "The purpose of this notebook is to explore the fiberassign algorithms and the consequences of introducing different target populations and realistic nominal positioner motions and exclusion zones.  Very often people see behavior of the fiberassign code which does not match their intuition and expectations.  Usually these tests also use real target data, which further complicates the situation and makes it impossible to determine what effects are due to the target data which effects are just due to geometry.\n",
    "\n",
    "The goal of this notebook is to examine the impacts of the just the DESI geometry applied to synthetic targets with different densities and priorities.\n",
    "\n",
    "For all these cells, we are using a single RA / DEC location on the sky and one or more tiles.  To give a more representative spread on our results we do 100 realizations of the simulation for each exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports and Definitions\n",
    "\n",
    "These are global imports and variables used throughout the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "from collections import OrderedDict\n",
    "\n",
    "import subprocess as sp\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from IPython.display import Image, display\n",
    "\n",
    "import fitsio\n",
    "\n",
    "from desitarget.targetmask import desi_mask\n",
    "\n",
    "from fiberassign import __version__ as fba_version\n",
    "\n",
    "from fiberassign.hardware import (\n",
    "    load_hardware,\n",
    ")\n",
    "\n",
    "from fiberassign.targets import (\n",
    "    Targets,\n",
    "    TargetTree,\n",
    "    TargetsAvailable,\n",
    "    LocationsAvailable,\n",
    "    load_target_table, \n",
    "    TARGET_TYPE_SCIENCE, \n",
    "    TARGET_TYPE_SKY,\n",
    "    TARGET_TYPE_SUPPSKY,\n",
    "    TARGET_TYPE_STANDARD\n",
    ")\n",
    "\n",
    "from fiberassign.tiles import (\n",
    "    load_tiles,\n",
    ")\n",
    "\n",
    "from fiberassign.assign import (\n",
    "    Assignment,\n",
    ")\n",
    "\n",
    "from fiberassign.vis import (\n",
    "    plot_assignment_tile,\n",
    ")\n",
    "\n",
    "# Date used for the focalplane model\n",
    "assign_date = \"2020-01-01T00:00:00\"\n",
    "\n",
    "# This is a small patch around one DESI pointing\n",
    "patch_ra_min = 148.0\n",
    "patch_ra_max = 152.0\n",
    "patch_dec_min = 29.0\n",
    "patch_dec_max = 33.0\n",
    "\n",
    "# Tile center (for our co-incident tiles)\n",
    "tile_ra = 150.0\n",
    "tile_dec = 31.0\n",
    "\n",
    "# Target densities (table 3.1 of DSR)\n",
    "target_density = {\n",
    "    \"ELG\": 2400,\n",
    "    \"LRG\": 350,\n",
    "    \"QSO-tracer\": 170,\n",
    "    \"QSO-lyman\": 90,\n",
    "    \"FAKE-high\": 4000,\n",
    "    \"standards\": 300,\n",
    "    \"sky\": 4000\n",
    "}\n",
    "\n",
    "# DESITARGET defs\n",
    "target_bitname = {\n",
    "    \"ELG\": \"ELG\",\n",
    "    \"LRG\": \"LRG\",\n",
    "    \"QSO-tracer\": \"QSO\",\n",
    "    \"QSO-lyman\": \"QSO\",\n",
    "    \"FAKE-high\": \"ELG\",\n",
    "    \"standards\": \"STD_WD\",\n",
    "    \"sky\": \"SKY\"\n",
    "}\n",
    "\n",
    "# Target priorities\n",
    "target_priority = {\n",
    "    \"ELG\": 3000,\n",
    "    \"LRG\": 3200,\n",
    "    \"QSO-tracer\": 3400,\n",
    "    \"QSO-lyman\": 3400,\n",
    "    \"FAKE-high\": 1000,\n",
    "    \"standards\": 0,\n",
    "    \"sky\": 0\n",
    "}\n",
    "\n",
    "# Target requested number of observations\n",
    "target_numobs = {\n",
    "    \"ELG\": 1,\n",
    "    \"LRG\": 1,\n",
    "    \"QSO-tracer\": 1,\n",
    "    \"QSO-lyman\": 4,\n",
    "    \"FAKE-high\": 1,\n",
    "    \"standards\": 0,\n",
    "    \"sky\": 0\n",
    "}\n",
    "\n",
    "# Plotting color\n",
    "target_pltcolor = {\n",
    "    \"ELG\": (0.12156862745098039, 0.4666666666666667, 0.7058823529411765),\n",
    "    \"LRG\": (1.0, 0.4980392156862745, 0.054901960784313725),\n",
    "    \"QSO-tracer\": (0.17254901960784313, 0.6274509803921569, 0.17254901960784313),\n",
    "    \"QSO-lyman\": (0.8392156862745098, 0.15294117647058825, 0.1568627450980392),\n",
    "    \"FAKE-high\": (0.5803921568627451, 0.403921568627451, 0.7411764705882353),\n",
    "    \"standards\": (1.0, 0.832, 0.0),\n",
    "    \"sky\": (0.0, 1.0, 1.0)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Working directory (where input and output files will be written).  Change this to wherever you like.\n",
    "workdir = os.path.join(os.environ[\"HOME\"], \"scratch\", \"desi\", \"tutorials\", \"fiberassign_part1\")\n",
    "# workdir = os.path.join(os.environ[\"SCRATCH\"], \"desi\", \"tutorials\", \"fiberassign_part1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper Functions\n",
    "\n",
    "Here we define some functions to simulate tiles and uniform random target distributions, as well as some functions to help with plotting later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim_tiles(path, ntile):\n",
    "    \"\"\"Function to generate some co-incident tiles.\n",
    "    \"\"\"\n",
    "    tile_dtype = np.dtype([\n",
    "        (\"TILEID\", \"i4\"),\n",
    "        (\"RA\", \"f8\"),\n",
    "        (\"DEC\", \"f8\"),\n",
    "        (\"IN_DESI\", \"i4\"),\n",
    "        (\"PROGRAM\", \"S6\"),\n",
    "        (\"OBSCONDITIONS\", \"i4\")\n",
    "    ])\n",
    "    fdata = np.zeros(ntile, dtype=tile_dtype)\n",
    "    for i in range(ntile):\n",
    "        fdata[i] = (1234+i, tile_ra, tile_dec, 1, \"DARK\", 1)\n",
    "\n",
    "    if os.path.isfile(path):\n",
    "        os.remove(path)\n",
    "    fd = fitsio.FITS(path, \"rw\")\n",
    "\n",
    "    header = dict()\n",
    "    header[\"FBAVER\"] = fba_version\n",
    "    fd.write(fdata, header=header)\n",
    "\n",
    "    return\n",
    "\n",
    "\n",
    "def sim_targets(tgtype, tgoffset, path=None, density=5000.0, priority=0, numobs=0, tgbits=0):\n",
    "    target_cols = OrderedDict([\n",
    "        (\"TARGETID\", \"i8\"),\n",
    "        (\"RA\", \"f8\"),\n",
    "        (\"DEC\", \"f8\"),\n",
    "        (\"DESI_TARGET\", \"i8\"),\n",
    "        (\"BRICKID\", \"i4\"),\n",
    "        (\"BRICK_OBJID\", \"i4\"),\n",
    "        (\"BRICKNAME\", \"a8\"),\n",
    "        (\"PRIORITY\", \"i4\"),\n",
    "        (\"SUBPRIORITY\", \"f8\"),\n",
    "        (\"OBSCONDITIONS\", \"i4\"),\n",
    "        (\"NUMOBS_MORE\", \"i4\"),\n",
    "        (\"FIBERFLUX_G\", \"f4\"),\n",
    "        (\"FIBERFLUX_R\", \"f4\"),\n",
    "        (\"FIBERFLUX_Z\", \"f4\"),\n",
    "        (\"FIBERFLUX_IVAR_G\", \"f4\"),\n",
    "        (\"FIBERFLUX_IVAR_R\", \"f4\"),\n",
    "        (\"FIBERFLUX_IVAR_Z\", \"f4\")\n",
    "    ])\n",
    "\n",
    "    target_dtype = np.dtype([(x, y) for x, y in target_cols.items()])\n",
    "\n",
    "    ndim = np.sqrt(density)\n",
    "    nra = int(ndim * (patch_ra_max - patch_ra_min))\n",
    "    ndec = int(ndim * (patch_dec_max - patch_dec_min))\n",
    "    ntarget = nra * ndec\n",
    "\n",
    "    fdata = np.zeros(ntarget, dtype=target_dtype)\n",
    "    fdata[\"TARGETID\"][:] = tgoffset + np.arange(ntarget)\n",
    "    fdata[\"RA\"][:] = np.random.uniform(low=patch_ra_min, high=patch_ra_max, size=ntarget)\n",
    "    fdata[\"DEC\"][:] = np.random.uniform(low=patch_dec_min, high=patch_dec_max, size=ntarget)\n",
    "    fdata[\"OBSCONDITIONS\"][:] = np.ones(ntarget, dtype=np.int32)\n",
    "    fdata[\"SUBPRIORITY\"][:] = np.random.uniform(low=0.0, high=1.0, size=ntarget)\n",
    "\n",
    "    sky_mask = desi_mask[\"SKY\"].mask\n",
    "    suppsky_mask = desi_mask[\"SUPP_SKY\"].mask\n",
    "    std_mask = desi_mask[\"STD_BRIGHT\"].mask\n",
    "\n",
    "    if tgtype == TARGET_TYPE_SKY:\n",
    "        fdata[\"PRIORITY\"][:] = np.zeros(ntarget, dtype=np.int32)\n",
    "        fdata[\"DESI_TARGET\"][:] |= sky_mask\n",
    "    elif tgtype == TARGET_TYPE_STANDARD:\n",
    "        fdata[\"PRIORITY\"][:] = priority * np.ones(ntarget, dtype=np.int32)\n",
    "        fdata[\"DESI_TARGET\"][:] |= std_mask\n",
    "    elif tgtype == TARGET_TYPE_SUPPSKY:\n",
    "        fdata[\"PRIORITY\"][:] = np.zeros(ntarget, dtype=np.int32)\n",
    "        fdata[\"DESI_TARGET\"][:] |= suppsky_mask\n",
    "    elif tgtype == TARGET_TYPE_SCIENCE:\n",
    "        fdata[\"PRIORITY\"][:] = priority * np.ones(ntarget, dtype=np.int32)\n",
    "        fdata[\"NUMOBS_MORE\"][:] = numobs\n",
    "        fdata[\"DESI_TARGET\"][:] |= desi_mask[tgbits].mask\n",
    "    else:\n",
    "        raise RuntimeError(\"unknown target type\")\n",
    "\n",
    "    if path is None:\n",
    "        # Just return the data table\n",
    "        return fdata\n",
    "    else:\n",
    "        # We are writing the output to a file.  Return the number of targets.\n",
    "        if os.path.isfile(path):\n",
    "            os.remove(path)\n",
    "        fd = fitsio.FITS(path, \"rw\")\n",
    "        header = dict()\n",
    "        header[\"FBAVER\"] = fba_version\n",
    "        fd.write(fdata, header=header)\n",
    "        return ntarget\n",
    "\n",
    "\n",
    "def plot_assignment_stats(\n",
    "        wd, outname, title, tile_ids, classes,\n",
    "        hist_tgassign, hist_tgavail, hist_tgconsid\n",
    "):\n",
    "    fiberbins = 0.01 * np.arange(101)\n",
    "    tgbins = {\n",
    "        \"ELG\": np.arange(0, 22220, 220),\n",
    "        \"LRG\": np.arange(0, 3232, 32),\n",
    "        \"QSO-tracer\": np.arange(0, 1515, 15),\n",
    "        \"QSO-lyman\": np.arange(0, 808, 8),\n",
    "        \"FAKE-high\": np.arange(0, 36360, 360),\n",
    "        \"standards\": np.arange(0, 3636, 36),\n",
    "        \"sky\": np.arange(0, 40400, 400)\n",
    "    }\n",
    "    \n",
    "    for tid in tile_ids:\n",
    "        # The plots\n",
    "        figfiber = plt.figure(figsize=(12, 6))\n",
    "        figtarget = plt.figure(figsize=(12, 12))\n",
    "        \n",
    "        # Plot the positioner assignment fractions\n",
    "        axfiber = figfiber.add_subplot(1, 1, 1)\n",
    "        \n",
    "        for tgclass in classes:\n",
    "            if np.sum(hist_tgassign[tid][tgclass]) > 0:\n",
    "                axfiber.hist(\n",
    "                    np.array(hist_tgassign[tid][tgclass]) / 5000, \n",
    "                    fiberbins, \n",
    "                    align=\"mid\", \n",
    "                    alpha=0.8, \n",
    "                    label=\"{}\".format(tgclass),\n",
    "                    color=target_pltcolor[tgclass]\n",
    "                )\n",
    "        axfiber.set_xlabel(\"Fraction of Assigned Positioners\", fontsize=\"large\")\n",
    "        axfiber.set_ylabel(\"Realization Counts\", fontsize=\"large\")\n",
    "        axfiber.set_ylim(0, 100)\n",
    "        axfiber.set_title(\n",
    "            \"Tile {}:  Positioner Assignment : {}\".format(tid, title)\n",
    "        )\n",
    "        axfiber.legend()\n",
    "        figfiber.savefig(\n",
    "            os.path.join(wd, \"fiberassign_{}_fibers_tile-{}.pdf\".format(outname, tid)),\n",
    "            dpi=300, \n",
    "            format=\"pdf\"\n",
    "        )\n",
    "        figfiber.show()\n",
    "            \n",
    "        nplot = len(classes)\n",
    "        pcols = 2\n",
    "        prows = (nplot + 1) // pcols\n",
    "        poff = 1\n",
    "        for tgclass in classes:\n",
    "            if np.sum(hist_tgavail[tid][tgclass]) == 0:\n",
    "                continue\n",
    "            axtarget = figtarget.add_subplot(prows, pcols, poff)\n",
    "            axtarget.hist(\n",
    "                np.array(hist_tgavail[tid][tgclass]), \n",
    "                tgbins[tgclass], \n",
    "                align=\"mid\", \n",
    "                alpha=0.4, \n",
    "                label=\"{} Reachable\".format(tgclass),\n",
    "                color=(0.7, 0.7, 0.7)\n",
    "            )\n",
    "            axtarget.hist(\n",
    "                np.array(hist_tgconsid[tid][tgclass]), \n",
    "                tgbins[tgclass], \n",
    "                align=\"mid\", \n",
    "                alpha=0.4, \n",
    "                label=\"{} Considered\".format(tgclass),\n",
    "                color=target_pltcolor[tgclass]\n",
    "            )\n",
    "            axtarget.hist(\n",
    "                np.array(hist_tgassign[tid][tgclass]), \n",
    "                tgbins[tgclass], \n",
    "                align=\"mid\", \n",
    "                alpha=0.8, \n",
    "                label=\"{} Assigned\".format(tgclass),\n",
    "                color=target_pltcolor[tgclass]\n",
    "            )\n",
    "            axtarget.set_xlabel(\"Targets Reachable, Considered, and Assigned\", fontsize=\"large\")\n",
    "            axtarget.set_ylabel(\"Realization Counts\", fontsize=\"large\")\n",
    "            axtarget.set_ylim(0, 100)\n",
    "            axtarget.legend()\n",
    "            poff += 1\n",
    "        \n",
    "        figtarget.suptitle(\n",
    "            \"Tile {}:  {}\".format(tid, title), fontsize=\"x-large\"\n",
    "        )\n",
    "        figtarget.tight_layout(rect=[0, 0, 1, 0.97])\n",
    "        figtarget.savefig(\n",
    "            os.path.join(wd, \"fiberassign_{}_targets_tile-{}.pdf\".format(outname, tid)),\n",
    "            dpi=300, \n",
    "            format=\"pdf\"\n",
    "        )\n",
    "        figtarget.show()\n",
    "\n",
    "        \n",
    "def create_histogram_dicts(tiles, classes):\n",
    "    hist_tgassign = dict()\n",
    "    hist_tgconsid = dict()\n",
    "    hist_tgavail = dict()\n",
    "    for tid in tiles:\n",
    "        hist_tgassign[tid] = dict()\n",
    "        hist_tgavail[tid] = dict()\n",
    "        hist_tgconsid[tid] = dict()\n",
    "        for tgclass in classes:\n",
    "            hist_tgassign[tid][tgclass] = list()\n",
    "            hist_tgavail[tid][tgclass] = list()\n",
    "            hist_tgconsid[tid][tgclass] = list()\n",
    "    return (hist_tgassign, hist_tgavail, hist_tgconsid)\n",
    "\n",
    "\n",
    "def accum_histogram_data(\n",
    "    tile, classes, tile_assign, tile_avail, id2class, tgobs, \n",
    "    hist_tgassign, hist_tgavail, hist_tgconsid, verbose=False\n",
    "):      \n",
    "    # The locations that were assigned to a target\n",
    "    locs_assigned = [x for x, y in tile_assign.items() if y >= 0]\n",
    "    \n",
    "    if verbose:\n",
    "        print(\n",
    "            \"    tile {}:  {} positioners ({:0.1f}%) assigned a target\".format(\n",
    "                tile, len(locs_assigned), 100*(len(locs_assigned)/5000)\n",
    "            )\n",
    "        )\n",
    "\n",
    "    # Compute the numbers of each target class that were assigned and available\n",
    "    assign_by_class = dict()\n",
    "    avail_by_class = dict()\n",
    "    consider_by_class = dict()\n",
    "\n",
    "    for tgclass in classes:\n",
    "        assign_by_class[tgclass] = list()\n",
    "        avail_by_class[tgclass] = list()\n",
    "        consider_by_class[tgclass] = list()\n",
    "\n",
    "    for loc in locs_assigned:\n",
    "        tgid = tile_assign[loc]\n",
    "        tgstr = id2class(tgid)\n",
    "        assign_by_class[tgstr].append(tgid)\n",
    "        for avtg in tile_avail[loc]:\n",
    "            avtgstr = id2class(avtg)\n",
    "            avail_by_class[avtgstr].append(avtg)\n",
    "            if avtgstr == \"standards\" or avtgstr == \"sky\" or tgobs[avtg] > 0:\n",
    "                # This target was actually considered for assignment\n",
    "                consider_by_class[avtgstr].append(avtg)\n",
    "    \n",
    "    # Decrement the obs remaining for the assigned targets\n",
    "    for loc in locs_assigned:\n",
    "        tgid = tile_assign[loc]\n",
    "        tgobs[tgid] -= 1\n",
    "\n",
    "    # Accumulate to the target histograms\n",
    "    for tgclass in classes:\n",
    "        hist_tgassign[tid][tgclass].append(len(assign_by_class[tgclass]))\n",
    "        uniq = np.unique(np.array(avail_by_class[tgclass]))\n",
    "        hist_tgavail[tid][tgclass].append(len(uniq))\n",
    "        uniq = np.unique(np.array(consider_by_class[tgclass]))\n",
    "        hist_tgconsid[tid][tgclass].append(len(uniq))\n",
    "    return\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Focalplane Coverage\n",
    "\n",
    "The hardware information that fiberassign gets from the `desimodel` package contains the positioner arm lengths and ranges of motion for all positioners.  We can use that along with a random sampling of locations to illustrate the regions that are not reachable by any positioner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directory for this section\n",
    "wdir = os.path.join(workdir, \"focalplane\")\n",
    "os.makedirs(wdir, exist_ok=True)\n",
    "\n",
    "# Set the random seed to ensure reproducibility of this cell\n",
    "np.random.seed(123456)\n",
    "\n",
    "# Read hardware properties\n",
    "test_date = assign_date\n",
    "# If you uncomment this line, you will get the commissioning\n",
    "# focalplane with restricted positioner reach.\n",
    "#test_date = \"2020-04-01T00:00:00\"\n",
    "hw = load_hardware(rundate=test_date)\n",
    "\n",
    "# Simulate a single tile and load\n",
    "tfile = os.path.join(wdir, \"footprint.fits\")\n",
    "sim_tiles(tfile, 1)\n",
    "tiles = load_tiles(tiles_file=tfile)\n",
    "\n",
    "# Generate target table at high density\n",
    "tgdensity = 10000\n",
    "tgdata = sim_targets(\n",
    "    TARGET_TYPE_SCIENCE, \n",
    "    0,\n",
    "    density=tgdensity,\n",
    "    priority=0,\n",
    "    numobs=0,\n",
    "    tgbits=\"SKY\"\n",
    ")\n",
    "print(\"{} fake targets\".format(len(tgdata)))\n",
    "tgs = Targets()\n",
    "load_target_table(tgs, tgdata)\n",
    "\n",
    "# Compute the targets available to all positioners\n",
    "tree = TargetTree(tgs, 0.01)\n",
    "tgsavail = TargetsAvailable(hw, tgs, tiles, tree)\n",
    "del tree\n",
    "\n",
    "# Availability for this single tile\n",
    "tid = tiles.id[0]\n",
    "tavail = tgsavail.tile_data(tid)\n",
    "\n",
    "# Get the unique set of target IDs reachable by *any* positioner\n",
    "reachable = np.unique([x for loc, tgl in tavail.items() for x in tgl])\n",
    "print(\"{} reachable targets\".format(len(reachable)))\n",
    "\n",
    "# Now get the rows of the original target table that are NOT in this set\n",
    "unreachable_rows = np.where(\n",
    "    np.isin(tgdata[\"TARGETID\"], reachable, invert=True)\n",
    ")[0]\n",
    "print(\"{} un-reachable targets\".format(len(unreachable_rows)))\n",
    "\n",
    "# Now we would like to plot the targets that were not reachable.  Also\n",
    "# plot the positioner angle ranges.\n",
    "\n",
    "locs = np.array(hw.locations)\n",
    "loc_theta_offset = hw.loc_theta_offset\n",
    "loc_theta_min = hw.loc_theta_min\n",
    "loc_theta_max = hw.loc_theta_max\n",
    "loc_phi_offset = hw.loc_phi_offset\n",
    "loc_phi_min = hw.loc_phi_min\n",
    "loc_phi_max = hw.loc_phi_max\n",
    "\n",
    "theta_off = np.array([loc_theta_offset[x]*180/np.pi for x in locs])\n",
    "theta_min = np.array([loc_theta_min[x]*180/np.pi for x in locs])\n",
    "theta_max = np.array([loc_theta_max[x]*180/np.pi for x in locs])\n",
    "\n",
    "phi_off = np.array([loc_phi_offset[x]*180/np.pi for x in locs])\n",
    "phi_min = np.array([loc_phi_min[x]*180/np.pi for x in locs])\n",
    "phi_max = np.array([loc_phi_max[x]*180/np.pi for x in locs])\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = fig.add_subplot(1, 2, 1)\n",
    "ax.scatter(locs, theta_off, marker=\"o\", s=0.5, label=\"theta off\")\n",
    "ax.scatter(locs, theta_min, marker=\"o\", s=0.5, label=\"theta min\")\n",
    "ax.scatter(locs, theta_max, marker=\"o\", s=0.5, label=\"theta max\")\n",
    "ax.set_xlabel(\"Positioner Location\", fontsize=\"large\")\n",
    "ax.set_ylabel(\"Degrees\", fontsize=\"large\")\n",
    "ax.legend()\n",
    "\n",
    "ax = fig.add_subplot(1, 2, 2)\n",
    "ax.scatter(locs, phi_off, marker=\"o\", s=0.5, label=\"phi off\")\n",
    "ax.scatter(locs, phi_min, marker=\"o\", s=0.5, label=\"phi min\")\n",
    "ax.scatter(locs, phi_max, marker=\"o\", s=0.5, label=\"phi max\")\n",
    "ax.set_xlabel(\"Positioner Location\", fontsize=\"large\")\n",
    "ax.set_ylabel(\"Degrees\", fontsize=\"large\")\n",
    "ax.legend()\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure(figsize=(10, 10))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.set_aspect(\"equal\")\n",
    "\n",
    "center_mm = hw.loc_pos_curved_mm\n",
    "theta_arm = hw.loc_theta_arm\n",
    "phi_arm = hw.loc_phi_arm\n",
    "patrol_buffer = hw.patrol_buffer_mm\n",
    "    \n",
    "for loc in hw.locations:\n",
    "    patrol_rad = theta_arm[loc] + phi_arm[loc] - patrol_buffer\n",
    "    patrol = plt.Circle(\n",
    "        (center_mm[loc][0], center_mm[loc][1]), \n",
    "        radius=patrol_rad, \n",
    "        fc=\"red\",\n",
    "        ec=\"none\", \n",
    "        alpha=0.1\n",
    "    )\n",
    "    ax.add_artist(patrol)\n",
    "\n",
    "tgxy = hw.radec2xy_multi(\n",
    "    tiles.ra[0], \n",
    "    tiles.dec[0],\n",
    "    tiles.obstheta[0],\n",
    "    tgdata[\"RA\"][unreachable_rows],\n",
    "    tgdata[\"DEC\"][unreachable_rows],\n",
    "    False,\n",
    "    1\n",
    ")\n",
    "tgx = np.array([x[0] for x in tgxy])\n",
    "tgy = np.array([x[1] for x in tgxy])\n",
    "\n",
    "ax.scatter(\n",
    "    tgx,\n",
    "    tgy,\n",
    "    color=(0.6, 0.6, 0.6),\n",
    "    marker=\".\",\n",
    "    s=0.01\n",
    ")\n",
    "\n",
    "ax.set_xlabel(\"Millimeters\", fontsize=\"large\")\n",
    "ax.set_ylabel(\"Millimeters\", fontsize=\"large\")\n",
    "pfile = os.path.join(wdir, \"coverage_{}.pdf\".format(tgdensity))\n",
    "plt.savefig(pfile, dpi=300, format=\"pdf\")\n",
    "plt.close()\n",
    "\n",
    "print(\"PDF written to {}\\n  Displaying low-res image inline...\".format(pfile))\n",
    "png_file = os.path.join(wdir, \"coverage_{}.png\".format(tgdensity))\n",
    "cmd = 'convert -density 200 {} {}'.format(pfile, png_file)\n",
    "sp.check_call(cmd, stderr=sp.STDOUT, shell=True)\n",
    "img = Image(filename=(png_file))\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effects of Positioner Geometry\n",
    "\n",
    "The exclusion zone around each positioner and its range of motion leads to an average limit on the fraction of targets that can be assigned for a given target density.  We can see this intrinsic limit by simulating realizations of a single population of random targets for several densities. The plots of the positioner assignments for the first realization can take a while- you can comment out that if you like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directory for this section\n",
    "wdir = os.path.join(workdir, \"geometry\")\n",
    "os.makedirs(wdir, exist_ok=True)\n",
    "\n",
    "# Set the random seed to ensure reproducibility of this cell\n",
    "np.random.seed(123456)\n",
    "\n",
    "# Read hardware properties\n",
    "hw = load_hardware(rundate=assign_date)\n",
    "\n",
    "# Simulate a single tile and load\n",
    "tfile = os.path.join(wdir, \"footprint.fits\")\n",
    "sim_tiles(tfile, 1)\n",
    "tiles = load_tiles(tiles_file=tfile)\n",
    "\n",
    "# Accumulation dictionaries\n",
    "(hist_tgassign, hist_tgavail, hist_tgconsid) = \\\n",
    "    create_histogram_dicts(tiles.id, list(target_density.keys()))\n",
    "\n",
    "# Consider these 5 classes.\n",
    "target_classes = [\"ELG\", \"LRG\", \"QSO-tracer\", \"QSO-lyman\", \"FAKE-high\"]\n",
    "\n",
    "for tgclass in target_classes:\n",
    "    tgdensity = target_density[tgclass]\n",
    "\n",
    "    print(\"Working on target density {}...\".format(tgdensity))\n",
    "    \n",
    "    # For each realization...\n",
    "    for mc in range(100):\n",
    "        # Generate target table\n",
    "        tgdata = sim_targets(\n",
    "            TARGET_TYPE_SCIENCE, \n",
    "            0,\n",
    "            density=tgdensity,\n",
    "            priority=target_priority[tgclass],\n",
    "            numobs=target_numobs[tgclass],\n",
    "            tgbits=target_bitname[tgclass]\n",
    "        )\n",
    "        tgs = Targets()\n",
    "        load_target_table(tgs, tgdata)\n",
    "        \n",
    "        # Make a working copy of the obs remaining for all targets, so we can decrement\n",
    "        # and track the targets that were considered (which had obsremain > 0)\n",
    "        tg_obs = dict()\n",
    "        for tgid in tgs.ids():\n",
    "            tgprops = tgs.get(tgid)\n",
    "            tg_obs[tgid] = tgprops.obsremain\n",
    "\n",
    "        # Create a hierarchical triangle mesh lookup of the targets positions\n",
    "        tree = TargetTree(tgs, 0.01)\n",
    "\n",
    "        # Compute the targets available to each fiber for each tile.\n",
    "        tgsavail = TargetsAvailable(hw, tgs, tiles, tree)\n",
    "\n",
    "        # Free the tree\n",
    "        del tree\n",
    "\n",
    "        # Compute the fibers on all tiles available for each target\n",
    "        favail = LocationsAvailable(tgsavail)\n",
    "\n",
    "        # Create assignment object\n",
    "        asgn = Assignment(tgs, tgsavail, favail)\n",
    "\n",
    "        # assignment of science targets\n",
    "        asgn.assign_unused(TARGET_TYPE_SCIENCE)\n",
    "        \n",
    "        # Targets available for all tile / locs\n",
    "        tgsavail = asgn.targets_avail()\n",
    "        \n",
    "        # Get the assignment for this one tile\n",
    "        tid = tiles.id[0]\n",
    "        tassign = asgn.tile_location_target(tid)\n",
    "        \n",
    "        # Targets available for this tile\n",
    "        tavail = tgsavail.tile_data(tid)\n",
    "        \n",
    "        # Plot the first realization\n",
    "        # THIS IS SLOW- after doing it once you can comment it out to save time.\n",
    "#         if mc == 0:\n",
    "#             pfile = os.path.join(wdir, \"assignment_density_{}.pdf\".format(tgdensity))\n",
    "#             #print(locs_assigned, flush=True)\n",
    "#             plot_assignment_tile(\n",
    "#                 hw, tgs, tiles.id[0], tiles.ra[0], tiles.dec[0], tiles.obstheta[0],\n",
    "#                 tassign, tile_avail=tgsavail.tile_data(tiles.id[0]), real_shapes=True, \n",
    "#                 outfile=pfile, figsize=8\n",
    "#             )\n",
    "#             print(\"Full resolution PDF written to {}, displaying low-res image inline...\".format(pfile))\n",
    "#             png_file = os.path.join(wdir, \"assignment_density_{}.png\".format(tgdensity))\n",
    "#             cmd = 'convert -density 200 {} {}'.format(pfile, png_file)\n",
    "#             sp.check_call(cmd, stderr=sp.STDOUT, shell=True)\n",
    "#             img = Image(filename=(png_file))\n",
    "#             display(img)\n",
    "        \n",
    "        # Accumulate results\n",
    "        \n",
    "        def target_id_class(target_id):\n",
    "            # We are doing one class at a time in this example\n",
    "            return tgclass\n",
    "        \n",
    "        accum_histogram_data(\n",
    "            tid, [tgclass], tassign, tavail, target_id_class, \n",
    "            tg_obs, hist_tgassign, hist_tgavail, hist_tgconsid\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot it\n",
    "\n",
    "plot_assignment_stats(\n",
    "    wdir,\n",
    "    \"geometry\", \n",
    "    \"One Target Class at a Time\", \n",
    "    [tiles.id[0]], \n",
    "    list(target_density.keys()), \n",
    "    hist_tgassign, \n",
    "    hist_tgavail, \n",
    "    hist_tgconsid\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous plots are not \"surprising\".  They simply show that for low (QSO-like) target densities there are fewer collisions and a larger fraction (80-90%) of targets can be assigned, even though most positioners do not get one of these targets.  For a medium (LRG-like) target density we have more collisions and around 70% of available targets can be assigned while slightly less than half of the positioners get a target.  Going to a higher (ELG-like) density we have far more targets than can be assigned, but there are sufficient numbers that nearly every positioner gets a target.\n",
    "\n",
    "**NOTE:  remember that these are uniform spatial target distributions with no physically realistic clustering.  Clustering of targets reduces the fraction that can be assigned, due to positioner collisions.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effects of Multiple Target Classes\n",
    "\n",
    "The previous section looked at different target densities but assigned each one in isolation.  Now we examine the consequences of having 3 target classes with different priorities (and densities).  The goal of this section is to demonstrate how priorities affect the resulting assignment.  The target locations are still spatially uniform with no clustering.  We will also demonstrate explicitly that low priority targets never \"take\" positioners from high priority targets.\n",
    "\n",
    "In the case of two QSO populations, these have the same `PRIORITY` value.  However the fiberassign code currently computes an effective \"total priority\" for each target that is a combination of the (integer) `PRIORITY`, the (random from 0.0-1.0) `SUBPRIORITY` and the number of remaining observations.  See the discussion here: https://github.com/desihub/fiberassign/issues/196 .  For now, targets with the same `PRIORITY` and a larger number of observations remaining are given a higher \"total priority\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directory for this section\n",
    "wdir = os.path.join(workdir, \"target_classes\")\n",
    "os.makedirs(wdir, exist_ok=True)\n",
    "\n",
    "# Set the random seed to ensure reproducibility of this cell\n",
    "np.random.seed(123456)\n",
    "\n",
    "# Read hardware properties\n",
    "hw = load_hardware(rundate=assign_date)\n",
    "\n",
    "# Simulate a single tile and load\n",
    "tfile = os.path.join(wdir, \"footprint.fits\")\n",
    "sim_tiles(tfile, 1)\n",
    "tiles = load_tiles(tiles_file=tfile)\n",
    "\n",
    "# Consider these 4 classes.\n",
    "target_classes = [\"ELG\", \"LRG\", \"QSO-tracer\", \"QSO-lyman\"]\n",
    "\n",
    "# Accumulation dictionaries\n",
    "(hist_tgassign, hist_tgavail, hist_tgconsid) = \\\n",
    "    create_histogram_dicts(tiles.id, target_classes)\n",
    "\n",
    "# For each realization...\n",
    "for mc in range(100):\n",
    "\n",
    "    # Generate target tables separately for each target class so we can combine them later.\n",
    "    tg_tables = dict()\n",
    "    tgid_offset = dict()\n",
    "    ntarget = dict()\n",
    "    toff = 0\n",
    "    for tgclass in target_classes:\n",
    "        tg_tables[tgclass] = sim_targets(\n",
    "            TARGET_TYPE_SCIENCE, \n",
    "            toff,\n",
    "            density=target_density[tgclass],\n",
    "            priority=target_priority[tgclass],\n",
    "            numobs=target_numobs[tgclass],\n",
    "            tgbits=target_bitname[tgclass]\n",
    "        )\n",
    "        tgid_offset[tgclass] = toff\n",
    "        ntarget[tgclass] = len(tg_tables[tgclass])\n",
    "        toff += ntarget[tgclass]\n",
    "        \n",
    "    def target_id_class(target_id):\n",
    "        \"\"\"Helper function to go from target ID back to the class\"\"\"\n",
    "        for cls in target_classes:\n",
    "            if (target_id >= tgid_offset[cls]) and (target_id < tgid_offset[cls] + ntarget[cls]):\n",
    "                return cls\n",
    "        return None\n",
    "\n",
    "    # Assign just the QSOs separately, and save for comparison.\n",
    "    \n",
    "    qso_tgs = Targets()\n",
    "    load_target_table(qso_tgs, tg_tables[\"QSO-tracer\"])\n",
    "    load_target_table(qso_tgs, tg_tables[\"QSO-lyman\"])\n",
    "\n",
    "    tree = TargetTree(qso_tgs, 0.01)\n",
    "    tgsavail = TargetsAvailable(hw, qso_tgs, tiles, tree)\n",
    "    del tree\n",
    "    favail = LocationsAvailable(tgsavail)\n",
    "    asgn = Assignment(qso_tgs, tgsavail, favail)\n",
    "    asgn.assign_unused(TARGET_TYPE_SCIENCE)\n",
    "    tid = tiles.id[0]\n",
    "    qso_assign = dict(asgn.tile_location_target(tid))\n",
    "\n",
    "    # Now assign all targets\n",
    "\n",
    "    all_tgs = Targets()\n",
    "    for tgclass in target_classes:\n",
    "        load_target_table(all_tgs, tg_tables[tgclass])\n",
    "        \n",
    "    # Make a working copy of the obs remaining for all targets, so we can decrement\n",
    "    # and track the targets that were considered (which had obsremain > 0)\n",
    "    tg_obs = dict()\n",
    "    for tgid in all_tgs.ids():\n",
    "        tgprops = all_tgs.get(tgid)\n",
    "        tg_obs[tgid] = tgprops.obsremain\n",
    "\n",
    "    tree = TargetTree(all_tgs, 0.01)\n",
    "    tgsavail = TargetsAvailable(hw, all_tgs, tiles, tree)\n",
    "    del tree\n",
    "    favail = LocationsAvailable(tgsavail)\n",
    "    asgn = Assignment(all_tgs, tgsavail, favail)\n",
    "    asgn.assign_unused(TARGET_TYPE_SCIENCE)\n",
    "    tid = tiles.id[0]\n",
    "    all_assign = asgn.tile_location_target(tid)\n",
    "\n",
    "    print(\n",
    "        \"Realization {:02d}:  assigning all targets:  {} positioners assigned\".format(\n",
    "            mc, \n",
    "            int(np.sum([1 for x, y in all_assign.items() if y >= 0]))\n",
    "        )\n",
    "    )\n",
    "\n",
    "    # Verify that QSO assignment for all targets is identical to the case of assigning only the QSOs\n",
    "    qso_pos = [x for x, y in qso_assign.items() if y >= 0]\n",
    "    mismatch = int(np.sum([1 for x, y in qso_assign.items() if all_assign[x] != y]))\n",
    "    if mismatch > 0:\n",
    "        print(\n",
    "            \"Realization {:02d}:    WARNING: {} positioners from QSO-only assignment changed!\".format(\n",
    "                mc, \n",
    "                mismatch\n",
    "            )\n",
    "        )\n",
    "    \n",
    "    # Targets available for all tile / locs\n",
    "    tgsavail = asgn.targets_avail()\n",
    "\n",
    "    # available for this tile\n",
    "    avail = tgsavail.tile_data(tid)\n",
    "    \n",
    "    # Accumulate results\n",
    "    accum_histogram_data(\n",
    "        tid, target_classes, all_assign, avail, target_id_class, \n",
    "        tg_obs, hist_tgassign, hist_tgavail, hist_tgconsid,\n",
    "        verbose=False\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results\n",
    "\n",
    "plot_assignment_stats(\n",
    "    wdir,\n",
    "    \"allscience\", \n",
    "    \"All Science Classes for One Tile\", \n",
    "    [tiles.id[0]], \n",
    "    target_classes, \n",
    "    hist_tgassign, \n",
    "    hist_tgavail, \n",
    "    hist_tgconsid\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exercise above shows several things:\n",
    "\n",
    "1.  The highest \"total priority\" targets (which also have a low density) are assigned first and nearly all get assigned to a fiber.\n",
    "\n",
    "2.  As the target \"total priority\" decreases (and the density increases), these targets take up larger fractions of the positioners.\n",
    "\n",
    "3.  During this first-pass assignment, about 98% of positioners get assigned to a science target (remember we are not assigning standards and sky in this excercise).\n",
    "\n",
    "4.  The introduction of lower priority targets **does not** affect the assignment of higher priority targets.\n",
    "\n",
    "Compare the \"Fraction of Assigned Positioners\" plot to the one from the previous section.  You can see that the highest priority class has the same assignment fraction, since it is assigned first.  The other target classes have smaller fractions of the positioners because they are forced to use only positioners leftover after the higher priority targets are assigned."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effects of Multiple Passes\n",
    "\n",
    "So far we have worked with one tile.  However the QSO-lyman population has 4 requested observations.  If we have multiple observations of the same point on the sky, we would expect these targets to continue to get preferential assignments until their observations remaining equals \"1\" at which point they have the same observations remaining as the other QSO-tracer population.  At that point the random subriority will determine the total priority ordering.  And of course positioners will continue to be assigned to all target classes as possible.\n",
    "\n",
    "In the real survey, tiles from multiple \"layers\" (passes over the footprint) are dithered so that they do not overlap perfectly.  For this study of one pointing on the sky, we will use 4 tiles (passes) centered on the same RA / DEC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directory for this section\n",
    "wdir = os.path.join(workdir, \"multipass\")\n",
    "os.makedirs(wdir, exist_ok=True)\n",
    "\n",
    "# Set the random seed to ensure reproducibility of this cell\n",
    "np.random.seed(123456)\n",
    "\n",
    "# Read hardware properties\n",
    "hw = load_hardware(rundate=assign_date)\n",
    "\n",
    "# Simulate 5 coincident tiles and load\n",
    "tfile = os.path.join(wdir, \"footprint.fits\")\n",
    "sim_tiles(tfile, 4)\n",
    "tiles = load_tiles(tiles_file=tfile)\n",
    "\n",
    "# Consider these 4 classes.\n",
    "target_classes = [\"ELG\", \"LRG\", \"QSO-tracer\", \"QSO-lyman\"]\n",
    "\n",
    "# Accumulation dictionaries\n",
    "(hist_tgassign, hist_tgavail, hist_tgconsid) = \\\n",
    "    create_histogram_dicts(tiles.id, target_classes)\n",
    "\n",
    "# For each realization...\n",
    "for mc in range(100):\n",
    "\n",
    "    # Generate target tables\n",
    "    tg_tables = dict()\n",
    "    tgid_offset = dict()\n",
    "    ntarget = dict()\n",
    "    toff = 0\n",
    "    for tgclass in target_classes:\n",
    "        tg_tables[tgclass] = sim_targets(\n",
    "            TARGET_TYPE_SCIENCE, \n",
    "            toff,\n",
    "            density=target_density[tgclass],\n",
    "            priority=target_priority[tgclass],\n",
    "            numobs=target_numobs[tgclass],\n",
    "            tgbits=target_bitname[tgclass]\n",
    "        )\n",
    "        tgid_offset[tgclass] = toff\n",
    "        ntarget[tgclass] = len(tg_tables[tgclass])\n",
    "        toff += ntarget[tgclass]\n",
    "        \n",
    "    def target_id_class(target_id):\n",
    "        for cls in target_classes:\n",
    "            if (target_id >= tgid_offset[cls]) and (target_id < tgid_offset[cls] + ntarget[cls]):\n",
    "                return cls\n",
    "        return None\n",
    "\n",
    "    # Load the targets\n",
    "\n",
    "    all_tgs = Targets()\n",
    "    for tgclass in target_classes:\n",
    "        load_target_table(all_tgs, tg_tables[tgclass])\n",
    "        \n",
    "    # Make a working copy of the obs remaining for all targets, so we can decrement\n",
    "    # and track the targets that were considered (which had obsremain > 0)\n",
    "    tg_obs = dict()\n",
    "    for tgid in all_tgs.ids():\n",
    "        tgprops = all_tgs.get(tgid)\n",
    "        tg_obs[tgid] = tgprops.obsremain\n",
    "\n",
    "    tree = TargetTree(all_tgs, 0.01)\n",
    "    tgsavail = TargetsAvailable(hw, all_tgs, tiles, tree)\n",
    "    del tree\n",
    "    favail = LocationsAvailable(tgsavail)\n",
    "    asgn = Assignment(all_tgs, tgsavail, favail)\n",
    "    asgn.assign_unused(TARGET_TYPE_SCIENCE)\n",
    "    \n",
    "    print(\"Realization {:02d}:\".format(mc))\n",
    "    \n",
    "    # Accumulate results for all tiles\n",
    "    for tid in tiles.id:\n",
    "        # Assignment for this tile\n",
    "        tile_assign = asgn.tile_location_target(tid)\n",
    "\n",
    "        # Targets available for all tile / locs\n",
    "        tgsavail = asgn.targets_avail()\n",
    "\n",
    "        # available for this tile\n",
    "        tile_avail = tgsavail.tile_data(tid)\n",
    "        \n",
    "        # Accumulate results\n",
    "        accum_histogram_data(\n",
    "            tid, target_classes, tile_assign, tile_avail, target_id_class, \n",
    "            tg_obs, hist_tgassign, hist_tgavail, hist_tgconsid,\n",
    "            verbose=True\n",
    "        )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results\n",
    "\n",
    "plot_assignment_stats(\n",
    "    wdir,\n",
    "    \"multipass\", \n",
    "    \"All Science Classes for {} Passes\".format(len(tiles.id)), \n",
    "    tiles.id, \n",
    "    target_classes, \n",
    "    hist_tgassign, \n",
    "    hist_tgavail, \n",
    "    hist_tgconsid\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the results from the first pass are identical to the previous exercise.  This is because all targets have at least one requested observation.  In the subsequent passes, we see the impact of some targets having met their required observations.  Although the numbers of \"reachable\" targets remains stable, we can see that the number of targets \"considered\" for assignment (i.e. that have remaining observations) goes down with subsequent tiles.  This is most notable with the QSO-lyman target class that has 4 requested observations per target."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sky and Standards\n",
    "\n",
    "After assigning science targets, we must free up a certain number of fibers for calibration standards and sky."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directory for this section\n",
    "wdir = os.path.join(workdir, \"sky-std\")\n",
    "os.makedirs(wdir, exist_ok=True)\n",
    "\n",
    "# Set the random seed to ensure reproducibility of this cell\n",
    "np.random.seed(123456)\n",
    "\n",
    "# Read hardware properties\n",
    "hw = load_hardware(rundate=assign_date)\n",
    "\n",
    "# Simulate 5 coincident tiles and load\n",
    "tfile = os.path.join(wdir, \"footprint.fits\")\n",
    "sim_tiles(tfile, 4)\n",
    "tiles = load_tiles(tiles_file=tfile)\n",
    "\n",
    "# Consider these 4 classes.\n",
    "target_classes = [\"ELG\", \"LRG\", \"QSO-tracer\", \"QSO-lyman\"]\n",
    "\n",
    "all_classes = list(target_classes)\n",
    "all_classes.extend([\"standards\", \"sky\"])\n",
    "    \n",
    "# Accumulation dictionaries\n",
    "(hist_tgassign, hist_tgavail, hist_tgconsid) = \\\n",
    "    create_histogram_dicts(tiles.id, all_classes)\n",
    "\n",
    "# For each realization...\n",
    "for mc in range(100):\n",
    "\n",
    "    # Generate target tables\n",
    "    tg_tables = dict()\n",
    "    tgid_offset = dict()\n",
    "    ntarget = dict()\n",
    "    toff = 0\n",
    "    for tgclass in target_classes:\n",
    "        tg_tables[tgclass] = sim_targets(\n",
    "            TARGET_TYPE_SCIENCE, \n",
    "            toff,\n",
    "            density=target_density[tgclass],\n",
    "            priority=target_priority[tgclass],\n",
    "            numobs=target_numobs[tgclass],\n",
    "            tgbits=target_bitname[tgclass]\n",
    "        )\n",
    "        tgid_offset[tgclass] = toff\n",
    "        ntarget[tgclass] = len(tg_tables[tgclass])\n",
    "        toff += ntarget[tgclass]\n",
    "    \n",
    "    tg_tables[\"standards\"] = sim_targets(\n",
    "        TARGET_TYPE_STANDARD, \n",
    "        toff,\n",
    "        density=target_density[\"standards\"],\n",
    "    )\n",
    "    tgid_offset[\"standards\"] = toff\n",
    "    ntarget[\"standards\"] = len(tg_tables[\"standards\"])\n",
    "    toff += ntarget[\"standards\"]\n",
    "    \n",
    "    tg_tables[\"sky\"] = sim_targets(\n",
    "        TARGET_TYPE_SKY, \n",
    "        toff,\n",
    "        density=target_density[\"sky\"],\n",
    "    )\n",
    "    tgid_offset[\"sky\"] = toff\n",
    "    ntarget[\"sky\"] = len(tg_tables[\"sky\"])\n",
    "    toff += ntarget[\"sky\"]\n",
    "    \n",
    "    def target_id_class(target_id):\n",
    "        for cls in all_classes:\n",
    "            if (target_id >= tgid_offset[cls]) and (target_id < tgid_offset[cls] + ntarget[cls]):\n",
    "                return cls\n",
    "        return None\n",
    "\n",
    "    # Load the targets\n",
    "\n",
    "    all_tgs = Targets()\n",
    "    for tgclass in all_classes:\n",
    "        load_target_table(all_tgs, tg_tables[tgclass])\n",
    "        \n",
    "    # Make a working copy of the obs remaining for all targets, so we can decrement\n",
    "    # and track the targets that were considered (which had obsremain > 0)\n",
    "    tg_obs = dict()\n",
    "    for tgid in all_tgs.ids():\n",
    "        tgprops = all_tgs.get(tgid)\n",
    "        tg_obs[tgid] = tgprops.obsremain\n",
    "\n",
    "    tree = TargetTree(all_tgs, 0.01)\n",
    "    tgsavail = TargetsAvailable(hw, all_tgs, tiles, tree)\n",
    "    del tree\n",
    "    favail = LocationsAvailable(tgsavail)\n",
    "    asgn = Assignment(all_tgs, tgsavail, favail)\n",
    "    asgn.assign_unused(TARGET_TYPE_SCIENCE)\n",
    "    \n",
    "    # Assign standards, up to some limit\n",
    "    asgn.assign_unused(TARGET_TYPE_STANDARD, 10)\n",
    "\n",
    "    # Assign sky to unused fibers, up to some limit\n",
    "    asgn.assign_unused(TARGET_TYPE_SKY, 40)\n",
    "\n",
    "    # Force assignment if needed\n",
    "    asgn.assign_force(TARGET_TYPE_STANDARD, 10)\n",
    "    asgn.assign_force(TARGET_TYPE_SKY, 40)\n",
    "    \n",
    "    print(\"Realization {:02d}:\".format(mc))\n",
    "    \n",
    "    # Accumulate results for all tiles\n",
    "    for tid in tiles.id:\n",
    "        # Assignment for this tile\n",
    "        tile_assign = asgn.tile_location_target(tid)\n",
    "\n",
    "        # Targets available for all tile / locs\n",
    "        tgsavail = asgn.targets_avail()\n",
    "\n",
    "        # available for this tile\n",
    "        tile_avail = tgsavail.tile_data(tid)\n",
    "        \n",
    "        # Accumulate results\n",
    "        accum_histogram_data(\n",
    "            tid, all_classes, tile_assign, tile_avail, target_id_class, \n",
    "            tg_obs, hist_tgassign, hist_tgavail, hist_tgconsid,\n",
    "            verbose=True\n",
    "        )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results\n",
    "\n",
    "plot_assignment_stats(\n",
    "    wdir,\n",
    "    \"sky-std\", \n",
    "    \"Science, Standards and Sky for {} Passes\".format(len(tiles.id)), \n",
    "    tiles.id, \n",
    "    all_classes, \n",
    "    hist_tgassign, \n",
    "    hist_tgavail, \n",
    "    hist_tgconsid\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When interpreting the above plots, first consider the \"fraction of assigned positioners\" for each tile.  The sharply peaked histogram for standards and sky is a result of us **forcing** 10 standards and 40 skies per petal.  When doing this \"bumping\" of science targets, we start with the lowest priority targets and this is why the overall number of ELGs drops for each tile compared to the previous section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "Hopefully this notebook has built some intuition about the limits of DESI positioner geometry and the assignment results that occur in the presence of multiple target classes, each with different priorities and densities, and observed over multiple passes / layers.  Part 2 of this study will bring in a \"real\" target sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "DESI 20.8",
   "language": "python",
   "name": "desi-20.8"
  },
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}