{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from collections import deque\n",
    "import scipy.stats as stats\n",
    "from statsmodels.tsa.stattools import adfuller\n",
    "from scipy.signal import savgol_filter\n",
    "from math import floor\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.lines import Line2D\n",
    "\n",
    "from BOCD import bocd_meanNstd, NormalUnKnownMeanPrecision\n",
    "\n",
    "\n",
    "def HotellingT2(window):\n",
    "    '''\n",
    "    Find indexes of anomaly using Hotelling T-square method\n",
    "    calculate T2 score for each point & compare with calculated Upper-confidence-level.\n",
    "    if exceed, anomaly\n",
    "\n",
    "    :param window: data list\n",
    "    :return: indexes of anomaly\n",
    "    '''\n",
    "    alpha = 0.01  #FIXME : must be larger value\n",
    "    p = 1\n",
    "    m = len(window)\n",
    "    q = 2 * (m - 1) ** 2 / (3 * m - 4)\n",
    "\n",
    "    UCL = (m - 1) ** 2 * stats.beta.isf(alpha, p / 2, (q - p - 1) / 2) / m\n",
    "\n",
    "    mean = np.mean(window)\n",
    "    V = np.array([])\n",
    "    T2_list = []\n",
    "\n",
    "    for ind in range(m - 1):\n",
    "        V = np.append(V, window[ind + 1] - window[ind])\n",
    "\n",
    "    S = np.array([0.5 * V.transpose() @ V / (m - 1)])\n",
    "\n",
    "    for item in window:\n",
    "        delta = np.array(item) - np.array(mean)\n",
    "        T2 = delta * np.linalg.inv(np.array([S])) * delta\n",
    "        T2_list.append(T2)\n",
    "\n",
    "    anomaly = []\n",
    "    for ind, value in enumerate(T2_list):\n",
    "        if value > UCL:\n",
    "            anomaly.append(ind)\n",
    "\n",
    "    return anomaly\n",
    "\n",
    "class Operator(object):    \n",
    "    def __init__(self, Nw=15, Nomin=3, Nomax=20, Nstd=15):\n",
    "        '''\n",
    "        :param Nw: window size for outlier\n",
    "        :param Nomin: min window size for waiting Mounter-FeedBack (MFB) reflection\n",
    "        :param Nomax: max window size for waiting Mounter-FeedBack (MFB) reflection\n",
    "        :param Nstd: window size for updating std.\n",
    "        '''\n",
    "        # status : two status considered i.e. 'noFBapplied' & 'FBapplied'\n",
    "        self.status = 'noFBapplied'\n",
    "\n",
    "        # window size\n",
    "        self.Nw = Nw  # mentioned\n",
    "        self.Nomin = Nomin  # mentioned\n",
    "        self.Nomax = Nomax  # mentioned\n",
    "        self.Nstd = Nstd  # mentioned\n",
    "\n",
    "        # data list\n",
    "        self.total_data = []  # total data list (FB applied list : collected data list if online production)\n",
    "        self.op_list = [0]  # indexes where operation were applied (not considered delay operation list)\n",
    "        self.CPD_list = [0]  # indexes of detected change-points\n",
    "        self.bocd_f_ind_list = [0]  # indexes of CPD fired points\n",
    "        self.sub_op_list = [0]  # indexes where the sub-operation was applied\n",
    "\n",
    "        # data window\n",
    "        self.window = deque([], maxlen=self.Nw)  # Raw data list for calculating offset & for smoothing (Nw data)\n",
    "        self.window_for_compare = deque([], maxlen=self.Nw)  # Smoothed data list to compare threshold & current state\n",
    "\n",
    "        # Target & Threshold\n",
    "        self.Target = 0\n",
    "        self.offset_threshold = None\n",
    "\n",
    "        # Operation\n",
    "        self.operation_cnt = 0  # number of operation\n",
    "        self.sub_operation = False  # where or not already applied sub-operation, as sub-op applied only once right after CPD\n",
    "\n",
    "        # waiting step for reflection\n",
    "        self.wait_num = 0  # counting number for monitoring Mounter-FeedBack delay\n",
    "        self.abrupt_num = 0  # number counting low p-value point for Abrupt CPD\n",
    "        self.p = 0.05  # respective lower & upper p-value for Abrupt CPD\n",
    "        self.estimator = None  # Gaussian KDE for Abrupt CPD => p-value calculated based on the Gaussian KDE regressed distribution\n",
    "\n",
    "    def _run_CPD(self):\n",
    "        '''\n",
    "        Check if change has occurred in the latest (No + Nw) data\n",
    "        :return: True/False, sequence point of CP\n",
    "        '''\n",
    "        # BOCD 가동\n",
    "        TimeSequence = self.Nw + len(self.total_data) - self.op_list[-1]  # window size for detecting change-point\n",
    "\n",
    "        std = np.std(np.array(self.total_data[-TimeSequence:][:self.Nw]))  # Std calculation for BOCPD FIXME (:self.Nw)\n",
    "\n",
    "        mu0 = 0  # Prior on Gaussian mean. (A parameter of normal-gamma prior distribution)\n",
    "        gamma0 = 1  # (A parameter of normal-gamma prior distribution)\n",
    "        alpha0 = 1  # (A parameter of normal-gamma prior distribution)\n",
    "        beta0 = alpha0 * std ** 2  # \"sqrt(beta/alpha) = std\" (A parameter of normal-gamma prior distribution)\n",
    "\n",
    "        hazard = 1 / 50.0  # Hazard survival function assumes probability to be a CP for each data point. #FIXME (1/20.0 ??)\n",
    "        message = np.array([1])  # Iterative message calculated using previously collected data.\n",
    "\n",
    "        Data_list = []\n",
    "        RL_dist_list = []\n",
    "        Temp_RL_dist = np.zeros((TimeSequence, TimeSequence))\n",
    "\n",
    "        model = NormalUnKnownMeanPrecision(mu0, gamma0, alpha0, beta0)  # Data assumed be to a normal distribution.\n",
    "        # Case : Both of mean and std unknown\n",
    "        # BOCPD class called.\n",
    "\n",
    "        RL = []  # Run-Length distribution initial list\n",
    "        Start = []  # Estimated Change-Point starting point\n",
    "        cp_list = [0]  # Detected Change-Point sequence saved list\n",
    "\n",
    "        for ind, cont in enumerate(range(TimeSequence)):\n",
    "            Data_list.append(self.total_data[-TimeSequence + ind])\n",
    "            RL_dist, new_joint = bocd_meanNstd(data_list=Data_list, model=model,\n",
    "                                               hazard=hazard, Message=message)\n",
    "\n",
    "            message = new_joint  # each point sequential likelihood, or numerator of RL posterior distribution.\n",
    "            RL_dist_list.append(RL_dist)\n",
    "            Temp_RL_dist[ind, :ind + 1] = RL_dist[1:]  # RL distribution temporary saved as a element of a list\n",
    "\n",
    "            RL.append(np.argmax(\n",
    "                Temp_RL_dist[ind]))  # The highest probability corresponding sequence of RL posterior distribution\n",
    "            Start.append((ind, ind - RL[-1]))\n",
    "\n",
    "            if ind >= 1:\n",
    "                if Start[-1][1] != Start[-2][1]:\n",
    "                    if max(cp_list) > ind - RL[-1]:\n",
    "                        cp_list = cp_list[:-1]  # Remove the lastly added Change-Point\n",
    "                    else:\n",
    "                        cp_list.append(ind - RL[-1])  # Change-point Detected and append into the list\n",
    "\n",
    "        if not cp_list:  # If no FeedBack candidate found\n",
    "            return False, None\n",
    "        elif cp_list[-1] < self.Nw:  # Exclude a CP ahead of FeedBack operation applied sequence\n",
    "            return False, None\n",
    "        else:\n",
    "            for cp in cp_list[1:]:  # Exclude initial point (the point '0')\n",
    "                self.CPD_list.append(len(self.total_data) - TimeSequence + cp_list[-1])\n",
    "            return True, cp_list[-1]\n",
    "\n",
    "    def _is_within_threshold(self):\n",
    "        '''\n",
    "        Check if the smoothed data is inside the threshold\n",
    "        :return: True/ False ('True' if data lies within threshold )\n",
    "        '''\n",
    "        is_withinThreshold = False\n",
    "\n",
    "        # Outlier : Last Nw smoothed data must be the same sign\n",
    "        '''\n",
    "        check = 0\n",
    "        for smoothed in self.window_for_compare:\n",
    "            check += np.sign(smoothed)\n",
    "\n",
    "        if np.abs(check) != self.Nw:\n",
    "            is_withinThreshold = True  \n",
    "        '''\n",
    "        check = np.sign(self.window_for_compare[0])\n",
    "        for smoothed in self.window_for_compare:\n",
    "            if(check != np.sign(smoothed)):\n",
    "                is_withinThreshold = True; break; #FIXME (Non-same-sign which is anomaly ??)\n",
    "\n",
    "        # Outlier : Last Nw smoothed data must be out of threshold.\n",
    "        for smoothed in self.window_for_compare:\n",
    "            if np.abs(smoothed - self.Target) <= self.offset_threshold:\n",
    "                is_withinThreshold = True\n",
    "\n",
    "        return is_withinThreshold\n",
    "\n",
    "    def _do_operation(self):\n",
    "        '''\n",
    "        Apply operation & Set operation delay\n",
    "        OUTPUT : operation\n",
    "        '''\n",
    "        self.operation_cnt += 1  # number of operations\n",
    "        self.sub_operation = False\n",
    "\n",
    "        # Exclude anomaly before calculating offset by using Hotelling's T2\n",
    "        windowForAnomaly = self.window\n",
    "        anomaly = HotellingT2(windowForAnomaly)\n",
    "        for delete_ind, anomaly_ind in enumerate(anomaly):\n",
    "            del windowForAnomaly[anomaly_ind - delete_ind] # windowForAnomaly : Anomaly excluded window \n",
    "\n",
    "        # Calculate offset & operation\n",
    "        offset = np.mean(windowForAnomaly)\n",
    "        \n",
    "        # Append operation index\n",
    "        self.op_list.append(len(self.total_data))\n",
    "        \n",
    "        return self.Target - offset\n",
    "\n",
    "\n",
    "    def _do_sub_operation(self):\n",
    "        '''\n",
    "        Apply sub-operation & Set sub-operation delay\n",
    "        OUTPUT : operation\n",
    "        '''\n",
    "        self.sub_operation = True  # number of sub operations\n",
    "\n",
    "        # Exclude anomaly before calculating offset by using Hotelling's T2\n",
    "        windowForAnomaly = self.total_data[max(self.CPD_list[-1], self.op_list[-1]):]\n",
    "        anomaly = HotellingT2(windowForAnomaly)\n",
    "        for delete_ind, anomaly_ind in enumerate(anomaly):\n",
    "            del windowForAnomaly[anomaly_ind - delete_ind] # windowForAnomaly : Anomaly excluded window \n",
    "\n",
    "        # Calculate offset & operation\n",
    "        offset = np.mean(windowForAnomaly)\n",
    "        \n",
    "        # Append sub-operation index\n",
    "        self.sub_op_list.append(len(self.total_data))\n",
    "        \n",
    "        return self.Target - offset\n",
    "\n",
    "\n",
    "    def step(self, feature):\n",
    "        '''\n",
    "        MFB Loop\n",
    "        :param feature: offset data\n",
    "        :return: None\n",
    "        '''\n",
    "        # START\n",
    "        if self.status == 'noFBapplied':\n",
    "            self._put_feature(feature)\n",
    "\n",
    "            # Enough data ?\n",
    "            if len(self.total_data) - self.CPD_list[-1] < self.Nw:\n",
    "                return 0, None\n",
    "\n",
    "            # Check if lies within calculated threshold. If not, do operation.\n",
    "            if self._is_within_threshold():\n",
    "                # If state is on target & stable, do one sseub-operation.\n",
    "                if not self.sub_operation: # \"False\" if no sub-operation from last CPD point. (As only one sub-op applied for one ADF-test steady-state)\n",
    "                    check_point = max(self.CPD_list[-1], self.op_list[-1]) #'op_list':not considered delay operation index list, 'CPD_list': BOCPD found CPD index list\n",
    "                    if (len(self.total_data) - check_point) >= self.Nw:  # There must be enough data after Op. or CPD.\n",
    "                        if adfuller(self.total_data[check_point:])[1] < 1e-6:  # Stable?\n",
    "                            return self._do_sub_operation(), 'SUB'\n",
    "                return 0, None\n",
    "            else:\n",
    "                self.estimator = stats.gaussian_kde(self.window)  # Set Gaussian KDE for Abrupt CPD.   # FIXME : minimum of delay given, then calculate the estimator\n",
    "                self.status = 'FBapplied'\n",
    "                return self._do_operation(), 'MAIN'  # SG-filter smoothing is not consider anomaly # FIXME\n",
    "                \n",
    "\n",
    "        # Waiting for FB to be applied.\n",
    "        if self.status == 'FBapplied':\n",
    "            self._put_feature(feature)\n",
    "            self.wait_num += 1\n",
    "\n",
    "            # Abrupt CPD : repeating low p-value point or waiting enough fires BOCD\n",
    "            if self.estimator.integrate_box_1d(-np.Inf, self.total_data[-1]) < self.p or \\\n",
    "                    self.estimator.integrate_box_1d(-np.Inf, self.total_data[-1]) > 1 - self.p:\n",
    "                self.abrupt_num += 1\n",
    "\n",
    "            if self.wait_num >= self.Nomax:\n",
    "                self._run_CPD()  # CPD fired\n",
    "                self.bocd_f_ind_list.append(len(self.total_data)-1)\n",
    "                self.wait_num = 0  # Reset self.wait_num\n",
    "                self.abrupt_num = 0  # Reset self.abrupt_num\n",
    "                self.status = 'noFBapplied'  # Reset self.status\n",
    "                return 0, None\n",
    "\n",
    "            elif self.abrupt_num >= self.Nomin:\n",
    "                self._run_CPD()  # CPD fired\n",
    "                self.bocd_f_ind_list.append(len(self.total_data)-1)\n",
    "                self.wait_num = 0  # Reset self.wait_num\n",
    "                self.abrupt_num = 0  # Reset self.abrupt_num\n",
    "                self.status = 'noFBapplied'  # Reset self.status\n",
    "                return 0, None\n",
    "\n",
    "            else:\n",
    "                return 0, None\n",
    "\n",
    "            \n",
    "    def show_param(self):\n",
    "        '''\n",
    "        Show parameters\n",
    "        :return: None\n",
    "        '''\n",
    "        print('\\n'+'='*40)\n",
    "        print('\\n>> Component :')\n",
    "        print('\\n>> Parameters :')\n",
    "        print('\\tNw =', self.Nw, ', Nomin =', self.Nomin, ', Nomax =', self.Nomax, ', Nstd =', self.Nstd)\n",
    "        print('\\tTarget Mean =', self.Target)\n",
    "        print('\\tOffset Threshold = 0.7 sigma')\n",
    "        print('\\n'+'='*40)\n",
    "\n",
    "    def draw(self):  # might be removed\n",
    "        '''\n",
    "        Draw figures\n",
    "        :return: None\n",
    "        '''\n",
    "        # Plot total data & operation index & cp index\n",
    "        fit, (ax1) = plt.subplots(nrows=2, ncols=1, figsize=(20, 8))\n",
    "        ax1[0].grid(True)\n",
    "        for ind in self.op_list[1:]:\n",
    "            ax1[0].axvline(ind, linestyle='--', color='g', linewidth=1)\n",
    "        for ind in self.sub_op_list[1:]:\n",
    "            ax1[0].axvline(ind, linestyle=':', color='r')\n",
    "        ax1[0].plot(self.total_data)\n",
    "        ax1[0].axhline(0, linestyle='--', color='grey')\n",
    "        legend_elements = [Line2D([0], [0], color='green', lw=2, label='Main Operation'),\n",
    "                           Line2D([0], [0], color='red', lw=2, label='Sub Operation'),\n",
    "                           Line2D([0], [0], color='blue', lw=2, label='After Delay')]\n",
    "        ax1[0].legend(handles=legend_elements)\n",
    "\n",
    "        # Plot total data & detected cp index\n",
    "        ax1[1].plot(self.total_data)\n",
    "        ax1[1].grid(True)\n",
    "        ax1[1].axhline(0, linestyle='--', color='grey')\n",
    "        for ind in self.op_list[1:]:\n",
    "            ax1[1].axvline(ind - self.Nw, linestyle='--', color='black', linewidth=1)\n",
    "        for ind in self.bocd_f_ind_list[1:]:\n",
    "            ax1[1].axvline(ind, linestyle='--', color='grey', linewidth=1)\n",
    "        for ind in self.CPD_list[1:]:\n",
    "            ax1[1].axvline(ind, linestyle='--', color='r', linewidth=1)\n",
    "        legend_elements = [Line2D([0], [0], color='red', lw=2, label='CP After Op.')]\n",
    "        ax1[1].legend(handles=legend_elements)\n",
    "\n",
    "        \n",
    "    def _put_feature(self, feature):\n",
    "        '''\n",
    "        Get offset data\n",
    "        :param feature: offset data\n",
    "        :return: None\n",
    "        '''\n",
    "        # Append feature to data lists\n",
    "        self.total_data.append(feature) # Collected total data\n",
    "        self.window.append(feature)     # Nw corresponding window size\n",
    "\n",
    "        # Initialize the threshold\n",
    "        if len(self.total_data) == self.Nstd:   # FIXME use once at the initial stage (put outside of the state stage)\n",
    "            self.offset_threshold = 0.7 * np.std(self.total_data)\n",
    "\n",
    "        # Data smoothing\n",
    "        if len(self.total_data) >= self.Nw:\n",
    "            savgol_result = savgol_filter(self.total_data[-self.Nw:], 2*floor((self.Nw-1)/2)+1, 3)\n",
    "            self.window_for_compare = savgol_result  ## compare threshold with sg-filter regressed line\n",
    "\n",
    "        # Reset the threshold\n",
    "        if (len(self.total_data) - self.CPD_list[-1]) % self.Nstd == 0:\n",
    "            if 0.7 * self.offset_threshold > np.std(buffer) or self.offset_threshold < 0.7 * np.std(buffer):\n",
    "                self.offset_threshold = 0.7 * np.std(buffer)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "class MounterFeedback(object):\n",
    "    \n",
    "    def __init__(self):\n",
    "        self.COMPONENT_INFO = []  # Tuple: (BLOCK_NO, REFERENCE)\n",
    "        self.operator_df = None  # DataFrame: rows: components / columns: Axis\n",
    "        self.operation_df = None  # DataFrame: rows: components / columns: Axis\n",
    "        self.feature_df = None  # DataFrame: rows: components / columns: Axis\n",
    "        \n",
    "        self.max_feedback_num = 20\n",
    "        \n",
    "    def ready(self, json_data):\n",
    "        '''\n",
    "        Initialize MFB algorithm & Ready for feedback\n",
    "        :param json_data\n",
    "        :return: None\n",
    "        '''\n",
    "        # Init operator_df & operation_df\n",
    "        operators = {}\n",
    "        operations = {}\n",
    "        features = {}\n",
    "        for ind in range(len(json_data['COMPONENT'][:5])):\n",
    "            BLOCK_NO = json_data['COMPONENT'][ind]['BLOCK_NO']\n",
    "            REFERENCE = json_data['COMPONENT'][ind]['REFERENCE']\n",
    "            \n",
    "            self.COMPONENT_INFO.append((BLOCK_NO,REFERENCE))\n",
    "            \n",
    "            operator_pair = {}\n",
    "            operation_pair = {}\n",
    "            feature_pair = {}\n",
    "            for AXIS in ['LENGTH','WIDTH']:\n",
    "                operator_pair[AXIS] = Operator(Nw=15, Nomin=3, Nomax=20, Nstd=15)\n",
    "                operation_pair[AXIS] = None\n",
    "                feature_pair[AXIS] = None\n",
    "                                       \n",
    "            operators[(BLOCK_NO,REFERENCE)] = operator_pair\n",
    "            operations[(BLOCK_NO,REFERENCE)] = operation_pair\n",
    "            features[(BLOCK_NO,REFERENCE)] = feature_pair\n",
    "            \n",
    "        self.operator_df = pd.DataFrame(operators).T\n",
    "        self.operation_df = pd.DataFrame(operations).T\n",
    "        self.feature_df = pd.DataFrame(features).T\n",
    "    \n",
    "    def _data_reader(self, json_data): # data form will be changed\n",
    "        '''\n",
    "        Save feature as dataframe from json data\n",
    "        :param json_data\n",
    "        :return: None   \n",
    "        '''\n",
    "        for ind, component_info in enumerate(self.COMPONENT_INFO[:5]):\n",
    "            BLOCK_NO = component_info[0]\n",
    "            REFERENCE = component_info[1]\n",
    "            \n",
    "            for AXIS in ['LENGTH','WIDTH']:      \n",
    "                self.feature_df[AXIS][(BLOCK_NO,REFERENCE)] = float(\n",
    "                    json_data['COMPONENT'][ind]['INSP_TYPE_VALUE'][0]['PADOVERHANG'][AXIS])\\\n",
    "                    +np.random.normal(10,10)  #FIXME: noize term should be deleted later\n",
    "                \n",
    "    def _ranking_discriminator(self):\n",
    "        '''\n",
    "        Determine the severity ranking of outlier\n",
    "        :return: (BLOCK_NO,REFERENCE) of high-ranked components\n",
    "        '''\n",
    "        return self.COMPONENT_INFO\n",
    "    \n",
    "    def step(self, json_data):\n",
    "        '''\n",
    "        param json_data\n",
    "        return: feedback information as dictionary\n",
    "        '''\n",
    "        self._data_reader(json_data)\n",
    "        \n",
    "        for ind, component_info in enumerate(self.COMPONENT_INFO):\n",
    "            for AXIS in ['LENGTH','WIDTH']:\n",
    "                self.operation_df[AXIS][component_info] \\\n",
    "                 = self.operator_df[AXIS][component_info].step(self.feature_df[AXIS][component_info])\n",
    "             \n",
    "        #############################\n",
    "        ### ranking_discriminator ###\n",
    "\n",
    "        high_ranked_component = self._ranking_discriminator()\n",
    "        \n",
    "        #############################\n",
    "        \n",
    "        FEEDBACK_INFO = {\"FeedbackType\": \"PLACEDLOCATIONOFFSET\",\n",
    "                         \"SeqID\": json_data[\"ASST_SEQNO\"],\n",
    "                         \"Time\": json_data[\"END_DATE_TIME\"],\n",
    "                         \"MachineName\": json_data[\"MACHINE_NM\"],\n",
    "                         \"LaneNo\": json_data[\"LANE_NO\"],\n",
    "                         \"Target\": '1',\n",
    "                         \"NumList\": len(json_data[\"COMPONENT\"]),\n",
    "                         \"OFFSETLIST\": []}\n",
    "        \n",
    "        for ind, component_info in enumerate(self.COMPONENT_INFO):\n",
    "            if component_info in high_ranked_component:\n",
    "                FEEDBACK_INFO['OFFSETLIST'].append({'BLOCK_NO': component_info[0],\n",
    "                                                    'REFERENCE': component_info[1],\n",
    "                                                    'AdjustX': self.operation_df['LENGTH'][component_info][0],\n",
    "                                                    'AdjustY': self.operation_df['WIDTH'][component_info][0]})\n",
    "            else:\n",
    "                FEEDBACK_INFO['OFFSETLIST'].append({'BLOCK_NO': component_info[0],\n",
    "                                                    'REFERENCE': component_info[1],\n",
    "                                                    'AdjustX': 0,\n",
    "                                                    'AdjustY': 0})\n",
    "        return FEEDBACK_INFO\n",
    "                "
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
