start_z_BDT.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy\n",
    "from iminuit import cost\n",
    "from iminuit import Minuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use([\"science\", \"notebook\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"font.size\"] = 14\n",
    "plt.rcParams[\"axes.formatter.limits\"] = -5, 4\n",
    "plt.rcParams[\"figure.figsize\"] = 6, 4\n",
    "colors = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"features_CNN_1d_99987.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.load(\"images_1d_99987.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = scaler.transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(X[:100, :])\n",
    "plotting.watermark()\n",
    "plt.xlabel(r\"$2\\times \\mathrm{station} + \\mathrm{plane}$\")\n",
    "plt.ylabel(\"Event\")\n",
    "plt.savefig(\"plots/1d_data.png\")\n",
    "plt.savefig(\"plots/1d_data.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import AdaBoostRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skopt.space import Real, Integer\n",
    "from skopt import BayesSearchCV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = df[\"start_z\"].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_scaler = StandardScaler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = y.reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_scaler.fit(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = y_scaler.transform(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = (y + 235) / 155"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "def y_to_z(y):\n",
    "    return (y * 155) - 235"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_to_z(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "metadata": {},
   "outputs": [],
   "source": [
    "z = y_scaler.inverse_transform(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "y.ravel().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y.ravel(), random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [],
   "source": [
    "estimator = AdaBoostRegressor(random_state=0, learning_rate=0.1, n_estimators=466)\n",
    "\n",
    "estimator.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {},
   "outputs": [],
   "source": [
    "bdt_opt = BayesSearchCV(\n",
    "    AdaBoostRegressor(),\n",
    "    {\n",
    "        \"learning_rate\": Real(0.005, 0.9, prior=\"log-uniform\"),\n",
    "        \"n_estimators\": Integer(1, 1000),\n",
    "    },\n",
    "    n_iter=100,\n",
    "    cv=5,\n",
    ")\n",
    "\n",
    "bdt_opt.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31",
   "metadata": {},
   "outputs": [],
   "source": [
    "bdt_opt.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = estimator.predict(X_test)\n",
    "\n",
    "estimator.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_pred = y_scaler.inverse_transform(y_pred.reshape(-1, 1)).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_test = y_scaler.inverse_transform(y_test.reshape(-1, 1)).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_pred = y_to_z(y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36",
   "metadata": {},
   "outputs": [],
   "source": [
    "h = hist.Hist.new.Regular(20, -10, +10, name=r\"𝛥z [cm]\").Double()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37",
   "metadata": {},
   "outputs": [],
   "source": [
    "h.fill(np.squeeze(z_pred) - np.squeeze(z_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38",
   "metadata": {},
   "outputs": [],
   "source": [
    "entries, edges = h.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39",
   "metadata": {},
   "outputs": [],
   "source": [
    "def residual_model(x, mu, sigma):\n",
    "    return scipy.stats.norm.cdf(x, mu, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Minuit(cost.BinnedNLL(entries, edges, residual_model), 0, 25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = m.migrad()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42",
   "metadata": {},
   "outputs": [],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43",
   "metadata": {},
   "outputs": [],
   "source": [
    "h.plot()\n",
    "plt.xlabel(r\"$\\Delta z\\;[\\mathrm{cm}]$\")\n",
    "plot_range = edges[0], edges[-1]\n",
    "x = np.linspace(*plot_range, 100)\n",
    "best_fit = scipy.stats.norm(res.params[0].value, res.params[1].value)\n",
    "# best_fit = scipy.stats.norm(0.044, 2.83) # TODO take from fit\n",
    "n_bins = len(entries)\n",
    "binsize = (plot_range[1] - plot_range[0]) / n_bins\n",
    "scale = h.sum() / (best_fit.cdf(plot_range[1]) - best_fit.cdf(plot_range[0])) * binsize\n",
    "plt.plot(x, scale * best_fit.pdf(x))\n",
    "ax = plt.gca()\n",
    "# plt.text(0.6, 0.9, r\"$\\mu = 0.044 $\\;cm\", transform=ax.transAxes, usetex=True)\n",
    "plt.text(\n",
    "    0.6,\n",
    "    0.9,\n",
    "    rf\"$\\mu = {res.params[0].value:.2f} \\pm {res.params[0].error:.2f}$\\;cm\",\n",
    "    transform=ax.transAxes,\n",
    "    usetex=True,\n",
    ")\n",
    "plt.text(\n",
    "    0.25,\n",
    "    0.1,\n",
    "    f\"Training dataset: {len(y_train)} events\\n\"\n",
    "    f\"Test dataset: {len(y_test)} events\\nAdaBoost BDT\",\n",
    "    transform=ax.transAxes,\n",
    "    usetex=True,\n",
    ")\n",
    "# plt.text(0.6, 0.81, r\"$\\sigma = 2.83 $\\;cm\", transform=ax.transAxes, usetex=True)\n",
    "plt.text(\n",
    "    0.6,\n",
    "    0.81,\n",
    "    rf\"$\\sigma = {res.params[1].value:.2f} \\pm {res.params[1].error:.2f}$\\;cm\",\n",
    "    transform=ax.transAxes,\n",
    "    usetex=True,\n",
    ")\n",
    "plotting.watermark()\n",
    "plt.savefig(f\"plots/h_dz_AdaBoost_n{len(y)}.pdf\")\n",
    "plt.savefig(f\"plots/h_dz_AdaBoost_n{len(y)}.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(y_to_z(y_train))\n",
    "plt.hist(y_to_z(y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45",
   "metadata": {},
   "outputs": [],
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}