Statistical_Change_Detection.ipynb

{
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.10.13",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "kaggle": {
      "accelerator": "none",
      "dataSources": [
        {
          "sourceId": 8344135,
          "sourceType": "datasetVersion",
          "datasetId": 4956394
        },
        {
          "sourceId": 8344150,
          "sourceType": "datasetVersion",
          "datasetId": 4956406
        }
      ],
      "dockerImageVersionId": 30698,
      "isInternetEnabled": true,
      "language": "python",
      "sourceType": "notebook",
      "isGpuEnabled": false
    },
    "colab": {
      "provenance": [],
      "include_colab_link": true
    }
  },
  "nbformat_minor": 5,
  "nbformat": 4,
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Mayurd9/KDE-test/blob/main/Statistical_Change_Detection.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "from scipy.stats import norm\n",
        "import matplotlib.pyplot as plt"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:53:59.772664Z",
          "iopub.execute_input": "2024-05-07T11:53:59.773077Z",
          "iopub.status.idle": "2024-05-07T11:53:59.778574Z",
          "shell.execute_reply.started": "2024-05-07T11:53:59.773046Z",
          "shell.execute_reply": "2024-05-07T11:53:59.777459Z"
        },
        "trusted": true,
        "id": "-UWnrccDmvhP"
      },
      "execution_count": 64,
      "outputs": [],
      "id": "-UWnrccDmvhP"
    },
    {
      "cell_type": "markdown",
      "source": [
        "# **$\\textbf{All the required functions}$**"
      ],
      "metadata": {
        "id": "B9wEdoVmmvhW"
      },
      "id": "B9wEdoVmmvhW"
    },
    {
      "cell_type": "code",
      "source": [
        "def small_p(xi, xj, cov):\n",
        "    k_dim = s1_matrix.shape[1]          #dimensions of s1 matrix\n",
        "    cov_1 = cov\n",
        "    det_cov = np.linalg.det(cov)\n",
        "    epsilon = 1e-2\n",
        "    det_cov = max(det_cov, epsilon)\n",
        "\n",
        "    if det_cov == 0.01:\n",
        "        alpha = 1e-1\n",
        "        cov_1 = cov_1 + alpha * np.identity(s1_matrix.shape[1])\n",
        "        det_cov = np.linalg.det(cov_1)        #new addition 4/5/2024\n",
        "\n",
        "    a = 1 / (((2 * np.pi) ** (k_dim / 2)) * det_cov ** 0.5)\n",
        "    b = -0.5 * ((xj - xi).T @ np.linalg.inv(cov_1) @ (xj - xi))\n",
        "\n",
        "    if b > 1000:\n",
        "        print('b value', b)\n",
        "    b = min(b, 50)\n",
        "    p = a * np.exp(b)\n",
        "\n",
        "    return p"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:00.395273Z",
          "iopub.execute_input": "2024-05-07T11:54:00.395687Z",
          "iopub.status.idle": "2024-05-07T11:54:00.406121Z",
          "shell.execute_reply.started": "2024-05-07T11:54:00.395656Z",
          "shell.execute_reply": "2024-05-07T11:54:00.404135Z"
        },
        "trusted": true,
        "id": "wsufGMqGmvhY"
      },
      "execution_count": 65,
      "outputs": [],
      "id": "wsufGMqGmvhY"
    },
    {
      "cell_type": "code",
      "source": [
        "def capital_p(i, j, cov_mat):\n",
        "    summation = 0\n",
        "    xi = s1_matrix[i, :].reshape(-1, 1)\n",
        "    xj = s1_matrix[j, :].reshape(-1, 1)\n",
        "\n",
        "    for k in range(len_s1):\n",
        "        if k != j:\n",
        "            xk = s1_matrix[k, :].reshape(-1, 1)\n",
        "            summation += small_p(xk,xj, cov_mat[k])\n",
        "\n",
        "    # Avoid division by zero\n",
        "    summation = np.clip(summation,1e-4,np.inf)\n",
        "    P = small_p(xi,xj,cov_mat[i]) / summation\n",
        "\n",
        "    return  P"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:00.608000Z",
          "iopub.execute_input": "2024-05-07T11:54:00.609069Z",
          "iopub.status.idle": "2024-05-07T11:54:00.616389Z",
          "shell.execute_reply.started": "2024-05-07T11:54:00.609032Z",
          "shell.execute_reply": "2024-05-07T11:54:00.615300Z"
        },
        "trusted": true,
        "id": "UT0nMQ92mvhZ"
      },
      "execution_count": 66,
      "outputs": [],
      "id": "UT0nMQ92mvhZ"
    },
    {
      "cell_type": "code",
      "source": [
        "def cov_xi(i,cov_mat):\n",
        "    denominator = 0\n",
        "    numerator = 0\n",
        "\n",
        "    xi = s1_matrix[i, :].reshape(-1, 1)\n",
        "    for j in range(len_s1):\n",
        "        if i!=j:\n",
        "            xj = s1_matrix[j, :].reshape(-1, 1)\n",
        "            P = capital_p(i, j, cov_mat)\n",
        "            denominator += P\n",
        "            numerator += P * ((xj - xi) @ (xj - xi).T)\n",
        "\n",
        "    return numerator / denominator"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:00.812170Z",
          "iopub.execute_input": "2024-05-07T11:54:00.812565Z",
          "iopub.status.idle": "2024-05-07T11:54:00.819730Z",
          "shell.execute_reply.started": "2024-05-07T11:54:00.812534Z",
          "shell.execute_reply": "2024-05-07T11:54:00.818440Z"
        },
        "trusted": true,
        "id": "5ExddSUQmvha"
      },
      "execution_count": 67,
      "outputs": [],
      "id": "5ExddSUQmvha"
    },
    {
      "cell_type": "code",
      "source": [
        "def L(cov_mat):\n",
        "    L_new = 0\n",
        "    for j in range(len_s1):\n",
        "        xj = s1_matrix[j, :].reshape(-1, 1)\n",
        "        summation = 0\n",
        "        for i in range(len_s1):\n",
        "            if i!=j:\n",
        "                xi = s1_matrix[i, :].reshape(-1, 1)\n",
        "                summation+=small_p(xi,xj,cov_mat[i])\n",
        "\n",
        "        summation=summation/(len_s1-1)\n",
        "        L_new+=np.log(summation)\n",
        "\n",
        "    return L_new[0][0]"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:01.021488Z",
          "iopub.execute_input": "2024-05-07T11:54:01.022128Z",
          "iopub.status.idle": "2024-05-07T11:54:01.029162Z",
          "shell.execute_reply.started": "2024-05-07T11:54:01.022070Z",
          "shell.execute_reply": "2024-05-07T11:54:01.027993Z"
        },
        "trusted": true,
        "id": "ZNKRoPbdmvha"
      },
      "execution_count": 68,
      "outputs": [],
      "id": "ZNKRoPbdmvha"
    },
    {
      "cell_type": "code",
      "source": [
        "def Algorithm1(cov_mat):\n",
        "    t = 0\n",
        "    phi = 0.01\n",
        "    L_prev = -250\n",
        "\n",
        "    while t < 41:\n",
        "        temp = []\n",
        "        for i in range(len_s1):\n",
        "            temp.append(cov_xi(i,cov_mat))\n",
        "\n",
        "        cov_mat = temp.copy()\n",
        "\n",
        "        for i in range(len_s1):\n",
        "            if np.linalg.det(cov_mat[i]) <= 0:\n",
        "                alpha = 1e-2\n",
        "                cov_mat[i] = cov_mat[i] + alpha * np.identity(s1_matrix.shape[1])\n",
        "\n",
        "        L_new = L(cov_mat)\n",
        "\n",
        "        print('L_prev=',L_prev)\n",
        "        print('L_new=',L_new)\n",
        "        print('t=',t)\n",
        "        print('-'*50)\n",
        "\n",
        "        if abs((L_new-L_prev)/L_prev) < phi:\n",
        "            break\n",
        "\n",
        "        else:\n",
        "            L_prev=L_new\n",
        "            t+=1\n",
        "\n",
        "    return cov_mat"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:01.234110Z",
          "iopub.execute_input": "2024-05-07T11:54:01.234651Z",
          "iopub.status.idle": "2024-05-07T11:54:01.245501Z",
          "shell.execute_reply.started": "2024-05-07T11:54:01.234610Z",
          "shell.execute_reply": "2024-05-07T11:54:01.244296Z"
        },
        "trusted": true,
        "id": "88-3G1i_mvhb"
      },
      "execution_count": 69,
      "outputs": [],
      "id": "88-3G1i_mvhb"
    },
    {
      "cell_type": "code",
      "source": [
        "def F(T):\n",
        "    s = 0\n",
        "    for i in range(len_s1):\n",
        "        xi = s1_matrix[i, :].reshape(-1, 1)\n",
        "        s += small_p(xi,T,cov_mat[i])\n",
        "    f = np.log(s/len_s1)\n",
        "    return f\n",
        "\n",
        "def EstVar(estsize,beta):\n",
        "    Est = []\n",
        "    for t in range(1,estsize):\n",
        "        #R = s2.sample(n=len(s2), replace=True, random_state=42)\n",
        "        R=s2_matrix       #Specified for image dataset, in case of non-image dataset comment this and uncomment above and below lines\n",
        "        #R = R.to_numpy()\n",
        "        F_R = []\n",
        "        for i in range(len_s2):\n",
        "            T = R[i, :].reshape(-1, 1)\n",
        "            F_R.append(F(T))\n",
        "\n",
        "        var_F_R = np.var(np.array(F_R))\n",
        "        Est.append((len_s2/(len_s2-1))*var_F_R)\n",
        "\n",
        "    V =  np.percentile(Est, estsize*(1 - beta))\n",
        "    var_delta = (len_sprime + ((len_sprime**2) / len_s2))*V\n",
        "    return var_delta"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:01.439085Z",
          "iopub.execute_input": "2024-05-07T11:54:01.439495Z",
          "iopub.status.idle": "2024-05-07T11:54:01.450568Z",
          "shell.execute_reply.started": "2024-05-07T11:54:01.439464Z",
          "shell.execute_reply": "2024-05-07T11:54:01.449189Z"
        },
        "trusted": true,
        "id": "cNvKy2YFmvhb"
      },
      "execution_count": 70,
      "outputs": [],
      "id": "cNvKy2YFmvhb"
    },
    {
      "cell_type": "code",
      "source": [
        "from scipy.stats import norm\n",
        "\n",
        "def critical_value(p,stepsize):\n",
        "    estsize = 100\n",
        "    M = int(np.ceil((p/stepsize)-1))\n",
        "    C = []\n",
        "    for i in range(1, M+1):\n",
        "        alpha_i = i*stepsize\n",
        "        beta_i = p-alpha_i\n",
        "        var_delta = EstVar(estsize,beta_i)\n",
        "        mean = 0\n",
        "        std_dev = np.sqrt(var_delta)\n",
        "        c = norm.ppf(alpha_i, loc = mean, scale = std_dev)\n",
        "        C.append(c)\n",
        "\n",
        "    C_max = max(C)\n",
        "    return C_max"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:01.709494Z",
          "iopub.execute_input": "2024-05-07T11:54:01.709920Z",
          "iopub.status.idle": "2024-05-07T11:54:01.718633Z",
          "shell.execute_reply.started": "2024-05-07T11:54:01.709886Z",
          "shell.execute_reply": "2024-05-07T11:54:01.717193Z"
        },
        "trusted": true,
        "id": "h2XJOcStmvhc"
      },
      "execution_count": 71,
      "outputs": [],
      "id": "h2XJOcStmvhc"
    },
    {
      "cell_type": "code",
      "source": [
        "def delta():\n",
        "\n",
        "    a = 0\n",
        "    for j in range(len_sprime):\n",
        "        y = sprime_matrix[j, :].reshape(-1, 1)\n",
        "        summation = 0\n",
        "        for i in range(len_s1):\n",
        "            x = s1_matrix[i, :].reshape(-1, 1)\n",
        "            summation += small_p(x,y,cov_mat[i])\n",
        "\n",
        "        a += np.log(summation/len_s1)\n",
        "\n",
        "    b = 0\n",
        "    for j in range(len_s1):\n",
        "        y = s1_matrix[j, :].reshape(-1, 1)\n",
        "        summation = 0\n",
        "        for i in range(len_s1):\n",
        "            x = s1_matrix[i, :].reshape(-1, 1)\n",
        "            summation += small_p(x,y,cov_mat[i])\n",
        "\n",
        "        b += (np.log(summation/len_s1))*len_sprime/len_s1\n",
        "\n",
        "    delta = a-b\n",
        "\n",
        "    return delta[0][0]"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:01.977480Z",
          "iopub.execute_input": "2024-05-07T11:54:01.978082Z",
          "iopub.status.idle": "2024-05-07T11:54:01.987442Z",
          "shell.execute_reply.started": "2024-05-07T11:54:01.978050Z",
          "shell.execute_reply": "2024-05-07T11:54:01.986392Z"
        },
        "trusted": true,
        "id": "nqHjpfDKmvhc"
      },
      "execution_count": 72,
      "outputs": [],
      "id": "nqHjpfDKmvhc"
    },
    {
      "cell_type": "markdown",
      "source": [
        "# $\\textbf{Data Pre-processing}$"
      ],
      "metadata": {
        "id": "DQ53u1Rjmvhd"
      },
      "id": "DQ53u1Rjmvhd"
    },
    {
      "cell_type": "code",
      "source": [
        "import tensorflow as tf\n",
        "from tensorflow.keras.datasets import mnist\n",
        "\n",
        "# Load the MNIST dataset\n",
        "(x_train, y_train), (x_test, y_test) = mnist.load_data()\n",
        "\n",
        "# Normalize the pixel values to be between 0 and 1\n",
        "x_train, x_test = x_train / 255.0, x_test / 255.0\n",
        "\n",
        "# Print the shape of the training and testing data\n",
        "print(\"Training data shape:\", x_train.shape)\n",
        "print(\"Training labels shape:\", y_train.shape)\n",
        "print(\"Testing data shape:\", x_test.shape)\n",
        "print(\"Testing labels shape:\", y_test.shape)"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:54:04.644523Z",
          "iopub.execute_input": "2024-05-07T11:54:04.644926Z",
          "iopub.status.idle": "2024-05-07T11:54:05.182229Z",
          "shell.execute_reply.started": "2024-05-07T11:54:04.644896Z",
          "shell.execute_reply": "2024-05-07T11:54:05.180825Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QHjlpaUTmvhd",
        "outputId": "8a27265e-1163-49b1-8740-445529534b63"
      },
      "execution_count": 73,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Training data shape: (60000, 28, 28)\n",
            "Training labels shape: (60000,)\n",
            "Testing data shape: (10000, 28, 28)\n",
            "Testing labels shape: (10000,)\n"
          ]
        }
      ],
      "id": "QHjlpaUTmvhd"
    },
    {
      "cell_type": "markdown",
      "source": [
        "# $\\textbf{Case 1: When S and S' are entirely different images}$"
      ],
      "metadata": {
        "id": "wcUY2zU92HGR"
      },
      "id": "wcUY2zU92HGR"
    },
    {
      "cell_type": "code",
      "source": [
        "image_index_s1=50100    #we can put any number here\n",
        "image_index_sprime= image_index_s1 + np.random.randint(1,10)   #Check whether these two images are different or not by running next cell. If they are not same, re-run this cell.\n",
        "\n",
        "s1_matrix=x_train[image_index_s1]\n",
        "len_s1=s1_matrix.shape[0]\n",
        "s2_matrix=s1_matrix.copy()\n",
        "len_s2=s2_matrix.shape[0]\n",
        "sprime_matrix=x_train[image_index_sprime]\n",
        "len_sprime=sprime_matrix.shape[0]"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:59:10.717871Z",
          "iopub.execute_input": "2024-05-07T11:59:10.718432Z",
          "iopub.status.idle": "2024-05-07T11:59:10.725755Z",
          "shell.execute_reply.started": "2024-05-07T11:59:10.718397Z",
          "shell.execute_reply": "2024-05-07T11:59:10.724360Z"
        },
        "trusted": true,
        "id": "0AXkrtaEmvhe"
      },
      "execution_count": 74,
      "outputs": [],
      "id": "0AXkrtaEmvhe"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Displaying images}$**"
      ],
      "metadata": {
        "id": "Pdilntsfmvhf"
      },
      "id": "Pdilntsfmvhf"
    },
    {
      "cell_type": "code",
      "source": [
        " #visualizing images\n",
        "s1_matrix_view=s1_matrix*255\n",
        "sprime_matrix_view=sprime_matrix*255\n",
        "\n",
        "s1_matrix_view=s1_matrix*255\n",
        "sprime_matrix_view=sprime_matrix*255\n",
        "\n",
        "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n",
        "axes[0].imshow(s1_matrix_view, cmap='gray')\n",
        "axes[0].set_title('S1')\n",
        "axes[0].axis('on')\n",
        "\n",
        "axes[1].imshow(sprime_matrix_view, cmap='gray')\n",
        "axes[1].set_title('S prime')\n",
        "axes[1].axis('on')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:59:14.048596Z",
          "iopub.execute_input": "2024-05-07T11:59:14.049795Z",
          "iopub.status.idle": "2024-05-07T11:59:14.990441Z",
          "shell.execute_reply.started": "2024-05-07T11:59:14.049755Z",
          "shell.execute_reply": "2024-05-07T11:59:14.989149Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 533
        },
        "id": "LEN0HcAVmvhf",
        "outputId": "3c65c328-702f-45ed-8cbc-5a389e7a3a62"
      },
      "execution_count": 75,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1000x500 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAIECAYAAAAEiE3VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtUklEQVR4nO3dfZSWdZ0/8M/Iw4jKDAHDDCMPogZjorQiTpYZriRgWSjuUSsPmkfNkLPKcStKJauNZLez9kBpbat5Now10bZybZME6wi4wRKxrSRkQfGgaM4oLGBw/f5om9+OPM587y/3PfB6nXOdM/d1X5/r/vD1xg/vue6HqqIoigAAAABK7qhyNwAAAACHK6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACAToRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAuriFCxdGVVVVLFy4sNytAK8jdMMR6Be/+EVceumlMXTo0Dj66KPj+OOPj3e+853xpS99qe2Yf//3f49rrrkmRo4cGd26dYsTTjihfA0DQIU7mNkKHJmqiqIoyt0EcOg89dRTcd5558WQIUNiypQp0dDQEOvXr48lS5bE2rVrY82aNRERcdVVV8W8efPijDPOiHXr1kW3bt3iN7/5TXmbB4AKdLCzNafdu3fHzp07o2fPnnHUUa6rQSURuuEI8653vSv+4z/+I371q19Fnz592t33/PPPx4ABAyIiYsOGDVFXVxc9evSId7/73bFq1SqhGwD24mBnaw7bt28XtKHC+dsJR5i1a9fGqaeeusc/CiKi3T8KGhsbo0ePHoewMwDomg52tu7L2LFjY+TIkbFs2bJ461vfGr169Yphw4bF3Xff3e64P79v+9vf/nbceuutcfzxx8cxxxwTra2te31P95/Pu3LlynjHO94RxxxzTJx88snxne98JyIiFi1aFM3NzdGrV68YMWJEPP7443v09vvf/z4++MEPRn19fVRXV8epp54a//RP/9SxBYIjnNANR5ihQ4fGsmXLYtWqVeVuBQAOC6WYrX/4wx/iwgsvjNGjR8fs2bNj0KBBccMNN+w14H7605+OH/zgB3HLLbfEZz/72ejZs+d+z/vud787mpubY/bs2VFdXR2XX355zJs3Ly6//PK48MIL43Of+1xs3bo1Lr300njllVfaajdv3hxvectb4vHHH48bb7wxvvCFL8TJJ58c11xzTdx1112d/rPCkaZ7uRsADq1bbrklJk6cGG9+85vjrLPOire//e1x/vnnx3nnnefKNgB0Qilm64YNG+Lzn/98TJ8+PSIirr/++mhubo4ZM2bElVde2e4827dvj5/97GfRq1evgzrv3Llz44orroiIiHe+853R1NQU73vf++Kpp56K5ubmiIg45ZRTYvz48fHQQw/FVVddFRERn/jEJ2LXrl3xi1/8Ivr16xcRER/60IfiiiuuiE9+8pNx/fXXH1QPcKRzpRuOMO985ztj8eLF8Z73vCd+/vOfx+zZs2P8+PFx/PHHx7/+67+Wuz0A6HJKMVu7d+8e119/fdvtnj17xvXXXx/PP/98LFu2rN2xU6ZMOeiwe9xxx8Xll1/ednvEiBHRp0+fOOWUU9oCd0S0/fzrX/86IiKKooiHHnooLrrooiiKIrZs2dK2jR8/PlpaWmL58uUH1QMc6YRuOAKNGTMm5s+fH3/4wx/i6aefjhkzZsQrr7wSl156afzyl78sd3sA0OWkztbGxsY49thj2+0bPnx4RMQeH2Q6bNiwg+5r0KBBUVVV1W5fbW1tDB48eI99EX96OXpExAsvvBAvv/xyfO1rX4u6urp229VXXx0Rf/qQOODAvLwcjmA9e/aMMWPGxJgxY2L48OFx9dVXx4MPPhgzZ84sd2sA0CUditnakZd0d+vWrUP7//zFRrt3746IiA984AMxZcqUvR57+umnH3QfcCQTuoGIiDjzzDMjImLjxo1l7gQADg8dma0bNmyIrVu3trva/atf/SoiIk444YQs/e1PXV1d9O7dO3bt2hXjxo075I8PhxMvL4cjzBNPPNH2W+z/69FHH42IP73XCwA4eKWYrX/84x/jnnvuabu9c+fOuOeee6Kuri5Gjx5dumYPUrdu3WLy5Mnx0EMP7fVT2V944YVD3hN0Va50wxFm2rRpsW3btrj44oujqakpdu7cGU899VTMmzcvTjjhhLb3aa1cubLtw1/WrFkTLS0t8ZnPfCYiIkaNGhUXXXRR2f4MAFBJDna27k9jY2Pceeed8Zvf/CaGDx8e8+bNixUrVsTXvva1sn27yOc+97l44oknorm5Oa699tp405veFC+99FIsX748Hn/88XjppZfK0hd0NUI3HGH+/u//Ph588MF49NFH42tf+1rs3LkzhgwZEh/+8Ifj1ltvjT59+kRExPLly+O2225rV/vn21OmTBG6AeB/Hexs3Z83vOEN8c1vfjOmTZsWX//616O+vj6+/OUvx7XXXpv/D7AP9fX18fTTT8enPvWpmD9/fnzlK1+Jfv36xamnnhp33nln2fqCrqaq2NtrYQAAgENi7NixsWXLlr2+jBvo+rynGwAAADIRugEAACAToRsAAAAy8Z5uAAAAyMSVbgAAAMhE6AYAAIBMKu57unfv3h0bNmyI3r17R1VVVbnbAYBDriiKeOWVV6KxsTGOOqpyfz9uZgNwJDvYeV1xoXvDhg0xePDgcrcBAGW3fv36GDRoULnb2CczGwAOPK8r7tfnvXv3LncLAFARKn0mVnp/AHAoHGgeVlzo9vI0APiTSp+Jld4fABwKB5qHFRe6AQAA4HCRLXTPmTMnTjjhhDj66KOjubk5nn766VwPBQB0knkNAHllCd3z5s2L6dOnx8yZM2P58uUxatSoGD9+fDz//PM5Hg4A6ATzGgDyqyqKoij1SZubm2PMmDHx5S9/OSL+9JUigwcPjmnTpsXHPvax/da2trZGbW1tqVsCgC6npaUlampqsp0/ZV5HmNkAEHHgeV3yK907d+6MZcuWxbhx4/7/gxx1VIwbNy4WL168x/E7duyI1tbWdhsAkFdH53WEmQ0AnVHy0L1ly5bYtWtX1NfXt9tfX18fmzZt2uP4WbNmRW1tbdvm+z4BIL+OzusIMxsAOqPsn14+Y8aMaGlpadvWr19f7pYAgL0wswGg47qX+oT9+/ePbt26xebNm9vt37x5czQ0NOxxfHV1dVRXV5e6DQBgPzo6ryPMbADojJJf6e7Zs2eMHj06FixY0LZv9+7dsWDBgjj77LNL/XAAQCeY1wBwaJT8SndExPTp02PKlClx5plnxllnnRV33XVXbN26Na6++uocDwcAdIJ5DQD5ZQndl112Wbzwwgtx++23x6ZNm+LNb35zPPbYY3t8WAsAUD7mNQDkl+V7ulP4zk8A+JPc39OdyswGgDJ8TzcAAADwJ0I3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCYlD92f/OQno6qqqt3W1NRU6ocBABKY1wBwaHTPcdJTTz01Hn/88f//IN2zPAwAkMC8BoD8skzX7t27R0NDQ45TAwAlYl4DQH5Z3tP97LPPRmNjY5x44onx/ve/P9atW5fjYQCABOY1AORXVRRFUcoT/tu//Vu8+uqrMWLEiNi4cWPccccd8fvf/z5WrVoVvXv33uP4HTt2xI4dO9put7a2xuDBg0vZEgB0SS0tLVFTU5Pl3B2d1xFmNgDszQHndZHZH/7wh6Kmpqb4x3/8x73eP3PmzCIibDabzWazvW5raWnJPaYPel6b2TabzWaz7X070LzO/pVhffr0ieHDh8eaNWv2ev+MGTOipaWlbVu/fn3ulgCA1znQvI4wswGgM7KH7ldffTXWrl0bAwcO3Ov91dXVUVNT024DAA6tA83rCDMbADqj5KH7lltuiUWLFsVvfvObeOqpp+Liiy+Obt26xRVXXFHqhwIAOsm8BoBDo+RfGfa73/0urrjiinjxxRejrq4uzjnnnFiyZEnU1dWV+qEAgE4yrwHg0Cj5p5enam1tjdra2nK3AQBll/PTy0vBzAaAA8/r7O/pBgAAgCNVyV9eTrrPf/7zSfU33HBDcg933nlnUv26deuS6n/1q18l1af64x//mFS/dOnSEnVy5Gpqakqq79evX1L9L37xi6T61tbWpHoAAA4PrnQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCbdy90AezrnnHOS6qurq5N7uP3225PPUU5VVVVJ9X/84x+T6tesWZNUX26p6xcRURRFUv3xxx+fVH/ccccl1a9bty6p/pRTTkmq3759e1I9AIdGXV1dUv0HPvCBpPpJkyYl1UdEvP3tb0+qT535c+fOTaq/8sork+ohN1e6AQAAIBOhGwAAADIRugEAACAToRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACAToRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACCT7uVugD1NmTIlqf5b3/pWcg+//e1vk+pPPfXUpPqTTz45qT5Vt27dkupHjBhRok7Ko6qqKvkcRVGUoJPyGTJkSFL9FVdckVR/7733JtUDcHDq6uqS6h999NGk+jPOOCOpvhTzNvUcqfVNTU1J9VDpXOkGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyqiqIoyt3E/9Xa2hq1tbXlbqNLO/roo5PPsX379qT6Y489Nqn+tNNOS6qvr69Pqh8+fHhSPREjR45Mqp80aVJS/XHHHZdUnyr170Dq30EODy0tLVFTU1PuNvbJzKYSnHvuuUn1CxcuTKpP/ad0VVVVUv38+fOT6iMiLr300qT6iy++OKn+oYceSqr/27/926T62267LakeDjSvXekGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyqiqIoyt3E/9Xa2hq1tbXlbgMos4cffjip/j3veU+JOumcbt26lfXxOTy0tLRETU1NudvYJzObUqirq0uqf/TRR5PqzzjjjKT6X/7yl0n1U6ZMSap/5plnkuojIrZt25Z8jhTf+c53kuonTZqUVP/AAw8k1V955ZVJ9XR9B5rXrnQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGTS4dD95JNPxkUXXRSNjY1RVVUVjzzySLv7i6KI22+/PQYOHBi9evWKcePGxbPPPluqfgGAg2BeA0Bl6HDo3rp1a4waNSrmzJmz1/tnz54dX/ziF+Puu++OpUuXxrHHHhvjx4+P7du3JzcLABwc8xoAKkP3jhZMnDgxJk6cuNf7iqKIu+66K2699dZ473vfGxER999/f9TX18cjjzwSl19+eVq3AMBBMa8BoDKU9D3dzz33XGzatCnGjRvXtq+2tjaam5tj8eLFe63ZsWNHtLa2ttsAgHw6M68jzGwA6IyShu5NmzZFRER9fX27/fX19W33vd6sWbOitra2bRs8eHApWwIAXqcz8zrCzAaAzij7p5fPmDEjWlpa2rb169eXuyUAYC/MbADouJKG7oaGhoiI2Lx5c7v9mzdvbrvv9aqrq6OmpqbdBgDk05l5HWFmA0BnlDR0Dxs2LBoaGmLBggVt+1pbW2Pp0qVx9tlnl/KhAIBOMq8B4NDp8KeXv/rqq7FmzZq2288991ysWLEi+vbtG0OGDImbbropPvOZz8Qb3/jGGDZsWNx2223R2NgYkyZNKmXfAMB+mNcAUBk6HLp/9rOfxXnnndd2e/r06RERMWXKlLjvvvviIx/5SGzdujWuu+66ePnll+Occ86Jxx57LI4++ujSdQ0A7Jd5DQCVocOhe+zYsVEUxT7vr6qqik996lPxqU99KqkxAKDzzGsAqAwdDt0AB6N///5J9e9973tL1EnnfOMb3yjr4wMcKT7wgQ8k1Z9xxhlJ9du2bUuqv/3225Pqly9fnlRfCmeeeWZSfVNTU1L9JZdcklS/v18wHozx48cn1Q8ZMiSpft26dUn1VL6yf2UYAAAAHK6EbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMupe7AeDwNHny5KT6oiiS6nft2pVUP2/evKR6AA7OiBEjkupT58VPfvKTstbX1dUl1d9///1J9RERF1xwQVJ96n+Dctf369cvqb5///5J9evWrUuqp/K50g0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmXQvdwPA4am2trasj//ss88m1T/++OMl6gSAnKqqqpLqJ0yYkFS/efPmpPrU/ouiSKovRQ9d/fFvuOGGpPrly5eXqBMOV650AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAm3cvdAFB5jjnmmORzTJkypQSdAHC4+/rXv55UXxRFUv11111X1sdPVe7Hr4QefvnLXybVz58/v0SdwN650g0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkEn3cjcAVJ7LLrss+RxNTU0l6KTzhgwZklT/k5/8pESddE1z585Nqp8/f35yD5s3b04+B1D5li1bllS/devWpPrjjjsuqT513qXOm1L8/3b8+PFJ9R//+MeT6rdt25ZU/1d/9VdJ9Vu2bEmqhwNxpRsAAAAyEboBAAAgE6EbAAAAMulw6H7yySfjoosuisbGxqiqqopHHnmk3f1XXXVVVFVVtdsmTJhQqn4BgINgXgNAZehw6N66dWuMGjUq5syZs89jJkyYEBs3bmzbHnjggaQmAYCOMa8BoDJ0+NPLJ06cGBMnTtzvMdXV1dHQ0NDppgCANOY1AFSGLO/pXrhwYQwYMCBGjBgRN9xwQ7z44ov7PHbHjh3R2trabgMA8uvIvI4wswGgM0oeuidMmBD3339/LFiwIO68885YtGhRTJw4MXbt2rXX42fNmhW1tbVt2+DBg0vdEgDwOh2d1xFmNgB0RodfXn4gl19+edvPp512Wpx++ulx0kknxcKFC+P888/f4/gZM2bE9OnT2263trYa4gCQWUfndYSZDQCdkf0rw0488cTo379/rFmzZq/3V1dXR01NTbsNADi0DjSvI8xsAOiM7KH7d7/7Xbz44osxcODA3A8FAHSSeQ0AeXT45eWvvvpqu9+CP/fcc7FixYro27dv9O3bN+64446YPHlyNDQ0xNq1a+MjH/lInHzyyTF+/PiSNg4A7Jt5DQCVocOh+2c/+1mcd955bbf//N6uKVOmxFe/+tVYuXJlfPOb34yXX345Ghsb44ILLohPf/rTUV1dXbquAYD9Mq8BoDJ0OHSPHTs2iqLY5/0//OEPkxoCANKZ1wBQGbK/pxsAAACOVFXF/n4NXgatra1RW1tb7jY4wvXo0SOp/tJLL02q39fX9Rysd73rXUn1DQ0NSfURsd8rbEeCqqqqpPqdO3cm1T/99NNJ9anuuuuu5HM89NBD6Y10cS0tLRX9CeFmNpRfU1NT8jn+67/+K6k+deb/53/+Z1L9mDFjkuoh1YHmtSvdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJBJ93I3AJWotrY2qf6f//mfS9RJeRRFURHnSPHrX/86qf473/lOUn19fX1S/dKlS5Pq77nnnqR6ALqGT3ziE8nnqKqqSqrfsmVLUv2VV16ZVA+VzpVuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyEbgAAAMike7kbAErvt7/9bVL9CSecUJpGEixdujSpfuzYsUn1O3fuTKoHgINRV1eXVH/OOeck91AURVL9Zz/72aT6Z555JqkeKp0r3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQSfdyNwCV6KWXXkqqHz16dFJ9nz59kuqHDx+eVH/33Xcn1ZfCD37wg6T6nTt3lqgTAMjn05/+dFL9kCFDkntYvnx5Uv0XvvCF5B7gcOZKNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABkInQDAABAJkI3AAAAZCJ0AwAAQCZCNwAAAGQidAMAAEAmQjcAAABk0r3cDUAl2r17d1L9ihUrStNIJ02bNq2sjx8RsX379qT6H/zgByXqBAAq18UXX5xUXxRFcg/z589PPgewb650AwAAQCZCNwAAAGQidAMAAEAmHQrds2bNijFjxkTv3r1jwIABMWnSpFi9enW7Y7Zv3x5Tp06Nfv36xXHHHReTJ0+OzZs3l7RpAGD/zGwAqAwdCt2LFi2KqVOnxpIlS+JHP/pRvPbaa3HBBRfE1q1b2465+eab43vf+148+OCDsWjRotiwYUNccsklJW8cANg3MxsAKkOHPr38sccea3f7vvvuiwEDBsSyZcvi3HPPjZaWlvjGN74Rc+fOjb/8y7+MiIh77703TjnllFiyZEm85S1vKV3nAMA+mdkAUBmS3tPd0tISERF9+/aNiIhly5bFa6+9FuPGjWs7pqmpKYYMGRKLFy9OeSgAIIGZDQDl0env6d69e3fcdNNN8ba3vS1GjhwZERGbNm2Knj17Rp8+fdodW19fH5s2bdrreXbs2BE7duxou93a2trZlgCAvTCzAaB8On2le+rUqbFq1ar49re/ndTArFmzora2tm0bPHhw0vkAgPbMbAAon06F7htvvDG+//3vxxNPPBGDBg1q29/Q0BA7d+6Ml19+ud3xmzdvjoaGhr2ea8aMGdHS0tK2rV+/vjMtAQB7YWYDQHl1KHQXRRE33nhjPPzww/HjH/84hg0b1u7+0aNHR48ePWLBggVt+1avXh3r1q2Ls88+e6/nrK6ujpqamnYbAJDGzAaAytCh93RPnTo15s6dG9/97nejd+/ebe/5qq2tjV69ekVtbW1cc801MX369Ojbt2/U1NTEtGnT4uyzz/YpqABwCJnZAFAZOhS6v/rVr0ZExNixY9vtv/fee+Oqq66KiIh/+Id/iKOOOiomT54cO3bsiPHjx8dXvvKVkjQLABwcMxsAKkOHQndRFAc85uijj445c+bEnDlzOt0UAJDGzAaAypD0Pd0AAADAvnX6e7qBfAYOHJhU/9a3vjWpvqqqKqk+IuKHP/xhUv2KFSuSewCA3D7xiU8k1dfV1SXVv/DCC0n1EREPP/xw8jmAfXOlGwAAADIRugEAACAToRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACAToRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACAToRsAAAAy6V7uBoA9ffCDH0yqHzBgQFJ9URRJ9RERq1evTj4HAOTW1NSUVP+xj30sqT515l555ZVJ9RERzzzzTPI5gH1zpRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACAToRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACAToRsAAAAyEboBAAAgE6EbAAAAMule7gaAPV188cXlbgEAjgjHHntsUv0xxxyTVL9+/fqk+uXLlyfVA/m50g0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkEn3cjcA7Gnz5s3lbiHZAw88UO4WAOCA7r///qT6oiiS6u+5556k+i1btiTVA/m50g0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmXQvdwPAnj7+8Y8n1f/FX/xFUv2zzz6bVB8R8fzzzyefAwAOpKmpKan+TW96U1L97t27k+pfeumlpHqg8rnSDQAAAJkI3QAAAJCJ0A0AAACZdCh0z5o1K8aMGRO9e/eOAQMGxKRJk2L16tXtjhk7dmxUVVW12z70oQ+VtGkAYP/MbACoDB0K3YsWLYqpU6fGkiVL4kc/+lG89tprccEFF8TWrVvbHXfttdfGxo0b27bZs2eXtGkAYP/MbACoDB369PLHHnus3e377rsvBgwYEMuWLYtzzz23bf8xxxwTDQ0NpekQAOgwMxsAKkPSe7pbWloiIqJv377t9n/rW9+K/v37x8iRI2PGjBmxbdu2fZ5jx44d0dra2m4DAErLzAaA8uj093Tv3r07brrppnjb294WI0eObNv/vve9L4YOHRqNjY2xcuXK+OhHPxqrV6+O+fPn7/U8s2bNijvuuKOzbQAAB2BmA0D5dDp0T506NVatWhU//elP2+2/7rrr2n4+7bTTYuDAgXH++efH2rVr46STTtrjPDNmzIjp06e33W5tbY3Bgwd3ti0A4HXMbAAon06F7htvvDG+//3vx5NPPhmDBg3a77HNzc0REbFmzZq9DvDq6uqorq7uTBsAwAGY2QBQXh0K3UVRxLRp0+Lhhx+OhQsXxrBhww5Ys2LFioiIGDhwYKcaBAA6zswGgMrQodA9derUmDt3bnz3u9+N3r17x6ZNmyIiora2Nnr16hVr166NuXPnxoUXXhj9+vWLlStXxs033xznnntunH766Vn+AADAnsxsAKgMHQrdX/3qVyMiYuzYse3233vvvXHVVVdFz5494/HHH4+77rortm7dGoMHD47JkyfHrbfeWrKGAYADM7MBoDJ0+OXl+zN48OBYtGhRUkMAQDozGwAqQ9L3dAMAAAD71umvDAPy+fnPf55U39jYWKJOAKCyPfPMM0n1d999d1L9KaecklT/0EMPJdUDlc+VbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIROgGAACATIRuAAAAyEToBgAAgEyEbgAAAMhE6AYAAIBMhG4AAADIpHu5GwAAgHK54YYbyt0CcJhzpRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACAToRsAAAAyEboBAAAgE6EbAAAAMhG6AQAAIBOhGwAAADIRugEAACAToRsAAAAyEboBAAAgk4oL3UVRlLsFAKgIlT4TK70/ADgUDjQPKy50v/LKK+VuAQAqQqXPxErvDwAOhQPNw6qiwn5NvXv37tiwYUP07t07qqqq9ri/tbU1Bg8eHOvXr4+ampoydNj1WcM01i+dNUxj/dJ0hfUriiJeeeWVaGxsjKOOqrjfj7cxs/OyfumsYRrrl8b6pav0NTzYed39EPZ0UI466qgYNGjQAY+rqampyIXvSqxhGuuXzhqmsX5pKn39amtry93CAZnZh4b1S2cN01i/NNYvXSWv4cHM68r99TkAAAB0cUI3AAAAZNLlQnd1dXXMnDkzqqury91Kl2UN01i/dNYwjfVLY/0OHWudxvqls4ZprF8a65fucFnDivsgNQAAADhcdLkr3QAAANBVCN0AAACQidANAAAAmQjdAAAAkEmXC91z5syJE044IY4++uhobm6Op59+utwtdQmf/OQno6qqqt3W1NRU7rYq2pNPPhkXXXRRNDY2RlVVVTzyyCPt7i+KIm6//fYYOHBg9OrVK8aNGxfPPvtseZqtQAdav6uuumqP5+SECRPK02wFmjVrVowZMyZ69+4dAwYMiEmTJsXq1avbHbN9+/aYOnVq9OvXL4477riYPHlybN68uUwdV56DWcOxY8fu8Tz80Ic+VKaODy/mdeeZ2R1jXqczs9OY2WmOhHndpUL3vHnzYvr06TFz5sxYvnx5jBo1KsaPHx/PP/98uVvrEk499dTYuHFj2/bTn/603C1VtK1bt8aoUaNizpw5e71/9uzZ8cUvfjHuvvvuWLp0aRx77LExfvz42L59+yHutDIdaP0iIiZMmNDuOfnAAw8cwg4r26JFi2Lq1KmxZMmS+NGPfhSvvfZaXHDBBbF169a2Y26++eb43ve+Fw8++GAsWrQoNmzYEJdcckkZu64sB7OGERHXXnttu+fh7Nmzy9Tx4cO8TmdmHzzzOp2ZncbMTnNEzOuiCznrrLOKqVOntt3etWtX0djYWMyaNauMXXUNM2fOLEaNGlXuNrqsiCgefvjhttu7d+8uGhoair/7u79r2/fyyy8X1dXVxQMPPFCGDivb69evKIpiypQpxXvf+96y9NMVPf/880VEFIsWLSqK4k/Ptx49ehQPPvhg2zH//d//XUREsXjx4nK1WdFev4ZFURTveMc7ir/+678uX1OHKfM6jZndeeZ1OjM7nZmd5nCc113mSvfOnTtj2bJlMW7cuLZ9Rx11VIwbNy4WL15cxs66jmeffTYaGxvjxBNPjPe///2xbt26crfUZT333HOxadOmds/H2traaG5u9nzsgIULF8aAAQNixIgRccMNN8SLL75Y7pYqVktLS0RE9O3bNyIili1bFq+99lq752BTU1MMGTLEc3AfXr+Gf/atb30r+vfvHyNHjowZM2bEtm3bytHeYcO8Lg0zuzTM69Ixsw+emZ3mcJzX3cvdwMHasmVL7Nq1K+rr69vtr6+vj2eeeaZMXXUdzc3Ncd9998WIESNi48aNcccdd8Tb3/72WLVqVfTu3bvc7XU5mzZtiojY6/Pxz/exfxMmTIhLLrkkhg0bFmvXro2Pf/zjMXHixFi8eHF069at3O1VlN27d8dNN90Ub3vb22LkyJER8afnYM+ePaNPnz7tjvUc3Lu9rWFExPve974YOnRoNDY2xsqVK+OjH/1orF69OubPn1/Gbrs28zqdmV065nVpmNkHz8xOc7jO6y4TukkzceLEtp9PP/30aG5ujqFDh8a//Mu/xDXXXFPGzjhSXX755W0/n3baaXH66afHSSedFAsXLozzzz+/jJ1VnqlTp8aqVau8pzPBvtbwuuuua/v5tNNOi4EDB8b5558fa9eujZNOOulQtwkRYWZTeczsg2dmpzlc53WXeXl5//79o1u3bnt8yt/mzZujoaGhTF11XX369Inhw4fHmjVryt1Kl/Tn55znY+mceOKJ0b9/f8/J17nxxhvj+9//fjzxxBMxaNCgtv0NDQ2xc+fOePnll9sd7zm4p32t4d40NzdHRHgeJjCvS8/M7jzzOg8ze+/M7DSH87zuMqG7Z8+eMXr06FiwYEHbvt27d8eCBQvi7LPPLmNnXdOrr74aa9eujYEDB5a7lS5p2LBh0dDQ0O752NraGkuXLvV87KTf/e538eKLL3pO/q+iKOLGG2+Mhx9+OH784x/HsGHD2t0/evTo6NGjR7vn4OrVq2PdunWeg//rQGu4NytWrIiI8DxMYF6XnpndeeZ1HmZ2e2Z2miNhXnepl5dPnz49pkyZEmeeeWacddZZcdddd8XWrVvj6quvLndrFe+WW26Jiy66KIYOHRobNmyImTNnRrdu3eKKK64od2sV69VXX23327PnnnsuVqxYEX379o0hQ4bETTfdFJ/5zGfijW98YwwbNixuu+22aGxsjEmTJpWv6Qqyv/Xr27dv3HHHHTF58uRoaGiItWvXxkc+8pE4+eSTY/z48WXsunJMnTo15s6dG9/97nejd+/ebe/5qq2tjV69ekVtbW1cc801MX369Ojbt2/U1NTEtGnT4uyzz463vOUtZe6+MhxoDdeuXRtz586NCy+8MPr16xcrV66Mm2++Oc4999w4/fTTy9x912ZepzGzO8a8TmdmpzGz0xwR87q8H57ecV/60peKIUOGFD179izOOuusYsmSJeVuqUu47LLLioEDBxY9e/Ysjj/++OKyyy4r1qxZU+62KtoTTzxRRMQe25QpU4qi+NPXkNx2221FfX19UV1dXZx//vnF6tWry9t0Bdnf+m3btq244IILirq6uqJHjx7F0KFDi2uvvbbYtGlTuduuGHtbu4go7r333rZj/ud//qf48Ic/XLzhDW8ojjnmmOLiiy8uNm7cWL6mK8yB1nDdunXFueeeW/Tt27eorq4uTj755OJv/uZvipaWlvI2fpgwrzvPzO4Y8zqdmZ3GzE5zJMzrqqIoijxxHgAAAI5sXeY93QAAANDVCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJCJ0A0AAACZCN0AAACQidANAAAAmQjdAAAAkInQDQAAAJkI3QAAAJDJ/wMOQDa5oWRxqgAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ],
      "id": "LEN0HcAVmvhf"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Training}$**"
      ],
      "metadata": {
        "id": "4S4P4LYXmvhf"
      },
      "id": "4S4P4LYXmvhf"
    },
    {
      "cell_type": "code",
      "source": [
        "# Initial guess\n",
        "cov_mat = []\n",
        "for i in range(len_s1):\n",
        "    xi = s1_matrix[i, :].reshape(-1, 1)\n",
        "    mat = 0\n",
        "    for j in range(len_s1):\n",
        "        xj = s1_matrix[j, :].reshape(-1, 1)\n",
        "        mat += (xj-xi) @ (xj-xi).T\n",
        "    cov_mat.append(mat/(len_s1-1))\n",
        "\n",
        "#Training\n",
        "cov_mat= Algorithm1(cov_mat)\n",
        "\n",
        "p = 0.08\n",
        "stepsize = 0.02\n",
        "C_max= critical_value(p,stepsize)"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:59:19.687328Z",
          "iopub.execute_input": "2024-05-07T11:59:19.688652Z",
          "iopub.status.idle": "2024-05-07T12:00:20.093467Z",
          "shell.execute_reply.started": "2024-05-07T11:59:19.688609Z",
          "shell.execute_reply": "2024-05-07T12:00:20.091898Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nP9qZOuymvhf",
        "outputId": "c191009a-48c1-4102-defa-265cceace393"
      },
      "execution_count": 76,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "L_prev= -250\n",
            "L_new= -24.979061419041678\n",
            "t= 0\n",
            "--------------------------------------------------\n",
            "L_prev= -24.979061419041678\n",
            "L_new= -2.7417944214373087\n",
            "t= 1\n",
            "--------------------------------------------------\n",
            "L_prev= -2.7417944214373087\n",
            "L_new= 12.923671832844727\n",
            "t= 2\n",
            "--------------------------------------------------\n",
            "L_prev= 12.923671832844727\n",
            "L_new= 26.498437651226823\n",
            "t= 3\n",
            "--------------------------------------------------\n",
            "L_prev= 26.498437651226823\n",
            "L_new= 36.745719528997235\n",
            "t= 4\n",
            "--------------------------------------------------\n",
            "L_prev= 36.745719528997235\n",
            "L_new= 42.59375909105305\n",
            "t= 5\n",
            "--------------------------------------------------\n",
            "L_prev= 42.59375909105305\n",
            "L_new= 44.80530555362608\n",
            "t= 6\n",
            "--------------------------------------------------\n",
            "L_prev= 44.80530555362608\n",
            "L_new= 45.70452806416698\n",
            "t= 7\n",
            "--------------------------------------------------\n",
            "L_prev= 45.70452806416698\n",
            "L_new= 46.07599296970359\n",
            "t= 8\n",
            "--------------------------------------------------\n"
          ]
        }
      ],
      "id": "nP9qZOuymvhf"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Change detection decision}$**"
      ],
      "metadata": {
        "id": "lT62Efrkmvhg"
      },
      "id": "lT62Efrkmvhg"
    },
    {
      "cell_type": "code",
      "source": [
        "Delta = delta()\n",
        "print('Delta:',Delta)\n",
        "\n",
        "if Delta < C_max:\n",
        "    print('THERE IS CHANGE. DISTRIBUTION OF S PRIME IS DIFFERENT THAN DISTRIBUTION OF S ')\n",
        "\n",
        "else:\n",
        "    print('THERE IS NO CHANGE. DISTRIBUTION OF S PRIME IS SAME AS DISTRIBUTION OF S ')"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T12:00:38.325910Z",
          "iopub.execute_input": "2024-05-07T12:00:38.326338Z",
          "iopub.status.idle": "2024-05-07T12:01:03.488229Z",
          "shell.execute_reply.started": "2024-05-07T12:00:38.326306Z",
          "shell.execute_reply": "2024-05-07T12:01:03.486998Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5pia1do0mvhh",
        "outputId": "e6d515a6-e66a-49e8-b681-f610f76fe73c"
      },
      "execution_count": 77,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Delta: -81.62426690012975\n",
            "THERE IS CHANGE. DISTRIBUTION OF S PRIME IS DIFFERENT THAN DISTRIBUTION OF S \n"
          ]
        }
      ],
      "id": "5pia1do0mvhh"
    },
    {
      "cell_type": "markdown",
      "source": [
        "# $\\textbf{Case 2: When S and S' are same}$"
      ],
      "metadata": {
        "id": "rJl9ZCD24BJc"
      },
      "id": "rJl9ZCD24BJc"
    },
    {
      "cell_type": "code",
      "source": [
        "image_index_s1 = image_index_sprime=14809    #we can put any number here\n",
        "s1_matrix=x_train[image_index_s1]\n",
        "len_s1=s1_matrix.shape[0]\n",
        "s2_matrix=s1_matrix.copy()\n",
        "len_s2=s2_matrix.shape[0]\n",
        "sprime_matrix=x_train[image_index_sprime]\n",
        "len_sprime=sprime_matrix.shape[0]"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:59:10.717871Z",
          "iopub.execute_input": "2024-05-07T11:59:10.718432Z",
          "iopub.status.idle": "2024-05-07T11:59:10.725755Z",
          "shell.execute_reply.started": "2024-05-07T11:59:10.718397Z",
          "shell.execute_reply": "2024-05-07T11:59:10.724360Z"
        },
        "trusted": true,
        "id": "2f2B1t6T4b7w"
      },
      "execution_count": 78,
      "outputs": [],
      "id": "2f2B1t6T4b7w"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Displaying images}$**"
      ],
      "metadata": {
        "id": "atYWeele4b7x"
      },
      "id": "atYWeele4b7x"
    },
    {
      "cell_type": "code",
      "source": [
        " #visualizing images\n",
        "s1_matrix_view=s1_matrix*255\n",
        "sprime_matrix_view=sprime_matrix*255\n",
        "\n",
        "s1_matrix_view=s1_matrix*255\n",
        "sprime_matrix_view=sprime_matrix*255\n",
        "\n",
        "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n",
        "axes[0].imshow(s1_matrix_view, cmap='gray')\n",
        "axes[0].set_title('S1')\n",
        "axes[0].axis('on')\n",
        "\n",
        "axes[1].imshow(sprime_matrix_view, cmap='gray')\n",
        "axes[1].set_title('S prime')\n",
        "axes[1].axis('on')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:59:14.048596Z",
          "iopub.execute_input": "2024-05-07T11:59:14.049795Z",
          "iopub.status.idle": "2024-05-07T11:59:14.990441Z",
          "shell.execute_reply.started": "2024-05-07T11:59:14.049755Z",
          "shell.execute_reply": "2024-05-07T11:59:14.989149Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 533
        },
        "outputId": "35a844ad-c7c0-4853-d749-b1982547806c",
        "id": "-1v4al4l4b7x"
      },
      "execution_count": 79,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1000x500 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAIECAYAAAAEiE3VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAss0lEQVR4nO3de5SU9XkH8GcAWTGyQ0DYZeUialArgRiDqzUxGFFEY7yeIya1aDnGGLRVjqYlVdGYlIS2J1Ylt140OQY11luj1lYR8OQIGrGGcNoSoXrAcvGSOIMYwcLbP2y2rlyW3Xd+vLPs53POew4z8/5mnnkz6zfffXdmSlmWZQEAAADUXK+iBwAAAIC9ldINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAA3dzChQujVCrFwoULix4F+AClG3qgX/7yl3HeeefFyJEjY999940DDzwwTj755Lj11lvb9vnXf/3XmDZtWowZMyZ69+4dBx10UHEDA0Cd251sBXqmUpZlWdFDAHvO008/HSeeeGKMGDEipk6dGs3NzbFmzZpYsmRJrFq1KlauXBkRERdddFHcc8898fGPfzxWr14dvXv3jpdffrnY4QGgDu1utqa0bdu22LJlS/Tt2zd69XJeDeqJ0g09zOmnnx4///nP41e/+lUMGDCg3W2vvvpqDBkyJCIi1q5dG4MHD4599tknPvvZz8by5cuVbgDYgd3N1hTeeecdRRvqnJ9O6GFWrVoVRx555Hb/pyAi2v2fgpaWlthnn3324GQA0D3tbrbuzIQJE2LMmDGxdOnS+P3f//3o169fjBo1Kr73ve+12+9379u+++6749prr40DDzww9ttvv6hWqzt8T/fv7nfZsmXx6U9/Ovbbb7849NBD4x//8R8jImLRokXR2toa/fr1i8MOOyyeeOKJ7Wb77//+7/ijP/qjaGpqioaGhjjyyCPjH/7hHzp3gKCHU7qhhxk5cmQsXbo0li9fXvQoALBXqEW2/uY3v4nTTjstjj766JgzZ04MGzYsLrvssh0W3JtuuikeeeSRuPrqq+Mv/uIvom/fvru8389+9rPR2toac+bMiYaGhpgyZUrcc889MWXKlDjttNPim9/8ZmzatCnOO++82LhxY9vaDRs2xLHHHhtPPPFEXH755fE3f/M3ceihh8a0adPi5ptv7vJzhZ6mT9EDAHvW1VdfHZMnT46Pfexjccwxx8SnPvWpOOmkk+LEE090ZhsAuqAW2bp27dr467/+65gxY0ZERFx66aXR2toaM2fOjAsvvLDd/bzzzjvx3HPPRb9+/XbrfufNmxcXXHBBREScfPLJcfjhh8fnP//5ePrpp6O1tTUiIo444oiYNGlS3HfffXHRRRdFRMSf//mfx9atW+OXv/xlDBo0KCIivvSlL8UFF1wQN9xwQ1x66aW7NQP0dM50Qw9z8sknx+LFi+Nzn/tc/OIXv4g5c+bEpEmT4sADD4x/+qd/Kno8AOh2apGtffr0iUsvvbTtct++fePSSy+NV199NZYuXdpu36lTp+522d1///1jypQpbZcPO+ywGDBgQBxxxBFthTsi2v79X//1XxERkWVZ3HfffXHGGWdElmXx+uuvt22TJk2KSqUSzz///G7NAD2d0g090Pjx4+P++++P3/zmN/Hss8/GzJkzY+PGjXHeeefFv//7vxc9HgB0O3mztaWlJT70oQ+1u2706NEREdt9kOmoUaN2e65hw4ZFqVRqd125XI7hw4dvd13Ee3+OHhHx2muvxZtvvhk/+MEPYvDgwe22iy++OCLe+5A4oGP+vBx6sL59+8b48eNj/PjxMXr06Lj44ovj3nvvjVmzZhU9GgB0S3siWzvzJ929e/fu1PW/+2Kjbdu2RUTEH/zBH8TUqVN3uO/YsWN3ew7oyZRuICIiPvGJT0RExLp16wqeBAD2Dp3J1rVr18amTZvane3+1a9+FRERBx10UJL5dmXw4MHRv3//2Lp1a0ycOHGPPz7sTfx5OfQwCxYsaPst9vs9+uijEfHee70AgN1Xi2z9n//5n/j+97/fdnnLli3x/e9/PwYPHhxHH3107YbdTb17945zzz037rvvvh1+Kvtrr722x2eC7sqZbuhhrrjiinj77bfj7LPPjsMPPzy2bNkSTz/9dNxzzz1x0EEHtb1Pa9myZW0f/rJy5cqoVCrx9a9/PSIixo0bF2eccUZhzwEA6snuZuuutLS0xLe+9a14+eWXY/To0XHPPffECy+8ED/4wQ8K+3aRb37zm7FgwYJobW2NSy65JH7v934vfv3rX8fzzz8fTzzxRPz6178uZC7obpRu6GH+6q/+Ku6999549NFH4wc/+EFs2bIlRowYEV/+8pfj2muvjQEDBkRExPPPPx/XXXddu7W/uzx16lSlGwD+z+5m6658+MMfjh/+8IdxxRVXxN/+7d9GU1NT3HbbbXHJJZekfwI70dTUFM8++2x87Wtfi/vvvz++853vxKBBg+LII4+Mb33rW4XNBd1NKdvR38IAAAB7xIQJE+L111/f4Z9xA92f93QDAABAIko3AAAAJKJ0AwAAQCLe0w0AAACJONMNAAAAiSjdAAAAkEjdfU/3tm3bYu3atdG/f/8olUpFjwMAe1yWZbFx48ZoaWmJXr3q9/fjMhuAnmx387ruSvfatWtj+PDhRY8BAIVbs2ZNDBs2rOgxdkpmA0DHeV13vz7v379/0SMAQF2o90ys9/kAYE/oKA/rrnT78zQAeE+9Z2K9zwcAe0JHeVh3pRsAAAD2FslK99y5c+Oggw6KfffdN1pbW+PZZ59N9VAAQBfJawBIK0npvueee2LGjBkxa9aseP7552PcuHExadKkePXVV1M8HADQBfIaANIrZVmW1fpOW1tbY/z48XHbbbdFxHtfKTJ8+PC44oor4s/+7M92ubZarUa5XK71SADQ7VQqlWhsbEx2/3nyOkJmA0BEx3ld8zPdW7ZsiaVLl8bEiRP//0F69YqJEyfG4sWLt9t/8+bNUa1W220AQFqdzesImQ0AXVHz0v3666/H1q1bo6mpqd31TU1NsX79+u32nz17dpTL5bbN930CQHqdzesImQ0AXVH4p5fPnDkzKpVK27ZmzZqiRwIAdkBmA0Dn9an1HR5wwAHRu3fv2LBhQ7vrN2zYEM3Nzdvt39DQEA0NDbUeAwDYhc7mdYTMBoCuqPmZ7r59+8bRRx8d8+fPb7tu27ZtMX/+/DjuuONq/XAAQBfIawDYM2p+pjsiYsaMGTF16tT4xCc+Ecccc0zcfPPNsWnTprj44otTPBwA0AXyGgDSS1K6zz///Hjttdfi+uuvj/Xr18fHPvaxeOyxx7b7sBYAoDjyGgDSS/I93Xn4zk8AeE/q7+nOS2YDQAHf0w0AAAC8R+kGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARPoUPQD1afjw4bnWv/zyy7UZpCC9euX7fdQrr7ySa/03vvGNXOvXrVuXa/3Pf/7zXOsjItauXZv7PgDomMyW2XnJbEjLmW4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASKRP0QNQnz7+8Y/nWp9lWY0mKca2bdtyrR86dGiu9bfddluu9Xm99tprue/j5JNPzrV++fLluWcA6AlktszOS2ZDWs50AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkUsqyLCt6iPerVqtRLpeLHqPHGzRoUK71S5cuzbV+2LBhudbnVSqVcq2vsx+rTsv7/CMiNm7cmGv9xRdfnGv9Aw88kGs91INKpRKNjY1Fj7FTMrs+yGyZnZfMhnw6ymtnugEAACARpRsAAAASUboBAAAgkZqX7htuuCFKpVK77fDDD6/1wwAAOchrANgz+qS40yOPPDKeeOKJ/3+QPkkeBgDIQV4DQHpJ0rVPnz7R3Nyc4q4BgBqR1wCQXpL3dL/44ovR0tISBx98cHzhC1+I1atXp3gYACAHeQ0A6dX8THdra2vccccdcdhhh8W6devixhtvjE996lOxfPny6N+//3b7b968OTZv3tx2uVqt1nokAOADOpvXETIbALqi5qV78uTJbf8eO3ZstLa2xsiRI+MnP/lJTJs2bbv9Z8+eHTfeeGOtxwAAdqGzeR0hswGgK5J/ZdiAAQNi9OjRsXLlyh3ePnPmzKhUKm3bmjVrUo8EAHxAR3kdIbMBoCuSl+633norVq1aFUOHDt3h7Q0NDdHY2NhuAwD2rI7yOkJmA0BX1Lx0X3311bFo0aJ4+eWX4+mnn46zzz47evfuHRdccEGtHwoA6CJ5DQB7Rs3f0/3KK6/EBRdcEG+88UYMHjw4PvnJT8aSJUti8ODBtX4oAKCL5DUA7Bk1L9133313re8SAKgxeQ0Ae0by93QDAABAT1XKsiwreoj3q1arUS6Xix6DnMaMGZNr/S233JJr/dKlS3OtL5VKudbn/bEaPXp0rvWnn356rvV5n39E/mPw29/+Ntf6KVOm5Fr/yCOP5FoPtVCpVOr6w8pk9t5BZstsmQ35dJTXznQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCKlLMuyood4v2q1GuVyuegxoEc788wzc9/H/fffX4NJuu7555/Ptf4zn/lMrvUbN27MtR4iIiqVSjQ2NhY9xk7JbCiezJbZFK+jvHamGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASKWVZlhU9xPtVq9Uol8tFjwE92n777Zf7PqZNm5Zr/U033ZRr/f77759r/cKFC3OtnzhxYq71EBFRqVSisbGx6DF2SmZD8WS2zKZ4HeW1M90AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkEifogcA6s/bb7+d+z5uvfXWXOs//vGP51p/4YUX5lr/kY98JNd6ANgTZLbMpv450w0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAifQpegCAHXn00Udzrb/wwgtzrR86dGiu9aeffnqu9Y888kiu9QCwp8hsmc2uOdMNAAAAiSjdAAAAkIjSDQAAAIko3QAAAJBIp0v3U089FWeccUa0tLREqVSKBx98sN3tWZbF9ddfH0OHDo1+/frFxIkT48UXX6zVvADAbpDXAFAfOl26N23aFOPGjYu5c+fu8PY5c+bELbfcEt/73vfimWeeiQ996EMxadKkeOedd3IPCwDsHnkNAPWh018ZNnny5Jg8efIOb8uyLG6++ea49tpr48wzz4yIiB/96EfR1NQUDz74YEyZMiXftADAbpHXAFAfavqe7pdeeinWr18fEydObLuuXC5Ha2trLF68eIdrNm/eHNVqtd0GAKTTlbyOkNkA0BU1Ld3r16+PiIimpqZ21zc1NbXd9kGzZ8+Ocrnctg0fPryWIwEAH9CVvI6Q2QDQFYV/evnMmTOjUqm0bWvWrCl6JABgB2Q2AHReTUt3c3NzRERs2LCh3fUbNmxou+2DGhoaorGxsd0GAKTTlbyOkNkA0BU1Ld2jRo2K5ubmmD9/ftt11Wo1nnnmmTjuuONq+VAAQBfJawDYczr96eVvvfVWrFy5su3ySy+9FC+88EIMHDgwRowYEVdeeWV8/etfj4985CMxatSouO6666KlpSXOOuusWs4NAOyCvAaA+tDp0v3cc8/FiSee2HZ5xowZERExderUuOOOO+IrX/lKbNq0Kb74xS/Gm2++GZ/85Cfjsccei3333bd2UwMAuySvAaA+dLp0T5gwIbIs2+ntpVIpvva1r8XXvva1XIMBAF0nrwGgPpSyXSVyAarVapTL5aLHALq51atX51o/bNiwXOtvu+22XOv/+I//ONd69g6VSqWuP6xMZgO1ILPp7jrK68K/MgwAAAD2Vko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJ9ih4AIIUsywpdf9RRR+VaDwA9hcxmb+dMNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAk0qfoAdg7HXvssbnWf+5zn8u1/t/+7d9yrd9///1zrW9ubs61fsGCBbnWL1myJNd6AHoOmS2zgbSc6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARPoUPQD1qaWlJdf6xx9/PNf6fv365VqfV6lUyrU+y7Jc63/729/mWv/rX/861/pNmzblWh8RcdZZZ+W+jzz233//Qh8fYE+R2TI7L5kNaTnTDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQSCnLsqzoId6vWq1GuVwueowe77nnnsu1/qijjqrRJMUolUq51tfZj1Wn5X3+EY5B3ud//PHH51q/ZMmSXOupD5VKJRobG4seY6dkdn2Q2TI7r55+DGQ2eXWU1850AwAAQCJKNwAAACSidAMAAEAinS7dTz31VJxxxhnR0tISpVIpHnzwwXa3X3TRRVEqldptp556aq3mBQB2g7wGgPrQ6dK9adOmGDduXMydO3en+5x66qmxbt26tu2uu+7KNSQA0DnyGgDqQ5/OLpg8eXJMnjx5l/s0NDREc3Nzl4cCAPKR1wBQH5K8p3vhwoUxZMiQOOyww+Kyyy6LN954Y6f7bt68OarVarsNAEivM3kdIbMBoCtqXrpPPfXU+NGPfhTz58+Pb33rW7Fo0aKYPHlybN26dYf7z549O8rlcts2fPjwWo8EAHxAZ/M6QmYDQFd0+s/LOzJlypS2f3/0ox+NsWPHxiGHHBILFy6Mk046abv9Z86cGTNmzGi7XK1WhTgAJNbZvI6Q2QDQFcm/Muzggw+OAw44IFauXLnD2xsaGqKxsbHdBgDsWR3ldYTMBoCuSF66X3nllXjjjTdi6NChqR8KAOgieQ0AaXT6z8vfeuutdr8Ff+mll+KFF16IgQMHxsCBA+PGG2+Mc889N5qbm2PVqlXxla98JQ499NCYNGlSTQcHAHZOXgNAfeh06X7uuefixBNPbLv8u/d2TZ06Nb773e/GsmXL4oc//GG8+eab0dLSEqecckrcdNNN0dDQULupAYBdktcAUB86XbonTJgQWZbt9PZ/+Zd/yTUQAJCfvAaA+pD8Pd0AAADQU9X8K8PYOxx11FG51u/q7EpP0NOff4RjkPf5X3jhhbnWr1+/Ptf6l19+Odd6YM+R2fn09Ocf4RjIbFJzphsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAEillWZYVPcT7VavVKJfLRY/R491www251l977bW1GaQgpVIp1/o6+7HqtLzPP8IxKPr5r1u3Ltf6O++8M9f6uXPn5lofEfHKK6/kvo/urlKpRGNjY9Fj7JTMrg8yu3v/9zovmd39XwMyu/vrKK+d6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACAREpZlmVFD/F+1Wo1yuVy0WP0eH369Mm1/tprry10fV6lUinX+jr7seq0vM8/wjHo6c//jTfeyD3DV7/61Vzr/+7v/i73DEWrVCrR2NhY9Bg7JbPrg8z23+u8evox6OnPX2bn11FeO9MNAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAImUsizLih7i/arVapTL5aLHoGDXX399rvVXXnllrvUDBgzItb7Ofqw6rVQq5b6Poo/BW2+9lWv9ggULcq1vbW3NtX7IkCG51ueV9zVQi//9X3zxxVzrJ06cmGv9K6+8kmt9LVQqlWhsbCx6jJ2S2UTI7KLJbJkts4vP7I7y2pluAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEiklGVZVvQQ71etVqNcLhc9Bt3c6NGjc60/8cQTc62fO3durvV5rVu3Ltf6lpaW3DPk/U/L/fffn2v9N77xjVzrf/GLX+RaP2jQoFzr874Gr7vuulzrx4wZk2t9PUTLt7/97Vzrr7nmmhpN0nWVSiUaGxuLHmOnZDa1ILNltsyW2d09szvKa2e6AQAAIBGlGwAAABJRugEAACCRTpXu2bNnx/jx46N///4xZMiQOOuss2LFihXt9nnnnXdi+vTpMWjQoNh///3j3HPPjQ0bNtR0aABg12Q2ANSHTpXuRYsWxfTp02PJkiXx+OOPx7vvvhunnHJKbNq0qW2fq666Kn7605/GvffeG4sWLYq1a9fGOeecU/PBAYCdk9kAUB/6dGbnxx57rN3lO+64I4YMGRJLly6NE044ISqVSvz93/99zJs3Lz7zmc9ERMTtt98eRxxxRCxZsiSOPfbY2k0OAOyUzAaA+pDrPd2VSiUiIgYOHBgREUuXLo133303Jk6c2LbP4YcfHiNGjIjFixfneSgAIAeZDQDF6NSZ7vfbtm1bXHnllXH88ce3fTfc+vXro2/fvjFgwIB2+zY1NcX69et3eD+bN2+OzZs3t12uVqtdHQkA2AGZDQDF6fKZ7unTp8fy5cvj7rvvzjXA7Nmzo1wut23Dhw/PdX8AQHsyGwCK06XSffnll8fDDz8cCxYsiGHDhrVd39zcHFu2bIk333yz3f4bNmyI5ubmHd7XzJkzo1KptG1r1qzpykgAwA7IbAAoVqdKd5Zlcfnll8cDDzwQTz75ZIwaNard7UcffXTss88+MX/+/LbrVqxYEatXr47jjjtuh/fZ0NAQjY2N7TYAIB+ZDQD1oVPv6Z4+fXrMmzcvHnrooejfv3/be77K5XL069cvyuVyTJs2LWbMmBEDBw6MxsbGuOKKK+K4447zKagAsAfJbACoD50q3d/97ncjImLChAntrr/99tvjoosuioiIb3/729GrV68499xzY/PmzTFp0qT4zne+U5NhAYDdI7MBoD50qnRnWdbhPvvuu2/MnTs35s6d2+WhAIB8ZDYA1Idc39MNAAAA7Fwp251fhe9B1Wo1yuVy0WPQwzU0NORaf9RRR9Vokq5Zt25drvVDhw7NPUPeY3DnnXfmWr9x48Zc67u7vP8b/uEf/mGu9Zdddlmu9RHR7pO2uyLvz0E9fB1WpVKp6w8rk9nUA5kts7s7md39M7ujvHamGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASKWVZlhU9xPtVq9Uol8tFjwFADoMGDcp9HxdffHGu9U1NTbnWX3PNNbnW10KlUonGxsaix9gpmQ3Q/cns/DrKa2e6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACCRUpZlWdFDvF+1Wo1yuVz0GABQuEqlEo2NjUWPsVMyGwA6zmtnugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABLpVOmePXt2jB8/Pvr37x9DhgyJs846K1asWNFunwkTJkSpVGq3felLX6rp0ADArslsAKgPnSrdixYtiunTp8eSJUvi8ccfj3fffTdOOeWU2LRpU7v9Lrnkkli3bl3bNmfOnJoODQDsmswGgPrQpzM7P/bYY+0u33HHHTFkyJBYunRpnHDCCW3X77ffftHc3FybCQGATpPZAFAfcr2nu1KpRETEwIED213/4x//OA444IAYM2ZMzJw5M95+++2d3sfmzZujWq222wCA2pLZAFCQrIu2bt2anX766dnxxx/f7vrvf//72WOPPZYtW7Ysu/POO7MDDzwwO/vss3d6P7Nmzcoiwmaz2Ww22we2SqXS1ZiW2TabzWaz7aGto7wuZVmWRRdcdtll8c///M/xs5/9LIYNG7bT/Z588sk46aSTYuXKlXHIIYdsd/vmzZtj8+bNbZer1WoMHz68KyMBwF6lUqlEY2Nj7vuR2QCQTkd53an3dP/O5ZdfHg8//HA89dRTuwzviIjW1taIiJ0GeENDQzQ0NHRlDACgAzIbAIrVqdKdZVlcccUV8cADD8TChQtj1KhRHa554YUXIiJi6NChXRoQAOg8mQ0A9aFTpXv69Okxb968eOihh6J///6xfv36iIgol8vRr1+/WLVqVcybNy9OO+20GDRoUCxbtiyuuuqqOOGEE2Ls2LFJngAAsD2ZDQB1ojMfxBI7eeP47bffnmVZlq1evTo74YQTsoEDB2YNDQ3ZoYceml1zzTWd+iCYSqVS+BvhbTabzWarhy3PB6nt7D5lts1ms9lstd2SfZBaKtVqNcrlctFjAEDhavVBaqnIbADoOK9zfU83AAAAsHNKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACRSd6U7y7KiRwCAulDvmVjv8wHAntBRHtZd6d64cWPRIwBAXaj3TKz3+QBgT+goD0tZnf2aetu2bbF27dro379/lEql7W6vVqsxfPjwWLNmTTQ2NhYwYffnGObj+OXnGObj+OXTHY5flmWxcePGaGlpiV696u73421kdlqOX36OYT6OXz6OX371fgx3N6/77MGZdkuvXr1i2LBhHe7X2NhYlwe+O3EM83H88nMM83H88qn341cul4seoUMye89w/PJzDPNx/PJx/PKr52O4O3ldv78+BwAAgG5O6QYAAIBEul3pbmhoiFmzZkVDQ0PRo3RbjmE+jl9+jmE+jl8+jt+e41jn4/jl5xjm4/jl4/jlt7ccw7r7IDUAAADYW3S7M90AAADQXSjdAAAAkIjSDQAAAIko3QAAAJBItyvdc+fOjYMOOij23XffaG1tjWeffbbokbqFG264IUqlUrvt8MMPL3qsuvbUU0/FGWecES0tLVEqleLBBx9sd3uWZXH99dfH0KFDo1+/fjFx4sR48cUXixm2DnV0/C666KLtXpOnnnpqMcPWodmzZ8f48eOjf//+MWTIkDjrrLNixYoV7fZ55513Yvr06TFo0KDYf//949xzz40NGzYUNHH92Z1jOGHChO1eh1/60pcKmnjvIq+7TmZ3jrzOT2bnI7Pz6Ql53a1K9z333BMzZsyIWbNmxfPPPx/jxo2LSZMmxauvvlr0aN3CkUceGevWrWvbfvaznxU9Ul3btGlTjBs3LubOnbvD2+fMmRO33HJLfO9734tnnnkmPvShD8WkSZPinXfe2cOT1qeOjl9ExKmnntruNXnXXXftwQnr26JFi2L69OmxZMmSePzxx+Pdd9+NU045JTZt2tS2z1VXXRU//elP4957741FixbF2rVr45xzzilw6vqyO8cwIuKSSy5p9zqcM2dOQRPvPeR1fjJ798nr/GR2PjI7nx6R11k3cswxx2TTp09vu7x169aspaUlmz17doFTdQ+zZs3Kxo0bV/QY3VZEZA888EDb5W3btmXNzc3ZX/7lX7Zd9+abb2YNDQ3ZXXfdVcCE9e2Dxy/Lsmzq1KnZmWeeWcg83dGrr76aRUS2aNGiLMvee73ts88+2b333tu2z3/8x39kEZEtXry4qDHr2gePYZZl2ac//ensT/7kT4obai8lr/OR2V0nr/OT2fnJ7Hz2xrzuNme6t2zZEkuXLo2JEye2XderV6+YOHFiLF68uMDJuo8XX3wxWlpa4uCDD44vfOELsXr16qJH6rZeeumlWL9+fbvXY7lcjtbWVq/HTli4cGEMGTIkDjvssLjsssvijTfeKHqkulWpVCIiYuDAgRERsXTp0nj33XfbvQYPP/zwGDFihNfgTnzwGP7Oj3/84zjggANizJgxMXPmzHj77beLGG+vIa9rQ2bXhryuHZm9+2R2PntjXvcpeoDd9frrr8fWrVujqamp3fVNTU3xn//5nwVN1X20trbGHXfcEYcddlisW7cubrzxxvjUpz4Vy5cvj/79+xc9Xrezfv36iIgdvh5/dxu7duqpp8Y555wTo0aNilWrVsVXv/rVmDx5cixevDh69+5d9Hh1Zdu2bXHllVfG8ccfH2PGjImI916Dffv2jQEDBrTb12twx3Z0DCMiPv/5z8fIkSOjpaUlli1bFn/6p38aK1asiPvvv7/Aabs3eZ2fzK4deV0bMnv3yex89ta87jalm3wmT57c9u+xY8dGa2trjBw5Mn7yk5/EtGnTCpyMnmrKlClt//7oRz8aY8eOjUMOOSQWLlwYJ510UoGT1Z/p06fH8uXLvaczh50dwy9+8Ytt//7oRz8aQ4cOjZNOOilWrVoVhxxyyJ4eEyJCZlN/ZPbuk9n57K153W3+vPyAAw6I3r17b/cpfxs2bIjm5uaCpuq+BgwYEKNHj46VK1cWPUq39LvXnNdj7Rx88MFxwAEHeE1+wOWXXx4PP/xwLFiwIIYNG9Z2fXNzc2zZsiXefPPNdvt7DW5vZ8dwR1pbWyMivA5zkNe1J7O7Tl6nIbN3TGbnszfndbcp3X379o2jjz465s+f33bdtm3bYv78+XHccccVOFn39NZbb8WqVati6NChRY/SLY0aNSqam5vbvR6r1Wo888wzXo9d9Morr8Qbb7zhNfl/siyLyy+/PB544IF48sknY9SoUe1uP/roo2OfffZp9xpcsWJFrF692mvw/3R0DHfkhRdeiIjwOsxBXteezO46eZ2GzG5PZufTE/K6W/15+YwZM2Lq1KnxiU98Io455pi4+eabY9OmTXHxxRcXPVrdu/rqq+OMM86IkSNHxtq1a2PWrFnRu3fvuOCCC4oerW699dZb7X579tJLL8ULL7wQAwcOjBEjRsSVV14ZX//61+MjH/lIjBo1Kq677rpoaWmJs846q7ih68iujt/AgQPjxhtvjHPPPTeam5tj1apV8ZWvfCUOPfTQmDRpUoFT14/p06fHvHnz4qGHHor+/fu3veerXC5Hv379olwux7Rp02LGjBkxcODAaGxsjCuuuCKOO+64OPbYYwuevj50dAxXrVoV8+bNi9NOOy0GDRoUy5Yti6uuuipOOOGEGDt2bMHTd2/yOh+Z3TnyOj+ZnY/MzqdH5HWxH57eebfeems2YsSIrG/fvtkxxxyTLVmypOiRuoXzzz8/Gzp0aNa3b9/swAMPzM4///xs5cqVRY9V1xYsWJBFxHbb1KlTsyx772tIrrvuuqypqSlraGjITjrppGzFihXFDl1HdnX83n777eyUU07JBg8enO2zzz7ZyJEjs0suuSRbv3590WPXjR0du4jIbr/99rZ9fvvb32Zf/vKXsw9/+MPZfvvtl5199tnZunXrihu6znR0DFevXp2dcMIJ2cCBA7OGhobs0EMPza655pqsUqkUO/heQl53nczuHHmdn8zOR2bn0xPyupRlWZamzgMAAEDP1m3e0w0AAADdjdINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJ/C+W5JsWqtPS4wAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ],
      "id": "-1v4al4l4b7x"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Training}$**"
      ],
      "metadata": {
        "id": "bMBJKJRA4b7y"
      },
      "id": "bMBJKJRA4b7y"
    },
    {
      "cell_type": "code",
      "source": [
        "# Initial guess\n",
        "cov_mat = []\n",
        "for i in range(len_s1):\n",
        "    xi = s1_matrix[i, :].reshape(-1, 1)\n",
        "    mat = 0\n",
        "    for j in range(len_s1):\n",
        "        xj = s1_matrix[j, :].reshape(-1, 1)\n",
        "        mat += (xj-xi) @ (xj-xi).T\n",
        "    cov_mat.append(mat/(len_s1-1))\n",
        "\n",
        "#Training\n",
        "cov_mat= Algorithm1(cov_mat)\n",
        "\n",
        "p = 0.08/2\n",
        "stepsize = 0.02\n",
        "C_max= critical_value(p,stepsize)"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:59:19.687328Z",
          "iopub.execute_input": "2024-05-07T11:59:19.688652Z",
          "iopub.status.idle": "2024-05-07T12:00:20.093467Z",
          "shell.execute_reply.started": "2024-05-07T11:59:19.688609Z",
          "shell.execute_reply": "2024-05-07T12:00:20.091898Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "fbf41039-db8b-4ff0-c7c8-a0141540df8b",
        "id": "NOJbSMJo4b7y"
      },
      "execution_count": 80,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "L_prev= -250\n",
            "L_new= -54.223803259768275\n",
            "t= 0\n",
            "--------------------------------------------------\n",
            "L_prev= -54.223803259768275\n",
            "L_new= -32.77853295616252\n",
            "t= 1\n",
            "--------------------------------------------------\n",
            "L_prev= -32.77853295616252\n",
            "L_new= -4.205711607345112\n",
            "t= 2\n",
            "--------------------------------------------------\n",
            "L_prev= -4.205711607345112\n",
            "L_new= 15.978974878133844\n",
            "t= 3\n",
            "--------------------------------------------------\n",
            "L_prev= 15.978974878133844\n",
            "L_new= 26.60367430738666\n",
            "t= 4\n",
            "--------------------------------------------------\n",
            "L_prev= 26.60367430738666\n",
            "L_new= 30.363463968687068\n",
            "t= 5\n",
            "--------------------------------------------------\n",
            "L_prev= 30.363463968687068\n",
            "L_new= 32.652175046006434\n",
            "t= 6\n",
            "--------------------------------------------------\n",
            "L_prev= 32.652175046006434\n",
            "L_new= 35.020940653487685\n",
            "t= 7\n",
            "--------------------------------------------------\n",
            "L_prev= 35.020940653487685\n",
            "L_new= 37.04688342424929\n",
            "t= 8\n",
            "--------------------------------------------------\n",
            "L_prev= 37.04688342424929\n",
            "L_new= 38.151664632589636\n",
            "t= 9\n",
            "--------------------------------------------------\n",
            "L_prev= 38.151664632589636\n",
            "L_new= 38.753004011115536\n",
            "t= 10\n",
            "--------------------------------------------------\n",
            "L_prev= 38.753004011115536\n",
            "L_new= 38.985286040124436\n",
            "t= 11\n",
            "--------------------------------------------------\n"
          ]
        }
      ],
      "id": "NOJbSMJo4b7y"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Change detection decision}$**"
      ],
      "metadata": {
        "id": "8Z5By3_q4b7z"
      },
      "id": "8Z5By3_q4b7z"
    },
    {
      "cell_type": "code",
      "source": [
        "Delta = delta()\n",
        "print('Delta:',Delta)\n",
        "\n",
        "if Delta < C_max:\n",
        "    print('THERE IS CHANGE. DISTRIBUTION OF S PRIME IS DIFFERENT THAN DISTRIBUTION OF S ')\n",
        "\n",
        "else:\n",
        "    print('THERE IS NO CHANGE. DISTRIBUTION OF S PRIME IS SAME AS DISTRIBUTION OF S ')"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T12:00:38.325910Z",
          "iopub.execute_input": "2024-05-07T12:00:38.326338Z",
          "iopub.status.idle": "2024-05-07T12:01:03.488229Z",
          "shell.execute_reply.started": "2024-05-07T12:00:38.326306Z",
          "shell.execute_reply": "2024-05-07T12:01:03.486998Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "1388b3e3-2511-4500-d355-c2fddbbf2ae1",
        "id": "n6e6JRAc4b7z"
      },
      "execution_count": 81,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Delta: 0.0\n",
            "THERE IS NO CHANGE. DISTRIBUTION OF S PRIME IS SAME AS DISTRIBUTION OF S \n"
          ]
        }
      ],
      "id": "n6e6JRAc4b7z"
    },
    {
      "cell_type": "markdown",
      "source": [
        "# $\\textbf{Case 3: Adding noise in the same image. (S' = noisy version of S)}$"
      ],
      "metadata": {
        "id": "5SrFr32Amvhh"
      },
      "id": "5SrFr32Amvhh"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Adding noise in the image}$**"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T08:56:26.545090Z",
          "iopub.execute_input": "2024-05-07T08:56:26.545484Z",
          "iopub.status.idle": "2024-05-07T08:56:26.552424Z",
          "shell.execute_reply.started": "2024-05-07T08:56:26.545456Z",
          "shell.execute_reply": "2024-05-07T08:56:26.550826Z"
        },
        "id": "9sAFb0Gsmvhi"
      },
      "id": "9sAFb0Gsmvhi"
    },
    {
      "cell_type": "code",
      "source": [
        "#adding noise\n",
        "\n",
        "image_index=53629    #We can put any number here\n",
        "\n",
        "s1_matrix=x_train[image_index]\n",
        "len_s1=s1_matrix.shape[0]\n",
        "s2_matrix=s1_matrix.copy()\n",
        "len_s2=s2_matrix.shape[0]\n",
        "\n",
        "sprime_matrix=x_train[image_index].copy()\n",
        "len_sprime=sprime_matrix.shape[0]\n",
        "\n",
        "\"\"\"Adding 3% noise, if it doesn't detect change, try slightly increasing noise percentage\"\"\"\n",
        "\n",
        "for i in range(int(0.03*28*28)):\n",
        "    sprime_matrix[np.random.randint(0,28), np.random.randint(0,28)]=np.random.uniform(0,1)\n"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:56:43.425370Z",
          "iopub.execute_input": "2024-05-07T11:56:43.425798Z",
          "iopub.status.idle": "2024-05-07T11:56:43.435244Z",
          "shell.execute_reply.started": "2024-05-07T11:56:43.425763Z",
          "shell.execute_reply": "2024-05-07T11:56:43.433953Z"
        },
        "trusted": true,
        "id": "3fokFPy-mvhi"
      },
      "execution_count": 82,
      "outputs": [],
      "id": "3fokFPy-mvhi"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Displaying images}$**"
      ],
      "metadata": {
        "id": "8WNoF3yOmvhi"
      },
      "id": "8WNoF3yOmvhi"
    },
    {
      "cell_type": "code",
      "source": [
        "#Visualizing images\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "s1_matrix_view=s1_matrix*255\n",
        "sprime_matrix_view=sprime_matrix*255\n",
        "\n",
        "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n",
        "axes[0].imshow(s1_matrix_view, cmap='gray')\n",
        "axes[0].set_title('S1')\n",
        "axes[0].axis('on')\n",
        "\n",
        "axes[1].imshow(sprime_matrix_view, cmap='gray')\n",
        "axes[1].set_title('S prime')\n",
        "axes[1].axis('on')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:56:49.361234Z",
          "iopub.execute_input": "2024-05-07T11:56:49.361627Z",
          "iopub.status.idle": "2024-05-07T11:56:49.908651Z",
          "shell.execute_reply.started": "2024-05-07T11:56:49.361597Z",
          "shell.execute_reply": "2024-05-07T11:56:49.907432Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 533
        },
        "id": "JTLY6kICmvhj",
        "outputId": "252951b7-a9a5-4944-c88b-4588f7b1013f"
      },
      "execution_count": 83,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1000x500 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAIECAYAAAAEiE3VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsT0lEQVR4nO3de5CV9XkH8GcFdgXdXQRhl5WrGtSIYOtltUkMBgQxoTEhHTUmA9bxlsWJUmtCqiFeRiJtvSTFGNtGk04w1sRLY61NQMA6glasIkmLQnUk4eZ1V9cIFt7+kWablcte3vPjPQufz8w7wznn/b3n4Z2jX7777jmnIsuyLAAAAICS26/oAQAAAGBvpXQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwBAD7dkyZKoqKiIJUuWFD0K8AFKN+yDnn/++fjc5z4XI0aMiP333z8OOeSQOO200+Lb3/522z4/+9nP4vzzz48xY8ZEr169YuTIkcUNDABlrjPZCuybKrIsy4oeAthznnjiiTj11FNj+PDhMX369Kivr49169bF8uXLY+3atbFmzZqIiJgxY0bcc8898Yd/+IfxyiuvRK9eveLll18udngAKEOdzdaUtm/fHlu3bo3KysrYbz/X1aCcKN2wj/nkJz8Z//7v/x4vvPBC9O/fv91jmzdvjsGDB0dExPr162PQoEHRp0+f+NSnPhWrVq1SugFgJzqbrSm89957ijaUOf91wj5m7dq1cfTRR+/wj4KIaPePgoaGhujTp88enAwAeqbOZuuujB8/PsaMGRMrVqyIP/qjP4q+ffvGqFGj4vbbb2+33+/et/2jH/0orrrqqjjkkEOiX79+0dLSstP3dP/uuCtXroyPf/zj0a9fvzj88MPjxz/+cURELF26NBobG6Nv375xxBFHxMKFC3eY7de//nX86Z/+adTV1UVVVVUcffTR8b3vfa9rJwj2cUo37GNGjBgRK1asiFWrVhU9CgDsFUqRrW+++WacccYZcdxxx8W8efNi6NChcckll+y04F533XXxz//8z3HFFVfEDTfcEJWVlbs97qc+9alobGyMefPmRVVVVZx99tlxzz33xNlnnx1nnHFGfPOb34zW1tb43Oc+F2+//Xbb2k2bNsVJJ50UCxcujJkzZ8att94ahx9+eJx//vlxyy23dPvvCvua3kUPAOxZV1xxRUyZMiWOPfbYOPHEE+NjH/tYTJgwIU499VRXtgGgG0qRrevXr4+//uu/jlmzZkVExEUXXRSNjY0xe/bs+OIXv9juOO+99148/fTT0bdv304dd8GCBXHOOedERMRpp50WRx55ZHz+85+PJ554IhobGyMi4qijjorJkyfHT37yk5gxY0ZERPzFX/xFbNu2LZ5//vkYOHBgRERcfPHFcc4558Q3vvGNuOiiizo1A+zrXOmGfcxpp50Wy5Ytiz/+4z+O5557LubNmxeTJ0+OQw45JP7pn/6p6PEAoMcpRbb27t07LrroorbblZWVcdFFF8XmzZtjxYoV7fadPn16p8vugQceGGeffXbb7SOOOCL69+8fRx11VFvhjoi2P//3f/93RERkWRY/+clPYurUqZFlWbz22mtt2+TJk6O5uTmeeeaZTs0A+zqlG/ZBJ5xwQtx3333x5ptvxlNPPRWzZ8+Ot99+Oz73uc/FL3/5y6LHA4AeJ2+2NjQ0xAEHHNDuvtGjR0dE7PBBpqNGjer0XEOHDo2Kiop299XW1sawYcN2uC/it7+OHhHx6quvxltvvRV33HFHDBo0qN123nnnRcRvPyQO6JhfL4d9WGVlZZxwwglxwgknxOjRo+O8886Le++9N+bMmVP0aADQI+2JbO3Kr3T36tWrS/f/7ouNtm/fHhERX/jCF2L69Ok73Xfs2LGdngP2ZUo3EBERxx9/fEREbNiwoeBJAGDv0JVsXb9+fbS2tra72v3CCy9ERMTIkSOTzLc7gwYNiurq6ti2bVtMnDhxjz8/7E38ejnsYxYvXtz2U+zf9/DDD0fEb9/rBQB0Ximy9X/+53/iu9/9btvtrVu3xne/+90YNGhQHHfccaUbtpN69eoV06ZNi5/85Cc7/VT2V199dY/PBD2VK92wj7n00kvj3Xffjc985jNx5JFHxtatW+OJJ56Ie+65J0aOHNn2Pq2VK1e2ffjLmjVrorm5Oa6//vqIiBg3blxMnTq1sL8DAJSTzmbr7jQ0NMSNN94YL7/8cowePTruueeeePbZZ+OOO+4o7NtFvvnNb8bixYujsbExLrjggvjwhz8cb7zxRjzzzDOxcOHCeOONNwqZC3oapRv2MX/1V38V9957bzz88MNxxx13xNatW2P48OHxpS99Ka666qro379/REQ888wzcfXVV7db+7vb06dPV7oB4P90Nlt356CDDorvf//7cemll8bf/u3fRl1dXfzN3/xNXHDBBen/ArtQV1cXTz31VFx77bVx3333xW233RYDBw6Mo48+Om688cbC5oKepiLb2e/CAAAAe8T48ePjtdde2+mvcQM9n/d0AwAAQCJKNwAAACSidAMAAEAi3tMNAAAAibjSDQAAAIko3QAAAJBI2X1P9/bt22P9+vVRXV0dFRUVRY8DAHtclmXx9ttvR0NDQ+y3X/n+fFxmA7Av62xel13pXr9+fQwbNqzoMQCgcOvWrYuhQ4cWPcYuyWwA6Divy650V1dXFz0CAJSFcs/Ecp+vMz7xiU/kWv/oo4+WaBIAeqqO8rDsSrdfTwOA3yr3TCz3+Tqjd++y+6cQAD1MR3lYvm8UAwAAgB4uWemeP39+jBw5Mvbff/9obGyMp556KtVTAQDdJK8BIK0kpfuee+6JWbNmxZw5c+KZZ56JcePGxeTJk2Pz5s0png4A6AZ5DQDpJSndN910U1xwwQVx3nnnxYc//OG4/fbbo1+/fvG9730vxdMBAN0grwEgvZKX7q1bt8aKFSti4sSJ//8k++0XEydOjGXLlu2w/5YtW6KlpaXdBgCk1dW8jpDZANAdJS/dr732Wmzbti3q6ura3V9XVxcbN27cYf+5c+dGbW1t2+b7PgEgva7mdYTMBoDuKPzTy2fPnh3Nzc1t27p164oeCQDYCZkNAF1X8i+nPPjgg6NXr16xadOmdvdv2rQp6uvrd9i/qqoqqqqqSj0GALAbXc3rCJkNAN1R8ivdlZWVcdxxx8WiRYva7tu+fXssWrQoTj755FI/HQDQDfIaAPaMkl/pjoiYNWtWTJ8+PY4//vg48cQT45ZbbonW1tY477zzUjwdANAN8hoA0ktSus8666x49dVX4+tf/3ps3Lgxjj322HjkkUd2+LAWAKA48hoA0qvIsiwreojf19LSErW1tUWPAQCFa25ujpqamqLH2KW9IbMnTZqUa/3PfvazEk0CQE/VUV4nudINANATKM0ApFb4V4YBAADA3krpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEikd9EDAADQfS0tLbnW19TUlGgS6JmefPLJXOuHDRuWa31DQ0Ou9ZQ/V7oBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIJHeRQ8AAED31dTUFPr8b7/9dq711dXVJZoEuqexsbHoEdjLudINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIn0LnoAAIB91XXXXZf7GFdffXUJJum+6urqQp8foNy50g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkEhFlmVZ0UP8vpaWlqitrS16DAAoXHNzc9TU1BQ9xi7JbADoOK9d6QYAAIBElG4AAABIROkGAACAREpeur/xjW9ERUVFu+3II48s9dMAADnIawDYM3qnOOjRRx8dCxcu/P8n6Z3kaQCAHOQ1AKSXJF179+4d9fX1KQ4NAJSIvAaA9JK8p/vFF1+MhoaGOPTQQ+Pcc8+NV155JcXTAAA5yGsASK/k39P9L//yL/HOO+/EEUccERs2bIhrrrkmfv3rX8eqVauiurp6h/23bNkSW7Zsabvd0tISw4YNK+VIANAjpfye7q7mdYTMBoCd6TCvs8TefPPNrKamJvu7v/u7nT4+Z86cLCJsNpvNZrN9YGtubk4d053Oa5lts9lsNtvOt47yOvlXhvXv3z9Gjx4da9as2enjs2fPjubm5rZt3bp1qUcCAD6go7yOkNkA0B3JS/c777wTa9eujSFDhuz08aqqqqipqWm3AQB7Vkd5HSGzAaA7Sl66r7jiili6dGm8/PLL8cQTT8RnPvOZ6NWrV5xzzjmlfioAoJvkNQDsGSX/yrBf/epXcc4558Trr78egwYNio9+9KOxfPnyGDRoUKmfCgDoJnkNAHtGyT+9PK+Wlpaora0tegwAKFzKTy8vBZkNAB3ndfL3dAMAAMC+quS/Xg4AQOfcdNNNuY8xa9asEkwCQCqudAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIhVZlmVFD/H7Wlpaora2tugxerSDDjoo9zEuv/zyXOunT5+ea/369etzrW9sbMy1nvwqKipyrS+z/zV12ebNm3Otv+uuu3KtX716da71P/jBD3Kt37ZtW671/FZzc3PU1NQUPcYuyez8ZLbMLgcyW2aTT0d57Uo3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACRSkWVZVvQQv6+lpSVqa2uLHqNHGzx4cO5jLFu2LNf6kSNH5p4BKM7o0aNzrV+7dm2JJtm3NTc3R01NTdFj7JLMzk9mA3nJ7OJ1lNeudAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIr2LHoDSq6mpyX2MkSNH5lrf3Nyca/0LL7yQa/3tt9+ea33e+aFol112WdEjAJ0gs2U2yOy9nyvdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJBI76IHoPR+/etf5z7GxIkTc61/4403cq1/7rnncq2Hnm7gwIG51l9yySUlmgRISWZDzyez6Ygr3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiXS5dD/22GMxderUaGhoiIqKinjggQfaPZ5lWXz961+PIUOGRN++fWPixInx4osvlmpeAKAT5DUAlIcul+7W1tYYN25czJ8/f6ePz5s3L771rW/F7bffHk8++WQccMABMXny5HjvvfdyDwsAdI68BoDy0OWvDJsyZUpMmTJlp49lWRa33HJLXHXVVfHpT386IiJ+8IMfRF1dXTzwwANx9tln55sWAOgUeQ0A5aGk7+l+6aWXYuPGje2+L7K2tjYaGxtj2bJlO12zZcuWaGlpabcBAOl0J68jZDYAdEdJS/fGjRsjIqKurq7d/XV1dW2PfdDcuXOjtra2bRs2bFgpRwIAPqA7eR0hswGgOwr/9PLZs2dHc3Nz27Zu3bqiRwIAdkJmA0DXlbR019fXR0TEpk2b2t2/adOmtsc+qKqqKmpqatptAEA63cnrCJkNAN1R0tI9atSoqK+vj0WLFrXd19LSEk8++WScfPLJpXwqAKCb5DUA7Dld/vTyd955J9asWdN2+6WXXopnn302BgwYEMOHD4/LLrssrr/++vjQhz4Uo0aNiquvvjoaGhrizDPPLOXcAMBuyGsAKA9dLt1PP/10nHrqqW23Z82aFRER06dPj7vuuiuuvPLKaG1tjQsvvDDeeuut+OhHPxqPPPJI7L///qWbGgDYLXkNAOWhy6V7/PjxkWXZLh+vqKiIa6+9Nq699tpcgwEA3SevAaA8dLl0U/5+85vf5D7G4sWLSzAJ7LsOOuigXOvvvvvuXOsnTJiQaz2wZ3zxi1/MfYw77rijBJNAz9XS0pJr/YgRI3Ktl9l0pPCvDAMAAIC9ldINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkEjvogcAKEfV1dW51v/4xz/OtX78+PG51gM9wx133FH0CNDjHXLIIbnWP/DAA7nWy2w64ko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACTSu+gBAFKoqKjItf7666/PtX78+PG51uf1i1/8Itf6N998s0STAMDuyWyZvbdzpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAEuld9AAAKdx888251s+cObNEk3TPc889l2v9lClTcq1/4403cq0HgM6S2TJ74sSJudYvXLiwRJOk4Uo3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAivYseAGBnzj333FzrZ86cWaJJivHVr3411/pNmzaVaBIA2D2Z3bMze+vWrbmPUVlZmWv9woULc89QzlzpBgAAgESUbgAAAEhE6QYAAIBEuly6H3vssZg6dWo0NDRERUVFPPDAA+0enzFjRlRUVLTbTj/99FLNCwB0grwGgPLQ5dLd2toa48aNi/nz5+9yn9NPPz02bNjQtt199925hgQAukZeA0B56PKnl0+ZMiWmTJmy232qqqqivr6+20MBAPnIawAoD0ne071kyZIYPHhwHHHEEXHJJZfE66+/vst9t2zZEi0tLe02ACC9ruR1hMwGgO4oeek+/fTT4wc/+EEsWrQobrzxxli6dGlMmTIltm3bttP9586dG7W1tW3bsGHDSj0SAPABXc3rCJkNAN3R5V8v78jZZ5/d9udjjjkmxo4dG4cddlgsWbIkJkyYsMP+s2fPjlmzZrXdbmlpEeIAkFhX8zpCZgNAdyT/yrBDDz00Dj744FizZs1OH6+qqoqampp2GwCwZ3WU1xEyGwC6I3np/tWvfhWvv/56DBkyJPVTAQDdJK8BII0u/3r5O++80+6n4C+99FI8++yzMWDAgBgwYEBcc801MW3atKivr4+1a9fGlVdeGYcffnhMnjy5pIMDALsmrwGgPHS5dD/99NNx6qmntt3+3Xu7pk+fHt/5zndi5cqV8f3vfz/eeuutaGhoiEmTJsV1110XVVVVpZsaANgteQ0A5aHLpXv8+PGRZdkuH//Xf/3XXAMBAPnJawAoD8nf0w0AAAD7qpJ/ZRhARERdXV2u9bfeemuu9RUVFbnW5/Vv//Zvudb/x3/8R4kmAYDdk9n7dmZXVlYWPcJez5VuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEikd9EDAOVnxowZuY8xe/bsXOsPOuigXOu3bduWa/2PfvSjXOv/7M/+LNf6V199Ndd6YN8xYMCAXOvfeOONEk1CEUqR2YsXL861/r777su1Pu/fIW9m19XV5Vovs+mIK90AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkEjvogcAdtS/f/9c65uamnKt/+pXv5prfUREv379ch8jjz/4gz/Itf4Xv/hFiSYBSOuNN97Itb6ysjLX+q1bt+Za39PJ7PzGjh2ba73Mpty50g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAifQuegBgR5/+9Kdzrb/22mtLNEn3/fKXv8y1/oYbbij0+QH2FVu3bi16hB5NZsts6Igr3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQSO+iB4ByVFFRkWv9zTffnGv9zJkzc63P69133819jJtuuinX+rvvvjv3DADs/WS2zIZy50o3AAAAJKJ0AwAAQCJKNwAAACTSpdI9d+7cOOGEE6K6ujoGDx4cZ555ZqxevbrdPu+99140NTXFwIED48ADD4xp06bFpk2bSjo0ALB7MhsAykOXSvfSpUujqakpli9fHj//+c/j/fffj0mTJkVra2vbPpdffnn89Kc/jXvvvTeWLl0a69evj89+9rMlHxwA2DWZDQDloUufXv7II4+0u33XXXfF4MGDY8WKFXHKKadEc3Nz/P3f/30sWLAgPvGJT0RExJ133hlHHXVULF++PE466aTSTQ4A7JLMBoDykOs93c3NzRERMWDAgIiIWLFiRbz//vsxceLEtn2OPPLIGD58eCxbtizPUwEAOchsAChGt7+ne/v27XHZZZfFRz7ykRgzZkxERGzcuDEqKyujf//+7fatq6uLjRs37vQ4W7ZsiS1btrTdbmlp6e5IAMBOyGwAKE63r3Q3NTXFqlWr4kc/+lGuAebOnRu1tbVt27Bhw3IdDwBoT2YDQHG6VbpnzpwZDz30UCxevDiGDh3adn99fX1s3bo13nrrrXb7b9q0Kerr63d6rNmzZ0dzc3Pbtm7duu6MBADshMwGgGJ1qXRnWRYzZ86M+++/Px599NEYNWpUu8ePO+646NOnTyxatKjtvtWrV8crr7wSJ5988k6PWVVVFTU1Ne02ACAfmQ0A5aFL7+luamqKBQsWxIMPPhjV1dVt7/mqra2Nvn37Rm1tbZx//vkxa9asGDBgQNTU1MSll14aJ598sk9BBYA9SGYDQHnoUun+zne+ExER48ePb3f/nXfeGTNmzIiIiJtvvjn222+/mDZtWmzZsiUmT54ct912W0mGBQA6R2YDQHnoUunOsqzDffbff/+YP39+zJ8/v9tDAQD5yGwAKA+5vqcbAAAA2LVuf0837M1uvvnmXOsvvfTSEk3SPb//Pbrdceqpp+ae4emnn859DIC93bvvvpv7GP369SvBJD2XzC4+s++9995c6//kT/4k13ryueaaa3IfY86cOSWYZO/lSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJFKRZVlW9BC/r6WlJWpra4segx5u0qRJudY/9NBDudb36tUr1/rXX3891/qLLroo1/r7778/13qgNJqbm6OmpqboMXZJZuf3/PPP5z7GMcccU4JJiiOzZTb0dB3ltSvdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJBI76IHgBS+9rWv5Vrfq1evEk3SPa2trbnWP/744yWaBICUjjnmmKJHKJzMltmQ1w033JBrfd7/D3XElW4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgER6Fz0A7MzkyZNzrW9sbCzRJMW4/fbbc61/9dVXSzQJAOyezJbZULSvfe1rRY+wW650AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAiSjcAAAAkonQDAABAIko3AAAAJKJ0AwAAQCJKNwAAACSidAMAAEAivYsegL3TKaeckmv9XXfdlWt9ZWVlrvXbt2/Ptf4LX/hCrvUPP/xwrvUA0Fkyu2dn9q233pr7GF/+8pdLMAmwK650AwAAQCJKNwAAACSidAMAAEAiXSrdc+fOjRNOOCGqq6tj8ODBceaZZ8bq1avb7TN+/PioqKhot1188cUlHRoA2D2ZDQDloUule+nSpdHU1BTLly+Pn//85/H+++/HpEmTorW1td1+F1xwQWzYsKFtmzdvXkmHBgB2T2YDQHno0qeXP/LII+1u33XXXTF48OBYsWJFu0++7NevX9TX15dmQgCgy2Q2AJSHXO/pbm5ujoiIAQMGtLv/hz/8YRx88MExZsyYmD17drz77ru7PMaWLVuipaWl3QYAlJbMBoBidPt7urdv3x6XXXZZfOQjH4kxY8a03f/5z38+RowYEQ0NDbFy5cr4yle+EqtXr4777rtvp8eZO3duXHPNNd0dAwDogMwGgOJ0u3Q3NTXFqlWr4vHHH293/4UXXtj252OOOSaGDBkSEyZMiLVr18Zhhx22w3Fmz54ds2bNarvd0tISw4YN6+5YAMAHyGwAKE63SvfMmTPjoYceisceeyyGDh26230bGxsjImLNmjU7DfCqqqqoqqrqzhgAQAdkNgAUq0ulO8uyuPTSS+P++++PJUuWxKhRozpc8+yzz0ZExJAhQ7o1IADQdTIbAMpDl0p3U1NTLFiwIB588MGorq6OjRs3RkREbW1t9O3bN9auXRsLFiyIM844IwYOHBgrV66Myy+/PE455ZQYO3Zskr8AALAjmQ0A5aFLpfs73/lORESMHz++3f133nlnzJgxIyorK2PhwoVxyy23RGtrawwbNiymTZsWV111VckGBgA6JrMBoDx0+dfLd2fYsGGxdOnSXAMBAPnJbAAoD7m+pxsAAADYtYqsox+F72EtLS1RW1tb9Bjk1KdPn1zrb7vttlzrjz322FzrP/nJT+Zav3nz5lzrASIimpubo6ampugxdklm7x1ktsyGfd0//MM/dGvdb37zm7jwwgs7zGtXugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgEaUbAAAAElG6AQAAIBGlGwAAABJRugEAACARpRsAAAASUboBAAAgkYosy7Kih/h9LS0tUVtbW/QYAFC45ubmqKmpKXqMXZLZANBxXrvSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJCI0g0AAACJKN0AAACQiNINAAAAiSjdAAAAkIjSDQAAAIko3QAAAJBI2ZXuLMuKHgEAykK5Z2K5zwcAe0JHeVh2pfvtt98uegQAKAvlnonlPh8A7Akd5WFFVmY/pt6+fXusX78+qquro6KiYofHW1paYtiwYbFu3bqoqakpYMKezznMx/nLzznMx/nLpyecvyzL4u23346GhobYb7+y+/l4G5mdlvOXn3OYj/OXj/OXX7mfw87mde89OFOn7LfffjF06NAO96upqSnLE9+TOIf5OH/5OYf5OH/5lPv5q62tLXqEDsnsPcP5y885zMf5y8f5y6+cz2Fn8rp8f3wOAAAAPZzSDQAAAIn0uNJdVVUVc+bMiaqqqqJH6bGcw3ycv/ycw3ycv3ycvz3Huc7H+cvPOczH+cvH+ctvbzmHZfdBagAAALC36HFXugEAAKCnULoBAAAgEaUbAAAAElG6AQAAIJEeV7rnz58fI0eOjP333z8aGxvjqaeeKnqkHuEb3/hGVFRUtNuOPPLIoscqa4899lhMnTo1GhoaoqKiIh544IF2j2dZFl//+tdjyJAh0bdv35g4cWK8+OKLxQxbhjo6fzNmzNjhNXn66acXM2wZmjt3bpxwwglRXV0dgwcPjjPPPDNWr17dbp/33nsvmpqaYuDAgXHggQfGtGnTYtOmTQVNXH46cw7Hjx+/w+vw4osvLmjivYu87j6Z3TXyOj+ZnY/MzmdfyOseVbrvueeemDVrVsyZMyeeeeaZGDduXEyePDk2b95c9Gg9wtFHHx0bNmxo2x5//PGiRyprra2tMW7cuJg/f/5OH583b15861vfittvvz2efPLJOOCAA2Ly5Mnx3nvv7eFJy1NH5y8i4vTTT2/3mrz77rv34ITlbenSpdHU1BTLly+Pn//85/H+++/HpEmTorW1tW2fyy+/PH7605/GvffeG0uXLo3169fHZz/72QKnLi+dOYcRERdccEG71+G8efMKmnjvIa/zk9mdJ6/zk9n5yOx89om8znqQE088MWtqamq7vW3btqyhoSGbO3dugVP1DHPmzMnGjRtX9Bg9VkRk999/f9vt7du3Z/X19dlf/uVftt331ltvZVVVVdndd99dwITl7YPnL8uybPr06dmnP/3pQubpiTZv3pxFRLZ06dIsy377euvTp0927733tu3zn//5n1lEZMuWLStqzLL2wXOYZVn28Y9/PPvyl79c3FB7KXmdj8zuPnmdn8zOT2bnszfmdY+50r1169ZYsWJFTJw4se2+/fbbLyZOnBjLli0rcLKe48UXX4yGhoY49NBD49xzz41XXnml6JF6rJdeeik2btzY7vVYW1sbjY2NXo9dsGTJkhg8eHAcccQRcckll8Trr79e9Ehlq7m5OSIiBgwYEBERK1asiPfff7/da/DII4+M4cOHew3uwgfP4e/88Ic/jIMPPjjGjBkTs2fPjnfffbeI8fYa8ro0ZHZpyOvSkdmdJ7Pz2RvzunfRA3TWa6+9Ftu2bYu6urp299fV1cV//dd/FTRVz9HY2Bh33XVXHHHEEbFhw4a45ppr4mMf+1isWrUqqqurix6vx9m4cWNExE5fj797jN07/fTT47Of/WyMGjUq1q5dG1/72tdiypQpsWzZsujVq1fR45WV7du3x2WXXRYf+chHYsyYMRHx29dgZWVl9O/fv92+XoM7t7NzGBHx+c9/PkaMGBENDQ2xcuXK+MpXvhKrV6+O++67r8BpezZ5nZ/MLh15XRoyu/Nkdj57a173mNJNPlOmTGn789ixY6OxsTFGjBgR//iP/xjnn39+gZOxrzr77LPb/nzMMcfE2LFj47DDDoslS5bEhAkTCpys/DQ1NcWqVau8pzOHXZ3DCy+8sO3PxxxzTAwZMiQmTJgQa9eujcMOO2xPjwkRIbMpPzK782R2PntrXveYXy8/+OCDo1evXjt8yt+mTZuivr6+oKl6rv79+8fo0aNjzZo1RY/SI/3uNef1WDqHHnpoHHzwwV6THzBz5sx46KGHYvHixTF06NC2++vr62Pr1q3x1ltvtdvfa3BHuzqHO9PY2BgR4XWYg7wuPZndffI6DZm9czI7n705r3tM6a6srIzjjjsuFi1a1Hbf9u3bY9GiRXHyyScXOFnP9M4778TatWtjyJAhRY/SI40aNSrq6+vbvR5bWlriySef9Hrspl/96lfx+uuve03+nyzLYubMmXH//ffHo48+GqNGjWr3+HHHHRd9+vRp9xpcvXp1vPLKK16D/6ejc7gzzz77bESE12EO8rr0ZHb3yes0ZHZ7MjuffSGve9Svl8+aNSumT58exx9/fJx44olxyy23RGtra5x33nlFj1b2rrjiipg6dWqMGDEi1q9fH3PmzIlevXrFOeecU/RoZeudd95p99Ozl156KZ599tkYMGBADB8+PC677LK4/vrr40Mf+lCMGjUqrr766mhoaIgzzzyzuKHLyO7O34ABA+Kaa66JadOmRX19faxduzauvPLKOPzww2Py5MkFTl0+mpqaYsGCBfHggw9GdXV123u+amtro2/fvlFbWxvnn39+zJo1KwYMGBA1NTVx6aWXxsknnxwnnXRSwdOXh47O4dq1a2PBggVxxhlnxMCBA2PlypVx+eWXxymnnBJjx44tePqeTV7nI7O7Rl7nJ7Pzkdn57BN5XeyHp3fdt7/97Wz48OFZZWVlduKJJ2bLly8veqQe4ayzzsqGDBmSVVZWZoccckh21llnZWvWrCl6rLK2ePHiLCJ22KZPn55l2W+/huTqq6/O6urqsqqqqmzChAnZ6tWrix26jOzu/L377rvZpEmTskGDBmV9+vTJRowYkV1wwQXZxo0bix67bOzs3EVEduedd7bt85vf/Cb70pe+lB100EFZv379ss985jPZhg0bihu6zHR0Dl955ZXslFNOyQYMGJBVVVVlhx9+ePbnf/7nWXNzc7GD7yXkdffJ7K6R1/nJ7Hxkdj77Ql5XZFmWpanzAAAAsG/rMe/pBgAAgJ5G6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgESUbgAAAEhE6QYAAIBElG4AAABIROkGAACARJRuAAAASETpBgAAgET+F5rNXg3NTM9fAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ],
      "id": "JTLY6kICmvhj"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Training for selected S1}$**"
      ],
      "metadata": {
        "id": "wcdI1729mvhj"
      },
      "id": "wcdI1729mvhj"
    },
    {
      "cell_type": "code",
      "source": [
        "#Initial guess\n",
        "cov_mat_initial = []\n",
        "\n",
        "for i in range(len_s1):\n",
        "    xi = s1_matrix[i, :].reshape(-1, 1)\n",
        "    mat = 0\n",
        "    for j in range(len_s1):\n",
        "        xj = s1_matrix[j, :].reshape(-1, 1)\n",
        "        mat += (xj-xi) @ (xj-xi).T\n",
        "\n",
        "    cov_mat_initial.append(mat/(len_s1-1))\n",
        "\n",
        "#Training\n",
        "cov_mat= Algorithm1(cov_mat_initial)\n",
        "\n",
        "p = 0.08\n",
        "stepsize = 0.02\n",
        "C_max= critical_value(p,stepsize)"
      ],
      "metadata": {
        "scrolled": true,
        "execution": {
          "iopub.status.busy": "2024-05-07T11:56:56.851756Z",
          "iopub.execute_input": "2024-05-07T11:56:56.852206Z",
          "iopub.status.idle": "2024-05-07T11:58:05.036199Z",
          "shell.execute_reply.started": "2024-05-07T11:56:56.852172Z",
          "shell.execute_reply": "2024-05-07T11:58:05.034590Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "glxRNfoYmvhj",
        "outputId": "cb3dbc7e-0ec8-4d28-f739-03aaffa6898e"
      },
      "execution_count": 84,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "L_prev= -250\n",
            "L_new= -43.842275419465764\n",
            "t= 0\n",
            "--------------------------------------------------\n",
            "L_prev= -43.842275419465764\n",
            "L_new= -20.14056892812331\n",
            "t= 1\n",
            "--------------------------------------------------\n",
            "L_prev= -20.14056892812331\n",
            "L_new= 7.962562871482346\n",
            "t= 2\n",
            "--------------------------------------------------\n",
            "L_prev= 7.962562871482346\n",
            "L_new= 26.080088303206274\n",
            "t= 3\n",
            "--------------------------------------------------\n",
            "L_prev= 26.080088303206274\n",
            "L_new= 31.95055718815513\n",
            "t= 4\n",
            "--------------------------------------------------\n",
            "L_prev= 31.95055718815513\n",
            "L_new= 33.85029680034366\n",
            "t= 5\n",
            "--------------------------------------------------\n",
            "L_prev= 33.85029680034366\n",
            "L_new= 34.76599292403555\n",
            "t= 6\n",
            "--------------------------------------------------\n",
            "L_prev= 34.76599292403555\n",
            "L_new= 35.2330975465774\n",
            "t= 7\n",
            "--------------------------------------------------\n",
            "L_prev= 35.2330975465774\n",
            "L_new= 35.529214655706554\n",
            "t= 8\n",
            "--------------------------------------------------\n"
          ]
        }
      ],
      "id": "glxRNfoYmvhj"
    },
    {
      "cell_type": "markdown",
      "source": [
        "**$\\textbf{Change detection decision}$**"
      ],
      "metadata": {
        "id": "cjgohxF9mvhk"
      },
      "id": "cjgohxF9mvhk"
    },
    {
      "cell_type": "code",
      "source": [
        "#Change detection test\n",
        "Delta = delta()\n",
        "print('delta:',Delta)\n",
        "\n",
        "if Delta < C_max:\n",
        "    print('THERE IS CHANGE. DISTRIBUTION OF S PRIME IS DIFFERENT THAN DISTRIBUTION OF S ')\n",
        "\n",
        "else:\n",
        "    print('THERE IS NO CHANGE. DISTRIBUTION OF S PRIME IS SAME AS DISTRIBUTION OF S ')"
      ],
      "metadata": {
        "execution": {
          "iopub.status.busy": "2024-05-07T11:58:15.144286Z",
          "iopub.execute_input": "2024-05-07T11:58:15.144702Z",
          "iopub.status.idle": "2024-05-07T11:58:40.351546Z",
          "shell.execute_reply.started": "2024-05-07T11:58:15.144669Z",
          "shell.execute_reply": "2024-05-07T11:58:40.349862Z"
        },
        "trusted": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "cFjpD81tmvhk",
        "outputId": "b8762636-b87f-47ff-eda4-b5dac5937178"
      },
      "execution_count": 85,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "delta: -38.79296070504994\n",
            "THERE IS CHANGE. DISTRIBUTION OF S PRIME IS DIFFERENT THAN DISTRIBUTION OF S \n"
          ]
        }
      ],
      "id": "cFjpD81tmvhk"
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "dBPoc86lDRhb"
      },
      "id": "dBPoc86lDRhb",
      "execution_count": null,
      "outputs": []
    }
  ]
}