variable_stars_starting_kit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"text-align: center;\">\n",
    "    <a href=\"http://www.datascience-paris-saclay.fr\">\n",
    "        <img border=\"0\" src=\"http://project.inria.fr/saclaycds/files/2017/02/logoUPSayPlusCDS_990.png\" width=\"90%\"></a>\n",
    "</div>\n",
    "\n",
    "# RAMP on variable star type prediction\n",
    "\n",
    "<i> Balázs Kégl (LAL/CNRS), Marc Moniez (LAL/CNRS), Alex Gramfort (Inria), Djalel Benbouzid (UPMC), Mehdi Cherti (LAL/CNRS) </i>\n",
    "\n",
    "This project was downloaded from https://github.com/ramp-kits/variable_stars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "### Variable stars\n",
    "\n",
    "Most stars emit light steadily in time, but a small fraction of them has a variable <i>light curve</i>: light emission versus time. We call them <a href = \"http://en.wikipedia.org/wiki/Variable_star\">variable stars</a>. The light curves are usually <a href=\"http://en.wikipedia.org/wiki/Periodic_function\">periodic</a> and highly regular. There are essentially two reasons why light emission can vary. First, the star itself can be <a href=\"https://www.youtube.com/watch?v=sXJBrRmHPj8\">oscillating</a>, so its light emission varies in time. Second, the star that seems a single point at Earth (because of our large distance) is actually a <a href=\"http://en.wikipedia.org/wiki/Binary_system\"><i>binary system</i></a>: two stars that <a href=\"https://www.youtube.com/watch?v=1kFFwHkxBiI\">orbit around their common center of gravity</a>. When the <a href=\"http://en.wikipedia.org/wiki/Orbital_plane_%28astronomy%29\">orbital plane</a> is parallel to our line of view, the stars <a href=\"http://csep10.phys.utk.edu/astr162/lect/binaries/eclipsing.html\">eclipse each other</a> periodically, creating a light curve with a <a href=\"http://www.dlr.de/en/desktopdefault.aspx/tabid-5170/8702_read-20474/\">charateristic signature</a>. Identifying, classifying, and analyzing variable stars are hugely important for calibrating distances, and making these analyses automatic will be crucial in the upcoming sky survey projects such as <a href=\"http://www.lsst.org/lsst/\">LSST</a>.\n",
    "\n",
    "### The EROS1 database, a catalog of light curves \n",
    "\n",
    "Our data is coming from the <a href=\"http://eros.in2p3.fr/\">EROS1 project</a> that took data between 1990 and 1994.\n",
    "400 <a href=\"http://www.atnf.csiro.au/outreach/education/senior/astrophysics/photometry_photographicastro.html\">photographic plates</a> were taken of the <a href=\"http://en.wikipedia.org/wiki/Large_Magellanic_Cloud\">Large Magellanic Cloud</a> with the <a href=\"http://www.eso.org/public/teles-instr/lasilla/1mschmidt/\">ESO-Schmidt 1m telescope</a>. Each plate covers the same $5^\\circ \\times 5^\\circ$ field, centered at $\\alpha  = 5h18m43s$, $\\delta = −69d42m17s$ in the <a href=\"http://en.wikipedia.org/wiki/Celestial_coordinate_system\">celestial coordinate system</a>. Photos were taken in two <a href = \"http://science-edu.larc.nasa.gov/EDDOCS/Wavelengths_for_Colors.html\">frequency bands</a> (red: 630nm; blue: 385nm), digitized at the <a href=\"http://en.wikipedia.org/wiki/Paris_Observatory\">Observatoire de Paris</a>, and analyzed at the <a href=\"http://cc.in2p3.fr/\">IN2P3 Computing Center</a>.\n",
    "\n",
    "### Selecting variable stars\n",
    "\n",
    "The full catalog contains 8 million objects. We estimated the probability of a star being stable using an in-house algorithm and selected the $\\simeq 1\\%$ least stable stars. The light curve of each selected star was then visually inspected, and the star was either declared stable or assigned a variability type and a quality index. The variability types are <a href=\"http://en.wikipedia.org/wiki/Binary_star#Eclipsing_binaries\">eclipsing binary</a>, <a href=\"http://en.wikipedia.org/wiki/Cepheid_variable\">Cepheid</a>, <a href=\"http://en.wikipedia.org/wiki/RR_Lyrae_variable\">RR-Lyrae</a>, <a href=\"http://en.wikipedia.org/wiki/Mira_variable\">Mira</a>, and other (the variability is clearly established, but the type is unclear). The quality index was visually estimated from 1 (lowest signal to background ratio) to 3 (best signal to background ratio). A total of 22802 variable objects were found that include 9046 RR-Lyrae, 2758 Cepheids, 1596 eclipsing binaries, 890 Miras, and 8512 unclassified objects. About 10% of the data was lost in the data archeology step that included converting the measurements from a native <a href=\"http://en.wikipedia.org/wiki/Physics_Analysis_Workstation\">PAW</a> format to csv, giving us a total of 19429 stars.\n",
    "\n",
    "### Selecting data\n",
    "\n",
    "For the RAMP, we decided to drop all instances in the \"unclassified object\" category since the interpretation of this type was unclear. We then randomly selected 30% of the data for training, giving us 3641 training instances.\n",
    "\n",
    "The data consists of two files, <a href=\"https://drive.google.com/file/d/0BzwKr6zuOkdRRWt5bnYyc2JNS0U/view?usp=sharing\">train.csv</a> contains \"static\" features in a classical row-wise csv table, and <a href=\"https://drive.google.com/file/d/0BzwKr6zuOkdRMXVuRE9FSXNjZTQ/view?usp=sharing\">train_varlength_features.csv.gz</a> is a table that contains all the time series. The contents is obtained with `get_train_data` function below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploratory data analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from itertools import chain\n",
    "pd.set_option('display.max_columns', None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get access to the training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import problem\n",
    "\n",
    "X_df, y = problem.get_train_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The static features\n",
    "\n",
    "<span style=\"color:red\">Red variables</span> are ids, labels, or other human-annotated features, so they should not be used as input in the classification. <span style=\"color:lightblue\">Light blue variables</span> are legal but not likely to contribute information to the classification. Black variables are definitely discriminative.\n",
    "\n",
    "<ul>\n",
    "    <li> <code><b><span style=\"color:red\">patch_id</span></b></code>: The catalogue is organized by tiles corresponding to $1$cm$^2$ patches on the plates, this is their id.\n",
    "    <li> <code><b><span style=\"color:red\">star_id_b, star_id_r</span></b></code>: The id of the star within the patch. <code>patch_id</code> and <code>star_id_b</code> or <code>star_id_r</code> identify the stars uniquely, and we use <code>[patch_id]_[star_id_b]</code> for indexing the stars in the pandas table.\n",
    "    <li> <code><b>magnitude_b, magnitude_r</b></code>: The average apparent luminosity of the star (in two frequency bands). <a href=\"http://www.astro-tom.com/technical_data/magnitude_scale.htm\">Magnitude</a> is a logarithmic measure, and the higher it is, the lower the apparent luminosity is.\n",
    "    <li> <code><b><span style=\"color:lightblue\">asc_d, asc_m, asc_s</span></b></code>: <a href=\"http://en.wikipedia.org/wiki/Celestial_coordinate_system\">Celestial</a> <a href=\"http://en.wikipedia.org/wiki/Right_ascension\">right ascension</a> (coordinate) of the star, measured in degrees, minutes, and seconds, respectively.\n",
    "    <li> <code><b><span style=\"color:lightblue\">dec_d, dec_m, dec_s</span></b></code>: <a href=\"http://en.wikipedia.org/wiki/Celestial_coordinate_system\">Celestial</a> <a href=\"http://en.wikipedia.org/wiki/Declination\">declination</a> (coordinate) of the star, measured in degrees, minutes, and seconds, respectively.\n",
    "    <li> <code><b><span style=\"color:black\">period</span></b></code>: the estimated period of the light curve. For the correct period, it should be divided by <code>div_period</code>\n",
    "    <li> <code><b><span style=\"color:lightblue\">frequency</span></b></code>: $1/$<code>period</code>, so it is redundant.\n",
    "    <li> <code><b><span style=\"color:lightblue\">num_points_good_b, num_points_good_r</span></b></code>: Number of good light curve measurements (some measurements can be corrupted). \n",
    "    <li> <code><b><span style=\"color:black\">asym_b, asym_r</span></b></code>: Unknown semantics.\n",
    "    <li> <code><b><span style=\"color:black\">log_p_not_variable</span></b></code>: Logarithm of the estimated probability that the star is stable.\n",
    "    <li> <code><b><span style=\"color:black\">sigma_flux_b, sigma_flux_r</span></b></code>: The square root of the total variance of the light measurements (indicating the amplitude of the variability).\n",
    "    <li> <code><b><span style=\"color:red\">type</span></b></code>: The label to predict.\n",
    "    <li> <code><b><span style=\"color:red \">quality</span></b></code>: Human-annotated measure of the quality of the time curve. The higher the better.\n",
    "    <li> <code><b><span style=\"color:black \">div_period</span></b></code>: The algorithm that estimates the <code>period</code> sometimes finds a multiple of the period. These cases were human-detected, and the divisor was recorded. In principle this variable is not available automatically as an observable, but we are confident that it could be obtained automatically, so we allow it as an input.\n",
    "</ul> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "label_names = {1: 'binary', 2: 'cepheid', 3: 'rr_lyrae', 4: 'mira'}\n",
    "labels = list(label_names.keys())\n",
    "y_series = pd.Series(y).replace(label_names)\n",
    "y_series.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = y_series.value_counts().plot(kind=\"bar\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some classwise histograms and scatterplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['r', 'b', 'g', 'm']\n",
    "def plot_classwise_normalized(feature, bins=None):\n",
    "    if bins is None:\n",
    "        bins = np.linspace(X_df[feature].min(), X_df[feature].max(), 15)\n",
    "    for label, color in zip(labels, colors):\n",
    "        plt.hist(X_df[y == label][feature].values, density=True, bins=bins, \n",
    "                 alpha=0.5, color=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('period')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('period', bins=np.linspace(0, 50, 15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe the aliasing below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('period', bins=np.linspace(0, 5, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_df['real_period'] = X_df['period'] / X_df['div_period']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('real_period')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('real_period', bins=np.linspace(0, 50, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('real_period', bins=np.linspace(0, 5, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('magnitude_b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('magnitude_r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('asym_b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('asym_r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('sigma_flux_b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('sigma_flux_b', bins=np.linspace(0, 1000, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('sigma_flux_b', bins=np.linspace(0, 100, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('sigma_flux_r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('sigma_flux_r', bins=np.linspace(0, 1000, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('sigma_flux_r', bins=np.linspace(0, 100, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_normalized('log_p_not_variable')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['r', 'b', 'g', 'm']\n",
    "def plot_classwise_scatter(feature1, feature2, range1=None, range2=None):\n",
    "    if range1 is None:\n",
    "        range1 = [X_df[feature1].min(), X_df[feature1].max()]\n",
    "    if range2 is None:\n",
    "        range2 = [X_df[feature2].min(), X_df[feature2].max()]\n",
    "    for label, color in zip(labels, colors):\n",
    "        plt.xlim(range1[0], range1[1])\n",
    "        plt.ylim(range2[0], range2[1])\n",
    "        plt.scatter(X_df[y == label][feature1], \n",
    "                    X_df[y == label][feature2],\n",
    "            alpha=0.3, s=80, c=color, marker='.');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_scatter('magnitude_b', 'magnitude_r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_classwise_scatter('magnitude_b', 'real_period', range1=None, range2=[0,10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The time series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each column contains a list of floating point numbers.\n",
    "\n",
    "<ul>\n",
    "    <li> <code><b><span style=\"color:black\">time_points_b, time_points_r</span></b></code>: The time (in unit of days) when the photos were taken. Note that the filters had to be changed so the time points of the blue and red frequency band are slighty different.\n",
    "    <li> <code><b><span style=\"color:black\">light_points_b, light_points_r</span></b></code>: The light points measured at the time points.\n",
    "    <li> <code><b><span style=\"color:black\">error_points_b, error_points_r</span></b></code>: Uncertainties (error bars) on the light measurements.\n",
    "    <li> <code><b><span style=\"color:lightblue\">bkg_points_, bkg_points_r</span></b></code>: Background noise measured at the time points.\n",
    "    <li> <code><b><span style=\"color:lightblue\">polltn_points_b, polltn_points_r</span></b></code>: Pollution noise measured at the time points.\n",
    "</ul> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_df.hist('num_points_good_r')\n",
    "print(min(X_df['num_points_good_b']))\n",
    "print(max(X_df['num_points_good_b']))\n",
    "print(min(X_df['num_points_good_r']))\n",
    "print(max(X_df['num_points_good_r']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting time curves\n",
    "\n",
    "Set the patch id and star id below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patch_id = 98\n",
    "star_id_b = 477\n",
    "\n",
    "def star_key(slab_id, star_id_b):\n",
    "    return str(slab_id) + '_' + str(star_id_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_df = X_df.set_index(X_df.apply(lambda row: \"%d_%d\" % (row['patch_id'], row['star_id_b']), axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_points = X_df.loc[star_key(patch_id, star_id_b)]['time_points_b']\n",
    "light_points = X_df.loc[star_key(patch_id, star_id_b)]['light_points_b']\n",
    "error_points = X_df.loc[star_key(patch_id, star_id_b)]['error_points_b']\n",
    "plt.errorbar(time_points, light_points, yerr=error_points, fmt='o');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The raw measurements seem rather messy. The scatter of the plots is visibly larger than the measurement uncertainty (which makes it, by definition, a variable star), but there is no visible periodicity. We can use the estimated period to overplot several periods of the curve (\"fold\" the time series) using the following function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fold_time_series(time_point, period, div_period):\n",
    "    real_period = period / div_period\n",
    "    return time_point % real_period  # modulo real_period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period = X_df.loc[star_key(patch_id, star_id_b)]['period']\n",
    "div_period = X_df.loc[star_key(patch_id, star_id_b)]['div_period']\n",
    "print(period, div_period)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_points_folded = [fold_time_series(time_point, period, div_period) \n",
    "                      for time_point in time_points]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting curve has a characteristic signature. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.gca().invert_yaxis()\n",
    "plt.errorbar(time_points_folded, light_points, yerr=error_points, fmt='o');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise\n",
    "\n",
    "The goal of the exercise is to classify the stars into one of the four types. In your code you will have access both to the static features and the time series. The submission site will have several examples that you can start from."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The feature extractor\n",
    "\n",
    "The input data are stored in a dataframe. To go from a dataframe to a numpy array we will\n",
    "use a scikit-learn column transformer. The first example we will write will just consist in\n",
    "selecting a subset of columns we want to work with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = [\n",
    "    'magnitude_b', \n",
    "    'magnitude_r',\n",
    "    'period',\n",
    "    'asym_b', \n",
    "    'asym_r', \n",
    "    'log_p_not_variable', \n",
    "    'sigma_flux_b', \n",
    "    'sigma_flux_r', \n",
    "    'quality', \n",
    "    'div_period',\n",
    "]\n",
    "\n",
    "from sklearn.compose import make_column_transformer\n",
    "\n",
    "transformer = make_column_transformer(\n",
    "    ('passthrough', cols)\n",
    ")\n",
    "\n",
    "X_array = transformer.fit_transform(X_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at how to transform line curves.\n",
    "\n",
    "The following feature extractor takes the light curve, bins it into <code>num_bins</code> bins, and return the bin means. It works with one band at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fold_time_series(time_point, period, div_period):\n",
    "    return (time_point -\n",
    "            (time_point // (period / div_period)) * period / div_period)\n",
    "\n",
    "\n",
    "def get_bin_means(X_df, num_bins, band):\n",
    "    feature_array = np.empty((len(X_df), num_bins))\n",
    "\n",
    "    for k, (_, x) in enumerate(X_df.iterrows()):\n",
    "        period = x['period']\n",
    "        div_period = x['div_period']\n",
    "        real_period = period / div_period\n",
    "        bins = [i * real_period / num_bins for i in range(num_bins + 1)]\n",
    "\n",
    "        time_points = np.array(x['time_points_' + band])\n",
    "        light_points = np.array(x['light_points_' + band])\n",
    "        time_points_folded = \\\n",
    "            np.array([fold_time_series(time_point, period, div_period)\n",
    "                      for time_point in time_points])\n",
    "        time_points_folded_digitized = \\\n",
    "            np.digitize(time_points_folded, bins) - 1\n",
    "\n",
    "        for i in range(num_bins):\n",
    "            this_light_points = light_points[time_points_folded_digitized == i]\n",
    "            if len(this_light_points) > 0:\n",
    "                feature_array[k, i] = np.mean(this_light_points)\n",
    "            else:\n",
    "                feature_array[k, i] = np.nan  # missing\n",
    "\n",
    "    return feature_array\n",
    "\n",
    "\n",
    "get_bin_means(X_df.iloc[:2], 5, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this we will use a funtion transformer that will get applied to both red and blues curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import FunctionTransformer\n",
    "\n",
    "transformer_r = FunctionTransformer(\n",
    "    lambda X_df: get_bin_means(X_df, 5, 'r')\n",
    ")\n",
    "\n",
    "transformer_b = FunctionTransformer(\n",
    "    lambda X_df: get_bin_means(X_df, 5, 'b')\n",
    ")\n",
    "\n",
    "transformer = make_column_transformer(\n",
    "    (transformer_r, ['period', 'div_period', 'time_points_r', 'light_points_r']),\n",
    "    (transformer_b, ['period', 'div_period', 'time_points_b', 'light_points_b']),\n",
    ")\n",
    "\n",
    "X_array = transformer.fit_transform(X_df)\n",
    "X_array.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combined with some static features and plugged into a random forest it reads;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sklearn.preprocessing import FunctionTransformer\n",
    "from sklearn.compose import make_column_transformer\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.impute import SimpleImputer\n",
    "\n",
    "\n",
    "def fold_time_series(time_point, period, div_period):\n",
    "    return (time_point -\n",
    "            (time_point // (period / div_period)) * period / div_period)\n",
    "\n",
    "\n",
    "def get_bin_means(X_df, num_bins, band):\n",
    "    feature_array = np.empty((len(X_df), num_bins))\n",
    "\n",
    "    for k, (_, x) in enumerate(X_df.iterrows()):\n",
    "        period = x['period']\n",
    "        div_period = x['div_period']\n",
    "        real_period = period / div_period\n",
    "        bins = [i * real_period / num_bins for i in range(num_bins + 1)]\n",
    "\n",
    "        time_points = np.array(x['time_points_' + band])\n",
    "        light_points = np.array(x['light_points_' + band])\n",
    "        time_points_folded = \\\n",
    "            np.array([fold_time_series(time_point, period, div_period)\n",
    "                      for time_point in time_points])\n",
    "        time_points_folded_digitized = \\\n",
    "            np.digitize(time_points_folded, bins) - 1\n",
    "\n",
    "        for i in range(num_bins):\n",
    "            this_light_points = light_points[time_points_folded_digitized == i]\n",
    "            if len(this_light_points) > 0:\n",
    "                feature_array[k, i] = np.mean(this_light_points)\n",
    "            else:\n",
    "                feature_array[k, i] = np.nan  # missing\n",
    "\n",
    "    return feature_array\n",
    "\n",
    "\n",
    "transformer_r = FunctionTransformer(\n",
    "    lambda X_df: get_bin_means(X_df, 5, 'r')\n",
    ")\n",
    "\n",
    "transformer_b = FunctionTransformer(\n",
    "    lambda X_df: get_bin_means(X_df, 5, 'b')\n",
    ")\n",
    "\n",
    "cols = [\n",
    "    'magnitude_b',\n",
    "    'magnitude_r',\n",
    "    'period',\n",
    "    'asym_b',\n",
    "    'asym_r',\n",
    "    'log_p_not_variable',\n",
    "    'sigma_flux_b',\n",
    "    'sigma_flux_r',\n",
    "    'quality',\n",
    "    'div_period',\n",
    "]\n",
    "\n",
    "common = ['period', 'div_period']\n",
    "transformer = make_column_transformer(\n",
    "    (transformer_r, common + ['time_points_r', 'light_points_r']),\n",
    "    (transformer_b, common + ['time_points_b', 'light_points_b']),\n",
    "    ('passthrough', cols)\n",
    ")\n",
    "\n",
    "pipe = make_pipeline(\n",
    "    transformer,\n",
    "    SimpleImputer(strategy='most_frequent'),\n",
    "    RandomForestClassifier(max_depth=5, n_estimators=10)\n",
    ")\n",
    "\n",
    "\n",
    "def get_estimator():\n",
    "    return pipe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing using a scikit-learn pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import problem\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "X_df, y = problem.get_train_data()\n",
    "\n",
    "scores = cross_val_score(get_estimator(), X_df, y, cv=2, scoring='accuracy')\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}