plasticc-eclipsing-binaries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5ee0a705-b529-45cc-afb9-382350e63fb6",
   "metadata": {},
   "source": [
    "# PLAsTiCC data exploration with TAPE\n",
    "\n",
    "Let's explore [PLAsTiCC](http://plasticc.org) data!\n",
    "\n",
    "It is publically avilable through [this Zenodo repository](https://zenodo.org/record/2539456)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a92caa3a0873700f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-19T14:50:37.291215Z",
     "start_time": "2023-10-19T14:50:37.289183Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Uncomment to install packages\n",
    "\n",
    "!pip install tape joblib requests pyarrow.csv pyarrow"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a8570b9-2a7b-4d79-be91-49a1d74d3f84",
   "metadata": {},
   "source": [
    "### Download the data\n",
    "\n",
    "Please get the data from [Zenodo](https://zenodo.org/record/2539456) and put them to the `./plasticc` folder (you may change the location bellow with `DATA_DIR`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4866d31cc059f9d7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-19T14:50:39.090278Z",
     "start_time": "2023-10-19T14:50:39.083605Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "DATA_DIR = Path(\"../../../../../datasets/plasticc\")\n",
    "\n",
    "META_FILENAME = \"plasticc_test_metadata.csv\"\n",
    "LC_FILENAMES = [f\"plasticc_test_lightcurves_{i:02d}.csv\" for i in range(1, 12)]\n",
    "\n",
    "N_PARTITIONS = len(LC_FILENAMES)\n",
    "\n",
    "N_PROCESSORS = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5cea8a584d359ad",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "You may skip the next cell if you already have the data downloaded in the `DATA_DIR`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61f1f943627a85da",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Read and analyse the data with TAPE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e023704e2d93eb13",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-18T21:26:16.650747Z",
     "start_time": "2023-10-18T21:24:42.075691Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "import light_curve as licu\n",
    "import dask.dataframe as dd\n",
    "import numpy as np\n",
    "from tape import Ensemble, ColumnMapper, TapeFrame\n",
    "\n",
    "# In TAPE's (and LSST's) terminology, sources are individual detections,\n",
    "# and objects are the underlying astrophysical objects.\n",
    "\n",
    "# We load object table first, from the metadata file.\n",
    "print(\"Loading object table...\")\n",
    "object_table = dd.read_csv(\n",
    "    DATA_DIR / META_FILENAME,\n",
    "    # Read data chunk by chunk, to avoid loading the whole file into memory.\n",
    "    blocksize=100e6,\n",
    ")\n",
    "\n",
    "# Then we load the sources:\n",
    "print(\"Loading source tables...\")\n",
    "source_table = dd.read_csv(\n",
    "    [DATA_DIR / filename for filename in LC_FILENAMES],\n",
    "    blocksize=100e6,\n",
    ")\n",
    "\n",
    "# Now we can make an Ensemble.\n",
    "# To make parallel processing work we need to partition the data.\n",
    "# After that, when we run analysis, TAPE will distribute the work:\n",
    "# one partition per worker.\n",
    "print(\"Building Ensemble...\")\n",
    "ens = Ensemble(\n",
    "    n_workers=N_PROCESSORS,\n",
    ")\n",
    "ens.from_dask_dataframe(\n",
    "    source_frame=source_table,\n",
    "    object_frame=object_table,\n",
    "    npartitions=None,\n",
    "    column_mapper=ColumnMapper(\n",
    "        id_col='object_id',\n",
    "        time_col='mjd',\n",
    "        flux_col='flux',\n",
    "        err_col='flux_err',\n",
    "        band_col='passband',\n",
    "    ),\n",
    "    # TODO: Doug, please put a comment here about the following line.\n",
    "    sync_tables=False,\n",
    "    sorted=True,\n",
    "    sort=False,\n",
    ")\n",
    "\n",
    "# Let's run some analysis!\n",
    "\n",
    "print(\"Starting analysis...\")\n",
    "# First, let's select only Galactic objects, by cutting on hostgal_photoz.\n",
    "print(\"First, filter by photoz\")\n",
    "ens = ens.query(\"hostgal_photoz < 1e-3\", table=\"object\")\n",
    "\n",
    "print(\"Extract durations\")\n",
    "duration = ens.batch(\n",
    "    lambda time, detected: np.ptp(time[np.asarray(detected, dtype=bool)]),\n",
    "    ens._time_col, 'detected_bool',\n",
    "    meta=('duration', 'float64'),\n",
    "    use_map=False,\n",
    "    compute=False,\n",
    ")\n",
    "\n",
    "print(\"Assign a column\")\n",
    "ens.assign(table=\"object\", duration=duration)\n",
    "print(\"Filter by duration\")\n",
    "ens = ens.query(\"duration > 366\", table=\"object\")\n",
    "\n",
    "# Next, we use Otsu's method to split light curves into two groups:\n",
    "# one with high flux, and one with low flux. Eclipsing binaries should have\n",
    "# lower flux group smaller than the higher flux group, but having larger \n",
    "# variability. We use light-curve package to extract these features.\n",
    "# (https://github.com/light-curve/light-curve-python)\n",
    "# For simplicity, we only calculate these features for the i (3) band.\n",
    "print(\"Extract Otsu features\")\n",
    "otsu_features = ens.batch(licu.OtsuSplit(), band_to_calc=3, use_map=False, compute=False)\n",
    "print(\"Assign columns\")\n",
    "ens = ens.assign(\n",
    "    table=\"object\",\n",
    "    otsu_lower_to_all_ratio=otsu_features['otsu_lower_to_all_ratio'],\n",
    "    otsu_std_lower=otsu_features['otsu_std_lower'],\n",
    "    otsu_std_upper=otsu_features['otsu_std_upper'],\n",
    ")\n",
    "\n",
    "print('Filter by Otsu features')\n",
    "ens = ens.query(\n",
    "    \"otsu_lower_to_all_ratio < 0.1 and otsu_std_lower > otsu_std_upper\",\n",
    "    table=\"object\",\n",
    ")\n",
    "print(\"Compute object table\")\n",
    "ens.compute(\"object\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afd58005",
   "metadata": {},
   "source": [
    "# Example with the Refactored TAPE API"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e71b40fc",
   "metadata": {},
   "source": [
    "The TAPE refactor currently is on branch [tape_ensemble_refactor](https://github.com/lincc-frameworks/tape/tree/tape_ensemble_refactor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4712ca61",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#%%memit\n",
    "\n",
    "import light_curve as licu\n",
    "import dask.dataframe as dd\n",
    "import numpy as np\n",
    "from tape import Ensemble, ColumnMapper\n",
    "\n",
    "# In TAPE's (and LSST's) terminology, sources are individual detections,\n",
    "# and objects are the underlying astrophysical objects.\n",
    "\n",
    "# We load object table first, from the metadata file.\n",
    "print(\"Loading object table...\")\n",
    "object_table = dd.read_csv(\n",
    "    DATA_DIR / META_FILENAME,\n",
    "    # Read data chunk by chunk, to avoid loading the whole file into memory.\n",
    "    blocksize=100e6,\n",
    ")\n",
    "\n",
    "# Then we load the sources:\n",
    "print(\"Loading source tables...\")\n",
    "source_table = dd.read_csv(\n",
    "    [DATA_DIR / filename for filename in LC_FILENAMES],\n",
    "    blocksize=100e6,\n",
    ")\n",
    "\n",
    "# Now we can make an Ensemble.\n",
    "# To make parallel processing work we need to partition the data.\n",
    "# After that, when we run analysis, TAPE will distribute the work:\n",
    "# one partition per worker.\n",
    "print(\"Building Ensemble...\")\n",
    "ens = Ensemble(\n",
    "    n_workers=N_PROCESSORS,\n",
    ")\n",
    "ens.from_dask_dataframe(\n",
    "    source_frame=source_table,\n",
    "    object_frame=object_table,\n",
    "    npartitions=None,\n",
    "    column_mapper=ColumnMapper(\n",
    "        id_col='object_id',\n",
    "        time_col='mjd',\n",
    "        flux_col='flux',\n",
    "        err_col='flux_err',\n",
    "        band_col='passband',\n",
    "    ),\n",
    "    # TODO: Doug, please put a comment here about the following line.\n",
    "    sync_tables=False,\n",
    "    sorted=True,\n",
    "    sort=False,\n",
    ")\n",
    "\n",
    "# Let's run some analysis!\n",
    "\n",
    "print(\"Starting analysis...\")\n",
    "# First, let's select only Galactic objects, by cutting on hostgal_photoz.\n",
    "print(\"First, filter by photoz\")\n",
    "ens = ens._object.query(\"hostgal_photoz < 1e-3\").update_ensemble()\n",
    "\n",
    "# Second, let's select persistent sources, by cutting on the duration of the light curve.\n",
    "print(\"Extract durations\")\n",
    "ens.batch(\n",
    "    lambda time, detected: np.ptp(time[np.asarray(detected, dtype=bool)]),\n",
    "    ens._time_col, 'detected_bool',\n",
    "    meta=('duration', 'float64'),\n",
    "    use_map=False,\n",
    "    compute=False,\n",
    "    label=\"duration\"\n",
    ")\n",
    "\n",
    "print(\"Filter by duration\")\n",
    "ens = ens._object.assign(duration=ens.select_frame(\"duration\")).query(\"duration > 366\").update_ensemble()\n",
    "\n",
    "# Next, we use Otsu's method to split light curves into two groups:\n",
    "# one with high flux, and one with low flux. Eclipsing binaries should have\n",
    "# lower flux group smaller than the higher flux group, but having larger \n",
    "# variability. We use light-curve package to extract these features.\n",
    "# (https://github.com/light-curve/light-curve-python)\n",
    "# For simplicity, we only calculate these features for the i (3) band.\n",
    "print(\"Extract Otsu features\")\n",
    "ens.batch(licu.OtsuSplit(), band_to_calc=3, use_map=False, compute=False, label=\"otsu_features\")\n",
    "\n",
    "print(\"Merge and filter by frame otsu_features, dropping unused 'otsu_mean_diff'\")\n",
    "ens = ens._object.merge(ens.select_frame(\"otsu_features\").drop(columns=\"otsu_mean_diff\")).query(\n",
    "    \"otsu_lower_to_all_ratio < 0.1 and otsu_std_lower > otsu_std_upper\",\n",
    ").update_ensemble()\n",
    "\n",
    "print(\"Compute object table\")\n",
    "ens.compute(\"object\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41476677",
   "metadata": {},
   "source": [
    "### Do the same, but with Pandas + joblib\n",
    "\n",
    "The workflow is challenger here: we don't like to load all the light curves into memory (~20 GBi), so we have to manually jump between object table processing and source table processing.\n",
    "\n",
    "NB: We also need `pyarrow.csv`, but only to read CSVs faster into pandas data-frames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50cd141e7891a5fe",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-19T14:51:47.541237Z",
     "start_time": "2023-10-19T14:50:43.446726Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "import light_curve as licu\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pyarrow\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "# Set number of threads to use when reading CSVs\n",
    "pyarrow.set_cpu_count(N_PROCESSORS)\n",
    "pyarrow.set_io_thread_count(N_PROCESSORS)\n",
    "\n",
    "# First we load object table, from the metadata file.\n",
    "print(\"Loading object table...\")\n",
    "object_table = pd.read_csv(DATA_DIR / META_FILENAME, engine='pyarrow', index_col='object_id')\n",
    "\n",
    "print(\"Pre-filter object table\")\n",
    "# Before doing light-curve analysis we select only Galactic objects, by cutting on hostgal_photoz.\n",
    "object_table = object_table[object_table['hostgal_photoz'] < 1e-3]\n",
    "\n",
    "\n",
    "def read_source_table(filename):\n",
    "    # We are going to read in parallel, so we set single thread \n",
    "    pyarrow.set_cpu_count(1)\n",
    "    pyarrow.set_io_thread_count(1)\n",
    "    return pd.read_csv(filename, engine='pyarrow', index_col='object_id')\n",
    "\n",
    "\n",
    "def process_source_table(filename):\n",
    "    source_table = read_source_table(filename)\n",
    "    \n",
    "    # Select sources for objects left\n",
    "    source_table = source_table[source_table.index.isin(object_table.index)]\n",
    "    \n",
    "    # Select persistent sources, by cutting on the duration of the light curve.\n",
    "    detections = source_table[source_table['detected_bool'].astype(bool)]\n",
    "    # We need to group detections by object_id, to get light curves:\n",
    "    durations = detections['mjd'].groupby(level=0).apply(np.ptp)\n",
    "    \n",
    "    # Then we filter the source table to keep only persistent sources:\n",
    "    source_table = source_table.loc[durations[durations > 366].index]\n",
    "    \n",
    "    # Next, we use Otsu's method to split light curves into two groups:\n",
    "    # one with high flux, and one with low flux. Eclipsing binaries should have\n",
    "    # lower flux group significantly smaller than the higher flux group,\n",
    "    # but having larger variability.\n",
    "    # We use light-curve package to extract these features.\n",
    "    # (https://github.com/light-curve/light-curve-python)\n",
    "    # For simplicity, we only calculate these features for the i (3) band.\n",
    "    \n",
    "    otsu_split = licu.OtsuSplit()\n",
    "    # i_band_sources = source_table[source_table['passband'] == 3]\n",
    "    otsu_features = source_table[['mjd', 'flux', 'flux_err', 'passband']].groupby(level=0).apply(lambda row: pd.Series(\n",
    "        otsu_split(\n",
    "            np.asarray(row['mjd']),\n",
    "            np.asarray(row['flux']),\n",
    "            np.asarray(row['flux_err']),\n",
    "            sorted=True,\n",
    "            check=False,\n",
    "        ),\n",
    "        index=otsu_split.names,\n",
    "    ))\n",
    "    \n",
    "    # And filter by Otsu features\n",
    "    otsu_idx = (otsu_features['otsu_lower_to_all_ratio'] < 0.1) & (otsu_features['otsu_std_lower'] > otsu_features['otsu_std_upper'])\n",
    "    \n",
    "    return source_table.loc[otsu_idx.index[otsu_idx]]\n",
    "\n",
    "\n",
    "# Running the main part of the analysis in parallel\n",
    "print(\"Run analysis in parallel...\")\n",
    "source_tables = Parallel(n_jobs=N_PROCESSORS)(\n",
    "    delayed(process_source_table)(DATA_DIR / filename) for filename in LC_FILENAMES\n",
    ")\n",
    "source_table = pd.concat(source_tables, ignore_index=False, sort=True)\n",
    "object_table = object_table.loc[np.unique(source_table.index)]\n",
    "object_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17869fda5032e67e",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# TODO\n",
    "### Redo without having everything in-memory\n",
    "### Add polars example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c92b01b629070426",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### Do the same, but with Pandas+PyArrow w/ nested arrays for dataframes, and joblib for parallelism"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e215607e9fd17fd",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-18T19:07:04.808063Z",
     "start_time": "2023-10-18T19:04:46.809682Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "import light_curve as licu\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pyarrow as pa\n",
    "import pyarrow.csv as pacsv\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "# Read the data\n",
    "# -------------\n",
    "\n",
    "# First we load object table, from the metadata file.\n",
    "print(\"Loading object table...\")\n",
    "object_table = pacsv.read_csv(\n",
    "    DATA_DIR / META_FILENAME,\n",
    "    # We'd like to load the whole file into a single partition\n",
    "    read_options=pacsv.ReadOptions(block_size=(1<<31)-1),\n",
    ")\n",
    "object_table = pd.DataFrame(\n",
    "    {\n",
    "        col: pd.Series(\n",
    "            object_table[col],\n",
    "            dtype=pd.ArrowDtype(object_table[col].type),\n",
    "            index=object_table['object_id'],\n",
    "            copy=False,\n",
    "        )\n",
    "        for col in object_table.column_names if col != 'object_id'\n",
    "    },\n",
    ")\n",
    "\n",
    "def read_source_table(filename):\n",
    "    table = pacsv.read_csv(\n",
    "        DATA_DIR / filename,\n",
    "        # We'd like to have a partition per an original file\n",
    "        read_options=pacsv.ReadOptions(block_size=(1<<31)-1),\n",
    "    )\n",
    "    return pd.DataFrame(\n",
    "        {\n",
    "            col: pd.Series(\n",
    "                table[col],\n",
    "                dtype=pd.ArrowDtype(table[col].type),\n",
    "                index=table['object_id'],\n",
    "                copy=False,\n",
    "            )\n",
    "            for col in table.column_names if col != 'object_id'\n",
    "        },\n",
    "    )\n",
    "\n",
    "# source_tables = []\n",
    "# for filename in LC_FILENAMES:\n",
    "#     source_table = pa.csv.read_csv(\n",
    "#         DATA_DIR / filename,\n",
    "#         # We'd like to have a partition per an original file\n",
    "#         read_options=pa.csv.ReadOptions(block_size=(1<<31)-1),\n",
    "#     )\n",
    "#     source_tables.append(pd.DataFrame(\n",
    "#         {\n",
    "#             col: pd.Series(\n",
    "#                 source_table[col],\n",
    "#                 dtype=pd.ArrowDtype(source_table[col].type),\n",
    "#                 index=source_table['object_id'],\n",
    "#                 copy=False,\n",
    "#             )\n",
    "#             for col in source_table.column_names if col != 'object_id'\n",
    "#         },\n",
    "#     ))\n",
    "source_tables = Parallel(backend='threading', n_jobs=N_PROCESSORS)(\n",
    "    delayed(read_source_table)(filename) for filename in LC_FILENAMES\n",
    ")\n",
    "source_table = pd.concat(source_tables, ignore_index=True, sort=False)\n",
    "\n",
    "\n",
    "# Add sources to the object table\n",
    "# -------------------------------\n",
    "\n",
    "# First, let's do some sanity checks\n",
    "print(\"Sanity checks...\")\n",
    "np.testing.assert_array_equal(\n",
    "    object_table.index.values,\n",
    "    np.unique(object_table.index.values),\n",
    "    err_msg=\"Object table has duplicate indices or is not sorted.\",\n",
    ")\n",
    "assert np.all(np.diff(source_table.index) >= 0), \"Source table index must be sorted.\"\n",
    "\n",
    "# We need an offsets array to know where each source light curve starts.\n",
    "source_offsets = []\n",
    "for table in source_tables:\n",
    "    offset = np.nonzero(np.diff(table.index, prepend=-1, append=-1))[0]\n",
    "    source_offsets.append(pa.array(offset))\n",
    "\n",
    "# Update the object table with list-arrays built from the source table\n",
    "\n",
    "print(\"Updating object table with list-arrays...\")\n",
    "for column in source_table.columns:\n",
    "    list_arrays = []\n",
    "    for table, offset in zip(source_tables, source_offsets):\n",
    "        list_arrays.append(pa.ListArray.from_arrays(\n",
    "            offset,\n",
    "            pa.array(table[column]),\n",
    "        ))\n",
    "    chunked_array = pa.chunked_array(list_arrays)\n",
    "    object_table[column] = pd.Series(\n",
    "        chunked_array,\n",
    "        dtype=pd.ArrowDtype(chunked_array.type),\n",
    "        index=object_table.index,\n",
    "    )\n",
    "    \n",
    "# Do analysis\n",
    "# -----------\n",
    "\n",
    "print(\"Starting analysis...\")\n",
    "# First, let's select only Galactic objects, by cutting on hostgal_photoz.\n",
    "df = object_table[object_table['hostgal_photoz'] < 1e-3]\n",
    "\n",
    "# Second, let's select persistent sources, by cutting on the duration of the light curve.\n",
    "df['duration'] = df[['mjd', 'detected_bool']].apply(\n",
    "    lambda row: np.ptp(row['mjd'][np.asarray(row['detected_bool'], dtype=bool)]),\n",
    "    axis=1\n",
    ")\n",
    "df = df[df['duration'] > 366]\n",
    "\n",
    "# Next, we use Otsu's method to split light curves into two groups:\n",
    "# one with high flux, and one with low flux. Eclipsing binaries should have\n",
    "# lower flux group significantly smaller than the higher flux group,\n",
    "# but having larger variability.\n",
    "# We use light-curve package to extract these features.\n",
    "# (https://github.com/light-curve/light-curve-python)\n",
    "# For simplicity, we only calculate these features for the i band.\n",
    "def extract_band(*arrays, bands, band_to_calc):\n",
    "    mask = np.asarray(bands) == band_to_calc\n",
    "    return [np.asarray(arr)[mask] for arr in arrays]\n",
    "\n",
    "otsu_split = licu.OtsuSplit()\n",
    "otsu_features = df[['mjd', 'flux', 'flux_err', 'passband']].apply(\n",
    "    lambda row: pd.Series(\n",
    "        otsu_split(\n",
    "            *extract_band(\n",
    "                row['mjd'],\n",
    "                row['flux'],\n",
    "                row['flux_err'],\n",
    "                bands=row['passband'],\n",
    "                band_to_calc=3,  # i band\n",
    "            ),\n",
    "            sorted=True,\n",
    "            check=False,\n",
    "        ),\n",
    "        index=otsu_split.names,\n",
    "    ),\n",
    "    axis=1\n",
    ")\n",
    "\n",
    "df = df[otsu_features['otsu_lower_to_all_ratio'] < 0.1]\n",
    "df = df[otsu_features['otsu_std_lower'] > otsu_features['otsu_std_upper']]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b336d01b8f74587",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-18T19:07:13.840204Z",
     "start_time": "2023-10-18T19:07:04.807562Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "LSST_BANDS = 'ugrizy'\n",
    "\n",
    "def plot(row):\n",
    "    plt.figure()\n",
    "    plt.title(f\"Object {row.Index}, true class {row.true_target}\")\n",
    "    plt.xlabel('MJD')\n",
    "    plt.ylabel('Flux, zp=27.5')\n",
    "    for band_idx, band_name in enumerate(LSST_BANDS):\n",
    "        mjd, flux, flux_err = extract_band(\n",
    "            row.mjd,\n",
    "            row.flux,\n",
    "            row.flux_err,\n",
    "            bands=row.passband,\n",
    "            band_to_calc=band_idx,\n",
    "        )\n",
    "        color = f'C{band_idx}'\n",
    "        plt.scatter(mjd, flux, c=color, label=band_name)\n",
    "        plt.errorbar(mjd, flux, yerr=flux_err, ls='none', c=color)\n",
    "        plt.legend()\n",
    "\n",
    "# Random objects from the selected sample\n",
    "for row in object_table[object_table['true_target'] == 16].sample(5, random_state=0).itertuples():\n",
    "    plot(row)\n",
    "    \n",
    "# Random objects from the selected sample\n",
    "for row in df.sample(5, random_state=0).itertuples():\n",
    "    plot(row)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}