Symulowanie-wizualne/sw_lab3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bf9393fa-7028-496b-9090-d7a1886364ce",
   "metadata": {},
   "source": [
    "### Regularyzacja przez mnożniki Lagrange'a. Algorytm SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "92fd0e85",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import subprocess\n",
    "import pkg_resources\n",
    "import numpy as np\n",
    "\n",
    "required = { 'scikit-image'}\n",
    "installed = {pkg.key for pkg in pkg_resources.working_set}\n",
    "missing = required - installed\n",
    "\n",
    "if missing: \n",
    "    python = sys.executable\n",
    "    subprocess.check_call([python, '-m', 'pip', 'install', *missing], stdout=subprocess.DEVNULL)\n",
    "\n",
    "def load_train_data(input_dir, newSize=(64,64)):\n",
    "    import numpy as np\n",
    "    import pandas as pd\n",
    "    import os\n",
    "    from skimage.io import imread\n",
    "    import cv2 as cv\n",
    "    from pathlib import Path\n",
    "    import random\n",
    "    from shutil import copyfile, rmtree\n",
    "    import json\n",
    "\n",
    "    import seaborn as sns\n",
    "    import matplotlib.pyplot as plt\n",
    "\n",
    "    import matplotlib\n",
    "    \n",
    "    image_dir = Path(input_dir)\n",
    "    categories_name = []\n",
    "    for file in os.listdir(image_dir):\n",
    "        d = os.path.join(image_dir, file)\n",
    "        if os.path.isdir(d):\n",
    "            categories_name.append(file)\n",
    "\n",
    "    folders = [directory for directory in image_dir.iterdir() if directory.is_dir()]\n",
    "\n",
    "    train_img = []\n",
    "    categories_count=[]\n",
    "    labels=[]\n",
    "    for i, direc in enumerate(folders):\n",
    "        count = 0\n",
    "        for obj in direc.iterdir():\n",
    "            if os.path.isfile(obj) and os.path.basename(os.path.normpath(obj)) != 'desktop.ini':\n",
    "                labels.append(os.path.basename(os.path.normpath(direc)))\n",
    "                count += 1\n",
    "                img = imread(obj)#zwraca ndarry postaci xSize x ySize x colorDepth\n",
    "                img = cv.resize(img, newSize, interpolation=cv.INTER_AREA)# zwraca ndarray\n",
    "                img = img / 255#normalizacja\n",
    "                train_img.append(img)\n",
    "        categories_count.append(count)\n",
    "    X={}\n",
    "    X[\"values\"] = np.array(train_img)\n",
    "    X[\"categories_name\"] = categories_name\n",
    "    X[\"categories_count\"] = categories_count\n",
    "    X[\"labels\"]=labels\n",
    "    return X\n",
    "\n",
    "def load_test_data(input_dir, newSize=(64,64)):\n",
    "    import numpy as np\n",
    "    import pandas as pd\n",
    "    import os\n",
    "    from skimage.io import imread\n",
    "    import cv2 as cv\n",
    "    from pathlib import Path\n",
    "    import random\n",
    "    from shutil import copyfile, rmtree\n",
    "    import json\n",
    "\n",
    "    import seaborn as sns\n",
    "    import matplotlib.pyplot as plt\n",
    "\n",
    "    import matplotlib\n",
    "\n",
    "    image_path = Path(input_dir)\n",
    "\n",
    "    labels_path = image_path.parents[0] / 'test_labels.json'\n",
    "\n",
    "    jsonString = labels_path.read_text()\n",
    "    objects = json.loads(jsonString)\n",
    "\n",
    "    categories_name = []\n",
    "    categories_count=[]\n",
    "    count = 0\n",
    "    c = objects[0]['value']\n",
    "    for e in  objects:\n",
    "        if e['value'] != c:\n",
    "            categories_count.append(count)\n",
    "            c = e['value']\n",
    "            count = 1\n",
    "        else:\n",
    "            count += 1\n",
    "        if not e['value'] in categories_name:\n",
    "            categories_name.append(e['value'])\n",
    "\n",
    "    categories_count.append(count)\n",
    "    \n",
    "    test_img = []\n",
    "\n",
    "    labels=[]\n",
    "    for e in objects:\n",
    "        p = image_path / e['filename']\n",
    "        img = imread(p)#zwraca ndarry postaci xSize x ySize x colorDepth\n",
    "        img = cv.resize(img, newSize, interpolation=cv.INTER_AREA)# zwraca ndarray\n",
    "        img = img / 255#normalizacja\n",
    "        test_img.append(img)\n",
    "        labels.append(e['value'])\n",
    "\n",
    "    X={}\n",
    "    X[\"values\"] = np.array(test_img)\n",
    "    X[\"categories_name\"] = categories_name\n",
    "    X[\"categories_count\"] = categories_count\n",
    "    X[\"labels\"]=labels\n",
    "    return X\n",
    "\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "\n",
    "# Data load\n",
    "data_train = load_train_data(\"train_test_sw/train_sw\", newSize=(16,16))\n",
    "X_train = data_train['values']\n",
    "y_train = data_train['labels']\n",
    "\n",
    "data_test = load_test_data(\"train_test_sw/test_sw\", newSize=(16,16))\n",
    "X_test = data_test['values']\n",
    "y_test = data_test['labels']\n",
    "\n",
    "class_le = LabelEncoder()\n",
    "y_train_enc = class_le.fit_transform(y_train)\n",
    "y_test_enc = class_le.fit_transform(y_test)\n",
    "\n",
    "X_train = X_train.flatten().reshape(X_train.shape[0], int(np.prod(X_train.shape) / X_train.shape[0]))\n",
    "X_test = X_test.flatten().reshape(X_test.shape[0], int(np.prod(X_test.shape) / X_test.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e63a5e75-7a26-4a97-b53d-67b858345126",
   "metadata": {},
   "source": [
    "#### 1.  Zadanie 1 (2pkt): \n",
    "\n",
    "Rozwiń algorytm regresji logistycznej z lab. 1, wprowadzając do niego człon regularyzacyjny"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ea300c45",
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogisticRegressionL2():\n",
    "    def __init__(self, l2=1):\n",
    "        self.l2 = l2\n",
    "\n",
    "    def mapY(self, y, cls):\n",
    "        m = len(y)\n",
    "        yBi = np.matrix(np.zeros(m)).reshape(m, 1)\n",
    "        yBi[y == cls] = 1.\n",
    "        return yBi\n",
    "\n",
    "    def indicatorMatrix(self, y):\n",
    "        classes = np.unique(y.tolist())\n",
    "        m = len(y)\n",
    "        k = len(classes)\n",
    "        Y = np.matrix(np.zeros((m, k)))\n",
    "        for i, cls in enumerate(classes):\n",
    "            Y[:, i] = self.mapY(y, cls)\n",
    "        return Y\n",
    "    \n",
    "    # Zapis macierzowy funkcji softmax\n",
    "    def softmax(self, X):\n",
    "        return np.exp(X) / np.sum(np.exp(X))\n",
    "    \n",
    "    # Funkcja regresji logistcznej\n",
    "    def h(self, theta, X):\n",
    "        return 1.0/(1.0 + np.exp(-X * theta))\n",
    "    \n",
    "    # Funkcja kosztu dla regresji logistycznej\n",
    "    def J(self, h, theta, X, y):\n",
    "        m = len(y)\n",
    "        h_val = h(theta, X)\n",
    "        s1 = np.multiply(y, np.log(h_val))\n",
    "        s2 = np.multiply((1 - y), np.log(1 - h_val))\n",
    "        s3 = np.sum(s1+s2, axis=0)/m\n",
    "        s4 = (self.l2 * np.sum(np.square(theta))) / 2*m\n",
    "        return -s3 + s4\n",
    "\n",
    "    # Gradient dla regresji logistycznej\n",
    "    def dJ(self, h, theta, X, y):\n",
    "        return 1.0 / (self.l2/len(y)) * (X.T * (h(theta, X) - y))\n",
    "\n",
    "    # Metoda gradientu prostego dla regresji logistycznej\n",
    "    def GD(self, h, fJ, fdJ, theta, X, y, alpha=0.01, eps=10**-3, maxSteps=10000):\n",
    "        errorCurr = fJ(h, theta, X, y)\n",
    "        errors = [[errorCurr, theta]]\n",
    "        while True:\n",
    "            # oblicz nowe theta\n",
    "            theta = theta - alpha * fdJ(h, theta, X, y)\n",
    "            # raportuj poziom błędu\n",
    "            errorCurr, errorPrev = fJ(h, theta, X, y), errorCurr\n",
    "            # kryteria stopu\n",
    "            if abs(errorPrev - errorCurr) <= eps:\n",
    "                break\n",
    "            if len(errors) > maxSteps:\n",
    "                break\n",
    "            errors.append([errorCurr, theta]) \n",
    "        return theta, errors\n",
    "\n",
    "    def trainMaxEnt(self, X, Y):\n",
    "        n = X.shape[1]\n",
    "        thetas = []\n",
    "        for c in range(Y.shape[1]):\n",
    "            YBi = Y[:,c]\n",
    "            theta = np.matrix(np.random.random(n)).reshape(n,1)\n",
    "            # Macierz parametrów theta obliczona dla każdej klasy osobno.\n",
    "            thetaBest, errors = self.GD(self.h, self.J, self.dJ, theta, \n",
    "                                X, YBi, alpha=0.1, eps=10**-4)\n",
    "            thetas.append(thetaBest)\n",
    "        return thetas\n",
    "\n",
    "    def classify(self, thetas, X):\n",
    "        regs = np.array([(X*theta).item() for theta in thetas])\n",
    "        probs = self.softmax(regs)\n",
    "        result = np.argmax(probs)\n",
    "        return result\n",
    "\n",
    "    def accuracy(self, expected, predicted):\n",
    "        return sum(1 for x, y in zip(expected, predicted) if x == y) / len(expected)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "84fc2187",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jonas\\AppData\\Local\\Temp\\ipykernel_16900\\514332287.py:33: RuntimeWarning: divide by zero encountered in log\n",
      "  s2 = np.multiply((1 - y), np.log(1 - h_val))\n",
      "C:\\Users\\jonas\\AppData\\Local\\Temp\\ipykernel_16900\\514332287.py:33: RuntimeWarning: invalid value encountered in multiply\n",
      "  s2 = np.multiply((1 - y), np.log(1 - h_val))\n",
      "C:\\Users\\jonas\\AppData\\Local\\Temp\\ipykernel_16900\\514332287.py:26: RuntimeWarning: overflow encountered in exp\n",
      "  return 1.0/(1.0 + np.exp(-X * theta))\n",
      "C:\\Users\\jonas\\AppData\\Local\\Temp\\ipykernel_16900\\514332287.py:32: RuntimeWarning: divide by zero encountered in log\n",
      "  s1 = np.multiply(y, np.log(h_val))\n",
      "C:\\Users\\jonas\\AppData\\Local\\Temp\\ipykernel_16900\\514332287.py:32: RuntimeWarning: invalid value encountered in multiply\n",
      "  s1 = np.multiply(y, np.log(h_val))\n",
      "C:\\Users\\jonas\\AppData\\Local\\Temp\\ipykernel_16900\\514332287.py:22: RuntimeWarning: overflow encountered in exp\n",
      "  return np.exp(X) / np.sum(np.exp(X))\n",
      "C:\\Users\\jonas\\AppData\\Local\\Temp\\ipykernel_16900\\514332287.py:22: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.exp(X) / np.sum(np.exp(X))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.5714285714285714"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_reg = LogisticRegressionL2(l2 = 0.1)\n",
    "\n",
    "Y = log_reg.indicatorMatrix(y_train_enc)\n",
    "thetas = log_reg.trainMaxEnt(X_train, Y)\n",
    "\n",
    "predicted = [log_reg.classify(thetas, x) for x in X_test]\n",
    "\n",
    "log_reg.accuracy(y_test_enc, np.array(predicted))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "945d1169-a13c-44a5-ad51-73ec62438487",
   "metadata": {},
   "source": [
    "#### Zadanie 2 (4pkt)\n",
    "\n",
    "Zaimplementuj algorytm SVM z miękkim marginesem (regularyzacją)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7a58dff1-f85d-4e07-b5e4-b639fc2af3ac",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "d0e86f0fc352ea3a95940ad21ba9cd6c",
     "grade": true,
     "grade_id": "cell-2e22217733fd5f34",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.base import BaseEstimator, ClassifierMixin\n",
    "from sklearn.utils import check_random_state\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "\n",
    "\n",
    "def projection_simplex(v, z=1):\n",
    "    n_features = v.shape[0]\n",
    "    u = np.sort(v)[::-1]\n",
    "    cssv = np.cumsum(u) - z\n",
    "    ind = np.arange(n_features) + 1\n",
    "    cond = u - cssv / ind > 0\n",
    "    rho = ind[cond][-1]\n",
    "    theta = cssv[cond][-1] / float(rho)\n",
    "    w = np.maximum(v - theta, 0)\n",
    "    return w\n",
    "\n",
    "\n",
    "class MulticlassSVM(BaseEstimator, ClassifierMixin):\n",
    "\n",
    "    def __init__(self, C=1, max_iter=50, tol=0.05,\n",
    "                 random_state=None, verbose=0):\n",
    "        self.C = C\n",
    "        self.max_iter = max_iter\n",
    "        self.tol = tol,\n",
    "        self.random_state = random_state\n",
    "        self.verbose = verbose\n",
    "\n",
    "    def _partial_gradient(self, X, y, i):\n",
    "        # Partial gradient for the ith sample.\n",
    "        g = np.dot(X[i], self.coef_.T) + 1\n",
    "        g[y[i]] -= 1\n",
    "        return g\n",
    "\n",
    "    def _violation(self, g, y, i):\n",
    "        # Optimality violation for the ith sample.\n",
    "        smallest = np.inf\n",
    "        for k in range(g.shape[0]):\n",
    "            if k == y[i] and self.dual_coef_[k, i] >= self.C:\n",
    "                continue\n",
    "            elif k != y[i] and self.dual_coef_[k, i] >= 0:\n",
    "                continue\n",
    "\n",
    "            smallest = min(smallest, g[k])\n",
    "\n",
    "        return g.max() - smallest\n",
    "\n",
    "    def _solve_subproblem(self, g, y, norms, i):\n",
    "        # Prepare inputs to the projection.\n",
    "        Ci = np.zeros(g.shape[0])\n",
    "        Ci[y[i]] = self.C\n",
    "        beta_hat = norms[i] * (Ci - self.dual_coef_[:, i]) + g / norms[i]\n",
    "        z = self.C * norms[i]\n",
    "\n",
    "        # Compute projection onto the simplex.\n",
    "        beta = projection_simplex(beta_hat, z)\n",
    "\n",
    "        return Ci - self.dual_coef_[:, i] - beta / norms[i]\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        n_samples, n_features = X.shape\n",
    "\n",
    "        # Normalize labels.\n",
    "        self._label_encoder = LabelEncoder()\n",
    "        y = self._label_encoder.fit_transform(y)\n",
    "\n",
    "        # Initialize primal and dual coefficients.\n",
    "        n_classes = len(self._label_encoder.classes_)\n",
    "        self.dual_coef_ = np.zeros((n_classes, n_samples), dtype=np.float64)\n",
    "        self.coef_ = np.zeros((n_classes, n_features))\n",
    "\n",
    "        # Pre-compute norms.\n",
    "        norms = np.sqrt(np.sum(X ** 2, axis=1))\n",
    "\n",
    "        # Shuffle sample indices.\n",
    "        rs = check_random_state(self.random_state)\n",
    "        ind = np.arange(n_samples)\n",
    "        rs.shuffle(ind)\n",
    "\n",
    "        violation_init = None\n",
    "        for it in range(self.max_iter):\n",
    "\n",
    "            for ii in range(n_samples):\n",
    "                i = ind[ii]\n",
    "\n",
    "                # All-zero samples can be safely ignored.\n",
    "                if norms[i] == 0:\n",
    "                    continue\n",
    "\n",
    "                g = self._partial_gradient(X, y, i)\n",
    "                v = self._violation(g, y, i)\n",
    "\n",
    "                if v < 1e-12:\n",
    "                    continue\n",
    "\n",
    "                # Solve subproblem for the ith sample.\n",
    "                delta = self._solve_subproblem(g, y, norms, i)\n",
    "\n",
    "                # Update primal and dual coefficients.\n",
    "                self.coef_ += (delta * X[i][:, np.newaxis]).T\n",
    "                self.dual_coef_[:, i] += delta\n",
    "\n",
    "                \n",
    "        return self\n",
    "\n",
    "    def predict(self, X):\n",
    "        decision = np.dot(X, self.coef_.T)\n",
    "        pred = decision.argmax(axis=1)\n",
    "        return self._label_encoder.inverse_transform(pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c694fbbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = data_train['values'][:,:,:,:3]\n",
    "X_train_flatten = np.array([x.flatten() for x in X_train])\n",
    "y_train = data_train['labels']\n",
    "\n",
    "X_test = data_test['values'][:,:,:,:3]\n",
    "X_test_flatten = np.array([x.flatten() for x in X_test])\n",
    "y_test = data_test['labels']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "62857aa5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6988416988416989\n"
     ]
    }
   ],
   "source": [
    "X, y = X_train_flatten, y_train\n",
    "\n",
    "clf = MulticlassSVM(C=0.1, tol=0.01, max_iter=100, random_state=0, verbose=1)\n",
    "clf.fit(X, y)\n",
    "print(clf.score(X_test_flatten, y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63ba4cfd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "7e1998ff7f8aa20ada591c520b972326324e5ea05489af9e422744c7c09f6dad"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}