From e4ad851e0d1a7c191d8ae98375c01b5e55ad851f Mon Sep 17 00:00:00 2001 From: Joshua Hoskins Date: Fri, 20 Dec 2024 12:14:51 -0500 Subject: [PATCH] Add encoder and antenna_map to from_visibility() function. Still need to work out how the prototype of the Jones matrix should work. Soon(TM) --- jones.ipynb | 269 ++++++++++++++++++++++++------------- src/calviper/jones.py | 29 ++-- src/calviper/math/tools.py | 2 +- testing.ipynb | 48 +++++-- 4 files changed, 234 insertions(+), 114 deletions(-) diff --git a/jones.ipynb b/jones.ipynb index 1c0e2b2..f5b2047 100644 --- a/jones.ipynb +++ b/jones.ipynb @@ -15,6 +15,7 @@ "import numpy as np\n", "import pandas as pd\n", "import calviper as cv\n", + "import matplotlib.pyplot as plt\n", "\n", "from xradio import measurement_set as ms" ] @@ -198,31 +199,40 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "id": "79e116f8-09d8-4151-b30a-ada716b467f3", "metadata": {}, "outputs": [], "source": [ - "V = sps.VISIBILITY.mean(dim=\"time\").data.compute()" + "V = dataset.VISIBILITY.mean(dim=\"time\").data.compute()" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "id": "4dc7448f-20ab-4eca-9e8c-11eba2e971ca", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "(957, 45, 8, 4)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "s = sps.VISIBILITY.shape\n", + "s = dataset.VISIBILITY.shape\n", "s" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "0aa8f93c-b00b-402f-9574-fe3f80a4c054", + "cell_type": "raw", + "id": "9caf29e9-a44f-412a-a260-dd61ff5533fa", "metadata": {}, - "outputs": [], "source": [ "def from_vis(dataset):\n", " shape = dataset.VISIBILITY.shape\n", @@ -239,7 +249,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "id": "0ce2493e-0c7e-4e6f-9711-5b63df1973ef", "metadata": {}, "outputs": [], @@ -259,17 +269,17 @@ }, { "cell_type": "code", - "execution_count": 10, - "id": "b590851a-7915-48d4-a3ab-9d38b8bdb284", + "execution_count": 11, + "id": "6403894b-a634-45a8-9219-05aed71457f1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "(957, 45, 8, 4, 9, 9)" + "(957, 9, 9)" ] }, - "execution_count": 10, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } @@ -280,172 +290,249 @@ }, { "cell_type": "code", - "execution_count": 66, - "id": "c942215e-72b3-4085-b1e8-66809ec6796b", + "execution_count": 24, + "id": "3ee0983f-5a19-4001-973c-c88216caa061", "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", - "\n", - "test = np.array([\n", - " [\n", - " [\n", - " [1, 0], \n", - " [0, 1]\n", - " ],\n", - " [\n", - " [2, 0], \n", - " [0, 2]\n", - " ]\n", - " ],\n", - " [\n", - " [\n", - " [3, 0], \n", - " [0, 3]\n", - " ],\n", - " [\n", - " [4, 0], \n", - " [0, 4]\n", - " ]\n", - " ]\n", - "])" + "G.matrix = G.matrix.mean(axis=0)" ] }, { "cell_type": "code", - "execution_count": 67, - "id": "3f5af5ee-1423-4f9a-a72f-b349d0941043", + "execution_count": 13, + "id": "bf6e3883-7923-4a5f-be89-ac0244ef8c17", "metadata": {}, "outputs": [], "source": [ - "def diagonal(x):\n", - " n_axis = len(x.shape)\n", - " axis1 = n_axis - 2\n", - " axis2 = n_axis - 1\n", - " \n", - " return test.diagonal(axis1=axis1, axis2=axis2)" + "# This should all be done within teh equivelent of VisEquation\n", + "v = V[:, 0, 0]\n", + "\n", + "index_a, ant = cv.math.tools.encode(dataset.baseline_antenna1_name.to_numpy())\n", + "index_b, _ = cv.math.tools.encode(dataset.baseline_antenna2_name.to_numpy())\n", + "\n", + "V = cv.math.tools.build_visibility_matrix(array=v, index_a=index_a, index_b=index_b)" ] }, { "cell_type": "code", - "execution_count": 68, - "id": "97ed10af-1230-4fe7-a0a4-d8bc1c182178", + "execution_count": 16, + "id": "467d3317-250f-4a98-8056-6615245700c3", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[[1, 1],\n", - " [2, 2]],\n", - "\n", - " [[3, 3],\n", - " [4, 4]]])" + "(9, 9)" ] }, - "execution_count": 68, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "diagonal(test)" + "V.shape" ] }, { "cell_type": "code", - "execution_count": 35, - "id": "f390f259-024f-462f-80b0-47183d770ca1", + "execution_count": 17, + "id": "2ff5cc79-02c8-49a0-9016-c7113a657226", "metadata": {}, "outputs": [], "source": [ - "param = np.array([1, 1, 2, 2, 3, 3, 4, 4])" + "solver = cv.math.solver.least_squares.LeastSquaresSolver()" ] }, { "cell_type": "code", - "execution_count": 64, - "id": "7a1277f3-611c-4e3e-88e0-f03c66f5d311", + "execution_count": 18, + "id": "99704e1f-e800-40bf-a0e4-5d4715e011b3", "metadata": {}, "outputs": [], "source": [ - "param.reshape((4, 2, 2))" + "gain_solutions = solver.solve(\n", + " vis=V,\n", + " iterations=40,\n", + " optimizer=cv.math.optimizer.MeanSquaredError(alpha=0.2),\n", + " stopping=1e-4\n", + ")" ] }, { "cell_type": "code", - "execution_count": 41, - "id": "4559a5c3-3734-47d0-8312-d8103a27bfb4", + "execution_count": 19, + "id": "4dba3ecf-0693-4462-9956-3a4e9541ba77", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKO1JREFUeJzt3XFs1HWe//HXF5BWpTMeYDttKF4Xb9FuLdqqONFtXCkCZzhZMdlVCWiIzTbFiN2985rs6XK7l+KauLi3WMXcyV1qZeNmUfFiiVtCu2aLQruNdPnZLKSG7o9Ocd0wU7rXgXTm90d/MzilLZ12+p1PP30+km9iv/Ptd96kf8zL7/f7eY0TjUajAgAAcMmcdA8AAABmF8IHAABwFeEDAAC4ivABAABcRfgAAACuInwAAABXET4AAICrCB8AAMBV89I9wEiRSERnzpxRVlaWHMdJ9zgAAGACotGo+vv7lZeXpzlzxr+2YVz4OHPmjPLz89M9BgAAmISenh4tWbJk3GOMCx9ZWVmShof3eDxpngYAAExEKBRSfn5+/HN8PMaFj9itFo/HQ/gAAGCGmcgjEzxwCgAAXEX4AAAAriJ8AAAAVxE+AACAqwgfAADAVYQPAADgKsIHAABwFeEDAAC4yriSsekyFInqk+6/6Gz/oLKzMnVnwULNncN3xwAA4LYpXfnYuXOnHMfR9u3b4/sGBwdVVVWlRYsWacGCBdq4caP6+vqmOueUNHb26p4XDumR14/o6X0deuT1I7rnhUNq7OxN61wAAMxGkw4fR48e1Wuvvabi4uKE/c8884wOHDigt99+W83NzTpz5oweeuihKQ86WY2dvaqsb1dvcDBhfyA4qMr6dgIIAAAum1T4OH/+vB577DG9/vrr+pu/+Zv4/mAwqP/4j//QSy+9pPvuu0+lpaV644039Lvf/U5HjhxJ2dATNRSJaseBE4qO8lps344DJzQUGe0IAAAwHSYVPqqqqvTAAw+ovLw8YX9bW5suXryYsP+mm27S0qVL1draOuq5wuGwQqFQwpYqn3T/5bIrHl8VldQbHNQn3X9J2XsCAIDxJf3A6b59+9Te3q6jR49e9logEND8+fN13XXXJezPyclRIBAY9Xy1tbXasWNHsmNMyNn+sYPHZI4DAABTl9SVj56eHj399NN68803lZmZmZIBampqFAwG41tPT09KzitJ2VkTm3GixwEAgKlLKny0tbXp7NmzKikp0bx58zRv3jw1Nzfr5z//uebNm6ecnBxduHBB586dS/i9vr4++Xy+Uc+ZkZEhj8eTsKXKnQULlevN1FgLah1Jud7hZbcAAMAdSYWPVatW6fjx4+ro6Ihvt99+ux577LH4f1911VVqamqK/05XV5dOnz4tv9+f8uGvZO4cR8+vL5SkywJI7Ofn1xfS9wEAgIuSeuYjKytLRUVFCfuuvfZaLVq0KL5/69atqq6u1sKFC+XxePTUU0/J7/frrrvuSt3USVhblKu6TSX60XsnFAhderbD583U8+sLtbYoNy1zAQAwW6W84fRnP/uZ5syZo40bNyocDmvNmjV65ZVXUv02k5C4nDYaZXktAADp4EQN+xQOhULyer0KBoMpef4jVjI28h8Zu9FSt6mEqx8AAExRMp/fVn+xHCVjAACYx+rwQckYAADmsTp8UDIGAIB5rA4flIwBAGAeq8MHJWMAAJjH6vBByRgAAOaxOnxIl0rGcjyJt1Z83kyW2QIAkAbWh49LKBkDAMAE1oePWMlYIBRO2N8XCquyvl2Nnb1pmgwAgNnJ6vBByRgAAOaxOnxQMgYAgHmsDh+UjAEAYB6rwwclYwAAmMfq8EHJGAAA5rE6fFAyBgCAeawOH9JwyVhFWYGcEfnCcaSKsgJKxgAAcJn14aOxs1d7Wro1cjVtJCrtaemm5wMAAJdZHT7G6/mIoecDAAB3WR0+6PkAAMA8VocPej4AADCP1eGDng8AAMxjdfig5wMAAPNYHT7o+QAAwDxWhw9puOejblOJcjyJt1Z83kzVbSqh5wMAAJdZHz4uSVxOG42yvBYAgHSwPnw0dvaqsr5dgVA4YX9fKKzK+nZKxgAAcJnV4WO8krHYPkrGAABwl9Xhg5IxAADMY3X4oGQMAADzWB0+KBkDAMA8VocPSsYAADCP1eGDkjEAAMyTVPioq6tTcXGxPB6PPB6P/H6/Pvjgg/jr9957rxzHSdi+973vpXzoZFAyBgCAWeYlc/CSJUu0c+dO/d3f/Z2i0aj+67/+Sw8++KB+//vf6xvf+IYk6cknn9S//uu/xn/nmmuuSe3Ek0bJGAAAJkgqfKxfvz7h53/7t39TXV2djhw5Eg8f11xzjXw+X+omnKJYydjIqBErGePqBwAA7pr0Mx9DQ0Pat2+fBgYG5Pf74/vffPNNLV68WEVFRaqpqdFf//rXlAw6qRkpGQMAwDhJXfmQpOPHj8vv92twcFALFizQ/v37VVg4/FDno48+qhtuuEF5eXn69NNP9eyzz6qrq0u//vWvxzxfOBxWOHyp+jwUCk3inzG6ZErG/MsWpex9AQDA2JIOH8uXL1dHR4eCwaB+9atfacuWLWpublZhYaEqKirix91yyy3Kzc3VqlWrdOrUKS1btmzU89XW1mrHjh2T/xeMg5IxAADMk/Rtl/nz5+vGG29UaWmpamtrtWLFCr388sujHrty5UpJ0smTJ8c8X01NjYLBYHzr6elJdqQxUTIGAIB5kr7yMVIkEkm4bfJVHR0dkqTc3LEf6MzIyFBGRsZUxxhVrGQsEBwc9bkPR8NLbikZAwDAPUmFj5qaGq1bt05Lly5Vf3+/GhoadPjwYR08eFCnTp1SQ0OD/v7v/16LFi3Sp59+qmeeeUZlZWUqLi6ervnHFSsZq6xvl6PExbaUjAEAkB5JhY+zZ89q8+bN6u3tldfrVXFxsQ4ePKjVq1erp6dHv/nNb7Rr1y4NDAwoPz9fGzdu1A9/+MPpmn1C1hblqqKsQK//tltfrfZwHOnJbxawzBYAAJc5UcPatkKhkLxer4LBoDwez5TPN1bPhzR89YOeDwAApi6Zz2+rv9tlvJ6PGHo+AABwl9XhI5meDwAA4A6rwwc9HwAAmMfq8EHPBwAA5rE6fMR6PsZaSOtIyqXnAwAAV1kdPmI9H5IuCyD0fAAAkB5Whw9puOejblOJcjyJt1Z83kyW2QIAkAbWh49LEpfTGlZvAgDArGF9+IiVjAVCid8/0xcKq7K+XY2dvWmaDACA2cnq8DFeyVhsHyVjAAC4y+rwQckYAADmsTp8UDIGAIB5rA4flIwBAGAeq8MHJWMAAJjH6vBByRgAAOaxOnxIwyVjFWUFckbkC8eRKsoKKBkDAMBl1oePxs5e7Wnp1sjVtJGotKelm54PAABcZnX4GK/nI4aeDwAA3GV1+KDnAwAA81gdPuj5AADAPFaHD3o+AAAwj9Xhg54PAADMY3X4oOcDAADzWB0+pOGej7pNJcrxJN5a8XkzVbephJ4PAABcZn34uCRxOW00yvJaAADSwfrw0djZq8r6dgVC4YT9faGwKuvbKRkDAMBlVoeP8UrGYvsoGQMAwF1Whw9KxgAAMI/V4YOSMQAAzGN1+KBkDAAA81gdPigZAwDAPFaHD0rGAAAwT1Lho66uTsXFxfJ4PPJ4PPL7/frggw/irw8ODqqqqkqLFi3SggULtHHjRvX19aV86GRQMgYAgFmSCh9LlizRzp071dbWpmPHjum+++7Tgw8+qD/84Q+SpGeeeUYHDhzQ22+/rebmZp05c0YPPfTQtAyePErGAAAwgROd4qfwwoUL9eKLL+rhhx/W9ddfr4aGBj388MOSpM8++0w333yzWltbddddd03ofKFQSF6vV8FgUB6PZyqjSbpUMjbyHxm70cLVDwAApi6Zz+9JP/MxNDSkffv2aWBgQH6/X21tbbp48aLKy8vjx9x0001aunSpWltbJ/s2U0LJGAAA5pmX7C8cP35cfr9fg4ODWrBggfbv36/CwkJ1dHRo/vz5uu666xKOz8nJUSAQGPN84XBY4fCl6vNQKJTsSGNKpmTMv2xRyt4XAACMLekrH8uXL1dHR4c+/vhjVVZWasuWLTpx4sSkB6itrZXX641v+fn5kz7XSJSMAQBgnqTDx/z583XjjTeqtLRUtbW1WrFihV5++WX5fD5duHBB586dSzi+r69PPp9vzPPV1NQoGAzGt56enqT/EWOhZAwAAPNMuecjEokoHA6rtLRUV111lZqamuKvdXV16fTp0/L7/WP+fkZGRnzpbmxLFUrGAAAwT1LPfNTU1GjdunVaunSp+vv71dDQoMOHD+vgwYPyer3aunWrqqurtXDhQnk8Hj311FPy+/0TXumSarGSscr6djlKXGxLyRgAAOmRVPg4e/asNm/erN7eXnm9XhUXF+vgwYNavXq1JOlnP/uZ5syZo40bNyocDmvNmjV65ZVXpmXwiVpblKuKsgK9/ttufXVRseNIT36zgGW2AAC4bMo9H6nmVs+HNHz1g54PAACmzpWej5lgvJ6PGHo+AABwl9XhI5meDwAA4A6rwwc9HwAAmMfq8EHPBwAA5rE6fNDzAQCAeawOH7GeD0mXBRB6PgAASA+rw4c03PNRt6lEOZ7EWys+bybLbAEASAPrw8clictpDas3AQBg1rA+fMRKxgKhcML+vlBYlfXtauzsTdNkAADMTlaHj/FKxmL7KBkDAMBdVocPSsYAADCP1eGDkjEAAMxjdfigZAwAAPNYHT4oGQMAwDxWhw9KxgAAMI/V4UMaLhmrKCuQMyJfOI5UUVZAyRgAAC6zPnw0dvZqT0u3Rq6mjUSlPS3d9HwAAOAyq8PHeD0fMfR8AADgLqvDBz0fAACYx+rwQc8HAADmsTp80PMBAIB5rA4f9HwAAGAeq8MHPR8AAJjH6vAhDfd81G0qUY4n8daKz5upuk0l9HwAAOAy68PHJYnLaaNRltcCAJAO1oePxs5eVda3KxAKJ+zvC4VVWd9OyRgAAC6zOnyMVzIW20fJGAAA7rI6fFAyBgCAeawOH5SMAQBgHqvDByVjAACYx+rwQckYAADmsTp8UDIGAIB5rA4fEiVjAACYJqnwUVtbqzvuuENZWVnKzs7Whg0b1NXVlXDMvffeK8dxErbvfe97KR16cigZAwDABEmFj+bmZlVVVenIkSP68MMPdfHiRd1///0aGBhIOO7JJ59Ub29vfPvpT3+a0qGTQckYAABmmZfMwY2NjQk/7927V9nZ2Wpra1NZWVl8/zXXXCOfz5eaCafgSiVjjoZLxlYX+njuAwAAl0zpmY9gMChJWrgwcbXIm2++qcWLF6uoqEg1NTX661//OuY5wuGwQqFQwpYqlIwBAGCepK58fFUkEtH27dt19913q6ioKL7/0Ucf1Q033KC8vDx9+umnevbZZ9XV1aVf//rXo56ntrZWO3bsmOwY46JkDAAA80w6fFRVVamzs1MfffRRwv6Kior4f99yyy3Kzc3VqlWrdOrUKS1btuyy89TU1Ki6ujr+cygUUn5+/mTHSkDJGAAA5plU+Ni2bZvef/99tbS0aMmSJeMeu3LlSknSyZMnRw0fGRkZysjImMwYVxQrGQsEB0d97sPR8JJbSsYAAHBPUs98RKNRbdu2Tfv379ehQ4dUUFBwxd/p6OiQJOXmut+nQckYAADmSSp8VFVVqb6+Xg0NDcrKylIgEFAgEND//u//SpJOnTqlH//4x2pra9Pnn3+u9957T5s3b1ZZWZmKi4un5R9wJWuLclVRViBnRL5wHKmirICSMQAAXOZEk2jbckZ+gv9/b7zxhh5//HH19PRo06ZN6uzs1MDAgPLz8/Xtb39bP/zhD+XxeCb0HqFQSF6vV8FgcMK/M55Yz8dYt11oOQUAYOqS+fxO6pmPK+WU/Px8NTc3J3PKaTVez0cMPR8AALjL6u92oecDAADzWB0+6PkAAMA8VocPej4AADCP1eEj1vMx1tMcjqRcej4AAHCV1eGDng8AAMxjdfiQhns+6jaVKMeTeGvF581kmS0AAGlgffi4JHHBbRL1JgAAIIWsDx+xkrFAKJywvy8UVmV9uxo7e9M0GQAAs5PV4WO8krHYvh0HTmgowlUQAADcYnX4oGQMAADzWB0+KBkDAMA8VocPSsYAADCP1eGDkjEAAMxjdfigZAwAAPNYHT4kSsYAADCN9eHjEkrGAAAwgfXhg5IxAADMYnX4oGQMAADzWB0+KBkDAMA8VocPSsYAADCP1eGDkjEAAMxjdfigZAwAAPNYHT4oGQMAwDxWhw9puGSsoqxAzoh84ThSRVkBJWMAALjM+vDR2NmrPS3dGrmaNhKV9rR00/MBAIDLrA4f4/V8xNDzAQCAu6wOH/R8AABgHqvDBz0fAACYx+rwQc8HAADmsTp80PMBAIB5rA4f9HwAAGAeq8OHNNzzUbepRDmexFsrPm+m6jaV0PMBAIDLkgoftbW1uuOOO5SVlaXs7Gxt2LBBXV1dCccMDg6qqqpKixYt0oIFC7Rx40b19fWldOjJSVxOG42yvBYAgHRIKnw0NzerqqpKR44c0YcffqiLFy/q/vvv18DAQPyYZ555RgcOHNDbb7+t5uZmnTlzRg899FDKB5+oxs5eVda3KxAKJ+zvC4VVWd9OyRgAAC5zolO4BPDFF18oOztbzc3NKisrUzAY1PXXX6+GhgY9/PDDkqTPPvtMN998s1pbW3XXXXdd8ZyhUEher1fBYFAej2eyo0kaLhm754VDY3Z9OBq+/fLRs/fx3AcAAFOQzOf3lJ75CAaDkqSFC4dXi7S1tenixYsqLy+PH3PTTTdp6dKlam1tHfUc4XBYoVAoYUsVSsYAADDPpMNHJBLR9u3bdffdd6uoqEiSFAgENH/+fF133XUJx+bk5CgQCIx6ntraWnm93viWn58/2ZEuQ8kYAADmmXT4qKqqUmdnp/bt2zelAWpqahQMBuNbT0/PlM73VZSMAQBgnnmT+aVt27bp/fffV0tLi5YsWRLf7/P5dOHCBZ07dy7h6kdfX598Pt+o58rIyFBGRsZkxriiWMlYIDg46pfLxZ75oGQMAAD3JHXlIxqNatu2bdq/f78OHTqkgoKChNdLS0t11VVXqampKb6vq6tLp0+flt/vT83ESaBkDAAA8yR15aOqqkoNDQ169913lZWVFX+Ow+v16uqrr5bX69XWrVtVXV2thQsXyuPx6KmnnpLf75/QSpfpsLYoVxVlBXr9t9366roex5Ge/GYBJWMAALgsqaW2jjP6FYI33nhDjz/+uKThkrHvf//7euuttxQOh7VmzRq98sorY952GSmVS22lSz0fY912oeUUAICpS+bze0o9H9OBng8AAGYe13o+TEfPBwAA5rE6fNDzAQCAeawOH/R8AABgHqvDR6znY6ynORxJufR8AADgKqvDBz0fAACYx+rwIQ33fNRtKlGOJ/HWis+byTJbAADSwPrwcUniimLDVhgDADBrWB8+YiVjgVA4YX9fKKzK+nY1dvamaTIAAGYnq8PHUCSqHQdOjNpuGtu348AJDUW4CgIAgFusDh+UjAEAYB6rwwclYwAAmMfq8EHJGAAA5rE6fFAyBgCAeawOH5SMAQBgHqvDh0TJGAAAprE+fFxCyRgAACawPnxQMgYAgFmsDh+UjAEAYB6rwwclYwAAmMfq8EHJGAAA5rE6fFAyBgCAeawOH5SMAQBgHqvDByVjAACYx+rwIQ2XjFWUFcgZkS8cR6ooK6BkDAAAl1kfPho7e7WnpVsjV9NGotKelm56PgAAcJnV4WO8no8Yej4AAHCX1eGDng8AAMxjdfig5wMAAPNYHT7o+QAAwDxWhw96PgAAMI/V4YOeDwAAzGN1+JCGez7qNpUox5N4a8XnzVTdphJ6PgAAcFnS4aOlpUXr169XXl6eHMfRO++8k/D6448/LsdxEra1a9emat4pSFxOG42yvBYAgHRIOnwMDAxoxYoV2r1795jHrF27Vr29vfHtrbfemtKQU9HY2avK+nYFQuGE/X2hsCrr2ykZAwDAZfOS/YV169Zp3bp14x6TkZEhn8836aFSZbySsaiGn/vYceCEVhf6eO4DAACXTMszH4cPH1Z2draWL1+uyspKffnll2MeGw6HFQqFErZUoWQMAADzpDx8rF27Vv/93/+tpqYmvfDCC2pubta6des0NDQ06vG1tbXyer3xLT8/P2WzUDIGAIB5kr7tciXf/e534/99yy23qLi4WMuWLdPhw4e1atWqy46vqalRdXV1/OdQKJSyAELJGAAA5pn2pbZf+9rXtHjxYp08eXLU1zMyMuTxeBK2VKFkDAAA80x7+PjTn/6kL7/8Urm57vdpUDIGAIB5kg4f58+fV0dHhzo6OiRJ3d3d6ujo0OnTp3X+/Hn94z/+o44cOaLPP/9cTU1NevDBB3XjjTdqzZo1qZ59QtYW5aqirEDOiHzhOFJFWQElYwAAuCzp8HHs2DHddtttuu222yRJ1dXVuu222/Tcc89p7ty5+vTTT/UP//AP+vrXv66tW7eqtLRUv/3tb5WRkZHy4SeisbNXe1q6FRmx3jYSlfa0dNPzAQCAy5yoYVWfoVBIXq9XwWBwys9/DEWiuueFQ2Mut3U0XLP+0bP3cesFAIApSObz2+rvdqHnAwAA81gdPuj5AADAPFaHD3o+AAAwj9Xhg54PAADMY3X4oOcDAADzWB0+pOGej7pNJcrxJN5a8XkzVbephJ4PAABcZn34uCRxRbFhK4wBAJg1rA8fjZ29qqxvVyAUTtjfFwqrsr6dkjEAAFxmdfgYikS148AJjXaNI7Zvx4ETGhpZfwoAAKaN1eGDkjEAAMxjdfigZAwAAPNYHT4oGQMAwDxWhw9KxgAAMI/V4YOSMQAAzGN1+JAoGQMAwDTWh49LKBkDAMAE1ocPSsYAADCL1eGDkjEAAMxjdfigZAwAAPNYHT4oGQMAwDxWhw9KxgAAMI/V4YOSMQAAzGN1+KBkDAAA81gdPqThkrGKsgI5I/KF40gVZQWUjAEA4DLrw0djZ6/2tHRr5GraSFTa09JNzwcAAC6zOnyM1/MRQ88HAADusjp80PMBAIB5rA4f9HwAAGAeq8MHPR8AAJjH6vBBzwcAAOaxOnzQ8wEAgHmsDh/ScM9H3aYS5XgSb634vJmq21RCzwcAAC5LOny0tLRo/fr1ysvLk+M4eueddxJej0ajeu6555Sbm6urr75a5eXl+uMf/5iqeacgcTltNMryWgAA0iHp8DEwMKAVK1Zo9+7do77+05/+VD//+c/16quv6uOPP9a1116rNWvWaHAwPStKGjt7VVnfrkAonLC/LxRWZX07JWMAALjMiU7hEoDjONq/f782bNggafhqQl5enr7//e/rBz/4gSQpGAwqJydHe/fu1Xe/+90rnjMUCsnr9SoYDMrj8Ux2NEnDJWP3vHBozK4PR8O3Xz569j6e+wAAYAqS+fxO6TMf3d3dCgQCKi8vj+/zer1auXKlWltbR/2dcDisUCiUsKUKJWMAAJgnpeEjEAhIknJychL25+TkxF8bqba2Vl6vN77l5+enbB5KxgAAME/aV7vU1NQoGAzGt56enpSdm5IxAADMk9Lw4fP5JEl9fX0J+/v6+uKvjZSRkSGPx5OwpQolYwAAmCel4aOgoEA+n09NTU3xfaFQSB9//LH8fn8q32pCKBkDAMA8SYeP8+fPq6OjQx0dHZKGHzLt6OjQ6dOn5TiOtm/frp/85Cd67733dPz4cW3evFl5eXnxFTFuo2QMAACzzEv2F44dO6Zvfetb8Z+rq6slSVu2bNHevXv1T//0TxoYGFBFRYXOnTune+65R42NjcrMTPdzFZSMAQBggin1fEyHVPZ8SJdKxkb+I2M3Wrj6AQDA1KWt58M0Q5Godhw4cVnwkC5dB9lx4ISGIkblLwAArGZ1+KBkDAAA81gdPigZAwDAPFaHD0rGAAAwj9Xhg5IxAADMY3X4oGQMAADzWB0+pOGSsYqyAjkj8oXjSBVlBSyzBQDAZdaHj8bOXu1p6dbI1bSRqLSnpVuNnb3pGQwAgFnK6vAxXs9HDD0fAAC4y+rwQc8HAADmsTp80PMBAIB5rA4f9HwAAGAeq8MHPR8AAJjH6vBBzwcAAOaxOnxIwz0fdZtK5PMm3lpZeO187X60hJ4PAABcZn34kIYDyL88cLMWXntVfN+XAxf04/85Qc8HAAAumxXho7GzV1UNv9dfBi4m7A8EB1VZ304AAQDARdaHj/GKxmL7KBoDAMA91ocPisYAADCL9eGDojEAAMxiffigaAwAALNYHz4oGgMAwCzWh4/xisak4Wc+/uUBisYAAHCL9eFDGrtoLIa+DwAA3DMrwod0qWhsNPR9AADgnlkTPoYiUf34f/7PqK/R9wEAgHtmTfig7wMAADPMmvBB3wcAAGaYNeGDvg8AAMwwa8IHfR8AAJhh1oSP8fo+Yj8/v56+DwAAptusCR/S2H0fPm+m6jaVaHWhT62nvtS7Hf9Xrae+ZOULAADTYF6qT/ijH/1IO3bsSNi3fPlyffbZZ6l+q0lZW5Sr1YU+fdL9F53tH1R21vCtlg9PBHTPC4cSVsTkejP1/PpCrS3KTePEAADYJeXhQ5K+8Y1v6De/+c2lN5k3LW8zaXPnOPIvWxT/ubGzV5X17Rp5nSNWPla3qYQAAgBAikxLKpg3b558Pt90nDrlhiJR7Thw4rLgIQ13fzgaLh9bXejjeRAAAFJgWp75+OMf/6i8vDx97Wtf02OPPabTp0+PeWw4HFYoFErY3ET5GAAA7kp5+Fi5cqX27t2rxsZG1dXVqbu7W9/85jfV398/6vG1tbXyer3xLT8/P9UjjYvyMQAA3OVEo9FpXdJx7tw53XDDDXrppZe0devWy14Ph8MKh8Pxn0OhkPLz8xUMBuXxeKZzNElS66kv9cjrR6543FtP3pXwnAgAALgkFArJ6/VO6PN72p8Eve666/T1r39dJ0+eHPX1jIwMZWRkTPcYY4qVjwWCg6M+9+FoeCku5WMAAKTGtPd8nD9/XqdOnVJurpmrRSgfAwDAXSkPHz/4wQ/U3Nyszz//XL/73e/07W9/W3PnztUjjzyS6rdKmSuVj010me1QJEpJGQAAV5Dy2y5/+tOf9Mgjj+jLL7/U9ddfr3vuuUdHjhzR9ddfn+q3SqmxyscmesWjsbNXOw6coKQMAIArmPYHTpOVzAMrphirpCwWWygpAwDYLpnP71n13S7T4UolZdJwSRm3YAAAGEb4mCJKygAASA7hY4ooKQMAIDmEjynKzsq88kFJHAcAgO0IH1MUKykba02Mo+FVL5SUAQAwjPAxRTOxpIw+EgBAOk17vfpsECspG9nz4TOw54M+EgBAutHzkUJDkeikS8rcQB8JAGC6GPXFcrPJ3DmOsd98e6U+EkfDfSSrC31GBSYAgH145mOWoI8EAGAKwscsQR8JAMAUhI9Zgj4SAIApCB+zBH0kAABTED5miZnYRwIAsBPhYxaJ9ZH4vIm3VnzeTJbZGo5iOAA2YantLLO2KFerC31G95EgEcVwAGxDyRhgMIrhAMwUyXx+c9sFMNSViuGk4WI4bsEAmGkIH4ChKIYDYCvCB2AoiuEA2IrwARiKYjgAtiJ8AIaiGA6ArQgfgKEohgOQaqZ0BtHzARgsVgw3sufDR88HgCSZ1BlEzwcwAwxFohTDAZg0NzqDkvn85soHMAPMnePIv2xRuscAMANdqTPI0XBn0OpCn2v/U8MzHwAAWMzEziDCBwAAFjOxM4jwAQCAxUzsDCJ8AABgMRM7gwgfAABYzMTOIMIHAACWi3UG+byJt1Z83syULLNN1rQttd29e7defPFFBQIBrVixQv/+7/+uO++8c7reDgAAjGNtUa5WF/qM6AyalvDxy1/+UtXV1Xr11Ve1cuVK7dq1S2vWrFFXV5eys7On4y0BAMAVmNIZNC23XV566SU9+eSTeuKJJ1RYWKhXX31V11xzjf7zP/9zOt4OAADMICkPHxcuXFBbW5vKy8svvcmcOSovL1dra+tlx4fDYYVCoYQNAADYK+Xh489//rOGhoaUk5OTsD8nJ0eBQOCy42tra+X1euNbfn5+qkcCAAAGSftql5qaGgWDwfjW09OT7pEAAMA0SvkDp4sXL9bcuXPV19eXsL+vr08+n++y4zMyMpSRkZHqMQAAgKFSfuVj/vz5Ki0tVVNTU3xfJBJRU1OT/H5/qt8OAADMMNOy1La6ulpbtmzR7bffrjvvvFO7du3SwMCAnnjiiel4OwAAMINMS/j4zne+oy+++ELPPfecAoGAbr31VjU2Nl72ECoAAJh9nGg0Gk33EF8VCoXk9XoVDAbl8XjSPQ4AAJiAZD6/p61efbJiWYi+DwAAZo7Y5/ZErmkYFz76+/slib4PAABmoP7+fnm93nGPMe62SyQS0ZkzZ5SVlSXHSe2X3YRCIeXn56unp4dbOjMAf6+Zhb/XzMLfa+Yx/W8WjUbV39+vvLw8zZkz/mJa4658zJkzR0uWLJnW9/B4PEb+4TA6/l4zC3+vmYW/18xj8t/sSlc8YtLecAoAAGYXwgcAAHDVrAofGRkZev7556lznyH4e80s/L1mFv5eM49NfzPjHjgFAAB2m1VXPgAAQPoRPgAAgKsIHwAAwFWEDwAA4KpZEz52796tv/3bv1VmZqZWrlypTz75JN0jYQwtLS1av3698vLy5DiO3nnnnXSPhHHU1tbqjjvuUFZWlrKzs7VhwwZ1dXWleyyMoa6uTsXFxfGiKr/frw8++CDdY2GCdu7cKcdxtH379nSPMiWzInz88pe/VHV1tZ5//nm1t7drxYoVWrNmjc6ePZvu0TCKgYEBrVixQrt37073KJiA5uZmVVVV6ciRI/rwww918eJF3X///RoYGEj3aBjFkiVLtHPnTrW1tenYsWO677779OCDD+oPf/hDukfDFRw9elSvvfaaiouL0z3KlM2KpbYrV67UHXfcoV/84heShr8/Jj8/X0899ZT++Z//Oc3TYTyO42j//v3asGFDukfBBH3xxRfKzs5Wc3OzysrK0j0OJmDhwoV68cUXtXXr1nSPgjGcP39eJSUleuWVV/STn/xEt956q3bt2pXusSbN+isfFy5cUFtbm8rLy+P75syZo/LycrW2tqZxMsBOwWBQ0vAHGsw2NDSkffv2aWBgQH6/P93jYBxVVVV64IEHEj7LZjLjvlgu1f785z9raGhIOTk5CftzcnL02WefpWkqwE6RSETbt2/X3XffraKionSPgzEcP35cfr9fg4ODWrBggfbv36/CwsJ0j4Ux7Nu3T+3t7Tp69Gi6R0kZ68MHAPdUVVWps7NTH330UbpHwTiWL1+ujo4OBYNB/epXv9KWLVvU3NxMADFQT0+Pnn76aX344YfKzMxM9zgpY334WLx4sebOnau+vr6E/X19ffL5fGmaCrDPtm3b9P7776ulpUVLlixJ9zgYx/z583XjjTdKkkpLS3X06FG9/PLLeu2119I8GUZqa2vT2bNnVVJSEt83NDSklpYW/eIXv1A4HNbcuXPTOOHkWP/Mx/z581VaWqqmpqb4vkgkoqamJu5xAikQjUa1bds27d+/X4cOHVJBQUG6R0KSIpGIwuFwusfAKFatWqXjx4+ro6Mjvt1+++167LHH1NHRMSODhzQLrnxIUnV1tbZs2aLbb79dd955p3bt2qWBgQE98cQT6R4Nozh//rxOnjwZ/7m7u1sdHR1auHChli5dmsbJMJqqqio1NDTo3XffVVZWlgKBgCTJ6/Xq6quvTvN0GKmmpkbr1q3T0qVL1d/fr4aGBh0+fFgHDx5M92gYRVZW1mXPT1177bVatGjRjH6ualaEj+985zv64osv9NxzzykQCOjWW29VY2PjZQ+hwgzHjh3Tt771rfjP1dXVkqQtW7Zo7969aZoKY6mrq5Mk3XvvvQn733jjDT3++OPuD4RxnT17Vps3b1Zvb6+8Xq+Ki4t18OBBrV69Ot2jYRaZFT0fAADAHNY/8wEAAMxC+AAAAK4ifAAAAFcRPgAAgKsIHwAAwFWEDwAA4CrCBwAAcBXhAwAAuIrwAQAAXEX4AAAAriJ8AAAAVxE+AACAq/4fYiCyn6/ys/wAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "t = np.linspace(1, len(solver.losses), len(solver.losses))\n", + "\n", + "plt.scatter(solver.losses, t)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "1afafaab-eef7-4960-b6b9-31de7e12579f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.98435739-0.07941222j, 0.97513866-0.07647791j,\n", + " 0.93544834-0.20615898j, 0.93922148-0.04046039j,\n", + " 0.88351005-0.11470297j, 0.88874281+0.00357627j,\n", + " 0.87630051+0.07383113j, 0.85103414+0.15219365j,\n", + " 0.81860585+0.27813677j])" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Gain solutions\n", + "\n", + "solver.parameter" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "9ff6f0c6-e349-45bf-ad2d-48eddfab3e65", "metadata": {}, "outputs": [], "source": [ - "z = np.zeros((4, 2, 2))" + "G.matrix = G.matrix * solver.parameter" ] }, { "cell_type": "code", - "execution_count": 48, - "id": "f5aff6da-b24c-4151-9aa3-1d9a956fee28", + "execution_count": 35, + "id": "3a47c317-419e-47b3-8743-ccfc45681723", "metadata": {}, "outputs": [], "source": [ - "I = np.identity(2)" + "G_inv = np.linalg.inv(G.matrix)\n", + "G_inv_H = np.linalg.inv(G.matrix.conj().T)" ] }, { "cell_type": "code", - "execution_count": 62, - "id": "ca2e6b88-1816-4698-aff2-801ddac71033", + "execution_count": 36, + "id": "1e33ec99-6df7-4e40-a23e-120a41635f40", "metadata": {}, "outputs": [], "source": [ - "Z = z+I" + "cache = np.matmul(V, G_inv_H)" ] }, { "cell_type": "code", - "execution_count": 61, - "id": "44cee9ec-454b-4da1-93c9-cf9d1351679d", + "execution_count": 37, + "id": "1772d75d-f507-4f73-ae95-b5bdb1aa7186", "metadata": {}, "outputs": [], "source": [ - "p = param.reshape(4, 2, 1)" + "T = np.matmul(G_inv, cache)" ] }, { "cell_type": "code", - "execution_count": 63, - "id": "48429ff8-95bb-4499-b9de-516ffeee3dd0", + "execution_count": 50, + "id": "3c9e9dad-d516-42f2-bcf5-99c7adaa97b2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[[1., 0.],\n", - " [0., 1.]],\n", - "\n", - " [[2., 0.],\n", - " [0., 2.]],\n", - "\n", - " [[3., 0.],\n", - " [0., 3.]],\n", - "\n", - " [[4., 0.],\n", - " [0., 4.]]])" + "" ] }, - "execution_count": 63, + "execution_count": 50, "metadata": {}, "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAGdCAYAAAAhXxuJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKsRJREFUeJzt3X94VPWZ///XZCCToElEIAnBYJAqv3+j+QZqhW2UUmTrXruKlBY2WrqlSQVy1ZVYISKFQLeyeCkSoSBcn0rBdiu1FeGj2Y0sCywQSD+yolRByaJJ4KpmIOrEzjnfP4TRlITJyZmZcybzfFzX+7qa43nPuZkqd+77/T7neEzTNAUAAByT5HQAAAAkOpIxAAAOIxkDAOAwkjEAAA4jGQMA4DCSMQAADiMZAwDgMJIxAAAO6xbrCxqGoffff19paWnyeDyxvjwAwAbTNHX+/Hnl5OQoKSl69dynn36qlpYW25+TnJyslJSUCEQUXTFPxu+//75yc3NjfVkAQATV1dXpuuuui8pnf/rppxpw/dWqbwza/qzs7GydOnXK9Qk55sk4LS1NkvRVfVPd1D3Wl++w+4+edDqELqFvt4+cDiGsbG/A6RA6JCUOOkndPO5f+eour9MhhNUYtF8RRsuFC4ZuzT8X+rs8GlpaWlTfGNSpmuuVntb5f6f85w0NGPeeWlpaSMZ/7VJrupu6q5vHvcm4R5r7/4ONB1d1c/9fzmle98coxUcy7h4Xydj9MX4SdH+MsVhmTE9LspWM40nMkzEAAB0RNA0FbbzKKGgakQsmykjGAABXMmTKUOezsZ25sUYyBgC4kiFDdmpbe7NjKzGa8QAAuBiVMQDAlYKmqaDZ+VaznbmxRjIGALhSIq0Z06YGAMBhVMYAAFcyZCqYIJUxyRgA4Eq0qQEAQMxQGQMAXInd1AAAOMy4OOzMjxe0qQEAcFinkvHatWuVl5enlJQU5efn6+DBg5GOCwCQ4IIXd1PbGfHCcjLevn27SktLVV5eriNHjmjUqFGaMmWKGhsboxEfACBBBU37I15YTsarV6/W3LlzVVRUpKFDh6qyslI9evTQpk2bohEfACBBGREY8cJSMm5paVFNTY0KCwu/+ICkJBUWFmr//v1tzgkEAvL7/a0GAAD4gqVkfO7cOQWDQWVlZbU6npWVpfr6+jbnVFRUKCMjIzRyc3M7Hy0AIGEY8ihoYxjyWL7mnj17NH36dOXk5Mjj8WjHjh1h51RXV2vs2LHy+Xz6yle+os2bN1u+btR3U5eVlampqSk06urqon1JAEAXYJj2h1XNzc0aNWqU1q5d26HzT506pWnTpmny5Mmqra3VggUL9L3vfU+7d++2dF1L9xn37t1bXq9XDQ0NrY43NDQoOzu7zTk+n08+n89SUAAAOGHq1KmaOnVqh8+vrKzUgAED9Pjjj0uShgwZor179+pf//VfNWXKlA5/jqXKODk5WePGjVNVVVXomGEYqqqqUkFBgZWPAgDgiuy0qC8NSZftWwoEAhGLcf/+/a32UUnSlClT2t1H1R7LberS0lJt2LBBW7Zs0fHjxzVv3jw1NzerqKjI6kcBANCuSCXj3NzcVnuXKioqIhZjfX19m/uo/H6/Pvnkkw5/juXHYc6YMUNnz57VkiVLVF9fr9GjR2vXrl2XBQMAgBvU1dUpPT099LMbl0479WzqkpISlZSURDoWAABCDNMjw7S+I/rL8yUpPT29VTKOpOzs7Db3UaWnpys1NbXDn8OLIgAArvTlVnNn50dbQUGBdu7c2erYK6+8YnkfFS+KAADgogsXLqi2tla1tbWSPr91qba2VqdPn5b0+e26s2fPDp3/gx/8QCdPntQ///M/680339TTTz+t559/XgsXLrR0XSpjAIArBZWkoI2aMdiJOYcPH9bkyZNDP5eWlkqS5syZo82bN+uDDz4IJWZJGjBggF566SUtXLhQTzzxhK677jr94he/sHRbk0QyBgC4lGlzzdjsxNxJkybJNNt/WkhbT9eaNGmSjh49avlaX0YyBgC4UjysGUcKa8YAADiMyhgA4EpBM0lB08aacRy9z5hkDABwJUMeGTYauIbiJxvTpgYAwGFUxgAAV0qkDVwkYwCAK9lfM6ZNDQAAOojKGADgSp9v4LLxogja1OHdf/SkeqR5nbp8WOtu/IrTIYR1/4lTTocQlmGjxRQrhtMBdCHeOPjLz+txf4zdXRxiLGMzbD4Ok93UAACgw2hTAwBcKZE2cJGMAQCuZCgpYR76QTIGALhS0PQoaOOtTXbmxhprxgAAOIzKGADgSkGbu6mDtKkBALDHMJNs3R5pxNEGLtrUAAA4jMoYAOBKtKkBAHCYIXs7ouPp6Xq0qQEAcBiVMQDAlew/9CN+6k2SMQDAlew/DjN+knH8RAoAQBdFZQwAcCXeZwwAgMNoU1/Bnj17NH36dOXk5Mjj8WjHjh1RCAsAkOgu3WdsZ8QLy5E2Nzdr1KhRWrt2bTTiAQAg4VhuU0+dOlVTp06NRiwAAIQYpkeGnYd+xNErFKO+ZhwIBBQIBEI/+/3+aF8SANAFGDZbzfF0n3HUI62oqFBGRkZo5ObmRvuSAADElagn47KyMjU1NYVGXV1dtC8JAOgCLr1C0c6IF1FvU/t8Pvl8vmhfBgDQxQTlUdDGvcJ25sZa/PzaAABAF2W5Mr5w4YLefvvt0M+nTp1SbW2trr32WvXv3z+iwQEAEpfdVnOXblMfPnxYkydPDv1cWloqSZozZ442b94cscAAAIktKHut5mDkQok6y8l40qRJMk0zGrEAAJCQeDY1AMCVaFMDAOCwRHpRBMkYAOBKps1XKJrc2gQAADqKyhgA4Eq0qQEAcFgivbUpfn5tAACgi6IyBgC4UtDmKxTtzI01kjEAwJVoUwMAgJihMgYAuJKhJBk2akY7c2ONZAwAcKWg6VHQRqvZztxYi59fGwAA6KKojNtx/4lTTocQ1sabBjgdQlhl75x3OoSwzhstTofQISlew+kQwkqKg9/vfZ7uTocQVh+v0xG0zxfDfw8TaQMXyRgA4Eqmzbc2mTyBCwAAe4LyKGjjZQ925sZa/PzaAABAF0VlDABwJcO0t+5rmBEMJspIxgAAVzJsrhnbmRtr8RMpAABdFJUxAMCVDHlk2NiEZWdurFEZAwBc6dITuOyMzli7dq3y8vKUkpKi/Px8HTx48Irnr1mzRoMGDVJqaqpyc3O1cOFCffrpp5auSTIGAOCi7du3q7S0VOXl5Tpy5IhGjRqlKVOmqLGxsc3zt27dqkWLFqm8vFzHjx/Xxo0btX37dj388MOWrksyBgC40qUNXHaGVatXr9bcuXNVVFSkoUOHqrKyUj169NCmTZvaPH/fvn2aOHGivv3tbysvL0933HGHZs6cGbaa/mskYwCAKxnyhB6J2alxcc3Y7/e3GoFAoM3rtbS0qKamRoWFhaFjSUlJKiws1P79+9ucM2HCBNXU1ISS78mTJ7Vz505985vftPRnJRkDALq03NxcZWRkhEZFRUWb5507d07BYFBZWVmtjmdlZam+vr7NOd/+9rf12GOP6atf/aq6d++ugQMHatKkSZbb1OymBgC4kmlzN7V5cW5dXZ3S09NDx30+n+3YLqmurtaKFSv09NNPKz8/X2+//bbmz5+vZcuWafHixR3+HJIxAMCVIvXWpvT09FbJuD29e/eW1+tVQ0NDq+MNDQ3Kzs5uc87ixYv13e9+V9/73vckSSNGjFBzc7O+//3v6yc/+YmSkjrWgKZNDQBwpVhv4EpOTta4ceNUVVX1RQyGoaqqKhUUFLQ55+OPP74s4Xq9n78D0zQ7/jxOKmMAAC4qLS3VnDlzNH78eN1yyy1as2aNmpubVVRUJEmaPXu2+vXrF1p3nj59ulavXq0xY8aE2tSLFy/W9OnTQ0m5Iywl44qKCv32t7/Vm2++qdTUVE2YMEGrVq3SoEGDrHwMAABhRapNbcWMGTN09uxZLVmyRPX19Ro9erR27doV2tR1+vTpVpXwI488Io/Ho0ceeURnzpxRnz59NH36dC1fvtzSdT2mhTr6G9/4hu69917dfPPN+stf/qKHH35Yx44d0xtvvKGrrrqqQ5/h9/uVkZGh/3N0hHqkdfy3hlhrMd0b2yUbbxrgdAhhlb3z/5wOIaw+3manQ+iQPl7D6RDCutrT3ekQwuqRlOx0CGEFzM+cDqFd/vOGsgfVqampqUPrsJ26xsU8Mf3/3q/uV3X+/6/Pmlv0+zs2RjXWSLFUGe/atavVz5s3b1ZmZqZqamr0ta99LaKBAQCQKGytGTc1NUmSrr322nbPCQQCrW6w9vv9di4JAEgQTrSpndLp3dSGYWjBggWaOHGihg8f3u55FRUVrW62zs3N7ewlAQAJxNbTt2wm8ljrdDIuLi7WsWPHtG3btiueV1ZWpqamptCoq6vr7CUBAOiSOtWmLikp0R/+8Aft2bNH11133RXP9fl8EX3aCQAgMSRSm9pSMjZNUz/60Y/0wgsvqLq6WgMGuH83LwAgPpGM21FcXKytW7fqd7/7ndLS0kIPzs7IyFBqampUAgQAoKuztGa8bt06NTU1adKkSerbt29obN++PVrxAQASlKmLr1Hs5Oj4wyidZ7lNDQBALNCmBgDAYYmUjHlrEwAADqMyBgC4UiJVxiRjAIArJVIypk0NAIDDqIwBAK5kmh6ZNqpbO3NjjWQMAHClS/cL25kfL2hTAwDgMCpjAIArJdIGLpIxAMCVEmnNmDY1AAAOozIGALgSbWoAAByWSG1qx5Jx324f6apu7u2SG6Z7Y7uk7J3zTocQVsXAkU6HENaSk0ecDqFDuns+cTqEsIKeFqdDCMuQ4XQIYX1qBp0OoV0XjNh9f6bNyjiekrH7Mw4AAF0cbWoAgCuZkkzT3vx4QTIGALiSIY88PIELAADEApUxAMCV2E0NAIDDDNMjT4LcZ0ybGgAAh1EZAwBcyTRt7qaOo+3UJGMAgCsl0poxbWoAABxGZQwAcKVEqoxJxgAAV0qk3dQkYwCAKyXSBi7WjAEAcBiVMQDAlT6vjO2sGUcwmCgjGQMAXCmRNnBZalOvW7dOI0eOVHp6utLT01VQUKCXX345WrEBAJAQLCXj6667TitXrlRNTY0OHz6sv/mbv9G3vvUt/c///E+04gMAJCgzAiNeWGpTT58+vdXPy5cv17p163TgwAENGzYsooEBABJbIrWpO71mHAwG9etf/1rNzc0qKCho97xAIKBAIBD62e/3d/aSAAB0SZZvbXr99dd19dVXy+fz6Qc/+IFeeOEFDR06tN3zKyoqlJGRERq5ubm2AgYAJIgE6lNbTsaDBg1SbW2t/vu//1vz5s3TnDlz9MYbb7R7fllZmZqamkKjrq7OVsAAgARxsU3d2aGu3KZOTk7WV77yFUnSuHHjdOjQIT3xxBN65pln2jzf5/PJ5/PZixIAkHB4ApcFhmG0WhMGAADWWKqMy8rKNHXqVPXv31/nz5/X1q1bVV1drd27d0crPgBAgmI3dTsaGxs1e/ZsffDBB8rIyNDIkSO1e/du3X777dGKDwCQqOyu+3bVZLxx48ZoxQEAQMLi2dQAAFdKpA1cJGMAgDvZvVc4jpIx7zMGAMBhVMYAAFdiNzUAAG4QR61mO2hTAwDgMCpjAIAr0aYGAMBp7KYGAMBpnggM69auXau8vDylpKQoPz9fBw8evOL5H330kYqLi9W3b1/5fD7ddNNN2rlzp6VrUhkDAHDR9u3bVVpaqsrKSuXn52vNmjWaMmWK3nrrLWVmZl52fktLi26//XZlZmbqN7/5jfr166f33ntP11xzjaXrkowBAO7kQJt69erVmjt3roqKiiRJlZWVeumll7Rp0yYtWrTosvM3bdqkP//5z9q3b5+6d+8uScrLy7N8XdrUAAB3MiMwJPn9/lajvdf+trS0qKamRoWFhaFjSUlJKiws1P79+9uc8+KLL6qgoEDFxcXKysrS8OHDtWLFCgWDQUt/VJIxAKBLy83NVUZGRmhUVFS0ed65c+cUDAaVlZXV6nhWVpbq6+vbnHPy5En95je/UTAY1M6dO7V48WI9/vjj+ulPf2opRsfa1NnegNK87v1dwHA6gA44b7Q4HUJYS04ecTqEsB67YazTIXTI90+cdDqEsPp1+9DpEMLK6faJ0yGElZHkdTqEdnWL5d1CEXqFYl1dndLT00OHfT6f3chCDMNQZmam1q9fL6/Xq3HjxunMmTP6l3/5F5WXl3f4c1gzBgC4UqTe2pSent4qGbend+/e8nq9amhoaHW8oaFB2dnZbc7p27evunfvLq/3i1+ghgwZovr6erW0tCg5OblDsbq3NAUAIIaSk5M1btw4VVVVhY4ZhqGqqioVFBS0OWfixIl6++23ZRhf9FNPnDihvn37djgRSyRjAIBbRWgDlxWlpaXasGGDtmzZouPHj2vevHlqbm4O7a6ePXu2ysrKQufPmzdPf/7znzV//nydOHFCL730klasWKHi4mJL16VNDQBwpwitGVsxY8YMnT17VkuWLFF9fb1Gjx6tXbt2hTZ1nT59WklJX9Sxubm52r17txYuXKiRI0eqX79+mj9/vh566CFL1yUZAwDwJSUlJSopKWnzn1VXV192rKCgQAcOHLB1TZIxAMCVPObnw878eEEyBgC4UwK9KIJkDABwJwfWjJ3CbmoAABxGZQwAcCfa1AAAOCyBkjFtagAAHEZlDABwpwSqjEnGAAB3Yjc1AACIFSpjAIArJdITuGxVxitXrpTH49GCBQsiFA4AABc58NYmp3Q6GR86dEjPPPOMRo4cGcl4AABIOJ1KxhcuXNCsWbO0YcMG9ezZM9IxAQCQUDqVjIuLizVt2jQVFhaGPTcQCMjv97caAACE49EX68adGk7/ASywvIFr27ZtOnLkiA4dOtSh8ysqKrR06VLLgQEAEhy3NrWtrq5O8+fP13PPPaeUlJQOzSkrK1NTU1No1NXVdSpQAAC6KkuVcU1NjRobGzV27NjQsWAwqD179uipp55SIBCQ1+ttNcfn88nn80UmWgBA4uAJXG37+te/rtdff73VsaKiIg0ePFgPPfTQZYkYAIBOIxm3LS0tTcOHD2917KqrrlKvXr0uOw4AADqGJ3ABAFwpkZ7AZTsZV1dXRyAMAAD+SgK1qXlRBAAADqNNDQBwpwSqjEnGAABXSqQ1Y9rUAAA4jMoYAOBOCfQ4TJIxAMCdWDMGAMBZrBkDAICYoTIGALgTbWoAABxms00dT8mYNjUAAA6jMgYAuBNtagAAHEYyjr4Uj0cpnvi5IduNUryG0yGE1d3zidMhhPX9EyedDqFD1t90g9MhhHX/iVNOhxBWkuec0yGEFTQDTofQrvNB9/+9E4+ojAEArsR9xgAAIGZIxgAAOIw2NQDAndjABQCAsxJpzZhkDABwrzhKqHawZgwAgMOojAEA7sSaMQAAzkqkNWPa1AAAOIzKGADgTrSpAQBwFm1qAAAQM1TGAAB3ok0NAIDDEigZ06YGAMBhlpLxo48+Ko/H02oMHjw4WrEBABLYpQ1cdka8sNymHjZsmF599dUvPqAbnW4AQBQkUJvacibt1q2bsrOzoxELAABfSKBkbHnN+E9/+pNycnJ0ww03aNasWTp9+nQ04gIAIGFYqozz8/O1efNmDRo0SB988IGWLl2qW2+9VceOHVNaWlqbcwKBgAKBQOhnv99vL2IAQEJIpId+WErGU6dODf3vkSNHKj8/X9dff72ef/553X///W3Oqaio0NKlS+1FCQBIPLSpO+aaa67RTTfdpLfffrvdc8rKytTU1BQadXV1di4JAEBUrV27Vnl5eUpJSVF+fr4OHjzYoXnbtm2Tx+PRXXfdZfmatpLxhQsX9M4776hv377tnuPz+ZSent5qAAAQjhO3Nm3fvl2lpaUqLy/XkSNHNGrUKE2ZMkWNjY1XnPfuu+/qxz/+sW699dZO/VktJeMf//jHeu211/Tuu+9q3759+ru/+zt5vV7NnDmzUxcHAKBdZgSGRatXr9bcuXNVVFSkoUOHqrKyUj169NCmTZvanRMMBjVr1iwtXbpUN9xwg/WLymIy/t///V/NnDlTgwYN0j333KNevXrpwIED6tOnT6cuDgBAtPn9/lbjy5uKv6ylpUU1NTUqLCwMHUtKSlJhYaH279/f7uc/9thjyszMbHfvVEdY2sC1bdu2Tl8IAABLIrSBKzc3t9Xh8vJyPfroo5edfu7cOQWDQWVlZbU6npWVpTfffLPNS+zdu1cbN25UbW2tjUB5UQQAwKU8F4ed+ZJUV1fXar+Sz+ezE1bI+fPn9d3vflcbNmxQ7969bX0WyRgA0KV1dPNw79695fV61dDQ0Op4Q0NDm0+efOedd/Tuu+9q+vTpoWOGYUj6/GmVb731lgYOHNihGHlrEwDAnWK8gSs5OVnjxo1TVVVV6JhhGKqqqlJBQcFl5w8ePFivv/66amtrQ+Nv//ZvNXnyZNXW1l7WHr8SKmMAgCs58QSu0tJSzZkzR+PHj9ctt9yiNWvWqLm5WUVFRZKk2bNnq1+/fqqoqFBKSoqGDx/eav4111wjSZcdD4dkDABwJweewDVjxgydPXtWS5YsUX19vUaPHq1du3aFNnWdPn1aSUmRbyqTjAEA+JKSkhKVlJS0+c+qq6uvOHfz5s2duibJGADgXnH0fGk7SMYAAFdKpLc2sZsaAACHURkDANwpgV6hSDIGALgSbWoAABAzVMYAAHeiTQ0AgLMSqU3tWDLu5klSd497u+ReW+8KiY2kOFhlCHpanA4hrH7dPnQ6hA65/8Qpp0MIa+NNA5wOIax7jqc6HUJYg33vOx1Cu5qDhtMhdElUxgAAd6JNDQCAw0jGAAA4K5HWjN2/6AgAQBdHZQwAcCfa1AAAOMtjmvKYnc+odubGGm1qAAAcRmUMAHAn2tQAADiL3dQAACBmqIwBAO5EmxoAAGfRpgYAADFDZQwAcCfa1AAAOIs29RWcOXNG3/nOd9SrVy+lpqZqxIgROnz4cDRiAwAkMjMCI05Yqow//PBDTZw4UZMnT9bLL7+sPn366E9/+pN69uwZrfgAAOjyLCXjVatWKTc3V88++2zo2IABAyIeFAAAUny1mu2w1KZ+8cUXNX78eN19993KzMzUmDFjtGHDhmjFBgBIZKZpf8QJS8n45MmTWrdunW688Ubt3r1b8+bN0wMPPKAtW7a0OycQCMjv97caAADgC5ba1IZhaPz48VqxYoUkacyYMTp27JgqKys1Z86cNudUVFRo6dKl9iMFACQUdlO3o2/fvho6dGirY0OGDNHp06fbnVNWVqampqbQqKur61ykAIDEwm7qtk2cOFFvvfVWq2MnTpzQ9ddf3+4cn88nn8/XuegAAEgAlpLxwoULNWHCBK1YsUL33HOPDh48qPXr12v9+vXRig8AkKA8xufDzvx4YalNffPNN+uFF17Qr371Kw0fPlzLli3TmjVrNGvWrGjFBwBIVLSp23fnnXfqzjvvjEYsAAAkJJ5NDQBwpUTaTU0yBgC4k90Hd8TRQz9IxgAAV0qkytjyW5sAAEBkURkDANzJ7o7oOKqMScYAAFeiTQ0AAGKGyhgA4E7spgYAwFm0qQEAQMxQGQMA3Ind1AAAOIs2NQAAiBkqYwCAOxnm58PO/DjhWDLuLq+6u7gw93o8TocQls/T3ekQwjLk/rd753T7xOkQOiTJc87pEMK653iq0yGE9fyQbKdDCCvnwI1Oh9Culgstkt6NzcVYMwYAwFke2Vwzjlgk0efe0hQAgARBZQwAcCeewAUAgLO4tQkAAMQMlTEAwJ3YTQ0AgLM8pimPjXVfO3NjjTY1AAAOIxkDANzJiMDohLVr1yovL08pKSnKz8/XwYMH2z13w4YNuvXWW9WzZ0/17NlThYWFVzy/PSRjAIArXWpT2xlWbd++XaWlpSovL9eRI0c0atQoTZkyRY2NjW2eX11drZkzZ+o//uM/tH//fuXm5uqOO+7QmTNnLF2XZAwAwEWrV6/W3LlzVVRUpKFDh6qyslI9evTQpk2b2jz/ueee0w9/+EONHj1agwcP1i9+8QsZhqGqqipL1yUZAwDcyYzAkOT3+1uNQCDQ5uVaWlpUU1OjwsLC0LGkpCQVFhZq//79HQr5448/1meffaZrr73W0h+VZAwAcKdLT+CyMyTl5uYqIyMjNCoqKtq83Llz5xQMBpWVldXqeFZWlurr6zsU8kMPPaScnJxWCb0juLUJAOBKkXoCV11dndLT00PHfT6fzcjatnLlSm3btk3V1dVKSUmxNJdkDADo0tLT01sl4/b07t1bXq9XDQ0NrY43NDQoO/vKr978+c9/rpUrV+rVV1/VyJEjLcdImxoA4E4RalN3VHJyssaNG9dq89WlzVgFBQXtzvvZz36mZcuWadeuXRo/fnyn/qiWknFeXp48Hs9lo7i4uFMXBwCgPR7D/rCqtLRUGzZs0JYtW3T8+HHNmzdPzc3NKioqkiTNnj1bZWVlofNXrVqlxYsXa9OmTcrLy1N9fb3q6+t14cIFS9e11KY+dOiQgsFg6Odjx47p9ttv1913323pogAAuNGMGTN09uxZLVmyRPX19Ro9erR27doV2tR1+vRpJSV9UceuW7dOLS0t+od/+IdWn1NeXq5HH320w9e1lIz79OnT6ueVK1dq4MCBuu2226x8DAAA4Tn0PuOSkhKVlJS0+c+qq6tb/fzuu+926hp/rdMbuFpaWvTLX/5SpaWl8ng87Z4XCARa3dPl9/s7e0kAQCJJoLc2dXoD144dO/TRRx/pH//xH694XkVFRav7u3Jzczt7SQAAuqROJ+ONGzdq6tSpysnJueJ5ZWVlampqCo26urrOXhIAkECceDa1UzrVpn7vvff06quv6re//W3Yc30+X9RusAYAdGEOrRk7oVOV8bPPPqvMzExNmzYt0vEAAJBwLFfGhmHo2Wef1Zw5c9StGw/wAgBEialOv5M4ND9OWM6mr776qk6fPq377rsvGvEAACBJttd9u/Sa8R133CEzjv6AAIA4ZcrmmnHEIok6nk0NAIDDWPQFALhTAu2mJhkDANzJkNT+Ax47Nj9O0KYGAMBhVMYAAFdiNzUAAE5LoDVj2tQAADiMyhgA4E4JVBmTjAEA7pRAyZg2NQAADqMyBgC4UwLdZ0wyBgC4Erc2AQDgtARaM3YsGTcGW/RJ0L1L1t3ttEZipI/X6QjC+9QMOh1CWBlJcfBFSgqaAadDCGuw732nQwgr58CNTocQ1vv/33mnQ2jXX8zPnA6hS6IyBgC4k2FKHhvVrUFlDACAPQnUpnZvnxgAgARBZQwAcCmblbHipzImGQMA3Ik2NQAAiBUqYwCAOxmmbLWa2U0NAIBNpvH5sDM/TtCmBgDAYVTGAAB3SqANXCRjAIA7sWYMAIDDEqgyZs0YAACHURkDANzJlM3KOGKRRB3JGADgTrSpAQBArFhKxsFgUIsXL9aAAQOUmpqqgQMHatmyZTLj6LcPAECcMAz7I05YalOvWrVK69at05YtWzRs2DAdPnxYRUVFysjI0AMPPBCtGAEAiSiB2tSWkvG+ffv0rW99S9OmTZMk5eXl6Ve/+pUOHjwYleAAAEgEltrUEyZMUFVVlU6cOCFJ+uMf/6i9e/dq6tSp7c4JBALy+/2tBgAAYV2qjO2MOGGpMl60aJH8fr8GDx4sr9erYDCo5cuXa9asWe3Oqaio0NKlS20HCgBIMAn0BC5LlfHzzz+v5557Tlu3btWRI0e0ZcsW/fznP9eWLVvanVNWVqampqbQqKursx00AABdiaXK+MEHH9SiRYt07733SpJGjBih9957TxUVFZozZ06bc3w+n3w+n/1IAQAJxTQNmTZeg2hnbqxZSsYff/yxkpJaF9Ner1dGHG0fBwDECdO012ruqmvG06dP1/Lly9W/f38NGzZMR48e1erVq3XfffdFKz4AQKIyba4Zd9Vk/OSTT2rx4sX64Q9/qMbGRuXk5Oif/umftGTJkmjFBwBAl2cpGaelpWnNmjVas2ZNlMIBAOAiw5A8NpZBu+qaMQAAMZNAbWpeFAEAgMOojAEArmQahkwbbeoue2sTAAAxQ5saAADECpUxAMCdDFPyJEZlTDIGALiTaUqyc2tT/CRj2tQAADiMyhgA4EqmYcq00aY2qYwBALDJNOyPTli7dq3y8vKUkpKi/Px8HTx48Irn//rXv9bgwYOVkpKiESNGaOfOnZavSTIGALiSaZi2h1Xbt29XaWmpysvLdeTIEY0aNUpTpkxRY2Njm+fv27dPM2fO1P3336+jR4/qrrvu0l133aVjx45Zui7JGACAi1avXq25c+eqqKhIQ4cOVWVlpXr06KFNmza1ef4TTzyhb3zjG3rwwQc1ZMgQLVu2TGPHjtVTTz1l6boxXzO+1MO/cMHdT0bp7nE6gvB8Xnd/h5J0IQ7edd0tDv6/lqTzQfd/l81xEGPLhRanQwjrL+ZnTofQrr/o89hisR77FzNg62UPl2L1+/2tjvt8Pvl8vsvOb2lpUU1NjcrKykLHkpKSVFhYqP3797d5jf3796u0tLTVsSlTpmjHjh2WYo15Mj5//rwk6db8c7G+NABIetfpALqE8+fPKyMjIyqfnZycrOzsbO2tt772+teuvvpq5ebmtjpWXl6uRx999LJzz507p2AwqKysrFbHs7Ky9Oabb7b5+fX19W2eX19fbynOmCfjnJwc1dXVKS0tTR6P/ZLE7/crNzdXdXV1Sk9Pj0CEiYnvMTL4HiOH7zIyIv09mqap8+fPKycnJwLRtS0lJUWnTp1SS4v9LoZpmpflmraqYqfFPBknJSXpuuuui/jnpqen8x9sBPA9RgbfY+TwXUZGJL/HaFXEX5aSkqKUlJSoX+fLevfuLa/Xq4aGhlbHGxoalJ2d3eac7OxsS+e3hw1cAADo8/b4uHHjVFVVFTpmGIaqqqpUUFDQ5pyCgoJW50vSK6+80u757eGhHwAAXFRaWqo5c+Zo/PjxuuWWW7RmzRo1NzerqKhIkjR79mz169dPFRUVkqT58+frtttu0+OPP65p06Zp27ZtOnz4sNavX2/punGfjH0+n8rLy125BhBP+B4jg+8xcvguI4Pv0ZoZM2bo7NmzWrJkierr6zV69Gjt2rUrtEnr9OnTSkr6oqk8YcIEbd26VY888ogefvhh3XjjjdqxY4eGDx9u6boeM56eFwYAQBfEmjEAAA4jGQMA4DCSMQAADiMZAwDgsLhPxlZfdYXWKioqdPPNNystLU2ZmZm666679NZbbzkdVtxbuXKlPB6PFixY4HQocefMmTP6zne+o169eik1NVUjRozQ4cOHnQ4rrgSDQS1evFgDBgxQamqqBg4cqGXLlsXV+30TTVwnY6uvusLlXnvtNRUXF+vAgQN65ZVX9Nlnn+mOO+5Qc3Oz06HFrUOHDumZZ57RyJEjnQ4l7nz44YeaOHGiunfvrpdffllvvPGGHn/8cfXs2dPp0OLKqlWrtG7dOj311FM6fvy4Vq1apZ/97Gd68sknnQ4N7YjrW5vy8/N18803h15VZRiGcnNz9aMf/UiLFi1yOLr4dPbsWWVmZuq1117T1772NafDiTsXLlzQ2LFj9fTTT+unP/2pRo8erTVr1jgdVtxYtGiR/uu//kv/+Z//6XQoce3OO+9UVlaWNm7cGDr293//90pNTdUvf/lLByNDe+K2Mr70qqvCwsLQsXCvukJ4TU1NkqRrr73W4UjiU3FxsaZNm9bq30t03Isvvqjx48fr7rvvVmZmpsaMGaMNGzY4HVbcmTBhgqqqqnTixAlJ0h//+Eft3btXU6dOdTgytCdun8DVmVdd4coMw9CCBQs0ceJEy0+PgbRt2zYdOXJEhw4dcjqUuHXy5EmtW7dOpaWlevjhh3Xo0CE98MADSk5O1pw5c5wOL24sWrRIfr9fgwcPltfrVTAY1PLlyzVr1iynQ0M74jYZI/KKi4t17Ngx7d271+lQ4k5dXZ3mz5+vV155JeZvmulKDMPQ+PHjtWLFCknSmDFjdOzYMVVWVpKMLXj++ef13HPPaevWrRo2bJhqa2u1YMEC5eTk8D26VNwm48686grtKykp0R/+8Aft2bMnKq+47OpqamrU2NiosWPHho4Fg0Ht2bNHTz31lAKBgLxer4MRxoe+fftq6NChrY4NGTJE//Zv/+ZQRPHpwQcf1KJFi3TvvfdKkkaMGKH33ntPFRUVJGOXits148686gqXM01TJSUleuGFF/Tv//7vGjBggNMhxaWvf/3rev3111VbWxsa48eP16xZs1RbW0si7qCJEydedmvdiRMndP311zsUUXz6+OOPW73MQJK8Xq8Mw3AoIoQTt5WxFP5VVwivuLhYW7du1e9+9zulpaWpvr5e0ucvD09NTXU4uviRlpZ22Tr7VVddpV69erH+bsHChQs1YcIErVixQvfcc48OHjyo9evXW34dXaKbPn26li9frv79+2vYsGE6evSoVq9erfvuu8/p0NAeM849+eSTZv/+/c3k5GTzlltuMQ8cOOB0SHFFUpvj2WefdTq0uHfbbbeZ8+fPdzqMuPP73//eHD58uOnz+czBgweb69evdzqkuOP3+8358+eb/fv3N1NSUswbbrjB/MlPfmIGAgGnQ0M74vo+YwAAuoK4XTMGAKCrIBkDAOAwkjEAAA4jGQMA4DCSMQAADiMZAwDgMJIxAAAOIxkDAOAwkjEAAA4jGQMA4DCSMQAADiMZAwDgsP8fBAQl3Rd07CYAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ - "Z*p" + "fig, ax = plt.subplots()\n", + "\n", + "rl = ax.imshow(T.real)\n", + "plt.colorbar(rl)" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "id": "bc987450-b6fd-4cb7-8b82-48a365bd985a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 51, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAGdCAYAAADUoZA5AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALoxJREFUeJzt3Xt0VOW9//FPJpAJaC4iJCEaDNBW7hcTyQrUS0vkInL0LI9VTzwicrDVRIFYV5OeyuVQCLTISUUOCAfUrkLB2mo9VuMPY5GjBgPBdIlF8AopdIgUSSDIJMzM7w9ldGTnMrNnMntP3q+1nrXMzn5mf2cEvvN8n2c/O87n8/kEAABswRHtAAAAQOeRuAEAsBESNwAANkLiBgDARkjcAADYCIkbAAAbIXEDAGAjJG4AAGykR1df0Ov16siRI0pKSlJcXFxXXx4AYILP59PJkyeVmZkphyNyY78zZ86opaXF9OskJCQoMTExDBFZR5cn7iNHjigrK6urLwsACKP6+npdeumlEXntM2fOaOBlF8rV4DH9WhkZGfr4449jKnl3eeJOSkqSJF268GdyWPiD7HnC+rMIF/zd+rvVnros2hF0LOEze1R+Pu/vjXYIHfKknI12CB1K6nM62iF06MqMg9EOoU0tza3acsMf/P+WR+QaLS1yNXj0ce1lSk4K/d/ippNeDcw5qJaWFhK3GefK447EREsn7vhE6yfu+ATrJ26Hdf8X+8U77ZG4HYnWT9y+XtZP3PG9zY/iIi3hwoRoh9ChrpjqTE5ymErcsarLEzcAAJ3h8XnlMTE+8fis/2U3FHyVAQBYklc+0y0Uq1evVnZ2thITE5WXl6eampp2zz9x4oSKiorUv39/OZ1Ofec739GLL74Y0rU7gxE3AMCSvPLKzJg5lN5bt25VSUmJ1q5dq7y8PFVUVGjy5Mnav3+/0tLSzju/paVF1113ndLS0vTMM8/okksu0cGDB5Wammoi8vaRuAEA+NLKlSs1e/ZszZw5U5K0du1a/elPf9LGjRtVWlp63vkbN27U8ePH9eabb6pnz56SpOzs7IjGSKkcAGBJHp/PdJOkpqamgOZ2uw2v19LSotraWhUUFPiPORwOFRQUqLq62rDP888/r/z8fBUVFSk9PV0jRozQ0qVL5fFEbhEkiRsAYEnhmuPOyspSSkqKv5WXlxte79ixY/J4PEpPTw84np6eLpfLZdjno48+0jPPPCOPx6MXX3xRDz/8sB555BH9/Oc/D++H8TWUygEAMa2+vl7Jycn+n51OZ9he2+v1Ki0tTevWrVN8fLxycnJ0+PBh/fKXv9SCBQvCdp2vI3EDACzJK588Ia4MP9dfkpKTkwMSd1v69u2r+Ph4HT16NOD40aNHlZGRYdinf//+6tmzp+Lj4/3Hhg4dKpfLpZaWFiUkhP+efErlAABL6urbwRISEpSTk6OqqqqvYvB6VVVVpfz8fMM+EyZM0AcffCCv96sV7AcOHFD//v0jkrQlEjcAAH4lJSVav369nnrqKe3bt0/33nuvmpub/avM77zzTpWVlfnPv/fee3X8+HHNmTNHBw4c0J/+9CctXbpURUVFEYuRUjkAwJK+vjI81P7BuvXWW/Xpp59q/vz5crlcGjNmjCorK/0L1g4dOhTwVLSsrCy9/PLLmjdvnkaNGqVLLrlEc+bM0U9+8pOQ4+4IiRsAYEneL5uZ/qEoLi5WcXGx4e+2b99+3rH8/Hzt3LkzxKsFj1I5AAA2ElLiDnYfVwAAguX5clW5mRaLgk7c5/ZxXbBggfbs2aPRo0dr8uTJamhoiER8AIBuyuMz32JR0In76/u4Dhs2TGvXrlXv3r21cePGSMQHAOimvGFosSioxB3KPq5ut/u8fWIBAEBogkrcoezjWl5eHrBHbFZWVujRAgC6Da/i5DHRvIqL9luIiIivKi8rK1NjY6O/1dfXR/qSAIAY4PWZb7EoqPu4Q9nH1el0hnVDdwAAurOgRtyh7OMKAEAozJTJz7VYFPTOaSUlJZoxY4Zyc3M1btw4VVRUBOzjCgBAOJhNviTuL3W0jysAAIickPYqb28fVwAAwsHri5PXF/qo2UxfK+MhIwAAS6JUboyHjAAAYCOMuAEAluSRQx4T40tPGGOxEhI3AMCSfCbnuH3McQMA0HWY4zbGHDcAADbCiBsAYEken0Men4k5bvYqBwCg63gVJ6+JwrBXsZm5KZUDAGAjjLgBAJbE4jRjJG4AgCWZn+OmVA4AAKKMETcAwJK+WJxm4iEjlMrDy5t8Vup1NlqX71CLDb7TJH1i/T+UjhYbxNga7Qg6x9Fq/c/S02r9Il5Lq/X/bp9sTYx2CG1qPdt1/4+9Jrc8ZVU5AACIOut/9QQAdEssTjNG4gYAWJJXDjZgMUDiBgBYkscXJ4+JJ3yZ6WtlzHEDAGAjjLgBAJbkMbmq3EOpHACAruP1OeQ1sTjNG6OL0yiVAwBgI4y4AQCWRKncGIkbAGBJXplbGe4NXyiWQqkcAAAbYcQNALAk8xuwxObYlMQNALAk81uexmbijs13BQBAjGLEDQCwJJ7HbYzEDQCwJErlxoJ+Vzt27ND06dOVmZmpuLg4PffccxEICwDQ3Z27j9tMi0VBv6vm5maNHj1aq1evjkQ8AACgHUGXyqdOnaqpU6dGIhYAAPy8vjh5zWzAEqOP9Yz4HLfb7Zbb7fb/3NTUFOlLAgBigNdkuTtW7+OO+LsqLy9XSkqKv2VlZUX6kgAAxKyIJ+6ysjI1Njb6W319faQvCQCIAece62mmxaKIvyun06nk5OSABgBARzyKM91CsXr1amVnZysxMVF5eXmqqanpVL8tW7YoLi5ON910U0jX7azY/DoCAEAItm7dqpKSEi1YsEB79uzR6NGjNXnyZDU0NLTb75NPPtGPf/xjXXXVVRGPMejEferUKdXV1amurk6S9PHHH6uurk6HDh0Kd2wAgG4sGqXylStXavbs2Zo5c6aGDRumtWvXqnfv3tq4cWObfTwejwoLC7Vo0SINGjTIzFvulKDf1e7duzV27FiNHTtWklRSUqKxY8dq/vz5YQ8OANB9eWS2XP6FpqamgPb1O52+rqWlRbW1tSooKPAfczgcKigoUHV1dZtx/ud//qfS0tI0a9asML77tgV9O9i1114rn88XiVgAAAi7b97NtGDBAi1cuPC8844dOyaPx6P09PSA4+np6XrvvfcMX/v111/Xhg0b/FXorsBe5QAASzK7Mvxc3/r6+oCF0U6n03RsknTy5En927/9m9avX6++ffuG5TU7g8QNALCkcD1kpLN3NPXt21fx8fE6evRowPGjR48qIyPjvPM//PBDffLJJ5o+fbr/mNfrlST16NFD+/fv1+DBg0OOvy2sKgcAWJLvy8d6htp8Qd4OlpCQoJycHFVVVfmPeb1eVVVVKT8//7zzhwwZonfeece/YLuurk7/9E//pO9973uqq6uL2IZjjLgBAPhSSUmJZsyYodzcXI0bN04VFRVqbm7WzJkzJUl33nmnLrnkEpWXlysxMVEjRowI6J+amipJ5x0PJxI3AMCSovE87ltvvVWffvqp5s+fL5fLpTFjxqiystK/YO3QoUNyOKJbrCZxAwAsKVpPBysuLlZxcbHh77Zv395u3yeffDKkawaDOW4AAGyEETcAwJI8Jh/raaavlZG4AQCWFK1SudXF5tcRAABiFCNuAIAleeWQ18T40kxfKyNxAwAsyeOLk8dEudtMXyuLza8jAADEqKiNuDMyPlOPC8Kz0XskJA44G+0QOvSJMqMdQodS9kc7go71+oc32iF0is8RH+0QOhR/xvpFPHdr72iH0KF3e56/L7ZVeE4bPxIzElicZsz6f8sAAN2Sz+TTwXwm+loZiRsAYEkexckT5INCvtk/FsXm1xEAAGIUI24AgCV5febmqb2+MAZjISRuAIAleU3OcZvpa2Wx+a4AAIhRjLgBAJbkVZy8JhaYmelrZSRuAIAlsXOaMUrlAADYCCNuAIAlsTjNGIkbAGBJXpnc8jRG57hj8+sIAAAxihE3AMCSfCZXlftidMRN4gYAWBJPBzNG4gYAWBKL04zF5rsCACBGBZW4y8vLdeWVVyopKUlpaWm66aabtH///kjFBgDoxs6Vys20WBRU4n7ttddUVFSknTt3atu2bWptbdWkSZPU3NwcqfgAAN3UuS1PzbRYFNQcd2VlZcDPTz75pNLS0lRbW6urr746rIEBAIDzmVqc1tjYKEnq06dPm+e43W653W7/z01NTWYuCQDoJlhVbizkxWler1dz587VhAkTNGLEiDbPKy8vV0pKir9lZWWFekkAQDfCHLexkBN3UVGR9u7dqy1btrR7XllZmRobG/2tvr4+1EsCANDthVQqLy4u1gsvvKAdO3bo0ksvbfdcp9Mpp9MZUnAAgO6LUrmxoBK3z+fT/fffr2effVbbt2/XwIEDIxUXAKCbI3EbCypxFxUVafPmzfrjH/+opKQkuVwuSVJKSop69eoVkQABAMBXgprjXrNmjRobG3Xttdeqf//+/rZ169ZIxQcA6KZ8Mncvty/abyBCgi6VAwDQFSiVG+MhIwAASyJxG+MhIwAA2AgjbgCAJTHiNkbiBgBYEonbGKVyAABshBE3AMCSfL44+UyMms30tTISNwDAksw+UztWn8dNqRwAABthxA0AsCQWpxkjcQMALIk5bmOUygEAsBFG3AAAS6JUbozEDQCwJErlxqKWuF2HLpajV2K0Lt+xBG+0I+hQ328fj3YIHTquPtEOoUO934h2BJ0TZ/0/kpINYozzWP8fc4/XurOYXRmbz+SIO9TEvXr1av3yl7+Uy+XS6NGjtWrVKo0bN87w3PXr1+vXv/619u7dK0nKycnR0qVL2zw/HKz7pwMAgC62detWlZSUaMGCBdqzZ49Gjx6tyZMnq6GhwfD87du36/bbb9ef//xnVVdXKysrS5MmTdLhw4cjFiOJGwBgST5JPp+JFsI1V65cqdmzZ2vmzJkaNmyY1q5dq969e2vjxo2G52/atEn33XefxowZoyFDhuh//ud/5PV6VVVVZeq9t4fEDQCwpHM7p5lpktTU1BTQ3G634fVaWlpUW1urgoIC/zGHw6GCggJVV1d3KubTp0+rtbVVffpEbpqQxA0AiGlZWVlKSUnxt/LycsPzjh07Jo/Ho/T09IDj6enpcrlcnbrWT37yE2VmZgYk/3BjVTkAwJLCtaq8vr5eycnJ/uNOp9N0bEaWLVumLVu2aPv27UpMjNziaxI3AMCSvL44xYXhPu7k5OSAxN2Wvn37Kj4+XkePHg04fvToUWVkZLTbd8WKFVq2bJleeeUVjRo1KuSYO4NSOQAAkhISEpSTkxOwsOzcQrP8/Pw2+/3iF7/Q4sWLVVlZqdzc3IjHyYgbAGBJ51aHm+kfrJKSEs2YMUO5ubkaN26cKioq1NzcrJkzZ0qS7rzzTl1yySX+efLly5dr/vz52rx5s7Kzs/1z4RdeeKEuvPDC0INvB4kbAGBJ0dg57dZbb9Wnn36q+fPny+VyacyYMaqsrPQvWDt06JAcjq+K1WvWrFFLS4v+5V/+JeB1FixYoIULF4Yce3tI3AAAfE1xcbGKi4sNf7d9+/aAnz/55JPIB/QNJG4AgCWxV7kxEjcAwJLCtao81pC4AQCWFI3FaXbA7WAAANgII24AgCV9MeI2M8cdxmAshMQNALAkFqcZC6pUvmbNGo0aNcq/fVx+fr5eeumlSMUGAAC+IajEfemll2rZsmWqra3V7t279f3vf1833nij3n333UjFBwDopnxhaLEoqFL59OnTA35esmSJ1qxZo507d2r48OFhDQwA0L1RKjcW8hy3x+PR7373OzU3N7e7+brb7Q54aHlTU1OolwQAoNsL+nawd955RxdeeKGcTqd+9KMf6dlnn9WwYcPaPL+8vDzgAeZZWVmmAgYAdBPUyg0Fnbgvv/xy1dXV6a233tK9996rGTNm6K9//Wub55eVlamxsdHf6uvrTQUMAOgmviyVh9pEqfwLCQkJ+ta3viVJysnJ0a5du/SrX/1Kjz/+uOH5TqdTTqfTXJQAgG6HndOMmd45zev1BsxhAwCAyAlqxF1WVqapU6dqwIABOnnypDZv3qzt27fr5ZdfjlR8AIBuilXlxoJK3A0NDbrzzjv197//XSkpKRo1apRefvllXXfddZGKDwDQXZmdpyZxSxs2bIhUHAAAoBPYqxwAYEksTjNG4gYAWJPZe7FjNHHzPG4AAGyEETcAwJJYVW6MxA0AsK4YLXebQakcAAAbYcQNALAkSuXGSNwAAGtiVbkhEjcAwKLivmxm+sce5rgBALARRtwAAGuiVG6IxA0AsCYStyFK5QAA2EjURtwJF32u+N7W/TrUs6cn2iF0aGDqP6IdQocGjzsW7RA6tEvfiXYInTL4d83RDqFDLRclRDuEDn1+sfULjSf/kRLtENrkdZ/puovxWE9D1v8TDADolng6mDFK5QAA2AgjbgCANbE4zRCJGwBgTcxxG6JUDgCAjTDiBgBYUpzvi2amfywicQMArIk5bkMkbgCANTHHbYg5bgAAbIQRNwDAmiiVGyJxAwCsicRtiFI5AAA2wogbAGBNjLgNkbgBANbEqnJDlMoBALARRtwAAEti5zRjpkbcy5YtU1xcnObOnRumcAAA+JIvDC0Eq1evVnZ2thITE5WXl6eampp2z//d736nIUOGKDExUSNHjtSLL74Y2oU7KeTEvWvXLj3++OMaNWpUOOMBACBqtm7dqpKSEi1YsEB79uzR6NGjNXnyZDU0NBie/+abb+r222/XrFmz9Pbbb+umm27STTfdpL1790YsxpAS96lTp1RYWKj169froosuCndMAABExcqVKzV79mzNnDlTw4YN09q1a9W7d29t3LjR8Pxf/epXmjJlih566CENHTpUixcv1hVXXKHHHnssYjGGlLiLioo0bdo0FRQUdHiu2+1WU1NTQAMAoCNx+mqeO6T25et8Mwe53W7D67W0tKi2tjYgtzkcDhUUFKi6utqwT3V19Xm5cPLkyW2eHw5BJ+4tW7Zoz549Ki8v79T55eXlSklJ8besrKyggwQAdEPnbgcz0yRlZWUF5KG28texY8fk8XiUnp4ecDw9PV0ul8uwj8vlCur8cAhqVXl9fb3mzJmjbdu2KTExsVN9ysrKVFJS4v+5qamJ5A0A6DL19fVKTk72/+x0OqMYjXlBJe7a2lo1NDToiiuu8B/zeDzasWOHHnvsMbndbsXHxwf0cTqdtv+QAABREKad05KTkwMSd1v69u2r+Ph4HT16NOD40aNHlZGRYdgnIyMjqPPDIahS+cSJE/XOO++orq7O33Jzc1VYWKi6urrzkjYAACHr4tvBEhISlJOTo6qqKv8xr9erqqoq5efnG/bJz88POF+Stm3b1ub54RDUiDspKUkjRowIOHbBBRfo4osvPu84AAB2U1JSohkzZig3N1fjxo1TRUWFmpubNXPmTEnSnXfeqUsuucQ/Tz5nzhxdc801euSRRzRt2jRt2bJFu3fv1rp16yIWIzunAQAsKRo7p91666369NNPNX/+fLlcLo0ZM0aVlZX+BWiHDh2Sw/FVsXr8+PHavHmzfvazn+mnP/2pvv3tb+u5556L6GDWdOLevn17GMIAAOAbovR0sOLiYhUXFxv+zijn3XLLLbrllltCu1gIeMgIAAA2QqkcAGBNPI/bEIkbAGBJPB3MGKVyAABshBE3AMCavrZtacj9YxCJGwBgTcxxGyJxAwAsiTluY8xxAwBgI4y4AQDWRKncEIkbAGBNJkvlsZq4KZUDAGAjjLgBANZEqdwQiRsAYE0kbkNRS9xnXRfIm5gYrct3yJ3gjXYIHXrXa/3NBRwO6//NKbru/0U7hE5ZrUnRDqFDg/7gjnYIHfIkWH+GsOfJ+GiH0CaH9f8XxzxG3AAAS+I+bmPW/+oJAAD8SNwAANgIpXIAgDWxOM0QiRsAYEnMcRsjcQMArCtGk68ZzHEDAGAjjLgBANbEHLchEjcAwJKY4zZGqRwAABthxA0AsCZK5YZI3AAAS6JUboxSOQAANsKIGwBgTZTKDZG4AQDWROI2RKkcAAAbCSpxL1y4UHFxcQFtyJAhkYoNANCNnVucZqbFoqBL5cOHD9crr7zy1Qv0oNoOAIgASuWGgs66PXr0UEZGRiRiAQDgKyRuQ0HPcb///vvKzMzUoEGDVFhYqEOHDkUiLgAAYCCoEXdeXp6efPJJXX755fr73/+uRYsW6aqrrtLevXuVlJRk2Mftdsvtdvt/bmpqMhcxAKBbYAMWY0El7qlTp/r/e9SoUcrLy9Nll12mp59+WrNmzTLsU15erkWLFpmLEgDQ/VAqN2TqdrDU1FR95zvf0QcffNDmOWVlZWpsbPS3+vp6M5cEAKBbM5W4T506pQ8//FD9+/dv8xyn06nk5OSABgBAR7gdzFhQifvHP/6xXnvtNX3yySd688039c///M+Kj4/X7bffHqn4AADdlS8MLQYFNcf9t7/9Tbfffrv+8Y9/qF+/fvrud7+rnTt3ql+/fpGKDwAAfE1QiXvLli2RigMAgEAsTjPEtmcAAEuK+7KZ6R+LeMgIAAA2wogbAGBNlMoNkbgBAJbEzmnGSNwAAGtixG2IOW4AAGyEETcAwLpidNRsBiNuAIAlWX3L0+PHj6uwsFDJyclKTU3VrFmzdOrUqXbPv//++3X55ZerV69eGjBggB544AE1NjYGdV0SNwAAISgsLNS7776rbdu26YUXXtCOHTt0zz33tHn+kSNHdOTIEa1YsUJ79+7Vk08+qcrKyjafrtkWSuUAAGuy8OK0ffv2qbKyUrt27VJubq4kadWqVbr++uu1YsUKZWZmntdnxIgR+v3vf+//efDgwVqyZInuuOMOnT17Vj16dC4lM+IGAFhSuErlTU1NAc3tdpuOrbq6Wqmpqf6kLUkFBQVyOBx66623Ov06jY2NSk5O7nTSlkjcAIAYl5WVpZSUFH8rLy83/Zoul0tpaWkBx3r06KE+ffrI5XJ16jWOHTumxYsXt1teN0KpHABgTWEqldfX1ys5Odl/2Ol0ttmltLRUy5cvb/dl9+3bZyKoLzQ1NWnatGkaNmyYFi5cGFRfEjcAwJLCtXNacnJyQOJuz4MPPqi77rqr3XMGDRqkjIwMNTQ0BBw/e/asjh8/royMjHb7nzx5UlOmTFFSUpKeffZZ9ezZs1OxnRO1xN3zhEPxidat1HsSrP9cmdPOXtEOoWMO69+E+cHnaR2fZAGzJv452iF0aIO+F+0QOnT5+uPRDqFDjtbUaIfQprOtZ6MdQkT169dP/fr16/C8/Px8nThxQrW1tcrJyZEkvfrqq/J6vcrLy2uzX1NTkyZPniyn06nnn39eiYmJQcdo3cwJAOjefGFoETJ06FBNmTJFs2fPVk1Njd544w0VFxfrtttu868oP3z4sIYMGaKamhpJXyTtSZMmqbm5WRs2bFBTU5NcLpdcLpc8Hk+nr02pHABgTRa+HUySNm3apOLiYk2cOFEOh0M333yzHn30Uf/vW1tbtX//fp0+fVqStGfPHv+K829961sBr/Xxxx8rOzu7U9clcQMALMnqTwfr06ePNm/e3Obvs7Oz5fN9FcS1114b8HOoKJUDAGAjjLgBANZk8VJ5tJC4AQCWFOfzKc5EadlMXyujVA4AgI0w4gYAWBOlckMkbgCAJVl9VXm0UCoHAMBGGHEDAKyJUrkhEjcAwJIolRujVA4AgI0w4gYAWBOlckMkbgCAJVEqNxZ0qfzw4cO64447dPHFF6tXr14aOXKkdu/eHYnYAADdmYUf6xlNQY24P/vsM02YMEHf+9739NJLL6lfv356//33ddFFF0UqPgAA8DVBJe7ly5crKytLTzzxhP/YwIEDwx4UAABS7Ja7zQiqVP78888rNzdXt9xyi9LS0jR27FitX78+UrEBALozn898i0FBJe6PPvpIa9as0be//W29/PLLuvfee/XAAw/oqaeearOP2+1WU1NTQAMAAKEJqlTu9XqVm5urpUuXSpLGjh2rvXv3au3atZoxY4Zhn/Lyci1atMh8pACAboVV5caCGnH3799fw4YNCzg2dOhQHTp0qM0+ZWVlamxs9Lf6+vrQIgUAdC+sKjcU1Ih7woQJ2r9/f8CxAwcO6LLLLmuzj9PplNPpDC06AAAQIKjEPW/ePI0fP15Lly7VD37wA9XU1GjdunVat25dpOIDAHRTcd4vmpn+sSioUvmVV16pZ599Vr/97W81YsQILV68WBUVFSosLIxUfACA7opSuaGgtzy94YYbdMMNN0QiFgAA0AH2KgcAWBKryo2RuAEA1mR2E5UY3YCFxA0AsCRG3MaCfjoYAACIHkbcAABrMrsyPEZH3CRuAIAlUSo3RqkcAAAbYcQNALAmVpUbInEDACyJUrkxSuUAANgII24AgDWxqtwQiRsAYEmUyo1RKgcAwEYYcQMArMnr+6KZ6R+Dopa4E49L8QnRunrHPAlx0Q6hQ74ePaMdQod88db/i/N+U79oh9App3pb+C/Mlx654TfRDqFDD+qOaIfQoW/N2xntENp01tfadRdjjtsQI24AgCXFyeQcd9gisRbmuAEAsBFG3AAAa2LnNEMkbgCAJXE7mDFK5QAA2AgjbgCANbGq3BCJGwBgSXE+n+JMzFOb6WtllMoBALARRtwAAGvyftnM9I9BJG4AgCVRKjdGqRwAgBAcP35chYWFSk5OVmpqqmbNmqVTp051qq/P59PUqVMVFxen5557LqjrkrgBANbkC0OLoMLCQr377rvatm2bXnjhBe3YsUP33HNPp/pWVFQoLi60TVkplQMArMnCO6ft27dPlZWV2rVrl3JzcyVJq1at0vXXX68VK1YoMzOzzb51dXV65JFHtHv3bvXv3z/oazPiBgBY0rmd08w0SWpqagpobrfbdGzV1dVKTU31J21JKigokMPh0FtvvdVmv9OnT+tf//VftXr1amVkZIR0bRI3ACCmZWVlKSUlxd/Ky8tNv6bL5VJaWlrAsR49eqhPnz5yuVxt9ps3b57Gjx+vG2+8MeRrUyoHAFhTmErl9fX1Sk5O9h92Op1tdiktLdXy5cvbfdl9+/aFFM7zzz+vV199VW+//XZI/c8JKnFnZ2fr4MGD5x2/7777tHr1alOBAADwdXHeL5qZ/pKUnJwckLjb8+CDD+quu+5q95xBgwYpIyNDDQ0NAcfPnj2r48ePt1kCf/XVV/Xhhx8qNTU14PjNN9+sq666Stu3b+9UjEEl7l27dsnj8fh/3rt3r6677jrdcsstwbwMAACW1K9fP/Xr16/D8/Lz83XixAnV1tYqJydH0heJ2ev1Ki8vz7BPaWmp/v3f/z3g2MiRI/Vf//Vfmj59eqdjDCpxf/PNLFu2TIMHD9Y111wTzMsAANAxC68qHzp0qKZMmaLZs2dr7dq1am1tVXFxsW677Tb/ivLDhw9r4sSJ+vWvf61x48YpIyPDcDQ+YMAADRw4sNPXDnlxWktLi37zm9/o7rvvbvdeNLfbfd6KPgAAOmTx+7g3bdqkIUOGaOLEibr++uv13e9+V+vWrfP/vrW1Vfv379fp06fDet2QF6c999xzOnHiRIdzAeXl5Vq0aFGolwEAwJL69OmjzZs3t/n77Oxs+ToY9Xf0eyMhj7g3bNigqVOntnuTuSSVlZWpsbHR3+rr60O9JACgGzm3V7mZFotCGnEfPHhQr7zyiv7whz90eK7T6Wx36T0AAIYsPMcdTSGNuJ944gmlpaVp2rRp4Y4HAAC0I+gRt9fr1RNPPKEZM2aoRw/2bwEARIhP5p6pHZsD7uAT9yuvvKJDhw7p7rvvjkQ8AABI4nncbQk6cU+aNCmkVXAAAATFJ5Nz3GGLxFJ4yAgAADbCJDUAwJpYVW6IxA0AsCavpLY35uxc/xhEqRwAABthxA0AsCRWlRsjcQMArIk5bkOUygEAsBFG3AAAa2LEbYjEDQCwJhK3IUrlAADYCCNuAIA1cR+3IRI3AMCSuB3MGIkbAGBNzHEbilriTjp4Vj16no3W5TvkcZqpz3SN+Jb4aIfQIW+89T/HDy9Oi3YInXIkKTnaIXSoV3xrtEPo0Ie3ro12CB0arB9FO4Q2ec+ckUr/GO0wujVG3AAAa/L6pDgTo2YvI24AALoOpXJD3A4GAICNMOIGAFiUyRG3YnPETeIGAFgTpXJDlMoBALARRtwAAGvy+mSq3M2qcgAAupDP+0Uz0z8GUSoHAMBGGHEDAKyJxWmGSNwAAGtijtsQiRsAYE2MuA0xxw0AgI0w4gYAWJNPJkfcYYvEUkjcAABrolRuiFI5AAA2ElTi9ng8evjhhzVw4ED16tVLgwcP1uLFi+WL0W81AIAo8nrNtxgUVKl8+fLlWrNmjZ566ikNHz5cu3fv1syZM5WSkqIHHnggUjECALojSuWGgkrcb775pm688UZNmzZNkpSdna3f/va3qqmpiUhwAAAgUFCl8vHjx6uqqkoHDhyQJP3lL3/R66+/rqlTp7bZx+12q6mpKaABANChcyNuMy0GBTXiLi0tVVNTk4YMGaL4+Hh5PB4tWbJEhYWFbfYpLy/XokWLTAcKAOhm2DnNUFAj7qefflqbNm3S5s2btWfPHj311FNasWKFnnrqqTb7lJWVqbGx0d/q6+tNBw0AQHcV1Ij7oYceUmlpqW677TZJ0siRI3Xw4EGVl5drxowZhn2cTqecTqf5SAEA3YrP55XPxKM5zfS1sqAS9+nTp+VwBA7S4+Pj5Y3RJfcAgCjy+cyVu5njlqZPn64lS5ZowIABGj58uN5++22tXLlSd999d6TiAwB0Vz6Tc9wkbmnVqlV6+OGHdd9996mhoUGZmZn64Q9/qPnz50cqPgAA8DVBJe6kpCRVVFSooqIiQuEAAPAlr1eKMzEVyxw3AABdiFK5IR4yAgBACI4fP67CwkIlJycrNTVVs2bN0qlTpzrsV11dre9///u64IILlJycrKuvvlqff/55p69L4gYAWJLP6zXdIqmwsFDvvvuutm3bphdeeEE7duzQPffc026f6upqTZkyRZMmTVJNTY127dql4uLi8+7Yag+lcgCANVm4VL5v3z5VVlZq165dys3NlfTFAu7rr79eK1asUGZmpmG/efPm6YEHHlBpaan/2OWXXx7UtRlxAwBi2jefl+F2u02/ZnV1tVJTU/1JW5IKCgrkcDj01ltvGfZpaGjQW2+9pbS0NI0fP17p6em65ppr9Prrrwd1bRI3AMCavD7zTVJWVpZSUlL8rby83HRoLpdLaWlpAcd69OihPn36yOVyGfb56KOPJEkLFy7U7NmzVVlZqSuuuEITJ07U+++/3+lrUyoHAFiTzyfJzO1gXyTu+vp6JScn+w+3tw13aWmpli9f3u7L7tu3L6Rwzu0y+sMf/lAzZ86UJI0dO1ZVVVXauHFjp79QkLgBADEtOTk5IHG358EHH9Rdd93V7jmDBg1SRkaGGhoaAo6fPXtWx48fV0ZGhmG//v37S5KGDRsWcHzo0KE6dOhQp+KTSNwAAIvyeX3yxYW+wMwXwuK0fv36qV+/fh2el5+frxMnTqi2tlY5OTmSpFdffVVer1d5eXmGfbKzs5WZman9+/cHHD9w4ICmTp3a6RiZ4wYAWJPPa75FyNChQzVlyhTNnj1bNTU1euONN1RcXKzbbrvNv6L88OHDGjJkiGpqaiRJcXFxeuihh/Too4/qmWee0QcffKCHH35Y7733nmbNmtXpazPiBgBYUjRG3MHYtGmTiouLNXHiRDkcDt1888169NFH/b9vbW3V/v37dfr0af+xuXPn6syZM5o3b56OHz+u0aNHa9u2bRo8eHCnr0viBgAgBH369NHmzZvb/H12drbhl4fS0tKA+7iD1eWJ+9ybOHv2TFdfOigeR1y0Q+iQpyU+2iF0yGuDyRjv5+bv6ewKnnjrx9lyqiXaIXSo6aT1HzzhPWPdfx/PxRbp0awknfW5TZW7z6o1jNFYR5yvKz79r/nb3/6mrKysrrwkACDM6uvrdemll0bktc+cOaOBAwe2eT90MDIyMvTxxx8rMTExDJFZQ5cnbq/XqyNHjigpKUlxceZHtU1NTcrKyjrvPj0Eh88xPPgcw4fPMjzC/Tn6fD6dPHlSmZmZQe2vHawzZ86opcV8BSchISGmkrYUhVK5w+GIyLe0YO7TQ9v4HMODzzF8+CzDI5yfY0pKSlhepz2JiYkxl3DDxQYzkAAA4BwSNwAANmL7xO10OrVgwYJ2955Fx/gcw4PPMXz4LMODzzH2dPniNAAAEDrbj7gBAOhOSNwAANgIiRsAABshcQMAYCO2T9yrV69Wdna2EhMTlZeX5398GjqnvLxcV155pZKSkpSWlqabbrrpvGfFInjLli1TXFyc5s6dG+1QbOfw4cO64447dPHFF6tXr14aOXKkdu/eHe2wbMXj8ejhhx/WwIED1atXLw0ePFiLFy/ukv3FEXm2Ttxbt25VSUmJFixYoD179mj06NGaPHmyGhoaoh2abbz22msqKirSzp07tW3bNrW2tmrSpElqbm6Odmi2tWvXLj3++OMaNWpUtEOxnc8++0wTJkxQz5499dJLL+mvf/2rHnnkEV100UXRDs1Wli9frjVr1uixxx7Tvn37tHz5cv3iF7/QqlWroh0awsDWt4Pl5eXpyiuv1GOPPSbpi33Qs7KydP/995t6ZFp39umnnyotLU2vvfaarr766miHYzunTp3SFVdcof/+7//Wz3/+c40ZM0YVFRXRDss2SktL9cYbb+j//u//oh2Krd1www1KT0/Xhg0b/Mduvvlm9erVS7/5zW+iGBnCwbYj7paWFtXW1qqgoMB/zOFwqKCgQNXV1VGMzN4aGxslffGcWQSvqKhI06ZNC/hzic57/vnnlZubq1tuuUVpaWkaO3as1q9fH+2wbGf8+PGqqqrSgQMHJEl/+ctf9Prrr2vq1KlRjgzh0OUPGQmXY8eOyePxKD09PeB4enq63nvvvShFZW9er1dz587VhAkTNGLEiGiHYztbtmzRnj17tGvXrmiHYlsfffSR1qxZo5KSEv30pz/Vrl279MADDyghIUEzZsyIdni2UVpaqqamJg0ZMkTx8fHyeDxasmSJCgsLox0awsC2iRvhV1RUpL179+r111+Pdii2U19frzlz5mjbtm080cgEr9er3NxcLV26VJI0duxY7d27V2vXriVxB+Hpp5/Wpk2btHnzZg0fPlx1dXWaO3euMjMz+RxjgG0Td9++fRUfH6+jR48GHD969KgyMjKiFJV9FRcX64UXXtCOHTsi8tjVWFdbW6uGhgZdccUV/mMej0c7duzQY489Jrfbrfj4+ChGaA/9+/fXsGHDAo4NHTpUv//976MUkT099NBDKi0t1W233SZJGjlypA4ePKjy8nISdwyw7Rx3QkKCcnJyVFVV5T/m9XpVVVWl/Pz8KEZmLz6fT8XFxXr22Wf16quvauDAgdEOyZYmTpyod955R3V1df6Wm5urwsJC1dXVkbQ7acKECefdjnjgwAFddtllUYrInk6fPi2HI/Cf9/j4eHm93ihFhHCy7YhbkkpKSjRjxgzl5uZq3LhxqqioUHNzs2bOnBnt0GyjqKhImzdv1h//+EclJSXJ5XJJklJSUtSrV68oR2cfSUlJ560LuOCCC3TxxRezXiAI8+bN0/jx47V06VL94Ac/UE1NjdatW6d169ZFOzRbmT59upYsWaIBAwZo+PDhevvtt7Vy5Urdfffd0Q4N4eCzuVWrVvkGDBjgS0hI8I0bN863c+fOaIdkK5IM2xNPPBHt0Gzvmmuu8c2ZMyfaYdjO//7v//pGjBjhczqdviFDhvjWrVsX7ZBsp6mpyTdnzhzfgAEDfImJib5Bgwb5/uM//sPndrujHRrCwNb3cQMA0N3Ydo4bAIDuiMQNAICNkLgBALAREjcAADZC4gYAwEZI3AAA2AiJGwAAGyFxAwBgIyRuAABshMQNAICNkLgBALAREjcAADby/wG0HglqJ6ciXAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "im = ax.imshow(T.imag)\n", + "plt.colorbar(im)" ] }, { "cell_type": "code", "execution_count": null, - "id": "c53f7616-bf7c-4b58-aa27-78841f1b7e87", + "id": "b56446bd-ac61-4e55-8099-ce23e52044be", "metadata": {}, "outputs": [], "source": [] diff --git a/src/calviper/jones.py b/src/calviper/jones.py index 9ab35c6..bdb0002 100644 --- a/src/calviper/jones.py +++ b/src/calviper/jones.py @@ -63,27 +63,40 @@ def calculate(self) -> None: @classmethod def from_visibility(cls: Type[T], dataset: xr.Dataset, time_dependence: bool = False) -> T: + """ + Build Jones matrix from visibility data. + :param dataset: + :param time_dependence: + :return: + """ + + # For the prototype we will pretend there is no polarization and one channel. import calviper.math.tools as tools shape = dataset.VISIBILITY.shape # This will encode the antenna values into an integer list. - index, antenna = tools.encode(dataset.baseline_antenna1_name.to_numpy()) + index, antennas = tools.encode(dataset.baseline_antenna1_name.to_numpy()) - cls._antenna_map = {antenna[i]: index[i] for i in range(len(index))} + instance = cls() # There should be a gain value for each independent antenna. Here we choose antenna_1 names but either # would work fine. - n_parameters = np.unique(dataset.baseline_antenna1_name).shape[0] + instance.n_antennas = np.unique(dataset.baseline_antenna1_name).shape[0] - identity = np.identity(n_parameters, dtype=np.complex64) - instance = cls() + # With no polarization and one channel, n_parameters = n_antennas + # instance.n_parameters = n_parameters + instance.n_parameters = instance.n_antennas + + identity = np.identity(instance.n_parameters, dtype=np.complex64) - instance.n_times, instance.n_antennas, instance.n_channels, instance.n_polarizations = shape + instance._antenna_map = {i: str(antenna) for i, antenna in enumerate(antennas)} - instance.n_parameters = n_parameters + instance.n_times, instance.n_baselines, instance.n_channels, instance.n_polarizations = shape - instance.matrix = np.tile(identity, reps=[*shape, 1, 1]) + # Build on the above idea ... wrong as they may be. Simplicity first. + # instance.matrix = np.tile(identity, reps=[*shape, 1, 1]) + instance.matrix = np.tile(identity, reps=[instance.n_times, 1, 1]) return instance \ No newline at end of file diff --git a/src/calviper/math/tools.py b/src/calviper/math/tools.py index 0fcddce..3983fd4 100644 --- a/src/calviper/math/tools.py +++ b/src/calviper/math/tools.py @@ -34,7 +34,7 @@ def build_visibility_matrix(array: np.ndarray, index_a: np.ndarray, index_b: np. size = index_a.shape[0] # Calculate the N X N matrix size needed - dimension = np.unique(index_a).shape[0] + 1 + dimension = np.unique(index_a).shape[0] # Build matrix matrix_ = np.zeros((dimension, dimension), dtype=np.complex64) diff --git a/testing.ipynb b/testing.ipynb index 9501ee2..09cbefc 100644 --- a/testing.ipynb +++ b/testing.ipynb @@ -57,21 +57,17 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 17, "id": "868b54cd-e0db-4243-bd26-aee45746b24a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([-0.63455564+0.47076233j, 1.15841049-0.81354628j,\n", - " 1.32849216+0.20964037j, 0.7710231 +0.13511301j,\n", - " 0.42394351-0.64036523j, 0.17546983-1.09119839j,\n", - " 0.41194325-0.57720941j, -0.00666993+0.52208018j,\n", - " 0.25330102-0.93767723j, 1.12460157+0.44485419j])" + "(10,)" ] }, - "execution_count": 4, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -79,7 +75,7 @@ "source": [ "df_gains = pd.read_csv(\"gains.csv\")\n", "gains = df_gains.gain.apply(complex).to_numpy()\n", - "gains" + "gains.shape" ] }, { @@ -126,12 +122,24 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 16, "id": "49ad6a73-c7c8-4cb4-b379-264bdcdd8ddb", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "(10, 10)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "X = build_vis(vis, ant_a, ant_b)" + "X = build_vis(vis, ant_a, ant_b)\n", + "X.shape" ] }, { @@ -154,7 +162,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "[\u001b[38;2;128;05;128m2024-12-19 18:54:42,480\u001b[0m] \u001b[38;2;50;50;205m INFO\u001b[0m\u001b[38;2;112;128;144m viperlog: \u001b[0m Iteration: (5)\tStopping criterion reached: 0.009207611105772 \n" + "[\u001b[38;2;128;05;128m2024-12-20 11:08:59,403\u001b[0m] \u001b[38;2;50;50;205m INFO\u001b[0m\u001b[38;2;112;128;144m viperlog: \u001b[0m Iteration: (5)\tStopping criterion reached: 0.009207611105772 \n" ] } ], @@ -418,14 +426,26 @@ "execution_count": 15, "id": "9c1f3593-91e2-49c1-b9b6-1b48e99c28bc", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "(10, 10)" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "vis_array = V[:, 0, 0]\n", "\n", "index_a, _ = cv.math.tools.encode(sps.baseline_antenna1_name.to_numpy())\n", "index_b, _ = cv.math.tools.encode(sps.baseline_antenna2_name.to_numpy())\n", "\n", - "V = cv.math.tools.build_visibility_matrix(array=vis_array, index_a=index_a, index_b=index_b)" + "V = cv.math.tools.build_visibility_matrix(array=vis_array, index_a=index_a, index_b=index_b)\n", + "V.shape" ] }, {