data_processing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data processing code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This page documents all the data processing code used in this project to convert `*.tdms` or `*.mat` acoustic signal files to Python readable format. There is also a miscelaneous section where codes on processing text file of voltage taps and other stuffs are presented, but might not be as well documented.\n",
    "\n",
    "Contents:\n",
    "\n",
    "    1. tdms processing\n",
    "    2. mat processing\n",
    "    3. Miscellaneous codes\n",
    "    \n",
    "Note that these are just documented codes, you can checkout other branches to see the results of the runs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. tdms processing"
   ]
  },
  {
   "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 nptdms import TdmsFile #Process ramping file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "These are codes that read tdms files and break them into .npy formats. There are also functions to \n",
    "support loading them.\n",
    "\n",
    "All the main codes and some run results are in `exploration` branch. These functions are copied from this file:\n",
    "\n",
    "https://github.com/Duchstf/quench-detector/blob/exploration/data-visualization-1.ipynb\n",
    "\n",
    "Additional documented and clean examples can also befound in signal-analysis branch:\n",
    "\n",
    "https://github.com/Duchstf/quench-detector/tree/signal-analysis\n",
    "\"\"\"\n",
    "\n",
    "def break_tdms_to_files(file_path):\n",
    "    \"\"\"Read the tdms file and break it into several files according to channels in a new directory\n",
    "    in ./data/ according to the ramp name.\n",
    "    \n",
    "    Example usage: break_tdms_to_files(some_file.tdms)\"\"\"\n",
    "    \n",
    "    tdms_file = TdmsFile.read(file_path)\n",
    "    \n",
    "    data_frame = tdms_file.as_dataframe()\n",
    "    \n",
    "    %reset -f in\n",
    "    \n",
    "    os.mkdir(\"./data/\" + file_path[-10:-5])\n",
    "    \n",
    "    for channel in data_frame.columns[:8]:\n",
    "        np.save(\"./data/{}/{}\".format(file_path[-10:-5], channel[-4:-1]), data_frame.loc[:, channel].to_numpy())\n",
    "\n",
    "def calculate_time(dir_path):\n",
    "    \"\"\"Take a tdms file directory path (after the file is broken up), and calculate and output a time file.\n",
    "    \n",
    "    Suppose that the tdms is already broken into npy file at ./data/Ramp4/\n",
    "    Example usage: calculate_time(\"./data/Ramp4/\")\n",
    "    \"\"\"\n",
    "    \n",
    "    #ai7 saved the \n",
    "    ai7 = np.load(dir_path + \"ai7.npy\")\n",
    "    \n",
    "    #Add time axis relative to the time when quench happens\n",
    "    time_range = np.asarray(range(ai7.shape[0]))\n",
    "    \n",
    "    #Center around the max value (quench happens at 0 time)\n",
    "    max_index = np.argmax(ai7)\n",
    "    time_range -= max_index\n",
    "    time_range = time_range.astype('float32')\n",
    "    \n",
    "    #Multiply by datarate\n",
    "    time_range = np.multiply(time_range, 1e-5, out=time_range, casting=\"unsafe\")\n",
    "    \n",
    "    np.save(dir_path + \"time\", time_range)\n",
    "\n",
    "def load_channel_and_time(dir_path, channel):\n",
    "    \"\"\"\n",
    "    Load the channel with time axis associated to it. Note that each tdms in this case has 5 main channels from\n",
    "    ai0 to ai4 according to 5 sensors.\n",
    "    \n",
    "    Suppose that the tdms is already broken into npy file at ./data/Ramp4/\n",
    "    Example usage: channel = load_channel_and_time(\"./data/Ramp4/\", \"ai2\",) \n",
    "    \"\"\"\n",
    "    data_frame = pd.DataFrame(data = {channel: np.load(dir_path + channel + \".npy\"),\n",
    "                                     \"time\": np.load(dir_path + \"time.npy\")})\n",
    "    return data_frame\n",
    "\n",
    "def plot_channel(data, channel, time_range = None):  \n",
    "    \"\"\"\n",
    "    Load a specific channel from the data and then plot it. \n",
    "    \n",
    "    Example usage:\n",
    "    \n",
    "    load_and_plot(\"./data/Ramp28/\", \"ai4\", time_range = (-20e-3, -1e-3))\n",
    "    \"\"\"\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",
    "    data = data[(data[\"time\"] > start) & (data[\"time\"] < end)]\n",
    "    \n",
    "    plt.figure(figsize=(20,2))\n",
    "    plt.plot(data[\"time\"], data[channel])\n",
    "    plt.xlabel(\"Time (s)\")\n",
    "    plt.ylabel(\"Volt\")\n",
    "    #plt.ylim([-1.5,1.5])\n",
    "    #plt.xlim([-650,100])\n",
    "    plt.title(\"Variable {}\".format(channel))\n",
    "    plt.show()\n",
    "    \n",
    "    %reset -f in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. mat processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mat file processing are mainly found in `15T_exploration` branch. The main processing codes can be found here:\n",
    "\n",
    "https://github.com/Duchstf/quench-detector/blob/15T-exploration/data-reading-dev.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "import seaborn as sns\n",
    "from nptdms import TdmsFile #Process ramping file\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_mat_time(ramp_num, dir_path):\n",
    "    \"\"\"\n",
    "    Take a ramp's mat (MATHLAB produced) file and calculate its time axis from the current data \n",
    "    (set maximum current point at 0 and multiply by sampling rate). \n",
    "    \n",
    "    Example usage: calculate_mat_time(\"Ramp_1\", \"15T_data\")\n",
    "    \n",
    "    \n",
    "    Args:\n",
    "        - ramp_num (str): ramp number\n",
    "        - dir_path (str): path to directory containing all the ramp's data\n",
    "        \n",
    "    Returns: None, save the calculated time axis to a npy file in dir_path + \"/time_axis/\"\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",
    "    current_data = data[3, :]\n",
    "    current_data = np.nan_to_num(current_data)\n",
    "    \n",
    "    #Add time axis relative to the time when quench happens\n",
    "    time_range = np.asarray(range(current_data.shape[0]))\n",
    "    \n",
    "    #Calculate its time axis, center around the max value (quench happens at 0 time)\n",
    "    max_index = np.argmax(current_data)\n",
    "    time_range -= max_index\n",
    "    time_range = time_range.astype('float32')\n",
    "    \n",
    "    #Multiply by datarate\n",
    "    time_range = np.multiply(time_range, 1e-5, out=time_range, casting=\"unsafe\")\n",
    "    \n",
    "    #Save the time axis\n",
    "    np.save(dir_path + \"/time_axis/\" + ramp_num + \"_time\", time_range)\n",
    "    \n",
    "    %reset -f in\n",
    "\n",
    "def plot_sensor(ramp_num, dir_path, sensor):\n",
    "    \n",
    "    \"\"\"\n",
    "    Plot a sensor from mat data\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",
    "    \n",
    "    Returns: None, just plot the sensor's data\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",
    "    plt.figure(figsize = (20,2))\n",
    "    plt.plot(time_axis, data[sensor, :])\n",
    "    plt.xlabel(\"Time [s]\")\n",
    "    plt.ylabel(\"Signal [V]\")\n",
    "    \n",
    "    %reset -f in\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",
    "    #Return the sensor's data accoring to the index \n",
    "    return data[sensor, :][selection_index]\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]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Miscellaneous"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These codes were not needed in the study. They are used to process voltage tap data files. Detailed usage can be found here:\n",
    "\n",
    "https://github.com/Duchstf/quench-detector/blob/exploration/Exploration-2.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_data_from_single_file(filepath):\n",
    "    \"\"\"Read in quench data from a given file, return a pandas dafa frame\"\"\"\n",
    "    \n",
    "    data_dict= {}\n",
    "    \n",
    "    with open(filepath) as f:\n",
    "        content = f.readlines()\n",
    "    #Remove`\\n` at the end of each line\n",
    "    content = [x.strip() for x in content]\n",
    "    \n",
    "    column_names = content[0].split(\" \")\n",
    "    data = [content[i].split(\"   \") for i in range(1, len(content))] \n",
    "    \n",
    "    for i in range(len(column_names)):\n",
    "        data_dict[column_names[i]] = [float(x[i]) for x in data]\n",
    "        \n",
    "    data_frame = pd.DataFrame(data_dict)\n",
    "    \n",
    "    return data_frame\n",
    "\n",
    "def read_quench_data(area_path, quench_name):\n",
    "    \"\"\"\n",
    "    Read the quench data from the provided file path, note that there are five quenches so we need to concatenate them together.\n",
    "    Just need to provide the quench's name, for e.g: \"./data/mqxfs1b.Quench.161011115654\"\n",
    "    \"\"\"\n",
    "    \n",
    "    data_list = []\n",
    "    \n",
    "    \n",
    "    num_file = 0 #Number of separate data files for this quench, assuming at least one\n",
    "    \n",
    "    for filename in os.listdir(area_path):\n",
    "        if filename.startswith(quench_name) and not filename.endswith(\".tar.gz\"):\n",
    "            #print(\"Reading file ... \" + filename)\n",
    "            if num_file == 0:\n",
    "                try:\n",
    "                    data_list.append(read_data_from_single_file(area_path + filename)) #read first file separately\n",
    "                except:\n",
    "                    pass\n",
    "            else:\n",
    "                try:\n",
    "                    data_list.append(read_data_from_single_file(area_path + filename).drop(columns = [\"time\"])) #Time is already saved in the first file\n",
    "                except:\n",
    "                    pass\n",
    "            num_file += 1\n",
    "                                 \n",
    "    #print(\"Total number of data files for {}: \".format(quench_name), num_file)\n",
    "    \n",
    "    if len(data_list) != 0:\n",
    "        data = pd.concat(data_list, axis = 1) #Concatenate all data files together\n",
    "        return data\n",
    "    else:\n",
    "        print(\"Returning None due to file errors in \" + quench_name)\n",
    "        return None\n",
    "\n",
    "def read_all_quench_in_area(area_path):\n",
    "    \"\"\"\n",
    "    Read all quench's files in the area and return a dictionary of different quench's data. name is the magnet's name\n",
    "    \"\"\"\n",
    "    data = {}\n",
    "    \n",
    "    index = 1 #For indexing different quench\n",
    "\n",
    "    for filename in os.listdir(area_path):\n",
    "        if filename.endswith(\".tar.gz\"): \n",
    "            #print(\"Reading data from ... \" + filename[0:-7])\n",
    "            data[filename[0:-7]] = read_quench_data(area_path, filename[0:-7])\n",
    "            index += 1\n",
    "            \n",
    "    #print(\"Completed. There are {} quenches in total.\".format(str(index - 1)))\n",
    "    \n",
    "    return data\n",
    "\n",
    "def plot_variables_with_time(data, variables = None, time_range = None):\n",
    "    \"\"\"\n",
    "    Take a data frame, and plot all other variables with time.\n",
    "    Optional argument: time_range to specify the range to plot, default is to plot all time. Example argument\n",
    "    is (start_time, end_time)\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",
    "    data = data[(data[\"time\"] > start) & (data[\"time\"] < end)]\n",
    "    \n",
    "    if variables:\n",
    "        for variable in variables:\n",
    "            plt.figure(figsize=(20,2))\n",
    "            plt.plot(data[\"time\"], data[variable])\n",
    "            plt.xlabel(\"Time (s)\")\n",
    "            plt.title(\"Variable {}\".format(variable))\n",
    "            plt.show()\n",
    "        \n",
    "    else:      \n",
    "        #If variables not specified then just plot all \n",
    "        for variable in data.columns[1:]:\n",
    "            plt.figure(figsize=(25,2))\n",
    "            plt.plot(data[\"time\"], data[variable])\n",
    "            plt.xlabel(\"Time (s)\")\n",
    "            plt.title(\"Variable {}\".format(variable))\n",
    "            plt.show()\n",
    "\n",
    "def plot_statistics(data, variables = None, time_range = None, window = 100):\n",
    "    \"\"\"Take a data frame and plot the variable with its moving average and real data in the specified range, if\n",
    "    variables are not specified then just plot all variables by default\"\"\"\n",
    "    #Pickout the data\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",
    "    data = data[(data[\"time\"] > start) & (data[\"time\"] < end)]\n",
    "    \n",
    "    def plot_signal_and_statistics(variable):\n",
    "        # Plotted by calculating Simple Moving Average (SMA)\n",
    "        plt.figure(figsize=(20,2))\n",
    "        plt.plot(data[\"time\"], data[variable], label = \"Signal with noise\")\n",
    "        plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).mean(),label = \"Moving average\")\n",
    "        #plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).min(),label = \"Moving min\")\n",
    "        #plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).max(),label = \"Moving max\")\n",
    "        plt.legend(loc = \"best\")\n",
    "        plt.xlabel(\"Time (s)\")\n",
    "        plt.title(\"Variable {}\".format(variable))\n",
    "        plt.show()\n",
    "        \n",
    "        \n",
    "        #Variance\n",
    "        plt.figure(figsize=(20,2))\n",
    "        plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).std(),\n",
    "                 label = \"Variance\",\n",
    "                 color = 'red')\n",
    "        plt.xlabel(\"Time (s)\")\n",
    "        plt.title(\"Variable {}'s Variance\".format(variable))\n",
    "        plt.show()\n",
    "        \n",
    "        #Kurtosis\n",
    "        plt.figure(figsize=(20,2))\n",
    "        plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).kurt(),\n",
    "                 label = \"Kurtosis\",\n",
    "                 color = 'green')\n",
    "        plt.xlabel(\"Time (s)\")\n",
    "        plt.title(\"Variable {}'s Kurtosis\".format(variable))\n",
    "        plt.show()\n",
    "        \n",
    "        #Skew\n",
    "        plt.figure(figsize=(20,2))\n",
    "        plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).skew(),\n",
    "                 label = \"Skew\",\n",
    "                 color = 'purple')\n",
    "        plt.xlabel(\"Time (s)\")\n",
    "        plt.title(\"Variable {}'s Skew\".format(variable))\n",
    "        plt.show()\n",
    "        \n",
    "        \"\"\"\n",
    "        #Quantile\n",
    "        plt.figure(figsize=(20,2))\n",
    "        plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).quantile(0.25), label = \"Quantile 25\",\n",
    "                 color = 'red')\n",
    "        plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).quantile(0.75), label = \"Quantile 75\",\n",
    "                 color = 'orange')\n",
    "        plt.plot(data[\"time\"], data.loc[:, variable].rolling(window=window).quantile(0.50), label = \"Quantile 50\",\n",
    "                 color = 'green')\n",
    "        plt.legend(loc = \"best\")\n",
    "        plt.title(\"Variable {}'s Quantile\".format(variable))\n",
    "        plt.show()\n",
    "        \"\"\"\n",
    "        \n",
    "    #Plot\n",
    "    if variables:\n",
    "        for variable in variables:\n",
    "            print(variable)\n",
    "            plot_signal_and_statistics(variable)         \n",
    "    else:\n",
    "        for variable in data.columns[1:]:\n",
    "            print(variable)\n",
    "            plot_signal_and_statistics(variable)\n",
    "    \n",
    "    %reset -f in\n",
    "\n",
    "def plot_variable(data_dict, variable, time_range = None, window = 100):\n",
    "    for quench_name in data_dict.keys():\n",
    "        print(\"Quench's index: \" + quench_name)\n",
    "        if data_dict[quench_name] is not None:\n",
    "            if variable[0] in list(data_dict[quench_name].columns):\n",
    "                plot_statistics(data_dict[quench_name], variables = variable, time_range=time_range, window=window)"
   ]
  }
 ],
 "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
}