EM_nonsparse.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3861c248-aff2-41ba-aeb1-f2bdcb847a9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "36bc47ac-906f-4173-8d55-cd7faebe6f52",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.load(\"simulated_large_A.npy\")\n",
    "x_flat = np.load(\"simulated_large_x.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1a0ae052-73c3-4f47-b3fe-ff5dc25c0a08",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3108, 2500), (2500,))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.shape, x_flat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d237a0ec-cb9d-4f3c-8f6d-da4ee91182be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image vector x is of shape (2500,)\n",
      "Mixing matrix A is of shape (3108, 2500)\n",
      "Observation matrix Y is of shape (3108,)\n",
      "M=2500, N=3108\n"
     ]
    }
   ],
   "source": [
    "Ax = A@x_flat\n",
    "y = np.random.poisson(Ax).reshape(len(A[:, 0]))\n",
    "N, M = A.shape\n",
    "\n",
    "print(\"Image vector x is of shape {}\".format(np.shape(x_flat)))\n",
    "print(\"Mixing matrix A is of shape {}\".format(np.shape(A)))\n",
    "print(\"Observation matrix Y is of shape {}\".format(np.shape(y)))\n",
    "print(f\"M={M}, N={N}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "33ac0e70-928f-4cc0-964a-2592cac0e784",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "step: 0, diff: 3852.9287154299577, mse: 1481.6489379363384\n",
      "step: 1, diff: 239.07995619466607, mse: 1362.6899482719655\n",
      "step: 2, diff: 146.54551973989464, mse: 1283.1207391362934\n",
      "step: 3, diff: 106.49011642760752, mse: 1219.8361733350876\n",
      "step: 4, diff: 85.7723614370583, mse: 1165.8964610276946\n",
      "step: 5, diff: 73.20595434376868, mse: 1118.5928358227402\n",
      "step: 6, diff: 64.50284708393981, mse: 1076.5811955527197\n",
      "step: 7, diff: 57.873579508267284, mse: 1039.0402543404894\n",
      "step: 8, diff: 52.501109420180214, mse: 1005.3851945898782\n",
      "step: 9, diff: 47.97305761367315, mse: 975.1570654541713\n",
      "step: 10, diff: 44.05922017874083, mse: 947.9744765699625\n",
      "step: 11, diff: 40.618859039065896, mse: 923.5098121450085\n",
      "step: 12, diff: 37.55915280163587, mse: 901.4761459262205\n",
      "step: 13, diff: 34.815036761619844, mse: 881.619308984794\n",
      "step: 14, diff: 32.33862116801013, mse: 863.7126599466859\n",
      "step: 15, diff: 30.09323142801531, mse: 847.5533714666008\n",
      "step: 16, diff: 28.04985008797984, mse: 832.9596128292952\n",
      "step: 17, diff: 26.184892208770727, mse: 819.7682863822465\n",
      "step: 18, diff: 24.478761752040285, mse: 807.8331227522772\n",
      "step: 19, diff: 22.91488444658536, mse: 797.0230220591967\n",
      "step: 20, diff: 21.47904041726303, mse: 787.2205755598516\n",
      "step: 21, diff: 20.158890122790357, mse: 778.3207294228454\n",
      "step: 22, diff: 18.943627781055405, mse: 770.2295679152705\n",
      "step: 23, diff: 17.82372081106832, mse: 762.8632019760743\n",
      "step: 24, diff: 16.790708754314082, mse: 756.146753824026\n",
      "step: 25, diff: 15.83704443041832, mse: 750.0134306297838\n",
      "step: 26, diff: 14.955965921760049, mse: 744.4036814056806\n",
      "step: 27, diff: 14.141391684196972, mse: 739.2644317420861\n",
      "step: 28, diff: 13.387833452861193, mse: 734.5483911938508\n",
      "step: 29, diff: 12.690323152204659, mse: 730.2134281857861\n",
      "step: 30, diff: 12.044351037913167, mse: 726.222007360894\n",
      "step: 31, diff: 11.445812987167542, mse: 722.5406843840849\n",
      "step: 32, diff: 10.890965332104098, mse: 719.1396533526192\n",
      "step: 33, diff: 10.376385973072585, mse: 715.9923421520091\n",
      "step: 34, diff: 9.898940759382665, mse: 713.0750513247743\n",
      "step: 35, diff: 9.455754314564313, mse: 710.3666322786273\n",
      "step: 36, diff: 9.04418462918564, mse: 707.8482009395149\n",
      "step: 37, diff: 8.661800859021122, mse: 705.5028832436216\n",
      "step: 38, diff: 8.306363857961887, mse: 703.3155891526726\n",
      "step: 39, diff: 7.97580904913696, mse: 701.2728121620682\n",
      "step: 40, diff: 7.668231298332341, mse: 699.3624515465875\n",
      "step: 41, diff: 7.381871503939859, mse: 697.5736548500541\n",
      "step: 42, diff: 7.11510465960066, mse: 695.8966783711556\n",
      "step: 43, diff: 6.866429181135342, mse: 694.3227636261688\n",
      "step: 44, diff: 6.634457319587875, mse: 692.8440279801175\n",
      "step: 45, diff: 6.417906508235136, mse: 691.4533678308551\n",
      "step: 46, diff: 6.215591513976527, mse: 690.1443729061746\n",
      "step: 47, diff: 6.026417283155676, mse: 688.9112503930395\n",
      "step: 48, diff: 5.849372389001341, mse: 687.7487577613447\n",
      "step: 49, diff: 5.6835230027889025, mse: 686.6521432733075\n",
      "step: 50, diff: 5.528007323755476, mse: 685.61709328479\n",
      "step: 51, diff: 5.382030413931512, mse: 684.6396855476804\n",
      "step: 52, diff: 5.244859393532418, mse: 683.7163478140418\n",
      "step: 53, diff: 5.115818960532304, mse: 682.8438211241187\n",
      "step: 54, diff: 4.994287204659604, mse: 682.0191272325108\n",
      "step: 55, diff: 4.879691691450918, mse: 681.2395396907942\n",
      "step: 56, diff: 4.771505796321856, mse: 680.5025581614849\n",
      "step: 57, diff: 4.669245272001155, mse: 679.8058855882806\n",
      "step: 58, diff: 4.572465035272433, mse: 679.1474078917212\n",
      "step: 59, diff: 4.4807561609071405, mse: 678.5251758984137\n",
      "step: 60, diff: 4.393743072081598, mse: 677.9373892463882\n",
      "step: 61, diff: 4.311080917562068, mse: 677.3823820394879\n",
      "step: 62, diff: 4.232453126615273, mse: 676.8586100504484\n",
      "step: 63, diff: 4.157569133046478, mse: 676.3646392958901\n",
      "step: 64, diff: 4.086162260049616, mse: 675.8991358272091\n",
      "step: 65, diff: 4.017987757737092, mse: 675.4608565996559\n",
      "step: 66, diff: 3.9528209853433256, mse: 675.048641298002\n",
      "step: 67, diff: 3.890455730198347, mse: 674.6614050113969\n",
      "step: 68, diff: 3.8307026556743393, mse: 674.2981316625272\n",
      "step: 69, diff: 3.7733878704283073, mse: 673.9578681072161\n",
      "step: 70, diff: 3.7183516114146444, mse: 673.6397188303092\n",
      "step: 71, diff: 3.665447033320601, mse: 673.3428411722665\n",
      "step: 72, diff: 3.6145390972913254, mse: 673.0664410284156\n",
      "step: 73, diff: 3.5655035520519296, mse: 672.8097689694947\n",
      "step: 74, diff: 3.5182260008036166, mse: 672.5721167379698\n",
      "step: 75, diff: 3.472601047562324, mse: 672.3528140798088\n",
      "step: 76, diff: 3.4285315169132695, mse: 672.1512258759502\n",
      "step: 77, diff: 3.3859277414755935, mse: 671.9667495417606\n",
      "step: 78, diff: 3.3447069116937884, mse: 671.7988126663255\n",
      "step: 79, diff: 3.3047924829009716, mse: 671.6468708665786\n",
      "step: 80, diff: 3.2661136349234314, mse: 671.5104058340567\n",
      "step: 81, diff: 3.228604779811904, mse: 671.3889235545274\n",
      "step: 82, diff: 3.1922051136005543, mse: 671.2819526829221\n",
      "step: 83, diff: 3.1568582082878853, mse: 671.1890430579258\n",
      "step: 84, diff: 3.1225116405265685, mse: 671.109764342289\n",
      "step: 85, diff: 3.089116653778183, mse: 671.043704776438\n",
      "step: 86, diff: 3.0566278509501026, mse: 670.9904700342938\n",
      "step: 87, diff: 3.0250029147745314, mse: 670.9496821714083\n",
      "step: 88, diff: 2.9942023534169526, mse: 670.9209786565705\n",
      "step: 89, diff: 2.9641892690170284, mse: 670.9040114789766\n",
      "step: 90, diff: 2.9349291470603447, mse: 670.898446323886\n",
      "step: 91, diff: 2.9063896646652885, mse: 670.9039618104267\n",
      "step: 92, diff: 2.8785405160380413, mse: 670.9202487858635\n",
      "step: 93, diff: 2.851353253505858, mse: 670.9470096712334\n",
      "step: 94, diff: 2.8248011426826105, mse: 670.9839578537675\n",
      "step: 95, diff: 2.798859030452407, mse: 671.0308171219825\n",
      "step: 96, diff: 2.7735032245788163, mse: 671.0873211397388\n",
      "step: 97, diff: 2.748711383856853, mse: 671.1532129559271\n",
      "step: 98, diff: 2.7244624178264747, mse: 671.2282445467778\n",
      "step: 99, diff: 2.7007363951587458, mse: 671.3121763880771\n"
     ]
    }
   ],
   "source": [
    "x_et, diff, mse, step = mle_em(100, A, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "13367f96-487b-477e-b873-ae04cbea400f",
   "metadata": {},
   "outputs": [],
   "source": [
    "size = int(np.sqrt(len(x_flat)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "bdd5ba94-6476-4e76-b6db-4a298a8adced",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAC2CAYAAAB6fF5CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2pklEQVR4nO2dd9RdVbnunxl6FRBFIIGAEGpI6C1AUIFQTyCgCEjHc49IUblw9XCAo+IQhl7wAnKVS0BAmgaCFGHQQmgJvQRI6EiXKiBBUdb9Y69vfr+1XM/ORpP15XjeZwwG77f32rPPlfm88y2pKAoFAoFAoB0MGugGBAKBwH8nxEs3EAgEWkS8dAOBQKBFxEs3EAgEWkS8dAOBQKBFxEs3EAgEWkS8dOdipJQmpZQOHuh2BOYepJS+k1L6f7P72R7KKlJKq8yOsv67Y7a/dFNKz6aUZqaU3kspvZJSOjeltOjsrmd2YE4upJTS0LL8eedQ+SeklC6YE2UH2kFKaf+U0sMppffLvXJmSmmJbr8piuIHRVH09A/xx3n2H0EcDj4e5tRJd+eiKBaVNFLSupK+PYfqmaOYUy/MQCCl9C1JJ0n6n5I+IWkTSStKuj6lNL/5TazHfwLMUfVCURSvSLpOnZevJCmltElK6Y6U0tsppQdTSqPx3VIppXNSSi+llN5KKU3Ed4eklJ5MKb2ZUvpNSmk5fFeklP5HSumJstwzUkqp/G6VlNItKaU/pJReTyldUn4+ufz5g+Wp/EsppdEppRdSSseklF6RdE55GrmN/eIJOaW0UErpxyml58o6bkspLSSpr/y3y/I3LZ8/MKX0WNm/61JKK6LcbVJK08tyTpeUeh3rsk1fK8fg3ZTS91JKny3H+p2U0qV9mzmltGRK6aqU0mtlO65KKQ1GWSullCaX5dxQjucF+N7OYWDWSCktLuk/JR1WFMW1RVF8WBTFs5K+KGmopH3K505IKf06pXRBSukdSfvXGU5Kad9y7b2RUvqP1GGaX8DvLyjlPua1X0rpd+Ve+HeUs1FK6c5yTl9OKZ3uXv6z6FvfHjo6pfT7sqyxKaUdUkqPl/v3O73Wm1LaNqU0o9wTPy338sH43u6nuRZFUczW/yQ9K+kLpTxY0sOSflL+vbykNyTtoM4Lf5vy70+V318t6RJJS0qaT9JW5eefk/S6pPUkLSDpNEmTUWch6SpJS0haQdJrksaU310k6d/L+haUNKr2u1Xw92hJf1HnBLKApIUk7S/ptlof8+8knSFpUtm3eSRtVv52aPncvPjdv0h6UtIakuaVdKykO8rvlpb0rqTdy75/o2zLwWacT5B0Qa1NV0haXNJakv4k6UZJK6tzknpU0n7ls5+UNE7SwpIWk/QrSRNR1p2SfiRpfkmjJL3TV9es5jD+62mPjCnndt6G734h6SLM8YeSxpZjvRDnXdKakt4r52j+cs4+VP/+47N96/GsspwR5RpZo/x+fXVO2/OWzz4m6Ui3V2ptntS3TtW/h44r1/Eh6uzHC8u1tpakmZJWmlW95Z54R9Ju5fdHlP3rq8vup7n5vzmxoJ4tF8K75UTdKGmJ8rtjJJ1fe/46SftJWlbSR5KWbCjzbEkn4+9Fy8EfigXBl+mlkv5XKZ8n6eeSBjeU2/TS/bOkBfHZ/jIv3XIjzJQ0oqHsvkXOl+5vJR2EvwdJel8dWrmvpCn4Lkl6QR/vpbs5/r5X0jH4+8eSTjVljZT0VimvUG6ahfH9BerfvHYOB3ox/1f5T52T7Cvmux9Kuh5zPLn2/QmYi+NUvqDLvxcu12+3l+5gPH+XpD1NO46UdHltffX60p0paZ7y78XK325cW5tjZ1VvuSfuxHdJ0vOoy+6ngZ7jbv/NKfXC2KIoFisnYHV1/sWSOi+XPUoq8XZK6W11/pVeVtIQSW8WRfFWQ3nLSXqu74+iKN5T53S1PJ55BfL76ryYJelodSbrrpTSIymlA2fR9teKovhg1l2U1OnXgpKe6vH5FSX9BH1/s2zb8ur08fm+B4vOKnq+qZAueBXyzIa/F5WklNLCKaWflbT0HXVUIUuklOYp2/FmURTv47dsR7c5DPSG1yUtnZp1tMuW3/eh2xqor5n31dkX3dC4T1JKw0o10yvlmviB+vftx8UbRVH8tZRnlv93a7FbvU174gWU020/zbWY0zrdWySdqw7tkToDeH5RFEvgv0WKovhh+d1Sqfn29iV1BliSlFJaRB2K/GIPbXilKIpDiqJYTtK/Svpp6m6xUA+79kd1ThB9dX8G370u6QNJn+2hHKnTx3+t9X+hoijukPSyOv/w9NWT+PdsxrckrabO6WNxSVv2VVu2Y6mU0sJ4nu3oNoeB3nCnOtR+N36YOlY+26vDDvvQLQzgy+qo8Pp+v5A6++LvwZmSpktatVwT39HHuFP4B9Ct3nr/Ev9W9/0016INO91TJW2TUhqhDk3dOaW0XUppnpTSgqXifXBRFC+rQxd+Wl70zJdS6nsZXCTpgJTSyJTSAur8azi16Fw+dEVKaQ9cEr2lziL+qPz7VXV0nt3woKS1yroXVIeySZKKovhI0nhJ/zultFzZp03LNr5W1sPy/6+kb6eU1irb9omU0h7ld1eX9exWnoAOl8QX/OzEYuqcNt5OKS0l6Xj06TlJ90g6IaU0f+pcAO6M39o5nENt/adDURR/UOci7bSU0phyrQ9VRy32gqTzeyzq1+rMxWbl5dMJ+vtflIupoz99L6W0uqR/+zvLmZ31Xi1peHkRN6+kQ1XdE93201yLOf7SLYriNXX0qscVRfG8Osrv76jzUnpeHZOZvnZ8RR1d7XRJv1dHv6OiKG6Q9B+SJqjzr99nJe3ZYxM2lDQ1pfSepN9IOqIoiqfL706Q9IuSnnzRtP9xSd+VdIOkJyTdVnvkKHUuC+9Wh96cJGlQSfVOlHR7Wf4mRVFcXn5/cUmlpqlzslFRFK9L2kMdnd4bklaVdHuPffy4OFWdy5TXJU2RdG3t+70lbVq24/vqXG7+qWznrOYw0AOKojhZnTH8kTovnanqjOXni6L4U49lPCLpMEkXq7Mv3lNn3/T0+xqOkrSXOncxZ6kz523A1os9cbI6a3FNdQ4EfWvR7qe5GalUQAcCFqljZje9KIrjZ/lwYMBQqifeVoeqPzPAzZntSCkNUocJ7F0Uxc0D3Z6/F3E6CfwNUkobpo6N76CU0hh1TrYTB7hZgQaklHYuL0YXUefU/LA6FkT/FCjVWEuUKrs+fe+UAW7WP4R46Qaa8Bl1zIDek/R/JP1bURT3D2iLAg7/os5F80vqqKT2LP656Oum6lgHva7O3cLYoihmdv/J3I1QLwQCgUCLiJNuIBAItIh46QYCgUCL6Bq1aJ555sm6B6ohOjbKf/t5L6qKhRZaKMvHH1+9DH/ppZey/KlPfSrLf/pTvwXMe++9l+Xnn+931hk8uN9M9LOf7fdVmHfe/i6yfSussEKW//KXv2T51Vf7HWcWW2yxLM+YMSPLgwb1/1vFti2wwAKV/iy6aH9ES/bni1/st0578skns3zPPfdk+fHHH8/yJz/Zb+/Oti6yyCKNdW+44YZZvvPOO7P81lv9zn7vvvtuljkuLF+Shg4dmuXLLrssy7ffPmtrNs71UUcdleUTTzwxyzNnzmzDAP9vMGjQoI+1tvm5e4Y4+uijK39zncw/f38cmT//+c+N5f71r3/N8jzzzJPlxRdfvPHzhRfu92Xh2n7ttdey/Prr/Y5uq666apYfe+yxLH/00UdqAte8JH3605/OMvfJaqutluXllssxqXTDDTdk+YUX+p3KPvig3/mT5bD/yy7b7+y43nrrZfnaa/stHf/4xz9mecEFF8zyO++8k+X336eTZbV999/ff2Vx+eWXZ5nzy3L5+Te/+c0s/+hHP8ryBx980Li246QbCAQCLaLrRVpKqfHL+r96eD7L7pRw6KGHNj4jVU+xa6+9dpZfeaXfXZynBJ4k55tvvkaZbeXJePPNN88yTyE/+MEPssxTwsyZ/RemLJ8nA542pOrpe7PNNsvyuuuum2WeOPn8G2/0u9DzBMB/kZdeut81nqdYnoaWWmqpxv7whM1TMk8GkvSZz/Q7AJE1nHbaaY194JzyeZbDU8wLL7wwICddrm2uEc4nT5Lsl3vm2GOPzTLnow6ON8uizLXt9hVZCBna8OHDs8z5P+yww7L8zDP9ZrxkN1zbrKu+53n63HvvvbO85ZZbZvmRRx7J8ttvv51l9m3kyJFZvvjii7PMtc3+E1xTTz/9dJb5HuEarJ90uV/Z7x//+MdZZr851x9++GGWV1mlP6oA+/nKK6/ESTcQCAQGGvHSDQQCgRbRVb3gLhsIHrlJj3iJ8r3vfa/xtw888EDlb1IhUhDShZVWWinLpLWkFw8//HCWSZcpk76wLtIuKudJcUhZSDM4FlJV8U6aQnUB5Q022CDLY8aMyfISSyyR5aee6o8iSZXHo48+mmWqF9gGglSLz//+97+vPLfiijm4m/7whz9k+eab+70wH3roocY2kbKxD6SgF1100YCrF2qfNz7v1v+3v92fiYpjWr+Q5NpwFzJcVwT3Esf97rvvznL9ErepTK5Vto+/peqH+5lqvW7f8WKQ87/LLrtk+YgjjsjyrbfemmVSc84D9y3HgmuNqgJ3qcaLeql6EU158uTJWWb72GeWy7buvHN/bCi3tuOkGwgEAi0iXrqBQCDQIrra6brbSx6zHT064IADsvy73/0uy6QvvK2Xqjeizz2XE0VUVAFTpvTHuqBdKykFrRGc+oNU5sUXX2x8xllpOJBOSdWxIe2iFQEp+7PPPpvlG2/sj2NNtcsaa6zR2D7aXXKMSf2cjS8pG9spVVUvtP5Ycskls0zVAeeBNpj8fPTo0fqvANJG0sljjjkmy7R95XzU1wIpPPcJ54fjResQlsX17NQRQ4b0x5zn3uO+YFu5TqkiYZvraheWxXXF/cNnaJnwq1/9Ksu09x07dmyWuV641mjjS9Ucx8JZe9Daod5uziOtQqZOnZpl7g33Xtl66601K8RJNxAIBFpEvHQDgUCgRXRVL/QCHq1pfM3bRFIW3kSSlkre3Zc383QO4O0lj/6ky3ShpcMBn3fG4c4t06la+Ez9b0fznVqEfaM6gmNJywK2iW7Qjjo664U6BZs+fXpjHauvvnqWaVHBfvK2nH2rz/tAw6mRSF95K835cDSda7BbfZdc0p+ggaoAgvNDGkwnCLqNP/HEE1nmfHA9O0N/rn/ukbp7MMt1zh5UnbnxoBPJWWedlWWq0bbYYossc+1wHznrkm5hCrjHuFeHDRuWZb4/Xn755cbfcizqaqUmxEk3EAgEWkS8dAOBQKBF9KxecMd3UjAeuT/xiU9kmUd0Hr9JISTp+uuvz/KkSZOyTGpKGsAbZUffqaZg+0h33M0nKbijzaRT3QziSSlZn1NPODrGebjvvvuyTGsEqgRWXrk/GTHVEe5GnHEupOrccbzffPPNLC+zzDJZ5u0y+8N2M2Lb3ADOgTPKpwUJ55/rwqmT6uWecsops6ybEbA4phxHxs+gKoB1s/xebvj5OddgXXXG+rguaIFBtYOLxsc28Z3BvlG9uPHGG2d5o402yjLXMNcm90V9Thj5j/uBah7GhqCVE/vAfd/L2o6TbiAQCLSIeOkGAoFAi+g5tCPpMQ3mDzrooMbPecvKoztv60eNGlWpjwbhLrCyoyOkNVQjEOwrqQb7RnpFSkSqxWdILeoqAf6ecNYSLsyfUzXwdpjlkAqTpjLmAQNd8wa5HpKQxuv8jhSb43HppZdmmbf8HHuqPKZNmzbgsRc4/6SjjBFAVYOj446uStK5556bZaoLWDcpPG/pnVOPW7cufoSzTOCeclY69f5wj5Fqc62S8nPfO5WaU3+wHc7S5JBDDmn83KkFpepc07KJddMK6+STT27sA+ug5cPdd98dsRcCgUBgoBEv3UAgEGgRXa0XHGXhbSVVCqSfzgB+0003zTLjKEjVsIKkcM7PmfSd6gVHRxxdcuEZHfV3xuROrSFVaRQpCz9nWfyc5XLsaTnC8aLageEveau7zjrrZJkqorvuuqvSbpbLMV5zzTWzTKNxWkjULSFm9XmbcGub65njSHD+GReENLOeRYSZHUijKXM9s25n4O/UC6T7nH/W5dQaLp5DnZpzn3BdONWby7zhnEucVQctivhbqiYPPPDALDO+Sz1sKdvN/jmrjV133TXLEyZMaGwTY8Y4xEk3EAgEWkS8dAOBQKBF9Kxe4HF/3333zTItCHh8p6ph+eWXb3yGYdOk6g0n6RmpDNtE2s3PneE/b2xJCRzVdA4N7va6bkBOOCpIOkuqSgNtUh/2gZYFvNUmTWVbaVFC43O2u56BgLEYGHqSDhik1euvv36WOb+9pC1vE27OeQtOiu9oMPtOKn/GGWdU6nMJDh39dzfofJ7gb11iRtbr5tw5N9TXNsfAxV5gqFbGA7njjjsay3XJHpnNheNIVRvX/JVXXpnlbbfdNsuMYSFVY4awbo4Z3zHuHeBULQ5x0g0EAoEWES/dQCAQaBHx0g0EAoEW0VWn6/RQ1MPQ64i6HZqMMSAKM5gyg6/UW6oQmnPwedc+6ndd0A2atjiPHGcW0y1jKttE3RDLctlKXbxfts9lK6YOsJ5+pw802zrnnHOyPG7cuMpz9B5jgB0XtMdlWKWesR4YaCDgYiI7jy+OO3XvHHfqDOnhVIczS3Rru5fYrU6X6LzTuC7YNz7PZziX9edc3Vyfd955Z5adORj7w73g9r/bz7xPYpogmjNK0h577JHls88+O8suqBb10s6ctBfESTcQCARaRLx0A4FAoEV0VS+Qpu6www5ZJtUgtaB5EYOu0HSE5h+kH5LPLEx61UtQGJbrPGwcPXDBP0iJSQ+J+ucuqAgpC5/hWNLbjNS87hnUBwZk6SXOMOtlcKF6+YwPSu8emoyR5tGMjV5v9HT7uHRsToDrZa+99sqyMyVjXGFmp6VK4dRTT238reTXLdc855Dr03mkcW65RlzGYGcOxnXB551aS6qqP5wHl5OdioDtdm3iGDnPPv72hhtuyPKxxx5b6QNNILlnnCqQc7rhhhtmmSrTXtZ2nHQDgUCgRcRLNxAIBFpEV/UCj+kjRozIMqkPY4OSctBTjd5V06ZNy3KdgjnvMRcIhrfLrNupJpznjFNBOO80R/Hqt7gukAjHle126W1ckBPCxUol/WW9LIeUrVu59BKitxb7TauVlVZaKcv0QnI36m2C64Iek/Tacylw6BVGFQq9/Opw88nx4pg670bOGy0tuL6cSom/7cV6h6irnXqJO+1i6Lr6+DlVOG599hKcim1gsC1Juu666xrLpYqM64FrhuompwpxiJNuIBAItIh46QYCgUCL6HoWdkd/qhRILUgJ+Py9996bZd6y1p0JSOccBXXG1Dz607qC9N1lLSVVcMbxLusrn6lTMLaPMukLM4yyDnczy/64VEFsB58hpeRY05j8sssuq/SBAXmGDh3a2B+m/nnggQfUBFJql36oTXAsSOudk4lziDnuuOOy7JyJ6r/nnHDs+RvW7QI7OXrt9gL3KlUkbLej+Ny39fY59Rf3CcfMWemwbj7D8l0gLNcfjtf+++9f6cMuu+ySZRf/m+PHPcPUTM7ZwyFOuoFAINAi4qUbCAQCLaKreoGUhbSb9JX+zCNHjswyaebVV1/d+Nu6ITGP8qQa/I2jMqTgvLF0KhKCdMxluSXcTXT95pJ/O1rkUpq49CHO+NoZx1NlwX46FQlT+kjSBhts0Fg36ZWLAcGbfUdhBwouziz7SNWKc4i5//77s8z5q/fRqbxcOik+49Y8KTj3KlUnbv6pamNb2X+qPuqqQK4ffsf1MmnSpMY2sT6uW7bJqchcfApH8Z0jilS17KHqhW1ifXwHMsu2S/3lECfdQCAQaBHx0g0EAoEW0VW9wBs90rE11lgjyzx+00/fpQ/pltKGVIA3hb2EZKQzBmmAUxGQprF8F2PBZQzu5rjA/lC94IzgXUhKR69IWdk+toP1uizGpHW8sZaqFIw3xJxTOr+wPufr30tKkzkNR9Md3XfgOnXqLsnHDOglyzTnmWuHY+rUdpwzt9ZcOiz3vOTH75577mksy6kwnKqNqkpaP/EZxnR54oknskzVF59nmFKpqj5y88C94SwznOrJIU66gUAg0CLipRsIBAItoqt6gZl7mWmAlJpHdFJTUiJHp+r0zWWMdTEQnJE9j/ukB+7GljL741QCzpqAt8P1cp1xOOGsF5xxONUiLhOAk51xO7P5SlWHiEcffTTLDNFJA3dndcG10S2rQluoz1UfuHac7NRA3eg4n3O36y6co1N58Rk63NACxTnZEM5hiWuK8Smkqgpv7bXXzjJDeLKt3RxHmupmObRSoMqT2SX4OdUOL7/8cpZXX331Sn1UzziLEpf1l3PlMqI7xEk3EAgEWkS8dAOBQKBFdFUv8GjtfI1J00jNeaPN2AYudkAdPOKzXEdfSalIi5zhO2WXfG677bbLMhNwkk7wc0apr7eJdJ718XMXJ4J9ZsT6V199NctM8ulC7blo/25upeoa6CVhHykfx2bIkCFZrjtgDARcGE13i+1UZFxfXOd1FUrdML8JLJdZWEhZXewFrgXX1m233baxfO4XqimI0047rfI3+0P1BOveYostsjx16tQsO0uLXXfdtfHza665JssuXCqfZ2wY7i/KUtWa4bnnnsuyU+FwzbNultOL6ixOuoFAINAi4qUbCAQCLaKresElfiR9XWaZZbLsjPIZ+o9JEOvWCy5mAOsjXMg7F/6R9NoZbpMqOWN3fs6b0rpjAa0UqJJwqgMXpZ4U58EHH1QTXPwDlsMx5XySHtGnvF4f55S/d04QtGDh2HRzkGkLHCPOrbOacfLw4cOzzDgMdcsaqh54o+7GgioF5zTBfcHPCe5hxkA55phjstyLNU3dyYjrme07/PDDs8xEnVwXHG+2+6qrrmpsh3OscWoAPj9u3LgsM/OJJF188cVZZpYTF8bVOZ1QPReZIwKBQGAuQ7x0A4FAoEX0rF7gUZ60kTfUNJInhWZCOMZnqIdBc4kj3U0zqQafd1TGUUQ+Q1UI6+Ltu6PpjFMgVWkk1TC8aV566aWz7GIyUHVAmbfOhKNdpPhrrbVWltdbb70s19ULDNU3ceLELFMl8eyzz2aZRuq8EV5ttdWyTMP3gYJzFGGMBOfQwvkfPXp0lh966KHG8qXqnLukqC48oYvv0YvzAWnw1772tVn2xyWyPPLIIyv9OffccxvbeuGFF2Z5n332yfIVV1zRWIcLsUnVCdWCfGdwn/O3LqRoPTvEAQcckGUmqewlBCz7zPdhL1YqcdINBAKBFhEv3UAgEGgRPTtHuFzy7haTGSVIg6l2qPu/82/Sn16S8bmMCvycNIrlUC1CquBuLhmHgtSHqpb6b0gvH3/88cZnXPg/gm0ivXJqBD5PVcZSSy2VZY4151aqZgAhPbvvvvuyTL98hthjm+gQQco7UHBG72ybs5rh2uEcuNB/klcXOEceZ43AclwsEZZJ9df06dOzTJUf2+1iSZxwwgmVdrAPTM555plnZnnChAmN7XZqBJeFhb9lf9i+3XffPct0puKa55qVqlYHLh5EL7ExXNwThzjpBgKBQIuIl24gEAi0iK7qBdJxxhUYNWpUlulM4CKuk+5uvPHGWb799ttt3S68GumFizpPGkAaTRrA23uqFzbaaKMs05rAGUyzLlo41NvH+hgi0d1YE6yb40r/efZ/5513biyHagCqSAhaXEhVdQP7QysPUkQ6i9DZg6qXjxtpf06Aa4o317TWoNqFlJXriPO37777ZvmXv/xlpT7OG/dML5lN+AxVPHyGqqmxY8c2tnXEiBGNnxMulOGBBx5YeW78+PFZZqwD1x+uT46Fy67CsefzbAfVA4z5wTYwuSZVKlJVdUa4rBDOWopjGdYLgUAgMJchXrqBQCDQIrqqF3jMZnI4hmwjxaF6gZSD6gXS95tvvrlSH4/mlOl04ELvEQxPx/pIRwhSaNJrl1yTFhu0dqCaQvKUx2UhoCqEKgU6LNDwm77kHBc+QxUOVSHDhg3LMp0b6iA9o4OHCxP51FNPZZkhH5mwcG6IvcA5mDJlSpaZOYNz5mKBuFgTdfrusohwTnrJTsE1v9dee2WZa4/ls61c56Tsbj7Yf2YQkaSvf/3rjfVxXbDdTqa65Mtf/nKW6XzjknROnjy58XPGUaAqj9Ybkg/dSQsW7kmCqhCq6rgeHOKkGwgEAi0iXrqBQCDQInqOvcCjOOk16Qs/d0bm3SiYC/nGIz5D3rlwi6TOW221lZrAm0hmXaBFAKm1S0RHOlXvD2/yL7300iw7tQgpFak81Qs77rhjlqdNm5Zl3q67bA90WKH1AdU/pGOST8JJywSn/uEzHJu6D/xAgOPCtcO+UKVEcHyZoJGUs+7c4KxUuLbpQMLPaUWz0047Zdk58hDsm1NxuOeJulUC+3rOOedkmWuJe4zqD4a5PProo7PMfUjVHNvKvcB9/rnPfS7LVCPstttuWWZoS6k6fs66gp+7eBVEOEcEAoHAXIZ46QYCgUCL6KpecBH1SZV32WWXLJM2k+6SjjNxW53K8CjPYzpVCi6OA2khjcCZINOFhaO1A60uSPF48+tiCrA9kk/C6SgLHQjWXXfdLDMsIimYyxZBH3Oqdvhb0iPS13rGA1pzkFIOHjw4y6ThHEuOMWln/RZ5IODCX9KihnPAsWZfuO6cw4HkaeeLL76YZY491RZcLxx3wmWXcODzHAuuIxd7RKrub8ZMcaEhaf1ACxGuSZd0knBxO84+++wsjxkzJsuMBVIfF7en2R8XD4NtdWPmECfdQCAQaBHx0g0EAoEW0bP1Ao/7NIAnPSZVIiXmEZ2Uu07HSYVJf/gbgrfgjOnAI75TWVAtwJgEpBak0yyHFg58pp51geO06qqrZtnd8LP/vO0lpWLdvAXmWFItQsNtzo+zKGGMCKlqUUFVDevg5wxb6TKMzA3OEc4Rh44SjMPAdUE1AFUKLtRm/W/SVGddQMrKmA5cO5SdxQrLd8+7xKJOBVMvi/EQaMnA35Dmc686iyWCn7OttC7iWmOcGKpj6tkvvv/972eZ2VM4v4zP4BLKujgZDnHSDQQCgRYRL91AIBBoEbNO0l7CJUrk56S7vBEnzeCtZ/1G1CWEY7mkI7RS4NGftIi39KQapLsuEwZv5UlrXALC+s0q1R9UCzjwJpxlsf9UF3Asnb/4mmuumWWG/xs5cmSW2Z+6ZQHVMATVIqR57DPnirfXpJoDBa49UlOqETi+tCZwRv+knHWa2Uu2CMp77713Y1tdQlXn7MDnXfhTt+9ccth6uXxu1113zfKvf/3rLFOd5eKNMLuKC//Jz+lARKsoqi/oZHH66adXyqJa4Mknn2xsK9UtTi3G9eOcJir1zvKJQCAQCMw2xEs3EAgEWkS8dAOBQKBF9Gwy5rL4UldJEylnVkWPr7rehnoi6pzoAbTJJptk2Zn9UOdGvQ3rpk6SXlQulQpjw1InyRi9dT0ex4PBbxjj98EHH8wydct8nv2njpU6R9ZF3TCD3LCtHC/qxurpeugNSN0Vx4nlUh9GvbnzHhwocK6cSaLTvXKe6jGUm34rVfWv1KHyc5qGuVRUXLfUh7pANdRDcg6cNyN1kk4HLFX7x/3N/hx++OGNn7MOmhjyzsV5hbns3scff3yWr7/++ixzr02cOLHSB65njofzqiPYf965uGBWRJx0A4FAoEXESzcQCARaRM8mY6SNBI/ipKx8njTAechIVfMMHt+pqmDaF9Jrmo/xc3pLsRzGk3UZPEmpSOvosUKaXR8jUjv2h2XxGXrGkbKTgjmzFddWmsyQ1lOmKVmdFpM+kyK6NCacX847TdpY30DBzTnbz3nm3DoPJI5PPWAL5410nCl3zj///CxTpbTddttl2WWnpamj85Zy2WxZJn/rzNCk6rrn2qMKg3W74DKs47nnnmt8hvuW5TPN1BFHHJFlBiq64447stxtTty4EnyGfaMqJDzSAoFAYC5DvHQDgUCgRfSsXiDcUZwyb/QYNILWB6RQUvVozrJ4o87bdFINBmrhbTypLG+pXVup4uDz9LqiJQM9YeoBfNyNNWWOBykLrRRIX0hzSc2ojqCaghYEbCtl0jeqI+r1USXz2GOPZZmBftgmjhmf6cVrZ06D9JXrjrKLj0uvO66vbgFvGDznkUceyfIvfvGLLDuaz/Gq75kmsG6uKa5Plsl+Uu4WJ5ZtZb+d+pBrmKoQqlq41lx8W8qrr756lhlDmp5tLHOzzTar9GHq1KlZpkqSfaWKiXuJ6jm+J5y1AxEn3UAgEGgR8dINBAKBFtFVveCcD3jk5tF6lVVWyTIDhNDQn+XUDcsZ79XFGaVBs4s/y2A7jGNL2sEbVxeLlOoBUijempKmkSrVf++cK2666aYsk6puu+22WWawGKpFOPZrr712lhnA5rbbbmtsH2nq9ttv39gfqTrevF1mm0hJncqDFhW8yR4ouFtmtp/9crTR3dzvs88+ledomcD1RtXWSiutlOUZM2ZkmWuPYB+cxQLh4uP2otaol+nq4PhxPS+99NKNdbMP7CfbwedZJtVXtKahOpLvpIcffrjSVqrL+L5iv7lWOb/OsiPS9QQCgcBchnjpBgKBQIvoql5wPs/Dhg3rLwAU2sWi5bGcNLNOUUgd+BvSaNJiqh1YN+mFsyCgioS3yewbrSAYt8GpOOqxZ3uJ00vLDlJVOnWQIvFzUrNbbrkly3Sm4PMcbxqQU21Qzzbr0sNwbVCd4+IKcG5dptc2wbEgJRw+fHiWqV6gxQL7SCrKZ04++eRKfVzPnH/S2nvvvTfLHGtSbY6dM+538XGdU4JLJeSel6r94Thxr5533nlZ5vo/5JBDssyxcDF0Caogr7nmmixzjI4++ugsU6W4ww47VMri2NNqh2PsYm4QtGQI54hAIBCYyxAv3UAgEGgRXXmeC7U4atSoLNNYmUd53pQ7SwamnpGqNIXUhlSGN5ak9meddVaWx40bl2VSfoY5ZDlMaUP1B+tl30hNSUc5FpI0ZMiQLE+bNq2xXI4Nx5hjw9CT66yzTuNvabj9zDPPZJmWFhxTPk91ydixYyt9YF/5HMeJahhH212ovoECx4Iqgt133z3LXCMuNQ77wvmoU3xanXD+naqC43XJJZdkmdTcjbWzLHDxJtg2toH9rFtQuBgbVC849ZLLhuziPvC9wNgrzpHjlFNOyfL48eOzXE8TRYsKZyHisonztx8303WcdAOBQKBFxEs3EAgEWkRX9YK7HeVtHSkLbyhJM2i4TFpa920nRXDqBd4g9hI6jn24//77s0yKs9pqq2WZaoSnn346y6SOjP9Ap4T6DS+tKG6++eYsjxkzJsukdmwTfclpQM92kILSAoFUjpTNWRYw/N1PfvKTSh9IbZ0BPueHMR2o/qGFQy/+6XMaLgMB162j487vnvPXLSwi542/4eecH6rIOLfcb87SyKk/SKdJjzlnzpJDqq4lZtk97LDDGvvDMfvhD3+Y5W984xtZ5l7lu4BqAWft4GI+HHTQQVmeMGFCpQ9U4bEsl7nYhXTl/ETshUAgEJjLEC/dQCAQaBE9x14gnaBlgqP1NL5/9NFHs8xjef2mz90gEqQOLnr7tdde2/hbGq+TFl5++eVZZog3OkSw/6RKtCaoG/2TsvA7qhQYM4IhFmnlQYrPMSINogrChZvkGHEcSTtZl+R94Pk5DdY5ruzDpEmTGp8ZKHBu2B5H8XsJ+ci9UF/LHC/OCeeQv3EZWbhWXdwTqsW4V2m9Q3WXU/+tscYaje2XqnuXqiNmCKHqhO2jNc5JJ52UZa5/5xDVS1YMWtnweVopSdXwpHyOfXNJKrlv6YwUiSkDgUBgLkO8dAOBQKBFdFUvkGqTjpJGkZpSdeDiFjj/cqlKkUkvSOd43KeBMm8QaYHgbm9ZJikB63W+1mwDnyFtkqrjwdCODHnI0JOUXbudesVF6Sd1dJYMvL2mk4kk7bTTTlkm/WUdpMucE1JyOhrUQ2AOBFzGAo6jo/gOV1xxRZbrNJPzyfmn8wp/Q4rsYilw7ZHict9S1eSclwhSee4RtlOqqgv23HPPxj5wzTNsJdUcV199dZa5Rthn9+7heqRMVQPXYD32ApOC8jmqSKiS5DuGa56WQ/XsMU2Ik24gEAi0iHjpBgKBQIvoql4g1ebNImkKVQJUKZCO8Bkexeu3wLwF5Q0in+ONKGkAn3dG1oTzMXfG7rzJZj+nTJmS5dGjR1fqoPE6k+hRDUE6wkR5vEUm2CZSOdIu9pnUmeoIZ03B+ZF89gCqLWiw7xxTSK978U+f06Dqh/NBWsuxdo47BMeOv5Wq68rFw+AYuVgFVAW5jA/cny4UJA39qSpi8kaqmuqZMFgu1WIcP5bFfXLVVVc1lsMxZjnc/26fO4cQFy61/hzVYnwPOQsu5/zSzSmmD3HSDQQCgRYRL91AIBBoEV3VCzwqL7fcclkmHXOWBTzWk8q4qOz1v11iyuWXXz7LPPqTUpGyk7K4WA0uASFl3o5usskmWaaBNW+KpSrtonqG1hzPP/98lp966qksO6N5jiWfoYrAJSDknPAZV6ZUTea31lprZZm335xT3qKTwruwhQMFrguqi7h2ejF0J7qFV3Rjz/HmmDIeCNcYx5pjSocdxhRwmT+4XlgmY4k4Z5r63y7E5AUXXJBlp3biePNzjlcvGTw4jizTJcGUqvNLZ4deknY++eSTWXYJNR3ipBsIBAItIl66gUAg0CJ6zhxBA33SABpTu+R4VAl0oyyk/M66wN3M8lhP2d1Ak165W31HA+m3znJ4s1xvBw3iGbaRCfFYh/P1p+xiO/C3pFruRpzP04hdqt5AUw1B/36GsORcM+aGU+HMDXC30hwXjq+zIHFODPXfuKStpLVU0zj6SpkqBaeCcJYMtC6iSoD1PvTQQ5X+HHrooVnmOmSiVTcenH9adTiV4sfNLsGx5vupDs6vy2zC/UnVhtu3vSRdnbtWfyAQCPyTI166gUAg0CK6noV5nKbBMI/iPPrzeVLOu+++O8ukBHUjc2c5wOec7NrBoz+pCakFn+ctsKMKpHJOBVMH+016yd/zBtVRH2coTprq1DnOIoLjuNFGG1XaTZUCb9Fp1bDCCitk2VFHl0RwoMCx5ryx/Wwn1w4pO8G1U4/bwTF2cR+IXgz/XSJXPs9yGJOBbeBaIG1m/Ie6s8cZZ5yRZY4Hy6WVk4sf4sLEOgcHrke2m887teARRxxR6QPHj9lgCPaN9XG/sU3hHBEIBAJzGeKlGwgEAi2iq3qBNJhHaNJxPkMKNmTIkCw7q4S6eoG/J+1wkdz5OandF77whca2Tps2Lct0SqAagVYGDE3329/+trGdBCPtS1Vfb9IcylQFuLHsJeEhx9IlhHRZEThGNMqvl0vrDN4K0/GDtMvdwJPaDxQ4Lu5GnDItP5x1gLsBl7y6gPuKzzhrlEMOOSTLHFOqfqZPn55ll9R1++23z/L48eMby6Szz1ZbbVXpDy14SNO59qjOcHuV4RKduoB73rXP7QtaF7FeqRo/5NZbb238vZurG264IcvsM1UyDnHSDQQCgRYRL91AIBBoEV3VC6RdPO6TBtNKgTSAlNMlUKyHWnMUgXU7v31SsMsuuyzLjuIQbPf666+fZcZCIDVzfut16wVSp15CzxEuPCWfJ9Vi+1z2A2eNMWLECPvMXXfdlWXOO6kZx4kWC3RkYR/qNG8gwDlwFjXuGYKfM5RnfS2w/5xb7jGX8JK/PfPMMxvrdiooqul23nnnLFM94JLOci5vu+22Spvc2FCdMWzYsCzT2qUXaxFnmcH3By0OXB/23nvvxnolacaMGY318Z3BdxfVa/ycv62r55oQJ91AIBBoEfHSDQQCgRYxa0fhEqREjLLOoz9v93jEZ9YEHsUZQV6qJqajFQFpuksi6MLCORrAZ6gKIb3mjbCjO6yXThOSN4Jn3fyNi5bvQmM6RxEX54C0iVYT7H/dMoN/kzouu+yyWeYN/m9+85ssc2zYvlGjRmmg4bJaONpMmWNCFc9RRx2V5e9+97uV+jh2DJ94++23Z9mFM+Q6crE3CJZDmXuK1jucG5eBpK6a41pyCV/r8Rr6wLHk+LkwrFRz0CLGWWasvPLKWabKh++Xen3O0YrjR3WE6wNjlTjESTcQCARaRLx0A4FAoEXESzcQCARaRM/xdKnnoD6DHhjUT/JzZu1lapR67Fb+nmZILo2PC7TBACz0JKJ5G3+74YYbZvnBBx/MMvVE1Ks6rziakknVvlLvw/GjCRzLdeYz1FFxXHvxnmJ76D3HMWXAE6nqrcYMxTSf4dg/88wzWXZZc2m6NFBwmWdpCsf+si80nWMaK+rtd9hhh0p9NGNkACiXTsbF7+Xa5jPMRMx27L///lk+7bTTssw+c02xby6ITr0O7g16G7ps3UQvKXe4d9zdDd83zMrNPVkPtOQ8ALnveXfBFD18nmNDk1OHOOkGAoFAi4iXbiAQCLSIruoFZ4Y0adKkLG+55ZZZprcIafOLL76YZWaUHTlyZKU+Zsl94oknsuyChfCIz/aRBjizKtIOBrMhHIVi30hr6hSKdfA3pE4u5Y5LXUJvLpZJqsV+0ryPQYjWW2+9LHcLHEKKSDMh1s3gMVRtODMcmisNFEg1KTNYDMea4+vWNk2VaCYpVdVZTo1EUzRSftJdqm9cuiquw9NPP11N4LrjHJPKu3bWn+Peo/mYM6106oI111wzy1QFcrxd7NqDDz5YTXBee1LVFI1t4r7nWqUHLcvi2u4WUzs/P8snAoFAIDDbEC/dQCAQaBHJ3TBL0nzzzZe/JNVgIIsxY8ZkmRSHR27SbN4s1rPn8mhOWuRS1JDau9QwLn6vs8zgDT+pI2/xnRceb34lH4TGpXtxwXlc0BqndmBbOd7s88Ybb5xl3nwzg69UtUwgfebc0dKEgVFcWifWfe211w5I7p555503TyLHmv396le/mmWnKuEc0OOr7i02dOjQLDOTrstE7FRNLtiSixvLdtMyg8FinKqA6XaogqrX57z4SMFdIB2CajfOCfvPZ8aNG9f4OeviO4V7Vap6qHHM+B7jO+DKK6/MMueNda+77rpZvummmxrXdpx0A4FAoEXESzcQCARaRFfrBWfQzBtU3vQxhQuNpGmVwPispFxSlRazLFIEUlbSZUfZefQnJSKl4k0xHTaoLqCqwFF/yt3a525K2TfSQrbP9ZNjSZm/5fg+/fTTWebtdf2WmhTW3XjfeOONs+wDnTGmTp3a2Ic24dY21wjXoIuNTOrL+MF0YpCqVhHOosYFsHE3/y69FVUE3IekvldffXVjGzjfVCnUY/06NRf3G8fGZVN2ViSUP//5z2eZwZk4Fuw/6+XndHbp1gfWzQBOzhqDab2mTJmiWSFOuoFAINAi4qUbCAQCLaLneLoEqQLjW2699dZZpsE8weM+DaClqp8zM/oy7i7pK29BeVNKiuNibpIu1W9m+0BrAvbHpRghfZeqN/xOvUCaznIpk1KRLpLW8HmqSHhDyzZQvcAx5RxIVcrLG2+Wy3FlO/jbXXbZJcsnnXSSBhrOeoW46KKLsrzNNts0PuPSMNUtWeiYcuSRR2b51FNPzbKjr46m0wHjkUceaWwf5+aWW27JMtdCL5YClKXquuV3Lgaxs8zgPPDd8JWvfKXxt7RYco4czhGlDqqPOPa04HHxV1g3LbjcPBBx0g0EAoEWES/dQCAQaBFd1QsufQjBdCM0enc3sYwFUPfzJ4VnShP68zPuA9ULzijbhYJz4Rld2g9Sa9e3XtOBsD6qC6jOoJPGFltskeUdd9wxy7wtZ6wKGumT5nIOBw8enGXSsbrvOMeAIQknT56cZedg4xwr3Pi1CbbBpRViDA9mleUtOMvhmuKalaoWP6SppNHjx4/PsnMUIK2lU4rLGO1Cj1J9d91112XZpaWq03SnknGZdJ2Fw4QJExrLYaofhhvlb6nioKqA6jKXqViq9o974Prrr2/8nGPMeede6gVx0g0EAoEWES/dQCAQaBFd1QsulgKP7KQj55xzTpYZsZ7HeBc6T6pSW6oXaNRNY/+JEydmmVSG7e7F/9tlPaWqwfmXk8rUM+m6rKykOYxjsd9++zWWxXEivZw2bVqWSd84XvQdp/qDt8C0ZKhTMMZlIJ2l2oJj4zIycK6cY0KbcDf2zvHlxBNPzDJjMrhy6k4mboy4nr/1rW9lmVYNbIfLYEI1FfeCy0ZBtaBTFTi1i+RDiXLdcm3/7Gc/yzKzaNx8881ZHj58eJa5Jp1VA+OzcD06K6C6ioTjwZghVMO4rNF871E90y2WTW7TLJ8IBAKBwGxDvHQDgUCgRXQN7Tho0KD8pbulrT2f5S996UtZXmWVVbLMm996HAF3O0jDcv6e9VEVcN555zW2j6oAUhYXvZ9tIOXgb13UfalK1RlikTfHtNhYccUVs8wxu/DCC7O8/PLLN/52xIgRje0jjXR+/iyHYyRVrSh408ybfZe1YPPNN8/yww8/nGWO2YcffjggoR25tnuh11wjBx54YJbr4QL7ULfQoEqKN/wu+4P7/Oc//3ljmxxc+EenLmR/qKbohu233z7LTMjJRLNck6NGjcry+eefn2WqvNg+rqMZM2ZkubaOskxVA9dzvT/cS1R5UAXq1GW01OLartUXoR0DgUBgoBEv3UAgEGgRXdULgUAgEJi9iJNuIBAItIh46QYCgUCLiJduIBAItIh46QYCgUCLiJduIBAItIh46QYCgUCL+P9vo0XO/ZgbhAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(x_et.reshape(size, size), cmap='gray')\n",
    "plt.title('Reconstructed Image')\n",
    "plt.axis('off')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(x_flat.reshape(size, size), cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.title('Original Image')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b8b4c69-2926-4516-9df3-117e4b419b3a",
   "metadata": {},
   "source": [
    "# functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8e5cd8ca-b79f-4b64-89f8-e059c6b3e203",
   "metadata": {},
   "outputs": [],
   "source": [
    "def intialize(y):\n",
    "    x_start = np.zeros(M)\n",
    "    for k, o in enumerate(y):\n",
    "        aver = o/sum(A[k])\n",
    "        x_start += A[k]*aver\n",
    "        # print(k, o, aver)\n",
    "    x_start = x_start/N\n",
    "    return x_start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d4c989c2-f8a1-4dc0-8ab3-d04d81039bfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def em(A, y, x_old):\n",
    "    z = np.zeros((N, M))\n",
    "    for i in range(N):\n",
    "        z[i] = y[i]*A[i]*x_old/sum(A[i]*x_old)\n",
    "    x_new = np.zeros(M)\n",
    "    for j in range(M):\n",
    "        x_new[j] = sum(z[:, j])/(sum(A[:, j]))\n",
    "    return x_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "44215791-168f-4c65-b817-6845645d5356",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mle_em(max_iter, A, y, x_true=x_flat):\n",
    "    x_old = intialize(y)\n",
    "    mse = []\n",
    "    for i in range(max_iter):\n",
    "        x_new = em(A, y, x_old)\n",
    "        mse = np.linalg.norm(x_new-x_true)\n",
    "        diff = np.linalg.norm(x_new-x_old)\n",
    "        if i%1 == 0:\n",
    "            print(f'step: {i}, diff: {diff}, mse: {mse}')\n",
    "        if diff < 1e-5:\n",
    "            return x_new, diff, mse, i\n",
    "        x_old = x_new\n",
    "    return x_new, diff, mse, max_iter"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nmf",
   "language": "python",
   "name": "nmf"
  },
  "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}