Latent_space_study.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Latent space study code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two separate latent space study for MQXFS1d (3D Latent space) and MDCPT1b (Simpler 2D latent space). So far, latent space of MDCPT1b's model is studied more intensively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. MDCPT1b latent space code (Latent space plot with heat map reconstruction loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More can be found at:\n",
    "\n",
    "https://github.com/Duchstf/quench-detector/blob/15T-exploration/latent-space.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",
    "import h5py\n",
    "from nptdms import TdmsFile #Process ramping file\n",
    "\n",
    "#For building ML models\n",
    "import keras\n",
    "import tensorflow as tf\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\n",
    "\n",
    "def load_sensor(ramp_num, dir_path, sensor, time_range = None):\n",
    "    \"\"\"\n",
    "    Load a sensor's data in the specified time range.\n",
    "    \n",
    "    Example usage: load_sensor(\"Ramp_7\", \"15T_data\", 0, time_range = (-20,0))\n",
    "    \n",
    "    Args:\n",
    "     - ramp_num (str): ramp number\n",
    "     - dir_path (str): path to directory containing all the ramp's data\n",
    "     - sensor (int): sensor index (0 or 1)\n",
    "     - time_range (tuple): time range in which the data should be loaded\n",
    "    \n",
    "    Returns: The specified sensor in specified time range\n",
    "    \"\"\"\n",
    "    \n",
    "    needed_file = None\n",
    "    \n",
    "    #Scan the directory for correct ramp file\n",
    "    for filename in os.listdir(dir_path):\n",
    "        if filename.startswith(ramp_num):\n",
    "            needed_file = dir_path + \"/\" + filename\n",
    "            \n",
    "    #Load the file\n",
    "    f = h5py.File(needed_file,'r')\n",
    "    data = np.array(f['chanvals'], dtype = 'float16')\n",
    "    \n",
    "    #Load time\n",
    "    time_axis = np.load(dir_path + \"/time_axis/\" + ramp_num + \"_time.npy\")\n",
    "    \n",
    "    #Calculate index for selecting data\n",
    "    selection_index = (time_axis > time_range[0]) & (time_axis < time_range[1])\n",
    "    \n",
    "    df_data = pd.DataFrame(data = {\"dummy\": data[sensor, :][selection_index]})\n",
    "    \n",
    "    #Return the sensor's data accoring to the index \n",
    "    return df_data[\"dummy\"]\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "def load_time_label(ramp_num, dir_path, time_range = None, window = 2000, step = 10):\n",
    "    \"\"\"\n",
    "    Load a ramp time axis in a specified range.\n",
    "    \n",
    "    Example usage: load_time_label(\"Ramp_7\", \"15T_data\", time_range = (-20,0))\n",
    "    \n",
    "    Args:\n",
    "     - ramp_num (str): ramp number\n",
    "     - dir_path (str): path to directory containing all the ramp's data\n",
    "     - time_range (tuple): time range in which the data should be loaded\n",
    "    \n",
    "    Returns: The specified time axis in specified time range\n",
    "    \"\"\"\n",
    "    \n",
    "    needed_file = None\n",
    "    \n",
    "    #Scan the directory for correct ramp file\n",
    "    for filename in os.listdir(dir_path):\n",
    "        if filename.startswith(ramp_num):\n",
    "            needed_file = dir_path + \"/\" + filename\n",
    "            \n",
    "    #Load the file\n",
    "    f = h5py.File(needed_file,'r')\n",
    "    data = np.array(f['chanvals'], dtype = 'float16')\n",
    "    \n",
    "    #Load time\n",
    "    time_axis = np.load(dir_path + \"/time_axis/\" + ramp_num + \"_time.npy\")\n",
    "    \n",
    "    #Calculate index for selecting data\n",
    "    selection_index = (time_axis > time_range[0]) & (time_axis < time_range[1])\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    return time_axis[selection_index][window-1::step]\n",
    "\n",
    "\n",
    "#### Mean of abs(signal)\n",
    "def generate_mean_data(ramp_num, dir_path, time_range, window = 2000, step = 10):\n",
    "    #Load the data\n",
    "    ai0 = load_sensor(ramp_num, dir_path, 0, time_range = time_range).abs().rolling(window=window).mean().to_numpy()[window-1::step]\n",
    "    ai1 = load_sensor(ramp_num, dir_path, 1, time_range = time_range).abs().rolling(window=window).mean().to_numpy()[window-1::step]\n",
    "    \n",
    "    #Calculate the product\n",
    "    product = ai0*ai1\n",
    "    \n",
    "    #Stack them together\n",
    "    all_mean = np.vstack((ai0,ai1,product)).transpose()\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    return all_mean\n",
    "\n",
    "#### SD of signal\n",
    "def generate_sd_data(ramp_num, dir_path, time_range, window = 2000, step = 10):\n",
    "    #Load the data\n",
    "    ai0 = load_sensor(ramp_num, dir_path, 0, time_range = time_range).rolling(window=window).std().to_numpy()[window-1::step]\n",
    "    ai1 = load_sensor(ramp_num, dir_path, 1, time_range = time_range).rolling(window=window).std().to_numpy()[window-1::step]\n",
    "    \n",
    "    #Calculate the product\n",
    "    product = ai0*ai1\n",
    "    \n",
    "    #Stack them together\n",
    "    all_sd = np.vstack((ai0,ai1, product)).transpose()\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    return all_sd\n",
    "\n",
    "def generate_data(ramp_num, dir_path, time_range, window = 2000, 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(ramp_num, dir_path, time_range, window = window, step = step)\n",
    "    moving_sd = generate_sd_data(ramp_num, dir_path, time_range, window = window, step = step)\n",
    "    time_label = load_time_label(ramp_num, dir_path, time_range, window = window, step = step)\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",
    "\n",
    "############DYNAMIC LEARNING IMPLEMENTATION############\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(4, activation = 'elu', kernel_initializer='glorot_uniform',\n",
    "                    kernel_regularizer=regularizers.l2(0.0),\n",
    "                    input_dim=6))\n",
    "\n",
    "    model.add(Dense(2,activation='elu',\n",
    "                    kernel_initializer='glorot_uniform', name = 'latent_space'))\n",
    "\n",
    "    model.add(Dense(4,activation='elu',\n",
    "                    kernel_initializer='glorot_uniform'))\n",
    "\n",
    "    model.add(Dense(6,\n",
    "                    kernel_initializer='glorot_uniform'))\n",
    "\n",
    "    model.compile(loss='mse',optimizer='adam')\n",
    "\n",
    "\n",
    "    return model\n",
    "\n",
    "def get_latent_output(model_weights, X):\n",
    "    \n",
    "    model = create_model()  # create the original model\n",
    "    model.set_weights(model_weights)\n",
    "\n",
    "    layer_name = 'latent_space'\n",
    "    intermediate_layer_model = keras.Model(inputs=model.input, outputs=model.get_layer(layer_name).output)\n",
    "    intermediate_layer_model.compile(loss='mse',optimizer='adam')\n",
    "    \n",
    "    intermediate_output = intermediate_layer_model.predict(X)\n",
    "    \n",
    "    return intermediate_output\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",
    "    #Get latent space output\n",
    "    latent_space = get_latent_output(current_weights, X_test)\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 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",
    "        #Get latent space of normal data points\n",
    "        #normal_latent = np.delete(latent_space, peaks, 0)\n",
    "        #normal_log_test_loss = np.delete(log_test_loss, peaks, 0)\n",
    "        \n",
    "        return list(time_test[peaks]), [list(x) for x in list(latent_space[peaks])], latent_space, log_test_loss, time_test\n",
    "    \n",
    "    else:\n",
    "        return None, None, None, None, None\n",
    "    \n",
    "\n",
    "def simulation(ramp_num, dir_path, section_size, time_range, reached_max):\n",
    "    \n",
    "    #Load data and create the model\n",
    "    X, time_axis = generate_data(ramp_num, 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",
    "    #Keep track of latent coordinate of trigger time\n",
    "    latent_coor = []\n",
    "    \n",
    "    #Keep track of normal latent space\n",
    "    latent_normal_coord = None\n",
    "    \n",
    "    #Test loss associated with latent space\n",
    "    all_test_loss = None\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(changing_range/section_size - (changing_range/section_size)%1)\n",
    "    \n",
    "    print(\"Number of triggering sections\", trigger_sections)\n",
    "    \n",
    "    trigger_thres = create_threshold_list(trigger_sections, start_thres = 3, end_thres = 1.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",
    "            #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_latent, normal_latent, log_test_loss, time_test = 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(\"Latent coordinates are: \", trigger_latent)\n",
    "                latent_coor.extend(trigger_latent)\n",
    "                \n",
    "                if latent_normal_coord is None:\n",
    "                    latent_normal_coord = normal_latent\n",
    "                else:\n",
    "                    latent_normal_coord = np.append(latent_normal_coord, normal_latent, 0)\n",
    "                    \n",
    "                if all_test_loss is None:\n",
    "                    all_test_loss = log_test_loss\n",
    "                else:\n",
    "                    all_test_loss = np.append(all_test_loss, log_test_loss, 0)\n",
    "        \n",
    "        print(\"<--------\")\n",
    "        \n",
    "    \n",
    "    #Create the animation here.\n",
    "    \n",
    "    #Plot latent space of normal data points\n",
    "    # Sort the points by density, so that the densest points are plotted last\n",
    "    idx = all_test_loss.argsort()\n",
    "    sorted_latent_normal_coord, sorted_all_test_loss  = latent_normal_coord[idx, :], all_test_loss[idx]\n",
    "    \n",
    "    fig, ax = plt.subplots()\n",
    "    norm_latent = ax.scatter(sorted_latent_normal_coord[:, 0], sorted_latent_normal_coord[:,1], c = sorted_all_test_loss, edgecolor='', cmap = 'inferno')\n",
    "    ax.scatter([x[0] for x in latent_coor], [x[1] for x in latent_coor], color = 'red', label = \"Anomaly events\")\n",
    "    plt.title(\"Latent space of normal points in {}\".format(ramp_num))\n",
    "    plt.legend(loc = 'best')\n",
    "    \n",
    "    #Set color bar\n",
    "    cbar = plt.colorbar(norm_latent)\n",
    "    cbar.set_label(r'$\\log_{10}$(loss)')\n",
    "    plt.show()\n",
    "    \n",
    "    return latent_normal_coord, all_test_loss, latent_coor, all_test_loss, time_test\n",
    "\n",
    "    print(\"Finished simulation, all the triggered times are: \", trigger_log)\n",
    "    print()\n",
    "    print(\"All latent coordinate: \", latent_coor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example usage\n",
    "#latent_normal_coord, all_test_loss, latent_coor, all_test_loss, time_test  = simulation(\"Ramp_57\", \"15T_data\", 10,  time_range = (-50,0), reached_max = -1.55)\n",
    "#Can use the outputs to obtain plots\n",
    "#More can be found in the results section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. MQXFS1d latent space code (Latent space plot in 3D)\n",
    "\n",
    "\n",
    "More can be found at:\n",
    "\n",
    "https://github.com/Duchstf/quench-detector/blob/MQFSX1D_latent/latent-space.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",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import seaborn as sns\n",
    "from nptdms import TdmsFile #Process ramping file\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\n",
    "\n",
    "######----PROCESSING DATA------########\n",
    "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",
    "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, step = 10):\n",
    "    \n",
    "    moving_mean = generate_mean_data(dir_path, time_range, window = window, step = step)\n",
    "    moving_sd = generate_sd_data(dir_path, time_range, window = window, step = step)\n",
    "    time_label = load_time_label(dir_path, time_range, window = window, step = step)\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', name = 'latent_space'))\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 get_latent_output(model_weights, X):\n",
    "    \n",
    "    model = create_model()  # create the original model\n",
    "    model.set_weights(model_weights)\n",
    "\n",
    "    layer_name = 'latent_space'\n",
    "    intermediate_layer_model = keras.Model(inputs=model.input, outputs=model.get_layer(layer_name).output)\n",
    "    intermediate_layer_model.compile(loss='mse',optimizer='adam')\n",
    "    \n",
    "    intermediate_output = intermediate_layer_model.predict(X)\n",
    "    \n",
    "    return intermediate_output\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",
    "    #Get latent space output\n",
    "    latent_space = get_latent_output(current_weights, X_test)\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 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",
    "        #Get latent space of normal data points\n",
    "        normal_latent = np.delete(latent_space, peaks, 0)\n",
    "        \n",
    "        return list(time_test[peaks]), latent_space[peaks], normal_latent\n",
    "    \n",
    "    else:\n",
    "        return None, None, None\n",
    "\n",
    "def plot_3d_latent(latent_normal_coord, latent_coor):\n",
    "    fig = plt.figure(figsize = (10,10))\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    \n",
    "    ax.scatter(latent_normal_coord[:,0], latent_normal_coord[:,1], latent_normal_coord[:,2],\n",
    "               color = 'deepskyblue', label = \"normal latent\")\n",
    "    ax.scatter(latent_coor[:,0], latent_coor[:,1], latent_coor[:,2],\n",
    "               color = 'red', label = \"anomaly latent\", alpha = 1, marker = 's', s = 50)\n",
    "    \n",
    "    plt.title(\"Latent space of normal vs abnormal points\")\n",
    "    plt.legend(loc = \"best\")\n",
    "    plt.show()\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",
    "    #Keep track of latent coordinate of trigger time\n",
    "    latent_coor = None\n",
    "    \n",
    "    #Keep track of normal latent space\n",
    "    latent_normal_coord = None\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(changing_range/section_size - (changing_range/section_size)%1)\n",
    "    \n",
    "    print(\"Number of triggering sections\", trigger_sections)\n",
    "    \n",
    "    trigger_thres = create_threshold_list(trigger_sections, start_thres = 5.0, 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",
    "            #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_latent, normal_latent = 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(\"Latent coordinates are: \", trigger_latent)\n",
    "                \n",
    "                if latent_coor is None:\n",
    "                    latent_coor = trigger_latent\n",
    "                else:\n",
    "                    latent_coor = np.append(latent_coor, trigger_latent, 0)\n",
    "                \n",
    "                print(\"latent_coor's shape: \", latent_coor.shape)\n",
    "                \n",
    "                if latent_normal_coord is None:\n",
    "                    latent_normal_coord = normal_latent\n",
    "                else:\n",
    "                    latent_normal_coord = np.append(latent_normal_coord, normal_latent, 0)\n",
    "                \n",
    "                print(\"Latent_normal_coord's shape: \", latent_normal_coord.shape)\n",
    "                \n",
    "                #Plot the latent space at this point\n",
    "                plot_3d_latent(latent_normal_coord, latent_coor)\n",
    "        \n",
    "        print(\"<--------\")\n",
    "    \n",
    "    #Plot the latent space when done\n",
    "    plot_3d_latent(latent_normal_coord, latent_coor)\n",
    "    \n",
    "    print(\"Finished simulation, all the triggered times are: \", trigger_log)\n",
    "    print()\n",
    "    print(\"All latent coordinate: \", latent_coor)"
   ]
  }
 ],
 "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
}