{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "91582ce2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4f6b472",
   "metadata": {},
   "source": [
    "# Gaussian case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "88895113",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_Y(d, lambd):\n",
    "    G = np.random.normal(0, 1, (d, d))\n",
    "    Z = (G + G.T) / np.sqrt(2)\n",
    "\n",
    "    x0 = np.random.normal(0, 1, d)\n",
    "    Y = np.sqrt(lambd/d) * np.outer(x0, x0) + Z\n",
    "    return Y, x0\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "65fc3267",
   "metadata": {},
   "outputs": [],
   "source": [
    "def AMP_Gaussian(d, lambd, T=500):\n",
    "    Y, x0 = generate_Y(d, lambd)\n",
    "\n",
    "    b = 0\n",
    "    x_est = np.random.normal(0, 1, d)\n",
    "    x_prev = np.zeros(d)\n",
    "    m_AMP = []\n",
    "    for t in range(T):\n",
    "        u = np.sqrt(lambd/d) * Y @ x_est - b * x_prev\n",
    "        sigma = lambd * x_est @ x_est / d\n",
    "        x_prev = x_est\n",
    "        x_est = u / (sigma + 1)\n",
    "        b = 1 / (sigma + 1)\n",
    "\n",
    "        m_AMP.append(np.abs(x_est @ x0) / d)\n",
    "        if t > 3 and np.linalg.norm(x_est - x_prev) < 1e-8:\n",
    "            break\n",
    "    \n",
    "\n",
    "    return m_AMP[-1], x_est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "72fd7783",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "16it [00:04,  3.34it/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG1CAYAAAAFuNXgAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAP8JJREFUeJzt3Ql4VNXBxvE3C0lIQsISSEgI+xpZBUFUCioWXEBtXapWEKvWqq2KrUpbRWsraJWilRYVLfLVDa1aF8QKgqigKIiCbLJDyMIiZIEsZOZ7zgmJCQRMwmTuzJ3/73nmmbl37syc3EsyL2cN83q9XgEAALhEuNMFAAAA8CXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcJVIhRiPx6OdO3eqSZMmCgsLc7o4AACgFsycw/n5+UpNTVV4+PHrZkIu3Jhgk56e7nQxAABAPWzfvl1t2rQ57jEhF25MjU3FyUlISHC6OAAAoBby8vJs5UTF9/jxhFy4qWiKMsGGcAMAQHCpTZcSOhQDAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAABXIdwAAADfKCmU7kssv5nHDiHcAADQUIryvv+y//Z9yVPGufYDwg0AAA1h9ZvStIHfbz9/iTS1Z/l+NCjCDQAAvmYCzOwxUn5W9f15WeX7CTgNinADAIAvmaanuXdJ8tbw5OF9c+92ZxOVp8rPtHWxYz8j4QYAAF8yX+p5O49zgFfKyyw/zk1WB04zHOEGAABfKsjx7XHBYHVgNcMRbgAA8KX4ZN8eF+g8gdcMR7gBAMCX2p0mJaRKCjvGAWFSQlr5cW6wNfCa4Qg3AAD49Js1Qhr50OGNIwPO4e2Rk8uPc8NkegWB1wxHuAEAwNcyRkuXzZKapFTfb2p0zH7zvFvEB14zXKTfPgkAgFBiAkz388ubY0ythflyN01RDVlj42QzXF7WMfrdmGa4VL82wxFuAABoKCbIdBgSGs1ws8ccbnbz+r8Z7sgi+e2TAACAO2UEVjOco+Fm0aJFGjVqlFJTUxUWFqY33njjB1+zcOFCnXzyyYqOjlbnzp01c+ZMv5QVAAAchwkwNy/9fvuqV6XbVjrSv8jRcFNYWKg+ffpo2rRptTp+8+bNOv/883XmmWdqxYoVuu2223Tdddfpvffea/CyAgCAH1C16cnB/kWO9rk599xz7a22pk+frg4dOujRRx+12z169NDHH3+sv/3tbxoxYkQDlhQAAASLoOpQvGTJEg0fPrzaPhNqTA3OsRQXF9tbhby8vAYtIwAAISsqTrpvv9OlCK4OxdnZ2UpOrj5O3mybwHLw4MEaXzNp0iQlJiZW3tLT0/1UWgBAwPDnpHZwXFCFm/qYMGGC9u/fX3nbvn2700UCAKBheaqs42Tm2fHTuk5FpWXKySvSzn01Vzj4S1A1S6WkpCgnp/r0zWY7ISFBjRs3rvE1ZlSVuQEAEBJWvym9e+f3289fUj4k28xFU4eRSyWHPNpbWKI9hcXaU1D1vkR7C0q072CJ9h0o1f6DpfbebBeVeuxrT+3YXC/dMFhOCapwM3jwYM2ZM6favvfff9/uBwAg5JlgM3vM0TMFm9mDzf7LZim/47m2diV7f7Gy84oOPy5Sbn6RDS8m0OwuKFZe0aF6nc7wMMlT00TFoRJuCgoKtGHDhmpDvc0Q7+bNm6tt27a2SSkzM1OzZs2yz99444164okndOedd+raa6/VBx98oNmzZ+udd95x8KcAACAAeMrkffcuG2yOXo/cK1OnkvPybTq9OEyeWvZKiQgPU/O4KLWIi1JSfLRaxJvH0Woe10hNY6PUNLaRmjYuv09s3EiJsY0UHxWpcJNwQjXcfPHFF3bOmgrjx4+392PHjrWT82VlZWnbtm2Vz5th4CbI3H777XrsscfUpk0bzZgxg2HgAICQUVB8SJt3FWrT7gJt3XNA2/ce0PbvDihp1+d6onTnMV9n4kzrsD0aGL5W30T1VkpCjFISY5Rs7hNi1CohujzAmDATbx5HKSGmkeNBJejCzbBhw+T1HrvuqqbZh81rvvzyywYuGQAAzjHfjZn7DmqjCTG7CrRxV4E27Sq09zl5309vUtXo8Bwp6off+7lL2yq6n7vnhguqPjcAALhNYfEhrcvJ15qsPK3NOnyfnW9raI7F1Kp0TIpXuxaxSm9ubo2VURQmvffED35edNNUuR3hBgAAP8kvKtXKzP36avt+fb1jnw0yW/ceUE2NGI0iwtS+RZw6toxTp5bx6mhvceqUFG/7thzF01paklreefjIDsVWWPmoKbMsgssRbgAAoTfvS6ezGnzdo0NlHlsDs2L7Pnv7avs+bdhVUGOQadUkWt1bJ6hH6ybqkWLuE2yQaRRRh+nozM9jhnvPNqOlTD+Zqh90uN/MyMmOrffkT4QbAIC7+Wjelx9SfKhMK3fs12eb99rbsi17VVhy9OR5aU0bq296U/VJT9RJqYnqntLEduD1CfPzXDar/OfNNzU4h9mfd7IjK3Q7Icx7vB69LmSWajDLMJjZis3kfwCAEJz3paImwwSBen7hl3m8+mrHPi1av0ufbdqr5du+U/Gh8knsKjSJibRBpuLWu01TtWzih4lli/KkyYeXG7rqVb/UVAXS9zc1NwAA9zZFzS2f9+VoZl+YNPduqfv5tf7iz9p/0IaZRet36+MNu+3svFWZYdQDOzTXoA7NNbBDC3VLaWLnivG78Co/j+ljE+TBpq4INwAAdzJ9a/KOPe+LDTh5meXHdRhS8xFer77esV//W52t91fnaH1OQbXnE2IidUaXJJ3eOUmDOrRQp5ZxCgsLvnlh3IZwAwBwp4Kceh1XWubR0s179d432frfNzl2iYIKphKmT3pT/ahLS/2oa0v1aZOoyLp0+oVfEG4AAO4Un1zr4zwer+0E/N8VmXp3VXa15qa4qAgN69ZKPz4pWcO6tqp5GDYCCuEGAOC/PjCmCcjUlJjg0dB9Qcz7Jxx73hezAlNpXIoeXd1M/33xg2o1NKbvzDkZyTbQnNYpSTGNQqvPSrAj3ABAqCkplB48PEvt73dKUXH+GbVkOvdW7QPTAMOxazvvS/nSkl79+rvL9d5HW8uLExOp83q11ui+qbb/jCMdgeEThBsAgDPDsU2Nitl/AsOxazvvi/fdOxVWZd6XLG9z3V96tRaEn6pzM1rpwr5pOrN7S0VHUkPjBoQbAEBQDceu67pNr+7rrVc8U/S2rrD7xpTcqb2tTtflp3bQw31SldiYPjRuQ7gBAAT0cOz6MPPRzFy8RS9+tk15RYfUWEVSTPlzv7t+nHq2b82QbRcj3AAAAm44dn2tytyvpxZt0pyVWTrkKa8tat8iVjec2kGaX35Mr7REibloXI1wAwChxp+LSNZhOPaJMAtT/n3+t5q/Nrdy36kdm+u6MzrqrO6tFH7oQGW4CQlRcdJ9+xWqCDcAEEr8tIhkbYdj2z435nlzXD18sWWvHv9gg10SwTADnEb1SdX1Qzqqp6mhQUgi3ABAqHBi1NJxhmNXLl5pVquuY83R6p15mjx3bWWoMcO2L+6XppvP7KwOSX4Y2o6ARrgBgFDg5Kilw8Oxa57nZnKdAtWO7w5oyv/W6/UVmfJ6pcjwMF06oI1+NbSz2raI9W25EbQINwAQChwatVTJBBgTnOo5Q7FZDuGJD77Vc0u2quSQx+67oHdr/W5EN7VrQU0NqiPcAEAo8POopRqZIFPH4GTWfHp12Q49NHet9hSW2H2DO7bQhPO6q3ebpg1UUAQ7wg0AhAI/jVrypZU79uveN1fpy2377HbnVvH6w/k9NKxrS+aowXERbgAgFDTwqCVfyi8q1cNz1+nfn221/WrMqty3De+qa05vr0YR4fV70xAfGh1qCDcAEAoaaNSSry1Ym6vfv75SWfvLV+g2I6AmnNtdrRIOTy8M1ALhBgBCRcWoJTPPTZVFJOszasnXviss0QNvr9ZrX2ba7XYtYjXpJ710Wqckx8qE4EW4AYBQYgJMx2HS5PTy7atebdgZimvhg7U5uvPVldpdUGwn4bv29A6648fd1DiKFbpRP4QbAAg1VYNMHYZj+1pRaZkenLNGs5ZstdtdWsXr4Ut6q1/bZo6UB+5BuAEA+N2arDz95sUv9W1ugd3+xRkd7Jw1MY2orcGJI9wAAPzG6/Xampq/vLNGJWUeJcVH69HL+mho15ZcBfgM4QYAQo1Dw6IPlBzShNdW6r8rymdKPrt7K9sM1SI+2u9lgbsRbgAADW7L7kLd+O9lWpudbxe5NMO7TVNUWNjhYeiADxFuAAAN6v3VORo/e4Xyiw7ZZqhpV/bToI4tOOtoMIQbAECD9a/554cb7WzDxoB2zTTtqpOVzIR8aGCEGwCAz5mVu//w+kq9smyH3R4zuJ3uuSCj/ssnAHVAuAEA+NS+AyW2f82nm/baSfkmjjpJY09rz1mG3xBuAAA+s+O7Axrz7FJt2lWo+OhI/f3KfjqzWyvOMPyKcAMA8In1Ofka88xSZecVKa1pYz1zzQB1T0ng7MLvCDcAgBP25bbvNG7m59p3oFRdk+M169pBSklkJW84g3ADADghi9bvsn1sDpSUqW96U80cd4qaxkZxVuEYwg0AoN7mr8nRr/693C6lMKRLkqb/vL/iovlqgbP4FwgAOOFgc27PFE39WV9FR7LwJZxHuAEAnFCwOa9Xih77WT/msEHAYDYlAECdEGwQ6Ag3AIBaW7xhNzU2CHiEGwBArXy9Y5+un/WFbYoacVIyTVEIWIQbAMAP2pBboGv+9bkKS8p0WqcWevwK+tggcBFuAADHtXPfQY155jPtLSxR7zaJemrMAEZFIaARbgAAx7T/QKldK2rn/iJ1bBmnf11zil0zCghkhBsAQI1Kyzz61fPLbJNUSkKM/u8Xg9QiPpqzhYBHuAEAHMXr9eqeN1Zp8cY9iouK0LPXnGIXwwSCAeEGAHCUpz/apJc+367wMNnOwxmprO6N4EG4AQBU89432Zr07lr7+I/nZ+jsHsmcIQQVx8PNtGnT1L59e8XExGjQoEFaunTpcY+fOnWqunXrpsaNGys9PV233367ioqK/FZeAHCzddn5uu2lFfJ6patPbadxp7d3ukhAcIWbl19+WePHj9fEiRO1fPly9enTRyNGjFBubm6Nx7/wwgu6++677fFr1qzRM888Y9/j97//vd/LDgBuk1dUqhv/vUwHS8t0eucWmjgqQ2FhYU4XCwiucDNlyhRdf/31GjdunDIyMjR9+nTFxsbq2WefrfH4xYsX6/TTT9eVV15pa3t+/OMf64orrvjB2h4AwPF5PF7dMfsrbd5daDsO//2KkxUZ4XjlPlAvjv3LLSkp0bJlyzR8+PDvCxMebreXLFlS42tOO+00+5qKMLNp0ybNmTNH55133jE/p7i4WHl5edVuAIDq/vnhRr2/OkdREeH6x1Unq3lcFKcIQcuxmZh2796tsrIyJSdX76hmtteuLe/IdiRTY2Ned8YZZ9hhiocOHdKNN9543GapSZMm6f777/d5+QHALRat36VH/rfOPv7ThSepT3pTp4sEnJCgqnNcuHChHnzwQf3jH/+wfXRee+01vfPOO3rggQeO+ZoJEyZo//79lbft27f7tcwAEMhy8op060tf2g7EPzslXT8b2NbpIgHBW3OTlJSkiIgI5eTkVNtvtlNSUmp8zT333KOrr75a1113nd3u1auXCgsLdcMNN+gPf/iDbdY6UnR0tL0BAKor83jtyKjvDpTqpNQE3Tf6JE4RXMGxmpuoqCj1799f8+fPr9zn8Xjs9uDBg2t8zYEDB44KMCYgGaaZCgBQe9M/3Kglm/YoNipCf7+in2Ialf89BYKdo6ufmWHgY8eO1YABAzRw4EA7h42piTGjp4wxY8YoLS3N9psxRo0aZUdY9evXz86Js2HDBlubY/ZXhBwAwA9bvu07TXl/vX18/+iT1LFlPKcNruFouLn88su1a9cu3XvvvcrOzlbfvn01d+7cyk7G27Ztq1ZT88c//tHOuWDuMzMz1bJlSxts/vKXvzj4UwBA8M1n85sXv7TNUqP6pOqS/m2cLhLgU2HeEGvPMUPBExMTbefihATWSgEQesbPXqHXlmeqTbPGmnPrECXENHK6SIBPv7+DarQUAODEmLlsTLAxC2I+9rO+BBu4EuEGAELEd4UlmvDaSvv4+iEd1b9dc6eLBDQIwg0AhIh7/rtKuwuK1aVVvG4/p6vTxQEaDOEGAELAO19n6e2vsxQRHqZHL+vDsG+4GuEGAFxub2GJrbUxbhrWSb3bsLwC3I1wAwAu95d31tiA0y25iX59VheniwM0OMINALjY4g279Z/lOxQWJk36aS9FRfJnH+7Hv3IAcKmi0jL94Y3y5qifD2qnk9s2c7pIgF8QbgDApf6xYIM27y5UqybR+t3Ibk4XB/Abwg0AuNCG3Hz988ONlWtHMQsxQgnhBgBcxqyqc+9/v1FpmVfDe7TSyJ4pThcJ8CvCDQC4zNxV2Vq8cY+iI8M1cdRJdsFhIJQQbgDARQ6WlOnP76yxj385tJPSm8c6XSTA7wg3AOAiTy7aqMx9B5WaGKNfDe3kdHEARxBuAMAldnx3QP9cWN6J+A/nZ6hxVITTRQIcQbgBAJeYNGetig95dGrH5jqvF52IEboINwDgAl9s2at3VmYpPEx0IkbII9wAgAuGfj84p7wT8eWntFWP1glOFwlwFOEGAJxSUijdl1h+M4/r6b1vsrV82z41bhSh24ezMCZAuAGAIFZa5tHDc9fZx9cN6aBWCTFOFwlwHOEGAILYS59v16bdhWoeF6UbftTR6eIAAYFwAwBBqqD4kB6bt94+/s1ZndUkppHTRQICAuEGAILUMx9t1u6CErVrEasrB7VzujhAwCDcAIBTPGXfP966uPr2D9h/oFQzPtpkH9/x426KiuTPOVCB3wYAcMLqN6VpA7/ffv4SaWrP8v218MzHm5RffEhdk+N1Qa/WDVdOIAgRbgDA30yAmT1Gys+qvj8vq3z/DwSc7wpL9OwnW+zj24d3VbiZuQ9AJcINAPiTaXqae5eZeq+GJw/vm3v3cZuonv5ok+1MbCbrG3ESyywARyLcAIA/mb41eTuPc4BXysssP64GewtLNHNxRa1NF2ptgBoQbgDAnwpyTui4Jxdt1IGSMvVMS9A5Gcm+LRvgEoQbAPCn+OR6H7enoFizFm+t7GsTFkZfG6AmhBsA8Kd2p0kJqZKOFUzCpIS08uOOYJqjDpaWqVdaos7q3qrBiwoEK8INAPj1r26ENPKhwxtHBpzD2yMnlx9XhelA/NzhvjY3DetErQ1wHIQbAPC3jNHSZbOkJkeMdDI1Oma/ef4IL3y2VXlFh9SxZRwjpIAfEPlDBwAAGoAJMB2HSZPTy7evelXqdNZRNTZG8aEyzfhos31849BOjJACfgA1NwDglKpBxvSxqSHYGK8tz1RufrFaJ8boor5p/isfEGo1N/v27dPSpUuVm5srj8dT7bkxY8b4omwAEPLKPF49+eFGex6uG9KRNaSAhgo3b731lq666ioVFBQoISGhWsc285hwAwC+MWdllrbsOaCmsY10xcDDTVgAfB9u7rjjDl177bV68MEHFRsbW5+3AABExUn37T/mefB6vXapBWPs4PaKjaKbJNBgfW4yMzP1m9/8hmADAA1o2dbv9PWO/bYpaszgdpxroCHDzYgRI/TFF1/U56UAgFp69pPyEVI/6ZemFvHRnDeglupVx3n++efrd7/7nVavXq1evXqpUaNG1Z4fPfroORoAALW3fe8BzV2VbR+PO70Dpw5o6HBz/fXX2/s//elPRz1nOhSXlZXV520BAIfNWrJFHq90RuckdUtpwnkBGjrcHDn0GwDgO2aphZc+324fX3tGe04tUEdM4gcAAeY/y3Yo3yy1kBSnYV1ZIBPwW7j58MMPNWrUKHXu3NneTD+bjz76qL5vBwCwNeNe/etwR+Jxp7dnqQXAX+Hm3//+t4YPH26Hgpsh4ebWuHFjnX322XrhhRfq85YAAPMfx2932Un7EmIi9ZOT23BOgHoI85pZouqoR48euuGGG3T77bdX2z9lyhQ9/fTTWrNmjQJVXl6eEhMTtX//fju7MgAEkuue+1zz1uTq2tM76N5RGU4XBwjK7+961dxs2rTJNkkdyTRNbd5cXp0KAKibzH0H9cHaXPv4qlPbcvqAeqpXuElPT9f8+fOP2j9v3jz7HACg7l5aus0O/x7csYU6tYznFAL+XlvK9LNZsWKFTjvtNLvvk08+0cyZM/XYY4/VtywAELJKyzyVw79/fipLLQB+Dze/+tWvlJKSokcffVSzZ8+u7Ifz8ssv68ILLzyhAgFAKHp/dY525RcrKT5a52QkO10cIKjVe4nZiy++2N4AACfu+c+22vufnZJuF8oEUH/8BgGAwzbtKtAnG/YoLEz62UD6LQJ+CzfNmzfX7t277eNmzZrZ7WPd6mLatGlq3769YmJiNGjQIC1duvS4x+/bt08333yzWrdurejoaHXt2lVz5syp02cCQCB54bNt9v6sbq3Uplms08UBQqdZ6m9/+5uaNGlS+dgskHmiTB+d8ePHa/r06TbYTJ06VSNGjNC6devUqtXRU46XlJTonHPOsc+9+uqrSktL09atW9W0adMTLgsAOKHkkEevfZlpH185iOHfgGOT+PmKCTSnnHKKnnjiicoFOc1Q8l//+te6++67jzrehKC//vWvWrt2rRo1alSvz2QSPwCBZO6qbN3472Vq1SRai+8+S5ER9BYAHJnELyIiQrm55RNNVbVnzx77XG2YWphly5bZZRwqCxMebreXLFlS42vefPNNDR482DZLJScnq2fPnnrwwQdVVlZ2zM8pLi62J6TqDQACxStflA//vvjkNIIN4CP1CjfHquwxQSIqKqpW72H675hQYkJKVWY7Ozv7mDMjm+Yo8zrTz+aee+6xw9H//Oc/H/NzJk2aZJNexY1JBgEEitz8Ii1cv8s+vrQ/HYkBR4aCP/744/be9LeZMWOG4uO/n0HTBI5Fixape/fuaiim2cr0t3nqqadsDVH//v2VmZlpm6omTpxY42smTJhg+/VUMDU3BBwAgeD15Zkq83h1ctum6tyKGYkBR8KN6UhcUXNj+r9UbYIyNTZm1JPZXxtJSUn29Tk5OdX2m20zQWBNzAgp09em6ueayQNNTY9p5qqp1siMqDI3AAgk5u/oK8t22MeXDqDWBnAs3FQsinnmmWfqtddes0PC68sEEVPzYtaouuiiiyprZsz2LbfcUuNrTj/9dL3wwgv2ONM/x1i/fr0NPbVtDgOAQLBi+z5tyC1QTKNwXdC7tdPFAVylXn1uFixYcELBpoJpLnr66af13HPPac2aNXZZh8LCQo0bN84+P2bMGNusVME8v3fvXt1666021Lzzzju2Q7HpYAwAwaSi1ua8nq3VJKZ+oz8B+Hj5hR07dtjRS9u2bbNNQlVNmTKlVu9x+eWXa9euXbr33ntt01Lfvn01d+7cyk7G5r0ramgM01fmvffe0+23367evXvbeW5M0Lnrrrvq+2MAgN8VlZbpra922seXDGjDFQACYZ4b03Q0evRodezY0c45Y4Zkb9myxbYhn3zyyfrggw8UqJjnBoDTTLD59Ytfqk2zxlr0uzMVHn7ik6ICbpfX0PPcmKai3/72t1q5cqVdNuE///mPtm/frqFDh+rSSy+tb7kBICT8d0V5rc1FfdMINkADqFe4Mf1jTH8YIzIyUgcPHrTDwv/0pz/poYce8nUZAcA19h0o0YfryydBvahfqtPFAVypXuEmLi6usp+NGam0cePGyucqFtcEABxtzspslZZ5ldE6QZ1bla/XByAAOhSfeuqp+vjjj+0cM+edd57uuOMO20Rlhoeb5wAANXtjRfkimRf2pdYGCKhwY0ZDFRQU2Mf333+/fWxW+O7SpUutR0oBQKjZue+glm7eq7AwaTThBgiscGNGSVVtoqrtrMQAEMrePDz8e2D75mqd2Njp4gCuVa8+N2ZklJnnpsLSpUt122232TWfAADHHyV1Yd80ThEQaOHmyiuvtLMUG2byveHDh9uA84c//MGOmAIAVLc+J19rsvLUKCJM5/Wqef08AA6Gm1WrVmngwIH28ezZs9WrVy8tXrxYzz//vGbOnOmjogGAe/z3cEfioV1bqWksa+EBARduSktLK1fanjdvnp2t2OjevbuysrJ8W0IACHJm9va3vir/20hHYiBAw81JJ51kOxF/9NFHev/99zVy5Ei7f+fOnWrRooWvywgAQe2bnXnatveAXQH87O6tnC4O4Hr1CjdmFuInn3xSw4YN0xVXXKE+ffrY/WYhzYrmKgBAuXdXldfaDOvaSnHR9V6vGEAt1eu3zIQaMxOxWcSqWbNmlftvuOEGxcbG1uctAcC1TVJmVmLjXDoSA35R7/9CREREVAs2Rvv27X1RJgDwr5JC6cHDMwb/fqcUFeezt16Xk6/NuwsVFRmus3sk++x9Afgg3Jx88smaP3++DTT9+vVTmJli8xiWL19e27cFAFeb83V5k9TQri0VT5MUEFjh5sILL6wcIWUeHy/cAADKzVlV3iTF3DZAAIabiRMnVj6+7777Gqo8AOAa3+bka0NugZ24jyYpIMBHS5m1pfbs2XPU/n379lVbdwoAQllFR+IhXVoqIaaR08UBQka9ws2WLVtUVlZ21P7i4uJqa04BQCirGAJ+bk+WWwACdrSUmcemwnvvvafExMTKbRN2TIfjDh06+LaEABCENu4q0NrsfEWGh+mcDEZJAQEbbi666CJ7bzoTjx07ttpzjRo1skPBH330Ud+WEACC0Purc+z94E4tWEsKCORw4/F47L2pnfn888+VlJTUUOUCAFeEmx9TawMExyR+mzdvrnxcVFSkmJgYX5YJAILa7oJiLd/2nX08nHADBEeHYlOD88ADDygtLU3x8fHatGmT3X/PPffomWee8XUZASCofLA2V16v1DMtQa0TGztdHCDk1Cvc/PnPf9bMmTP18MMPKyoqqnJ/z549NWPGDF+WDwCCtklqOMstAMETbmbNmqWnnnpKV111lV1jqoJZHXzt2rW+LB8ABJWi0jJ99O0u+5hRUkAQhZvMzEx17ty5xuaq0tJSX5QLAILSJxt2q6jUo9TEGGW0TnC6OEBIqle4ycjI0EcffXTU/ldffdUuqgkACvUmqYxk1uADgmm01L333mvnuTE1OKa25rXXXtO6detsc9Xbb7/t+1ICQBDweLyatybXPqa/DRBkNTdmVfC33npL8+bNU1xcnA07a9assfvOOecc35cSAILAVzv22WHg8dGROrVjC6eLA4SsetXcGEOGDNH777/v29IAgAuapIZ2a6moyHr93xGAk+HGKCkpUW5ubuXMxRXatm17ouUCgKAz/3CT1DkMAQeCL9x8++23uvbaa7V48eJq+71er+1AV9OK4QDgZjv3HdS6nHyFh0lDu7Z0ujhASKtXuLnmmmsUGRlpOw+3bt2aEQEAQt6H68vntumT3lTN4r6f3BRAkISbFStWaNmyZerevbvvSwQAQWjhuvImqWFdWzldFCDk1Xuem927d4f8yQMAo+SQR59s2GMfD+tGkxQQlOHmoYce0p133qmFCxdqz549ysvLq3YDgFBiVgAvKD6kFnFR6pWW6HRxgJBXr2ap4cOH2/uzzz672n46FAMIRQvXlfe3+VHXlgo3PYoBBF+4WbBgge9LAgDB3t+GJikgeMPN0KFDfV8SAAhC2fuLtDY7X2Fh0pAu9LcBgjbcfP311zXuN3PcxMTE2En8oqOjT7RsABDwPlxfXmvTu01TNWcIOBC84aZv377HndumUaNGuvzyy/Xkk0/asAMAbu9vM4yJ+4DgHi31+uuvq0uXLnrqqafsnDfmZh5369ZNL7zwgp555hl98MEH+uMf/+j7EgNAgCgt8+jjb8unxaC/DRDkNTd/+ctf9Nhjj2nEiBGV+3r16qU2bdronnvu0dKlS+1q4XfccYceeeQRX5YXAALGl9v2Kb/4kJrFNrLNUgCCuOZm5cqVateu3VH7zT7zXEXTVVZW1omXEAAC1MffljdJndGlpSIYAg4Ed7gxyy5MnjzZrgpeobS01O6rWJIhMzNTycnJvispAASYjzeUN0kN6ZzkdFEAnGiz1LRp0zR69GjbDNW7d2+7z9TYmNXAzWKaxqZNm3TTTTfV5+0BIODlFZXqqx377ePTuxBugKAPN6eddpo2b96s559/XuvXr7f7Lr30Ul155ZVq0qSJ3b766qt9W1IACCCfbdqrMo9XHZLilNa0sdPFAXCi4cYwIebGG2+s78sBIKh9crhJ6rROLZwuCgBfhRtj9erV2rZtW7W+N4ZpsgKAUAg3Z9DfBnBHuDH9aS6++GLbz8ZM5mcWzDQqJvYzfW8AwK1y8or0bW6BXXJhMDU3gDtGS916663q0KGDcnNzFRsbq2+++UaLFi3SgAEDtHDhQt+XEgACsNamV1qimsZGOV0cAL6ouVmyZImdgTgpKUnh4eH2dsYZZ2jSpEn6zW9+oy+//LI+bwsAQTUE/HSapICAVK+aG9PsVDEqygScnTt3Vk7it27dunoNLW/fvr1dh2rQoEF2huPaeOmll2xT2EUXXVTnzwSA+jDN8Is37LGPT+/EEHDANeGmZ8+e+uqrr+xjE0YefvhhffLJJ/rTn/6kjh071um9Xn75ZY0fP14TJ07U8uXL1adPH7usg2nyOp4tW7bot7/9rYYMGVKfHwEA6mXjrkJl5xUpKjJcA9o34ywCbgk3ZkFMj8djH99///12zhsTMubMmWPXnKqLKVOm6Prrr9e4ceOUkZGh6dOn2348zz777HFrjq666ir72XUNUwDgi/42p7RvpphGEZxMwC19bqoumGlWB1+7dq327t2rZs2aVY6Yqg0zhHzZsmWaMGFC5T7Tf2f48OG2X8+xmBqiVq1a6Re/+IU++uij435GcXGxvVXIy8urdfkA4Ej0twFcFm6uvfbaWh13vFqXqnbv3m1rYY5cg8psm8BUk48//ljPPPOMVqxYUavPMJ2cTQ0PAJwoMyPxZ5vK+9ucRn8bwB3hZubMmbbTcL9+/SrntvGn/Px8u6zD008/bTsy14apFTJ9eqrW3KSnpzdgKQG41ZqsPOUVHVJ8dKR6piY4XRwAvgg3v/rVr/Tiiy/aPjamj8zPf/5zNW/eXPVlAkpERIRycnKq7TfbKSkpRx2/ceNG25F41KhRlfsq+v5ERkbakVqdOnWq9pro6Gh7A4AT9dnmvfbedCSOjKhXl0UAfhBe1yHbWVlZuvPOO/XWW2/ZGpDLLrtM7733Xr1qcqKiotS/f3/Nnz+/Wlgx24MHDz7q+O7du9tZkU2TVMXNLPVw5pln2sfUyABoSBVNUoM6sJ4U4KoOxaYW5IorrrC3rVu32qaqm266SYcOHbIzFcfHx9fp/UyT0dixY+3sxgMHDtTUqVNVWFhoa4aMMWPGKC0tzfadMfPgmGHoVTVt2tTeH7kfAHzJ4/Fq6ZbymptBHetfYw0gwBfONCObKtaWqu96Updffrl27dqle++9V9nZ2erbt6/mzp1b2cnYLMxpPgcAnLQ+N1/7DpQqNirCLrsAwEXhxgyrfu211+yIKDNy6YILLtATTzyhkSNH1juE3HLLLfZWkx9aq8rUHAFAQ/t0Y3mTVP92zdSI/jaAe8KNaX4ySx6Yvi1mWLjpXFzbUUsA4IbOxIM60CQFuCrcmNmD27Zta2cF/vDDD+2tJqZmBwDcwjS9L60INx3pTAy4KtyYzr11mYEYANxgQ26B9hSWKDoyXL3b0N8GcN0kfgDgOp4qAyK2LpY6nSWFf79u1KeHa21Mf5voSNaTAgIdw5AAhLbVb0rTBn6//fwl0tSe5fsPY34bILgQbgCELhNgZo+R8rOq78/LKt+/+k3b36ayMzHz2wBBgXADIHSboubeZboL1/Dk4X1z79bm3Dztyi9WVGS4+qaXTxoKILARbgCEJtO3Jm/ncQ7wSnmZ2rx8nt0ywSamEf1tgGBAuAEQmgqqL9h7LFk7ttj7U9o3a+ACAfAVwg2A0BRfvsTLD/liT5S9H9COyfuAYEG4ARCa2p0mJaRKOtbcXWEqa5KqN/e1t1snt6XmBggWhBsAocnMYzPyocMbRwac8u2ve06QR+HqmhyvxNhGfi8igPoh3AAIXRmjpctmSU1Squ83NTqXzdLcslMqJ+8D4OJVwQHAdQGn4zBpcnr59lWvVs5Q/MWHi+2u/vS3AYIKNTcAUGWpBdsXJzxCRaVlWrljv901gJobIKgQbgCgBt/s3K+SMo+S4qPUrkUs5wgIIoQbAKjBF1u+qxwlFRZ2rBFVAAIR4QYAavDF1vJwM4DJ+4CgQ7gBgCOYxTKXHw43dCYGgg/hBgCOsGXPAe0pLLGLZfZMS+D8AEGGcAMAR/hiy1573zstUdGRLJYJBBvCDQAcYfm2w01S9LcBghLhBgCOMVKKxTKB4ES4AYAq8opKtWFXgX3cr21Tzg0QhAg3AFDFNzvz5PVK6c0bKyk+mnMDBCHCDQBU8fXhJRf6tKHWBghWhBsAqOKr7fvsfd90wg0QrAg3AFBDzQ39bYDgRbgBgCr2FBYrMjxMJ6Umcl6AIEW4AYAjdG/dRDGNmLwPCFaEGwA4Ap2JgeBGuAGAI9CZGAhuhBsAOALhBghuhBsAqKJJdKQ6tYznnABBjHADAFX0TEtUeHgY5wQIYoQbAKiiVxpDwIFgR7gBgCp6tyHcAMGOcAMg5BUUH6o8B70IN0DQI9wACHnfZOZVnoNWTWJC/nwAwY5wAyDkrdxZvp4UAHcg3AAIed8QbgBXIdwACHmrd37fLAUg+BFuAIS0/QdLtW3vAaeLAcCHCDcAQto3mfS3AdyGcAMgpK0k3ACuQ7gBENJW0d8GcB3CDYCQtoqaG8B1CDcAQlZeUak27y50uhgAfIxwA0ChPjNxamJjp4sCwIcINwAU6k1SGalNnC4KAB8i3AAIWasOz0yc0ZqVwAE3IdwAUKgPA89IS3C6KAB8iHADICQVFB+q7Ex8UmvCDeAmARFupk2bpvbt2ysmJkaDBg3S0qVLj3ns008/rSFDhqhZs2b2Nnz48OMeDwDHmpnY65VaJ8YoKT6akwS4iOPh5uWXX9b48eM1ceJELV++XH369NGIESOUm5tb4/ELFy7UFVdcoQULFmjJkiVKT0/Xj3/8Y2VmZvq97ACCv0mqZxr9bQC3cTzcTJkyRddff73GjRunjIwMTZ8+XbGxsXr22WdrPP7555/XTTfdpL59+6p79+6aMWOGPB6P5s+f7/eyAwj+kVK9CDeA6zgabkpKSrRs2TLbtFRZoPBwu21qZWrjwIEDKi0tVfPmzWt8vri4WHl5edVuAFCx7EJPOhMDruNouNm9e7fKysqUnJxcbb/Zzs7OrtV73HXXXUpNTa0WkKqaNGmSEhMTK2+mGQtAaDtQckgbdxXYxz1TaZYC3MbxZqkTMXnyZL300kt6/fXXbWfkmkyYMEH79++vvG3fvt3v5QQQWNZl59vOxKYjcauEmv92AAhekU5+eFJSkiIiIpSTk1Ntv9lOSUk57msfeeQRG27mzZun3r17H/O46OhoewOACmuy8u19j9bMTAy4kaM1N1FRUerfv3+1zsAVnYMHDx58zNc9/PDDeuCBBzR37lwNGDDAT6UF4BarsypmJmZ+G8CNHK25Mcww8LFjx9qQMnDgQE2dOlWFhYV29JQxZswYpaWl2b4zxkMPPaR7771XL7zwgp0bp6JvTnx8vL0BQG1rbjJSCTeAGzkebi6//HLt2rXLBhYTVMwQb1MjU9HJeNu2bXYEVYV//vOfdpTVJZdcUu19zDw59913n9/LDyC4eDxerc0qHynVg5obwJUcDzfGLbfcYm/HmrSvqi1btvipVADcaNveAyosKVNUZLg6JsU5XRwADSCoR0sBQF2tOVxr0y25iSIj+BMIuBG/2QBCMtwwUgpwr4BolgIAf1ldU3+bqDjpvvIRVACCHzU3AEJzpBSdiQHXItwACBn7D5Qqc99B+7g74QZwLcINgJBrkkpr2liJjRs5XRwADYRwAyDkOhMzeR/gboQbACE4UoqZiQE3I9wACLlmqQwWzARcjXADICSUlnn0bU6BfZzROtHp4gBoQIQbACFh065ClZR5FB8dqTbNGjtdHAANiHADICSsztpfOTNxeHiY08UB0IAINwBCwtrs8sn7uqfQmRhwO8INgJCw/nC46ZbSxOmiAGhghBsAIWH94c7EhBvA/Qg3AFwvv+j7ZRe6tqLmBnA7wg2AkKm1SUmIUWIsyy4Abke4AeB663PK+9t0pb8NEBIINwBcb11FZ+LkeKeLAsAPCDcAQibcdE2mvw0QCgg3AEKmWYqRUkBoINwAcLXdBcXaU1iisDCpcyuapYBQQLgBEBKT97VtHqvYqEiniwPADwg3AFxtXcVIKfrbACGDcAMgNPrbEG6AkEG4ARAaI6WY4wYIGYQbAK7l9Xq/X1OKmhsgZBBuALjWzv1FKig+pMjwMHVIinO6OAD8hHADwPUjpTq2jFNUJH/ugFDBbzsA12KkFBCaCDcAXF9zQ38bILQQbgC4v+aGkVJASCHcAHClMo9XG3IZKQWEIsINAFfavveAig95bEfi9OaxThcHgB8RbgC4UkWtTcekOEWEhzldHAB+RLgB4Eobd5WHG1YCB0IP4QaAq2tuCDdA6CHcAHClDdTcACGLcAPAlWtKUXMDhC7CDQDX2ZVfrPyiQzL9iFlTCgg9hBsArlNRa9O2eayiIyOcLg4APyPcAHAd+tsAoY1wA8C1NTedWsU7XRQADiDcAHCdys7ELQk3QCgi3ABwHUZKAaGNcAPAVfKKSpWbX2wf0ywFhCbCDQBX1tokJ0QrIaaR08UB4ADCDQBXoUkKQCSnwEc8ZdLWxVJBjhSfLLU7TQpnfg3A3zbSmRgIeYQbX1j9pjT3Lilv5/f7ElKlkQ9JGaMb9h8ZoQqohpobAIQbXwSb2WPMajbV9+dlle+/bFbDBZxQDFWEOdRyAj86EwOhi3Bzol+0JlwcGWwssy9Mmnu31P1833/xh2KoIsy5P0Se4OcWlZZp+94D9jFz3AChKyA6FE+bNk3t27dXTEyMBg0apKVLlx73+FdeeUXdu3e3x/fq1Utz5syRI8wf4apf8EfxSnmZ5cf5NVSpPFSZ43ytIlQd+XNXhCrzfENw6nMrPntqT+m5C6T//KL83mw35GfyufU6z5t3F8rjlZrERKplk+gGvTwAApfj4ebll1/W+PHjNXHiRC1fvlx9+vTRiBEjlJubW+Pxixcv1hVXXKFf/OIX+vLLL3XRRRfZ26pVq/xedvu/S18eV1uhFqoIc+4PkT763Kr9bcLCwhqipACCgOPhZsqUKbr++us1btw4ZWRkaPr06YqNjdWzzz5b4/GPPfaYRo4cqd/97nfq0aOHHnjgAZ188sl64okn/F52W23uy+NqK9RCFWGu+jl2W4j04ed+m5Nv72mSAkKbo31uSkpKtGzZMk2YMKFyX3h4uIYPH64lS5bU+Bqz39T0VGVqet54440ajy8uLra3Cnl5eT4rv+0PYPqamP9d1vCH2SPpu4iWGv9BpLxhx29qq4vuRd/p97U47sGPvtPaz333uacWfqKbanHcP97+RJ/GRQf953Yv+kq/r0WYe3D6s1ob04fPDYDzvHrnfnvfr20zn10PAMHH0XCze/dulZWVKTm5es2G2V67dm2Nr8nOzq7xeLO/JpMmTdL999+vBmE6OppOtLPHyKswhVUJOKbd3/j9wav04bd7ffqxHytF10Q3V4r2KryGmnfz2dlqoRnbUuTRLp99bkl4hG6K+uHjFmVF6FNP8H9u0/AtUi0+N3vHFi3ypPK5AXKeoyLDNbxHK59dDwDBx/WjpUytUNWaHlNzk56e7rsPMKODLpslz5y7FFHw/f8+ixqn6Oued2tE6jkaId/btvNetf7i9sNjsr4PVTZkhUnbBtyjR1L7+fZDvT114P2n1bgot9pnVv3sgzHJumzUZbosLCLoPzdpd7FUixa2S4b117Ak39Xc8Lkndp67JjdRq4QYn10PAMHH0XCTlJSkiIgI5eRU7xtitlNSUmp8jdlfl+Ojo6PtrUFljFaEGe5dZQhrbLvTdGpDDp09eZzUscVRQ7LD7JDsyTq1oYZkN37k8BB0U2VUNWiYmispdvRf9ZOMdu74XE9radWxmx1tWRJS9aPhF/p2mDSf65/zDMC1HO1QHBUVpf79+2v+/PmV+zwej90ePHhwja8x+6seb7z//vvHPN5vzB/dDkOkXpeU3/vjj7AJMLetksa+Lf30mfL721Y27Fwzh2uqlNC6+n4Tqhpybh0nPrei2dE6sv3v8PbIyb6/1nyuf84zANcK83q9Nf2X1K9DwceOHasnn3xSAwcO1NSpUzV79mzb58b0pRkzZozS0tJs35mKoeBDhw7V5MmTdf755+ull17Sgw8+aIeR9+zZ8wc/zzRLJSYmav/+/UpISPDDT+hSQTrJm+8mD0wr/8L1+6SFfC6A0JRXh+9vx8ONYYZx//Wvf7Wdgvv27avHH3/cTuZnDBs2zE7wN3PmzGqT+P3xj3/Uli1b1KVLFz388MM677zzavVZhBvUSyiFuVD8XAABL+jCjT8RbgAAcPf3t+OT+AEAAPgS4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALiKo6uCO6FiQmYz0yEAAAgOFd/btVlYIeTCTX5+vr1PT093uigAAKAe3+NmGYbjCbm1pTwej3bu3KkmTZooLCyscv8pp5yizz//vMbXHOu5I/ebVGlC0/bt2x1fcfx4P48/368ur6vNsfW5Tsd6rqZ9gXINQ/H6He95fge5hk4Ixd/DUwL4u9DEFRNsUlNTFR5+/F41IVdzY05ImzZtjtofERFxzItwrOeOtd/sczrcHO/n8ef71eV1tTm2PtfpWM8d73inr2EoXr/jPc/vINfQCaH4exgR4N+FP1RjU4EOxYfdfPPNdX7ueK9xmq/LVt/3q8vranNsfa7TsZ7j+vnmmtTl2B86pq6/a1xD31yXuhzLNaw/N/4dDdTfw5BrlgqU5dgRmLiGwY3rF/y4hsEvLwC+C6m58aHo6GhNnDjR3iM4cQ2DG9cv+HENg190AHwXUnMDAABchZobAADgKoQbAADgKoQbAADgKoQbAADgKoQbAADgKoQbB5gpqYcNG6aMjAz17t1br7zyihPFwAm6+OKL1axZM11yySWcyyDx9ttvq1u3burSpYtmzJjhdHFQR/zOBbftfvzuYyi4A7KyspSTk6O+ffsqOztb/fv31/r16xUXF+dEcVBPCxcutOucPPfcc3r11Vc5jwHu0KFD9o/qggUL7ARj5vdu8eLFatGihdNFQy3xOxfcsvz43UfNjQNat25tL66RkpKipKQk7d2714mi4ASY/4GYBVgRHJYuXaqTTjpJaWlpio+P17nnnqv//e9/ThcLdcDvXHBr7cfvPsJNDRYtWqRRo0bZlUfNyuFvvPHGUcdMmzZN7du3V0xMjAYNGmT/cNbHsmXLVFZWZldQRXBeQwTHNd25c6cNNhXM48zMTC5fkFw/uOsaLmvg7z7CTQ0KCwvVp08fe5Fq8vLLL2v8+PF2eunly5fbY0eMGKHc3NzKY0w67dmz51E38we2gkmsY8aM0VNPPdUQ1zak+esaIriuKZzD9Qt+hT76HfTLd59ZOBPHZk7R66+/Xm3fwIEDvTfffHPldllZmTc1NdU7adKkWp/KoqIi75AhQ7yzZs3i9AfpNTQWLFjg/elPf+qzsqLhruknn3ziveiiiyqfv/XWW73PP/88pzzIfif5nQvua1jkp+8+am7qqKSkxFanDR8+vHJfeHi43V6yZEltA6WuueYanXXWWbr66qvrWgQEwDVE8F3TgQMHatWqVbYpqqCgQO+++679XyWcx+9kaFxDrx+/+wg3dbR7927bTpicnFxtv9k2vb9r45NPPrHVd6a90jR9mNvKlSvrWhQ4eA0N80t76aWXas6cOWrTpg3BKMCvaWRkpB599FGdeeaZ9nfujjvuYKRUkP1O8jsX3NfwEz9+90U2yLviuM444wx5PB7OUpCbN2+e00VAHY0ePdreEJz4nQtuZ/jxu4+amzoyQ9ciIiLsWP2qzLYZ2obAxzV0H65pcOP6Bb+kAPtuJNzUUVRUlJ14aP78+ZX7TBI124MHD/b19UED4Bq6D9c0uHH9gl9UgH030ixVA9PZcMOGDZXbmzdv1ooVK9S8eXO1bdvWDnUbO3asBgwYYDspTp061Q6RGzdunD+vHY6Da+g+XNPgxvULfgXB9N3YoGOxgpQZamhOzZG3sWPHVh7z97//3du2bVtvVFSUHf726aefOlpmVMc1dB+uaXDj+gW/BUH03cjaUgAAwFXocwMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAPAb4YNG6bbbrstID+jffv2drp4AMGPcAMAAFyFcAMAAFyFcAPAEf/3f/9nVw9u0qSJUlJSdOWVVyo3N7fy+YULFyosLEzvvfee+vXrp8aNG+uss86yx7z77rvq0aOHEhIS7OsOHDhQ7b0PHTqkW265RYmJiUpKStI999xjFgmufN68x6hRo+x7dujQQc8///xR5ZsyZYp69eqluLg4paen66abbrKrIgMIfIQbAI4oLS3VAw88oK+++kpvvPGGtmzZomuuueao4+677z498cQTWrx4sbZv367LLrvM9o154YUX9M477+h///uf/v73v1d7zXPPPafIyEgtXbpUjz32mA0qM2bMqHzefI55rwULFujVV1/VP/7xj2rByggPD9fjjz+ub775xr7fBx98oDvvvLMBzwgAn3FkLXIAIWno0KHeW2+9tcbnPv/8c1O14s3Pz7fbCxYssNvz5s2rPGbSpEl238aNGyv3/fKXv/SOGDGi2mf06NHD6/F4Kvfddddddp+xbt06+x5Lly6tfH7NmjV239/+9rdjlv2VV17xtmjRot4/OwD/oeYGgCOWLVtmm4batm1rm6aGDh1q92/btq3acb179658nJycrNjYWHXs2LHaviNrXU499VTbpFVh8ODB+vbbb1VWVqY1a9bYWp3+/ftXPt+9e3c1bdq02nvMmzdPZ599ttLS0mz5rr76au3Zs+eoJjAAgYdwA8DvCgsLNWLECNtnxvR3+fzzz/X666/b50pKSqod26hRo8rHJrBU3a7Y5/F4fFo+00R2wQUX2GD1n//8xwaxadOm1Vg+AIEn0ukCAAg9a9eutbUgkydPtp11jS+++MJn7//ZZ59V2/7000/VpUsXRURE2Foa0+HYBJZTTjnFPr9u3Trt27ev8njznAlMjz76qO17Y8yePdtn5QPQsKi5AeB3pikqKirKdgTetGmT3nzzTdu52FdM09b48eNtaHnxxRft59x66632uW7dumnkyJH65S9/aUOQCTLXXXedHTlVoXPnzrbDc0X5zMiu6dOn+6x8ABoW4QaA37Vs2VIzZ87UK6+8ooyMDFuD88gjj/js/ceMGaODBw9q4MCBuvnmm22wueGGGyqf/9e//qXU1FTbz+cnP/mJfa5Vq1aVz/fp08eOsHrooYfUs2dP23Q2adIkn5UPQMMKM72KG/gzAAAA/IaaGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAAIDf5fyjgGpmbwiUsAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lambda_SE_list = np.logspace(-2, 2, 1024)\n",
    "m_SE_list = 1 - 1/lambda_SE_list\n",
    "m_SE_list[m_SE_list < 0] = 0\n",
    "plt.plot(lambda_SE_list, m_SE_list, label='SE')\n",
    "\n",
    "d = 1000\n",
    "samples = 4\n",
    "lambda_AMP_list = np.logspace(-2, 2, 16)\n",
    "m_AMP_list = np.zeros((len(lambda_AMP_list), samples))\n",
    "for i, lambd in tqdm(enumerate(lambda_AMP_list)):\n",
    "    for s in range(samples):\n",
    "        m_AMP_list[i, s], _ = AMP_Gaussian(d, lambd)\n",
    "plt.errorbar(lambda_AMP_list, m_AMP_list.mean(axis=1), m_AMP_list.std(axis=1), fmt='o', label=f'AMP (d={d})')\n",
    "\n",
    "\n",
    "plt.xlabel('lambda')\n",
    "plt.ylabel('Magnetisation')\n",
    "plt.xscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b05829f1",
   "metadata": {},
   "source": [
    "# Checking that PCA and AMP estimate two alligned vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d32d4e28",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cosine similarity PCA-AMP:  0.9998200332471496\n"
     ]
    }
   ],
   "source": [
    "d = 1000\n",
    "lambd = 2.0\n",
    "T = 500\n",
    "\n",
    "Y, x0 = generate_Y(d, lambd)\n",
    "\n",
    "b = 0\n",
    "x_est = np.random.normal(0, 1, d)\n",
    "x_prev = np.zeros(d)\n",
    "m_AMP = []\n",
    "for t in range(T):\n",
    "    u = np.sqrt(lambd/d) * Y @ x_est - b * x_prev\n",
    "    sigma = lambd * x_est @ x_est / d\n",
    "    x_prev = x_est\n",
    "    x_est = u / (sigma + 1)\n",
    "    b = 1 / (sigma + 1)\n",
    "\n",
    "    m_AMP.append(np.abs(x_est @ x0) / d)\n",
    "    if t > 2 and np.abs(m_AMP[-1] - m_AMP[-2]) < 1e-6:\n",
    "        break\n",
    "\n",
    "# diagonalize\n",
    "eigenvalues, eigenvectors = np.linalg.eigh(Y)\n",
    "\n",
    "print(\"Cosine similarity PCA-AMP: \", np.abs(eigenvectors[:,-1] @ x_est / (np.linalg.norm(eigenvectors[:,-1]) * np.linalg.norm(x_est))))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af3bd9e2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "global (3.14.0)",
   "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.14.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}