TP optimization gradient descent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP de SD 211\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.optimize\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_matrix_from_faces(folder='orl_faces', minidata=False):\n",
    "    # load images\n",
    "    # 400 images of size (112, 92)\n",
    "    M = []\n",
    "    if minidata is True:\n",
    "        nb_subjects = 1\n",
    "    else:\n",
    "        nb_subjects = 40\n",
    "    for subject in range(1, nb_subjects + 1\n",
    "                        ):\n",
    "        for image in range(1, 11):\n",
    "            face = plt.imread(folder + '/s' + str(subject)\n",
    "                              + '/' + str(image) + '.pgm')\n",
    "            M.append(face.ravel())\n",
    "\n",
    "    return np.array(M, dtype=float)\n",
    "\n",
    "def vectorize(W, H):\n",
    "    return np.concatenate((W.ravel(), H.ravel()))\n",
    "\n",
    "def unvectorize_M(W_H, M):\n",
    "    # number of elements in W_H is (n+p)*k where M is of size n x m\n",
    "    # W has the nk first elements\n",
    "    # H has the kp last elements\n",
    "    n, p = M.shape\n",
    "    k = W_H.shape[0] // (n + p)\n",
    "    W = W_H[:n * k].reshape((n, k))\n",
    "    H = W_H[n * k:].reshape((k, p))\n",
    "    return W, H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Small data to test the algorithm\n",
    "M = build_matrix_from_faces(folder=r'C:\\Users\\Antoine\\Downloads\\att_faces', minidata=True)\n",
    "def unvectorize(W_H): return unvectorize_M(W_H, M)\n",
    "k = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN4AAAD8CAYAAAAYAxqKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXuMpFd55p/TVT3dM+0Zjwdf8fgyThwDgXCRwx1CuGizbGITRCKyaBciEv+TXUgUKZisomSlRCFSsoTVroiskKyzQhBuweDlKjAQEtaLbSyDGYMNxvbYZjw2jJkZd89Md5/9o+v56tTv+15XT489VWOfRxr1VNV3Oed8Ve9z3nvKOauiouLEYmbSA6ioeDKi/vAqKiaA+sOrqJgA6g+vomICqD+8iooJoP7wKiomgPrDq6iYAB6XH15K6ZdSSt9JKd2RUrry8bhHRcXJjPRYO9BTSj1J35X0Gkl7JH1d0m/knL/9mN6oouIkRv9xuObzJd2Rc/6+JKWUPijpcknhD6/f7+fZ2dnmtYVBSmnkuJmZNYJeXl4e+dzv55xb7/V6Pd+jdWx5Lb9eXV2VJK2srIy8T3CMvO/MzExzb89tbm5u5NyjR4+O3Mt/N23aNPLaY/TxXA9f32PnXFJKzf997OHDh0de+xzO1597DIbn6+uW9yrBMZVj59rxHF7T6+m/PN7ws/Z6+Xivpz8vn8O4e/JeXnvD11pcXHww53yGxuDx+OGdK+me4vUeSS/gQSmlKyRdIa093IsuuqiZLH9YnvzmzZslSQ8++ODI+1u3bpW0toD+oszPz0uSTjvtNEnSjh07JElbtmxpjpWkH//4x5KkI0eOSJIeeeQRSdLDDz8sSVpaWhq5l8fk9/0QfT//XVhY0Pbt2yVJZ599tiRp165dI/Pbu3evJOlHP/rRyD3PP/98SdJPfvKTkTHee++9koYP3/M+66yzJA1/TP578ODBZkye1znnnCNJuuOOO0Zee0z8u3PnzpGxUpj4uocOHRpZD38x/b7HVK6X/++/vrbPKccvSaeeeqokNetqeJ18/lOe8pSR9fJ5vp6/C37me/fubT7zex6/z/UPy59biBr+Ud988813aR14PH54qeO9Fm3knK+SdJUkzc/P55WVlWaBvfD8Mh84cGDkGn74/NFI0imnnCJp+CM977zzJA0XlOziL73vTab0l92vfX3/9Y/AD6Tf7zcPzfe4++67JUlPfepTJUmnn36612Jknnfdtfbs+EXzPPft2ydpKFQ8b6+Hx2ghMzs723yxPD8LMQsSH+t7+HiP3WPkTsDrRkZcWFjoXEcf5+uX8/NnvqfhY73WHqOv7e+F7+nPvb6+nq/vHyqZsjzG35OHHnpo5N6ch9fjWPF4GFf2SDqveL1T0n2Pw30qKk5aPB6M93VJF6eUdkm6V9IbJf379Zxo6UEp6W0KYalkCbeystKSxGagBx54YOSalvxmFd/DkoxbSUvlbdu2SRqyi8dqqeptzPz8fPOZz/G2zozmsflaZh+/9hwMM6LZ1XPyazOpGdEsdvjw4Waevoev5Xl6/GRpsynHZiag3k19lM/QrDQ7OxvuaDwmw+tohvLY/NzPPffckeM9b3+fvM4ekxnP5y8vL7e2jv5scXFxZNyeH8dCxh+Hx/yHl3NeTin9J0mfldST9Hc551sf6/tUVJzMeDwYTznnT0n61HqPTymp3+83UoOWNEtNGlto9Zubm2v29YalqPUhSyrrOJauPs+SzRLdn1sC2hBiKU0r4f79+yWt6ZRmv+9+97uShuxhY4gNFtZh7r//fklD1vHnZnGzi+9tCU7mNzt7fbosxl4HGhP82vfgeTQycZdituHnHgOtgdJwzQxaDnmO7+Ex+nOzFtnHc/J5Pt7rt7Cw0DxXP0eOwbsI716MrvmsBzVypaJiAnhcGG8j6LIOWZpY+tDN4HNsrl9aWmrY44wz1lwpllCWZGSPcp8vDVnGry2NzVJmTl/PUtPvWwrv37+/0Q/Mrma0PXv2SJIuueQSSUNp6nnZ+mkJbaY0q5D5KHWtd3in8PDDDzfz9jWN0vJZzsfX8Bxs1aMv0WPyWDxXvy6tvOX9Dhw40KyHnxl3MmQyfw9o/eaYyOaGdTvrxOW6lc+tfO219/h9TbpCeK9xqIxXUTEBTA3jSW09whLJepalryWgpasl3pYtWxom8jUYHWLJbv+epaAlnF8b1pd8LzMhoyjot5qdnW0cuL63dTfP65vf/KYk6eKLL5Y0ZFsfZ6smX9NqZ3byvelzeuihh5p7ej4eryW8Gdv38Hw5T+vMlPi+55lnnimprSOZzT3WhYWFFsv42n7tMXI35OMYffTDH/5Q0lBf507BY/Kz93dhcXGxZcXma+6UaM2l73EcKuNVVEwAU8F4OWflnFuSiWFHtG7St7S6utrodGVkhDRkPp9jaeh70c9n3xBDqnyc9/TWFyyFLcUfeuih5toei+/tkCZLVfv1zBZmdsOvPTdLX7OUx2YrqtnKoWabN29uhV2Vn5Vj818zmOdr1qEV1MdZz6b10mNnOODi4uJIeJ00fEbWz8myfk290+vImFXDOwmPyfcz5ubmmmswioixvD7OO4eq41VUnESYCsZLKTX/SnjfbJ3OLGXGsBQq4+csFf2Z/XEG4/1oAbzwwgtHjjMzWNIx1pEBymVEC6UgpamPtU7jyBNKZPrK/NrHWS/z+jFSo9/vtyS44Wt4ba0f0Y/ne/p9skpkTaZl2my1srLSrKVhxva9aIn2PTwv+u38vudCJqQ13J8fPHiwOcbP1eA1fZxtAdyFrReV8SoqJoCpYDxpTTKSwZhJYB3IVj1HG1gqzc/Pt5gp8ukw745+J+tPZhlbwiyVqYcYHtPZZ5/dsKWlodnFfjv78TyG8txyDmYPz81j8LypM9K31uv1Gmb3MdZNfU+vNWM0zfweu/Utn0/dyK/JiPTJraystJidPlVaRJln5/XxjoF+XlrJzaC+rhl069atLZ+rX1vHt8WUsa2MSV0vKuNVVEwAU8F4KSVt2rSpkSbUC5jcailtC52l7759+1rZBZTklnqM+7Rkp2+QvjPqWYykL6Mqfvqnf3rkGtQDfA9mSFh/MJP5NTPYmSHh14zKP+2005p7My6S1jmvg2EWYja3GcHPyAxJHxyj+o1er9ewJ7M0/My8Y+C59M36Hn4WXi8fx3w+j9XnHzx4sJUT6M88H+u+vjZ1vmPFVPzwpLUvE8syGDTLe6EvuugiSdItt9wiae0B+mHROetr+iHTHO4F9/n8gtFZy7AthiGtrq42x3pbynAtX9tfEI/FD9f39PkURPxyU6j4fjnn5lz/yBkMzS0lf1C+NtfV86bhg2Z4v/b9H3nkkWaN/RkD1z0mz9d/o7Qh//U9vBXleQxR6/f7LbcHS0F4vlQtPN9j/QHWrWZFxQQwFYyXc9bKykrLIW5pG6UHfe1rXxu5ztatW1tBrJH5m5L4ggsukDQ0qlAqsp6HmZOGDDPA5s2bm7GQbY3S+FHe05Ke6T/cipVJpeX1GVTepfiTkcwSvidrjJiFGUrle7EOireR3rJx676wsNAwEg0z0XwMzz8ywvleTAFjOJjvv7Ky0tru2qBlpvM53kWQwekaGYfKeBUVE8BUMJ4TYS1VLKEsRah3eD/OcgYHDx5spQExwNgSnYxIfYJBwGYyS3ZLSBpdfPymTZuaa9uww1A4Fg4ii9AM7nsyKNpz9li4TimlVrnDaBdhIwJdIXZh+Np+JkzNYWCCX/P4w4cPt0IE/fyZ7sR1o4uIJQR9Ha8LjTQ0pPT7/eY75xA573zs2rELiM/Mz6YGSVdUnASYCsaTNBIkzb+Wyiz1YOemndtS21FuBrPzneURzBJ+bT3L59miyIBeMmNXYVi6Hjgfpj9RctMVwrAmWvF4nfL6TK2hhCbbml2pdzFx1sd53Rg0TLb3Oh4+fLjRA3lPOr49P1qoGQxAndHHM0DbY/czn5+fbz6jldZrzbAz6sC0WI9DZbyKiglgKhgvpaRer9cZaFyCvjMzXamvsfCqz7HkotOdAde+pwvg0mpJtjJYRKjf77f8kbROsvQdHeE+n4HcDECmTsS0mfn5+VahWoP+J5bcs27MsvNeTzKpJT8L6DLgvfRzuhQGE1+j3QtD58xsXkemLnkOtAp7Ttu3b2+uwXt4jGZnJkRzl7FeVMarqJgApoLxVldX9cgjj7SiIix9LLn9vqWMX7u+/44dOxrpRv3B1zZb+C+TKZ3mQwshkzVZ0p2vSx3PIJOxjERpES3nSV8Zk1fN4mY6n2dWX11dHdGtynt6XizNwGRb6p8MmeL1zWx8pqX+6edWlt4v58Fy77QSe8zRvJ0A6+MY/mYL7uzsbCvYOWos410Wk5Ir41VUnASYCsaT1iQgrZmWiOVeXJJ+8IMfjJzrlI3FxcWWfmA9wNc0s1ki00/F5MqoaQfH2tViiz4y6pH0T5l9eU2PiREtLPxL9imDpslEHBOtlCwZyNhFj8WS32MwI/i6DLou07V8LT9nl/ujrkZ/rVnajGZWLos7SUNLtlnLvjhbqo0f/ehHzfzNlp6HA/JZPIu7mWNFZbyKiglgahgvpdRIT6ZzmNGe+cxnShpKJTLC4cOHW341Fii1fkQdL2pESH8XC7RSHyuZMGpqSf8di72ygaIlPPVJlhunHluuC/Ui34OvaVmOfJA8z3/JcGXpxXK9FhcXW/oyU2+suzEKibsZtgpjcxKvh5mOWRwrKythojT1Q4I64XpRGa+iYgKYGsaT2omgjJcsm5NIQwtkGYfJeDxLOf81W1JvICNYwkfRFF3lFaTRKIyobW/kn2S0PcdAK2hXK+qu9UoptSypBqNkohbUtjx6/fyMGJnChFqzFpuI9nq9ljXT93RpRfoU/XlUkpCRLFHrLJZmvPvuu1sJwZ6XdVJG8rBA8gkrdpRSOi+ldF1KaXdK6daU0tsH7+9IKX0+pXT74O9pG71HRcUTFcfDeMuSfj/nfFNKaaukG1NKn5f0FklfyDm/K6V0paQrJb3j0S7k7ARbkCxd6HexVGY2dJkLZv+dJa6ZzvqApRylKPf0jJMsxyrFJc6NXq/XiuqIYjUJsid9aBwj14cFYKV2+2JmE0QlBMmeLPnAPDRaf42uDG0WM/JztZWS1QF4Tcbu0gpK35x3SDfeeKOkoeVybm6usQU4ksl/aSllefxj9d81c9/QWZJyzvfnnG8a/P+ApN2SzpV0uaSrB4ddLel1G71HRcUTFY+JjpdSulDScyVdL+msnPP90tqPM6V05nqvU5Z+k9ol+JwrVTYUlEZbTDl/isVQWafEYNnvqO0vs5nHHTczM9PylVFHpa5GvYrsQX8dM/ap05TSmIxOXY6tlBnR7+PNeJE1j/omY2ZLvd0M5TX1Gpp92EiSeXu2ejO6xjo/W4u5/ZmLFu/evbuZk8frRjP0qdKvy/kfq1XzuH94KaVTJH1U0u/mnH+yXiUzpXSFpCukjVdqqqg4WXFcP7yU0qzWfnTvzzl/bPD23pTSOQO2O0fSA13n5pyvknSVJJ1yyil5y5YtjZWL+VaWOo48YP6VXw/GJClmBYMRGvRrsYQ5/YMGpXEpCcsolvLaZIGoESMZjDoudSSyb5e+FVn6jEjn5c6A+mbEfBxLWfeEpeb9HL0riYrk0tLM1mpcJ5Z6ZFHesuQhi+NyVxWVtj/WlszHY9VMkt4naXfO+b8VH31C0psH/3+zpGs2eo+KiicqjofxXiLpP0j6Zkrp5sF7fyjpXZI+lFJ6q6S7Jf3aei5W6iOs40EpS6ugdYGjR4+2/EisT0Ldhn68SCchq5D5GI3S9RmtkVGpcrJO1OSR60DLZXl/VuAyyMbU7ajjcC6s3UK9nLpQqadS/+Z6GSyXTh+bI1M8f1vHWWaejFnmWLLEvVmU7Mv5bjRyZcM/vJzzVyVFCt2rNnrdioonA6YicmVmZkabN29uNXekj8n+GWYYlBZMSzlamxjvR8ltsPYlWSbKseuqhOX/ex68BvVDWi8N6o+Gj6MviXGFvV6vZa2k9ZH6IBneiOrBRGzr85nXVsayRr4wvs+x25pp3Z96GBuqsH5KmefHqCHWE/U8vY6+B62860WN1ayomACmgvFWV1e1uLjYii5n5SdLHet0tFqdccYZrZ4AtBRSYlPPYMYD9anIisWm9L1er9VAsYuJusYSjTFqE8zoEkbTHD16tDU/jsFrGDG736eOaERjYIxryaRRpgTbcXGtWU+TbMO2ZvSx0soptXc6bIZKWwF3CLV3QkXFSYCpYDxpTYKw7a0lu/fylFSsObJly5ZWpa5xkfyMQGEFMIKda4yuqse8d+TX45go+QlaDBlnGXVdKo8xWPks8stFY6K1lNZjrlO5PmRH7ibKiKRyTNSnOVbrZf5euJI0OyCVlasd7ULfqttjO67TY7ZFdaMxm5XxKiomgKlhvJRSS6qymhij+2m9KuP/GJnCqlGsSRkxIyM4aAWlblfqY2QPnhvVPaH05PvMVCdTcA7luTyHPlIj8s9F/eDItrS4khFLq68RZZ/QnxuxKqtX00ru75FZ7L777hu5fnlvWi/tzzMzOiOGY1svKuNVVEwAU8F4KysrOnDgQMvPYgnGDqHebztnyu8//PDDoR+J0pK6ClmJ7EtpyixxWh5XV1dbrEAW4T2jvnbUCdl3L4qQL/UVMjcZKbJCRrqeQT2Mkp/+w5Jhox0AcwZ5HKNIvD7shej1s7+OES9mvqWlpZaF3Oe6yjWZL7KYrhdT8cOTRhfXRW/uvPNOSW1XgBMavfBetM2bNzcLQ6NJ9KXktrbLLVAiCiEzyi9y1LRk3I88cqzzIUfhWfxx5aLNdVSQ1+C1/COn8YVbVaYmccvZFTxA1wSfUbTN57OLCgP7B+exsnyFvzeHDh1q5mlhbqe8x+Qtpg020bNaL+pWs6JiApgaxlteXm6kSJmmIQ3dCVSmmXIitRtWeIsQOY75mkxGVwDfN1gWL+fcChWjY5xbS0p8bms5XxonGABeGlv4WdfadY2JbMOg6XGJoFHqUtcx3PYSTNalmyF6pp6rd0M2uvj57Nu3r5Xoa1Y080VqQLQdHofKeBUVE8BUMJ6bljjUxw7QrlLk0lCxNawTli11GSpmCcaEVkpbn0/9iRKbf7uSUWnmp95JfYhsxPPpdqCORHbpCv+KknB574itx5WbiIxbDHdbWloKxxuF8TEtjLsbBs/7td0vPt7GOBcwWlxc7NSLy/WgTkzmP1ZUxquomACmgvFSSpqdnW2VArC0YTNF78NZDGdubq5VDt1WKl+bxVDZXiqyTo1L8e9ikqj0HyV3FBrGMgNs8xUFMjMx9tEsbtRVIvM4JTstqNRxIotuqSONK1/IBjR0M5Cd/WxZAt6gHse0ofLaUWAESwZutHlJZbyKiglgKhiv1+vp1FNPbQXcspQaLUpdBW/ZUJLOd7IMywEwLSiyuNGPxQJFMzMzrfCsqBwC3zeop0bBwbSoUn/t9/uhbkffF/XGKHk38tdFOwKOeXl5uTN9qTyGfjs/S+q43tW4fF+XPikN04UMWzfLa3oto+BnsypR/XgVFScBpoLxUkratGlTI6ncUJAhYy7HbvZiMO2mTZtajMa2WtRNosgDhlRFkpzSlwmVXccyLCkaEy2lkU5nSU5Ws3Q+cuRIi12iCBRGuBiMnqGvkIxIn2sUZSO1dzIem0O7eA/qWdbZzHwudmQwaNrH2xp+5MiRRi/k82NUDN83aumHioqTAFPDePPz843UYfl1SzhLQEsbM2BpsYzaZ0UtsIxxuklUWo8+p5IJonQgRtNElrFIv2IUjtfDDEedcnZ2trEMU3f1mtOPxXlT56FP0ogCv7t8iixJH0XeRP5M6sAs3e65ke19nNHv95trOt3HPj6uaXlOOYeoJEaEyngVFRPA1DBeSqmJHqd1y5ZK7uUtRUurld9jQVsWPxpXgHRcs3n60rrSX8h01G+ieMmoIQabZ9KqGxVL8o6i/IyRGEw1YtGoqIgRdb8ofYr3W11dbbEm9U6mVUVrzsJU9NXyO+DGl04vO3jwYHNNxgkb1Cu5M6hWzYqKkwBTwXgu70dpYqbbuXOnpGF+nvfoXZY2WqXGFQwyxkWXU7eLdKFSAlKCR0mzkTWTxVWpnxlkM//1DqHL0jqupMW4cn5dTVq6XkfnlZE91OWjgkvR2A1bd/394I7IRY+e8YxnSJL27t0rac1W4LXifKnDRQW3jrWEe2W8iooJYCoYL+es5eXlJrv3rrvukjRsIGhJxkK3UfxleY6vSatjV/Gd8vOukuzlvaPjSoaNjomsdFHpPTIWYxcZ2cH20eU9aTmMokU4BoPlKCL/nhGxeZkj6PdYkr+rZGIXmCXvnZLHZKbzd+L73/++pOGz3L9/f3OOrb/c0XDnxN0Md1DjUBmvomICmArGW11d1cGDB1sMZl8KGc86Dcu2nXrqqc01GOdJK10UBcHydWS6rphMqduiGOXq0SoX5QYaUf2TcTVJjJRSa770nZllWGae+imjSaIxMOa1S68dl4dHfcpjK5tbluezwaXPc8TKeeedJ0m64447Ro733MvxUpf3dzDS5U545EpKqZdS+kZK6drB610ppetTSrenlP4xpdTe91RUPMnxWDDe2yXtluTQ77+Q9O6c8wdTSn8j6a2S3rueC9myZFZy5IEZzTUzWGzUUvjIkSPasWOHpKGEYjYBs5ajKlhRZEtk5fP5ZXT7uEK2kU5Hvx3PY1k7si/bm5W5b7T08Z4G9S++H9Vqiax91EtnZ2dDXyn/RjsHxmoajNP1Dop5eOUzY9k++4FpV/D7vhZtB+vFcTFeSmmnpH8n6W8Hr5OkV0r6yOCQqyW97njuUVHxRMTxMt5fS/oDSVsHr58iaX/O2aHbeySdu66BFDlj3MM/8MADkoZ7d0sfR7qYZY4ePRrWXKTVKSqXTv2LelkUTUG9Kufc0msYfc9rRCxC6RrpilHVseXl5bA2CC2vBtePkSvjmmUaj2YV5rOKrLJRfZfIn0fdMcqtLK8f6fBRxBJ13xMWuZJS+mVJD+Scbyzf7ji0k4NTSleklG5IKd1wrKbYioqTHcfDeC+RdFlK6bWS5rWm4/21pO0ppf6A9XZKuq/r5JzzVZKukqQtW7bkfr/f0qMYq2n9zflWzK8699xzR3LQukAms55I/cHRD/QNGV01RMrrHD16tMUaUawh4xrHVe6K9K+o1src3FzoUxwXsRNVDePcoigdrne5Q6Bu52PL9lklWC1sXK0Wf259/5577pE03CmVFej8nXIcsL9T0Xx5z2OtvbJhxss5vzPnvDPnfKGkN0r6Ys75TZKuk/SGwWFvlnTNRu9RUfFExePhx3uHpA+mlP5U0jckvW/cCXnQlNISyBYmWzOdKez3LSkdZWBL06FDh0LdLaoexXhBSzoznlmW7abYPJOWtSNHjoSZ5ZbYluBkBbISM+0jto2qO6+srITtjqPcPn7OWEXmTNqyzAgWWhhLlqNVmow9ztrJHRJzBs10jvH1M/PYzHKl/9efeefEGq5RLPCx5uM9Jj+8nPOXJH1p8P/vS3r+Y3HdioonKqYickUajd2jz8RWTcdumuns5zMjltKJkjxq08VYO0pfdpphNSrrmR5zmddG/cDs6XEyKsKg5DdYGdmI2imXFalpfWQ1rcgX6PVwDwGvPZmAETxRha9yh+EdDCu8+Z5G5L/z/LhevrezyX1vR6x4d1Lqt45qsR7IHD4+X0bXnDAdr6KiYuOYCsZLKanf77eYznpWJPnZm02KMwCi2EFemzGeBvUnM4ZZlzpgzrlhJp/j2FPPj1XUIguiwRbM1KPox/NYjh492speZwQLx0CdjpZW+lKjfgeMovH6z83NteJEqbv5fX5upmR/B/rUrKdxh0Br6LZt21o1VxjBwp3NuOyFcaiMV1ExAUwF40nDnDypzTqWRpS+tDCVHXqoo1Gn8z4/qt9hKWlGo7XPUrOMmin/lvOhpLYOQ+slK0aTnSKrZ1StzOxTWnsj3ZedmKKq1b6H9SavE628ZhvPyfcv2yL7XFqi6Yv1+9whjMv9c6ST5+ZKBrfeequkIZullBoLOhnMdVl+6qd+StJQT6TFuubjVVScBJgaxrOeJw2ljiWhJZ1rZDAPrytuMNKXuEePMs3Nupaalugc2znnnDNyH4/hkUceaUlBWsRoOSNL+32zhlnakjrKbuiq3UI/E6PqPZYyN608zl15o5hP6ttd/Ruk0Vw5Zg+wWpjX3Czp16xGx1hM/vU9rb/ZCu65HDp0qFlLz5+7Kd+T9oeNdg2amh+eNBy8v2BPfepTJQ1NvH6Y/tzbQC/G6upq8zDYYJI/NH+5/VAYQsaGlpFxhQ5m3/fUU09txhltQ1hykKXXmZTprarXwesVbb3KsCZuR+lIN1wqw/PjD40GD26xmR5FtcHr2ev1WltDugdoTDKidCg+a6+T3VB2pHt9u4Qkgw9+7ud+buTevofLSWwUdatZUTEBTBXj7dq1S9KwaYmZ7t5775U0lITeYtJ83Ov1WkYEJo9aytk4YGlp57alLMPTaB5nOlDXdsjuA27fzHRmFzbY5Na5y4hUzo2NUxgskFJqGSjIFmRAJoDanG5HuncbNj4wCODGG9eSVrxrcRBE6U4w81x00UWShttZByqzqSjZl4zIbZ/n6lL/LKZ19tlnN9cyg/nYV77ylZKkF73oRZKGriB/9z7wgQ9Ikn74wx9Kaj/jcaiMV1ExAUwF423btk2vec1rGslnaXvBBRdIGko8S1ezjVnJr3fs2BE2tLceaFB/YngRQ6IYpuYAW0s8l4yzFD948GBzrCWyDTFmHyr7lrYsTeA50XFMfZXsXir8UbkErhd1Ns/TLEE98vzzz5ckPfe5z5Uk3X///SN/XTz2ZS97mSTp3e9+tyTp4osv1ne+8x1JQ8azDnb99ddLkp7znOeMfG6m9xpTjyTz0f3i9TZKPdZl3V/84hdLGroPWFjLz/Atb3mLJOnP/uzPJNWQsYqKkwJTwXhbt27Vy1/+8tY+maFCX/7yl0c+ZyDu5s2bW+k+/mspyRJwTLa0c9uuix/84AeSho5TFkU13Ab4bW+wJlUZAAAgAElEQVR7myTpwx/+cKNDmC1+/ud/XtKQDXyP++5byxW2HmUrHCW4x8agYoIl38twPIbIkfkYhuZ7Wnf1enpu1pu4UzDL3HzzzZKGbGO97DnPeU7zfJ/1rGeNHPNP//RPkobM553D05/+9JExUC+Nmo36c4/VNgTP6cwzz9Qll1wycg/r/AwW9/fnmc98piTppS99qSTp85//vI4FlfEqKiaAqWC8mZkZLSwsNNLSbGSrlvUz6zj/+q//KmmoE91+++2S1iQafUF+benPtCD6wqyzWdpaKl988cWSpGc/+9mShkz5L//yLyNj+fM//3NJa5LT12SRXes9lsi23lriewxOVWGTFrO3rZnjmiT2er3m3gxX8xjZGovlNRxu5ff9jHhdj9k6kvVyX9dM+ZrXvKbRl3xtWxb92tf2/L1O1vn4jMnmLKVhpvN5nsNZZ52ln/3Zn5U0/K55Tb078b3Nzrfddpsk6YorrpAkfeYzn9GxoDJeRcUEMBWMt7i4qG9+85uNz8e+Iu+3md7hv2aG733ve837TBGJoiG8v/fnlmjW1fz65S9/uaSh3mUp+drXvlbSkM12794tadSiaBax9DejWTLbmkk99Lvf/a6koVXXuol1HYNz9ZyY3Ds/P98KS7M/imXMPTYznd9nKo2ZwHNikeHLLrtMkvSVr3xF0vBZnXXWWc0crVfZT2td1zsbMxSDoz0mJtIaTAnz2DxWr6dZ7bzzzmv5DD0Wn+tn5Wt6/XzNF7zgBZKka65ZX4mhyngVFRPAVDDe6uqqDh061Egq60u2BjqixdLFUtVS6qabbpK0xkaM+6Mu588Z7MrgWEcsWLJ5LK973VphbPv1vOd/+9vfLmno51leXm6kpPUeM7n/0q9k9rFFjQwWNTmJWoeVOrP/73v6mizlYFY1y3gXwoK+ZsQy9UgaMr/ZyqxkP5mf2b59+8LCUQZfe0x+VmXcZzkWg8nOXi/vpPz5rl27mnXxd4kWUF+brHzLLbdIGlqzK+NVVEwxpobxlpaWWuXcLHWsJ1jfMltZCr3kJS+RNBrLF7WsooWRqUWMo3SMnqMs3vOe90gaRjhYAlpfMYNu3769FVnCsnOer614TEliiXrrp0zqjSJWSn+XmY0pSdYrzXQek/17vqaZnM0jzSa8nnVA+tLK5GaPxbqsWbTMGijnzTGZsaIyDPRJGoxG2b59e2Ot9fM2Q9MK7mdnH6z/+ju4XlTGq6iYAKaG8cqSbizy87SnPU2S9NnPflaSdPnll498XhaooUWLia4sX8eCOZZoluy2ZjJJ9atf/erIeSxbsLCw0CoFyGRaWlypo/GePI6tsKL2wCmlVvFXroulP9fPTGi91LqNj2ecpI+nj9HHOQOh3+83OhqjiHwNNoOkfkowWyVq3czcwIWFhSbChqVB6DP1bstxpdaNzfDrRWW8iooJYCoYT+qOrmBOnKWOfW2OjDfKlsMs0WApaNZhxIalK3WgMrZQGuohLLzE7OeZmZmG6cpYUo9Tamezs00VLbJ8zYYjbNlc+vm8HmxvzCgOr4fHwDw+7hwiXZpjNkpdkHl24xrO8BpR5A7jTn2cnz3HtLS01NzbkTUsr+ExMSPG93L01HpRGa+iYgKYCsbr9/vasWNHI5m8l/drsxEzDLraVbGpIQu0RrodfWPMRGf8qK2hZkBaA8smmSxgy8JCLCnOmiC0xHKstJ5Sj92/f38r345FiaIWYp6v51mW55OGepnPZwMRsk9ZoIjZJ/RTsswhWYbW3Ih1WTXA8Dp/7Wtf0/Ofv9buw9Zs+kIJ+/eYp7leVMarqJgApoLxZmdndfbZZzeswggEM4alrP16tiSx6pbU9vGwhS73/9SbGJtIa6DByPhSryIj07LIa5glWIOFPrCupiTlGBjR0mUxjix+zPXzcYzJ9DUZL0o2jnLkDh061Pzf86UvlQWOo2LFkTWT62H4eo4RveWWW5p8S0e1OFrK8/GY/AxYueyEZqCnlLanlD6SUrotpbQ7pfSilNKOlNLnU0q3D/6edjz3qKh4IuJ4Ge89kj6Tc35DSmmTpC2S/lDSF3LO70opXSnpSq01qwzR6/W0bdu2RjoyXtBZ4LYcWSey7+XVr3712mT6/RYjUS+gPhFFP/h9WjFZpcvoanPFtlmUxAYZ3WDdTepNjIjpqi5m0KpJiyj1UV6TuYDUnWlNHrfeZVY81y5qjkkWJcuw1opBv5/XucwV/Jmf+RlJw8ijb3/725KGETvW/XxvZjg4vni92DDjpZS2SXq5Bh1fc85Hcs77JV0u6erBYVdLet1G71FR8UTF8TDeRZL2Sfr7lNKzJd0o6e2Szso53y9JOef7U0pnrudivV6vkZ6WOo6DY7suSyH7ZZy5bAvk4N4jfynhqReQLWidY2wi9Qxe/5RTTmk1f2RVLOpFtN5F7bc4Zlokqb+W5dJ5DYLsyuYuHCt1Qq4j9fXSh0adlc+Kuw8eR+sn9UiCWR6Ox7zpppsa662/Q2ZDW68N1jT1HFjKfxyOR8frS3qepPfmnJ8r6ZDWtpXrQkrpipTSDSmlG47VFFtRcbLjeBhvj6Q9OefrB68/orUf3t6U0jkDtjtH0gNdJ+ecr5J0lSRdcskleWFhoZGu/uvqw4ajCixNnStWsg7ZIGI6NviImJHRJPyc+oeRUmr563xvsk+UXxf5pei/onWvK07V82NNSl7LY46stQaZkX5CWpGZpVDqyGQu9mFgtoHHRv2d68VdT9TW7KKLLmrmaeutd1McE5+RdWI+03HYMOPlnH8o6Z6U0iWDt14l6duSPiHpzYP33ixpfZmBFRVPIhyvVfM/S3r/wKL5fUm/qbUf84dSSm+VdLekX1vPhVJKjRRhNInh2EznxnmPboZcXV1t7e9ZUZrREaxdyQwBn2eJSL2LPidKyK75WGKT2chYBhksstSyUlip01DaG4ypZE6g58eanozljKyYZLpyXXjsuJZXUYwqmY3W0Mh66s/PP//8Jv6X1ei4PmwdZ2smdcFxOK4fXs75ZkmXdnz0quO5bkXFEx1TEbnitsXRXt3SxlnClprW+RyBUF6DESjMVWPVLcZRGqzKzKgTVmgu4wEjCd3VSLNrvnwdxRwy/47xmP1+v8UG9M+Ny3iwLhP1kYsqeFP3K31wUWYD9cfIahnp57Q0c6fQ5Tf0e+wLyOwN18Pxd9GvXatlvaixmhUVE8BUMN7KyooOHDjQSA9LKPtU2FnVr32862QcOnSokazsDkSdjHt1VrSifunr0lJpKUrdYGlpqXVP6k+Gr8Exc75kW1+f3Vl5/dLiRmYiSxB8FmQTWj3X23l2ZmYm9MMR1H1pYaT/koh8umVlAEYVmfFc9duRKszts2XdkS7rRWW8iooJYCoYL+esxcXFllSxlHTkuKWRGcBxlM7XW1xcDGuuRF1UI/9SFCVh8DUl/sLCQkv6G9RNzCasfs3oCepf1DtteeyK7YwqkVFnY5SNr8l40kh/iqJneFy5frTCGowuIiK2JgvzO8Bq0KeccoouvfTSkXP913YEn0PdzlXmbr311s4xRpiKH97s7Kx27typPXv2SBo2pfDk/GVwaJgf1LXXXitpWC583759ra0CzeL8YRE0RRv8IdOIwC9wSikMiWJZCia82mnv+TIZk0YoNtdk2lAJpjlFbhCPjQ53/wCjUoJODPV6M/SuXBP+MGgYM6IfWLT1pkvIYHCFvytbt25tlXrwVtKgQ91r/clPflLSsET9elG3mhUVE8BUMF6v19PWrVsbGnfTRzvMLY0sdZyk6JAyN/nYtGlTmDRJdjCiVsSRM9vgFooOZamdiGlGixiLQeI+n6xBljFoyPAWtmQQbiW7HNvl2LyN99ozkIAMyCReomQvBlxHwQxR2UOyN8+3KhIlQXusZ5xxRmgUsuHO83PTFaequeTkCQsZq6io2DimgvH279+vT37yk/rVX/1VScMWWTaqWF+w9PVrNxa54YYbJK0pwtSnaGKmTsbwLRpbmPbCAF1K39IRHwUSRyzLvwZdFda7aBDh9cvwOIZPkfmiorrUlW3YYsMQH+exeixROcEuHZvhWWSy6PiuML3ynjS4eQ5lKwA617nbsE7ne37605+WNNTt6Aoah8p4FRUTwFQw3t69e/VXf/VX+oVf+AVJw+Bnh4KxiYlh5nN75LvuuqvZt7PUgCWX9R5aN2km72KN8n06cS1FS8lJJjJ4bZr66eilPhW5SpgmVCIqbUErLsPOqFex3TPdKpxjFMidUmqVFozSfKL3x5X54Pt8Hj5/bm6uxfRsFOOxfu5zn5MkffGLX5TUfibrRWW8iooJYCoYr9/v6/TTT9cf/dEfSZLe9a53SWonZUYhUeVen/oTEzGtczDwOEpJoY5HyU2pXI6N4VW0rvE4MpalLRNooyK0htnXKSsLCwtNKUQWz+UYOOZIN6T1N2qNZX20qz0y9W+yJ+85rqRD5AeMSjT6uEOHDjX3YtsAv//lL39Z0pDpuObR2CJUxquomACmgvEc5WHfyDvesVYN8Ld+67ckDQNRLfFZ/q+0qNHCFZUBJ1tEZcTpG4r8gT7femcZJE2dIgpro07HkntRypLPY+SKI1+Wl5dbraZ9D/tCGSzOMTCUyvDxEYNyfcoiSdQfo1IXZEayC8t5RMWg+B1wgaO5ublW8Pv3vvc9SdL1169VNrGl3WF8fP4s9zgOlfEqKiaAqWG8ubm5Rno6Df9P/uRPJElvfOMbJUm//uu/LqnNFJbiO3bsaKI+omiFKAia1kwW+zED0rfE5h9lDGNU8Ifj5vuGme7000+XNNRP3SSSljhLcEdX+DjHT5bjdCQKI1I8fjfloBWT5Q2juEpaGLvSr6hXU49mEDgtylEakBEFT7OV9+zsbBPl8q1vfUvS0HppPx3bN7OcxLGiMl5FxQQwNYw3OzvbSGZLJEu0D33oQ5KGEutXfuVXJA3b+pqdzjjjDP3zP/+zpHajRTYWpJRlSYioMQhjNGkF7CpfF+katAwa1mE9Px932223jfy17usGGvZ7OpXFFsXSgmjrppmOxYOtZ3t3wRbNjqd1BgnTrKJGIV1lLuhDZGGqKOqIOwPuYqImLrSKW2+77777mnhfr6H146hIVDSX9aIyXkXFBDAVjJdz1tLSUrOnZ9yjpcnHP/5xSdLXv/51ScNmJS972cskrek4tPxFbbqigkJROcCoaA+tn6We0uW7Kl/TmufX1jccg+ry4PbB7dy5s/N9j8WsVEbMey2jgr9mSa/T7t27JQ31RevdjJelrkdLZGTZ7ff7LQvpestRGFxP5tlFRaU+9rGPSRqy/crKSrPbos/P4E6GVvBq1ayoOAkwFYxnWNpY/6AFyZEc3oe///3vlzTMidq2bVuj97ziFa8YuYZBCyNzwqJmiBxjVFa8vE/UTirKx/O1HKtqlnGJAs+F2dKGpTavPzMz08zLUt4S3dc0S9iC6mdgy6h9hGZA5hZyfTj3sgWz14fsGDXe9FisszI73Ogqoy8NnzHzN8t1Z0QOI5KikorHqts1Y93QWRUVFceFqWC8nLOOHj3aav3EfTOLyzrKwBaoBx98ULfccosk6Rd/8RdHrsH8M0tPg1bJKP+Ox0e6Y8655WfjvKjj8Fq2zNoSG+XddbWiLj9fXl5umM7w7oFZGSzA5OMMrxv18aiGDedm9Hq90ErJa9JHSlDHYwSL5+hoFPtDvXMo2wdE2RRR4SpGOq0XlfEqKiaAqWA8aXSvzFwvShWyWCltbdFzxMEFF1wgqd1uixI4ioJgxApfMzq/KyeO16ZUtW5mNqFFzWNnMV5a7fw+fZErKyvNMfYNMgLFYDQJC/ZG5fw4lii7oVw3MhqbYFJPjHYhtFwzBtR+UVex47Pq9/thrGlUDtEYx/QRKuNVVEwAU8F4KSX1er1G/2K0wLhirMby8nKjF33pS1+SJP32b//2yLFRlTEiaiFMC1w5hxK9Xq8leRk3ymsz0oJsxBwx6yjUFQ3rZysrKy1/Gv1PUS1K1nsh1ps9bjYv1ylq6hndgzsBg7od9VRbMV18lvV0lpeXWzuaaLdFPy5Zd704LsZLKf1eSunWlNK3UkofSCnNp5R2pZSuTyndnlL6x7TWO6+ioqLAhhkvpXSupLdJekbOeTGl9CFJb5T0Wknvzjl/MKX0N5LeKum94663uro6dr8cSZVSB7IUc0lt6hpRFnsUZW/QKhplN3TlkJWNTB5tPlEWg+HXkW8yssStrq4247c+GVVvjqJqGKvqdYjy8QyuS3k/rwcbdnK+3Akw6oY7AWYhuJEpGbbcMbFiduSfG8eM68Xx6nh9SZtTSn1JWyTdL+mVWuuHLklXS3rdcd6jouIJhw0zXs753pTSX2qt3fKipM9JulHS/pyzf/57JJ27nuv1er2W1KBeZVBal/qEpaCjW2zldIWycblczLyOPo/qp5SWOrbyooRnLCrvYTYhK0VVyfzXem7XNVmRjNdwhArZw+c7aiTKV4yyNrg+OefWOZFlkfOOfGqGx2Z2d31W+mrL7wAjedysJWoXzWcQ6dkRNsx4KaXTJF0uaZekp0pakPRvOw7t/IanlK5IKd2QUrphXDBsRcUTDcdj1Xy1pDtzzvskKaX0MUkvlrQ9pdQfsN5OSfd1nZxzvkrSVZK0efPmXPqaIktR1C65jNCwjue/N998syQ1NTstwSzRGPVB6yel67gI+tJySSZif4VxOWucJ32HXB9/br9VaZklazKrgv0ZDPvWGH1j9o6iTsiItG52jYlgTRnqi7RmkpXuuusuScNdD8daslYULRXVMKVv9URWGbtb0gtTSlvS2gxeJenbkq6T9IbBMW+WdM1x3KOi4gmJ49Hxrk8pfUTSTZKWJX1Dawz2fyR9MKX0p4P33ne8g4wqJjN7QRpKSUdo2KJlxqN1Loq1i3S3KDua7FVWkqYljNdmleooYiOqH0NWot41MzPTqh5GBvPacgcQxbgyEshjYJ0UZid01dWMjuGOIeoORLbxmO+8805Jbb20y2rMe5OxGS3TVWngWHBcDvSc8x9L+mO8/X1Jzz+e61ZUPNExFZEr0poEoZ/LiOr9d4F60U033SSp7Z/hcVG+1biIDUrrMsKB2RQ+h+wT1QClb43HRbUdPXbPuewWZMaKWIJ6ksEeAWRIWnmjamOlZZdryMxxg/66qKYKq4W7zyJ9iHwuR44caa45rnI4z2XO4HoxNT+8krK5JeAXKkrDl9pbIL922NCznvUsSXH5v3HjixpfdDXMYMOLqDgu70EnPH/c3HrafcCHX/4oOBaGPkVJpNEWi64cnhcZTMqe6r6Wn1FX2F15byMqM2+hYveBXUr8EUVBD1J768jnThdIVNxpHGqQdEXFBDA1jFcmRjIlg0YVbv/K5oqUPN4ifeUrX5EkPe95z5MUG0co4boKsZbH8Xz/nZ+fb223okBrGjzMAD7O7gHPk41HWGy2awtqZ3IUeM1WYGQfJsyut7gRUa4rx8JAbe5GoqRbGnJcwtDJvw6WNhj2VpZ+YHhfFJ5njJtvhMp4FRUTwNQwXln8hiDrsHBOaWygdHSomJtP2Jnqoj1sRWx9iU7uyOltUOebnZ1ttSF22T4GalNBd4Eh6mxeB+osNOmTtcu1pY5KtjRLeG2tG5oJrUfRIBIxYcSMpUHDYGCAX/sZ8V5+7Xl7btYjGdYVlf5fWlpq2Q38XaNRjjsj41ijryrjVVRMAFPDeGUpgKgVFovLPlriKFmBut7rX//6znPpFoja/0YWyS53giWy2daMRhaypLZOx7bAvo4/9/k7duwYOZ7m+Pn5+Ya5yMK+tkHmNpuaHdzMxDsDrosZ0X9Zjr8rQJm6ncfAIHNaM6mPkYX5PWLpCJ+3sLDQ2l1ELiuypxHt1iJUxquomACmhvGk2HlL31xk3Sz32dQxXHr8mmvWQkcvv/xySW0HsME0kSgYmC2lPIalpaVGXzLT2a9EKyUtZ5bYUSAyfWy+HlN9zKCSWowXldXwWDgfNj0hE/r6bhXmv76u040ejemikh9R+J51ZrMV2yhHOmRXaf2obB8tzlFpi2rVrKg4CTAVjJdz1vLycsM+DKWydOW+mp+7MG6JKMrl05/+tCTpsssuGzmekqwcY3l+xEJGabXzX5eXdwC32YdpLVGjRkt2hozR7+fPrRNu2rSpJf3NhrwGk3CpA5Lp/MzY4tnw+dTPZ2ZmQqsjg585Br8me9OnFrFVFA5WjoUB+FFQdC3hXlFxEmEqGC+lNMIYlL4suU0JVloiyUS0bFkP+OhHPypp2NzEEntcGXnqAlE8Ya/Xa9jAvsTIB+aoEussZAeDTSL9Ocv8sSjSzMxMY/k0a5ItOE82M6GV1+voOXl9be1k8aBHYwr677yWPpfWTs8hKoAbNVJ5tMD3KGKJUS5RatGJTIStqKjYIKaC8aTRqJMoioBMx1ZaZQp/VCqP/qtrr71WkvSmN71JUjuFh/65KB7SKFuKsb2zwSJGtgCa8WyltGRnOlE0N8Ykeo6HDx8OE19pCTUbsGyf//oe/pw6n6/j4zimMtqI1keOhVkpXgfrp9QNae2lvs/WYuWzYkpa5AOMSi7WVswVFScBpobxynJvjK6ICs50lYWLolu4V3es5qc+9SlJQ+smrXSM5KDe9WgFi+y78rnW5Qiyh32OPp66W6SnUc/09fr9frMuZEkWZvVY2J6L/j0mwpLdGZUTJciWY/E8OV/GqHr+ZFmW76Almvcrv1++thmaTVZ8Lz8T69W0uK4XlfEqKiaAqWG8lFJYuJb6SJQ9vbq62mqtbGnpcyklLenuueceSUOJzfIMHhslPPf8ZbtgZg+wlF6kR7LpBhuHmI0Ytc/rlWUWeM8oG9734Fg4Bj4Tsg11w6gJp9SuNEDG81jdkprPnSUL+R2Iyk/4WR04cKBlWY0s51FpiBqrWVFxEmBqGE+KY+vsr2Hjxi7dj1Kf/ihKaBaxefDBB5trlfekxI9i9cqseVpAfY1ovlFTDUpT6oRRFL8/n5+fD6/NeEjuBMhwUTn5qIR7xF6OVpKGz5fFnPh5lBVOXyoz+6OGpl11dvheWTCqHJt3HV3Fg9eDyngVFRPA1DBel+ffks1R/uMiFMr3WIiV7GkJZUajb9D6RGQZ83mPFsvHEuORha8sM1fei8VhmadG1iFjlOxPXY4FW6OWzAb9mzw+KkVIVi91TVsnqbv7L/MQWR6e5fq4S2Esb2QVn52dDQvTssoa9cQaq1lRcRJhahhvZmamtUdnvhr1lC4pRGnIOi30T1lyPfTQQ5KkCy+8UFK77LfPo3UzitFbWlpqZQvQ2hbVkeRr6rJRiy0fF1Uhk9pWzKiUfdQghO2lOX+vG62+Xfo7W2/TqunPHanimFezL/Uvg7mUvHdXLczoGflYxp76XJa0Xy8q41VUTABTw3ilD47SknGS9DmVFjiyCCNZ6K9hhApLeVuCW3pawkUxm6UuRP0yasJB/SLKsCBLkU14XMmY1PuihpI8N6qDE+X3GcyJY+3LQ4cOtRic84lyBKPmLszn9PmMtukq/c5xeodDK2+0PjUDvaLiJMDUMF4pMcgQrCjNqPOS1RjvSclECRZVYSYj2MrJmE1GQZR+PdawjOI8qQvyHmQX+pqiFsWlxTXyz/FeUTa4QWYwmN3A6xmMMy2vxUra/uuMfVaCY3VrI6oOQP20yyoe6WxRn4do1zIOYxkvpfR3KaUHUkrfKt7bkVL6fErp9sHf0wbvp5TSf08p3ZFSuiWl9LxjGk1FxZME62G8/yXpf0j6h+K9KyV9Ief8rpTSlYPX79BaD/SLB/9eIOm9g79j0dUInjofpXbXvpqSh5nDkeWQsXhRHU1LbFq1WMV58+bNLeusEV2bDEhdlyxMkDmNo0ePtu4ZZVuwxidjW6PKYFG9FI69HFukJzoaxLGqbJbJMTFShTpjxHBdvleeyxqcUWU7WsvHYSzj5Zy/IulHePtySVcP/n+1pNcV7/9DXsP/1Vo/9HOOaUQVFU8CbFTHOyvnfL8k5ZzvTymdOXj/XEn3FMftGbx3Py+QUrpC0hXSmkQrJW9kmaRPjlEUXRkOUc17dofhNS1lWbPRzOZoGmZ/l4xA31jUBSeKHokyIJifRj21qx5IxJrUaamPRrGtXZkh5fGRLl3OgZ+xT4Hrt5R5hSWo40U9DI2o+9Tc3Fzr3qz8Ni5CZdK9E7psqp0dH3POV+WcL805X3qsKRUVFSc7Nsp4e1NK5wzY7hxJDwze3yPpvOK4nZLuW88Fc85htD4jwZ392+X/omSm/4kSzHpFFOHue1uiO6t87969ktrWzjK+0ixIyye727D+I5mRdf2pZ5HpqI+trKy0rLbMYCDrRvVveDwjX6JIFWbPHz16tJU3aIZi91ruVsyE3DlwrFwHRvKU1uCuinXlMXw2jPc8VmyU8T4h6c2D/79Z0jXF+/9xYN18oaSHvSWtqKgYYizjpZQ+IOkVkk5PKe2R9MeS3iXpQymlt0q6W9KvDQ7/lKTXSrpD0iOSfnO9Ayn9XlEfNEYuUNps2rSp5Yey9GSdDlf0ouSnXkArp48//fTTJUkPPLBG9swRLCU+LV6UuIy6L9dEanfFoUUxyrQus74jxidbmo05j8ifF7EK2dg7i1I/jfQnWjMZsUP9ks+2q3tUOVY/D+ZoluP2MXxWzCTZKPON/eHlnH8j+OhVHcdmSb9zTCMYoDSMRK2iaNLlF/DIkSNh8G95H2m4baUboSs5shwDtzFnn322pOEPsCypwIBag+UCfBwbbRoskUeXhhGVKyh/NJFRhClMTADmNSPhR4OJf3As9LSystIK12PBJQa6c74s58dtLgUS1ZByK+41ZlFhFnGKttzHihoyVlExAUxNyNjq6mqYFhQ18+AWpTSbM8THWwe3zGJKiT/3tS3xud3jdX1+V0k+H+uUIxtmmFzJ4rF0kHtMTouJCukx6nYAAAcoSURBVMGSEUr27woMLv+yqBEZj+vCnYTnzfXxWL2F9fXm5uZaRXK5E/Az4XOn+4DPnEHoDEzoWi8GVEcJzvzOMURuvaiMV1ExAUwF49mVQKaLnJZ+n2b51dXVsIyCJRols48zYxl+366LrgJCHnt5fnldljawhLU5fJzuwrn4fZbEoC5slMwZuQXoIC8ZqQtkTJYYtI7kdaYOWBaN4j1ZWIlFovwMadChs5tGk3H6aNnsht85vua1fG8/0/WiMl5FxQQwFYznNl00F0fl2KLmgWWzQzpdaaVimXinnnCvzoRO6mF+3xLe+kq/328+s5Pd7GAp6dZZUfmFKISMeiWLBDFooAwsYPhaNC+6cthazM/C1mHDc6W53TuHcm5MUKUDneF8ka5Hay7n5J0H3QjcIZSfeZxRwH5U3Gm9qIxXUTEBTAXjSWvSLEq9iJpOGGVyJ61tdDpTUlkKRvoRWcag85oB3CVz2tlunYx+ODY3oV8vAgs70cnN8oDlMUyejRJ8yXgsVGs92/46j93WT+s+fGb9fr85hgxH1qX1NiqNwbF6bN6FRClMUjvJmrod/ZssCVKtmhUVJwGmgvFcztvSw9LTkipqUsISB6VllEw3LrA4KtDKe0RNPfzakvPQoUONXmPfIa2ZtgCyNDmLOzGqhilMbLDSlSjMZFC2+IoCrX2ex0pGZJNIszd3BGTvfr/f0vGiIkZRoHtUFIvPmgWMaJksA/TZYKbLV1zey3jMSz9UVFQ89pgKxrNFk4VI6c+jlZOpP6urq2GpdqbWMOLAsZa7du2S1JaalHxRQ5FSN2SbLY+JhVn9vlmFEp5l5tlkk7oddZ7l5eWW9Oc8DLIJ/XS8J9uaMT4yKn84NzcXtvyiRZo6HXcxfN8WZkasMDqljF7iPelb5S4kKha1XlTGq6iYAKaC8XLOOnr0aKvZBNFVnKZ8f3V1tZFqUZaBYQllX9g3vvENSdIll1zSXKtEFJVuvY33KzMKzHCUioxc8bXIgPaNMa6UvkgyQWkVZBKoj/G9qMtSB6TvlNZcWn/pe+P1t27d2kqHiooU0T9pcL4+7sc//vHI+4YZzutcxrxG0UOcN9d2o6iMV1ExAUwF41nHY7l1FjOKdJmyaM2+fftGrvFo95SGPrbrrrtOkvT6179+5J6MJqGlMYqCmJuba/mEokiTKDuBuq6Ptw7DoquMYTTKtY2supGOE0XRsCgQiyRF7dLM/vPz8y39u4upy3uTRbkLsfXYzUUZpUM2Li2Z0a7Bz4CNYDbansuojFdRMQGkcbrQCRlESvskHZL04KTHEuB01bFtBNM6tsdzXBfknM8Yd9BU/PAkKaV0Q8750kmPowt1bBvDtI5tGsZVt5oVFRNA/eFVVEwA0/TDu2rSA3gU1LFtDNM6tomPa2p0vIqKJxOmifEqKp40mIofXkrpl1JK3xk0tLxyguM4L6V0XUppd0rp1pTS2wfvdzbinNAYeymlb6SUrh283pVSun4wtn9MKbXrGZyYcW1PKX0kpXTbYP1eNC3rllL6vcHz/FZK6QMppflJr9vEf3gppZ6k/6m1ppbPkPQbKaVnTGg4y5J+P+f8dEkvlPQ7g7G4EefFkr4weD0pvF3S7uL1X0h692BsP5b01omMSnqPpM/knJ8m6dlaG+PE1y2ldK6kt0m6NOf8TEk9SW/UpNct5zzRf5JeJOmzxet3SnrnpMc1GMs1kl4j6TuSzhm8d46k70xoPDu19gV+paRrtdYW7UFJ/a61PIHj2ibpTg1sBsX7E183DXs27tBaiOS1kv7NpNdt4oynuJnlRJFSulDScyVdLzTilHRmfObjir+W9AeSnDrxFEn7c85O15jU2l0kaZ+kvx9sg/82pbSgKVi3nPO9kv5Sa8117pf0sKQbNeF1m4Yf3rqbWZ4opJROkfRRSb+bc/7JuONPBFJKvyzpgZzzjeXbHYdOYu36kp4n6b055+dqLfxvktvxBgO98nJJuyQ9VdKC1tQa4oSu2zT88DbczPLxQEppVms/uvfnnD82eHtvGvRyT6ONOE8kXiLpspTSDyR9UGvbzb/WWp95h8pPau32SNqTc75+8PojWvshTsO6vVrSnTnnfTnno5I+JunFmvC6TcMP7+uSLh5YmTZpTfH9xCQGktZyPt4naXfO+b8VH0WNOE8Ycs7vzDnvzDlfqLU1+mLO+U2SrpP0hgmP7YeS7kkpXTJ461WSvq0pWDetbTFfmFLaMni+Httk1+1EK7uBAvxaSd+V9D1J/2WC43ip1rYct0i6efDvtVrTpb4g6fbB3x0TXq9XSLp28P+LJP0/rTUD/bCkuQmN6TmSbhis3cclnTYt6ybpv0q6TdK3JP1vSXOTXrcauVJRMQFMw1azouJJh/rDq6iYAOoPr6JiAqg/vIqKCaD+8CoqJoD6w6uomADqD6+iYgKoP7yKigng/wNVBOwxOVqrsgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# To see the first face:\n",
    "plt.imshow(M[0].reshape((112, 92)), cmap='gray'); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il y a 400 images et ces images ont 92*112 pixels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Propriété 1:On remarque tout d'abord que si une fonction f $\\mathbb{R}^{n}\\rightarrow\\mathbb{R}$ est convexe alors toute fonction $\\tilde{f}$ définit par $\\forall(x,y)\\in \\mathbb{R}^{2}, \\tilde{f}(x,y)=f(0,...,0,x,0,...,0,y,0,..)$ est aussi convexe.\n",
    "Considérons la fonction g:$(x,y)\\rightarrow (1-xy)^{2}$\n",
    "\n",
    "On a: $g(0,0)=1$, $g(1,1)=0$ or $g(1/2,1/2)=9/16$\n",
    "\n",
    "Ainsi, $1/2\\times g(0,0)+1/2\\times g(1,1)\\leq g(1/2,1/2)$\n",
    "\n",
    "Donc g n'est pas convexe.\n",
    "\n",
    "Or si g n'est pas convexe alors par contraposé de la propriété 1, la focntion objectif n'est elle aussi pas convexe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Claculons le gradient de la fonction par la methode de la différentielle.\n",
    "Commençons par calculer la première partie du gradient en différenciant par rapprot à H:\n",
    "On a: $\\frac{1}{2np}||M-W(H+\\delta H)||^{2}_{2}=\\frac{1}{2np}<M-W(H+\\delta H);M-W(H+\\delta H)>$\n",
    "                                         $=\\frac{1}{2np}(||M-WH||^{2}_{2}-2<M-WH;W\\delta H> +o(H))$\n",
    "                                         $=\\frac{1}{2np}(||M-WH||^{2}_{2}-2<^{t}W(M-WH);\\delta H>) +o(H)$\n",
    "                                         \n",
    "                                         \n",
    "                                         \n",
    "On en déduit par identification que la première partie du gradient est $\\frac{-1}{np}^{t}W(M-WH)$\n",
    "\n",
    "\n",
    "Maintenant calculons la deuxième parite du gradient par rapport à W:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{1}{2np}||M-(W+\\delta W)H||^{2}_{2}&=\\frac{1}{2np}<M-(W+\\delta W)H;M-(W+\\delta W)H> \\\\\n",
    "                                        &=\\frac{1}{2np}(||M-WH||^{2}_{2}-2<M-WH;\\delta WH> +o(W))\\\\\n",
    "                                        &=\\frac{1}{2np}(||M-WH||^{2}_{2}-2<(M-WH)^{t}H;\\delta W> +o(W))\n",
    "                                        \\end{align}\n",
    "\n",
    "\n",
    "Donc la deuxième partie du gradient est $\\frac{-1}{np}(M-WH)^{t}H$\n",
    "\n",
    "Si on appelle g la fonction objectif alors son gradient est donné par:\n",
    "$\\nabla g=\\begin{pmatrix} \n",
    "\\frac{-1}{np}(M-WH)^{t}H  \\\\\n",
    "\\frac{-1}{np}^{t}W(M-WH)\n",
    "\\end{pmatrix}$\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "Montrons maintenant que la gradient n'est  pas Lipschitz continue,\n",
    "et pour cela nous allons nous placer dans la cas où n=k=p=1:\n",
    "\n",
    "On rappelle que le gradient de G est lipschitzien s'il existe L $\\in \\mathbb{R}^{+}$ tel que:\n",
    "$|| \\nabla G (W_{1},H_{1})- \\nabla G(W_{2},H_{2})|| \\leq L||(W_{1},H_{1})-(W_{2},H_{2})||$\n",
    "\n",
    "\n",
    "prenons maintenant $W_{1} et H_{1}$ les fonctions identités et $W_{2}, H_{2} $nulles.\n",
    "\n",
    "\n",
    "Dans notre cas:\n",
    "$|| \\nabla G (W_{1}(x),H_{1}(x))- \\nabla G(W_{2}(x),H_{2}(x)||^{2}=(-x(m-x^{2})+m^{2})^{2}+(-x(m-x^{2})+m^{2})^{2}=2x^{2}(x^{2}-m)^{2}$ \n",
    "\n",
    "Or:\n",
    "\n",
    "$||(W_1(x)-W_2(x),H_1(x)-H_2(x))||^2 = (x-0)^2 + (x-0)^2 = 2x^2$\n",
    "\n",
    "Ainsi si on écrit l'inégalité de lichitz on obtient que $\\forall x \\in \\mathbb{R}$\n",
    "$2x^{2}(x^{2}-m)^{2}\\leq L2x^{2}$\n",
    "en divisant par $x^{2}$, on obtient qu'un polynome de degrès 4 est inférieur à une constante  $\\forall x \\in \\mathbb{R}$ Ce qui est ilmpossible donc le gradient n'est pas lipschitzien."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'avantage de ce choix est que l'on est pas loin du resultat à obtenir à savoir la factorisation par deux matrices W et H.\n",
    "En effet la décomposition svd permet de transformer M en un produit de trois matrices U, $\\Sigma$ et V $M=U\\Sigma V$ avec $\\Sigma$ une matrice qui ressemble à une homothètie mais pas totalement car elle est de dimentinon $p\\times k$ avec $p\\ne k$. On se rend bien compte que l'on est proche de la factorisation WH c'est donc un bonne initialisation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour ce qui est du gradient nous l'avons déjà calculé précédemment car il s'agit de la deuxième composante du gradient de g précédemment définit soit:\n",
    "On sait que la fonction carré préserve la convexité donc montrer que g est convexe revient à montrer que:\n",
    "\n",
    "$h:W \\rightarrow ||M-WH_{0}||$\n",
    "est convexe:\n",
    "\n",
    "\n",
    "soit \n",
    "$(W_{1},W_{2})\\in M_{n,k}$ et  $\\lambda \\in [0,1]$\n",
    "\n",
    "\\begin{align}\n",
    "||M-(\\lambda W_{1}+(1-\\lambda)W_{2})H_{0}|| &=||\\lambda M+(1-\\lambda) \\\\\n",
    "    &=M-(\\lambda W_{1}+(1-\\lambda)W_{2})H_{0}|| \\\\\n",
    "    &=||\\lambda M-\\lambda W_{1}H_{0}+(1-\\lambda)M-(1-\\lambda)W_{2}H_0|| \\\\\n",
    "    &=\\lambda|| M- W_{1}H_{0}||+(1-\\lambda)||M-W_{2}H_0||\n",
    "\\end{align}\n",
    "\n",
    "Par inégalité triangulaire on obtient bien que g est convexe\n",
    "    \n",
    "\n",
    "\n",
    "$ \\nabla g =\\frac{-1}{np}^{t}W(M-WH_{0})$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_g_w(W):\n",
    "    n,p=M.shape\n",
    "    L=M-np.dot(W,H0)\n",
    "    return (1/(2*n*p))*np.linalg.norm(L,\"fro\")**(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_grad_g(W):\n",
    "    n,p=M.shape\n",
    "    return -1/(n*p)*np.dot(M-np.dot(W,H0),H0.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$prox_{\\gamma\\iota_{\\mathbb{R}^{+}}}=argmin({\\gamma\\iota_{\\mathbb{R}^{+}}(y)+\\frac{1}{2}||y-x||^{2}})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, on remarque que si $x<0$ alors y ne peut pas être négatif car sinon le prox vaudrait $+\\infty$ pour minimiser il faut donc que y soit nul.\n",
    "\n",
    "Dans le cas où $0\\leq x$ alors en prenant $y=x$ on remarque ${\\gamma\\iota_{\\mathbb{R}^{+}}(y)+\\frac{1}{2}||y-x||^{2}}=0$ et $y=x$ est le seul point pour lequel $||x-y||=0$ \n",
    "\n",
    "Donc:\n",
    "Pour $x<0$, le prox vaut 0 et pour $ 0 \\leq x$ le prox vaut x:\n",
    "\n",
    "Donc le prox correspond donc bien à la projection sur $\\mathbb{R}^{+}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def projected_gradient_method(val_g, grad_g,W0, gamma,N):\n",
    "    W=W0\n",
    "    for i in range (N):\n",
    "        W=W-gamma*grad_g(W)\n",
    "        W=np.maximum(W,0)\n",
    "    return val_g(W)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.222664966607934\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "442.2878543960478"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W0, S, H0 = scipy.sparse.linalg.svds(M,k)\n",
    "W0 = np.maximum(0, W0 * np.sqrt(S)) \n",
    "H0 = np.maximum(0,(H0.T * np.sqrt(S)).T)\n",
    "n,p=M.shape\n",
    "K=np.dot(np.transpose(H0),H0)\n",
    "L0=np.linalg.norm(K,\"fro\")\n",
    "print((n*p)/L0)\n",
    "projected_gradient_method(compute_g_w,compute_grad_g,W0,(n*p)/L0,100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 4.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons ici utiliser le line search:\n",
    "Cette methode consiste à trouver à chaque étape une valeur de $\\gamma$ idéale permettant d'éviter certains problèmes.\n",
    "Comme le gradient est lipschitz continue, comme il est indiqué à la question 3.2 cela est possible.\n",
    "à chaque étape on doit trouver $\\gamma_{k}$ tel que:\n",
    "$\\gamma_{k}=\\underset{\\gamma \\in \\underset{+}{\\mathbb{R}}}{argmin} \\ G(x_{k}-\\gamma \\nabla f(x_{k}))$\n",
    "Pour trouver ce $\\gamma_k$ on va dériver cette fonction pour en trouver le minimum.\n",
    "Si on nomme $F_k$ cette fonction alors on a:\n",
    "$\\frac{d F_k}{d \\gamma}(\\gamma)=<-\\nabla g(x_k), \\nabla g(x_k - \\gamma \\nabla g(x_k))>$\n",
    "\n",
    "\n",
    "Le gamma que l'on cherche est donc celui qui annule cette expression.\n",
    "\n",
    "Pour le trouver nous allons utiliser la mathode de newton."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def line_search_projected_gradient(val_g, grad_g, W0, gamma0,N):\n",
    "    W=W0\n",
    "    gamma = gamma0\n",
    "    for i in range (1,N):\n",
    "        F= lambda gamma: np.trace(-np.dot(grad_g(W).T,grad_g(W - gamma*grad_g(W))))\n",
    "        W = ((W - gamma*grad_g(W))>0)*(W - gamma*grad_g(W))\n",
    "        gamma = scipy.optimize.newton(F, gamma)\n",
    "    return val_g(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "442.6972009318645"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "line_search_projected_gradient(compute_g_w, compute_grad_g, W0, n*p/L0,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 4.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On remarque que le second algorithme met moins de temps à s'executer que le premier. C'est tout à fait normal car la procédure du line search est plus efficace."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 5.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_g_w_2(W,H):\n",
    "    n,p=M.shape\n",
    "    L=M-np.dot(W,H0)\n",
    "    W=W.reshape((n,k))\n",
    "    return (1/(2*n*p))*np.linalg.norm(L,\"fro\")**(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_grad_g_2(W,H):\n",
    "    n,p=M.shape\n",
    "    W=W.reshape((n,k))\n",
    "    return vectorize((-1/n*p)*np.dot((M-np.dot(W,H)),H.T),-1/(n*p)*np.dot(W.T,M-np.dot(W,H)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def projected_gradient_method_2(valg,gradg,W0,H0,gamma0,N):\n",
    "    W=W0\n",
    "    H=H0\n",
    "    gamma=gamma0\n",
    "    for i in range (N):\n",
    "        W,H=unvectorize_M(np.maximum(0,vectorize(W,H)-gamma*gradg(W,H)),M)\n",
    "        L=unvectorize_M(vectorize(W,H)-gamma*gradg(W,H),M)\n",
    "        F=lambda gamma: -np.dot(gradg(W,H).T,gradg(L[0],L[1]))\n",
    "        gamma=scipy.optimize.newton(F,gamma)\n",
    "\n",
    "    return vectorize(W,H), valg(W,H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00019999999999997797 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.0002000099999999172 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020002000049990087 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.000200030001499929 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020004000300000158 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.0002000500050001186 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020006000750050212 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.0002000700105009301 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020008001400140252 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020009001800214143 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.0002001000225029248 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.0002001100275041967 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020012003300551306 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020013003900709592 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.0002001400455091673 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020015005251128315 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020016006001388753 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020017006801698045 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020018007652033987 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020019008552418782 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020020009502852432 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.0002002101050331273 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020022011553844088 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.000200230126544243 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020024013805053364 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020025015005753488 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020026016256502466 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py:195: RuntimeWarning: Tolerance of 0.00020027017557322502 reached\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\Users\\Antoine\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:4: RuntimeWarning: overflow encountered in multiply\n",
      "  after removing the cwd from sys.path.\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "Failed to converge after 50 iterations, value is nan",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-14-df8d6a06707d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mprojected_gradient_method_2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcompute_g_w_2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcompute_grad_g_2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mW0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mH0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-13-f11dcd4c3992>\u001b[0m in \u001b[0;36mprojected_gradient_method_2\u001b[1;34m(valg, gradg, W0, H0, gamma0, N)\u001b[0m\n\u001b[0;32m      7\u001b[0m         \u001b[0mL\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0munvectorize_M\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvectorize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mW\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mH\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mgamma\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mgradg\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mW\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mH\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mM\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m         \u001b[0mF\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mlambda\u001b[0m \u001b[0mgamma\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;33m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgradg\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mW\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mH\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mgradg\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mL\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mL\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m         \u001b[0mgamma\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mscipy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0moptimize\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnewton\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mF\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mgamma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     10\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     11\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mvectorize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mW\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mH\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvalg\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mW\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mH\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\zeros.py\u001b[0m in \u001b[0;36mnewton\u001b[1;34m(func, x0, fprime, args, tol, maxiter, fprime2)\u001b[0m\n\u001b[0;32m    204\u001b[0m             \u001b[0mq1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    205\u001b[0m     \u001b[0mmsg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m\"Failed to converge after %d iterations, value is %s\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mmaxiter\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 206\u001b[1;33m     \u001b[1;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    207\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    208\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mRuntimeError\u001b[0m: Failed to converge after 50 iterations, value is nan"
     ]
    }
   ],
   "source": [
    "projected_gradient_method_2(compute_g_w_2,compute_grad_g_2,W0,H0,1,1000)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 5.2 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons la tième itération on a par définition:\n",
    "\n",
    "$W_t = \\underset{W}{argmin} \\frac{1}{2np} ||M - WH_{t-1}||^2$\n",
    "\n",
    "$H_t = \\underset{H}{argmin} \\frac{1}{2np} ||M - W_t H||^2$\n",
    "\n",
    "On a alors les inégalités suivantes:\n",
    "\n",
    "$\\frac{1}{2np} ||M - W_t H_t||^2 \\le \\frac{1}{2np} ||M - W_t H_{t-1}||^2$\n",
    "\n",
    "$\\frac{1}{2np} ||M - W_t H_{t-1}||^2 \\le \\frac{1}{2np} ||M - W_{t-1} H_{t-1}||^2$\n",
    "\n",
    "Par définition\n",
    "\n",
    "On en déduit donc que la suite est décroissante.\n",
    "Puis comme elle est minorée par 0 alors elle est convergente"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 5.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour coder cette methode, il  suffit de déterminer $W_{t}$ et $H_{t}$ à chaque étape. Pour déterminer les minimums des fonctions nous allons encore utiliser la methode de la descente du gradient. Sauf que cette fois si il s'agit de déterminer les gradients en fixant une des deux variables.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_H(W,H):\n",
    "    return -(1/(n*p))*np.dot(W.T,M - np.dot(W,H))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_W(W,H):\n",
    "    return -(1/(n*p))*np.dot((M-np.dot(W,H)),H.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def determ_W_t(g, gradw, W0, H0, gamma, N):\n",
    "    W,H = W0, H0\n",
    "    for i in range(1,N):\n",
    "        W = np.maximum(0,W - gamma*gradw(W,H))\n",
    "    return W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def determ_H_t(g, gradh, W0, H0, gamma, N):\n",
    "    W,H = W0, H0\n",
    "    for i in range(1,N):\n",
    "        H = np.maximum(0,H - gamma*gradh(W,H))\n",
    "    return H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def alternate_minimization_method(g, gradh, gradw, W0, H0, gamma, N):\n",
    "    W, H = W0, H0\n",
    "    for i in range(1,N):\n",
    "        H = determ_W_t(g, gradh, W, H, gamma, N)\n",
    "        W = determ_H_t(g, gradw, W, H, gamma, N)\n",
    "    return g(vectorize(W,H))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def g_modified(W):\n",
    "    W, H = unvectorize_M(W,M)\n",
    "    return (1/(2*n*p))*np.linalg.norm(M-np.dot(W,H),ord='fro')**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alternate_minimization_method(g_modified,gradient_H, gradient_W, W0, H0, n*p/L0,1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 5.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On remarque que la methode alternate  est très longue en temps d'execution, beaucoup plus longue que la methode précédemment utilisée"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 5.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il s'agit ici d'un problème de convergence et de limite. Introduire un $\\epsilon$ suffisamment petit de manière à arreter le processus lorque $||M-WH||<\\epsilon$ est la solution la intuitive. Elle permet aussi de limiter des itérations parfois inutile lorsque la précision est déjà bonne."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}