2.3_GroupLASSOMultinomialRegression.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Part 0: Import Libraries."
      ],
      "metadata": {
        "id": "nw0iwePutBgc"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "t8AsmGprplSG"
      },
      "outputs": [],
      "source": [
        "!pip install group-lasso --quiet"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.datasets import fetch_openml\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.patches as mpatches\n",
        "import seaborn as sns\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "from sklearn.linear_model import LogisticRegression, LogisticRegressionCV\n",
        "from group_lasso import GroupLasso\n",
        "from sklearn.metrics import accuracy_score\n",
        "from group_lasso import LogisticGroupLasso\n",
        "from sklearn.metrics import accuracy_score, classification_report"
      ],
      "metadata": {
        "id": "g7--ZF34p1kQ"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Part 1: Import & Process Data."
      ],
      "metadata": {
        "id": "H5umRaCvtFL8"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Run necessary code from notebook 1."
      ],
      "metadata": {
        "id": "HEUsm7XZg0kJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Load Dataset\n",
        "mnist = fetch_openml(\"mnist_784\")\n",
        "\n",
        "# Filter and Prepare Dataset\n",
        "keys = list(mnist.keys())\n",
        "df = pd.concat([mnist[keys[0]],mnist[keys[1]]],axis=1)\n",
        "df['class'] = df['class'].apply(lambda x: int(x))\n",
        "bool_musk = df['class'].isin([3,5,8])\n",
        "df = df.loc[bool_musk].reset_index(drop=True)\n",
        "\n",
        "# Split data into training and testing\n",
        "cutoff = int(df.shape[0]*(0.8))\n",
        "X, y = df.iloc[:,:-1], df['class']\n",
        "X_train, X_test = X[:cutoff], X[cutoff:]\n",
        "y_train, y_test = y[:cutoff], y[cutoff:]\n",
        "\n",
        "# Randomize the split\n",
        "n = X_train.shape[0]\n",
        "rand = np.random.permutation(n)\n",
        "X_train = X_train.iloc[list(rand),:].reset_index(drop = True)\n",
        "y_train = y_train[rand].reset_index(drop = True)\n",
        "\n",
        "assert X_train.shape[0] + X_test.shape[0] == df.shape[0], 'ERROR'\n",
        "assert X_train.shape[1] == X_test.shape[1] == df.shape[1] - 1, 'ERROR' #class column == y\n",
        "assert y_train.shape[0] + y_test.shape[0] == df.shape[0], 'ERROR'\n",
        "\n",
        "# Formats\n",
        "print(f'X_train: {X_train.shape},\\nX_test: {X_test.shape},\\ny_train: {y_train.shape},\\ny_test: {y_test.shape}')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "cvnnJwGsqSp0",
        "outputId": "7e767706-58f3-4aa4-e6aa-0cc4817d42a3"
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "X_train: (16223, 784),\n",
            "X_test: (4056, 784),\n",
            "y_train: (16223,),\n",
            "y_test: (4056,)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "df = pd.concat([mnist[keys[0]],mnist[keys[1]]],axis=1)\n",
        "i = 7 #select any row; i = 0 -> row 0\n",
        "obs = df.iloc[i,:-1].to_numpy()\n",
        "obs_image = obs.astype(np.uint8).reshape(28, 28)\n",
        "plt.imshow(obs_image, plt.cm.binary)\n",
        "plt.show()\n",
        "print(f'\\n*** Image is classified as a {df.loc[i,\"class\"]} ***')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 465
        },
        "id": "isovz7ZHGTGA",
        "outputId": "8b5eea88-d313-4968-ed55-2a055d89e16f"
      },
      "execution_count": 31,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGdCAYAAABU0qcqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAG39JREFUeJzt3X9s1PUdx/HX8aMHaHus1vbaUVgBlU2gbgy6RmU4OkpJlApZQF0CxkDEYoboNF1UdFtWxcQxDcMsUVAnoCQCwTEcFlt0FhaqBMm2jjZ1VKFl4rgrRQqhn/1BvHnQAt/jru9e+3wk34Te3af39rvv7smXu37rc845AQDQzfpZDwAA6JsIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMDHAeoBzdXR06NChQ0pNTZXP57MeBwDgkXNOra2tysnJUb9+XZ/n9LgAHTp0SLm5udZjAAAuU1NTk4YNG9bl/T0uQKmpqZLODp6WlmY8DQDAq3A4rNzc3MjreVcSFqCVK1fqmWeeUXNzs/Lz8/X8889r0qRJF1331T+7paWlESAASGIXexslIR9CeP3117V06VItW7ZMH374ofLz81VcXKwjR44k4ukAAEkoIQF69tlntWDBAt199936zne+oxdeeEFDhgzRSy+9lIinAwAkobgH6NSpU6qtrVVRUdH/n6RfPxUVFammpua8x7e3tyscDkdtAIDeL+4B+vzzz3XmzBllZWVF3Z6VlaXm5ubzHl9RUaFAIBDZ+AQcAPQN5j+IWl5erlAoFNmampqsRwIAdIO4fwouIyND/fv3V0tLS9TtLS0tCgaD5z3e7/fL7/fHewwAQA8X9zOglJQUTZgwQZWVlZHbOjo6VFlZqcLCwng/HQAgSSXk54CWLl2qefPm6fvf/74mTZqkFStWqK2tTXfffXcing4AkIQSEqA5c+boP//5jx5//HE1Nzfrhhtu0LZt2877YAIAoO/yOeec9RBfFw6HFQgEFAqFuBICACShS30dN/8UHACgbyJAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMDLAeAEiEf/3rXzGtO3XqlOc17733nuc19913n+c1Pp/P85reqLS01POa9evXx/RcKSkpMa3DpeEMCABgggABAEzEPUBPPPGEfD5f1DZmzJh4Pw0AIMkl5D2g66+/Xu+8887/n2QAbzUBAKIlpAwDBgxQMBhMxLcGAPQSCXkP6MCBA8rJydHIkSN111136eDBg10+tr29XeFwOGoDAPR+cQ9QQUGB1qxZo23btmnVqlVqbGzUzTffrNbW1k4fX1FRoUAgENlyc3PjPRIAoAeKe4BKSkr0k5/8ROPHj1dxcbG2bt2qY8eO6Y033uj08eXl5QqFQpGtqakp3iMBAHqghH86YOjQobr22mtVX1/f6f1+v19+vz/RYwAAepiE/xzQ8ePH1dDQoOzs7EQ/FQAgicQ9QA899JCqq6v1ySef6IMPPtDtt9+u/v3764477oj3UwEAkljc/wnu008/1R133KGjR4/q6quv1k033aRdu3bp6quvjvdTAQCSmM8556yH+LpwOKxAIKBQKKS0tDTrcRBn+/fv97zm5Zdf9rxmw4YNntdIUkdHh+c1n332mec1sfzfjouRxm7evHkxrVuxYoXnNbxuXfrrONeCAwCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMcDFSdKvbbrvN85o//elPCZjEFhcjTQ7V1dWe19x0000JmCS5cDFSAECPRoAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMDrAdA3/LjH//Y85ruvBp2Zmam5zX33HOP5zUdHR2e1/Tr131/X/zggw88r4nlytHo2zgDAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMcDFSdKtFixZ5XlNaWhr/QbowcOBAz2uCwWACJrEVDoc9rxk7dqznNZ999pnnNbGI9RiaOHFifAdBFM6AAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATXIwU3WrAAO+HXG5ubgImwYW8/fbbntf897//TcAk8RHrMeT3++M8Cb6OMyAAgAkCBAAw4TlAO3fu1K233qqcnBz5fD5t2rQp6n7nnB5//HFlZ2dr8ODBKioq0oEDB+I1LwCgl/AcoLa2NuXn52vlypWd3r98+XI999xzeuGFF7R7925dccUVKi4u1smTJy97WABA7+H5HeGSkhKVlJR0ep9zTitWrNCjjz6qmTNnSpJeeeUVZWVladOmTZo7d+7lTQsA6DXi+h5QY2OjmpubVVRUFLktEAiooKBANTU1na5pb29XOByO2gAAvV9cA9Tc3CxJysrKiro9Kysrct+5KioqFAgEIhsfuQWAvsH8U3Dl5eUKhUKRrampyXokAEA3iGuAgsGgJKmlpSXq9paWlsh95/L7/UpLS4vaAAC9X1wDlJeXp2AwqMrKysht4XBYu3fvVmFhYTyfCgCQ5Dx/Cu748eOqr6+PfN3Y2Ki9e/cqPT1dw4cP15IlS/TrX/9a11xzjfLy8vTYY48pJydHpaWl8ZwbAJDkPAdoz549uuWWWyJfL126VJI0b948rVmzRg8//LDa2tq0cOFCHTt2TDfddJO2bdumQYMGxW9qAEDS8znnnPUQXxcOhxUIBBQKhXg/CLhM69evj2ndH/7wB89rqqurY3qu7hDrhVJ5DYrNpb6Om38KDgDQNxEgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMCE51/HAODy/fGPf/S85qmnnvK8pqGhwfMaSTp16lRM67rDDTfc4HnNwIED4z8ILhtnQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACS5Gim71ySefeF7z6quvel7zzjvveF7Tnd577z3Pa3w+XwImiZ+0tDTPa55++mnPa2bMmOF5zeDBgz2vQeJxBgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmOBipIjZxx9/7HnNbbfd5nnNwYMHPa9B95s8ebLnNQsXLkzAJEgWnAEBAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACa4GCl6POec9Qhx1xv/m7Zs2eJ5zdatWz2vmTFjhuc16Jk4AwIAmCBAAAATngO0c+dO3XrrrcrJyZHP59OmTZui7p8/f758Pl/UNn369HjNCwDoJTwHqK2tTfn5+Vq5cmWXj5k+fboOHz4c2datW3dZQwIAeh/PH0IoKSlRSUnJBR/j9/sVDAZjHgoA0Psl5D2gqqoqZWZm6rrrrtOiRYt09OjRLh/b3t6ucDgctQEAer+4B2j69Ol65ZVXVFlZqaefflrV1dUqKSnRmTNnOn18RUWFAoFAZMvNzY33SACAHijuPwc0d+7cyJ/HjRun8ePHa9SoUaqqqtLUqVPPe3x5ebmWLl0a+TocDhMhAOgDEv4x7JEjRyojI0P19fWd3u/3+5WWlha1AQB6v4QH6NNPP9XRo0eVnZ2d6KcCACQRz/8Ed/z48aizmcbGRu3du1fp6elKT0/Xk08+qdmzZysYDKqhoUEPP/ywRo8ereLi4rgODgBIbp4DtGfPHt1yyy2Rr796/2bevHlatWqV9u3bp5dfflnHjh1TTk6Opk2bpl/96lfy+/3xmxoAkPQ8B2jKlCkXvJDi22+/fVkDIXmMGzfO85qqqirPa1599VXPa2K9+sagQYNiWtdTvfjiizGte+655+I8CXA+rgUHADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEz53oUtbGwiHwwoEAgqFQvx2VOAyhUKhmNalp6fHeZLObdmyxfOaGTNmJGASxNOlvo5zBgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmBhgPQCAxHn77betRwC6xBkQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCi5H2MqdPn/a8JtYLVk6dOtXzmsGDB8f0XJBeeuklz2uWLFkS/0GAOOEMCABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwwcVIe7D33nvP85rf/OY3ntf85S9/8bxGkj755BPPa3Jzc2N6rp7siy++8Lxm69atntc8+OCDnte0tbV5XhOrIUOGeF7DxWn7Ns6AAAAmCBAAwISnAFVUVGjixIlKTU1VZmamSktLVVdXF/WYkydPqqysTFdddZWuvPJKzZ49Wy0tLXEdGgCQ/DwFqLq6WmVlZdq1a5e2b9+u06dPa9q0aVH/zvzAAw9oy5Yt2rBhg6qrq3Xo0CHNmjUr7oMDAJKbpw8hbNu2LerrNWvWKDMzU7W1tZo8ebJCoZBefPFFrV27Vj/60Y8kSatXr9a3v/1t7dq1Sz/4wQ/iNzkAIKld1ntAoVBIkpSeni5Jqq2t1enTp1VUVBR5zJgxYzR8+HDV1NR0+j3a29sVDoejNgBA7xdzgDo6OrRkyRLdeOONGjt2rCSpublZKSkpGjp0aNRjs7Ky1Nzc3On3qaioUCAQiGy98WO6AIDzxRygsrIy7d+/X+vXr7+sAcrLyxUKhSJbU1PTZX0/AEByiOkHURcvXqy33npLO3fu1LBhwyK3B4NBnTp1SseOHYs6C2ppaVEwGOz0e/n9fvn9/ljGAAAkMU9nQM45LV68WBs3btSOHTuUl5cXdf+ECRM0cOBAVVZWRm6rq6vTwYMHVVhYGJ+JAQC9gqczoLKyMq1du1abN29Wampq5H2dQCCgwYMHKxAI6J577tHSpUuVnp6utLQ03X///SosLOQTcACAKJ4CtGrVKknSlClTom5fvXq15s+fL0n67W9/q379+mn27Nlqb29XcXGxfv/738dlWABA7+FzzjnrIb4uHA4rEAgoFAopLS3NehxTN9xwg+c1H3/8cfwH6cJ9993neU1qamoCJrG1fft2z2tqa2s9r/H5fJ7XxOrcv2ReiliOh9mzZ3teg57vUl/HuRYcAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATMT0G1EBSfyajW6WmZnpec1tt90W03P97ne/87xm0KBBMT0X+i7OgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAE1yMtAdbvXq15zXPP/+85zUvv/yy5zW91ejRoz2vGTJkiOc1N998s+c1CxYs8Lxm3LhxntcA3YUzIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABBcj7cG++93vel6zatUqz2sKCgo8r5GkRx991POaL774wvOa0tJSz2umTZvmeY0kzZw50/OaYDAY03MBfR1nQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACZ9zzlkP8XXhcFiBQEChUEhpaWnW4wAAPLrU13HOgAAAJggQAMCEpwBVVFRo4sSJSk1NVWZmpkpLS1VXVxf1mClTpsjn80Vt9957b1yHBgAkP08Bqq6uVllZmXbt2qXt27fr9OnTmjZtmtra2qIet2DBAh0+fDiyLV++PK5DAwCSn6ffiLpt27aor9esWaPMzEzV1tZq8uTJkduHDBnCb4kEAFzQZb0HFAqFJEnp6elRt7/22mvKyMjQ2LFjVV5erhMnTnT5Pdrb2xUOh6M2AEDv5+kM6Os6Ojq0ZMkS3XjjjRo7dmzk9jvvvFMjRoxQTk6O9u3bp0ceeUR1dXV68803O/0+FRUVevLJJ2MdAwCQpGL+OaBFixbpz3/+s95//30NGzasy8ft2LFDU6dOVX19vUaNGnXe/e3t7Wpvb498HQ6HlZuby88BAUCSutSfA4rpDGjx4sV66623tHPnzgvGR5IKCgokqcsA+f1++f3+WMYAACQxTwFyzun+++/Xxo0bVVVVpby8vIuu2bt3ryQpOzs7pgEBAL2TpwCVlZVp7dq12rx5s1JTU9Xc3CxJCgQCGjx4sBoaGrR27VrNmDFDV111lfbt26cHHnhAkydP1vjx4xPyHwAASE6e3gPy+Xyd3r569WrNnz9fTU1N+ulPf6r9+/erra1Nubm5uv322/Xoo49e8vs5XAsOAJJbQt4DulircnNzVV1d7eVbAgD6KK4FBwAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwMcB6gHM55yRJ4XDYeBIAQCy+ev3+6vW8Kz0uQK2trZKk3Nxc40kAAJejtbVVgUCgy/t97mKJ6mYdHR06dOiQUlNT5fP5ou4Lh8PKzc1VU1OT0tLSjCa0x344i/1wFvvhLPbDWT1hPzjn1NraqpycHPXr1/U7PT3uDKhfv34aNmzYBR+TlpbWpw+wr7AfzmI/nMV+OIv9cJb1frjQmc9X+BACAMAEAQIAmEiqAPn9fi1btkx+v996FFPsh7PYD2exH85iP5yVTPuhx30IAQDQNyTVGRAAoPcgQAAAEwQIAGCCAAEATCRNgFauXKlvfetbGjRokAoKCvS3v/3NeqRu98QTT8jn80VtY8aMsR4r4Xbu3Klbb71VOTk58vl82rRpU9T9zjk9/vjjys7O1uDBg1VUVKQDBw7YDJtAF9sP8+fPP+/4mD59us2wCVJRUaGJEycqNTVVmZmZKi0tVV1dXdRjTp48qbKyMl111VW68sorNXv2bLW0tBhNnBiXsh+mTJly3vFw7733Gk3cuaQI0Ouvv66lS5dq2bJl+vDDD5Wfn6/i4mIdOXLEerRud/311+vw4cOR7f3337ceKeHa2tqUn5+vlStXdnr/8uXL9dxzz+mFF17Q7t27dcUVV6i4uFgnT57s5kkT62L7QZKmT58edXysW7euGydMvOrqapWVlWnXrl3avn27Tp8+rWnTpqmtrS3ymAceeEBbtmzRhg0bVF1drUOHDmnWrFmGU8ffpewHSVqwYEHU8bB8+XKjibvgksCkSZNcWVlZ5OszZ864nJwcV1FRYThV91u2bJnLz8+3HsOUJLdx48bI1x0dHS4YDLpnnnkmctuxY8ec3+9369atM5iwe5y7H5xzbt68eW7mzJkm81g5cuSIk+Sqq6udc2f/tx84cKDbsGFD5DH/+Mc/nCRXU1NjNWbCnbsfnHPuhz/8ofvZz35mN9Ql6PFnQKdOnVJtba2Kiooit/Xr109FRUWqqakxnMzGgQMHlJOTo5EjR+quu+7SwYMHrUcy1djYqObm5qjjIxAIqKCgoE8eH1VVVcrMzNR1112nRYsW6ejRo9YjJVQoFJIkpaenS5Jqa2t1+vTpqONhzJgxGj58eK8+Hs7dD1957bXXlJGRobFjx6q8vFwnTpywGK9LPe5ipOf6/PPPdebMGWVlZUXdnpWVpX/+859GU9koKCjQmjVrdN111+nw4cN68skndfPNN2v//v1KTU21Hs9Ec3OzJHV6fHx1X18xffp0zZo1S3l5eWpoaNAvfvELlZSUqKamRv3797ceL+46Ojq0ZMkS3XjjjRo7dqyks8dDSkqKhg4dGvXY3nw8dLYfJOnOO+/UiBEjlJOTo3379umRRx5RXV2d3nzzTcNpo/X4AOH/SkpKIn8eP368CgoKNGLECL3xxhu65557DCdDTzB37tzIn8eNG6fx48dr1KhRqqqq0tSpUw0nS4yysjLt37+/T7wPeiFd7YeFCxdG/jxu3DhlZ2dr6tSpamho0KhRo7p7zE71+H+Cy8jIUP/+/c/7FEtLS4uCwaDRVD3D0KFDde2116q+vt56FDNfHQMcH+cbOXKkMjIyeuXxsXjxYr311lt69913o359SzAY1KlTp3Ts2LGox/fW46Gr/dCZgoICSepRx0OPD1BKSoomTJigysrKyG0dHR2qrKxUYWGh4WT2jh8/roaGBmVnZ1uPYiYvL0/BYDDq+AiHw9q9e3efPz4+/fRTHT16tFcdH845LV68WBs3btSOHTuUl5cXdf+ECRM0cODAqOOhrq5OBw8e7FXHw8X2Q2f27t0rST3reLD+FMSlWL9+vfP7/W7NmjXu73//u1u4cKEbOnSoa25uth6tWz344IOuqqrKNTY2ur/+9a+uqKjIZWRkuCNHjliPllCtra3uo48+ch999JGT5J599ln30UcfuX//+9/OOeeeeuopN3ToULd582a3b98+N3PmTJeXl+e+/PJL48nj60L7obW11T300EOupqbGNTY2unfeecd973vfc9dcc407efKk9ehxs2jRIhcIBFxVVZU7fPhwZDtx4kTkMffee68bPny427Fjh9uzZ48rLCx0hYWFhlPH38X2Q319vfvlL3/p9uzZ4xobG93mzZvdyJEj3eTJk40nj5YUAXLOueeff94NHz7cpaSkuEmTJrldu3ZZj9Tt5syZ47Kzs11KSor75je/6ebMmePq6+utx0q4d99910k6b5s3b55z7uxHsR977DGXlZXl/H6/mzp1qqurq7MdOgEutB9OnDjhpk2b5q6++mo3cOBAN2LECLdgwYJe95e0zv77JbnVq1dHHvPll1+6++67z33jG99wQ4YMcbfffrs7fPiw3dAJcLH9cPDgQTd58mSXnp7u/H6/Gz16tPv5z3/uQqGQ7eDn4NcxAABM9Pj3gAAAvRMBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYOJ/Oa3sLCcwMkMAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "*** Image is classified as a 3 ***\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We will perform group LASSO on a multinomial regression in order to determin which of these pixel features are significant in identifying if the digit is 3, 5, or 8.  Intuitively it would make sense if the corners were not relevant for classifying the image, etc.  We will see the results of this process below."
      ],
      "metadata": {
        "id": "h2SS7hEsGhX0"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Part 2: Fit Group LASSO on Multinomial Regression"
      ],
      "metadata": {
        "id": "ZQikYR-d9QBN"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Here we perform a Multinomial Regression Group LASSO for groups. First we need to scale the data."
      ],
      "metadata": {
        "id": "M1_rtxQ9tQvR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Scale the data\n",
        "scaler = StandardScaler()\n",
        "X_train_scaled = scaler.fit_transform(X_train)\n",
        "X_test_scaled = scaler.transform(X_test)"
      ],
      "metadata": {
        "id": "AQd8WLiIre2r"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 1x1 Chunk Grouping"
      ],
      "metadata": {
        "id": "le__1AKlkT7H"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now initialize the model with #groups and l1 regularization.  Here we use the number of groups = number of features. This is done in order to determine the relevant features (or pixels) that are used across all digits 3, 5, and 8.  If we fit a regular LASSO, the features for each digit that shrink to 0 would differ across classes.  Now with group LASSO the features that shrink to 0 must do so across all classes which is a desirable feature.  Note, we could also use different groups such as 2x2 pixels which would likely increase accuracy; but I value interpretability and simplicity more in this fit so I just keep it simple and fit each feature (or pixel) as its own group."
      ],
      "metadata": {
        "id": "yaZrq1g-iRVf"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "print(f'# of groups = {X_train_scaled.shape[1]}')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QeOadi-xkgVF",
        "outputId": "4d10cd32-4f5c-4f51-bded-55afd0ffa5eb"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "# of groups = 784\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Initialize and train the model\n",
        "model = LogisticGroupLasso(groups=np.arange(X_train_scaled.shape[1]),\n",
        "                           group_reg=0.05,  # Group regularization strength\n",
        "                           l1_reg=0.01,     # L1 regularization within groups\n",
        "                           n_iter=1000,\n",
        "                           tol=1e-4,\n",
        "                           supress_warning=True\n",
        ")"
      ],
      "metadata": {
        "id": "3OdlKpkzh0P2"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Fit the model and calculate the accuracy (this takes 8 minutes to run; I will make more groups to improve the computational costs in the next part)."
      ],
      "metadata": {
        "id": "Gy6XC4RPisP7"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "model.fit(X_train_scaled, y_train) # Fit the model using scaled data\n",
        "y_pred = model.predict(X_test_scaled)\n",
        "y_prob = model.predict_proba(X_test_scaled)\n",
        "accuracy = accuracy_score(y_test, y_pred)\n",
        "report = classification_report(y_test, y_pred)\n",
        "\n",
        "print(f\"Accuracy: {accuracy:.4f}\")\n",
        "print(\"Classification Report:\")\n",
        "print(report)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "tCt-20snhxyp",
        "outputId": "d4fa0459-89f2-4c9d-dcb5-7b7dc44a2963"
      },
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Accuracy: 0.8084\n",
            "Classification Report:\n",
            "              precision    recall  f1-score   support\n",
            "\n",
            "           3       0.77      0.87      0.82      1426\n",
            "           5       0.80      0.76      0.78      1257\n",
            "           8       0.86      0.79      0.82      1373\n",
            "\n",
            "    accuracy                           0.81      4056\n",
            "   macro avg       0.81      0.81      0.81      4056\n",
            "weighted avg       0.81      0.81      0.81      4056\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Visualize the results of the group LASSO."
      ],
      "metadata": {
        "id": "IHEE-8OUjGrj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "coefficients = model.coef_ # Get the coefficients to see which groups were selected\n",
        "coef_df = pd.DataFrame(coefficients)\n",
        "coef_df.index = X_train.columns\n",
        "coef_df['sum'] = coef_df[0] + coef_df[1] + coef_df[2]\n",
        "bool_musk = coef_df['sum'] != 0\n",
        "nonzero_coef_df = coef_df.loc[bool_musk] #non-zero features\n",
        "relevant_features = list(nonzero_coef_df.index)\n",
        "print(f'There are {len(relevant_features)} nonzero features after applying group LASSO.')\n",
        "print('These are the pixels group LASSO found to be significant in classification of the digit to either 3, 5, or 8.')\n",
        "print(f'\\nPixel List = {relevant_features}')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "K0YnT5jqrbWZ",
        "outputId": "48355c97-59e0-47d0-eae2-ea105af6e5e3"
      },
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "There are 23 nonzero features after applying group LASSO.\n",
            "These are the pixels group LASSO found to be significant in classification of the digit to either 3, 5, or 8.\n",
            "\n",
            "Pixel List = ['pixel151', 'pixel152', 'pixel153', 'pixel177', 'pixel220', 'pixel291', 'pixel292', 'pixel297', 'pixel318', 'pixel319', 'pixel325', 'pixel327', 'pixel353', 'pixel380', 'pixel381', 'pixel408', 'pixel435', 'pixel462', 'pixel487', 'pixel488', 'pixel489', 'pixel515', 'pixel516']\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "coef_df['binary'] = coef_df['sum'].apply(lambda x: 0 if x == 0 else 1)\n",
        "image_matrix = coef_df['binary'].values.reshape((28, 28))\n",
        "plt.imshow(image_matrix, cmap='binary')\n",
        "plt.title(f'# of groups = {X_train_scaled.shape[1]}')\n",
        "\n",
        "black_patch = mpatches.Patch(color='black', label='Pixel On')\n",
        "white_patch = mpatches.Patch(color='white', label='Pixel Off')\n",
        "plt.legend(handles=[black_patch, white_patch], loc='upper right') # Add legend\n",
        "\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 452
        },
        "id": "tcqH-lHIa39w",
        "outputId": "7390034a-33a8-4a27-9370-9c9241773c69"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGzCAYAAABpdMNsAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALkhJREFUeJzt3XlYVfW+x/HPRmE7AJsQFEgknI+ZZg5kg2aSouWUlqmZNlgZWiR2yiaHBs5VKdNjds85T1Le1HMb1I6n7BCKmqEdPZUNNwcOTlcUtQTFGILf/cPHfduBA7jhB/h+Pc96Hvdv/fZa371Yz/641vrttRzGGCMAAKqZj+0CAACXJgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIItd7hw4c1YsQINWnSRA6HQ/PmzbNdEoALQAChynz99ddyOBzasWOHJOnVV1/VFVdc4fX1PP744/rkk080bdo0LVmyRHFxcV5fB/6fw+E463TLLbd49M3OztaDDz6o6OhoNWzYUK1atdKUKVN07Nixsy6/uLhYHTp0kMPh0Ny5c6v648Ci+rYLQN21ZcsWBQcHq23btpKkjIwMXXvttV5fz9q1azVkyBBNnTrV68tGWUuWLCnTtnXrVr322mvq16+fu+3kyZPq2bOn8vPz9cgjjygyMlJff/21/vjHP2rdunXatm2bfHzK/h94wYIF2rdvX5V+BtQMBBCqzBdffKEePXrI4XBIOh1AU6ZM8fp6cnJyFBQU5PXl/lZ+fr4aN25c5eup6e6+++4ybenp6XI4HBo1apS77cMPP9TevXu1evVq3Xrrre724OBgzZo1S19//bW6dOnisZycnBzNmjVLTz75pJ5//vmq+xCoETgFB6/66aefdPToUR09elRbtmxRx44ddfToUX333Xc6cOCA2rRpo6NHj+rkyZPnXda///1v3XHHHQoODlajRo107bXX6u9//7t7fkpKihwOh4wxWrhwofs00LkcO3ZMY8eOVWBgoIKCgjRu3Dj3qcKUlBR3v/Hjx8vf31+ZmZkaOHCgAgICNGbMGEmngygxMVGRkZFyOp1q166d5s6dq1/fWH7Pnj1llnmGw+HQjBkz3K9nzJghh8OhH374QXfeeacCAwPVpEkTPfbYYyooKPB4b2pqqm644QYFBQXJ399f7dq109NPP33ebVmVCgsL9f7776t3795q3ry5uz0vL0+S1KxZM4/+4eHhkqSGDRuWWdZTTz2ldu3alRtyqHs4AoJXdenSRXv37nW//vbbbz3O4w8aNEiSNG7cuHK/nM84fPiwrrvuOp06dUqPPvqomjRporfeekuDBw/We++9p2HDhqlXr15asmSJxo4dq1tuuUX33HPPOWsrLS3VoEGD9MUXX2jixIlq3769Vq1apXHjxpXb/5dfflH//v11ww03aO7cuWrUqJGMMRo8eLDWrVun+++/X1dffbU++eQTPfHEE/rf//1fvfrqqxXYWp7uvPNOXXHFFUpKStLmzZs1f/58/fTTT3r77bclSd99951uu+02derUSbNmzZLT6dTu3bu1adOm8y47NzdXxcXF5+3XoEED+fv7V6jujz76SMePH3cH9Bm9evWSj4+PHnvsMSUnJ6t58+bavn27XnrpJQ0dOlTt27f36P/FF1/orbfe0meffXbe/0igjjCAF3322WcmNTXVPPfcc6Z+/frm448/NqmpqWbAgAGmW7duJjU11aSmpprvvvvunMtJSEgwkszGjRvdbSdOnDDR0dHmiiuuMCUlJe52SSY+Pv68tb3//vtGkpk3b567raSkxNx8881Gklm8eLG7fdy4cUaSeeqppzyWsXLlSiPJvPjiix7tI0aMMA6Hw+zevdsYY0xWVlaZZf663unTp7tfT58+3UgygwcP9uj3yCOPGEnm66+/NsYY8+qrrxpJ5siRI+f9rL/Vu3dvI+m807hx4yq87OHDhxun02l++umnMvP+8pe/mKCgoDLrKC4u9uhXWlpqevToYUaNGmWM+f/tN2fOnArXg9qDIyB41fXXXy/p9P+Ku3fv7h6RlpCQoDvuuEOxsbEXtJyPPvpIPXr00A033OBu8/f314MPPqhp06bp+++/V8eOHStU25o1a+Tr66sJEya423x8fBQfH6+1a9eW+56JEyeWqatevXp69NFHPdoTExP13nvv6eOPP9akSZMqVNcZ8fHxHq8nT56s119/XR999JE6derkvs61atUq3XvvveVewD+b5ORk/fTTT+ftFxERUaGa8/Ly9Pe//10DBw4s9zrc5Zdfrh49emjgwIGKiorSxo0bNX/+fIWEhHgcGaekpOibb77Re++9V6H1o3YjgOA1vz7Nk5aWpptvvllHjx7Vjz/+qO+++04vvviijh49Kl9fX7lcrnMua+/evYqJiSnT/rvf/c49v6IBtHfvXoWHh6tRo0Ye7a1bty63f/369T2uaZxZRkREhAICAs5aV2W1adPG43WrVq3k4+OjPXv2SJJGjhypv/zlL3rggQf01FNPqW/fvrr99ts1YsSI84ZR165dK13Xubz//vsqKCgoc/pNkjZt2qTbbrtNmzdvVrdu3SRJQ4cOVWBgoGbOnKn77rtPHTp0UF5enqZNm6YnnnhCkZGRVVInaiYGIcBrhgwZotDQUIWGhmr79u2aN2+eQkND1a5dO0nSsGHDFBoaqiFDhliu9MI4nc4KHWX82tmuYZSUlFR6GQ0bNtSGDRv06aefauzYsdq+fbtGjhypW2655bzL/fHHH3Xo0KHzTrm5uRdcnyS98847crlcuu2228rM+8///E81a9bMHT5nDB48WMYYff7555KkuXPnqqioSCNHjtSePXu0Z88eHThwQNLpQS179uxRUVFRhepC7UAAwWuSk5OVmpqqWbNmqV69evr444+VmpqqwYMHq2vXrkpNTVVqaqqSk5PPu6yoqCj3D1h/7YcffnDPr6ioqChlZ2fr1KlTHu27d++u0DIOHjyoEydOnLOuyy67TJJ0/Phxj37nOkLatWtXmbpKS0s9frzr4+Ojvn376pVXXtH333+vl156SWvXrtW6devOWfftt9+u8PDw806PPfbYOZfza9nZ2Vq3bp2GDx8up9NZZv7hw4fLDcYzR8m//PKLJGnfvn366aefdOWVVyo6OlrR0dG68cYbJUkvv/yyoqOj9f33319wXag9OAUHrzlzmmfjxo3q2LGj+/rPs88+q9jY2Au+/iNJAwcO1Lx585SRkaGePXtKOj38+U9/+pOuuOIKdejQocL19e/fX3/+85/15z//2f1FW1paqoULF1aorj/96U/64x//qGnTprnbX331VTkcDg0YMECSFBgYqJCQEG3YsEEJCQnufq+//vpZl71w4UKPH3IuWLBAktzL/PHHHxUcHOzxnquvvlrS6aHQ51IV14CWL1+u0tLSck+/SVLbtm31j3/8Q+np6brpppvc7cuWLZMk92+AHn30UQ0dOtTjvTk5OXrooYc0fvx4DRkyRNHR0RdcF2oPAghet2nTJl133XWSpIKCAn355ZcV/q3KU089pWXLlmnAgAF69NFHFRwcrLfeektZWVl6//33K3VqbOjQoerRo4cSExO1e/dutW/fXh9++KF+/PFHSWc/bfZrgwYNUp8+ffTMM89oz5496ty5s/7xj39o1apVSkhIUKtWrdx9H3jgAf3hD3/QAw88oG7dumnDhg3auXPnWZedlZWlwYMHKy4uThkZGfqv//ovjR49Wp07d5YkzZo1Sxs2bNCtt96qqKgo5eTk6PXXX1fz5s09BmuUpyquAb3zzjuKiIjwCJdfmzRpkhYvXqxBgwZp8uTJioqK0vr167Vs2TLdcsst7mt811xzja655hqP95657nXllVeWCSfUIbaH4aFu+eWXX4y/v79ZsmSJMeb0sGxJJicnp8LLyszMNCNGjDBBQUGmQYMGpkePHmb16tVl+ukCh2EbY8yRI0fM6NGjTUBAgHG5XGb8+PFm06ZNRpJZvny5u9+4ceNM48aNy13GiRMnzOOPP24iIiKMr6+vadOmjZkzZ44pLS316Hfq1Clz//33G5fLZQICAsydd95pcnJyzjoM+/vvvzcjRowwAQEB5rLLLjOTJk0yP//8s7tfWlqaGTJkiImIiDB+fn4mIiLCjBo1yuzcufOCPrs3/fDDD0aSmTJlynn7jRgxwkRGRhpfX18TFRVlpk6davLz88/5PoZhXxocxvzq59vAJWjlypUaNmyYPvvsM/cw8uo0Y8YMzZw5U0eOHFFISEi1rx+whUEIuKT8/PPPHq9LSkq0YMECBQYGljkNBKBqcQ0Il5TJkyfr559/Vs+ePVVYWKgPPvhAn3/+uV5++eVy700GoOoQQLik3HzzzUpOTtbq1atVUFCg1q1ba8GCBZW+ewGAyuMaEADACq4BAQCsIIAAAFbUuGtApaWlOnjwoAICAngmCADUQsYYnThxQhEREef80XiNC6CDBw9yR1wAqAP2799f5o7yv1bjAujMbe7379+vwMBAy9UAACoqLy9PkZGRZR5b8ltVFkALFy7UnDlzdOjQIXXu3FkLFixQjx49zvu+M6fdAgMDCSAAqMXOdxmlSgYh/PWvf9WUKVM0ffp0/etf/1Lnzp3Vv39/5eTkVMXqAAC1UJUE0CuvvKIJEybo3nvvVYcOHfTGG2+oUaNGevPNN6tidQCAWsjrAVRUVKRt27Z5PPvFx8dHsbGxysjIKNO/sLBQeXl5HhMAoO7z+jWgo0ePqqSkRM2aNfNob9asmfupkb+WlJSkmTNnersMALVYSUmJ+8mpqHnq1aun+vXrX/RPZayPgps2bZqmTJnifn1m9ASAS9PJkyd14MABcZewmq1Ro0YKDw+Xn59fpZfh9QAKCQlRvXr1dPjwYY/2w4cPKywsrEx/p9NZ7vPkAVx6SkpKdODAATVq1EihoaH8GL0GMsaoqKhIR44cUVZWltq0aVOpJxRLVRBAfn5+6tq1q9LS0tyP0i0tLVVaWhp3HAZwTsXFxTLGKDQ0lMdj1GANGzaUr6+v9u7dq6KiIjVo0KBSy6mSU3BTpkzRuHHj1K1bN/Xo0UPz5s1Tfn6+7r333qpYHYA6hiOfmq+yRz2/ViUBNHLkSB05ckTPP/+8Dh06pKuvvlpr1qwpMzABAHDpqrJBCJMmTeKUGwCv2Ldvn44ePVpt6wsJCVGLFi2qbX2XKuuj4ADgXPbt26d27dqpoKCg2tbZoEED7dixw2shNH78eB0/flwrV670yvJSUlKUkJCg48ePe2V5tvA8IAA12tGjR6s1fCSpoKCgQkdc48ePl8PhkMPhkJ+fn1q3bq1Zs2bpl19+kSS99tprSklJqaJqz2716tXq3bu3AgIC1KhRI3Xv3t1KHWdDAAGAF8TFxSk7O1u7du1SYmKiZsyYoTlz5kiSXC6XgoKCqrWeBQsWaMiQIbr++uu1ZcsWbd++XXfddZcefvhhTZ06tVprORsCCAC8wOl0KiwsTFFRUZo4caJiY2P14YcfSjp9hHTmZylHjhxRWFiYXn75Zfd7P//8c/n5+SktLU3S6VuUTZ06VZdffrkaN26smJgYpaenX3At+/fvV2JiohISEvTyyy+rQ4cOat26tRITEzVnzhwlJydry5YtkqT09HQ5HA6lpaWpW7duatSoka677jrt2LHDOxvmHAggAKgCDRs2VFFRUZn20NBQvfnmm5oxY4a2bt2qEydOaOzYsZo0aZL69u0r6fQgroyMDC1fvlzbt2/XHXfcobi4OO3ateuC1v3ee++puLi43COdhx56SP7+/lq2bJlH+zPPPKPk5GRt3bpV9evX13333VeJT10xDEIAAC8yxigtLU2ffPKJJk+eXG6fgQMHasKECRozZoy6deumxo0bKykpSdLpQReLFy/Wvn37FBERIUmaOnWq1qxZo8WLF3scOZ3Nzp075XK5FB4eXmaen5+fWrZsqZ07d3q0v/TSS+rdu7ck6amnntKtt96qgoKCSv/I9EIQQADgBatXr5a/v7+Ki4tVWlqq0aNHa8aMGWftP3fuXHXs2FHvvvuutm3b5r4l2TfffKOSkhK1bdvWo39hYaGaNGlSZfV36tTJ/e8zwZWTk1Olw9EJIADwgj59+mjRokXy8/NTRESE6tc/99drZmamDh48qNLSUu3Zs0dXXXWVpNM3Y61Xr562bdumevXqebzH39//gmpp27atcnNzdfDgQfdR1BlFRUXKzMxUnz59PNp9fX3d/z5zJ4rS0tILWl9lcQ0IALygcePGat26tVq0aHHe8CkqKtLdd9+tkSNH6oUXXtADDzzgfmJ0ly5dVFJSopycHLVu3dpjKu+GzuUZPny4fH19lZycXGbeG2+8ofz8fI0aNariH9LLOAICgGr2zDPPKDc3V/Pnz5e/v78++ugj3XfffVq9erXatm2rMWPG6J577lFycrK6dOmiI0eOKC0tTZ06ddKtt9563uW3aNFCs2fPVmJioho0aKCxY8fK19dXq1at0tNPP63ExETFxMRUwyc9N46AANRoISEhVXohvDwNGjRQSEhIlSw7PT1d8+bN05IlSxQYGCgfHx8tWbJEGzdu1KJFiyRJixcv1j333KPExES1a9dOQ4cO1T//+c8KXY9JSEjQihUrtHHjRnXr1k0dO3bU0qVLtWjRIs2dO7dKPltFOUwNe+pTXl6eXC6XcnNzFRgYaLscANWooKBAWVlZio6O9ggd7gVX85ztbyVd+Pc4p+AA1HgtWrQgEOogTsEBAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABQBX79RNRvSElJcWrj/jetGmTrrrqKvn6+rrrLK/N2wggALhI48ePl8PhkMPhkJ+fn1q3bq1Zs2bpl19+kSS99tprSklJqfa6Vq9erd69eysgIECNGjVS9+7dy61jypQpuvrqq5WVleWeX16btxFAAOAFcXFxys7O1q5du5SYmKgZM2Zozpw5kiSXy+XVI5YLsWDBAg0ZMkTXX3+9tmzZou3bt+uuu+7Sww8/XOZR3ZmZmbr55pvVvHlzd53ltXkbAQQAXuB0OhUWFqaoqChNnDhRsbGx+vDDDyV5noI7cuSIwsLCPB6t/fnnn8vPz09paWmSTj/9dOrUqbr88svVuHFjxcTEKD09/YJr2b9/vxITE5WQkKCXX35ZHTp0UOvWrZWYmKg5c+YoOTlZW7Zs0Z49e+RwOHTs2DHdd999cjgcSklJKbetKhBAAFAFGjZsqKKiojLtoaGhevPNNzVjxgxt3bpVJ06c0NixYzVp0iT17dtXkjRp0iRlZGRo+fLl2r59u+644w7FxcVp165dF7Tu9957T8XFxWWOdCTpoYcekr+/v5YtW6bIyEhlZ2crMDBQ8+bNU3Z2tu64444ybSNHjry4jXEW3A0bALzIGKO0tDR98sknmjx5crl9Bg4cqAkTJmjMmDHq1q2bGjdurKSkJEmnHz2xePFi7du3z/047alTp2rNmjVavHixx5HT2ezcuVMul0vh4eFl5vn5+ally5bauXOn6tWrp7CwMDkcDrlcLvcTVxs3blymrSoQQADgBatXr5a/v7+Ki4tVWlqq0aNHa8aMGWftP3fuXHXs2FHvvvuutm3bJqfTKUn65ptvVFJSorZt23r0LywsVJMmTaryI1Q7AggAvKBPnz5atGiR/Pz8FBERofr1z/31mpmZqYMHD6q0tFR79uzRVVddJUk6efKk6tWrp23btqlevXoe7/H397+gWtq2bavc3FwdPHjQfRR1RlFRkTIzM9WnT58KfLqqwTUgAPCCxo0bq3Xr1mrRosV5w6eoqEh33323Ro4cqRdeeEEPPPCAcnJyJEldunRRSUmJcnJy1Lp1a4/pQk+HDR8+XL6+vkpOTi4z74033lB+fr5GjRpV8Q/pZRwBAUA1e+aZZ5Sbm6v58+fL399fH330ke677z6tXr1abdu21ZgxY3TPPfcoOTlZXbp00ZEjR5SWlqZOnTrp1ltvPe/yW7RoodmzZysxMVENGjTQ2LFj5evrq1WrVunpp59WYmKiYmJiquGTnhtHQABQjdLT0zVv3jwtWbJEgYGB8vHx0ZIlS7Rx40YtWrRIkrR48WLdc889SkxMVLt27TR06FD985//rNBjyRMSErRixQpt3LhR3bp1U8eOHbV06VItWrRIc+fOraqPVyEOY4yxXcSv5eXlyeVyKTc3V4GBgbbLAVCNCgoKlJWVpejoaDVo0MB2OTiHc/2tLvR7nCMgAIAVBBAAwAoCCABgBQEEALCCAAJQ49SwsVEohzf+RgQQgBrjzC//y7uJJ2qWU6dOSZJ8fX0rvQx+iAqgxqhfv74aNWqkI0eOyNfXVz4+/B+5pjHG6NSpU8rJyVFQUFCZ2wVVBAEEoMZwOBwKDw9XVlaW9u7da7scnENQUNBF3ymbAAJQo/j5+alNmzachqvBfH19L+rI5wwCCECN4+Pjw50QLgGcYAUAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWMGdEFDjORwO2yWcE48OACqHIyAAgBUEEADACq8H0IwZM+RwODym9u3be3s1AIBarkquAV155ZX69NNP/38l9bnUBADwVCXJUL9+/Yt+UBEAoG6rkmtAu3btUkREhFq2bKkxY8Zo3759Z+1bWFiovLw8jwkAUPd5PYBiYmKUkpKiNWvWaNGiRcrKytKNN96oEydOlNs/KSlJLpfLPUVGRnq7JABADeQwVfwjhuPHjysqKkqvvPKK7r///jLzCwsLVVhY6H6dl5enyMhI5ebmKjAwsCpLQy3B74CA2iUvL08ul+u83+NVPjogKChIbdu21e7du8ud73Q65XQ6q7oMAEANU+W/Azp58qQyMzMVHh5e1asCANQiXg+gqVOnav369dqzZ48+//xzDRs2TPXq1dOoUaO8vSoAQC3m9VNwBw4c0KhRo3Ts2DGFhobqhhtu0ObNmxUaGurtVQEAajGvB9Dy5cu9vUjUUJUZHFCZC/Zc5AfqJu4FBwCwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWVPkD6VB3cZNQoKzquklvXcAREADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKzgbtgA4EWX6p2tK4MjIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFfVtFwCcj8PhqLZ1GWOqbV3VobLbrq5tB9RMHAEBAKwggAAAVlQ4gDZs2KBBgwYpIiJCDodDK1eu9JhvjNHzzz+v8PBwNWzYULGxsdq1a5e36gUA1BEVDqD8/Hx17txZCxcuLHf+7NmzNX/+fL3xxhvasmWLGjdurP79+6ugoOCiiwUA1B0VHoQwYMAADRgwoNx5xhjNmzdPzz77rIYMGSJJevvtt9WsWTOtXLlSd91118VVCwCoM7x6DSgrK0uHDh1SbGysu83lcikmJkYZGRnlvqewsFB5eXkeEwCg7vNqAB06dEiS1KxZM4/2Zs2auef9VlJSklwul3uKjIz0ZkkAgBrK+ii4adOmKTc31z3t37/fdkkAgGrg1QAKCwuTJB0+fNij/fDhw+55v+V0OhUYGOgxAQDqPq8GUHR0tMLCwpSWluZuy8vL05YtW9SzZ09vrgoAUMtVeBTcyZMntXv3bvfrrKwsffXVVwoODlaLFi2UkJCgF198UW3atFF0dLSee+45RUREaOjQod6sGwBQy1U4gLZu3ao+ffq4X0+ZMkWSNG7cOKWkpOj3v/+98vPz9eCDD+r48eO64YYbtGbNGjVo0MB7VQMAaj2HqWF3HczLy5PL5VJubi7Xg+qg6rqxaA3bra2pzpuRVmZd/J3qpgv9Hrc+Cg4AcGkigAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADAigo/jgG4GNz9uO7ib4uK4ggIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKzgZqRALeFwOGyXAHgVR0AAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAU3IwXqMGOM7RKAs+IICABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCs4GakgAUOh6PC7+HGoqhrOAICAFhBAAEArKhwAG3YsEGDBg1SRESEHA6HVq5c6TF//PjxcjgcHlNcXJy36gUA1BEVDqD8/Hx17txZCxcuPGufuLg4ZWdnu6dly5ZdVJEAgLqnwoMQBgwYoAEDBpyzj9PpVFhYWKWLAgDUfVVyDSg9PV1NmzZVu3btNHHiRB07duysfQsLC5WXl+cxAQDqPq8HUFxcnN5++22lpaXpP/7jP7R+/XoNGDBAJSUl5fZPSkqSy+VyT5GRkd4uCQBQAznMRfy4wOFwaMWKFRo6dOhZ+/z73/9Wq1at9Omnn6pv375l5hcWFqqwsND9Oi8vT5GRkcrNzVVgYGBlSwNqNH4HhLosLy9PLpfrvN/jVT4Mu2XLlgoJCdHu3bvLne90OhUYGOgxAQDqvioPoAMHDujYsWMKDw+v6lUBAGqRCo+CO3nypMfRTFZWlr766isFBwcrODhYM2fO1PDhwxUWFqbMzEz9/ve/V+vWrdW/f3+vFg4AqN0qHEBbt25Vnz593K+nTJkiSRo3bpwWLVqk7du366233tLx48cVERGhfv366YUXXpDT6fRe1QCAWu+iBiFUhQu9eAXUFAwoADzVmEEIAACUhwACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsq/DgGoLpV5m7T1Yk7WwOVwxEQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFjBzUiBX+HGokD14QgIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKzgZqSo8bhBKFA3cQQEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwIoKBVBSUpK6d++ugIAANW3aVEOHDtWOHTs8+hQUFCg+Pl5NmjSRv7+/hg8frsOHD3u1aABA7VehAFq/fr3i4+O1efNmpaamqri4WP369VN+fr67z+OPP66//e1vevfdd7V+/XodPHhQt99+u9cLBwDUbg5zEY+bPHLkiJo2bar169erV69eys3NVWhoqJYuXaoRI0ZIkn744Qf97ne/U0ZGhq699trzLjMvL08ul0u5ubkKDAysbGkAAEsu9Hv8oq4B5ebmSpKCg4MlSdu2bVNxcbFiY2Pdfdq3b68WLVooIyOj3GUUFhYqLy/PYwIA1H2VDqDS0lIlJCTo+uuvV8eOHSVJhw4dkp+fn4KCgjz6NmvWTIcOHSp3OUlJSXK5XO4pMjKysiUBAGqRSgdQfHy8vv32Wy1fvvyiCpg2bZpyc3Pd0/79+y9qeQCA2qF+Zd40adIkrV69Whs2bFDz5s3d7WFhYSoqKtLx48c9joIOHz6ssLCwcpfldDrldDorUwYAoBar0BGQMUaTJk3SihUrtHbtWkVHR3vM79q1q3x9fZWWluZu27Fjh/bt26eePXt6p2IAQJ1QoSOg+Ph4LV26VKtWrVJAQID7uo7L5VLDhg3lcrl0//33a8qUKQoODlZgYKAmT56snj17XtAIOADApaNCw7AdDke57YsXL9b48eMlnf4hamJiopYtW6bCwkL1799fr7/++llPwf0Ww7ABoHa70O/xi/odUFUggACgdquW3wEBAFBZBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFZUKICSkpLUvXt3BQQEqGnTpho6dKh27Njh0eemm26Sw+HwmB5++GGvFg0AqP0qFEDr169XfHy8Nm/erNTUVBUXF6tfv37Kz8/36DdhwgRlZ2e7p9mzZ3u1aABA7Ve/Ip3XrFnj8TolJUVNmzbVtm3b1KtXL3d7o0aNFBYW5p0KAQB10kVdA8rNzZUkBQcHe7S/8847CgkJUceOHTVt2jSdOnXqrMsoLCxUXl6exwQAqPsqdAT0a6WlpUpISND111+vjh07uttHjx6tqKgoRUREaPv27XryySe1Y8cOffDBB+UuJykpSTNnzqxsGQCAWsphjDGVeePEiRP18ccf67PPPlPz5s3P2m/t2rXq27evdu/erVatWpWZX1hYqMLCQvfrvLw8RUZGKjc3V4GBgZUpDQBgUV5enlwu13m/xyt1BDRp0iStXr1aGzZsOGf4SFJMTIwknTWAnE6nnE5nZcoAANRiFQogY4wmT56sFStWKD09XdHR0ed9z1dffSVJCg8Pr1SBAIC6qUIBFB8fr6VLl2rVqlUKCAjQoUOHJEkul0sNGzZUZmamli5dqoEDB6pJkybavn27Hn/8cfXq1UudOnWqkg8AAKidKnQNyOFwlNu+ePFijR8/Xvv379fdd9+tb7/9Vvn5+YqMjNSwYcP07LPPXvD1nAs9dwgAqJmq5BrQ+bIqMjJS69evr8giAQCXKO4FBwCwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwor7tAn7LGCNJysvLs1wJAKAyznx/n/k+P5saF0AnTpyQJEVGRlquBABwMU6cOCGXy3XW+Q5zvoiqZqWlpTp48KACAgLkcDg85uXl5SkyMlL79+9XYGCgpQrtYzucxnY4je1wGtvhtJqwHYwxOnHihCIiIuTjc/YrPTXuCMjHx0fNmzc/Z5/AwMBLegc7g+1wGtvhNLbDaWyH02xvh3Md+ZzBIAQAgBUEEADAiloVQE6nU9OnT5fT6bRdilVsh9PYDqexHU5jO5xWm7ZDjRuEAAC4NNSqIyAAQN1BAAEArCCAAABWEEAAACsIIACAFbUmgBYuXKgrrrhCDRo0UExMjL744gvbJVW7GTNmyOFweEzt27e3XVaV27BhgwYNGqSIiAg5HA6tXLnSY74xRs8//7zCw8PVsGFDxcbGateuXXaKrULn2w7jx48vs3/ExcXZKbaKJCUlqXv37goICFDTpk01dOhQ7dixw6NPQUGB4uPj1aRJE/n7+2v48OE6fPiwpYqrxoVsh5tuuqnM/vDwww9bqrh8tSKA/vrXv2rKlCmaPn26/vWvf6lz587q37+/cnJybJdW7a688kplZ2e7p88++8x2SVUuPz9fnTt31sKFC8udP3v2bM2fP19vvPGGtmzZosaNG6t///4qKCio5kqr1vm2gyTFxcV57B/Lli2rxgqr3vr16xUfH6/NmzcrNTVVxcXF6tevn/Lz8919Hn/8cf3tb3/Tu+++q/Xr1+vgwYO6/fbbLVbtfReyHSRpwoQJHvvD7NmzLVV8FqYW6NGjh4mPj3e/LikpMRERESYpKcliVdVv+vTppnPnzrbLsEqSWbFihft1aWmpCQsLM3PmzHG3HT9+3DidTrNs2TILFVaP324HY4wZN26cGTJkiJV6bMnJyTGSzPr1640xp//2vr6+5t1333X3+Z//+R8jyWRkZNgqs8r9djsYY0zv3r3NY489Zq+oC1Djj4CKioq0bds2xcbGutt8fHwUGxurjIwMi5XZsWvXLkVERKhly5YaM2aM9u3bZ7skq7KysnTo0CGP/cPlcikmJuaS3D/S09PVtGlTtWvXThMnTtSxY8dsl1SlcnNzJUnBwcGSpG3btqm4uNhjf2jfvr1atGhRp/eH326HM9555x2FhISoY8eOmjZtmk6dOmWjvLOqcXfD/q2jR4+qpKREzZo182hv1qyZfvjhB0tV2RETE6OUlBS1a9dO2dnZmjlzpm688UZ9++23CggIsF2eFYcOHZKkcvePM/MuFXFxcbr99tsVHR2tzMxMPf300xowYIAyMjJUr1492+V5XWlpqRISEnT99derY8eOkk7vD35+fgoKCvLoW5f3h/K2gySNHj1aUVFRioiI0Pbt2/Xkk09qx44d+uCDDyxW66nGBxD+34ABA9z/7tSpk2JiYhQVFaX//u//1v3332+xMtQEd911l/vfV111lTp16qRWrVopPT1dffv2tVhZ1YiPj9e33357SVwHPZezbYcHH3zQ/e+rrrpK4eHh6tu3rzIzM9WqVavqLrNcNf4UXEhIiOrVq1dmFMvhw4cVFhZmqaqaISgoSG3bttXu3bttl2LNmX2A/aOsli1bKiQkpE7uH5MmTdLq1au1bt06j+eHhYWFqaioSMePH/foX1f3h7Nth/LExMRIUo3aH2p8APn5+alr165KS0tzt5WWliotLU09e/a0WJl9J0+eVGZmpsLDw22XYk10dLTCwsI89o+8vDxt2bLlkt8/Dhw4oGPHjtWp/cMYo0mTJmnFihVau3atoqOjPeZ37dpVvr6+HvvDjh07tG/fvjq1P5xvO5Tnq6++kqSatT/YHgVxIZYvX26cTqdJSUkx33//vXnwwQdNUFCQOXTokO3SqlViYqJJT083WVlZZtOmTSY2NtaEhISYnJwc26VVqRMnTpgvv/zSfPnll0aSeeWVV8yXX35p9u7da4wx5g9/+IMJCgoyq1atMtu3bzdDhgwx0dHR5ueff7ZcuXedazucOHHCTJ061WRkZJisrCzz6aefmmuuuca0adPGFBQU2C7dayZOnGhcLpdJT0832dnZ7unUqVPuPg8//LBp0aKFWbt2rdm6davp2bOn6dmzp8Wqve9822H37t1m1qxZZuvWrSYrK8usWrXKtGzZ0vTq1cty5Z5qRQAZY8yCBQtMixYtjJ+fn+nRo4fZvHmz7ZKq3ciRI014eLjx8/Mzl19+uRk5cqTZvXu37bKq3Lp164ykMtO4ceOMMaeHYj/33HOmWbNmxul0mr59+5odO3bYLboKnGs7nDp1yvTr18+EhoYaX19fExUVZSZMmFDn/pNW3ueXZBYvXuzu8/PPP5tHHnnEXHbZZaZRo0Zm2LBhJjs7217RVeB822Hfvn2mV69eJjg42DidTtO6dWvzxBNPmNzcXLuF/wbPAwIAWFHjrwEBAOomAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACw4v8AOnolCzLUbAkAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Define Functions (Testing Groups)"
      ],
      "metadata": {
        "id": "VSgbZeJf-__K"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def fit_summary(groups):\n",
        "    # Initialize and train the model\n",
        "  model = LogisticGroupLasso(groups=groups,\n",
        "                            group_reg=0.05,  # Group regularization strength\n",
        "                            l1_reg=0.01,     # L1 regularization within groups\n",
        "                            n_iter=1000,\n",
        "                            tol=1e-4,\n",
        "                            supress_warning=True\n",
        "  )\n",
        "  model.fit(X_train_scaled, y_train) # Fit the model using scaled data\n",
        "  y_pred = model.predict(X_test_scaled)\n",
        "  y_prob = model.predict_proba(X_test_scaled)\n",
        "  accuracy = accuracy_score(y_test, y_pred)\n",
        "  report = classification_report(y_test, y_pred)\n",
        "\n",
        "  print(f\"Accuracy: {accuracy:.4f}\")\n",
        "  print(\"Classification Report:\")\n",
        "  print(report)\n",
        "\n",
        "  coefficients = model.coef_ # Get the coefficients to see which groups were selected\n",
        "  coef_df = pd.DataFrame(coefficients)\n",
        "  coef_df.index = X_train.columns\n",
        "  coef_df['sum'] = coef_df[0] + coef_df[1] + coef_df[2]\n",
        "  bool_musk = coef_df['sum'] != 0\n",
        "  nonzero_coef_df = coef_df.loc[bool_musk] #non-zero features\n",
        "  relevant_features = list(nonzero_coef_df.index)\n",
        "  print(f'There are {len(relevant_features)} nonzero features after applying group LASSO.')\n",
        "  print('These are the pixels group LASSO found to be significant in classification of the digit to either 3, 5, or 8.')\n",
        "  print(f'\\nPixel List = {relevant_features}')\n",
        "  coef_df['binary'] = coef_df['sum'].apply(lambda x: 0 if x == 0 else 1)\n",
        "  image_matrix = coef_df['binary'].values.reshape((28, 28))\n",
        "  plt.imshow(image_matrix, cmap='binary')\n",
        "  plt.title(f'# of groups = {len(np.unique(groups))}')\n",
        "\n",
        "  black_patch = mpatches.Patch(color='black', label='Pixel On')\n",
        "  white_patch = mpatches.Patch(color='white', label='Pixel Off')\n",
        "  plt.legend(handles=[black_patch, white_patch], loc='upper right') # Add legend\n",
        "\n",
        "  plt.show()"
      ],
      "metadata": {
        "id": "k0XIKPOi_Erb"
      },
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now we try to improve the results by using fewer groups:\n",
        "\n",
        "- Encourages more aggressive feature selection (more features dropped together).\n",
        "- Can improve interpretability by selecting entire feature sets.\n",
        "- Helps when features within a group are highly correlated."
      ],
      "metadata": {
        "id": "4X5oVlmAnBYP"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 2x1 Chunk Grouping"
      ],
      "metadata": {
        "id": "aOhf3E-vsxa8"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Create a 28x28 matrix with paired numbers that continue across rows\n",
        "rows = 28\n",
        "cols = 28\n",
        "pairs_per_row = cols // 2  # How many pairs in each row\n",
        "\n",
        "# Initialize the matrix\n",
        "groups = np.zeros((rows, cols), dtype=int)\n",
        "\n",
        "# Fill the matrix with pairs of continuing numbers\n",
        "for i in range(rows):\n",
        "    row_numbers = []\n",
        "    start_number = i * pairs_per_row + 1\n",
        "    for j in range(pairs_per_row):\n",
        "        current_number = start_number + j\n",
        "        row_numbers.extend([current_number, current_number])\n",
        "    groups[i] = row_numbers\n"
      ],
      "metadata": {
        "id": "-e4PZXKds7xN"
      },
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The groups can be visualized by 1x2 chunks in a 28x28 image.  As given below:"
      ],
      "metadata": {
        "id": "zgR2HbEVwM2d"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Print the matrix with proper alignment\n",
        "for row in groups:\n",
        "    formatted_row = []\n",
        "    for num in row:\n",
        "        # Format each individual number, not the entire array\n",
        "        formatted_row.append(f\"{num:3d}\")\n",
        "    print(\"[\" + \" \".join(formatted_row) + \"]\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BzqFmqQuvikA",
        "outputId": "39e4e25c-dd47-4cfd-ef21-01ed92e1d83d"
      },
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[  1   1   2   2   3   3   4   4   5   5   6   6   7   7   8   8   9   9  10  10  11  11  12  12  13  13  14  14]\n",
            "[ 15  15  16  16  17  17  18  18  19  19  20  20  21  21  22  22  23  23  24  24  25  25  26  26  27  27  28  28]\n",
            "[ 29  29  30  30  31  31  32  32  33  33  34  34  35  35  36  36  37  37  38  38  39  39  40  40  41  41  42  42]\n",
            "[ 43  43  44  44  45  45  46  46  47  47  48  48  49  49  50  50  51  51  52  52  53  53  54  54  55  55  56  56]\n",
            "[ 57  57  58  58  59  59  60  60  61  61  62  62  63  63  64  64  65  65  66  66  67  67  68  68  69  69  70  70]\n",
            "[ 71  71  72  72  73  73  74  74  75  75  76  76  77  77  78  78  79  79  80  80  81  81  82  82  83  83  84  84]\n",
            "[ 85  85  86  86  87  87  88  88  89  89  90  90  91  91  92  92  93  93  94  94  95  95  96  96  97  97  98  98]\n",
            "[ 99  99 100 100 101 101 102 102 103 103 104 104 105 105 106 106 107 107 108 108 109 109 110 110 111 111 112 112]\n",
            "[113 113 114 114 115 115 116 116 117 117 118 118 119 119 120 120 121 121 122 122 123 123 124 124 125 125 126 126]\n",
            "[127 127 128 128 129 129 130 130 131 131 132 132 133 133 134 134 135 135 136 136 137 137 138 138 139 139 140 140]\n",
            "[141 141 142 142 143 143 144 144 145 145 146 146 147 147 148 148 149 149 150 150 151 151 152 152 153 153 154 154]\n",
            "[155 155 156 156 157 157 158 158 159 159 160 160 161 161 162 162 163 163 164 164 165 165 166 166 167 167 168 168]\n",
            "[169 169 170 170 171 171 172 172 173 173 174 174 175 175 176 176 177 177 178 178 179 179 180 180 181 181 182 182]\n",
            "[183 183 184 184 185 185 186 186 187 187 188 188 189 189 190 190 191 191 192 192 193 193 194 194 195 195 196 196]\n",
            "[197 197 198 198 199 199 200 200 201 201 202 202 203 203 204 204 205 205 206 206 207 207 208 208 209 209 210 210]\n",
            "[211 211 212 212 213 213 214 214 215 215 216 216 217 217 218 218 219 219 220 220 221 221 222 222 223 223 224 224]\n",
            "[225 225 226 226 227 227 228 228 229 229 230 230 231 231 232 232 233 233 234 234 235 235 236 236 237 237 238 238]\n",
            "[239 239 240 240 241 241 242 242 243 243 244 244 245 245 246 246 247 247 248 248 249 249 250 250 251 251 252 252]\n",
            "[253 253 254 254 255 255 256 256 257 257 258 258 259 259 260 260 261 261 262 262 263 263 264 264 265 265 266 266]\n",
            "[267 267 268 268 269 269 270 270 271 271 272 272 273 273 274 274 275 275 276 276 277 277 278 278 279 279 280 280]\n",
            "[281 281 282 282 283 283 284 284 285 285 286 286 287 287 288 288 289 289 290 290 291 291 292 292 293 293 294 294]\n",
            "[295 295 296 296 297 297 298 298 299 299 300 300 301 301 302 302 303 303 304 304 305 305 306 306 307 307 308 308]\n",
            "[309 309 310 310 311 311 312 312 313 313 314 314 315 315 316 316 317 317 318 318 319 319 320 320 321 321 322 322]\n",
            "[323 323 324 324 325 325 326 326 327 327 328 328 329 329 330 330 331 331 332 332 333 333 334 334 335 335 336 336]\n",
            "[337 337 338 338 339 339 340 340 341 341 342 342 343 343 344 344 345 345 346 346 347 347 348 348 349 349 350 350]\n",
            "[351 351 352 352 353 353 354 354 355 355 356 356 357 357 358 358 359 359 360 360 361 361 362 362 363 363 364 364]\n",
            "[365 365 366 366 367 367 368 368 369 369 370 370 371 371 372 372 373 373 374 374 375 375 376 376 377 377 378 378]\n",
            "[379 379 380 380 381 381 382 382 383 383 384 384 385 385 386 386 387 387 388 388 389 389 390 390 391 391 392 392]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Reshape into array to feed it into the model."
      ],
      "metadata": {
        "id": "c--QdeCfxVKj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "groups = groups.reshape((1,784))[0]"
      ],
      "metadata": {
        "id": "QKCWjwdwxYkU"
      },
      "execution_count": 16,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now to fit Group LASSO Multinomial Regression on this group schematic!"
      ],
      "metadata": {
        "id": "47UoTQUnwjfR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fit_summary(groups)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 750
        },
        "id": "Pthp-WLZ-XGa",
        "outputId": "36898fed-ebd8-4bc8-c0ee-5a3dad7f12f0"
      },
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Accuracy: 0.8151\n",
            "Classification Report:\n",
            "              precision    recall  f1-score   support\n",
            "\n",
            "           3       0.78      0.88      0.83      1426\n",
            "           5       0.81      0.76      0.78      1257\n",
            "           8       0.86      0.81      0.83      1373\n",
            "\n",
            "    accuracy                           0.82      4056\n",
            "   macro avg       0.82      0.81      0.81      4056\n",
            "weighted avg       0.82      0.82      0.81      4056\n",
            "\n",
            "There are 38 nonzero features after applying group LASSO.\n",
            "These are the pixels group LASSO found to be significant in classification of the digit to either 3, 5, or 8.\n",
            "\n",
            "Pixel List = ['pixel151', 'pixel152', 'pixel153', 'pixel154', 'pixel177', 'pixel178', 'pixel191', 'pixel192', 'pixel219', 'pixel220', 'pixel291', 'pixel292', 'pixel297', 'pixel298', 'pixel317', 'pixel318', 'pixel319', 'pixel320', 'pixel325', 'pixel326', 'pixel327', 'pixel328', 'pixel345', 'pixel346', 'pixel351', 'pixel352', 'pixel353', 'pixel354', 'pixel407', 'pixel408', 'pixel461', 'pixel462', 'pixel487', 'pixel488', 'pixel489', 'pixel490', 'pixel515', 'pixel516']\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGzCAYAAABpdMNsAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALnVJREFUeJzt3XtYVXW+x/HPRmErAptUFEgwFC9TqJkXxmwykxRtSi1vaaalNhpaJk6TNeWlJjpeytFj+kznpGWpk01qwzE7Sl4LbbSLXb0Q3kbBSwmKggi/84eP+7QDFHTDD/D9ep71PO61fvu3vnu53B/XWr+9lsMYYwQAQAXzsV0AAODaRAABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQChSsvMzFS/fv1Ur149ORwOzZ4923ZJAEqJAEK5+Oqrr+RwOLRr1y5J0quvvqobbrjB6+t58skn9dFHH2nSpElavHix4uPjvb4O/L/XX39dXbp0UcOGDeV0OhUVFaWHH35Y+/btK9I2MzNTDz/8sBo0aKDatWvrlltu0fLly4u0e//99zVw4EA1adJE/v7+atGihRITE3Xy5Mny/0CwysG94FAe/va3v2nSpEk6fvy4HA6HBgwYIB8fHy1btsyr6wkNDVVcXJzefvttr/aL4j322GM6c+aMWrVqpeuuu07p6el6/fXXVVBQoK+++krh4eGSpOzsbLVr106ZmZl64oknFBoaqnfffVebNm3SO++8o8GDB7v7rF+/vsLDw9WnTx9FRkbq66+/1oIFC9SkSRN9/vnnql27tq2Pi/JmgHIwYsQIEx8f737dqFEj88orr3h9PQ6HwyQkJHi93187ffp0ua+jqtq+fbuRZJKSktzzpk+fbiSZlJQU97yCggLToUMHExoaavLy8tzz169fX6TPN99800gyr7/+ernWDrs4BQev+fnnn3X8+HEdP35c27ZtU0xMjI4fP65vv/1Whw4dUrNmzXT8+HGdPn36sn39+OOP6t+/v+rWrSt/f3/99re/1f/8z/+4ly9atEgOh0PGGM2bN08Oh0MOh+OSfZ44cUJDhw5VUFCQgoODNWzYMPepwkWLFrnbDR8+XAEBAUpLS1OvXr0UGBioIUOGSJJycnKUmJioiIgIOZ1OtWjRQjNnzpT5xYmEffv2FenzIofDoSlTprhfT5kyRQ6HQz/88IMGDBigoKAg1atXT0888YRyc3M93rt27VrddtttCg4OVkBAgFq0aKFnnnnmstuyvF08tfrLU2abN29WSEiI7rzzTvc8Hx8fDRgwQBkZGdq4caN7/h133FGkz759+0qSvv/++3KpGZVDTdsFoPpo27at9u/f7379zTffaObMme7X99xzjyRp2LBhxX45X5SZmalbb71VZ86c0eOPP6569erpzTff1L333qv33ntPffv21e23367Fixdr6NChuuuuu/TQQw9dsrbCwkLdc889+uyzzzRmzBi1bNlSq1at0rBhw4ptf/78efXo0UO33XabZs6cKX9/fxljdO+992r9+vUaMWKEbr75Zn300Uf64x//qH//+9969dVXy7C1PA0YMEA33HCDkpKStHXrVs2ZM0c///yz3nrrLUnSt99+q9///vdq3bq1pk2bJqfTqb179+qTTz65bN9ZWVnKz8+/bLtatWopICCgVPWeOHFCBQUFOnDggKZNmyZJ6tatm3t5Xl5esafO/P39JUk7duzQXXfdVWL/GRkZki6cnkM1ZvkIDNXIli1bzNq1a81zzz1natasaT788EOzdu1a07NnT9O+fXuzdu1as3btWvPtt99esp/x48cbSWbz5s3ueadOnTJRUVHmhhtuMAUFBe75kkp1Cu4f//iHkWRmz57tnldQUGDuvPNOI8ksXLjQPX/YsGFGknn66ac9+li5cqWRZF588UWP+f369TMOh8Ps3bvXGGNMenp6kT5/We/kyZPdrydPnmwkmXvvvdej3WOPPWYkma+++soYY8yrr75qJJljx45d9rP+WpcuXYyky07Dhg0rdZ9Op9P9vnr16pk5c+Z4LB83bpzx8fEx+/bt85g/aNAgI8mMHTv2kv2PGDHC1KhRw+zevbvUNaHq4QgIXtO5c2dJ0urVq9WhQwf3iLTx48erf//+iouLK1U/q1evVseOHXXbbbe55wUEBOjRRx/VpEmT9N133ykmJqZMta1Zs0a+vr4aNWqUe56Pj48SEhL08ccfF/ueMWPGFKmrRo0aevzxxz3mJyYm6r333tOHH36osWPHlqmuixISEjxejxs3Tq+99ppWr16t1q1bKzg4WJK0atUqPfzww/LxKf3Z81mzZunnn3++bLuLAwhK48MPP1Rubq6+//57vf3228rJyfFYPnLkSC1YsEADBgzQq6++qoYNG+rdd9/VihUrJElnz54tse8lS5bov//7v/XUU0+pWbNmpa4JVQ8BBK/45WmelJQU3XnnnTp+/Lh++uknffvtt3rxxRd1/Phx+fr6yuVyXbKv/fv3KzY2tsj83/zmN+7lZQ2g/fv3KywszH0K6KLo6Ohi29esWVONGjUq0kd4eLgCAwNLrOtK/fqLtmnTpvLx8XEPbx44cKD+67/+SyNHjtTTTz+tbt266b777lO/fv0uG0bt2rW74rpK0rVrV0lSz5491bt3b8XExCggIMAdwK1bt9aSJUs0evRo939MQkNDNXv2bI0ZM6bEU32bN2/WiBEj1KNHD/3lL3/xet2oXBiEAK/o3bu3QkJCFBISop07d2r27NkKCQlRixYtJF24qBwSEqLevXtbrrR0nE5nmY4yfqmkwRAFBQVX3Eft2rW1adMmrVu3TkOHDtXOnTs1cOBA3XXXXZft96efflJGRsZlp6ysrFLX90tNmzZV27Zt9c4773jM79evnw4fPqzPPvtMqamp2r9/v5o0aSJJat68eZF+vvrqK917772KiYnRe++9p5o1+f9xdcffMLzi4mme1NRUTZ06VcnJyapZs6bmzp2rf//733r55ZclSdddd91l+2rcuLH7B6y/9MMPP7iXl1Xjxo21fv16nTlzxuMoaO/evWXqY926dTp16pTHUdCv67r4GX/9Q8pLHSHt2bNHUVFRHnUVFhZ6/HjXx8dH3bp1U7du3fTKK6/opZde0rPPPqv169df8vTmfffd5zHqrCSXGxxyKWfPnlVeXl6R+X5+furQoYP79bp16ySpSL1paWmKj49XgwYNtHr16lIPhkDVxhEQvKJdu3aKi4vT+fPnFRMTo/j4eMXFxSkzM1NxcXHuqTSng3r16uX+X/NFOTk5+tvf/qYbbrhBN954Y5nr69Gjh/Lz8/X666+75xUWFmrevHml7qNXr14qKCjQf/7nf3rMf/XVV+VwONSzZ09JUlBQkOrXr69NmzZ5tHvttddK7PvXdcydO1eS3H3+9NNPRd5z8803S1KxX/y/NGvWLK1du/ay01NPPXXJfs6fP1/staTPPvtMX3/9tdq3b3/J9+/Zs0cLFizQ73//e48joIyMDHXv3l0+Pj766KOPFBIScsl+UH1wBASv+uSTT3TrrbdKknJzc/XFF1+U+bcqTz/9tJYuXaqePXvq8ccfV926dfXmm28qPT1d//jHP67o1FifPn3UsWNHJSYmau/evWrZsqU++OAD9xf75X5DJF0YRt61a1c9++yz2rdvn9q0aaP//d//1apVqzR+/Hg1bdrU3XbkyJF6+eWXNXLkSLVv316bNm3S7t27S+w7PT1d9957r+Lj45Wamqq3335bgwcPVps2bSRJ06ZN06ZNm3T33XercePGOnr0qF577TU1atTIY7BGcbx1Dej06dOKiIjQwIEDddNNN6lOnTr6+uuvtXDhQrlcLj333HMe7W+88Ub1799fkZGRSk9P1/z581W3bl0tWLDAo118fLx+/PFHPfXUU9qyZYu2bNniXtawYcNLDtdGFWd7GB6qj/Pnz5uAgACzePFiY8yFYdmSzNGjR8vcV1pamunXr58JDg42tWrVMh07djTJyclF2qmUw7CNMebYsWNm8ODBJjAw0LhcLjN8+HDzySefGElm2bJl7nbDhg0zderUKbaPU6dOmSeffNKEh4cbX19f06xZMzNjxgxTWFjo0e7MmTNmxIgRxuVymcDAQDNgwABz9OjREodhf/fdd6Zfv34mMDDQXHfddWbs2LHm7Nmz7nYpKSmmd+/eJjw83Pj5+Znw8HDzwAMPVOgw5by8PPPEE0+Y1q1bm6CgIOPr62saN25sRowYYdLT04u0HzRokImIiHDXO3r0aJOZmVmknS4xNLxLly7l/8FgDfeCwzVt5cqV6tu3r7Zs2eIerVWRpkyZoqlTp+rYsWP86BLXHK4B4Zrx69+eFBQUaO7cuQoKCtItt9xiqSrg2sU1IFwzxo0bp7Nnz6pTp07Ky8vT+++/r08//VQvvfQSd1wGLCCAcM248847NWvWLCUnJys3N1fR0dGaO3fuFd+9AMDV4RoQAMAKrgEBAKwggAAAVlS6a0CFhYU6fPiwAgMDS/XjQABA5WKM0alTpxQeHn7JH45XugA6fPiwIiIibJcBALhKBw8eLHJX+V+qdAF08SaPBw8eVFBQkOVqAABllZ2drYiIiCKPLvm1cgugefPmacaMGcrIyFCbNm00d+5cdezY8bLvu3jaLSgoiAACgCrscpdRymUQwt///ndNmDBBkydP1ueff642bdqoR48eOnr0aHmsDgBQBZVLAL3yyisaNWqUHn74Yd14441asGCB/P399cYbb5TH6gAAVZDXA+jcuXPasWOHxwOnfHx8FBcX5/F8l4vy8vKUnZ3tMQEAqj+vXwM6fvy4CgoK1LBhQ4/5DRs2dD858peSkpI0depUb5cBoAorKChQfn6+7TJQgho1aqhmzZpX/VMZ66PgJk2apAkTJrhfXxw9AeDadPr0aR06dEjcJaxy8/f3V1hYmPz8/K64D68HUP369VWjRg1lZmZ6zM/MzFRoaGiR9k6nU06n09tlAKiCCgoKdOjQIfn7+yskJIQfo1dCxhidO3dOx44dU3p6upo1a3ZFTymWyiGA/Pz81K5dO6WkpKhPnz6SLtzdICUlhbsOA7ik/Px8GWMUEhLCIzIqsdq1a8vX11f79+/XuXPnVKtWrSvqp1xOwU2YMEHDhg1T+/bt1bFjR82ePVs5OTl6+OGHy2N1AKoZjnwqvys96vmlcgmggQMH6tixY3r++eeVkZGhm2++WWvWrCkyMAEAcO0qt0EIY8eO5ZQbAK84cOCAjh8/XmHrq1+/viIjIytsfdcq66PgAOBSDhw4oBYtWig3N7fC1lmrVi3t2rXLayE0fPhwnTx5UitXrvRKf4sWLdL48eN18uRJr/RnC88DAlCpHT9+vELDR5Jyc3PLdMQ1fPhwORwOORwO+fn5KTo6WtOmTdP58+clSX/961+1aNGicqq2ZMnJyerSpYsCAwPl7++vDh06WKmjJAQQAHhBfHy8jhw5oj179igxMVFTpkzRjBkzJEkul0vBwcEVWs/cuXPVu3dvde7cWdu2bdPOnTs1aNAgjR49WhMnTqzQWkpCAAGAFzidToWGhqpx48YaM2aM4uLi9MEHH0i6cIR08Wcpx44dU2hoqF566SX3ez/99FP5+fkpJSVF0oVblE2cOFHXX3+96tSpo9jYWG3YsKHUtRw8eFCJiYkaP368XnrpJd14442Kjo5WYmKiZsyYoVmzZmnbtm2SpA0bNsjhcCglJUXt27eXv7+/br31Vu3atcs7G+YSCCAAKAe1a9fWuXPniswPCQnRG2+8oSlTpmj79u06deqUhg4dqrFjx6pbt26SLgziSk1N1bJly7Rz5071799f8fHx2rNnT6nW/d577yk/P7/YI50//OEPCggI0NKlSz3mP/vss5o1a5a2b9+umjVr6pFHHrmCT102DEIAAC8yxiglJUUfffSRxo0bV2ybXr16adSoURoyZIjat2+vOnXqKCkpSdKFQRcLFy7UgQMHFB4eLkmaOHGi1qxZo4ULF3ocOZVk9+7dcrlcCgsLK7LMz89PTZo00e7duz3m/+Uvf1GXLl0kSU8//bTuvvtu5ebmXvGPTEuDAAIAL0hOTlZAQIDy8/NVWFiowYMHa8qUKSW2nzlzpmJiYrR8+XLt2LHDfUuyr7/+WgUFBWrevLlH+7y8PNWrV6/c6m/durX7zxeD6+jRo+U6HJ0AAgAv6Nq1q+bPny8/Pz+Fh4erZs1Lf72mpaXp8OHDKiws1L59+9SqVStJF27GWqNGDe3YsUM1atTweE9AQECpamnevLmysrJ0+PBh91HURefOnVNaWpq6du3qMd/X19f954t3oigsLCzV+q4U14AAwAvq1Kmj6OhoRUZGXjZ8zp07pwcffFADBw7UCy+8oJEjR7qfGN22bVsVFBTo6NGjio6O9piKu6Fzce6//375+vpq1qxZRZYtWLBAOTk5euCBB8r+Ib2MIyAAqGDPPvussrKyNGfOHAUEBGj16tV65JFHlJycrObNm2vIkCF66KGHNGvWLLVt21bHjh1TSkqKWrdurbvvvvuy/UdGRmr69OlKTExUrVq1NHToUPn6+mrVqlV65plnlJiYqNjY2Ar4pJfGERCASq1+/frleiG8OLVq1VL9+vXLpe8NGzZo9uzZWrx4sYKCguTj46PFixdr8+bNmj9/viRp4cKFeuihh5SYmKgWLVqoT58++te//lWm6zHjx4/XihUrtHnzZrVv314xMTFasmSJ5s+fr5kzZ5bLZysrh6lkT33Kzs6Wy+VSVlaWgoKCbJcDoALl5uYqPT1dUVFRHqHDveAqn5L+rqTSf49zCg5ApRcZGUkgVEOcggMAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAKAcvbLJ6J6w6JFi7z6iO9PPvlErVq1kq+vr7vO4uZ5GwEEAFdp+PDhcjgccjgc8vPzU3R0tKZNm6bz589Lkv76179q0aJFFV5XcnKyunTposDAQPn7+6tDhw7F1jFhwgTdfPPNSk9Pdy8vbp63EUAA4AXx8fE6cuSI9uzZo8TERE2ZMkUzZsyQJLlcLq8esZTG3Llz1bt3b3Xu3Fnbtm3Tzp07NWjQII0ePbrIo7rT0tJ05513qlGjRu46i5vnbQQQAHiB0+lUaGioGjdurDFjxiguLk4ffPCBJM9TcMeOHVNoaKjHo7U//fRT+fn5KSUlRdKFp59OnDhR119/verUqaPY2Fht2LCh1LUcPHhQiYmJGj9+vF566SXdeOONio6OVmJiombMmKFZs2Zp27Zt2rdvnxwOh06cOKFHHnlEDodDixYtKnZeeSCAAKAc1K5dW+fOnSsyPyQkRG+88YamTJmi7du369SpUxo6dKjGjh2rbt26SZLGjh2r1NRULVu2TDt37lT//v0VHx+vPXv2lGrd7733nvLz84sc6UjSH/7wBwUEBGjp0qWKiIjQkSNHFBQUpNmzZ+vIkSPq379/kXkDBw68uo1RAu6GDQBeZIxRSkqKPvroI40bN67YNr169dKoUaM0ZMgQtW/fXnXq1FFSUpKkC4+eWLhwoQ4cOOB+nPbEiRO1Zs0aLVy40OPIqSS7d++Wy+VSWFhYkWV+fn5q0qSJdu/erRo1aig0NFQOh0Mul8v9xNU6deoUmVceCCAA8ILk5GQFBAQoPz9fhYWFGjx4sKZMmVJi+5kzZyomJkbLly/Xjh075HQ6JUlff/21CgoK1Lx5c4/2eXl5qlevXnl+hApHAAGAF3Tt2lXz58+Xn5+fwsPDVbPmpb9e09LSdPjwYRUWFmrfvn1q1aqVJOn06dOqUaOGduzYoRo1ani8JyAgoFS1NG/eXFlZWTp8+LD7KOqic+fOKS0tTV27di3DpysfXAMCAC+oU6eOoqOjFRkZednwOXfunB588EENHDhQL7zwgkaOHKmjR49Kktq2bauCggIdPXpU0dHRHlNpT4fdf//98vX11axZs4osW7BggXJycvTAAw+U/UN6GUdAAFDBnn32WWVlZWnOnDkKCAjQ6tWr9cgjjyg5OVnNmzfXkCFD9NBDD2nWrFlq27atjh07ppSUFLVu3Vp33333ZfuPjIzU9OnTlZiYqFq1amno0KHy9fXVqlWr9MwzzygxMVGxsbEV8EkvjSMgAKhAGzZs0OzZs7V48WIFBQXJx8dHixcv1ubNmzV//nxJ0sKFC/XQQw8pMTFRLVq0UJ8+ffSvf/2rTI8lHz9+vFasWKHNmzerffv2iomJ0ZIlSzR//nzNnDmzvD5emTiMMcZ2Eb+UnZ0tl8ulrKwsBQUF2S4HQAXKzc1Venq6oqKiVKtWLdvl4BIu9XdV2u9xjoAAAFYQQAAAKwggAIAVBBAAwAoCCEClU8nGRqEY3vg7IoAAVBoXf/lf3E08UbmcOXNGkuTr63vFffBDVACVRs2aNeXv769jx47J19dXPj78H7myMcbozJkzOnr0qIKDg4vcLqgsCCAAlYbD4VBYWJjS09O1f/9+2+XgEoKDg6/6TtkEEIBKxc/PT82aNeM0XCXm6+t7VUc+FxFAACodHx8f7oRwDeAEKwDACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACu6EgErP4XDYLsHreNwAwBEQAMASAggAYIXXA2jKlClyOBweU8uWLb29GgBAFVcu14BuuukmrVu37v9XUpNLTQAAT+WSDDVr1rzqBxUBAKq3crkGtGfPHoWHh6tJkyYaMmSIDhw4UGLbvLw8ZWdne0wAgOrP6wEUGxurRYsWac2aNZo/f77S09P1u9/9TqdOnSq2fVJSklwul3uKiIjwdkkAgErIYcr5BwknT55U48aN9corr2jEiBFFlufl5SkvL8/9Ojs7WxEREcrKylJQUFB5loYqgt8BAVVLdna2XC7XZb/Hy310QHBwsJo3b669e/cWu9zpdMrpdJZ3GQCASqbcfwd0+vRppaWlKSwsrLxXBQCoQrweQBMnTtTGjRu1b98+ffrpp+rbt69q1KihBx54wNurAgBUYV4/BXfo0CE98MADOnHihEJCQnTbbbdp69atCgkJ8faqAABVmNcDaNmyZd7uEpVURQ0O4II9bGEfL1/cCw4AYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCj3B9Kh+rpWb6AIwDs4AgIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAV3A0bAErAHd/LF0dAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAAratouALgch8NRYesyxlTYuioC2w6VGUdAAAArCCAAgBVlDqBNmzbpnnvuUXh4uBwOh1auXOmx3Bij559/XmFhYapdu7bi4uK0Z88eb9ULAKgmyhxAOTk5atOmjebNm1fs8unTp2vOnDlasGCBtm3bpjp16qhHjx7Kzc296mIBANVHmQch9OzZUz179ix2mTFGs2fP1p///Gf17t1bkvTWW2+pYcOGWrlypQYNGnR11QIAqg2vXgNKT09XRkaG4uLi3PNcLpdiY2OVmppa7Hvy8vKUnZ3tMQEAqj+vBlBGRoYkqWHDhh7zGzZs6F72a0lJSXK5XO4pIiLCmyUBACop66PgJk2apKysLPd08OBB2yUBACqAVwMoNDRUkpSZmekxPzMz073s15xOp4KCgjwmAED159UAioqKUmhoqFJSUtzzsrOztW3bNnXq1MmbqwIAVHFlHgV3+vRp7d271/06PT1dX375perWravIyEiNHz9eL774opo1a6aoqCg999xzCg8PV58+fbxZNwCgiitzAG3fvl1du3Z1v54wYYIkadiwYVq0aJGeeuop5eTk6NFHH9XJkyd12223ac2aNapVq5b3qgYAVHkOU8nuIJidnS2Xy6WsrCyuB1VyFXmjy4pSyf45XDX+jmBDab/HrY+CAwBcmwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCizI9jACpadbz7MXepBjgCAgBYQgABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAAruBkprhg3n7xybDuAIyAAgCUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIKbkQIWOByOClkPNz1FZcYREADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYwc1IAQu4SSjAERAAwBICCABgRZkDaNOmTbrnnnsUHh4uh8OhlStXeiwfPny4HA6HxxQfH++tegEA1USZAygnJ0dt2rTRvHnzSmwTHx+vI0eOuKelS5deVZEAgOqnzIMQevbsqZ49e16yjdPpVGho6BUXBQCo/srlGtCGDRvUoEEDtWjRQmPGjNGJEydKbJuXl6fs7GyPCQBQ/Xk9gOLj4/XWW28pJSVF//Ef/6GNGzeqZ8+eKigoKLZ9UlKSXC6Xe4qIiPB2SQCASshhruIHCQ6HQytWrFCfPn1KbPPjjz+qadOmWrdunbp161ZkeV5envLy8tyvs7OzFRERoaysLAUFBV1paQAAS7Kzs+VyuS77PV7uw7CbNGmi+vXra+/evcUudzqdCgoK8pgAANVfuQfQoUOHdOLECYWFhZX3qgAAVUiZR8GdPn3a42gmPT1dX375perWrau6detq6tSpuv/++xUaGqq0tDQ99dRTio6OVo8ePbxaOACgaitzAG3fvl1du3Z1v54wYYIkadiwYZo/f7527typN998UydPnlR4eLi6d++uF154QU6n03tVAwCqvKsahFAeSnvxCigPDoejQtZTyf7ZAV5VaQYhAABQHAIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwo8+MYgIpWUXeoBlCxOAICAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACu4GSnwC8YY2yUA1wyOgAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACm5GikqPG4QC1RNHQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsKFMAJSUlqUOHDgoMDFSDBg3Up08f7dq1y6NNbm6uEhISVK9ePQUEBOj+++9XZmamV4sGAFR9ZQqgjRs3KiEhQVu3btXatWuVn5+v7t27Kycnx93mySef1D//+U8tX75cGzdu1OHDh3Xfffd5vXAAQNXmMFfxuMljx46pQYMG2rhxo26//XZlZWUpJCRES5YsUb9+/SRJP/zwg37zm98oNTVVv/3tby/bZ3Z2tlwul7KyshQUFHSlpQEALCnt9/hVXQPKysqSJNWtW1eStGPHDuXn5ysuLs7dpmXLloqMjFRqamqxfeTl5Sk7O9tjAgBUf1ccQIWFhRo/frw6d+6smJgYSVJGRob8/PwUHBzs0bZhw4bKyMgotp+kpCS5XC73FBERcaUlAQCqkCsOoISEBH3zzTdatmzZVRUwadIkZWVluaeDBw9eVX8AgKqh5pW8aezYsUpOTtamTZvUqFEj9/zQ0FCdO3dOJ0+e9DgKyszMVGhoaLF9OZ1OOZ3OKykDAFCFlekIyBijsWPHasWKFfr4448VFRXlsbxdu3by9fVVSkqKe96uXbt04MABderUyTsVAwCqhTIdASUkJGjJkiVatWqVAgMD3dd1XC6XateuLZfLpREjRmjChAmqW7eugoKCNG7cOHXq1KlUI+AAANeOMg3Ddjgcxc5fuHChhg8fLunCD1ETExO1dOlS5eXlqUePHnrttddKPAX3awzDBoCqrbTf41f1O6DyQAABQNVWIb8DAgDgShFAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYUaYASkpKUocOHRQYGKgGDRqoT58+2rVrl0ebO+64Qw6Hw2MaPXq0V4sGAFR9ZQqgjRs3KiEhQVu3btXatWuVn5+v7t27Kycnx6PdqFGjdOTIEfc0ffp0rxYNAKj6apal8Zo1azxeL1q0SA0aNNCOHTt0++23u+f7+/srNDTUOxUCAKqlq7oGlJWVJUmqW7eux/x33nlH9evXV0xMjCZNmqQzZ86U2EdeXp6ys7M9JgBA9VemI6BfKiws1Pjx49W5c2fFxMS45w8ePFiNGzdWeHi4du7cqT/96U/atWuX3n///WL7SUpK0tSpU6+0DABAFeUwxpgreeOYMWP04YcfasuWLWrUqFGJ7T7++GN169ZNe/fuVdOmTYssz8vLU15envt1dna2IiIilJWVpaCgoCspDQBgUXZ2tlwu12W/x6/oCGjs2LFKTk7Wpk2bLhk+khQbGytJJQaQ0+mU0+m8kjIAAFVYmQLIGKNx48ZpxYoV2rBhg6Kioi77ni+//FKSFBYWdkUFAgCqpzIFUEJCgpYsWaJVq1YpMDBQGRkZkiSXy6XatWsrLS1NS5YsUa9evVSvXj3t3LlTTz75pG6//Xa1bt26XD4AAKBqKtM1IIfDUez8hQsXavjw4Tp48KAefPBBffPNN8rJyVFERIT69u2rP//5z6W+nlPac4cAgMqpXK4BXS6rIiIitHHjxrJ0CQC4RnEvOACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFTVtF/BrxhhJUnZ2tuVKAABX4uL398Xv85JUugA6deqUJCkiIsJyJQCAq3Hq1Cm5XK4SlzvM5SKqghUWFurw4cMKDAyUw+HwWJadna2IiAgdPHhQQUFBliq0j+1wAdvhArbDBWyHCyrDdjDG6NSpUwoPD5ePT8lXeirdEZCPj48aNWp0yTZBQUHX9A52EdvhArbDBWyHC9gOF9jeDpc68rmIQQgAACsIIACAFVUqgJxOpyZPniyn02m7FKvYDhewHS5gO1zAdrigKm2HSjcIAQBwbahSR0AAgOqDAAIAWEEAAQCsIIAAAFYQQAAAK6pMAM2bN0833HCDatWqpdjYWH322We2S6pwU6ZMkcPh8Jhatmxpu6xyt2nTJt1zzz0KDw+Xw+HQypUrPZYbY/T8888rLCxMtWvXVlxcnPbs2WOn2HJ0ue0wfPjwIvtHfHy8nWLLSVJSkjp06KDAwEA1aNBAffr00a5duzza5ObmKiEhQfXq1VNAQIDuv/9+ZWZmWqq4fJRmO9xxxx1F9ofRo0dbqrh4VSKA/v73v2vChAmaPHmyPv/8c7Vp00Y9evTQ0aNHbZdW4W666SYdOXLEPW3ZssV2SeUuJydHbdq00bx584pdPn36dM2ZM0cLFizQtm3bVKdOHfXo0UO5ubkVXGn5utx2kKT4+HiP/WPp0qUVWGH527hxoxISErR161atXbtW+fn56t69u3JyctxtnnzySf3zn//U8uXLtXHjRh0+fFj33Xefxaq9rzTbQZJGjRrlsT9Mnz7dUsUlMFVAx44dTUJCgvt1QUGBCQ8PN0lJSRarqniTJ082bdq0sV2GVZLMihUr3K8LCwtNaGiomTFjhnveyZMnjdPpNEuXLrVQYcX49XYwxphhw4aZ3r17W6nHlqNHjxpJZuPGjcaYC3/3vr6+Zvny5e4233//vZFkUlNTbZVZ7n69HYwxpkuXLuaJJ56wV1QpVPojoHPnzmnHjh2Ki4tzz/Px8VFcXJxSU1MtVmbHnj17FB4eriZNmmjIkCE6cOCA7ZKsSk9PV0ZGhsf+4XK5FBsbe03uHxs2bFCDBg3UokULjRkzRidOnLBdUrnKysqSJNWtW1eStGPHDuXn53vsDy1btlRkZGS13h9+vR0ueuedd1S/fn3FxMRo0qRJOnPmjI3ySlTp7ob9a8ePH1dBQYEaNmzoMb9hw4b64YcfLFVlR2xsrBYtWqQWLVroyJEjmjp1qn73u9/pm2++UWBgoO3yrMjIyJCkYvePi8uuFfHx8brvvvsUFRWltLQ0PfPMM+rZs6dSU1NVo0YN2+V5XWFhocaPH6/OnTsrJiZG0oX9wc/PT8HBwR5tq/P+UNx2kKTBgwercePGCg8P186dO/WnP/1Ju3bt0vvvv2+xWk+VPoDw/3r27On+c+vWrRUbG6vGjRvr3Xff1YgRIyxWhspg0KBB7j+3atVKrVu3VtOmTbVhwwZ169bNYmXlIyEhQd988801cR30UkraDo8++qj7z61atVJYWJi6deumtLQ0NW3atKLLLFalPwVXv3591ahRo8golszMTIWGhlqqqnIIDg5W8+bNtXfvXtulWHNxH2D/KKpJkyaqX79+tdw/xo4dq+TkZK1fv97j+WGhoaE6d+6cTp486dG+uu4PJW2H4sTGxkpSpdofKn0A+fn5qV27dkpJSXHPKywsVEpKijp16mSxMvtOnz6ttLQ0hYWF2S7FmqioKIWGhnrsH9nZ2dq2bds1v38cOnRIJ06cqFb7hzFGY8eO1YoVK/Txxx8rKirKY3m7du3k6+vrsT/s2rVLBw4cqFb7w+W2Q3G+/PJLSapc+4PtURClsWzZMuN0Os2iRYvMd999Zx599FETHBxsMjIybJdWoRITE82GDRtMenq6+eSTT0xcXJypX7++OXr0qO3SytWpU6fMF198Yb744gsjybzyyivmiy++MPv37zfGGPPyyy+b4OBgs2rVKrNz507Tu3dvExUVZc6ePWu5cu+61HY4deqUmThxoklNTTXp6elm3bp15pZbbjHNmjUzubm5tkv3mjFjxhiXy2U2bNhgjhw54p7OnDnjbjN69GgTGRlpPv74Y7N9+3bTqVMn06lTJ4tVe9/ltsPevXvNtGnTzPbt2016erpZtWqVadKkibn99tstV+6pSgSQMcbMnTvXREZGGj8/P9OxY0ezdetW2yVVuIEDB5qwsDDj5+dnrr/+ejNw4ECzd+9e22WVu/Xr1xtJRaZhw4YZYy4MxX7uuedMw4YNjdPpNN26dTO7du2yW3Q5uNR2OHPmjOnevbsJCQkxvr6+pnHjxmbUqFHV7j9pxX1+SWbhwoXuNmfPnjWPPfaYue6664y/v7/p27evOXLkiL2iy8HltsOBAwfM7bffburWrWucTqeJjo42f/zjH01WVpbdwn+F5wEBAKyo9NeAAADVEwEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWPF/cyBXF+0ihKsAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "This looks a bit nicer then the single pixel grouping schema and the accuracy has improved.  I bet the accuracy would improve if we implemented 4x4 chunks of pixels as groups!  Lets try this."
      ],
      "metadata": {
        "id": "tJtbFAmIzHwq"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 2x2 Chunk Grouping"
      ],
      "metadata": {
        "id": "e1nGHv8kzXYV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Create a 28x28 matrix filled with zeros\n",
        "groups = np.zeros((28, 28))\n",
        "\n",
        "# Calculate how many 2x2 chunks we'll have in total\n",
        "chunks_per_row = 28 // 2\n",
        "chunks_per_col = 28 // 2\n",
        "total_chunks = chunks_per_row * chunks_per_col  # This should equal 196 (28×28/4)\n",
        "\n",
        "\n",
        "counter = 1 # Fill the matrix with 2x2 chunks counting up from 1\n",
        "for i in range(0, 28, 2):\n",
        "    for j in range(0, 28, 2):\n",
        "        groups[i:i+2, j:j+2] = counter\n",
        "        counter += 1"
      ],
      "metadata": {
        "id": "qjEscX8J20Si"
      },
      "execution_count": 18,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "groups"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "zvIfCJjg3kdZ",
        "outputId": "5a75f0db-e376-4c36-914d-d59bf4e771fc"
      },
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[  1.,   1.,   2.,   2.,   3.,   3.,   4.,   4.,   5.,   5.,   6.,   6.,   7.,   7.,   8.,   8.,   9.,   9.,  10.,  10.,  11.,  11.,  12.,  12.,  13.,  13.,  14.,  14.],\n",
              "       [  1.,   1.,   2.,   2.,   3.,   3.,   4.,   4.,   5.,   5.,   6.,   6.,   7.,   7.,   8.,   8.,   9.,   9.,  10.,  10.,  11.,  11.,  12.,  12.,  13.,  13.,  14.,  14.],\n",
              "       [ 15.,  15.,  16.,  16.,  17.,  17.,  18.,  18.,  19.,  19.,  20.,  20.,  21.,  21.,  22.,  22.,  23.,  23.,  24.,  24.,  25.,  25.,  26.,  26.,  27.,  27.,  28.,  28.],\n",
              "       [ 15.,  15.,  16.,  16.,  17.,  17.,  18.,  18.,  19.,  19.,  20.,  20.,  21.,  21.,  22.,  22.,  23.,  23.,  24.,  24.,  25.,  25.,  26.,  26.,  27.,  27.,  28.,  28.],\n",
              "       [ 29.,  29.,  30.,  30.,  31.,  31.,  32.,  32.,  33.,  33.,  34.,  34.,  35.,  35.,  36.,  36.,  37.,  37.,  38.,  38.,  39.,  39.,  40.,  40.,  41.,  41.,  42.,  42.],\n",
              "       [ 29.,  29.,  30.,  30.,  31.,  31.,  32.,  32.,  33.,  33.,  34.,  34.,  35.,  35.,  36.,  36.,  37.,  37.,  38.,  38.,  39.,  39.,  40.,  40.,  41.,  41.,  42.,  42.],\n",
              "       [ 43.,  43.,  44.,  44.,  45.,  45.,  46.,  46.,  47.,  47.,  48.,  48.,  49.,  49.,  50.,  50.,  51.,  51.,  52.,  52.,  53.,  53.,  54.,  54.,  55.,  55.,  56.,  56.],\n",
              "       [ 43.,  43.,  44.,  44.,  45.,  45.,  46.,  46.,  47.,  47.,  48.,  48.,  49.,  49.,  50.,  50.,  51.,  51.,  52.,  52.,  53.,  53.,  54.,  54.,  55.,  55.,  56.,  56.],\n",
              "       [ 57.,  57.,  58.,  58.,  59.,  59.,  60.,  60.,  61.,  61.,  62.,  62.,  63.,  63.,  64.,  64.,  65.,  65.,  66.,  66.,  67.,  67.,  68.,  68.,  69.,  69.,  70.,  70.],\n",
              "       [ 57.,  57.,  58.,  58.,  59.,  59.,  60.,  60.,  61.,  61.,  62.,  62.,  63.,  63.,  64.,  64.,  65.,  65.,  66.,  66.,  67.,  67.,  68.,  68.,  69.,  69.,  70.,  70.],\n",
              "       [ 71.,  71.,  72.,  72.,  73.,  73.,  74.,  74.,  75.,  75.,  76.,  76.,  77.,  77.,  78.,  78.,  79.,  79.,  80.,  80.,  81.,  81.,  82.,  82.,  83.,  83.,  84.,  84.],\n",
              "       [ 71.,  71.,  72.,  72.,  73.,  73.,  74.,  74.,  75.,  75.,  76.,  76.,  77.,  77.,  78.,  78.,  79.,  79.,  80.,  80.,  81.,  81.,  82.,  82.,  83.,  83.,  84.,  84.],\n",
              "       [ 85.,  85.,  86.,  86.,  87.,  87.,  88.,  88.,  89.,  89.,  90.,  90.,  91.,  91.,  92.,  92.,  93.,  93.,  94.,  94.,  95.,  95.,  96.,  96.,  97.,  97.,  98.,  98.],\n",
              "       [ 85.,  85.,  86.,  86.,  87.,  87.,  88.,  88.,  89.,  89.,  90.,  90.,  91.,  91.,  92.,  92.,  93.,  93.,  94.,  94.,  95.,  95.,  96.,  96.,  97.,  97.,  98.,  98.],\n",
              "       [ 99.,  99., 100., 100., 101., 101., 102., 102., 103., 103., 104., 104., 105., 105., 106., 106., 107., 107., 108., 108., 109., 109., 110., 110., 111., 111., 112., 112.],\n",
              "       [ 99.,  99., 100., 100., 101., 101., 102., 102., 103., 103., 104., 104., 105., 105., 106., 106., 107., 107., 108., 108., 109., 109., 110., 110., 111., 111., 112., 112.],\n",
              "       [113., 113., 114., 114., 115., 115., 116., 116., 117., 117., 118., 118., 119., 119., 120., 120., 121., 121., 122., 122., 123., 123., 124., 124., 125., 125., 126., 126.],\n",
              "       [113., 113., 114., 114., 115., 115., 116., 116., 117., 117., 118., 118., 119., 119., 120., 120., 121., 121., 122., 122., 123., 123., 124., 124., 125., 125., 126., 126.],\n",
              "       [127., 127., 128., 128., 129., 129., 130., 130., 131., 131., 132., 132., 133., 133., 134., 134., 135., 135., 136., 136., 137., 137., 138., 138., 139., 139., 140., 140.],\n",
              "       [127., 127., 128., 128., 129., 129., 130., 130., 131., 131., 132., 132., 133., 133., 134., 134., 135., 135., 136., 136., 137., 137., 138., 138., 139., 139., 140., 140.],\n",
              "       [141., 141., 142., 142., 143., 143., 144., 144., 145., 145., 146., 146., 147., 147., 148., 148., 149., 149., 150., 150., 151., 151., 152., 152., 153., 153., 154., 154.],\n",
              "       [141., 141., 142., 142., 143., 143., 144., 144., 145., 145., 146., 146., 147., 147., 148., 148., 149., 149., 150., 150., 151., 151., 152., 152., 153., 153., 154., 154.],\n",
              "       [155., 155., 156., 156., 157., 157., 158., 158., 159., 159., 160., 160., 161., 161., 162., 162., 163., 163., 164., 164., 165., 165., 166., 166., 167., 167., 168., 168.],\n",
              "       [155., 155., 156., 156., 157., 157., 158., 158., 159., 159., 160., 160., 161., 161., 162., 162., 163., 163., 164., 164., 165., 165., 166., 166., 167., 167., 168., 168.],\n",
              "       [169., 169., 170., 170., 171., 171., 172., 172., 173., 173., 174., 174., 175., 175., 176., 176., 177., 177., 178., 178., 179., 179., 180., 180., 181., 181., 182., 182.],\n",
              "       [169., 169., 170., 170., 171., 171., 172., 172., 173., 173., 174., 174., 175., 175., 176., 176., 177., 177., 178., 178., 179., 179., 180., 180., 181., 181., 182., 182.],\n",
              "       [183., 183., 184., 184., 185., 185., 186., 186., 187., 187., 188., 188., 189., 189., 190., 190., 191., 191., 192., 192., 193., 193., 194., 194., 195., 195., 196., 196.],\n",
              "       [183., 183., 184., 184., 185., 185., 186., 186., 187., 187., 188., 188., 189., 189., 190., 190., 191., 191., 192., 192., 193., 193., 194., 194., 195., 195., 196., 196.]])"
            ]
          },
          "metadata": {},
          "execution_count": 19
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now let us verify if this grouping method improves the accuracy."
      ],
      "metadata": {
        "id": "7XA-gY_D3xPs"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "groups = groups.reshape((1,784))[0]\n",
        "\n",
        "# Initialize and train the model\n",
        "model = LogisticGroupLasso(groups=groups,\n",
        "                           group_reg=0.05,  # Group regularization strength\n",
        "                           l1_reg=0.01,     # L1 regularization within groups\n",
        "                           n_iter=1000,\n",
        "                           tol=1e-4,\n",
        "                           supress_warning=True\n",
        ")"
      ],
      "metadata": {
        "id": "QuScjJeq303P"
      },
      "execution_count": 20,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "fit_summary(groups)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 750
        },
        "id": "WfBn2hu__XQX",
        "outputId": "52aa4df9-2508-4a15-91c1-9d6cbad74082"
      },
      "execution_count": 21,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Accuracy: 0.8074\n",
            "Classification Report:\n",
            "              precision    recall  f1-score   support\n",
            "\n",
            "           3       0.77      0.87      0.82      1426\n",
            "           5       0.82      0.74      0.78      1257\n",
            "           8       0.84      0.81      0.82      1373\n",
            "\n",
            "    accuracy                           0.81      4056\n",
            "   macro avg       0.81      0.80      0.81      4056\n",
            "weighted avg       0.81      0.81      0.81      4056\n",
            "\n",
            "There are 44 nonzero features after applying group LASSO.\n",
            "These are the pixels group LASSO found to be significant in classification of the digit to either 3, 5, or 8.\n",
            "\n",
            "Pixel List = ['pixel123', 'pixel124', 'pixel151', 'pixel152', 'pixel191', 'pixel192', 'pixel219', 'pixel220', 'pixel291', 'pixel292', 'pixel297', 'pixel298', 'pixel319', 'pixel320', 'pixel325', 'pixel326', 'pixel345', 'pixel346', 'pixel351', 'pixel352', 'pixel353', 'pixel354', 'pixel373', 'pixel374', 'pixel379', 'pixel380', 'pixel381', 'pixel382', 'pixel407', 'pixel408', 'pixel435', 'pixel436', 'pixel459', 'pixel460', 'pixel461', 'pixel462', 'pixel487', 'pixel488', 'pixel489', 'pixel490', 'pixel515', 'pixel516', 'pixel543', 'pixel544']\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGzCAYAAABpdMNsAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALjdJREFUeJzt3XtYVXW+x/HPBmErApsQFUgwFS/jeKm8kOnkqEyiTV7KW5phqaWhReI0WVOSOTHHSzl6TJ9zetLypJZNajllYyhqhTY6lWOdvDB4GwUvJSDGJfidP3zcp50ooht+oO/X86zncf/Wj7W+e7HcH9Zav72WwxhjBABANfOxXQAA4PpEAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAKHWy8nJ0ZAhQ9SgQQM5HA7NmzfPdkkALgMBhCrz1VdfyeFwaM+ePZKkl19+WTfddJPX1/PEE0/oo48+0rRp07Rs2TLFx8d7fR34f59//rkeffRRderUSX5+fnI4HBftm5OTowcffFCNGjVSvXr1dOutt2rVqlUX7f/WW2+pW7duql+/vkJCQnT77bdr48aNVfE2UAPUsV0Arl3bt29XaGioWrVqJUnKyMjQbbfd5vX1bNy4UQMHDtTUqVO9vmxc6IMPPtCrr76qDh06qHnz5tq7d2+5/fLy8tSjRw/l5OTo8ccfV3h4uN5++20NGzZMb775pkaOHOnRPyUlRTNmzNCQIUM0ZswYlZSUaPfu3fr3v/9dHW8LNhigiowdO9bEx8e7Xzdp0sS89NJLXl+Pw+EwiYmJXl/uz505c6bK11EbZGdnm7NnzxpjjElMTDQX+xiZNWuWkWTS0tLcbaWlpaZLly4mPDzcFBUVudszMjKMw+Gokv0DNRen4OBV33//vU6ePKmTJ09q+/btateunU6ePKmvv/5aR44cUcuWLXXy5EmdOXOmwmX961//0tChQxUaGqqAgADddttt+utf/+qev3TpUjkcDhljtHDhQjkcjkueDpKkU6dOafTo0QoODlZISIgSEhLcpwqXLl3q7jdmzBgFBgYqMzNT/fv3V1BQkEaNGiVJKigoUHJysqKiouR0OtW6dWvNmTNH5ic3lj9w4MAFyzzP4XAoJSXF/TolJUUOh0Pffvuthg0bpuDgYDVo0ECPP/64CgsLPX52w4YN6tGjh0JCQhQYGKjWrVvr6aefrnBbelPjxo1Vr169Cvtt3bpVDRs2VO/evd1tPj4+GjZsmLKzs7V582Z3+7x58xQeHq7HH39cxpjL2j9Q+3EKDl51yy236ODBg+7Xu3fv1pw5c9yv7777bklSQkJCuR/O5+Xk5Oj222/X2bNn9dhjj6lBgwZ6/fXXNWDAAL3zzjsaPHiw7rjjDi1btkyjR4/Wb37zGz3wwAOXrK2srEx33323Pv/8c02cOFFt2rTR2rVrlZCQUG7/H3/8UX379lWPHj00Z84cBQQEyBijAQMGaNOmTRo7dqxuvvlmffTRR/rd736nf//733r55ZcrsbU8DRs2TDfddJNSU1O1bds2zZ8/X99//73eeOMNSdLXX3+t3/72t+rQoYNmzJghp9Op/fv369NPP61w2bm5uSopKamwX926dRUYGHjF7+GnioqKyg2qgIAASdLOnTv1m9/8RpKUlpam22+/XfPnz9fMmTN16tQphYeH65lnntGkSZO8Ug9qILsHYLjWfPLJJ2bDhg3m2WefNXXq1DEffvih2bBhg+nXr5/p3Lmz2bBhg9mwYYP5+uuvL7mcpKQkI8ls3brV3Zafn2+aNWtmbrrpJlNaWupul3RZp+D+8pe/GElm3rx57rbS0lLTu3dvI8ksWbLE3Z6QkGAkmaeeespjGWvWrDGSzMyZMz3ahwwZYhwOh9m/f78xxpisrKwLlvnTeqdPn+5+PX36dCPJDBgwwKPfo48+aiSZr776yhhjzMsvv2wkmRMnTlT4Xn+uZ8+eRlKFU0JCQqWWe6lTcJMnTzY+Pj7mwIEDHu0jRowwksykSZOMMcZ89913RpJp0KCBCQwMNLNnzzZvvfWWiY+PN5LM4sWLK/1+UTtwBASv6t69u6RzF6q7dOniHpGWlJSkoUOHKi4u7rKW88EHH6hr167q0aOHuy0wMFAPP/ywpk2bpm+++Ubt2rWrVG3r16+Xn5+fxo8f727z8fFRYmLiRUdaTZw48YK6fH199dhjj3m0Jycn65133tGHH354xX+xJyYmeryePHmyXnnlFX3wwQfq0KGDQkJCJElr167Vgw8+KB+fyz+DPnfuXH3//fcV9ouMjKxUzZcybtw4LV68WMOGDdPLL7+sxo0b6+2339bq1aslST/88IMkuU+3nTp1SitXrtTw4cMlSUOGDFH79u01c+ZMPfLII16rCzUHAQSv+elpnrS0NPXu3VsnT57Ud999p6+//lozZ87UyZMn5efnJ5fLdcllHTx4ULGxsRe0/+IXv3DPr2wAHTx4UBEREe5TQOfFxMSU279OnTpq0qTJBcuIjIxUUFDQReu6Ui1btvR43aJFC/n4+OjAgQOSpOHDh+vVV1/VuHHj9NRTT6lPnz665557NGTIkArDqFOnTldc15Xq0KGDli9frgkTJrj/MAkPD9e8efM0ceJE96m+86fp/Pz8NGTIEPfP+/j4aPjw4Zo+fboOHTqk6Ojoan8PqFoEELxm4MCBHheWd+3a5fGl0MGDB0uSevbsqfT09GqurvKcTmeljjJ+6mKDIUpLS694GfXq1dOWLVu0adMm/fWvf9X69ev11ltvqXfv3vrb3/4mX1/fiy7ru+++U3FxcYXrrFevXoV/HFTGkCFDNGDAAH311VcqLS3Vrbfe6v7dnx+eHxoaqrp16yokJOSC99CoUSNJ5wa3EEDXHkbBwWvmzp2rDRs2aMaMGfL19dWHH36oDRs2aMCAAerUqZM2bNigDRs2aO7cuRUuq2nTpu4vsP7Ut99+655fWU2bNtWxY8d09uxZj/b9+/dXahlHjx5Vfn7+Jeu64YYbJEmnT5/26HepI6R9+/ZdUFdZWZnHl3d9fHzUp08fvfTSS/rmm2/0xz/+URs3btSmTZsuWfc999yjiIiICqfHH3/8ksu5Ev7+/urSpYtuu+02+fv76+OPP5Yk9+lYHx8f3XzzzTpx4sQFIXn06FFJUsOGDb1eF+wjgOA1nTp1UlxcnH788Ue1a9dO8fHxiouLU05OjuLi4tzT5ZwO6t+/vz7//HNlZGS42woKCvRf//Vfuummm9S2bdtK19e3b1+VlJTov//7v91tZWVlWrhw4WUvo3///iotLdV//ud/erS//PLLcjgc6tevnyQpODhYYWFh2rJli0e/V1555aLL/nkdCxYskCT3Mr/77rsLfubmm2+WdG7E2aWc/+OgounJJ5+85HKu1r59+7R48WL99re/dR8BSedOL5aWlur11193txUWFurNN99U27ZtvXptCjUHp+DgdZ9++qluv/12Sec+RL744otKf1flqaee0ooVK9SvXz899thjCg0N1euvv66srCz95S9/uaJTY4MGDVLXrl2VnJys/fv3q02bNnrvvffcH+wVfYdIOjeMvFevXnrmmWd04MABdezYUX/729+0du1aJSUlqUWLFu6+48aN05/+9CeNGzdOnTt31pYtWy561wBJysrK0oABAxQfH6+MjAz9z//8j0aOHKmOHTtKkmbMmKEtW7borrvuUtOmTXX8+HG98soratKkicdgjfJ48xrQwYMHtWzZMknSjh07JEkzZ86UdO4IcPTo0e6+bdu21dChQxUdHa2srCwtWrRIoaGhWrx4sccyH3nkEb366qtKTEzU3r17FR0drWXLlungwYN6//33vVY7ahjbw/Bwbfnxxx9NYGCgWbZsmTHm3LBsSeb48eOVXlZmZqYZMmSICQkJMXXr1jVdu3Y169atu6CfLnMYtjHGnDhxwowcOdIEBQUZl8tlxowZYz799FMjyaxcudLdLyEhwdSvX7/cZeTn55snnnjCREZGGj8/P9OyZUsze/ZsU1ZW5tHv7NmzZuzYscblcpmgoCAzbNgwc/z48YsOw/7mm2/MkCFDTFBQkLnhhhvMpEmTzA8//ODul5aWZgYOHGgiIyONv7+/iYyMNPfdd5/Zu3fvZb13b9m0adNFh3H37NnTo++IESNMVFSUu94JEyaYnJyccpebk5NjEhISTGhoqHE6nSY2NtasX7++Gt4RbHEY85OvbwPXoTVr1mjw4MH65JNP3KO1qlNKSoqef/55nThxQmFhYdW+fsAWrgHhunL+uyfnlZaWasGCBQoODtatt95qqSrg+sQ1IFxXJk+erB9++EHdunVTUVGR3n33XX322Wd68cUXL+v+ZgC8hwDCdaV3796aO3eu1q1bp8LCQsXExGjBggXcbwywgGtAAAAruAYEALCCAAIAWFHjrgGVlZXp6NGjCgoKuqwvBgIAahZjjPLz8xUZGXnJL43XuAA6evSooqKibJcBALhKhw8fvuCO8j9V4wLo/G3uDx8+rODgYMvVAAAqKy8vT1FRURc8tuTnqiyAFi5cqNmzZys7O1sdO3bUggUL1LVr1wp/7vxpt+DgYAIIAGqxii6jVMkghLfeektTpkzR9OnT9Y9//EMdO3ZU3759dfz48apYHQCgFqqSAHrppZc0fvx4Pfjgg2rbtq0WL16sgIAAvfbaa1WxOgBALeT1ACouLtbOnTvdD5uSzj1wKi4uzuPZLucVFRUpLy/PYwIAXPu8fg3o5MmTKi0tVePGjT3aGzdu7H5q5E+lpqbq+eef93YZAGqx0tJSlZSU2C4DF+Hr66s6depc9VdlrI+CmzZtmqZMmeJ+fX70BIDr05kzZ3TkyBFxl7CaLSAgQBEREfL397/iZXg9gMLCwuTr66ucnByP9pycHIWHh1/Q3+l0yul0ersMALVQaWmpjhw5ooCAADVs2JAvo9dAxhgVFxfrxIkTysrKUsuWLa/oCcVSFQSQv7+/OnXqpLS0NA0aNEjSubsbpKWlccdhAJdUUlIiY4waNmzI4zFqsHr16snPz08HDx5UcXGx6tate0XLqZJTcFOmTFFCQoI6d+6srl27at68eSooKNCDDz5YFasDcI3hyKfmu9Kjnp+qkgAaPny4Tpw4oeeee07Z2dm6+eabtX79+gsGJgAArl9VNghh0qRJnHID4BWHDh3SyZMnq219YWFhio6Orrb1Xa+sj4IDgEs5dOiQWrdurcLCwmpbZ926dbVnzx6vhdCYMWN0+vRprVmzxivLW7p0qZKSknT69GmvLM8WngcEoEY7efJktYaPJBUWFlbqiGvMmDFyOBxyOBzy9/dXTEyMZsyYoR9//FGS9Oc//1lLly6tomovbt26derZs6eCgoIUEBCgLl26WKnjYgggAPCC+Ph4HTt2TPv27VNycrJSUlI0e/ZsSZLL5VJISEi11rNgwQINHDhQ3bt31/bt27Vr1y6NGDFCEyZM0NSpU6u1loshgADAC5xOp8LDw9W0aVNNnDhRcXFxeu+99ySdO0I6/7WUEydOKDw8XC+++KL7Zz/77DP5+/srLS1N0rlblE2dOlU33nij6tevr9jYWKWnp192LYcPH1ZycrKSkpL04osvqm3btoqJiVFycrJmz56tuXPnavv27ZKk9PR0ORwOpaWlqXPnzgoICNDtt9+uPXv2eGfDXAIBBABVoF69eiouLr6gvWHDhnrttdeUkpKiHTt2KD8/X6NHj9akSZPUp08fSecGcWVkZGjlypXatWuXhg4dqvj4eO3bt++y1v3OO++opKSk3COdRx55RIGBgVqxYoVH+zPPPKO5c+dqx44dqlOnjh566KEreNeVwyAEAPAiY4zS0tL00UcfafLkyeX26d+/v8aPH69Ro0apc+fOql+/vlJTUyWdG3SxZMkSHTp0SJGRkZKkqVOnav369VqyZInHkdPF7N27Vy6XSxERERfM8/f3V/PmzbV3716P9j/+8Y/q2bOnJOmpp57SXXfdpcLCwiv+kunlIIAAwAvWrVunwMBAlZSUqKysTCNHjlRKSspF+8+ZM0ft2rXTqlWrtHPnTvctyf75z3+qtLRUrVq18uhfVFSkBg0aVFn9HTp0cP/7fHAdP368SoejE0AA4AW9evXSokWL5O/vr8jISNWpc+mP18zMTB09elRlZWU6cOCA2rdvL+nczVh9fX21c+dO+fr6evxMYGDgZdXSqlUr5ebm6ujRo+6jqPOKi4uVmZmpXr16ebT7+fm5/33+ThRlZWWXtb4rxTUgAPCC+vXrKyYmRtHR0RWGT3Fxse6//34NHz5cL7zwgsaNG+d+YvQtt9yi0tJSHT9+XDExMR5TeTd0Ls+9994rPz8/zZ0794J5ixcvVkFBge67777Kv0kv4wgIAKrZM888o9zcXM2fP1+BgYH64IMP9NBDD2ndunVq1aqVRo0apQceeEBz587VLbfcohMnTigtLU0dOnTQXXfdVeHyo6OjNWvWLCUnJ6tu3boaPXq0/Pz8tHbtWj399NNKTk5WbGxsNbzTS+MICECNFhYWVqUXwstTt25dhYWFVcmy09PTNW/ePC1btkzBwcHy8fHRsmXLtHXrVi1atEiStGTJEj3wwANKTk5W69atNWjQIP3973+v1PWYpKQkrV69Wlu3blXnzp3Vrl07LV++XIsWLdKcOXOq5L1VlsPUsKc+5eXlyeVyKTc3V8HBwbbLAVCNCgsLlZWVpWbNmnmEDveCq3ku9ruSLv9znFNwAGq86OhoAuEaxCk4AIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAKhiP30iqjcsXbrUq4/4/vTTT9W+fXv5+fm56yyvzdsIIAC4SmPGjJHD4ZDD4ZC/v79iYmI0Y8YM/fjjj5KkP//5z1q6dGm117Vu3Tr17NlTQUFBCggIUJcuXcqtY8qUKbr55puVlZXlnl9em7cRQADgBfHx8Tp27Jj27dun5ORkpaSkaPbs2ZIkl8vl1SOWy7FgwQINHDhQ3bt31/bt27Vr1y6NGDFCEyZMuOBR3ZmZmerdu7eaNGnirrO8Nm8jgADAC5xOp8LDw9W0aVNNnDhRcXFxeu+99yR5noI7ceKEwsPDPR6t/dlnn8nf319paWmSzj39dOrUqbrxxhtVv359xcbGKj09/bJrOXz4sJKTk5WUlKQXX3xRbdu2VUxMjJKTkzV79mzNnTtX27dv14EDB+RwOHTq1Ck99NBDcjgcWrp0abltVYEAAoAqUK9ePRUXF1/Q3rBhQ7322mtKSUnRjh07lJ+fr9GjR2vSpEnq06ePJGnSpEnKyMjQypUrtWvXLg0dOlTx8fHat2/fZa37nXfeUUlJyQVHOpL0yCOPKDAwUCtWrFBUVJSOHTum4OBgzZs3T8eOHdPQoUMvaBs+fPjVbYyL4G7YAOBFxhilpaXpo48+0uTJk8vt079/f40fP16jRo1S586dVb9+faWmpko69+iJJUuW6NChQ+7HaU+dOlXr16/XkiVLPI6cLmbv3r1yuVyKiIi4YJ6/v7+aN2+uvXv3ytfXV+Hh4XI4HHK5XO4nrtavX/+CtqpAAAGAF6xbt06BgYEqKSlRWVmZRo4cqZSUlIv2nzNnjtq1a6dVq1Zp586dcjqdkqR//vOfKi0tVatWrTz6FxUVqUGDBlX5FqodAQQAXtCrVy8tWrRI/v7+ioyMVJ06l/54zczM1NGjR1VWVqYDBw6offv2kqQzZ87I19dXO3fulK+vr8fPBAYGXlYtrVq1Um5uro4ePeo+ijqvuLhYmZmZ6tWrVyXeXdXgGhAAeEH9+vUVExOj6OjoCsOnuLhY999/v4YPH64XXnhB48aN0/HjxyVJt9xyi0pLS3X8+HHFxMR4TJd7Ouzee++Vn5+f5s6de8G8xYsXq6CgQPfdd1/l36SXcQQEANXsmWeeUW5urubPn6/AwEB98MEHeuihh7Ru3Tq1atVKo0aN0gMPPKC5c+fqlltu0YkTJ5SWlqYOHTrorrvuqnD50dHRmjVrlpKTk1W3bl2NHj1afn5+Wrt2rZ5++mklJycrNja2Gt7ppXEEBADVKD09XfPmzdOyZcsUHBwsHx8fLVu2TFu3btWiRYskSUuWLNEDDzyg5ORktW7dWoMGDdLf//73Sj2WPCkpSatXr9bWrVvVuXNntWvXTsuXL9eiRYs0Z86cqnp7leIwxhjbRfxUXl6eXC6XcnNzFRwcbLscANWosLBQWVlZatasmerWrWu7HFzCpX5Xl/s5zhEQAMAKAggAYAUBBACwggACAFhBAAGocWrY2CiUwxu/IwIIQI1x/pv/5d3EEzXL2bNnJUl+fn5XvAy+iIoaz+FwVNu6+Mvbrjp16iggIEAnTpyQn5+ffHz4G7mmMcbo7NmzOn78uEJCQi64XVBlEEAAagyHw6GIiAhlZWXp4MGDtsvBJYSEhFz1nbIJIAA1ir+/v1q2bMlpuBrMz8/vqo58ziOAANQ4Pj4+3AnhOsAJVgCAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwwusBlJKSIofD4TG1adPG26sBANRyVXI37F/+8pf6+OOP/38ldbjpNgDAU5UkQ506da76QUUAgGtblVwD2rdvnyIjI9W8eXONGjVKhw4dumjfoqIi5eXleUwAgGuf1wMoNjZWS5cu1fr167Vo0SJlZWXpV7/6lfLz88vtn5qaKpfL5Z6ioqK8XRIAoAZyGGNMVa7g9OnTatq0qV566SWNHTv2gvlFRUUqKipyv87Ly1NUVJRyc3MVHBxclaWhlnA4HNW2rir+7wBcF/Ly8uRyuSr8HK/y0QEhISFq1aqV9u/fX+58p9Mpp9NZ1WUAAGqYKv8e0JkzZ5SZmamIiIiqXhUAoBbxegBNnTpVmzdv1oEDB/TZZ59p8ODB8vX11X333eftVQEAajGvn4I7cuSI7rvvPp06dUoNGzZUjx49tG3bNjVs2NDbqwIA1GJeD6CVK1d6e5G4zjEwALZU1wCY63Uf515wAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBV1bBcAADWVMcZ2Cdc0joAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFbUsV0AUBGHw1Ft6zLGVNu6qgPbDjUZR0AAACsIIACAFZUOoC1btujuu+9WZGSkHA6H1qxZ4zHfGKPnnntOERERqlevnuLi4rRv3z5v1QsAuEZUOoAKCgrUsWNHLVy4sNz5s2bN0vz587V48WJt375d9evXV9++fVVYWHjVxQIArh2VHoTQr18/9evXr9x5xhjNmzdPf/jDHzRw4EBJ0htvvKHGjRtrzZo1GjFixNVVCwC4Znj1GlBWVpays7MVFxfnbnO5XIqNjVVGRka5P1NUVKS8vDyPCQBw7fNqAGVnZ0uSGjdu7NHeuHFj97yfS01Nlcvlck9RUVHeLAkAUENZHwU3bdo05ebmuqfDhw/bLgkAUA28GkDh4eGSpJycHI/2nJwc97yfczqdCg4O9pgAANc+rwZQs2bNFB4errS0NHdbXl6etm/frm7dunlzVQCAWq7So+DOnDmj/fv3u19nZWXpyy+/VGhoqKKjo5WUlKSZM2eqZcuWatasmZ599llFRkZq0KBB3qwbAFDLVTqAduzYoV69erlfT5kyRZKUkJCgpUuX6sknn1RBQYEefvhhnT59Wj169ND69etVt25d71UNAKj1HKaG3UEwLy9PLpdLubm5XA+CJG6oeTXYdrDhcj/HrY+CAwBcnwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCi0o9jAM6rrjstX4t3Wa7Ou1QDNRVHQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAK+rYLgC1lzHGdgmoAL8j1GQcAQEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwotIBtGXLFt19992KjIyUw+HQmjVrPOaPGTNGDofDY4qPj/dWvQCAa0SlA6igoEAdO3bUwoULL9onPj5ex44dc08rVqy4qiIBANeeSj8RtV+/furXr98l+zidToWHh19xUQCAa1+VXANKT09Xo0aN1Lp1a02cOFGnTp26aN+ioiLl5eV5TACAa5/XAyg+Pl5vvPGG0tLS9B//8R/avHmz+vXrp9LS0nL7p6amyuVyuaeoqChvlwQAqIEcxhhzxT/scGj16tUaNGjQRfv861//UosWLfTxxx+rT58+F8wvKipSUVGR+3VeXp6ioqKUm5ur4ODgKy0NqNEcDke1rOcq/nsDVywvL08ul6vCz/EqH4bdvHlzhYWFaf/+/eXOdzqdCg4O9pgAANe+Kg+gI0eO6NSpU4qIiKjqVQEAapFKj4I7c+aMx9FMVlaWvvzyS4WGhio0NFTPP/+87r33XoWHhyszM1NPPvmkYmJi1LdvX68WDgCo3SodQDt27FCvXr3cr6dMmSJJSkhI0KJFi7Rr1y69/vrrOn36tCIjI3XnnXfqhRdekNPp9F7VAIBa76oGIVSFy714hetHdV2wr0417L8d4FU1ZhACAADlIYAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFbUsV0AUJMYY2yXAFw3OAICAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhRx3YBQEWMMbZLAFAFOAICAFhBAAEArKhUAKWmpqpLly4KCgpSo0aNNGjQIO3Zs8ejT2FhoRITE9WgQQMFBgbq3nvvVU5OjleLBgDUfpUKoM2bNysxMVHbtm3Thg0bVFJSojvvvFMFBQXuPk888YTef/99rVq1Sps3b9bRo0d1zz33eL1wAEDt5jBXcYX3xIkTatSokTZv3qw77rhDubm5atiwoZYvX64hQ4ZIkr799lv94he/UEZGhm677bYKl5mXlyeXy6Xc3FwFBwdfaWkAAEsu93P8qq4B5ebmSpJCQ0MlSTt37lRJSYni4uLcfdq0aaPo6GhlZGSUu4yioiLl5eV5TACAa98VB1BZWZmSkpLUvXt3tWvXTpKUnZ0tf39/hYSEePRt3LixsrOzy11OamqqXC6Xe4qKirrSkgAAtcgVB1BiYqJ2796tlStXXlUB06ZNU25urns6fPjwVS0PAFA7XNEXUSdNmqR169Zpy5YtatKkibs9PDxcxcXFOn36tMdRUE5OjsLDw8tdltPplNPpvJIyAAC1WKWOgIwxmjRpklavXq2NGzeqWbNmHvM7deokPz8/paWludv27NmjQ4cOqVu3bt6pGABwTajUEVBiYqKWL1+utWvXKigoyH1dx+VyqV69enK5XBo7dqymTJmi0NBQBQcHa/LkyerWrdtljYADAFw/KjUM2+FwlNu+ZMkSjRkzRtK5L6ImJydrxYoVKioqUt++ffXKK69c9BTczzEMGwBqt8v9HL+q7wFVBQIIAGq3avkeEAAAV4oAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADACgIIAGAFAQQAsIIAAgBYQQABAKwggAAAVhBAAAArCCAAgBUEEADAikoFUGpqqrp06aKgoCA1atRIgwYN0p49ezz6/PrXv5bD4fCYJkyY4NWiAQC1X6UCaPPmzUpMTNS2bdu0YcMGlZSU6M4771RBQYFHv/Hjx+vYsWPuadasWV4tGgBQ+9WpTOf169d7vF66dKkaNWqknTt36o477nC3BwQEKDw83DsVAgCuSVd1DSg3N1eSFBoa6tH+5ptvKiwsTO3atdO0adN09uzZiy6jqKhIeXl5HhMA4NpXqSOgnyorK1NSUpK6d++udu3audtHjhyppk2bKjIyUrt27dLvf/977dmzR++++265y0lNTdXzzz9/pWUAAGophzHGXMkPTpw4UR9++KE++eQTNWnS5KL9Nm7cqD59+mj//v1q0aLFBfOLiopUVFTkfp2Xl6eoqCjl5uYqODj4SkoDAFiUl5cnl8tV4ef4FR0BTZo0SevWrdOWLVsuGT6SFBsbK0kXDSCn0ymn03klZQAAarFKBZAxRpMnT9bq1auVnp6uZs2aVfgzX375pSQpIiLiigoEAFybKhVAiYmJWr58udauXaugoCBlZ2dLklwul+rVq6fMzEwtX75c/fv3V4MGDbRr1y498cQTuuOOO9ShQ4cqeQMAgNqpUteAHA5Hue1LlizRmDFjdPjwYd1///3avXu3CgoKFBUVpcGDB+sPf/jDZV/PudxzhwCAmqlKrgFVlFVRUVHavHlzZRYJALhOcS84AIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVBBAAwAoCCABgBQEEALCCAAIAWEEAAQCsIIAAAFYQQAAAKwggAIAVdWwX8HPGGElSXl6e5UoAAFfi/Of3+c/zi6lxAZSfny9JioqKslwJAOBq5Ofny+VyXXS+w1QUUdWsrKxMR48eVVBQkBwOh8e8vLw8RUVF6fDhwwoODrZUoX1sh3PYDuewHc5hO5xTE7aDMUb5+fmKjIyUj8/Fr/TUuCMgHx8fNWnS5JJ9goODr+sd7Dy2wzlsh3PYDuewHc6xvR0udeRzHoMQAABWEEAAACtqVQA5nU5Nnz5dTqfTdilWsR3OYTucw3Y4h+1wTm3aDjVuEAIA4PpQq46AAADXDgIIAGAFAQQAsIIAAgBYQQABAKyoNQG0cOFC3XTTTapbt65iY2P1+eef2y6p2qWkpMjhcHhMbdq0sV1WlduyZYvuvvtuRUZGyuFwaM2aNR7zjTF67rnnFBERoXr16ikuLk779u2zU2wVqmg7jBkz5oL9Iz4+3k6xVSQ1NVVdunRRUFCQGjVqpEGDBmnPnj0efQoLC5WYmKgGDRooMDBQ9957r3JycixVXDUuZzv8+te/vmB/mDBhgqWKy1crAuitt97SlClTNH36dP3jH/9Qx44d1bdvXx0/ftx2adXul7/8pY4dO+aePvnkE9slVbmCggJ17NhRCxcuLHf+rFmzNH/+fC1evFjbt29X/fr11bdvXxUWFlZzpVWrou0gSfHx8R77x4oVK6qxwqq3efNmJSYmatu2bdqwYYNKSkp05513qqCgwN3niSee0Pvvv69Vq1Zp8+bNOnr0qO655x6LVXvf5WwHSRo/frzH/jBr1ixLFV+EqQW6du1qEhMT3a9LS0tNZGSkSU1NtVhV9Zs+fbrp2LGj7TKskmRWr17tfl1WVmbCw8PN7Nmz3W2nT582TqfTrFixwkKF1ePn28EYYxISEszAgQOt1GPL8ePHjSSzefNmY8y5372fn59ZtWqVu8///u//GkkmIyPDVplV7ufbwRhjevbsaR5//HF7RV2GGn8EVFxcrJ07dyouLs7d5uPjo7i4OGVkZFiszI59+/YpMjJSzZs316hRo3To0CHbJVmVlZWl7Oxsj/3D5XIpNjb2utw/0tPT1ahRI7Vu3VoTJ07UqVOnbJdUpXJzcyVJoaGhkqSdO3eqpKTEY39o06aNoqOjr+n94efb4bw333xTYWFhateunaZNm6azZ8/aKO+iatzdsH/u5MmTKi0tVePGjT3aGzdurG+//dZSVXbExsZq6dKlat26tY4dO6bnn39ev/rVr7R7924FBQXZLs+K7OxsSSp3/zg/73oRHx+ve+65R82aNVNmZqaefvpp9evXTxkZGfL19bVdnteVlZUpKSlJ3bt3V7t27SSd2x/8/f0VEhLi0fda3h/K2w6SNHLkSDVt2lSRkZHatWuXfv/732vPnj169913LVbrqcYHEP5fv3793P/u0KGDYmNj1bRpU7399tsaO3asxcpQE4wYMcL97/bt26tDhw5q0aKF0tPT1adPH4uVVY3ExETt3r37urgOeikX2w4PP/yw+9/t27dXRESE+vTpo8zMTLVo0aK6yyxXjT8FFxYWJl9f3wtGseTk5Cg8PNxSVTVDSEiIWrVqpf3799suxZrz+wD7x4WaN2+usLCwa3L/mDRpktatW6dNmzZ5PD8sPDxcxcXFOn36tEf/a3V/uNh2KE9sbKwk1aj9ocYHkL+/vzp16qS0tDR3W1lZmdLS0tStWzeLldl35swZZWZmKiIiwnYp1jRr1kzh4eEe+0deXp62b99+3e8fR44c0alTp66p/cMYo0mTJmn16tXauHGjmjVr5jG/U6dO8vPz89gf9uzZo0OHDl1T+0NF26E8X375pSTVrP3B9iiIy7Fy5UrjdDrN0qVLzTfffGMefvhhExISYrKzs22XVq2Sk5NNenq6ycrKMp9++qmJi4szYWFh5vjx47ZLq1L5+fnmiy++MF988YWRZF566SXzxRdfmIMHDxpjjPnTn/5kQkJCzNq1a82uXbvMwIEDTbNmzcwPP/xguXLvutR2yM/PN1OnTjUZGRkmKyvLfPzxx+bWW281LVu2NIWFhbZL95qJEycal8tl0tPTzbFjx9zT2bNn3X0mTJhgoqOjzcaNG82OHTtMt27dTLdu3SxW7X0VbYf9+/ebGTNmmB07dpisrCyzdu1a07x5c3PHHXdYrtxTrQggY4xZsGCBiY6ONv7+/qZr165m27ZttkuqdsOHDzcRERHG39/f3HjjjWb48OFm//79tsuqcps2bTKSLpgSEhKMMeeGYj/77LOmcePGxul0mj59+pg9e/bYLboKXGo7nD171tx5552mYcOGxs/PzzRt2tSMHz/+mvsjrbz3L8ksWbLE3eeHH34wjz76qLnhhhtMQECAGTx4sDl27Ji9oqtARdvh0KFD5o477jChoaHG6XSamJgY87vf/c7k5ubaLfxneB4QAMCKGn8NCABwbSKAAABWEEAAACsIIACAFQQQAMAKAggAYAUBBACwggACAFhBAAEArCCAAABWEEAAACv+D+tWG2pXFCrEAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Looks good.  I think all of these methods are acceptable.  Hyperparameter tuning can be done to optimize the parameters group_regularization and l1_regularization (as well as group selection like what we have done here on a naive level)."
      ],
      "metadata": {
        "id": "dEbIlG73F5zG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "*End of Notebook 2*"
      ],
      "metadata": {
        "id": "pqd_bbZ-GM_T"
      }
    }
  ]
}