Waterfilling_examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f184e1d0-29a8-4719-b5ce-16d0a094b7eb",
   "metadata": {},
   "source": [
    "# Waterfilling examples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33405aae-4058-4fbf-95b2-aaac775be642",
   "metadata": {},
   "source": [
    "We consider two parallel channels, one with gain $h_1$ and another with gain $h_2$, both fixed numbers known to both the transmitter and receiver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "676309d0-96f1-4d27-beed-0ebffa2c851c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Channel 1 absolute gain is 1.0770329614269007\n",
      "Channel 2 absolute gain is 0.3605551275463989\n",
      "[503.37832812 496.62166477]\n",
      "50.34% in channel 1, 49.66% in channel 2\n",
      "10.55445134463954\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import cvxpy as cp\n",
    "\n",
    "h1 = 1 + 0.4j\n",
    "h2 = -0.3 + 0.2j\n",
    "P_total = 1000\n",
    "P = cp.Variable(2)\n",
    "objective = cp.Maximize(cp.log(1 + np.abs(h1)**2 * P[0]) + \\\n",
    "                        cp.log(1 + np.abs(h2)**2 * P[1]))\n",
    "constraints = [P >= 0, P[0] + P[1] <= P_total]\n",
    "prob = cp.Problem(objective, constraints)\n",
    "result = prob.solve()\n",
    "print(f\"Channel 1 absolute gain is {np.abs(h1)}\")\n",
    "print(f\"Channel 2 absolute gain is {np.abs(h2)}\")\n",
    "print(P.value)\n",
    "print(f\"{P.value[0] / P_total * 100:.2f}% in channel 1, {P.value[1] / P_total * 100:.2f}% in channel 2\")\n",
    "print(objective.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db5c88d7-78ca-4850-89eb-c3434a0c958e",
   "metadata": {},
   "source": [
    "## OFDM case\n",
    "We will consider an FIR channel and perform waterfilling across the OFDM subcarriers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "832e799d-8837-4a93-9cd3-3d9b18df2583",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUdJJREFUeJzt3Xt8VNW5P/7PnnsScgFCEgKBAKJIQRDENNaqp6YGSy+0llKOp3ioR79a0qON5VisBa09P6xVqy0cOfX8rPZrEUpbKVpOWgiCrQQityKC3CQEIZMQQu7J3Pb6/pHsyYVJyCQzs/be83m/XvNqmdmz55ktzDyz1rOepQghBIiIiIh0zCI7ACIiIqIrYcJCREREuseEhYiIiHSPCQsRERHpHhMWIiIi0j0mLERERKR7TFiIiIhI95iwEBERke7ZZAcQCaqq4vz580hOToaiKLLDISIiogEQQqCpqQnZ2dmwWPofQzFFwnL+/Hnk5OTIDoOIiIgG4ezZsxg7dmy/x5giYUlOTgbQ8YZTUlIkR0NEREQD0djYiJycnOD3eH9MkbBo00ApKSlMWIiIiAxmIOUcLLolIiIi3WPCQkRERLrHhIWIiIh0jwkLERER6R4TFiIiItI9JixERESke0xYiIiISPeYsBAREZHumaJxXLQEVIHy03WoaWpHRrILN04YAauFexURERHFGhOWPpQcrsKTbx1BVUN78L5Ulw0F12YgKy0BQgDDEx0YkeRAfasXaYkd/ztimBNZKUxuiIhI/7Qf5u6GNtS1dHyX1bV4UN/mu+x7Tvb3GxOWEEoOV+HB1/dD9Lq/od2PPxw4P6BzJLus+PqssbjjU6OZvBARkW4EVIHdpy7i9T0V+NuJWjR7AmE9f3SqCyu/NBVzp42OUoShKUKI3t/LhtPY2IjU1FQ0NDQMeS+hgCpw80+39xhZGapUlw2fn5qJz0wexdEXIiKKqe6jKO+drMWWw260esNLUnpTALz0L7OGnLSE8/3NEZZeyk/XBZMVIQTqtr4Ei90FxZEAW8oo2EfmwJE5EYrVPuBzNrT78fv95/D7/ecAcPSFiIiip3eCsvVoDRrafIM6l9reDE/1KfgvfgLV24rUTy8IPvbkW0fw+alZMfsOY8LSS01Tt5GVgA/NB7Zcdoxid8E1/jokXXsrEiZ/Gha7M6zXaGoP4Ne7zuDXu85w9IWIiIZMS1L++mEVfr//HJra/YM+l7+5Dq1HdqDl2HvwVp0AhAoAUGxOpOTdBUWxQACoamhH+ek65E8aGaF30T8mLL1kJLt6/Dn15rshvG1QPS3w17vhrTkNta0RbSfL0XayHBbXMCTP+iKSZ38J1sTUsF+v9+hLWoIdSz6Ti6LPTWbiQkRE/QqoAqu3n8Sv3zuN+kGOomi8FyrQuPv3aDn6bjBJAQBbWhbs6eNgHzEWwu+FYu/6nuzxIz/KWMPSy5VqWIRQ4aupQOux99D84TsINNYAABRHItJuvhvJs+ZBsQ49D0xyWvGzu67DF67LHvK5iIjIXLRE5b/fPTXkepRASz3q3/0Nmg9tBTqXmzizpyBp2ueQMGkObCmj+nzuG/d9ekgjLOF8fzNhCaHkcBUeeH3/FY8TagCtJ3ajYdcG+Go+BgA4Midh5LzvwTEqd8hxAMCXrsvCC9+cxdEWIiICAGw5VIX/+MMhNHsGP+2jaTn6Lur++hLU9iYAQOI1n0Fq/jfgyJzU7/MUAFmpLvz90c8N6fuJCUsElByuwqN/OISGtiv/hRBqAM0fbEP9jlc7/qNbbRhx+/0YNvNOKMrQE41EhxX/55aJnCYiIopjAVXgofUH8PahqiGfS/W1o+6v/4WWw9sBAPaMCRjx+QfhGjt1QM+XsUqICUs/mtp9mP7EXwd8vL+5DnUlv0TbqfcBAEnTP4+Rhd8Ja0VRf5i4EBHFH236Z+3Ok2jzqVd+whX4G2tQ84efdMwMKBak5i9E6k0LB1zOEMk+LExYIqTV68fUFX8J6zlCCDSW/wH1O38DCBWu3Osx6quPweJIiFhcrG8hIjK/SNapaLwXzqDmdysQaL4IS2IqRn3lUbjGXTfg57+6ZA4+O3lUxH40h/P9zc0PI0xRFKTmfR0ZX18Jxe5Ee8UBVK//IdT25oi9RosngO+sO4DvrtuHgGr4fJOIiHrZcqgKM578K36+7XjEkhVP1QlUr3sUgeaLsI8ch9H3vBBWsgJAausNJixRkjBxNjK/+f/B4kqGt+o4qn+3AqqnJaKv8dYhN6Y/8Re8uO04ExciIhMIqAJF6/bjO+v2R6SoVuOpOoGaDY9DbW+GI/saZN79035X/+gRE5YocmZfg8xF/xlMWmo2PgnV54noa7R6A/j5thNMXIiIDCygCry47QSmrSyJSFFtd97aStT87kdQPS1wjpmKzG88BWtCckRfIxYGlbCsWbMGubm5cLlcyMvLQ3l5eZ/Hfvjhh7jrrruQm5sLRVHwwgsvXHbMqlWrMGfOHCQnJyMjIwPz58/HsWPHBhOa7jgyJnYkLc4keM4dQe3bz0KokRne605LXK578i/YcmhgGzQSEZFcWqIy/Ym/4OfbjkekqLY7f1Mtan63MjiykrHgCViciRF9jVgJO2HZsGEDiouLsXLlSuzfvx8zZsxAYWEhampqQh7f2tqKiRMn4umnn0ZWVlbIY3bu3ImlS5di9+7d2Lp1K3w+H+644w60tER2CkUWR8ZEjLrrR4DVhrbjZbi049dRey3WtxAR6V/vRCVSdSrdqb52XPjDUwg0XYBtxBhkfH2lYZMVYBCrhPLy8jBnzhysXr0aAKCqKnJycvDd734XP/jBD/p9bm5uLh5++GE8/PDD/R534cIFZGRkYOfOnbjllluuGJOeVgn1p+Xo31C7+acAgJHzijFs2ucidu5QuAyaiEh/Itn4rS9CCNS+9Sxaj+6EJSEFWYufhz0t9KBBOI78uBCJjsjt6hO1VUJerxf79u1DQUFB1wksFhQUFKCsrGxw0YbQ0NAAABgxYkTIxz0eDxobG3vcjCDp2s8iNX8hAOBiyS/hrT4V1dfTpolm/2QrSg5Hdk6UiIjCt2rLkYgX1IbStPdPaD26E7BYMWr+8ogkK7KFlbDU1tYiEAggMzOzx/2ZmZlwu90RCUhVVTz88MP4zGc+g2nTpoU8ZtWqVUhNTQ3ecnJyIvLasZD62buRMGkOEPDhwp+ehuppjfpr1rf68MDr+1nbQkQk0dsHz+O/3z0d9dfxnD8WLD0Y/rl/g2vc9Ki/ZizobpXQ0qVLcfjwYaxfv77PY5YvX46Ghobg7ezZszGMcGgUxYKR84phTR4F/6UqXPzrmpi99tJ1B/DCVq4kIiKKpYAq8POtx1G0/kDUX0ttb8aFzc8AagCJ19yM5FlfjPprxkpYCUt6ejqsViuqq6t73F9dXd1nQW04ioqK8Pbbb+Odd97B2LFj+zzO6XQiJSWlx81IrAnJSP/yMkCxoPXIzo6tvGNAAHihlFNERESxUnK4CrN/shUvlp6IyevVbftvBBqqYUvNxMg7vxuR/ez0IqyExeFwYPbs2SgtLQ3ep6oqSktLkZ+fP+gghBAoKirCm2++ie3bt2PChAmDPpdRuMZODdaz1P31Jfib62L22toUEZMWIqLo2XKoCg+8vh/1rb6YvF7r8V1o+fAdQLEg/Uvfh8WZFJPXjZWwp4SKi4vx8ssv47XXXsPRo0fx4IMPoqWlBUuWLAEALF68GMuXLw8e7/V6cfDgQRw8eBBerxfnzp3DwYMHcfLkyeAxS5cuxeuvv45169YhOTkZbrcbbrcbbW1tEXiL+pV600I4MidBbW9C3V9iNzWkeWTjP+D1R3bNPxERddSrLH1jf8xeL9DWiIud3yMpeV+Dc8y1MXvtWAk7YVm4cCGeffZZrFixAjNnzsTBgwdRUlISLMStrKxEVVXXL/fz58/j+uuvx/XXX4+qqio8++yzuP766/Fv//ZvwWNeeuklNDQ04LbbbsPo0aODtw0bNkTgLeqXYrVh5LxiwGJD28k9aD2xJ6av3+IJsEMuEVEEaa31i9YfQCy3Fq7f+RuorQ2wjxyHtM/cHbsXjiHu1tyPSPdh6culna+hcfdGWFMykH3vf8HicEX9NXtLS7Tj6a9Nj8h24URE8ajkcBUe/cMhNLRFd8lyb55zH8H9+jIAApn//DRcOaFX2EaCYfqwUHSk3rQQ1pQMBBpr0FAmZ1SJS5+JiAZPq1eJdbIi1ADqtr4EQCBp2u1RTVZkY8KiAxa7CyMK/g8AoLH8Tfhq5S3TLnrjALZEeOMtIiIz23LoPIpiWK/SXdOBLfBWn4LFmYThty2REkOsMGHRicTJeUi46kZA9ePi1pcga6ZOFcB31nEFERHRQGw5VIXvrDsAGWWA/uY61L/7fwEAabfeA2tSWuyDiCEmLDoy/Pb7odic8FQeQuux96TGwhVERET9i/VKoN7qd74G4W2FY/RkDJtRKC2OWGHCoiP2tCyk5N0FAKh/9zcQgdjOhXbX4glg1lNsMEdEFMqqLUdivhKoO++FCrQc3g4AGFHwABSLVU4gMcSERWdS5syHJTEV/kvn0fzBNqmxNHv8LMQlIuolVnsC9af+b68DEEi8+iY4s6+RGkusMGHRGYszMdgBt+G9dVB9HskRsRCXiEiz5dB5fDcGewL1x3PuKNpO7AYUC9Ju+ZbUWGKJCYsOJc+8s2OZc3Mdmva/LTucYCEuR1qIKJ5pBbYym5cJIXBp52sAgGHTC2AfmSMxmthiwqJDis2OtM92dCps3L0Ranuz5Ig6cKSFiOKVzKXL3bWf3g/P2cOA1Y7UzyySHU5MMWHRqaSpt8GePg5qezMa9vxBdjgAuOSZiOKTzKXL3QmhBkdXUmZ9EbaUUXIDijEmLDqlWKxIu+UeAEDTvrcQaGuSHFEXLnkmonghe+lyd20ndsNX8zEURwJS8hfIDifmmLDoWMJVN8KeMQHC146mA3+WHU4QlzwTUTyQvXS5OyEEGsp+BwBInv1lWBMit2+eUTBh0TFFUZCa93UAQNPezVB97ZIj6sIlz0RkZnpYutxde8UBeN0nodidSLnhy7LDkYIJi84lTrkZttRMqG2NaJHclyUUFuISkdnoYelyb9royrAZc2FNTJUcjRxMWHROsViRcuPXAAAN5W9CqAHJEfXEJc9EZCZ6WLrcm6fqRMfKIIsNKTd+VXY40jBhMYCk6QWwJKYi0FCN1o/+JjuckDjSQkRGp5ely7017t0EAEi69rOwJafLDUYiJiwGYLE7kTK7Y86yYffvpe3k3B+OtBCRkell6XJv/qZatH70dwBA8g1fkRyNXExYDGLYrHlQHAnwXahA+8f7ZIfTJ460EJHR6HVkBUBHt3M1AGfONDizrpIdjlRMWAzC6hqGYdfdAQC6aNffFzaXIyIj0evICgCo3nY0HywBAKTE+egKwITFUJJnzQMAtH28D756t+Ro+veDP36AgB4/AYiIOul5ZAUAWj7cDrW9Gba0LCRcdaPscKRjwmIg9uHZcE2YBUCg+cAW2eH0q77Vh1+UnpAdBhFRSCWH9TuyAnS04W/cuxlAR6M4xWKVHJF8TFgMRhtlaT60FarPIzma/v1i+wnWsxCR7gRUgR/88QPZYfSr/eN98Nd9AsWRiGHTC2SHowtMWAwmYeINsKZkQG1vClaO65XgyiEi0qFflp5AfatPdhj9ajr4vwCAYdd9HhZnouRo9IEJi8EoFiuSr78TANB0QL/Ft91x5RAR6cWWQ+fxos6nq/1NtWg7tRcAkDxzruRo9IMJiwENu+4OwGqDt+oEPFXHZYdzRezRQkR6oMcutqG0fFAKCBXOsZ+CfWSO7HB0gwmLAVkTU5E05bMAgKb9+i6+7Y4jLUQki95XBGmEUNF86K8AgGEz7pAcjb4wYTGo5Ou/AABo/ehdqJ4WydEMDHu0EJEMel8R1F17xT/gb6iG4kxC4jWfkR2OrjBhMShH9hTYR46D8HvRclSf+wv15cm3jrBHCxHFREAVeGLzEdlhDFjzP/4CABj2qdtgsbskR6MvTFgMSlEUJE3/HACg5fB2ydGEp6qhHeWn62SHQURxYPX2E3A3tssOY0ACrQ1oPbEbADDsukLJ0egPExYDS5r6T4BigefcEfjqzskOJyx/+ZDTQkQUXSWHq/DzbfpeEdRdy+FSQPXDkTUZjsyJssPRHSYsBmZLHglX7vUAjDfK8puyMyzAJaKoMUJzuO6EEGj6h1Zsy9GVUJiwGNywaR3TQs0fbocQquRoBo5LnYkomozQHK47r/tER2dbmxNJ194iOxxdYsJicAmTPw3FmYRA4wW0Vxrn14SGS52JKNKM0Byut5YPdwAAEibnsbNtH5iwGJzF7gz2ZGk5XCo5mvBxpIWIIskozeG6E2oALR+9CwBImnqb3GB0jAmLCSRNux0A0HpsF1Rvm+RoBocjLUQ0VEZpDtdb+5lDUFvqYXElI2HC9bLD0S0mLCbgHDMFtuHZEL52tB7bJTucQWFTOSIaCiM1h+ut5chOAEDilJuhWO2So9EvJiwmoCgKkj71TwCAlqPvSo5maNhUjojCZbTmcN2pPg9aj78HAEiaeqvkaPSNCYtJaHUs7WcOItDWKDmawWNTOSIKl5Gaw/XWdup9CG8brCmj4Bw7VXY4usaExSTsI8fCnjEBUANoPb5bdjhDwqZyRDRQRmsO11vLkR0AgKRrb4Wi8Cu5P7w6JqKNsrQe+7vkSIaGTeWIaCCM1hyut0B7M9o+3gsASPrUbXKDMQAmLCai7ezZXmHsaSEW4BLRQKzebqzmcL21HnsPCPhhTx8Px6hc2eHoHhMWE7GPGANH5iRAqGg9XiY7nCF7YvOHLMAlopACqsAr71XIDmNIWj/qGA1nse3AMGExGW2URfuHYGTuRg9Wbz8pOwwi0qHV20+goc24oytqezPaKw8B6Prcpv4NKmFZs2YNcnNz4XK5kJeXh/Ly8j6P/fDDD3HXXXchNzcXiqLghRdeGPI5qW+JU24GALSf+QcCrQ2Soxm6n287zi64RNTDlkPGLrQFgNZT7wNqAPb0cbCPGCM7HEMIO2HZsGEDiouLsXLlSuzfvx8zZsxAYWEhampqQh7f2tqKiRMn4umnn0ZWVlZEzkl9sw/PNtW0EMAuuETUxajdbHtrPd7R5DPx6pskR2IcYScszz//PO677z4sWbIEU6dOxdq1a5GYmIhXXnkl5PFz5szBz372M3zzm9+E0+mMyDmpf4naaqGP/iY5kshgES4RAcbuZtud6m1H+8cdSRcTloELK2Hxer3Yt28fCgoKuk5gsaCgoABlZYP7NT+Yc3o8HjQ2Nva4UZfgtFDlB6aYFtKwCJcofhm5m21v7af3Q/g9sKZmdvTPogEJK2Gpra1FIBBAZmZmj/szMzPhdrsHFcBgzrlq1SqkpqYGbzk5OYN6bbOyp2V1/CMQanCNvxmwCJcofhm5m21vXdNB+VAURXI0xmHIVULLly9HQ0ND8Hb27FnZIelO4lV5AIC2E3skRxJZP992nFNDRHHGDEW2GhHwdRTcgtNB4QorYUlPT4fVakV1dXWP+6urq/ssqI3GOZ1OJ1JSUnrcqKcELWE5vR/C75UcTWT94I8fcGqIKE6YpchW037mEISnBdak4XCOmSI7HEMJK2FxOByYPXs2SktLg/epqorS0lLk5+cPKoBonJMAR9YkWIeNgPC1o73SuK2rQ6lv9XFqiCgOmKXItjtt9WbC5E9z76AwhX21iouL8fLLL+O1117D0aNH8eCDD6KlpQVLliwBACxevBjLly8PHu/1enHw4EEcPHgQXq8X586dw8GDB3Hy5MkBn5PCpygWJFx1IwCg9aS5poUA4JX3PuYoC5GJBVSBJ98yR5GtRqgBtJ7o2JyW00Hhs4X7hIULF+LChQtYsWIF3G43Zs6ciZKSkmDRbGVlJSyWrjzo/PnzuP7664N/fvbZZ/Hss8/i1ltvxY4dOwZ0ThqcxKvy0HywBG0n9kB8/kFTFXc1tPmxevtJPFQwWXYoRBQF5afrUNVgjiJbjbfqBNTWeijOJLjGTZcdjuGEnbAAQFFREYqKikI+piUhmtzcXAhx5V/C/Z2TBsc1fgYUuxOB5ovwVp+CM+sq2SFF1M+3HcfkjCR84bps2aEQUYT99cPBrTzVs7bOYtuECbOgWAf19RvXOIFmYorNAdeEWQDMt1pIwy64ROaz5dB5vLqrQnYYEae1mUiYdIPkSIyJCYvJacubzVjHArALLpHZaIW2ZqtQ8zfXwVt9CgCQMGG25GiMiQmLySVMmgNAga/mY/gbzbs3E7vgEhmfmbrZ9tbeObriGD0Z1qQ0ucEYFBMWk7MmpsI55loAQKtJp4UAdsElMgMzdbPtre1U53TQxDmSIzEuJixxQFve3P7xPsmRRBe74BIZV8lh83Sz7U0EfGirOACA9StDwYQlDiRM7Ci8bT/7AYTfJzma6HryrSOcGiIyGDNPBQGA55OjEN42WBLT4DDZas1YYsISB+yjcmFJSoPwedB+zrwfCgBQ1dCO8tN1ssMgojCYeSoI6LaceeJsdrcdAl65OKAoFiTkdjTvaz9tnj05+rL1iPn6NxCZlZmngjTB5cwTOR00FExY4kSwH0scJCybDpzjtBCRAZix/X5vvno3fBfPAooFCROuv/ITqE9MWOKENsLiqzmNQPMlydFEVx03RyQyBDO23+9NW87sHDsVFtcwydEYGxOWOGFNSoMjcxIAoK3C/KMsP992HFsOnZcdBhH1w4zt93tr61ydyemgoWPCEkfiaVoIYNt+Ij0za/v97kTAj/azhwGA00ERwIQljiR0Jiztpw9ACFVyNNHHtv1E+mTW9vu9eapOdCxnTkiBPWOC7HAMjwlLHHGOmQLFkQC1rRHe6o9lhxMzbNtPpB9m77nSXfuZgwAA17jruJw5AngF44hitcM17joA8bG8WcO2/UT6YfaeK921VxwEALhyZ0qNwyyYsMSZhDirY9GwbT+RfPHQc0WjetvgOX8MABOWSGHCEmdcnYVfnnNHoXpaJUcTW2zbTyRPPE0FAYDn7IeA6oc1NRP2tCzZ4ZgCE5Y4Yx+eDWtqJqAG4Pkkfj48ALbtJ5IpnqaCAKCts34lYfwMuYGYCBOWOOQaNx1Ax2aI8YZt+4liL56mgjSsX4k8JixxKFh4W3lIciSxx7b9RLEVD+33ewu01MN3oQIA4OIIS8QwYYlDWsLidZ+C6mmRHE1ssW0/UWzFQ/v93trP/AMAYM+YAGtiquRozIMJSxyypaTDNnw0IFS0x1kdC8AVQ0SxtC0Op2G1hCVh/Ey5gZgME5Y45crpqGPxnIm/aSGAzeSIYiGgCvzxwDnZYcSUEAJtFQcAsH4l0piwxCnX+PitYwHYTI4oFlZvP4FLrT7ZYcSUv96NQOMFwGKDc+ynZIdjKkxY4pRTq2Op/hiB9mbJ0cjBqSGi6InHlUEA4Onc7NCZfTUsDpfkaMyFCUucsg0bAduIsQBER4OjOMWpIaLIi7cmcd21d36ecnQl8piwxDHXuGkA4ndaCODUEFE0xFuTuO48n3SMsLhypkmOxHyYsMSxeO7H0h2nhogiJ16nggDA31QLf70bUCxwjrlWdjimw4Qljmkdb301FQi0NUmORi7uM0Q0dPE8FQQgOL3uyJwIizNRcjTmw4QljlmThsM+chw66ljir01/d9xniGjo4nkqCADaP2H9SjQxYYlzTm1focr4TlgA7jNENBTxPBWk0VYIsX4lOpiwxDntH1a87dwcCvcZIhqceJ8KAoBAawN8tZUAAOfYqZKjMScmLHHOObajMMxbcxqqp1VyNHJxnyGiwYn3qSCg60effeQ47h8UJUxY4pwtOR3W1ExAqPBUHZcdjnRcMUQUHk4FdQjWr+SwfiVamLAQXJ3L7zyfxG8Due7YTI5oYAKqwJNvxfdUkEZbIcT6lehhwkLB+VbPJ0clR6IPbCZHNDDlp+tQ1RDfU0EAoHpa4a0+BYArhKKJCQsFGxx5qo5BqAHJ0egDp4aIrmwbV9YBADznjgJChS0tC7aUdNnhmBYTFoI9fRwURyKEtw2+CxWyw9ENTg0R9S2gCvzxwDnZYegC+6/EBhMWgmKxwjlmCgCgncubgzg1RNS31dtP4FKrT3YYuqCtEGLCEl1MWAhA9zoWJizdcWqI6HJcGdRFqAF4qzquhYv7B0UVExYCALjGdCYs51h42xunhoi6sElcT96a0xB+DyzOJNhGjpEdjqkxYSEAgGP01YBiQaCpFv7GGtnh6Aqnhoi6sElcT97zHwEAHNnXQFH4lRpNg7q6a9asQW5uLlwuF/Ly8lBeXt7v8Rs3bsSUKVPgcrkwffp0bNmypcfjzc3NKCoqwtixY5GQkICpU6di7dq1gwmNBsnicMGROQkAlzeHwqkhIk4FheI5fwwA4MyeIjkS8ws7YdmwYQOKi4uxcuVK7N+/HzNmzEBhYSFqakL/Kt+1axcWLVqEe++9FwcOHMD8+fMxf/58HD58OHhMcXExSkpK8Prrr+Po0aN4+OGHUVRUhM2bNw/+nVHYtOXNLLwN7cm3jnBqiOIWp4JC85zrGGFxZl8jORLzCzthef7553HfffdhyZIlwZGQxMREvPLKKyGPf/HFFzF37lwsW7YM1157LZ566inMmjULq1evDh6za9cu3HPPPbjtttuQm5uL+++/HzNmzLjiyA1FVrDw9hw/lEKpamhH+ek62WEQScGpoMsFWhvgr+8YeWXCEn1hJSxerxf79u1DQUFB1wksFhQUFKCsrCzkc8rKynocDwCFhYU9jr/pppuwefNmnDt3DkIIvPPOOzh+/DjuuOOOkOf0eDxobGzscaOh00ZYfBfOxP1GiH3ZykZZFIc4FRSap7N+xT4yBxbXMMnRmF9YCUttbS0CgQAyMzN73J+ZmQm3O/QHudvtvuLxv/zlLzF16lSMHTsWDocDc+fOxZo1a3DLLbeEPOeqVauQmpoavOXk5ITzNqgPtuSR3AjxCv508DynhSiucL+gvmnTQQ7Wr8SELkqaf/nLX2L37t3YvHkz9u3bh+eeew5Lly7Ftm3bQh6/fPlyNDQ0BG9nz56NccTm5Rx9NQDAy4QlpIstXk4LUVzhfkF9CxbcjmHCEgu2cA5OT0+H1WpFdXV1j/urq6uRlZUV8jlZWVn9Ht/W1obHHnsMb775JubNmwcAuO6663Dw4EE8++yzl00nAYDT6YTT6QwndBog5+ir0frR3zjC0o+tR9zInzRSdhhEMVHTxGQllI6GcR2fk6xfiY2wRlgcDgdmz56N0tLS4H2qqqK0tBT5+fkhn5Ofn9/jeADYunVr8HifzwefzweLpWcoVqsVqqqGEx5FgCO7c4Tl/DEIwamPUF55r4JLnCluVNS2yA5Bl3wXzkD42qE4EmFPHyc7nLgQ1ggL0LEE+Z577sENN9yAG2+8ES+88AJaWlqwZMkSAMDixYsxZswYrFq1CgDw0EMP4dZbb8Vzzz2HefPmYf369di7dy9+9atfAQBSUlJw6623YtmyZUhISMD48eOxc+dO/OY3v8Hzzz8fwbdKA+HInNTRQK7lEgJNF7nzaB+e2PwhPj81C1aLIjsUoqgJqALr9lTKDkOXtIJb5+ir2TAuRsJOWBYuXIgLFy5gxYoVcLvdmDlzJkpKSoKFtZWVlT1GS2666SasW7cOjz/+OB577DFMnjwZmzZtwrRp04LHrF+/HsuXL8fdd9+Nuro6jB8/Hv/5n/+JBx54IAJvkcJhsbtgH5ULX83H8FYdZ8LSB6377UMFk2WHQhQ1q7efQHWTR3YYuhRMWFi/EjOKMMG4f2NjI1JTU9HQ0ICUlJSInbfV68fUFX+J2PmM4mLJajT/owQpeV/H8Nv+VXY4urb2X2Zh7rTRssMgiriSw1V44PX9ssPQrXMv/x/4684h4+srkTBpjuxwYubIjwuR6Ah7rKNP4Xx/cxyLLuPoXCnEwtsr48aIZEbsatu/QFsj/HXnAHTsIUSxwYSFLuPUCm/dJyDUgORo9I0bI5IZsatt/7znO37M2YZnw5oQuVF96h8TFrqMfWQOFLsLwtsG38VPZIeje9wYkcyEXW2vzOPuuD5a3yqKDSYsdBnFYoUj6yoAbCA3UNwYkcyAU0ED460+BQBwZLHoPpaYsFBITtaxhIUbI5IZcCpoYLxVHSMsjtFXSY4kvjBhoZAcbNEfNm6MSEbGqaCB8TfXIdB8EYACR8ZE2eHEFSYsFFKw8PZCBVQf+zAMBDdGJKPiBocD53V3FNnbR+bA4kiQHE18YcJCIVmTR8GaNBxQA/DVfCw7HEPgxohkVNzgcOC8bk4HycKEhUJSFIX9WAaB00JkRNzgcOC0ERYW3MYeExbqU7Dw9jwTloHixohkRNzgcGCEEMGExZnFEZZYY8JCfQoube5cwkcDw+63ZCTc4HDgAs0XEWi5BCgW2DMmyA4n7jBhoT45MicBAPx156B6WiVHYxzsfktGwg0OBy5YcJs+Dha7S3I08YcJC/XJmpQGa3I6AAHvhdOywzEUdr8lI+BS5vAE+69wOkgKJizUL0dmR58BbzVXCoWLU0OkZ+xqGz5PtVa/woJbGZiwUL+0aSGvm3Us4eLUEOkZu9qGp3vBLUdY5GDCQv0KJiw1TFgGg1NDpEecCgpfoOkC1NYGwGKFfVSu7HDiEhMW6peWsPhqKyH8XsnRGBM3RiQ9YVfbwfFWaQW342GxOyVHE5+YsFC/rMnpsCSkAGoA3gtnZIdjSNwYkfSEXW0Hx9PZ4Zb9V+RhwkL9UhSla1qI/VgGjZ1ESS+2sRvzoATrV0az4FYWJix0RUxYhq6iln1sSL6Sw1X4/9+rkB2G4Qgh4O3cU407NMvDhIWuiAnL0L1RfoZ1LCQVlzEPXqDlUkfBrWKBfdR42eHELSYsdEVaLxbfhQoINSA5GmPiEmeSjcuYB8/X2YfKPmIMO9xKxISFrsg2fDQURwKE3wvfxbOywzEsLnEmWbiMeWi06SDuHyQXExa6IkWxBOdt2fF2aLjEmWKNy5iHzlvTsTUJ61fkYsJCA8I6lsjgEmeKNS5jHrquhIUjLDIxYaEBcWQxYYmUrVxWSjHEJfVDo3rb4a87B4AjLLIxYaEB6Rph+RhCqJKjMbZX3qtgLQvFTEVti+wQDM1XewaAgCUpDdZhw2WHE9eYsNCA2EfmQLE5ILyt8F/il+1QcSdnioWAKrBuT6XsMAwtOB00itNBsjFhoQFRLNZg/wHtHzANHpc5Uyys3n4C1U0e2WEYWrBhXCang2RjwkIDpv3C8HFPoYjgMmeKJi5ljgwfC251gwkLDVhwhOUCR1gihcucKRq4lDkyhFCDI8p2FtxKx4SFBowjLJHHZc4UDVzKHBn+S1UQvnbAaod9xBjZ4cQ9Jiw0YNoIi7++Cqq3TXI05sFlzhRp3JE5MroKbnOhWKySoyEmLDRg1sRUWIeNAMBRlkjiMmeKJO7IHDlsGKcvTFgoLPZRuQAA74UKqXGYDWtZKBK4I3Nk+biHkK4wYaGwODoTFh8TlohiLQtFAndkjqzgCAuXNOsCExYKC0dYooe1LDQUXMYcWYG2RgSaagGwaZxeMGGhsDgycgF0jLAIwSmMSPrTwfOcFqJB4TLmyNNGV2xpWbA4EyVHQwATFgqTfUQOoFigtjcj0HRRdjimcrHFy2khGhQuY448bWGBNqpM8jFhobAoNjvsI8cCYB1LNHBaiAaDOzJHXsemh4AjfbzkSEjDhIXCxjqW6OESZxoM7sgcecERlvRxkiMhDRMWCpsjmLCwRX80cCdnCgd3ZI48IQS8tR3XVGuYSfIxYaGwaT0J2DwuOriTM4WDOzJHXqCpFsLbClisbMmvI4NKWNasWYPc3Fy4XC7k5eWhvLy83+M3btyIKVOmwOVyYfr06diyZctlxxw9ehRf/vKXkZqaiqSkJMyZMweVlfzVoEeOzl8cvotnIQI+ydGYE3dypoHgUuboCE4HjRgDxWqXHA1pwk5YNmzYgOLiYqxcuRL79+/HjBkzUFhYiJqampDH79q1C4sWLcK9996LAwcOYP78+Zg/fz4OHz4cPObUqVO4+eabMWXKFOzYsQOHDh3Cj370I7hcrsG/M4oaa/IoKM4kQA3AV3dOdjimxe631B8uZY4eb61Wv8LpID0JO2F5/vnncd9992HJkiWYOnUq1q5di8TERLzyyishj3/xxRcxd+5cLFu2DNdeey2eeuopzJo1C6tXrw4e88Mf/hBf+MIX8Mwzz+D666/HpEmT8OUvfxkZGRmDf2cUNYqidNWx1LCOJVrY/Zb6w6XM0ePT6ldYcKsrYSUsXq8X+/btQ0FBQdcJLBYUFBSgrKws5HPKysp6HA8AhYWFweNVVcWf//xnXH311SgsLERGRgby8vKwadOmPuPweDxobGzscaPYsrNFf0xwmTP1hTsyR09wSTMLbnUlrISltrYWgUAAmZmZPe7PzMyE2x36H4/b7e73+JqaGjQ3N+Ppp5/G3Llz8de//hVf/epX8bWvfQ07d+4Mec5Vq1YhNTU1eMvJyQnnbVAEaB1vvTUVUuMwOy5zplC4I3P0CKHCV3sWAKeE9Eb6KiFVVQEAX/nKV/C9730PM2fOxA9+8AN88YtfxNq1a0M+Z/ny5WhoaAjezp49G8uQCV3/kH0XWRgdTQpYy0I9sXYluvz11RB+D2C1w5aWJTsc6iashCU9PR1WqxXV1dU97q+urkZWVuj/sFlZWf0en56eDpvNhqlTp/Y45tprr+1zlZDT6URKSkqPG8WWfWTHqFag8QJUT6vkaMxLgLUs1BNrV6Krq8PtOCgWq+RoqLuwEhaHw4HZs2ejtLQ0eJ+qqigtLUV+fn7I5+Tn5/c4HgC2bt0aPN7hcGDOnDk4duxYj2OOHz+O8eM5HKdX1oRkWJOGAwB8dZ9Ijsb82HqdNKxdiS52uNUvW7hPKC4uxj333IMbbrgBN954I1544QW0tLRgyZIlAIDFixdjzJgxWLVqFQDgoYcewq233ornnnsO8+bNw/r167F371786le/Cp5z2bJlWLhwIW655Rb80z/9E0pKSvDWW29hx44dkXmXFBX29BwEWi7BV3sWztFXyw7H1CpqOYpFrF2JhWCHW9av6E7YNSwLFy7Es88+ixUrVmDmzJk4ePAgSkpKgoW1lZWVqKrqKhK86aabsG7dOvzqV7/CjBkz8Pvf/x6bNm3CtGnTgsd89atfxdq1a/HMM89g+vTp+J//+R/84Q9/wM033xyBt0jRYh/Z8QuEdSzR90b5GdaxxDnWrsSGNiXElvz6E/YICwAUFRWhqKgo5GOhRkUWLFiABQsW9HvOb3/72/j2t789mHBIEm3I1HeRRc/RprXrf6hgsuxQSBLWrkSfCPjhu9gxxc1dmvVH+iohMi77yLEAEFwCSNHFdv3xjbUr0ee7dB5Q/VAcCbCmjJIdDvXChIUGTZsS8te7ofq4+VoscIlzfAqoAm8e5DYY0RbscDtyHBRFkRwN9caEhQbNkpgKS0IKAAE/9xSKCS5xjk/lp+tQ18KNRqNNWyHEDrf6xISFBk1RlGA/Fhbexg7b9ccfTgfFho+bHuoaExYaEnt6Z8LCOpaYYbv++MKlzLHj5aaHusaEhYaka2kzE5ZYYbv++MGlzLEjAn746zt+CGg/xEhfmLDQkGhTQtovE4o+tuuPH1zKHDv++ipADXSsEBo2UnY4FAITFhoS7ZeI/9J5iACLAmOJtSzmx9qV2NH6r9hHjOUKIZ1iwkJDYh02EoojERAqfHXnZYcTV1jLYm6sXYktbU80rb8U6Q8TFhoSRVG6Cm9ZxxJTrGUxL9auxF73ERbSJyYsNGQsvJWDtSzmxdqV2NM+v2wcYdEtJiw0ZI7g0mYW3spQ08QvNrNh7UpsCSE4wmIATFhoyLqax3GERYaK2lbZIVAEsXYl9gItlyC8rYBigX14tuxwqA9MWGjIgrs2130CoQYkRxN/3ig/wzoWkwioAk9sZu1KrPk7R1dsaZlQbHbJ0VBfmLDQkFlTRkGxO4GAH/56DmXHmrvRg9XbT8oOgyJg9fYTcDdyii/WgiuEOB2ka0xYaMgUxdI1LcQ6Fil+vu04lzgbXMnhKvx82wnZYcQlbTpb+xwjfWLCQhFhGzEGAOC7xF2bZeESZ+PiMma5tIJbG0dYdI0JC0WENpTqu8iERRYucTYuLmOWi03jjIEJC0WEVlnv5wiLVGzXb0xcmi6P6m1HoPECACYseseEhSJC+4fO9vxysV2/MVXUtsgOIW5pP7IsCSmwJqRIjob6w4SFIsLWOcKittZDbW+WHE38Yrt+4wmoAm+Us1hdlmDDOI6u6B4TFooIiyMB1mEjAAC+Ok4LycJ2/cZTfroO7kaP7DDiVnCFEAtudY8JC0WMVmHPhEU+1rIYB9vwy8URFuNgwkIRYx/RWXjLhEU61rIYA9vwy6etELKxB4vuMWGhiLEP7+zFwoRFOtay6B97r8gn1EDw84pTQvrHhIUiRtuWnc3j5GMti/6x94p8/sYLQMAHWO2wpWbIDoeugAkLRUxXL5bzEEKVHA0BrGXRM9auyOcPFtyOgWKxSo6GroQJC0WMLS0LsFghfB4Emi7KDofAWha9Yu2KPgSngzp/bJG+MWGhiFEsVtjSRgNgHYtesJZFf1i7oh++Sx2NLrW90EjfmLBQRHGlkL6wlkV/WLuiH/7Oztz24aMlR0IDwYSFIsrOXiy6xL1q9IO1K/rhq++YLrVxSsgQmLBQRGn/8Jmw6EtFbavsEAisXdET4fd1bXo4nFNCRsCEhSJK6xbJXZv15Y3yM6xjkYy1K/rir3cDQoXiSIAlKU12ODQATFgoorRfKv6GGgi/T3I0pHE3erB6+0nZYcQ11q7oi1Zwax+eDUVRJEdDA8GEhSLKkpQGxZEACDU4P0z68PNtx7nEWSLWruiLNgqsrWwk/WPCQhGlKEqw8JYrhfSHS5zlCKgCbx7kvwc98V3qLLgdwYJbo2DCQhGnfQCw8FZ/uMRZjvLTdahr4RSpnvi7TQmRMTBhoYjjJoj6xnb9scfpIP0JNo1jwmIYTFgo4uydXSM5JaRPbNcfW1zKrD/C70WgsRYAR1iMhAkLRZzW5lr7BUP6wnb9scOlzPrUUb8iOpY0J6bKDocGiAkLRZzW5lptrYfqYcMyvWG7/tjhUmZ98neuYOSSZmNhwkIRZ3EmwZKQAgDwN3DuXq9YyxJ9rF3RJ18d61eMiAkLRUVw1+ZLrJXQK9ayRBdrV/SLK4SMaVAJy5o1a5CbmwuXy4W8vDyUl5f3e/zGjRsxZcoUuFwuTJ8+HVu2bOnz2AceeACKouCFF14YTGikE7bhWQC6hl5Jf1jLEj2sXdG34Aoh9mAxlLATlg0bNqC4uBgrV67E/v37MWPGDBQWFqKmpibk8bt27cKiRYtw77334sCBA5g/fz7mz5+Pw4cPX3bsm2++id27dyM7m3+JjM7eOcLiv8Qhcb1iLUv0sHZF37QRFlsav2uMJOyE5fnnn8d9992HJUuWYOrUqVi7di0SExPxyiuvhDz+xRdfxNy5c7Fs2TJce+21eOqppzBr1iysXr26x3Hnzp3Dd7/7Xfz2t7+F3W4f3Lsh3dCmhPz1XCmkdzVN/GKNNNau6Jfqa0egqXNJM0dYDCWshMXr9WLfvn0oKCjoOoHFgoKCApSVlYV8TllZWY/jAaCwsLDH8aqq4lvf+haWLVuGT33qU1eMw+PxoLGxsceN9MU2XKth4Qe33lXUciVXJLF2Rd/89R2fSUq3xQFkDGElLLW1tQgEAsjMzOxxf2ZmJtzu0F9Mbrf7isf/9Kc/hc1mw7//+78PKI5Vq1YhNTU1eMvJyQnnbVAMaFNCgaZa7tqsc+vfr2QdS4SwdkX//Nyl2bCkrxLat28fXnzxRbz66qsD/suzfPlyNDQ0BG9nz56NcpQULktSGhS7CxAq/A3VssOhfrCOJXJYu6J/bMlvXGElLOnp6bBaraiu7vkFVF1djaysrJDPycrK6vf4v/3tb6ipqcG4ceNgs9lgs9lw5swZPPLII8jNzQ15TqfTiZSUlB430hdFUWBL40oho2BPlshgPZD++eu4pNmowkpYHA4HZs+ejdLS0uB9qqqitLQU+fn5IZ+Tn5/f43gA2Lp1a/D4b33rWzh06BAOHjwYvGVnZ2PZsmX4y1/+Eu77IR0J1rHU88tQ79iTJTIqaltkh0BX0DXCMlpyJBQuW7hPKC4uxj333IMbbrgBN954I1544QW0tLRgyZIlAIDFixdjzJgxWLVqFQDgoYcewq233ornnnsO8+bNw/r167F371786le/AgCMHDkSI0eO7PEadrsdWVlZuOaaa4b6/kgie9potKFrzpj0S+vJ8vmpWbBaOK8/GAFV4I3yStlh0BX4L3W15SdjCTthWbhwIS5cuIAVK1bA7XZj5syZKCkpCRbWVlZWwmLpGri56aabsG7dOjz++ON47LHHMHnyZGzatAnTpk2L3LsgXdJ+wfg5wqJ73Xuy5E8aecXj6XLlp+vgbvTIDoP6IfxeBJovAuAIixGFnbAAQFFREYqKikI+tmPHjsvuW7BgARYsWDDg81dUVAwmLNIZtuc3nq1H3ExYBom9V/RPWwCgOBK4pNmApK8SIvMKFt02uCHUgORoaCBYyzI47L1iDNpory01k0uaDYgJC0WNLWUUYLECAX9wGJb0jfsLhY+9V4xDWwCg/ZgiY2HCQlGjWKywpXbUNnFPIWPg/kLhY+8V4/AzYTE0JiwUVcE6FvZiMRT2ZRk41q4Yh1bDYmfCYkhMWCiq7MM761hYeGsorGUZGNauGEtXDQsTFiNiwkJRpW3fzm63xsJalitj7YqxCCGCIyycEjImJiwUVTZthIW9WAyFtSxXxtoVY1HbGiG8bQAAW2qG5GhoMJiwUFR19WI5DyH4a91ouDdO31i7YizajybrsJFQbA7J0dBgMGGhqNJWCQlvG9S2RsnRULgqaltlh6BLrF0xHq4QMj4mLBRVFrsT1uR0ACy8NaL171eyjqUX1q4YE+tXjI8JC0Wd9gHBpc3GwzqWy7F2xZg4wmJ8TFgo6rpa9FdLjoQGgz1ZemLtijGxy63xMWGhqAt2u61nwmJE7MnShbUrxsUeLMbHhIWiTktYAo1MWIyIPVk6sHbFuETAh0BTLQDAlpYpORoaLCYsFHVazwOOsBgTe7J0YO2KcfkbLwBChWJzwJo0XHY4NEhMWCjqtCFYf1MthBqQHA0NVrzXsrAnjXFpP5ZsqVlQFEVyNDRYTFgo6qzDhgMWG6AGEGi6KDscGqR4r2WpqG2RHQINkr9BK7jldJCRMWGhqFMsVthSRwHgSiEji+daloAq8EZ5pewwaJC4pNkcmLBQTNhSOlcKNdRIjoQGK55rWcpP18Hd6JEdBg0SExZzYMJCMREsvG2I7zoIM4jHWhb2XjE2JizmwISFYqKreRxHWIwu3mpZ2HvF+Lp6sLCGxciYsFBMdI2wsIbF6OKploW9V4wv0N4M1dNRMM2mccbGhIViItjtlgmL4cVTLQt7rxifNrpiSUqDxeGSHA0NBRMWigmr1u226SJEwC85GoqEeKhlYe2K8WkJi52jK4bHhIViwpo0HIrNAQgV/s4W2WRsZq9lYe2KOWiF/lb2YDE8JiwUE4qiwJrCOhYzMXMtC2tXzMPfcAEAC27NgAkLxQz3FDIXM9eysHbFPAKNHSsTbZ0/mMi4mLBQzAR3beYIi6mYcY8d1q6Yh9ZKwZYySnIkNFRMWChmuFLInCpqW2WHEFGsXTEPIQT82ghLKkdYjI4JC8VMV8LC5nFmsv79StPUsbB2xVxUTwuEtw0AYOUIi+ExYaGYYXt+czJTHQtrV8wl0PnjyJKYCoudPViMjgkLxYzWnj/QXAfh90mOhiLJLD1ZzFiPE8/8LLg1FSYsFDOWhBQodieArg8SMgez9GSpqG2RHQJFEAtuzYUJC8WMoiiwpbDw1ozM0JMloAq8UV4pOwyKIO1zxsqCW1NgwkIxZUtj4a0ZmaEnS/npOrgbPbLDoAgKNHY2jeOUkCkwYaGYYuGtuRm5loW9V8yHS5rNhQkLxRSXNpubUWtZ2HvFnLpqWJiwmAETFoopbddmtuc3JyPWsrD3ijmp3naobY0AWMNiFkxYKKbY7dbcjFjLwt4r5qTVryiORFicSZKjoUhgwkIxpS0vVFvrIfxeydFQtBiploW1K+bUvX5FURTJ0VAkMGGhmLIkpECxdfZiaaqVHA1Fi1FqWVi7Yl5dTePYg8UsmLBQTCmKEtzTw985ZEvmY4RaFtaumFuw4Jb1K6bBhIViTvvFE2jkCItZGaGWhbUr5qaNsFi5Qsg0BpWwrFmzBrm5uXC5XMjLy0N5eXm/x2/cuBFTpkyBy+XC9OnTsWXLluBjPp8Pjz76KKZPn46kpCRkZ2dj8eLFOH/+/GBCIwOwJqcDYHv+eKDnvXlYu2JugQY2jTObsBOWDRs2oLi4GCtXrsT+/fsxY8YMFBYWoqYm9JfPrl27sGjRItx77704cOAA5s+fj/nz5+Pw4cMAgNbWVuzfvx8/+tGPsH//fvzxj3/EsWPH8OUvf3lo74x0q2uEhVNCZldR2yo7hJBYu2J+bBpnPmEnLM8//zzuu+8+LFmyBFOnTsXatWuRmJiIV155JeTxL774IubOnYtly5bh2muvxVNPPYVZs2Zh9erVAIDU1FRs3boV3/jGN3DNNdfg05/+NFavXo19+/ahspL7epiRlrCw6Nb81r9fqbs6FtaumJ8I+BFo7piO5AiLeYSVsHi9Xuzbtw8FBQVdJ7BYUFBQgLKyspDPKSsr63E8ABQWFvZ5PAA0NDRAURSkpaWFfNzj8aCxsbHHjYzDyhGWuKHHOhbWrpifv6kWECpgtcOSlCo7HIqQsBKW2tpaBAIBZGZm9rg/MzMTbnfo+WC32x3W8e3t7Xj00UexaNEipKSkhDxm1apVSE1NDd5ycnLCeRskma3bKiEh9PXrmyJPbz1Z9FxXQ5ER6LakWVG4tsQsdPVf0ufz4Rvf+AaEEHjppZf6PG758uVoaGgI3s6ePRvDKGmotKJb4WuH2t4sORqKNr31ZKmobZEdAkWZnwW3phRWwpKeng6r1Yrq6p5t1aurq5GVlRXyOVlZWQM6XktWzpw5g61bt/Y5ugIATqcTKSkpPW5kHBa7E5bEjmHaQBOnhcxOTz1ZAqrAG+WsjTM7FtyaU1gJi8PhwOzZs1FaWhq8T1VVlJaWIj8/P+Rz8vPzexwPAFu3bu1xvJasnDhxAtu2bcPIkSPDCYsMyBZc2syExez01JOl/HQd3I0e2WFQlGlN46zscmsqYU8JFRcX4+WXX8Zrr72Go0eP4sEHH0RLSwuWLFkCAFi8eDGWL18ePP6hhx5CSUkJnnvuOXz00Ud44oknsHfvXhQVFQHoSFa+/vWvY+/evfjtb3+LQCAAt9sNt9sNr5d7zZgVC2/jjx5qWdh7JT4EGjtG9TnCYi62cJ+wcOFCXLhwAStWrIDb7cbMmTNRUlISLKytrKyExdKVB910001Yt24dHn/8cTz22GOYPHkyNm3ahGnTpgEAzp07h82bNwMAZs6c2eO13nnnHdx2222DfGukZ12Ft1zaHC9eea8CN04YgbnTRkt5ffZeiR+sYTGnsBMWACgqKgqOkPS2Y8eOy+5bsGABFixYEPL43NxcrhSJQzbuJxR3tFqWz0/NgtUS291z2Xslfgihwt9ZG8cRFnPR1Sohih/WZE4JxRuZtSzsvRI/1NYGIOAHoMA6jPWQZsKEhaTgCEv8klHLwtqV+KFNM1uHDYdiHdQkAukUExaSIlh023wRQg1IjoZiKdZ9WVi7El8CnVt+aP2eyDyYsJAU1mHDAYsNEGpwzw+KD7Hsy8Lalfij7VFmY8JiOkxYSApFscCa3DG/zGmh+BLLWhbWrsQfjrCYFxMWksbGXixxLRZ7+rB2Jf4Ea1iYsJgOExaSJlh4y/b8camitjWq52ftSnzSRlhsKUxYzIYJC0mj/QLiCEt8Wv9+ZdTqWFi7Er/8nBIyLSYsJA2XNse3aNaxsHYlPgmhItB0EQBHWMyICQtJw/b8FK06lljUx5D+qC0NgNrZNC5phOxwKMKYsJA03ACRolXHkpHsisp5Sd+C00FsGmdKTFhIGm2ERW1vgurlL+J49MK241FpInepxRPxc5L+BToL+Fm/Yk5MWEgaizMJiiMRQNcHDcWfSDeRC6gCT/35aMTOR8bh1+pXmLCYEhMWkoqFt/EtGk3kWHAbv7TpZY6wmBMTFpKKCQsBkd0Qkc3i4pefK4RMjQkLSWVNYS8WityGiGwWF99Yw2JuTFhIKu2DRftlRPEpEhsislkcaZ8j1uRRkiOhaGDCQlJpxXFaO22KT5GoZWHtSnzr2TRupORoKBqYsJBUwfb8HGEhDK2WhbUr8S3YNE6xsGmcSTFhIalsyR2/hPzNTFho8LUsrF0hbRNVa1Iam8aZFBMWkso6rCNhEZ4WqJ7o7t5L+jeYWhbWrhDQNa3M+hXzYsJCUlmciVCcSQCAAEdZ4t5gallYu0JA155k2qgtmQ8TFpIuOC3EOhbqFE4tC2tXCOg2wpLCERazYsJC0mnTQlwpRJqB1rKwdoU02saHHGExLyYsJF1XLxYmLNRhILUsrF2h7ljDYn5MWEg69mKh3gZSy8LaFepOq2Fhl1vzYsJC0lmTtSkh1rBQTzVNfSckrF0hjRBqsGif+wiZFxMWks7GKSHqQ0Vt6KXurF2h7gIt9YAa6GgaN4xN48yKCQtJF9wAkSMs1Mv69ysvq2Nh7Qr1FqxfSRoOxWKVHA1FCxMWkk6bc1bbGqH6PJKjIT0JVcfC2hXqLcD6lbjAhIWksziToNidAIBA8+A3vyNz6t2Tpb+6FopPwSXNrF8xNSYsJJ2iKF2bIDZekBwN6U3vniwVtS0SoyE96lrSzITFzJiwkC7YOnsnaBuYEWm692QJqAJvlFfKDol0xt/5Q8eWkiE5EoombmlJuqC10/ZzhIV66d2Txd3IOifqyd9YAwCwsS2/qTFhIV3QPmg4JUR92XrEDUV2EKRL2ucG9xEyNyYspAtawuJvqJEcCenVK+y7QiGIgA+B5ksAOMJidqxhIV3glBARDUbHLu8CsNphSUyVHQ5FERMW0gVbakexXKDpAoToe8M7IqLuAsGC21FQFE4amhkTFtIFrT2/8HmgtjdJjoaIjIIFt/GDCQvpgmJzwJKUBoCFt0Q0cH4W3MYNJiykGyy8JaJwBdiDJW4wYSHdCDaP4wgLEQ2Qv1sNC5nboBKWNWvWIDc3Fy6XC3l5eSgvL+/3+I0bN2LKlClwuVyYPn06tmzZ0uNxIQRWrFiB0aNHIyEhAQUFBThx4sRgQiMDs2qFt0xYiGiA2IMlfoSdsGzYsAHFxcVYuXIl9u/fjxkzZqCwsBA1NaGH8Xft2oVFixbh3nvvxYEDBzB//nzMnz8fhw8fDh7zzDPP4Be/+AXWrl2LPXv2ICkpCYWFhWhv5yZn8cTGpc1EFAYhBEdY4ogiwlxDmpeXhzlz5mD16tUAAFVVkZOTg+9+97v4wQ9+cNnxCxcuREtLC95+++3gfZ/+9Kcxc+ZMrF27FkIIZGdn45FHHsH3v/99AEBDQwMyMzPx6quv4pvf/OYVY2psbERqaioaGhqQkpISztvpV4vHh1k/fCti56P+NZ3Yg6rNz8A1ejLG/fPTssMhIp0LtDXh1H/9KwDgqn9fB0vnru8UPfv/80tIctojdr5wvr/D6nTr9Xqxb98+LF++PHifxWJBQUEBysrKQj6nrKwMxcXFPe4rLCzEpk2bAACnT5+G2+1GQUFB8PHU1FTk5eWhrKwsZMLi8Xjg8XTtJ9LY2BjO2xgw0daGTW//MCrnpst92N6OBQCs1aeQ98q/yQ6HiHSuIRDAKQAjrVZs/suPZYcTF8TjBUAEE5ZwhJWw1NbWIhAIIDMzs8f9mZmZ+Oijj0I+x+12hzze7XYHH9fu6+uY3latWoUnn3wynNAHJdHBnQtiaYzdDiuAFlXF/710SXY4RGQQ4x0O2SHEDZnfi4b8Rl6+fHmPUZvGxkbk5ORE/HWUhARcs39fxM9Lffv9n/+M8vfflx0GSXauvg1vHzp/xeO+eF02xqQlxCAi0iur1Ypv3HUXrpk6VXYocUFJkPfvLayEJT09HVarFdXV1T3ur66uRlZWVsjnZGVl9Xu89r/V1dUYPXp0j2NmzpwZ8pxOpxNOZ/TnKhVFgZKYGPXXoS7zFyzA/AULZIdBkgVUgZt/uh3uhnaEKrJTAGSluvDKo5+D1cJ27ETxIKxVQg6HA7Nnz0ZpaWnwPlVVUVpaivz8/JDPyc/P73E8AGzdujV4/IQJE5CVldXjmMbGRuzZs6fPcxKRuVktClZ+qeMXc+90RPvzyi9NZbJCFEfCXtZcXFyMl19+Ga+99hqOHj2KBx98EC0tLViyZAkAYPHixT2Kch966CGUlJTgueeew0cffYQnnngCe/fuRVFREYCOUYyHH34YP/nJT7B582Z88MEHWLx4MbKzszF//vzIvEsiMpy500bjpX+ZhaxUV4/7s1JdeOlfZmHutNF9PJOIzCjsGpaFCxfiwoULWLFiBdxuN2bOnImSkpJg0WxlZSUslq486KabbsK6devw+OOP47HHHsPkyZOxadMmTJs2LXjMf/zHf6ClpQX3338/6uvrcfPNN6OkpAQul+uy1yei+DF32mh8fmoWyk/XoaapHRnJLtw4YQRHVojiUNh9WPQoWn1YiIiIKHrC+f7mXkJERESke0xYiIiISPeYsBAREZHuMWEhIiIi3WPCQkRERLrHhIWIiIh0jwkLERER6R4TFiIiItI9JixERESke2G35tcjrVlvY2Oj5EiIiIhooLTv7YE03TdFwtLU1AQAyMnJkRwJERERhaupqQmpqan9HmOKvYRUVcX58+eRnJwMRYnspmiNjY3IycnB2bNnuU9RFPD6Rh+vcXTx+kYfr3F0yby+Qgg0NTUhOzu7x8bJoZhihMVisWDs2LFRfY2UlBT+Q4kiXt/o4zWOLl7f6OM1ji5Z1/dKIysaFt0SERGR7jFhISIiIt1jwnIFTqcTK1euhNPplB2KKfH6Rh+vcXTx+kYfr3F0GeX6mqLoloiIiMyNIyxERESke0xYiIiISPeYsBAREZHuMWEhIiIi3WPCcgVr1qxBbm4uXC4X8vLyUF5eLjskQ3riiSegKEqP25QpU4KPt7e3Y+nSpRg5ciSGDRuGu+66C9XV1RIj1rd3330XX/rSl5CdnQ1FUbBp06YejwshsGLFCowePRoJCQkoKCjAiRMnehxTV1eHu+++GykpKUhLS8O9996L5ubmGL4LfbvSNf7Xf/3Xy/5Oz507t8cxvMZ9W7VqFebMmYPk5GRkZGRg/vz5OHbsWI9jBvK5UFlZiXnz5iExMREZGRlYtmwZ/H5/LN+KLg3k+t52222X/R1+4IEHehyjp+vLhKUfGzZsQHFxMVauXIn9+/djxowZKCwsRE1NjezQDOlTn/oUqqqqgre///3vwce+973v4a233sLGjRuxc+dOnD9/Hl/72tckRqtvLS0tmDFjBtasWRPy8WeeeQa/+MUvsHbtWuzZswdJSUkoLCxEe3t78Ji7774bH374IbZu3Yq3334b7777Lu6///5YvQXdu9I1BoC5c+f2+Dv9xhtv9Hic17hvO3fuxNKlS7F7925s3boVPp8Pd9xxB1paWoLHXOlzIRAIYN68efB6vdi1axdee+01vPrqq1ixYoWMt6QrA7m+AHDffff1+Dv8zDPPBB/T3fUV1Kcbb7xRLF26NPjnQCAgsrOzxapVqyRGZUwrV64UM2bMCPlYfX29sNvtYuPGjcH7jh49KgCIsrKyGEVoXADEm2++GfyzqqoiKytL/OxnPwveV19fL5xOp3jjjTeEEEIcOXJEABDvv/9+8Jj//d//FYqiiHPnzsUsdqPofY2FEOKee+4RX/nKV/p8Dq9xeGpqagQAsXPnTiHEwD4XtmzZIiwWi3C73cFjXnrpJZGSkiI8Hk9s34DO9b6+Qghx6623ioceeqjP5+jt+nKEpQ9erxf79u1DQUFB8D6LxYKCggKUlZVJjMy4Tpw4gezsbEycOBF33303KisrAQD79u2Dz+frca2nTJmCcePG8VoPwunTp+F2u3tcz9TUVOTl5QWvZ1lZGdLS0nDDDTcEjykoKIDFYsGePXtiHrNR7dixAxkZGbjmmmvw4IMP4uLFi8HHeI3D09DQAAAYMWIEgIF9LpSVlWH69OnIzMwMHlNYWIjGxkZ8+OGHMYxe/3pfX81vf/tbpKenY9q0aVi+fDlaW1uDj+nt+ppi88NoqK2tRSAQ6PEfCgAyMzPx0UcfSYrKuPLy8vDqq6/immuuQVVVFZ588kl89rOfxeHDh+F2u+FwOJCWltbjOZmZmXC73XICNjDtmoX6u6s95na7kZGR0eNxm82GESNG8JoP0Ny5c/G1r30NEyZMwKlTp/DYY4/hzjvvRFlZGaxWK69xGFRVxcMPP4zPfOYzmDZtGgAM6HPB7XaH/HuuPUYdQl1fAPjnf/5njB8/HtnZ2Th06BAeffRRHDt2DH/84x8B6O/6MmGhmLjzzjuD//+6665DXl4exo8fj9/97ndISEiQGBnR4Hzzm98M/v/p06fjuuuuw6RJk7Bjxw7cfvvtEiMznqVLl+Lw4cM96toocvq6vt3rqaZPn47Ro0fj9ttvx6lTpzBp0qRYh3lFnBLqQ3p6OqxW62UV6dXV1cjKypIUlXmkpaXh6quvxsmTJ5GVlQWv14v6+voex/BaD452zfr7u5uVlXVZ8bjf70ddXR2v+SBNnDgR6enpOHnyJABe44EqKirC22+/jXfeeQdjx44N3j+Qz4WsrKyQf8+1x6jv6xtKXl4eAPT4O6yn68uEpQ8OhwOzZ89GaWlp8D5VVVFaWor8/HyJkZlDc3MzTp06hdGjR2P27Nmw2+09rvWxY8dQWVnJaz0IEyZMQFZWVo/r2djYiD179gSvZ35+Purr67Fv377gMdu3b4eqqsEPLQrPJ598gosXL2L06NEAeI2vRAiBoqIivPnmm9i+fTsmTJjQ4/GBfC7k5+fjgw8+6JEYbt26FSkpKZg6dWps3ohOXen6hnLw4EEA6PF3WFfXN+Zlvgayfv164XQ6xauvviqOHDki7r//fpGWltajYpoG5pFHHhE7duwQp0+fFu+9954oKCgQ6enpoqamRgghxAMPPCDGjRsntm/fLvbu3Svy8/NFfn6+5Kj1q6mpSRw4cEAcOHBAABDPP/+8OHDggDhz5owQQoinn35apKWliT/96U/i0KFD4itf+YqYMGGCaGtrC55j7ty54vrrrxd79uwRf//738XkyZPFokWLZL0l3envGjc1NYnvf//7oqysTJw+fVps27ZNzJo1S0yePFm0t7cHz8Fr3LcHH3xQpKamih07doiqqqrgrbW1NXjMlT4X/H6/mDZtmrjjjjvEwYMHRUlJiRg1apRYvny5jLekK1e6vidPnhQ//vGPxd69e8Xp06fFn/70JzFx4kRxyy23BM+ht+vLhOUKfvnLX4px48YJh8MhbrzxRrF7927ZIRnSwoULxejRo4XD4RBjxowRCxcuFCdPngw+3tbWJr7zne+I4cOHi8TERPHVr35VVFVVSYxY39555x0B4LLbPffcI4ToWNr8ox/9SGRmZgqn0yluv/12cezYsR7nuHjxoli0aJEYNmyYSElJEUuWLBFNTU0S3o0+9XeNW1tbxR133CFGjRol7Ha7GD9+vLjvvvsu+zHDa9y3UNcWgPj1r38dPGYgnwsVFRXizjvvFAkJCSI9PV088sgjwufzxfjd6M+Vrm9lZaW45ZZbxIgRI4TT6RRXXXWVWLZsmWhoaOhxHj1dX0UIIWI3nkNEREQUPtawEBERke4xYSEiIiLdY8JCREREuseEhYiIiHSPCQsRERHpHhMWIiIi0j0mLERERKR7TFiIiIhI95iwEBERke4xYSEiIiLdY8JCREREuseEhYiIiHTv/wEBuQdeFlMMQgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.49990098e+01 3.49898358e+01 3.49599073e+01 3.49098768e+01\n",
      " 3.48420544e+01 3.47559647e+01 3.46481716e+01 3.45228280e+01\n",
      " 3.43770294e+01 3.42094074e+01 3.40247067e+01 3.38176291e+01\n",
      " 3.35900995e+01 3.33414894e+01 3.30715494e+01 3.27795947e+01\n",
      " 3.24665168e+01 3.21289285e+01 3.17733466e+01 3.13885187e+01\n",
      " 3.09836277e+01 3.05534429e+01 3.00984042e+01 2.96160807e+01\n",
      " 2.91092297e+01 2.85752058e+01 2.80117451e+01 2.74242599e+01\n",
      " 2.68033707e+01 2.61554202e+01 2.54726307e+01 2.47626374e+01\n",
      " 2.40161943e+01 2.32384652e+01 2.24229193e+01 2.15717943e+01\n",
      " 2.06813953e+01 1.97517383e+01 1.87818786e+01 1.77698262e+01\n",
      " 1.67143898e+01 1.56112436e+01 1.44597158e+01 1.32584108e+01\n",
      " 1.20054142e+01 1.06984422e+01 9.33324429e+00 7.91059180e+00\n",
      " 6.42447592e+00 4.87298040e+00 3.25290708e+00 1.56204736e+00\n",
      " 1.29983117e-05 1.31838488e-06 6.91402967e-07 4.67559183e-07\n",
      " 3.52722771e-07 2.82890614e-07 2.35968101e-07 2.02291855e-07\n",
      " 1.76957429e-07 1.57218712e-07 1.41416081e-07 1.28484286e-07\n",
      " 1.17715518e-07 1.08611142e-07 1.00818167e-07 9.40763880e-08\n",
      " 8.81898424e-08 8.30084283e-08 7.84152790e-08 7.43182559e-08\n",
      " 7.06428504e-08 6.73295907e-08 6.43291347e-08 6.16008376e-08\n",
      " 5.91109729e-08 5.68311429e-08 5.47371260e-08 5.28085956e-08\n",
      " 5.10277648e-08 4.93795581e-08 4.78509755e-08 4.64305241e-08\n",
      " 4.51082841e-08 4.38754231e-08 4.27241094e-08 4.16476098e-08\n",
      " 4.06395096e-08 3.96941847e-08 3.88068309e-08 3.79730376e-08\n",
      " 3.71889739e-08 3.64511158e-08 3.57563820e-08 3.51019635e-08\n",
      " 3.44852689e-08 3.39039708e-08 3.33559677e-08 3.28393165e-08\n",
      " 3.23522224e-08 3.18930705e-08 3.14604253e-08 3.10528815e-08\n",
      " 3.06692400e-08 3.03083445e-08 2.99691594e-08 2.96507151e-08\n",
      " 2.93521383e-08 2.90725898e-08 2.88113335e-08 2.85676491e-08\n",
      " 2.83408881e-08 2.81304407e-08 2.79357592e-08 2.77563639e-08\n",
      " 2.75918388e-08 2.74418159e-08 2.73061720e-08 2.71848234e-08\n",
      " 2.70772760e-08 2.69827382e-08 2.69003407e-08 2.68314700e-08\n",
      " 2.67758238e-08 2.67329386e-08 2.67024796e-08 2.66842577e-08\n",
      " 2.66781910e-08 2.66842577e-08 2.67024796e-08 2.67329386e-08\n",
      " 2.67758234e-08 2.68314700e-08 2.69003407e-08 2.69827382e-08\n",
      " 2.70772753e-08 2.71848234e-08 2.73061722e-08 2.74418163e-08\n",
      " 2.75918405e-08 2.77563645e-08 2.79357585e-08 2.81304402e-08\n",
      " 2.83408881e-08 2.85676496e-08 2.88113318e-08 2.90725926e-08\n",
      " 2.93521383e-08 2.96507151e-08 2.99691589e-08 3.03083445e-08\n",
      " 3.06692297e-08 3.10528815e-08 3.14604171e-08 3.18930715e-08\n",
      " 3.23522224e-08 3.28393165e-08 3.33560064e-08 3.39039813e-08\n",
      " 3.44852574e-08 3.51019548e-08 3.57563891e-08 3.64511323e-08\n",
      " 3.71889853e-08 3.79729916e-08 3.88068429e-08 3.96942255e-08\n",
      " 4.06395169e-08 4.16476098e-08 4.27241742e-08 4.38754231e-08\n",
      " 4.51082841e-08 4.64305323e-08 4.78509755e-08 4.93796049e-08\n",
      " 5.10277648e-08 5.28086076e-08 5.47371260e-08 5.68311429e-08\n",
      " 5.91109729e-08 6.16008376e-08 6.43291347e-08 6.73295907e-08\n",
      " 7.06428504e-08 7.43182559e-08 7.84152790e-08 8.30084283e-08\n",
      " 8.81898424e-08 9.40763880e-08 1.00818167e-07 1.08611142e-07\n",
      " 1.17715518e-07 1.28484286e-07 1.41416081e-07 1.57218712e-07\n",
      " 1.76957429e-07 2.02291855e-07 2.35968101e-07 2.82890614e-07\n",
      " 3.52722771e-07 4.67559183e-07 6.91410567e-07 1.31838488e-06\n",
      " 1.30042372e-05 1.56204736e+00 3.25290708e+00 4.87298040e+00\n",
      " 6.42447592e+00 7.91049007e+00 9.33324429e+00 1.06962698e+01\n",
      " 1.20054142e+01 1.32584108e+01 1.44590878e+01 1.56112436e+01\n",
      " 1.67143898e+01 1.77690554e+01 1.87825929e+01 1.97517383e+01\n",
      " 2.06813953e+01 2.15715509e+01 2.24229193e+01 2.32384652e+01\n",
      " 2.40161943e+01 2.47635320e+01 2.54757849e+01 2.61554202e+01\n",
      " 2.68033707e+01 2.74237791e+01 2.80118690e+01 2.85752740e+01\n",
      " 2.91092297e+01 2.96160807e+01 3.00972422e+01 3.05534429e+01\n",
      " 3.09836277e+01 3.13885187e+01 3.17733466e+01 3.21289285e+01\n",
      " 3.24665168e+01 3.27795947e+01 3.30715494e+01 3.33414894e+01\n",
      " 3.35896537e+01 3.38176291e+01 3.40247067e+01 3.42094074e+01\n",
      " 3.43770294e+01 3.45228280e+01 3.46481716e+01 3.47559647e+01\n",
      " 3.48420544e+01 3.49098768e+01 3.49598119e+01 3.49887619e+01]\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "NFFT = 256\n",
    "\n",
    "# Channel\n",
    "#h = [1.0 + 0.5j, 0.5 - 0.3j, -0.2 + 0.3j]\n",
    "h = [1, 1]\n",
    "\n",
    "# Subcarrier gains\n",
    "H = np.fft.fft(h, NFFT) / np.sqrt(NFFT)\n",
    "\n",
    "# Per symbol power constraint\n",
    "P_total = 10\n",
    "P = cp.Variable(NFFT)\n",
    "objective = cp.Maximize(cp.sum(cp.log(1 + cp.multiply(cp.abs(H)**2,  P))))\n",
    "constraints = [P >= 0, cp.sum(P) <= NFFT * P_total]\n",
    "prob = cp.Problem(objective, constraints)\n",
    "result = prob.solve()\n",
    "\n",
    "plt.stem(np.abs(H))\n",
    "plt.plot(P.value / np.max(P.value) * np.max(abs(H)), 'k')\n",
    "plt.show()\n",
    "print(P.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf502e03-525a-462e-a918-1b43f63f8f65",
   "metadata": {},
   "source": [
    "## Ergodic capacity with CSIT\n",
    "We consider the SISO case. $h$ varies rapidly allowing us to allocate power effectively based on its current value. Let us first find $(P, \\lambda)$ pairs using\n",
    "\\begin{equation}\n",
    "{\\bf E}\\left[\\left(\\frac{1}{\\lambda} - \\frac{N_0}{|h|^2}\\right)^+\\right] = P\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "28727eca-10ca-4248-bb26-b2c50ab2ce45",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.2668008383164935\n"
     ]
    }
   ],
   "source": [
    "N_ITER = 10000\n",
    "LAMBDA = 0.1\n",
    "N0 = 1\n",
    "P = 0\n",
    "for i in range(N_ITER):\n",
    "    h = (np.random.randn() + np.random.randn() * 1j) / np.sqrt(2)\n",
    "    P_alloc = 1 / LAMBDA - N0 / np.abs(h)**2\n",
    "    if P_alloc > 0:\n",
    "        P = P + P_alloc\n",
    "P = P / N_ITER\n",
    "print(P)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}