{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plasma detection from $U_{loop}$\n",
    "\n",
    "This is a quick plasma start and end detection from the peaks in the derivative of the $U_{loop}$ signal.\n",
    "This procedure is run as soon as possible in order to provide the information on plasma existance to other routines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal  # for peak detection\n",
    "import pandas as pd\n",
    "import requests\n",
    "import subprocess #workarround for db operation\n",
    "import os\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shot_no = 38422"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_db_current_shot_alchemy(field_name, value):\n",
    "    try:\n",
    "        # local UNIX socket connection\n",
    "        engine = sqlalchemy.create_engine('postgresql://golem@/golem_database?host=/var/run/postgresql')\n",
    "        #engine = sqlalchemy.create_engine('postgresql://golem:rabijosille@192.168.2.116/golem_database')\n",
    "    except:\n",
    "        return\n",
    "    engine.execute(f\"\"\"UPDATE shots SET \"{field_name}\"={value} WHERE shot_no IN(SELECT max(shot_no) FROM shots)\"\"\")  \n",
    "\n",
    "def update_db_current_shot(field_name, value):\n",
    "    os.system('export PGPASSWORD=\"rabijosille\";psql -c \"UPDATE shots SET '+field_name+'='+str(value)+'WHERE shot_no IN(SELECT max(shot_no) FROM shots)\" -q -U golem golem_database')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs('Results', exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_scalar(phys_quant, value, format_str='%.3f'):\n",
    "    with open(\"Results/\"+phys_quant, 'w') as f:\n",
    "        f.write(format_str % value)\n",
    "    update_db_current_shot(phys_quant,value) \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_CD = requests.get(\n",
    "    f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Production/Parameters/Tcd')\n",
    "t_CD = float(t_CD.content) * 1e-3  # from us to ms\n",
    "t_CD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U_loop = pd.read_csv(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/U_Loop.csv',\n",
    "                    names=['Time', 'U_loop'], index_col='Time', squeeze=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The processing is done only beyond the CD system activation time (to prevent thyristor noise from e.g. the $B_t$ system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U_loop.index *= 1e3  # s to ms "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U_loop = U_loop.loc[t_CD:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = np.diff(U_loop.index[[0,-1]]).item() / U_loop.index.size\n",
    "dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_points(time, must_be_odd=False):\n",
    "    n = int(np.rint(time / dt))\n",
    "    if n % 2 == 0 and must_be_odd:\n",
    "        n += 1\n",
    "    return n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate $\\frac{d U_{loop}}{dt}$ by a second order Savitzky-Golay 1. derivative filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Dt = 0.1   # window length in ms\n",
    "dU = pd.Series(signal.savgol_filter(U_loop.loc[t_CD:], to_points(Dt, must_be_odd=True),\n",
    "                          2, 1, delta=dt), index=U_loop.index, name='dU_dt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def argrelmax_time(sig, Dt, res_index, threshold_value=None, \n",
    "                   comparator=np.greater):\n",
    "    idxs = signal.argrelextrema(sig.values, comparator, order=to_points(Dt))\n",
    "    sig_sel = sig.iloc[idxs]\n",
    "    if threshold_value is not None:\n",
    "        sig_sel = sig_sel[comparator(sig_sel, threshold_value)]\n",
    "    return sig_sel.index[res_index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plasma starts when the $U_{loop}$ decreases significantly for the first  (index 0) time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    t_plasma_start = argrelmax_time(dU, Dt/2, 0, threshold_value=-3,\n",
    "                                    comparator=np.less)\n",
    "except IndexError:\n",
    "    t_plasma_start = -1.0    # not found, nan convention\n",
    "t_plasma_start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plasma ends when the $U_{loop}$ increases significantly for the last  (index -1) time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    t_plasma_end = argrelmax_time(dU.loc[t_plasma_start:], Dt/2, -1,\n",
    "                                  threshold_value=10,\n",
    "                                  comparator=np.greater)\n",
    "except IndexError:\n",
    "    t_plasma_end = -1.0    # not found, nan convention\n",
    "t_plasma_end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "is_plasma = int(t_plasma_start > 0 and t_plasma_end > 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if is_plasma:\n",
    "    t_plasma_duration = t_plasma_end - t_plasma_start\n",
    "else:\n",
    "    t_plasma_duration = -1.0  # convention instead of nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'plasma: {bool(is_plasma)}, from {t_plasma_start:.1f} to {t_plasma_end:.1f} ms (lifetime {t_plasma_duration:.1f} ms)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_scalar(\"is_plasma\", is_plasma)\n",
    "save_scalar(\"t_plasma_start\", t_plasma_start)\n",
    "save_scalar(\"t_plasma_end\", t_plasma_end)\n",
    "save_scalar(\"t_plasma_duration\", t_plasma_duration)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}