.ipynb_checkpoints/anomaly-detection-md-conversion-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C Conversion\n",
    "---\n",
    "All we need to do to convert our Mahalanobis Distance model to C is to save the mean vector and inverse of the covariance matrix as constant arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os import listdir\n",
    "from os.path import join\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python 3.7.6\n",
      "Numpy 1.18.1\n",
      "SciPy 1.4.1\n"
     ]
    }
   ],
   "source": [
    "# Print versions\n",
    "!python --version\n",
    "print('Numpy ' + np.__version__)\n",
    "print('SciPy ' + sp.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Settings (define file locations)\n",
    "models_path = 'models'  # Where we can find the model files (relative path location)\n",
    "md_file_name = 'md_model'   # Mahalanobis Distance model arrays (.npz will be added)\n",
    "c_model_name = 'md_model'   # Will be given .h suffix\n",
    "\n",
    "c_mean_name = 'model_mu'\n",
    "c_inv_cov_name = 'model_inv_cov'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.01204995 0.00709583 0.00497362]\n",
      "[[5.36191246e-07 8.13786731e-08 3.46141514e-08]\n",
      " [8.13786731e-08 4.76039917e-07 1.12038078e-08]\n",
      " [3.46141514e-08 1.12038078e-08 3.15563520e-07]]\n"
     ]
    }
   ],
   "source": [
    "# Load model\n",
    "with np.load(join(models_path, md_file_name) + '.npz') as data:\n",
    "    model_mu = data['model_mu']\n",
    "    model_cov = data['model_cov']\n",
    "print(model_mu)\n",
    "print(model_cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1927196.25782373 -324748.88969416 -199864.163168  ]\n",
      " [-324748.88969416 2157143.90610544  -40965.82067725]\n",
      " [-199864.163168    -40965.82067725 3192311.64978127]]\n"
     ]
    }
   ],
   "source": [
    "# Calculate inverse of covariance matrix (as it's constant)\n",
    "inv_cov = sp.linalg.inv(model_cov)\n",
    "print(inv_cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to convert an array into a C string (requires Numpy)\n",
    "def create_c_array(np_array, var_type, var_name, line_limit=80, indent=4):\n",
    "    \n",
    "    c_str = ''\n",
    "    \n",
    "    # Add array shape\n",
    "    for i, dim in enumerate(np_array.shape):\n",
    "        c_str += 'const unsigned int ' + var_name + '_dim' + str(i + 1) + ' = ' + str(dim) + ';\\n'\n",
    "    c_str += '\\n'\n",
    "        \n",
    "    # Declare C variable\n",
    "    c_str += 'const ' + var_type + ' ' + var_name\n",
    "    for dim in np_array.shape:\n",
    "        c_str += '[' + str(dim) + ']'\n",
    "    c_str += ' = {\\n'\n",
    "    \n",
    "    # Create string for the array\n",
    "    indent = ' ' * indent\n",
    "    array_str = indent\n",
    "    line_len = len(indent)\n",
    "    val_sep = ', '\n",
    "    array = np_array.flatten()\n",
    "    for i, val in enumerate(array):\n",
    "        \n",
    "        # Create a new line if string is over line limit\n",
    "        val_str = str(val)\n",
    "        if line_len + len(val_str) + len(val_sep) > line_limit:\n",
    "            array_str += '\\n' + indent\n",
    "            line_len = len(indent)\n",
    "\n",
    "        # Add value and separator\n",
    "        array_str += val_str\n",
    "        line_len += len(val_str)\n",
    "        if (i + 1) < len(array):\n",
    "            array_str += val_sep\n",
    "            line_len += len(val_sep)\n",
    "\n",
    "    # Add closing brace\n",
    "    c_str += array_str + '\\n};\\n'\n",
    "        \n",
    "    return c_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to create a header file with given C code as a string\n",
    "def create_c_header(c_code, name):\n",
    "    \n",
    "    c_str = ''\n",
    "    \n",
    "    # Create header guard\n",
    "    c_str += '#ifndef ' + name.upper() + '_H\\n'\n",
    "    c_str += '#define ' + name.upper() + '_H\\n\\n'\n",
    "    \n",
    "    # Add provided code\n",
    "    c_str += c_code\n",
    "    \n",
    "    # Close out header guard\n",
    "    c_str += '\\n#endif //' + name.upper() + '_H'\n",
    "    \n",
    "    return c_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create constant C arrays for model (mu and inverse covariance)\n",
    "c_mu = create_c_array(model_mu, 'float', c_mean_name)\n",
    "c_inv_cov = create_c_array(inv_cov, 'float', c_inv_cov_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#ifndef MD_MODEL_H\n",
      "#define MD_MODEL_H\n",
      "\n",
      "const unsigned int model_mu_dim1 = 3;\n",
      "\n",
      "const float model_mu[3] = {\n",
      "    0.012049953994017095, 0.00709583131025641, 0.004973618240170936\n",
      "};\n",
      "\n",
      "const unsigned int model_inv_cov_dim1 = 3;\n",
      "const unsigned int model_inv_cov_dim2 = 3;\n",
      "\n",
      "const float model_inv_cov[3][3] = {\n",
      "    1927196.2578237287, -324748.88969416154, -199864.16316800224, \n",
      "    -324748.88969416154, 2157143.906105444, -40965.82067724897, \n",
      "    -199864.16316800224, -40965.82067724897, 3192311.649781269\n",
      "};\n",
      "\n",
      "#endif //MD_MODEL_H\n"
     ]
    }
   ],
   "source": [
    "# Construct header file\n",
    "header_str = create_c_header(c_mu + '\\n' + c_inv_cov, c_model_name)\n",
    "print(header_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save C header file\n",
    "with open(join(models_path, c_model_name) + '.h', 'w') as file:\n",
    "    file.write(header_str)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Develop C functions\n",
    "---\n",
    "We don't have access to an easy/efficient library for calculating the Mahalanobis Distance in C, so we need to write a few functions from scratch. It's easiest to get the algorithms down in Python first so we can compare the output from a few samples to known good algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set some settings to test out our functions\n",
    "sample_file_path = 'test-samples'\n",
    "sample_file_name = 'normal_anomaly_samples'  # Will be given .npz suffix\n",
    "sensor_sample_rate = 200    # Hz\n",
    "sample_time = 0.64           # Time (sec) length of each sample\n",
    "max_measurements = int(sample_time * sensor_sample_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(200, 3)\n",
      "(200, 3)\n"
     ]
    }
   ],
   "source": [
    "# Load test samples\n",
    "with np.load(join(sample_file_path, sample_file_name) + '.npz') as data:\n",
    "    normal_sample = data['normal_sample']\n",
    "    anomaly_sample = data['anomaly_sample']\n",
    "print(normal_sample.shape)\n",
    "print(anomaly_sample.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.01158355 0.00651454 0.00434402]\n"
     ]
    }
   ],
   "source": [
    "# Test extracting features (median absolute deviation) using SciPy\n",
    "sample = normal_sample[0:max_measurements]                  # Truncate to 128 measurements\n",
    "sp_test_features = stats.median_absolute_deviation(sample)  # Calculate MAD\n",
    "print(sp_test_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find median in an array (we'll use qsort in C)\n",
    "# We want to calculate MAD in a C-like fashion to make it easier to port\n",
    "def median(arr):\n",
    "    n = len(arr)\n",
    "    sorted_arr = np.sort(arr)\n",
    "    if n % 2 == 0:\n",
    "        median = (sorted_arr[int((n - 1) / 2)] + sorted_arr[int(n / 2)]) / 2.0\n",
    "    else:\n",
    "        median = devs[n / 2]\n",
    "    return median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate MAD manually\n",
    "def calc_mad(arr):\n",
    "    \n",
    "    # Find number of elements\n",
    "    n = len(arr)\n",
    "    \n",
    "    # Get median of array\n",
    "    med = median(arr)\n",
    "\n",
    "    # Calculate absolute deviation from median for each element\n",
    "    devs = [0] * n\n",
    "    for i in range(n):\n",
    "        devs[i] = abs(arr[i] - med)\n",
    "    \n",
    "    # Find median of deviations\n",
    "    return median(devs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wrapper to extract MAD features for all axis\n",
    "def extract_mad_features(sample, scale=1.4826):\n",
    "    mads = []\n",
    "    \n",
    "    # Calculate MAD for each axis (note: SciPy scales by 1.4826)\n",
    "    # https://github.com/scipy/scipy/issues/11090\n",
    "    for axis in range(sample.shape[1]):\n",
    "        mads.append(scale * calc_mad(sample[:, axis]))\n",
    "        \n",
    "    return mads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0115835538, 0.006514544399999972, 0.0043440179999999825]\n"
     ]
    }
   ],
   "source": [
    "# Test MAD algorithm and compare to SciPy answer\n",
    "sample = normal_sample[0:max_measurements]                  # Truncate to 128 measurements\n",
    "normal_x = extract_mad_features(sample)\n",
    "print(normal_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0900788347195736\n"
     ]
    }
   ],
   "source": [
    "# Test: calculate Mahalanobis distance using Numpy\n",
    "x_minus_mu = normal_x - model_mu\n",
    "left_term = np.dot(x_minus_mu, inv_cov)\n",
    "mahal = np.dot(left_term, x_minus_mu.T)\n",
    "print(mahal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the dot product of two vectors\n",
    "def dot_product(a, b):\n",
    "    sum = 0\n",
    "    for i in range(len(a)):\n",
    "        sum += a[i] * b[i]\n",
    "        \n",
    "    return sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do matrix multiplication\n",
    "def matrix_multiply(a, b):\n",
    "    \n",
    "    # Find number of rows and columns\n",
    "    a_rows = len(a)\n",
    "    a_cols = len(a[0])\n",
    "    b_rows = len(b)\n",
    "    b_cols = len(b[0])\n",
    "    \n",
    "    # Check to make sure we can multiply the matrices\n",
    "    if a_cols != b_rows:\n",
    "        print('Dimension mismatch')\n",
    "        return [[]]\n",
    "    \n",
    "    # Create return matrix\n",
    "    prod = [[0 for i in range(b_cols)] for j in range(a_rows)]\n",
    "    \n",
    "    # Calculate the dot product for each element in the return matrix\n",
    "    for i in range(a_rows):\n",
    "        for j in range(b_cols):\n",
    "            for k in range(a_cols):\n",
    "                prod[i][j] += a[i][k] * b[k][j]\n",
    "    \n",
    "    return prod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Mahalanobis distance the C way (with const inverse covariance)\n",
    "def mahalanobis(x, mu, inv_cov):\n",
    "    \n",
    "    # Subtract each element in X from the mean\n",
    "    x_minus_mu = [0 for i in range(len(x))]\n",
    "    for i in range(len(x)):\n",
    "        x_minus_mu[i] = x[i] - mu[i]\n",
    "        \n",
    "    # Compute dot product of prev term and inverse covariance\n",
    "    left_term = [dot_product(x_minus_mu, inv_cov)]\n",
    "    \n",
    "    # Transpose difference matrix\n",
    "    x_minus_mu_t = [[i] for i in x_minus_mu]\n",
    "    \n",
    "    # Matrix multiply prev term and difference\n",
    "    mahal = matrix_multiply(left_term, x_minus_mu_t)\n",
    "    \n",
    "    return mahal[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0900788347195736"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test Mahalanobis function with normal sample\n",
    "mahalanobis(normal_x, model_mu, inv_cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4844.40534566512"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Find MAD for anomaly sample\n",
    "sample = anomaly_sample[0:max_measurements]  # Truncate to 128 measurements\n",
    "anomaly_x = extract_mad_features(sample)\n",
    "mahalanobis(anomaly_x, model_mu, inv_cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}