MQXFS1d_dynamic_learning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MQXFS1d dynamic learning code documentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This page documents the code for the dynamic learning algorithm and dynamic threshold setting as presented in paper. More can be found in the trigger branch.\n",
    "\n",
    "The results can be viewed in this page:\n",
    "\n",
    "https://github.com/Duchstf/quench-detector/blob/trigger/dynamic-thres-1.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Data processing\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from nptdms import TdmsFile #Process ramping file\n",
    "import math\n",
    "\n",
    "#For building ML models\n",
    "import keras\n",
    "import keras.models as models\n",
    "from keras.layers.core import Dense\n",
    "from keras import regularizers\n",
    "from sklearn.preprocessing import MinMaxScaler, StandardScaler\n",
    "from scipy.signal import find_peaks\n",
    "\n",
    "np.random.seed(1337) # for reproducibility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_sensor(dir_path, sensor, time_range = None):\n",
    "    \n",
    "    data = pd.DataFrame(data = {sensor: np.load(dir_path + sensor + \".npy\"),\n",
    "                                \"time\": np.load(dir_path + \"time.npy\")})\n",
    "    \n",
    "    start = min(data[\"time\"])\n",
    "    end = max(data[\"time\"])\n",
    "    \n",
    "    if time_range:\n",
    "        start = time_range[0]\n",
    "        end = time_range[1]\n",
    "    \n",
    "    \n",
    "    data = data[(data[\"time\"] > start) & (data[\"time\"] < end)]\n",
    "   \n",
    "    %reset -f in\n",
    "    \n",
    "    return data[sensor]\n",
    "\n",
    "def load_time_label(dir_path, time_range, window = 2000, step = 10):\n",
    "    \n",
    "    time_label =  np.load(dir_path + \"time.npy\")\n",
    "    \n",
    "    start = min(time_label)\n",
    "    end = max(time_label)\n",
    "    \n",
    "    if time_range:\n",
    "        start = time_range[0]\n",
    "        end = time_range[1]\n",
    "    \n",
    "    \n",
    "    time_label = time_label[(time_label > start) & (time_label < end)][window-1::step]\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    return time_label\n",
    "\n",
    "#### Mean of abs(signal)\n",
    "def generate_mean_data(dir_path, time_range, window = 2000, step = 10):\n",
    "    #Load the data\n",
    "    ai0 = load_sensor(dir_path, \"ai0\", time_range = time_range).abs().rolling(window=window).mean().to_numpy()[window-1::step]\n",
    "    ai1 = load_sensor(dir_path, \"ai1\", time_range = time_range).abs().rolling(window=window).mean().to_numpy()[window-1::step]\n",
    "    ai2 = load_sensor(dir_path, \"ai2\", time_range = time_range).abs().rolling(window=window).mean().to_numpy()[window-1::step]\n",
    "    ai3 = load_sensor(dir_path, \"ai3\", time_range = time_range).abs().rolling(window=window).mean().to_numpy()[window-1::step]\n",
    "    ai4 = load_sensor(dir_path, \"ai4\", time_range = time_range).abs().rolling(window=window).mean().to_numpy()[window-1::step]\n",
    "    \n",
    "    #Calculate the product\n",
    "    product = ai0*ai1*ai2*ai3*ai4\n",
    "    \n",
    "    #Stack them together\n",
    "    all_mean = np.vstack((ai0,ai1,ai2,ai3,ai4, product)).transpose()\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    return all_mean\n",
    "\n",
    "#### SD of signal\n",
    "def generate_sd_data(dir_path, time_range, window = 2000, step = 10):\n",
    "    #Load the data\n",
    "    ai0 = load_sensor(dir_path, \"ai0\", time_range = time_range).rolling(window=window).std().to_numpy()[window-1::step]\n",
    "    ai1 = load_sensor(dir_path, \"ai1\", time_range = time_range).rolling(window=window).std().to_numpy()[window-1::step]\n",
    "    ai2 = load_sensor(dir_path, \"ai2\", time_range = time_range).rolling(window=window).std().to_numpy()[window-1::step]\n",
    "    ai3 = load_sensor(dir_path, \"ai3\", time_range = time_range).rolling(window=window).std().to_numpy()[window-1::step]\n",
    "    ai4 = load_sensor(dir_path, \"ai4\", time_range = time_range).rolling(window=window).std().to_numpy()[window-1::step]\n",
    "    \n",
    "    #Calculate the product\n",
    "    product = ai0*ai1*ai2*ai3*ai4\n",
    "    \n",
    "    #Stack them together\n",
    "    all_sd = np.vstack((ai0,ai1,ai2,ai3,ai4, product)).transpose()\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    return all_sd\n",
    "\n",
    "#### Thresholded zero crossings\n",
    "def zero_crossings(array):\n",
    "    \n",
    "    #Set values outside of range (0.001, 0.1) = 0\n",
    "    array[abs(array) > 0.1] = 0\n",
    "    array[abs(array) < 0.005] = 0\n",
    "    \n",
    "    #Calculate number of zero-crossing points, normalized by the window size\n",
    "    zero_crossings = ((array[:-1] * array[1:]) < 0).sum()/array.size\n",
    "    \n",
    "    return zero_crossings\n",
    "\n",
    "\n",
    "def generate_crossings_data(dir_path, time_range = None, window = 5000, step = 10):\n",
    "    \n",
    "    #Load the data\n",
    "    ai0 = load_sensor(dir_path, \"ai0\", time_range = time_range).rolling(window=window).apply(zero_crossings, raw = True).to_numpy()[window-1::step]\n",
    "    ai1 = load_sensor(dir_path, \"ai1\", time_range = time_range).rolling(window=window).apply(zero_crossings, raw = True).to_numpy()[window-1::step]\n",
    "    ai2 = load_sensor(dir_path, \"ai2\", time_range = time_range).rolling(window=window).apply(zero_crossings, raw = True).to_numpy()[window-1::step]\n",
    "    ai3 = load_sensor(dir_path, \"ai3\", time_range = time_range).rolling(window=window).apply(zero_crossings, raw = True).to_numpy()[window-1::step]\n",
    "    ai4 = load_sensor(dir_path, \"ai4\", time_range = time_range).rolling(window=window).apply(zero_crossings, raw = True).to_numpy()[window-1::step]\n",
    "    \n",
    "    product = ai0*ai1*ai2*ai3*ai4\n",
    "    \n",
    "    #Stack them together\n",
    "    all_crossings = np.vstack((ai0,ai1,ai2,ai3,ai4, product)).transpose()\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    return all_crossings\n",
    "\n",
    "\n",
    "def generate_data(dir_path, time_range, window = 2000, thres_win = 25000, step = 10):\n",
    "    \n",
    "    #Selection index due to using different window\n",
    "    selection_index = int((thres_win - window)/step)\n",
    "    \n",
    "    moving_mean = generate_mean_data(dir_path, time_range, window = window, step = step)[selection_index:]\n",
    "    moving_sd = generate_sd_data(dir_path, time_range, window = window, step = step)[selection_index:]\n",
    "    #moving_thres_crossings = generate_crossings_data(dir_path, time_range, window = thres_win, step = step)\n",
    "    time_label = load_time_label(dir_path, time_range, window = window, step = step)[selection_index:]\n",
    "    \n",
    "    all_data = np.concatenate((moving_mean, moving_sd), axis = 1)\n",
    "    \n",
    "    print(\"All data's shape: \", all_data.shape)\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    return all_data, time_label\n",
    "\n",
    "#####-------INCREMENTAL LEARNING IMPLEMENTATION-------######\n",
    "\n",
    "def create_threshold_list(trigger_sections, start_thres = 4.0, end_thres = 2.6):\n",
    "    \n",
    "    increment = (start_thres - end_thres)/(trigger_sections-3)\n",
    "    \n",
    "    trigger_thres = [0,0] #No trigger in the first two sections\n",
    "    \n",
    "    trigger_thres.extend([start_thres - i*increment for i in range(trigger_sections - 2)])\n",
    "    \n",
    "    return trigger_thres\n",
    "\n",
    "def create_model():\n",
    "    \n",
    "    #================BUILD THE MODEL====================\n",
    "    # Simple model\n",
    "    # Simple model\n",
    "    model = models.Sequential()\n",
    "\n",
    "    model.add(Dense(6, activation = 'elu', kernel_initializer='glorot_uniform',\n",
    "                    kernel_regularizer=regularizers.l2(0.0),\n",
    "                    input_dim=12))\n",
    "\n",
    "    model.add(Dense(3,activation='elu',\n",
    "                    kernel_initializer='glorot_uniform'))\n",
    "\n",
    "    model.add(Dense(6,activation='elu',\n",
    "                    kernel_initializer='glorot_uniform'))\n",
    "\n",
    "    model.add(Dense(12,\n",
    "                    kernel_initializer='glorot_uniform'))\n",
    "\n",
    "    model.compile(loss='mse',optimizer='adam')\n",
    "\n",
    "\n",
    "    return model\n",
    "\n",
    "def train(current_weights, X, time_axis, time_range):\n",
    "    \"\"\"Take a model, train it in the appropritate time range, return the new weights\n",
    "    and its reconstruction loss distribution\"\"\"\n",
    "    \n",
    "    print(\"Training for section: \", time_range)\n",
    "    #Select the appropriate time\n",
    "    start_index = None\n",
    "    \n",
    "    if time_range[0] is not None:\n",
    "        start_index = np.argmax(time_axis > time_range[0])\n",
    "    else:\n",
    "        start_index = 0\n",
    "    \n",
    "    if time_range[1] == 0:\n",
    "        end_index = None\n",
    "    else:\n",
    "        end_index = np.argmax(time_axis > time_range[1])\n",
    "    \n",
    "    X_train = np.copy(X[start_index:end_index])\n",
    "    time_train = time_axis[start_index:end_index]\n",
    "    \n",
    "    #Create model\n",
    "    model = create_model()\n",
    "    \n",
    "    if current_weights:\n",
    "        model.set_weights(current_weights)\n",
    "    \n",
    "    #Scaler initiation\n",
    "    scaler = MinMaxScaler()\n",
    "    scaler.fit(X_train)\n",
    "    X_train = scaler.transform(X_train)\n",
    "    \n",
    "    #Train the model and extract new weights\n",
    "    model.fit(X_train, X_train,\n",
    "              batch_size=2000,\n",
    "              epochs=30, verbose = 0)\n",
    "    \n",
    "    new_weights = model.get_weights()\n",
    "    \n",
    "    #Evaluate the losses and update the distribution\n",
    "    X_pred = model.predict(X_train)\n",
    "    loss = np.mean(np.abs(X_pred-X_train), axis = 1)\n",
    "    \n",
    "    #Plot the reconstruction loss\n",
    "    plt.figure(figsize=(20,2))\n",
    "    plt.plot(time_train, loss)\n",
    "    plt.title(\"TRAINING reconstruction loss on {}\".format(time_range))\n",
    "    plt.xlabel(\"Time [s]\")\n",
    "    plt.show()\n",
    "    \n",
    "    #Take the mean of the log loss and update the mean threshold\n",
    "    current_median = np.median(np.log10(loss))\n",
    "    \n",
    "    print(\"Finished training, median of the log(loss) is: \", current_median)\n",
    "    \n",
    "    return new_weights, current_median, scaler\n",
    "\n",
    "def trigger(trigger_threshold, current_weights, X, time_axis, time_range, prev_median, scaler):\n",
    "    \n",
    "    print(\"Triggering for section: \", time_range)\n",
    "    print(\"Using this median for triggering:\", prev_median)\n",
    "    print(\"Triggering threshold is: \", trigger_threshold)\n",
    "    \n",
    "    #Select the appropriate time\n",
    "    start_index = None\n",
    "    \n",
    "    if time_range[0] is not None:\n",
    "        start_index = np.argmax(time_axis > time_range[0])\n",
    "    else:\n",
    "        start_index = 0\n",
    "    \n",
    "    if time_range[1] == 0:\n",
    "        end_index = None\n",
    "    else:\n",
    "        end_index = np.argmax(time_axis > time_range[1])\n",
    "    \n",
    "    time_test = time_axis[start_index:end_index]\n",
    "    X_test = np.copy(X[start_index:end_index])\n",
    "    \n",
    "    X_test = scaler.transform(X_test) #Scale the current data\n",
    "    \n",
    "    #Create model\n",
    "    model = create_model()\n",
    "    model.set_weights(current_weights)\n",
    "        \n",
    "    #Find the test loss (doing this in parallel for the sake of simulation), but in a practical system\n",
    "    #Ideally we need to do inference every single data point come in\n",
    "    X_pred_test = model.predict(X_test)\n",
    "    test_loss = np.mean(np.abs(X_pred_test-X_test), axis = 1)\n",
    "    \n",
    "    \n",
    "    #Plot the reconstruction loss just to check\n",
    "    plt.figure(figsize=(20,2))\n",
    "    plt.plot(time_test, test_loss)\n",
    "    plt.title(\"Reconstruction loss on {}\".format(time_range))\n",
    "    plt.xlabel(\"Time [s]\")\n",
    "    plt.show()\n",
    "    \n",
    "    #Quantify how much each data point is far away from the distribution\n",
    "    log_test_loss = np.log10(test_loss) - prev_median\n",
    "    \n",
    "    #Find the trigger time\n",
    "    trigger_index = np.argmax(log_test_loss > trigger_threshold)\n",
    "    \n",
    "    if  trigger_index != 0:\n",
    "        #Plot the triggered time\n",
    "        plt.figure(figsize=(20,2))\n",
    "        plt.plot(time_test, log_test_loss, color = \"firebrick\")\n",
    "        plt.title(\"Log reconstruction loss on {}\".format(time_range))\n",
    "        plt.xlabel(\"Time [s]\")\n",
    "        \n",
    "        #Trigger time\n",
    "        plt.vlines(time_test[trigger_index], log_test_loss.min(), log_test_loss.max(), color = \"blue\", linestyles = \"dashed\", label = \"Trigger time\")\n",
    "        plt.legend(loc = \"upper left\")\n",
    "        \n",
    "        \n",
    "        ### Find trigger times\n",
    "        peaks, _ = find_peaks(log_test_loss, height=trigger_threshold, distance = 5000)\n",
    "        print(\"Other trigger time: \", time_test[peaks])\n",
    "        print(\"Their log loss values are: \", log_test_loss[peaks])\n",
    "        \n",
    "        for i in range(len(peaks)):\n",
    "            plt.vlines(time_test[peaks[i]], log_test_loss.min(), log_test_loss.max(), color = \"blue\", linestyles = \"dashed\")\n",
    "            \n",
    "        plt.show()\n",
    "        \n",
    "        return list(time_test[peaks])\n",
    "    \n",
    "    else:\n",
    "        return None\n",
    "    \n",
    "\n",
    "def simulation(dir_path, section_size, time_range, reached_max):\n",
    "    \n",
    "    #Load data and create the model\n",
    "    X, time_axis = generate_data(dir_path, time_range = time_range, window = 2000, step = 10)\n",
    "    \n",
    "    #For updating model\n",
    "    current_weights = None\n",
    "    new_weights = None\n",
    "    \n",
    "    #For keeping track of distribution history\n",
    "    prev_median = None\n",
    "    updated_median = None\n",
    "    \n",
    "    #Keeping track of scaler after training\n",
    "    past_scaler = None\n",
    "    \n",
    "    #Calculate number of sections\n",
    "    num_section = int(abs(time_range[0])/section_size)\n",
    "    print(\"Number of sections:\", num_section)\n",
    "    \n",
    "    #Keeping log of all the trigger time\n",
    "    trigger_log = []\n",
    "    \n",
    "    #Setting dynamic trigger threshold\n",
    "    start_time = time_range[0]\n",
    "    \n",
    "    changing_range = abs(start_time - reached_max)\n",
    "    \n",
    "    trigger_sections = int(math.ceil(changing_range/section_size))\n",
    "    \n",
    "    print(\"Number of triggering sections\", trigger_sections)\n",
    "    \n",
    "    trigger_thres = create_threshold_list(trigger_sections, start_thres = 5, end_thres = 2.6)\n",
    "    \n",
    "    print(\"List of trigger thresholds: \", trigger_thres)\n",
    "    \n",
    "    #Sequentially go through the data\n",
    "    for i in range(num_section):\n",
    "        \n",
    "        print(\"Entering sections: \", (-section_size*(num_section - i), -section_size*(num_section - i - 1)))\n",
    "        print(\"-------->\")\n",
    "        \n",
    "        if i == 0:\n",
    "            continue\n",
    "        elif i == 1:\n",
    "            #Train on the first segment\n",
    "            time_range = (None, -section_size*(num_section - i))\n",
    "            new_weights, updated_median, new_scaler = train(current_weights, X, time_axis, time_range)\n",
    "        else:\n",
    "            \n",
    "            #Set previous distrubution to the newly trained disttribution\n",
    "            prev_median = updated_median\n",
    "            current_scaler = new_scaler\n",
    "            current_weights = new_weights\n",
    "            \n",
    "            #Selecting appropriate time range for training and triggering\n",
    "            train_time_range = (-section_size*(num_section - i + 1), -section_size*(num_section - i))\n",
    "            trigger_time_range = (-section_size*(num_section - i), -section_size*(num_section - i - 1))\n",
    "            \n",
    "            #Train and predict simutaneously\n",
    "            new_weights, updated_median, new_scaler = train(current_weights, X, time_axis, train_time_range)\n",
    "            \n",
    "            \n",
    "            #Select the trigger threshold\n",
    "            \n",
    "            if i >= len(trigger_thres):\n",
    "                trigger_threshold = trigger_thres[-1]\n",
    "            else:\n",
    "                trigger_threshold = trigger_thres[i]\n",
    "                \n",
    "            \n",
    "            trigger_time = trigger(trigger_threshold, current_weights, X, time_axis, trigger_time_range, prev_median, current_scaler)\n",
    "            \n",
    "            if trigger_time is not None:\n",
    "                print(\"Triggered at: \", trigger_time)\n",
    "                trigger_log.extend(trigger_time)\n",
    "        \n",
    "        print(\"<--------\")\n",
    "        \n",
    "    print(\"Finished simulation, all the triggered times are: \", trigger_log)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}