diff --git a/frank/Examples/TestingMmatrixSymmetry.ipynb b/frank/Examples/TestingMmatrixSymmetry.ipynb new file mode 100644 index 00000000..d90dcb25 --- /dev/null +++ b/frank/Examples/TestingMmatrixSymmetry.ipynb @@ -0,0 +1,231 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np \n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "from frank.geometry import SourceGeometry\n", + "from frank.io import load_uvtable\n", + "from frank.radial_fitters import FrankFitter, FourierBesselFitter" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Upload gridded data" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Data source AS209 1mm visibilities gridded\n", + "\n", + "rad_to_arcsec = 3600 * 180 / np.pi\n", + "\n", + "# Huang 2018 \n", + "inc = 34.97\n", + "pa = 85.76\n", + "dra = 1.9e-3\n", + "ddec = -2.5e-3\n", + "r_out = 1.9\n", + "\n", + "# Frank Parameters\n", + "n_pts = 300\n", + "alpha = 1.3\n", + "w_smooth = 1e-1\n", + "\n", + "# UVtable used for the simulations on AS209 at 1mm.\n", + "dir = \"../data/\"\n", + "data_file = dir +'AS209_continuum_prom_1chan_30s_keepflagsFalse_gridded.txt'\n", + "\n", + " # load data\n", + "u, v, Vis, Weights = np.loadtxt(data_file, unpack = True)\n", + "\n", + "geom = SourceGeometry(inc= inc, PA= pa, dRA= dra, dDec= ddec)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Test exploring 2D-DFT with Frank" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "N = 25" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "#FF = FrankFitter(1.9, 300, geom, alpha = alpha, weights_smooth = w_smooth)\n", + "import time\n", + "start_time = time.time()\n", + "FF = FourierBesselFitter(r_out, N, geom)\n", + "sol_new = FF.fit(u, v, Vis, Weights)\n", + "print(\"--- %s minutes ---\" % (time.time()/60 - start_time/60))" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "625" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(sol_new.mean)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "dx = dy = 2*r_out/(N**2*rad_to_arcsec)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "I = sol_new.mean.real/(dx*dy)\n", + "#I" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "I_reshape = I.reshape((N,N))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "r_out_rad = r_out" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGaCAYAAACSWkBBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVrElEQVR4nO3deXxTVd4/8M9NSpuyJGVtWii0LAMKSKVIbQcEhj4UBMYODIIbiwyIigOWkR8gFhhx+ggiiKAVRi2gPDAMiohMx1JBEDplLSgCsrS2A6ZstoFCt+T+/qiNxm45ubdZyOft674wt+fec26aJt+cc+75SrIsyyAiIiKfpXF3A4iIiMi9GAwQERH5OAYDREREPo7BABERkY9jMEBEROTjGAwQERH5OAYDREREPs7P3Q0gIiLvVlJSgrKyMtXO5+/vD51Op9r5qH4MBoiIyGklJSWIiGgLk+m6auc0Go3IyclhQOBCDAaIiMhpZWVlMJmuI/f7f0Cvb6z4fGbzLYR3eBhlZWUMBlyIwQARESmmb6qDvmmg8hNZrcrPQcIYDBARkXJWqzof5AwG3IJ3ExAREfk49gwQEZFy7BnwagwGiIhIOVmu3NQ4D7kchwmIiIh8HHsGiIhIOaus0jABewbcgcEAEREpxzkDXo3DBERERD6OPQNERKQcewa8GoMBIiJSjsGAV+MwAamqtLQUTz75JNq3bw+9Xo/7778fmZmZ7m4WERHVgT0DpKqKigqEh4fjq6++Qrt27fCPf/wDI0eORG5uLpo2beru5hFRQ5FV6hmQ2TPgDuwZIFU1adIESUlJaN++PTQaDcaNGwd/f3+cOXPG5W1JTU2FJEnIzc11ed0NTeTaFi5cCEmScPXq1YZvGPksSbaqtpHrMRigev373/+GJEmQJAnffvtttZ+PHDkS7dq1q/HYs2fP4vr16+jcuXNDN1NVBw4cwMKFC1FYWOjupjikIdtbFXjodDpcvHix2s8HDhyIHj16qF5vbQ4dOoTp06eje/fuaNKkCdq3b4+HH34Y3333XbWye/bssb12f7395z//cVmbiTwdgwGq1/HjxwEAGo0GO3bsqPHn99xzT7X9t2/fxuOPP465c+fCYDA0eDt/7YknnsDt27fRoUMH4WMPHDiARYsWeWww8Otrc0V7S0tL8b//+78Ndn5Hvfrqq9i6dSsGDx6MN954A1OnTsXevXvRu3dvfPPNNzUe8+c//xkbNmyw27wtQPV4VRMI1djI5ThngOp14sQJ6PV6/Pa3v8Wnn36K2bNn2372448/Ij8/H4899pjdMeXl5RgzZgw6d+6MpKQkVzcZAKDVaqHVat1Sd0Nzx7VFRkZi7dq1mDt3LkJDQ11a9y8lJiZi48aN8Pf3t+0bO3Ysevbsif/93//FBx98UO2Y/v37449//KMrm+l7rLI6qwdyBUK3YM8A1ev48ePo2bMnRowYgczMTFy7ds3uZwDsegasViueeOIJSJKEdevWQZKkeuuoGtc+ffo0Hn74Yej1erRs2RIzZsxASUlJtfLHjh3DsGHDoNfr0bRpUwwePLhat29N4+pV9Zw7dw4TJ05EUFAQDAYDJk2ahFu3btnKvPDCCwCAiIgIW7dy1Xlu3LiBmTNnIjw8HAEBAWjTpg3+53/+B0ePHnXsCf3JiRMnIEkStm/fbtt35MgRSJKE3r1725UdNmwYoqOja7y2+tpbpbCwsNZrdsS8efNgsVjc3jsQGxtrFwgAQJcuXdC9e3ecOnWq1uNu3LiBioqKhm4ekVdiMEB1Kisrw5kzZ9CrVy+MGDECFosFO3futP38xIkTAIBevXrZ9j311FP44YcfsGXLFvj5iXU+PfzwwygpKUFycjIefPBBrFy5ElOnTrUrc/LkSfTv3x/Hjx/H7Nmz8dJLLyEnJwcDBw5EVlaWw/XcuHEDycnJePjhh5GamopFixYBAEaNGoVHHnkEALB8+XJbt3Lr1q0BANOmTcPbb7+N0aNH46233sJf/vIXBAYG1vlBVJMePXogKCgIe/fute3bt28fNBoNjh8/DrPZDKAyuDpw4AAeeOCBGs9TX3sduWZHREREYPz48Vi7di0uXbokdK3l5eW4evWqQ5vViW5iWZZRUFCAVq1a1fjzSZMmQa/XQ6fTYdCgQTh8+LBwHVQPNw4TrF69GuHh4dDpdIiOjsbBgwfrLL9lyxZ069YNOp0OPXv2tHtPA4CPPvoIQ4YMQcuWLSFJErKzs+1+fv36dTz33HPo2rUrAgMD0b59e/z5z39GUVGRcNs9hkxUh2PHjskA5JSUFFmWZblnz57ymDFjbD9/8skn5YCAALmiokKWZVnOzc2VAcg6nU5u0qSJbdu7d2+d9SxYsEAGIP/+97+32//MM8/IAOTjx4/b9iUkJMj+/v7y+fPnbfsuXbokN2vWTH7ggQds+95//30ZgJyTk1OtnieffNKunj/84Q9yy5YtbY+XLl1a7dgqBoNBfvbZZ+u8HkcNHz5c7tu3r+3xqFGj5FGjRslarVb+17/+JcuyLB89elQGIH/yySe1Xltd7XX0mmtTVdehQ4fk8+fPy35+fvKf//xn288HDBggd+/evc5z7N69Wwbg0FbTNdRnw4YNMgD53Xfftdu/f/9+efTo0fK7774rf/LJJ3JycrLcsmVLWafTyUePHhWuh6orKiqSAcg/frtWtuR/qHj78du1MgC5qKjIofo3bdok+/v7y++995588uRJecqUKXJQUJBcUFBQY/n9+/fLWq1WXrJkifztt9/K8+fPlxs1aiR//fXXtjLr16+XFy1aJK9dW9mWY8eO2Z3j66+/lkeNGiVv375dPnfunJyRkSF36dJFHj16tNPPo7txzgDVqeqbf9UwwIgRI7B69WqUl5ejUaNGOH78OLp3724bv+7QoQNkBfnIn332WbvHzz33HN566y3s3LkT99xzDywWCz7//HMkJCSgY8eOtnIhISF49NFHsXbtWpjNZuj1+jrrmTZtmt3j/v374+OPP3bo2KCgIGRlZeHSpUuKx8779++P+fPno7i4GE2aNMFXX32Fv/3tb/j++++xb98+DB06FPv27YMkSejXr5+iupRcc5WOHTviiSeewJo1azBnzhyEhIQ4dFyvXr2Qnp7uUFmj0ehQuSqnT5/Gs88+i5iYGEyYMMHuZ7GxsYiNjbU9/v3vf48//vGPuOeeezB37lykpaUJ1UWe5/XXX8eUKVMwadIkAEBKSgo+++wzvPfee5gzZ0618m+88QaGDh1qG1p7+eWXkZ6ejlWrViElJQVA5QRdALXeutujRw9s3brV9rhTp0545ZVX8Pjjj6OiokK4R9QTcJiA6nT8+HFIkoSePXsCqAwGzGYz9u7dC4vFgpMnT9Z4J4GzunTpYve4U6dO0Gg0tj/KK1eu4NatW+jatWu1Y++66y5YrVbk5+fXW0/79u3tHjdv3hxA5YTI+ixZsgTffPMNwsLC0LdvXyxcuBAXLlyo97ia9O/fHxUVFcjMzMSZM2dw+fJl9O/fHw888AD27dsHoHLo4O6770aLFi2cqqOKkmv+pfnz56OiokJo7kDz5s0RFxfn0KbT6Rw+r8lkwvDhw2EwGPDPf/7ToUmVnTt3xkMPPYTdu3fDYrE4XBfVQ5YrFwxSvFV+mTCbzXZbaWlptSrLyspw5MgRxMXF2fZpNBrExcXVuvJpZmamXXkAiI+PV7xSalFREfR6vVcGAgCDAarHiRMn0LFjR9vqgffffz9atWqFTz/9FN999x1KSkrs5guozZHJh86o7UPDkV6Nhx9+GBcuXMCbb76J0NBQLF26FN27d8e//vUv4Xb06dMHOp0Oe/fuxb59+9CmTRv85je/Qf/+/XHw4EGUlpZi37596N+/v/C5f03JNf9Sx44d8fjjj2PNmjX44YcfHDqmrKwMJpPJoc3RD+iioiIMGzYMhYWFSEtLE+qlCQsLQ1lZGYqLix0+huqh8pyBsLAwGAwG25acnFytyqtXr8JisSA4ONhuf3BwMEwmU43NNJlMQuUdcfXqVbz88svV5jd5E+8MYchlTpw4gd/+9re2xxqNBsOGDcOnn36K+++/HwBU7Rk4e/YsIiIibI/PnTsHq9WK8PBwAEDr1q3RuHHjGlc0PH36NDQaDcLCwhS3o74gJCQkBM888wyeeeYZXL58Gb1798Yrr7yCYcOGCdXj7++Pvn37Yt++fWjfvr3tQ79///4oLS3Fhx9+iIKCglonDzraXrXNnz8fH3zwAV599VWHyh84cACDBg1yqGxOTo7t912bkpISjBw5Et999x127dqFu+++26FzV7lw4QJ0Oh2XyPZg+fn5dsNXAQEBbmxN7cxmM4YPH467774bCxcudHdznMZggGplMplw+fLlat/8R4wYgQ0bNuD//u//AEDVnoHVq1djyJAhtsdvvvkmANg+ZLVaLYYMGYJPPvkEubm5tg+NgoICbNy4Ef369XN4/LsuTZo0AYBqi/hYLBbcvHnTbhGlNm3aIDQ0tMZuTEf0798fr7/+Os6fP49Zs2YBAFq1aoW77rrL9mFbX89Abe1tKJ06dcLjjz+Od955Bx06dKi3a1TNOQMWiwVjx45FZmYmPvnkE8TExNRa9sqVK9Xuqjh+/Di2b9+OYcOGQaNh56hqVF5nQK/X1/u33KpVK2i1WhQUFNjtLygoqPV1ZDQahcrX5caNGxg6dCiaNWuGjz/+GI0aNRI+h6dgMEC1qmkNAaByfK1Ro0b49NNPERoaipYtW6pWZ05ODn7/+99j6NChyMzMxAcffIBHH33ULuBYvHgx0tPT0a9fPzzzzDPw8/PDO++8g9LSUixZskSVdkRFRQEAXnzxRYwbNw6NGjXCyJEjUV5ejnbt2uGPf/wjevXqhaZNm2LXrl04dOgQli1bZjtekiQMGDAAe/bsqbeu/v3745VXXkF+fr7dh/4DDzyAd955B+Hh4bUu91xfe6uChIbw4osvYsOGDThz5gy6d+9eZ9mqOQNqmDVrFrZv346RI0fi+vXr1RYZevzxx23/P3bsWAQGBiI2NhZt2rTBt99+izVr1qBx48ZuXy/hjuOGFMb+/v6IiopCRkYGEhISfjrcioyMDEyfPr3GY2JiYpCRkYGZM2fa9qWnp9cZVNbEbDYjPj4eAQEB2L59u9BcF0/EYIBq9es7CaoYDAb069cPu3fvVnWIAAA2b96MpKQkzJkzB35+fpg+fTqWLl1qV6Z79+7Yt28f5s6di+TkZFitVkRHR+ODDz6wW5hHifvuuw8vv/wyUlJSkJaWBqvVipycHISGhuKZZ57B559/jo8++ghWqxWdO3fGW2+9haeffhoAcPPmTQBweKZ9bGwstFotGjdubBf09O/fH++8845D8wVqa29DBgOdO3fG448/jnXr1jVYHTWpuuf7008/xaefflrt578MBhISEvDhhx/i9ddfh9lsRuvWrTFq1CgsWLCAyxHfIRITEzFhwgT06dMHffv2xYoVK1BcXGy7u2D8+PFo27atbc7BjBkzMGDAACxbtgzDhw/Hpk2bcPjwYaxZs8Z2zuvXryMvL8+2nkbVsKTRaITRaITZbMaQIUNw69YtfPDBB7ZJjkDlUKY3rnwqyUruAyNSycKFC7Fo0SJcuXKl1oVjvMXOnTsxYsQI28qNRHcys9kMg8GAH4++CX3TQOXnu3kbzXs/Z5ud74hVq1Zh6dKlMJlMiIyMxMqVK21fDAYOHIjw8HCkpqbaym/ZsgXz589Hbm4uunTpgiVLluDBBx+0/Tw1NdUWTPzSggULsHDhQuzZs6fWOTCOzHnxRAwGyCPcScHACy+8gIsXL2Ljxo3ubgpRg6sKBgoPv6FaMBDUZ4ZQMEDKcZiASGW/HtYgIvJ0DAaIiEg5WbYtGKT4PORyvK+GPMLChQshy7LXDxEQ+Sw3Jioi5RgMEBER+TgOExARkXJuWGeA1MNggIiIlFN5BUJyLQYD9bBarbh06RKaNWvm8vXfiYgagizLuHHjBkJDQ7kkMwFgMFCvS5cuqZL4hojI0+Tn59e71LXDOEzg1RgM1KNZs2Y//V8Tj+oZkGXxPxhJEvsGIEme+fKQIPZ7kCHW7SjLFULlK48R/X0484bnijdJ0W+JouVdcd2e2s0s+v4h/vfn6N945Vpzxb94f1OBVVYpGPDU39+dzTPf7T1IVQAgSZJHBQPibyziaW4963p/JjX4B5Az1y16jBN1yC74fQj/zj3vuj31o0Q0iBX/XYj9zcqy5/6Nk+sxGCAiIuW46JBXYzBARETKcc6AV+M0UiIiIh/HngEiIlJOVmmdAQ4TuAWDASIiUo7DBF6NwwREREQ+jj0DRESkHHsGvBqDASIiUo65CbwagwGHaeDogimiK/05w5m1QoQX63HFdTgxUiX6/IquDig70SYrRFctFP/2I4vW4cQqleLHiD63zrzR3xkrEMrQCpWXnFqtkSO/5BwGA0REpJxsdS4Arek85HIMBoiISDkOE3g19ikRERH5OPYMEBGRcrybwKsxGCAiIuU4TODVOExARETk49gzQEREyllllYYJ2DPgDgwGiIhIOQ4TeDUOExAREfk49gwQEZEKVFp0yKmVF0kpBgNERKQchwm8GocJiIiIfBx7Bhyk0TRyOEGORmr4p1U0+Y6riCYRcua5auhEUFaraNIhQCPYtWmVxeuQBY8RLe/UMR76OvRMos+VWGKjSo7+bTTAt2/2DHg1BgNERKQcVyD0ahwmICIi8nHsGSAiIuU4TODVGAwQEZFyDAa8mlcNE+zduxcjR45EaGgoJEnCtm3b6iy/Z88eSJJUbTOZTK5pMBERkRfwqp6B4uJi9OrVC08++SRGjRrl8HFnzpyBXq+3PW7Tpk1DNI+IyHdxAqFX86pgYNiwYRg2bJjwcW3atEFQUJBDZUtLS1FaWmp7bDabhesjIvI5sly5qXEecjmvGiZwVmRkJEJCQvA///M/2L9/f51lk5OTYTAYbFtYWJiLWklEROQed3QwEBISgpSUFGzduhVbt25FWFgYBg4ciKNHj9Z6zNy5c1FUVGTb8vPzXdhiIiIvVTWBUI2NXM6rhglEde3aFV27drU9jo2Nxfnz57F8+XJs2LChxmMCAgIQEBDgqiYSEd0ZeDeBV7ujewZq0rdvX5w7d87dzSAiIvIYd3TPQE2ys7MREhIifJxW8nc8N4FG/GkVzTUgO5HmU/QYyYlYUTw3QSPhOkSfX9H8B7LWiedW8PfnTP4Dq1wuVL7CWiJcR4PnP3Aml4Fom2ARr8MFJMFcA5ITeTs0Gn+HysmyFRa1nyZZpbsJmO/CLbwqGLh586bdt/qcnBxkZ2ejRYsWaN++PebOnYuLFy9i/fr1AIAVK1YgIiIC3bt3R0lJCf7+97/jiy++wOeff+6uSyAiujNxmMCreVUwcPjwYQwaNMj2ODExEQAwYcIEpKam4ocffkBeXp7t52VlZZg1axYuXryIxo0b45577sGuXbvszkFEROTrvGrOwMCBAyHLcrUtNTUVAJCamoo9e/bYys+ePRvnzp3D7du3ce3aNezevZuBABFRQ7BCpbsJxKtevXo1wsPDodPpEB0djYMHD9ZZfsuWLejWrRt0Oh169uyJnTt32v38o48+wpAhQ9CyZUtIkoTs7Oxq5ygpKcGzzz6Lli1bomnTphg9ejQKCgrEG+8hvCoYICIiD+WmWws3b96MxMRELFiwAEePHkWvXr0QHx+Py5cv11j+wIEDeOSRRzB58mQcO3YMCQkJSEhIwDfffGMrU1xcjH79+uHVV1+ttd7nn38en376KbZs2YIvv/wSly5dEloZ19NIsszlnupiNpthMBjg7xfiWRMInZhkwwmEjnFqciYnEDpYnhMIHS7v4GTAX9JodA6Vq5xAeA1FRUV2S7U7o+o9snDVVOgDxdtc7Xy3yxA0fY3DbYuOjsZ9992HVatWAQCsVivCwsLw3HPPYc6cOdXKjx07FsXFxdixY4dt3/3334/IyEikpKTYlc3NzUVERASOHTuGyMhI2/6ioiK0bt0aGzduxB//+EcAwOnTp3HXXXchMzMT999/vzOX7lbsGSAiIsVkq6zaBlQGGb/cfrlMfJWysjIcOXIEcXFxtn0ajQZxcXHIzMyssZ2ZmZl25QEgPj6+1vI1OXLkCMrLy+3O061bN7Rv317oPJ6EwQARESlXlZtAjQ1AWFiY3dLwycnJ1aq8evUqLBYLgoOD7fYHBwfXmp3WZDIJla/tHP7+/tVy3oiex5N41d0ERETkG/Lz8+2GCbgybMNiMEBERMqpvM6AXq+vd85Aq1atoNVqq83iLygogNForPEYo9EoVL62c5SVlaGwsNCud0D0PJ6EwwRERKScG+4m8Pf3R1RUFDIyMn5uhtWKjIwMxMTE1HhMTEyMXXkASE9Pr7V8TaKiotCoUSO785w5cwZ5eXlC5/Ek7BkgIiKvlZiYiAkTJqBPnz7o27cvVqxYgeLiYkyaNAkAMH78eLRt29Y252DGjBkYMGAAli1bhuHDh2PTpk04fPgw1qxZYzvn9evXkZeXh0uXLgGo/KAHKnsEjEYjDAYDJk+ejMTERLRo0QJ6vR7PPfccYmJivPJOAoDBABERqcFNyxGPHTsWV65cQVJSEkwmEyIjI5GWlmabJJiXlweN5udO8NjYWGzcuBHz58/HvHnz0KVLF2zbtg09evSwldm+fbstmACAcePGAQAWLFiAhQsXAgCWL18OjUaD0aNHo7S0FPHx8XjrrbecvWq34zoD9ai6hzbQP1z4HnoRovdGO7XOgOAxzlyv6LoBWq34fclajdhEItF1BlzBmd+f8DoDFvF1BkTrsFjLxM4vWB4ArA2dPKmyFieOESOaeEgjOfO34fg6A2UVP6i6zsCPr05UbZ2B5v8vVZW2keM4Z4CIiMjHed5XJiIi8jqy/POCQUrPQ67HYICIiJRjCmOvxmECIiIiH8eeASIiUo49A16NwQARESnHYMCrcZiAiIjIx7FngIiIlPtFxkHF5yGXYzBARESKydbKTY3zkOtxmICIiMjHsWeAiIiU4wRCr8ZgwEs5lZtAcP11Z9b012jEjhHNMwAAfoJrtovWoRXMrwAAErRC5UVzUQCAVfB3btGWCtdRYRU7psJ6W6y8Raw8AFitYrkGLLJ4/gPRvmnRfAmAeK4P0VwGztShKgYDXo3DBERERD6OPQNERKQYJxB6NwYDRESknKzSMAFvLXQLDhMQERH5OPYMEBGRctafNjXOQy7HYICIiBSTrTJkFYYJ1DgHieMwARERkY/zqmBg7969GDlyJEJDQyFJErZt21bvMXv27EHv3r0REBCAzp07IzU1tcHbSUTkc6wqbuRyXhUMFBcXo1evXli9erVD5XNycjB8+HAMGjQI2dnZmDlzJv70pz/h3//+dwO3lIiIyHt41ZyBYcOGYdiwYQ6XT0lJQUREBJYtWwYAuOuuu/DVV19h+fLliI+Pb6hmEhH5HvmnTY3zkMt5VTAgKjMzE3FxcXb74uPjMXPmzFqPKS0tRWnpz0uyms3mhmoeEdEdgxMIvZtXDROIMplMCA4OttsXHBwMs9mM27drXiM9OTkZBoPBtoWFhbmiqURERG5zR/cMOGPu3LlITEy0PTabzQgLC/spsUzDRayiyVhEkw4BgCQY+4km3wHEkxuJJh0CAD9NY6Hy/pJYeT/JieRJEDtG40QcbpXEfucVEE9UVKIpEipfbhW77nInElNZBJMnWSziiYpE65CcSFQkyq1Jh5zBdQa82h0dDBiNRhQUFNjtKygogF6vR2BgYI3HBAQEICBA/A2LiMiXMTeBd/Oy0FNMTEwMMjIy7Palp6cjJibGTS0iIiLyPF4VDNy8eRPZ2dnIzs4GUHnrYHZ2NvLy8gBUdvGPHz/eVn7atGm4cOECZs+ejdOnT+Ott97CP/7xDzz//PPuaD4R0Z2L6wx4Na8aJjh8+DAGDRpke1w1tj9hwgSkpqbihx9+sAUGABAREYHPPvsMzz//PN544w20a9cOf//733lbIRGRyjhM4N28KhgYOHAg5DrSW9a0uuDAgQNx7NixBmwVERGRd/OqYICIiDyUDHW6+LnMgFswGCAiIsVkuXJT4zzkel41gZCIiIjUx54BIiJSjBMIvRuDASIiUo4rEHo1DhMQERH5OPYMOMhiLXN4rXDRHACAc7kGRImuda7RiL88tIJrz4uWBwCdRi9UvjEMQuX95ZqXqq5LE1ks/0EjJ/70KgRfIyVSzcm46nJTErv2Eu1NsfKSeBbQCumWUPlyJ65bsoj9bYjmMgAAq2A+A9mJ/nIrHKvDmXPXe04OE3g1BgNERKQY7ybwbhwmICIi8nHsGSAiIuWsUuWmxnnI5RgMEBGRYpwz4N04TEBEROTj2DNARESKybIEWVbexa/GOUgcgwEiIlKMwwTejcMEREREPo49A0REpJgsq9QzwHUG3ILBABERKcY5A96NwwREROTVVq9ejfDwcOh0OkRHR+PgwYN1lt+yZQu6desGnU6Hnj17YufOnXY/l2UZSUlJCAkJQWBgIOLi4nD27Fm7Mt999x0eeughtGrVCnq9Hv369cPu3btVvzZXYc+Ag6xyhVM5BxqKM22RNIK5CSTxl4efYK6BAKmZcB2iuQaaW1sKlW+haSJUHgCaBzQSKt/YT/zbj0Ww+/RGufhze61MLO/DNalIqHyRRux5AoASSfwYTyRbxPrQHc0zYFeHg/30DZGbAFYJshsWHdq8eTMSExORkpKC6OhorFixAvHx8Thz5gzatGlTrfyBAwfwyCOPIDk5GSNGjMDGjRuRkJCAo0ePokePHgCAJUuWYOXKlVi3bh0iIiLw0ksvIT4+Ht9++y10Oh0AYMSIEejSpQu++OILBAYGYsWKFRgxYgTOnz8Po9Go/HlwMUmWOUJTF7PZDIPBAK22pXCin4bkTDDgp9UJlW+kFf9QDPAT+zAJ1DQXrqOpJPbh7pJgwN8TgwHxN/xrZSVi5UWDAemKUHkAKJHFkhuVWsTaBADlVrHkRhUWseep8hixOixymXAdjpJlKyyWaygqKoJeL/Y3+2tV75EXRj2KZo38FbftRnkZOn60Efn5+XZtCwgIQEBA9S8b0dHRuO+++7Bq1SoAgNVqRVhYGJ577jnMmTOnWvmxY8eiuLgYO3bssO27//77ERkZiZSUFMiyjNDQUMyaNQt/+ctfAABFRUUIDg5Gamoqxo0bh6tXr6J169bYu3cv+vfvX9nuGzeg1+uRnp6OuLg4xc+Dq3nOpxsREdFPwsLCYDAYbFtycnK1MmVlZThy5Ijdh69Go0FcXBwyMzNrPG9mZma1D+v4+Hhb+ZycHJhMJrsyBoMB0dHRtjItW7ZE165dsX79ehQXF6OiogLvvPMO2rRpg6ioKMXX7g4cJiAiIsXUnkBYU8/Ar129ehUWiwXBwcF2+4ODg3H69Okaz28ymWosbzKZbD+v2ldbGUmSsGvXLiQkJKBZs2bQaDRo06YN0tLS0Ly5eG+nJ2AwQEREiskqzRmoOoder1c8hNFQZFnGs88+izZt2mDfvn0IDAzE3//+d4wcORKHDh1CSEiIu5sojMMERETklVq1agWtVouCggK7/QUFBbVO4jMajXWWr/q3rjJffPEFduzYgU2bNuG3v/0tevfujbfeeguBgYFYt26dKtfmagwGiIhIMVlWb3OUv78/oqKikJGRYdtntVqRkZGBmJiYGo+JiYmxKw8A6enptvIREREwGo12ZcxmM7Kysmxlbt26BaByfsIvaTQaWK3euZ4yhwmIiEgxdy06lJiYiAkTJqBPnz7o27cvVqxYgeLiYkyaNAkAMH78eLRt29Y2AXHGjBkYMGAAli1bhuHDh2PTpk04fPgw1qxZA6ByPsDMmTOxePFidOnSxXZrYWhoKBISEgBUBhTNmzfHhAkTkJSUhMDAQKxduxY5OTkYPny44ufAHRgMEBGR1xo7diyuXLmCpKQkmEwmREZGIi0tzTYBMC8vz+4bfGxsLDZu3Ij58+dj3rx56NKlC7Zt22ZbYwAAZs+ejeLiYkydOhWFhYXo168f0tLSbGsMtGrVCmlpaXjxxRfxu9/9DuXl5ejevTs++eQT9OrVy7VPgEq4zkA9uM6AGK4z4BiuM+A4rjOgvoZYZ+D08PGqrTPQ7bP1qrSNHMeeASIiUkx0vL+u85DrMRggIiJyk/Xr16t2rvHjxzt9rNcFA6tXr8bSpUthMpnQq1cvvPnmm+jbt2+NZVNTU22TSKoEBASgpES8i4+IiGrHrIXOmThxIiRJnWv2mWBANCEFULlwxZkzZ2yPnX3SZdkCoOFuGRGeA+DE/AUJWqHyGieSyogmNwqQmgrX0cwaJFQ+WCuWsKejXvy6w5uI9W229LcI1yHqSpnY7xsAzt8Qmy+Re0PwuXKmC1jwT9aiKReuwiKLJQXSaMSTCElWwb9ZJ54r2cHraIipYgwGnNerVy889NBDTh+/bds2nDhxQlEbvCoYeP311zFlyhTbt/2UlBR89tlneO+992pMSAFUfvh7YwYpIiLyDZGRkViwYIHTx+fm5ioOBjxnenw9nElIAQA3b95Ehw4dEBYWhoceeggnT56ss57S0lKYzWa7jYiI6maVJdU2X6LX69G4cWNF5wgMDFR854XXBAN1JaSoSh7xa127dsV7772HTz75BB988AGsVitiY2Px3//+t9Z6kpOT7TJlhYWFqXodRER3oqrcBGpsvqSwsNCWftlZb731Fn788UdF5/CaYMAZMTExGD9+PCIjIzFgwAB89NFHaN26Nd55551aj5k7dy6KiopsW35+vgtbTERE5HpeM2fAmYQUv9aoUSPce++9OHfuXK1lAgICakyVSUREteM6A97Na3oGnElI8WsWiwVff/21V6aXJCLyZFaoNGdA9PYRUoXXBANAZUKKtWvXYt26dTh16hSefvrpagkp5s6dayv/17/+FZ9//jkuXLiAo0eP4vHHH8f333+PP/3pT+66BCIionqlpaWhY8eOLqvPa4YJAPGEFD/++COmTJkCk8mE5s2bIyoqCgcOHMDdd9/trksgIrojcZ0BdRUXF+P77793WX1eFQwAwPTp0zF9+vQaf7Znzx67x8uXL8fy5ctd0CoiIt8mq3Rb4J0eDCQlJTlU7tSpUw3cEnteFwwQERF5q8WLFyMoKAgGg6HOcrdu3XJRiyoxGCAiIsU4TOCYTp06oX///njvvffqLPfPf/4TY8eOdVGrGAw4rHLNb0dfpE7kDRDMNSCaAwAANBqxY4TzJQDwk3RC5f0RKFxHc8F8Bh2aia2fHxkkvu78fW2uCZU3thVf2VL0TbLgklhOBgA4dLml4BFi+etLzc0Fzw+U4rZQ+XJJrDwAWDSlYuWtYuVdRZYdy5/SELkJrFAne0vDZYDxDDExMThw4EC95SRJapDfU20YDBAREbnIjBkzsH///nrLDRgwALt373ZBiyoxGCAiIsU4TOCYqKgoREVF1VuuVatWGDBggAtaVInBABERKWaVocrdBFauQOgWXrXoEBER0Z0oPz8fv/vd79xWP3sGiIhIMQ4TKHPr1i18+eWXbqufwQARESlWOUygznnI9ThMQERE5OPYM0BERIpxmMC7MRggIiLFrFAn/TBTGLsHhwmIiIh8HHsGiIhIMVmu3NQ4D7kegwEiIlLMqlIKYzXO4a1cmYvg1xgMOKgyAYhjL1LJBa9l0cRGzhzjpwkQrsNPEjsmQBZPVBTUSCw5TvsmYn9gvVv9KFQeADrGlQiV1/a9S7gOWMVSuDQ5nCNeRfp1ofLXyloLlb9yW+x3V1mHWGKq2xrxJFCiib9c8fcHB5MO2XP0GH799jRGoxFvv/222+pnMEBERIrJKk0glH10AqHBYMBTTz3ltvo5gZCIiBSrmjOgxubLnn/+eWzcuNHl9TIYICIi8hCrVq3Chx9+6PJ6OUxARESKcQKhOtq1a4dmzZq5vF6HgoETJ04In/juu++Gnx9jDSIiXyBDUmW831fnDFT5wx/+gE2bNqGsrAz+/uITbp3l0Kd1ZGQkJEly+LYHjUaD7777Dh07dlTUOCIiIl/y0ksvYcuWLZgxY4ZL7y5w+Kt7VlYWWreu/zYiWZbRo0cPRY0iIiLvwqyF6hg1ahR69eqFNWvW4OTJk1i6dCmio6MbvF6HgoEBAwagc+fOCAoKcuikDzzwAAIDxe8fJyIi78Q5A+r48ssvbf//1VdfITY2FkajEQ888ADuvfde29aqVStV63UoGNi9e7fQSXfu3OlUY4iIiHxZTk4OsrOzcfz4cWRnZyM7Oxu5ubnYvHkzNm/eDOmnVe1CQ0PRu3dvfPLJJ6rUyxl+RESkGCcQqqNDhw7o0KEDHnroIds+s9lsCw6q/j158iR27NihWr0OBwN//etf7R4nJSWp1ggiIvJunDOgjm+//RYFBQXo27cvmjRpAgDQ6/Xo378/+vfvbytnsVhw+vRp1ep1OBjIyfl5nXPJFYvvk+okwTWmRNdrBwA/iOUmaCQ3Eq4jUCv2+mvlL7bGe0hokVB5AND2+Y1QeUu/WOE6RGmdOCbk2++Eyre80kKofKCf+Dpn/mVirxEtxF9TriA7lWuAfM3SpUuxfv16/Oc//8F9991n219QUIB3330XVqsVv//973HPPfege/fuqtXr8Lv9+++/r1qlRER0Z+EwgToyMzPRuXNnu0CgtLQUMTEx+P777yHLMhYuXIjk5GS88MILqtXL5YiJiEixqmECNTZf9sMPP+A3v7Hvady0aRNyc3PRp08fLF++HJ06dcKcOXOwf/9+1ep1OBiYNGkSnnzySdvmLqtXr0Z4eDh0Oh2io6Nx8ODBOstv2bIF3bp1g06nQ8+ePXmnAxEReazS0lI0bWqftnvr1q3QarXYvHkzZsyYgV27dsHPzw9vvPGGavU6HAyEh4fbZjl26NBBtQaI2Lx5MxITE7FgwQIcPXoUvXr1Qnx8PC5fvlxj+QMHDuCRRx7B5MmTcezYMSQkJCAhIQHffPONi1tORHRnq1pnQI3Nl7Vt2xa5ubm2x7du3UJGRgZiY2MRHh4OAAgLC0P//v1V7RlweM7AggULVKvUWa+//jqmTJmCSZMmAQBSUlLw2Wef4b333sOcOXOqlX/jjTcwdOhQ27jKyy+/jPT0dKxatQopKSkubTsR0Z1M/mlT4zy+bODAgVi3bh1OnDiBe+65B+vXr8ft27cxbNgwu3JGoxH79u1TrV6vmTNQVlaGI0eOIC4uzrZPo9EgLi4OmZmZNR6TmZlpVx4A4uPjay0PVHbRmM1mu42IiDyX2sPHsiwjKSkJISEhCAwMRFxcHM6ePVvtPJ999hmio6MRGBiI5s2bIyEhQfG1vPDCC2jUqBF+97vf4Q9/+AMSExOh1WoxduxYu3LXrl2DXq9XXF8Vh4KB3r1748cff3T4pP369cPFixedblRNrl69CovFguDgYLv9wcHBMJlMNR5jMpmEygNAcnIyDAaDbQsLC1PeeCKiO5wMdYYIRO8maIjh4yVLlmDlypVISUlBVlYWmjRpgvj4eJSUlNjKbN26FU888QQmTZqE48ePY//+/Xj00Uede/J+oVu3bvj444+h0+nwySefoLS0FIsWLUJERIStjNVqxaFDh9CuXTvF9VVxaJigatWjFi0cu6c4OzsbpaWlihrmLnPnzkViYqLtsdlsZkBARFQP60+bGucRofbwsSzLWLFiBebPn29bBXD9+vUIDg7Gtm3bMG7cOFRUVGDGjBlYunQpJk+ebDv33Xff7dxF/8rQoUORl5eHs2fPwmAwwGg02v38888/x/Xr16v1Fijh8JyBwYMHO5zCuCEWJWrVqhW0Wi0KCgrs9hcUFFR7oqoYjUah8gAQEBCAgACxhXOIiEhdvx6irem9uWr4eO7cubZ9jgwf//ILH1A5fLxt2zYAlQvsmUwmuyFmg8GA6OhoZGZmYty4cTh69CguXrwIjUaDe++9FyaTCZGRkVi6dKlqWXs1Gg26du1a688mTZqEUaNGqVIX4GAw8MvVBx2lZvcFAPj7+yMqKgoZGRm2cRmr1YqMjAxMnz69xmNiYmKQkZGBmTNn2valp6cjJiZG1bYREfk6WZYgq3AnQNU5ft0ju2DBAixcuNBuX13Dx7Ut1Vvf8HHVv3WVuXDhAgBg4cKFeP311xEeHo5ly5Zh4MCB+O677xzuRXfWkCFDMGTIEFXP6VAw4K5bCX8tMTEREyZMQJ8+fdC3b1+sWLECxcXFtu6h8ePHo23btkhOTgYAzJgxAwMGDMCyZcswfPhwbNq0CYcPH8aaNWvceRlERHcctYcJ8vPz7SbIeVKPrdVa2coXX3wRo0ePBlC5Sm+7du2wZcsWPPXUU+5snlO8Kmvh2LFjceXKFSQlJdm6ZdLS0mwRXF5eHjSan+dExsbGYuPGjZg/fz7mzZuHLl26YNu2bap14xARUcPQ6/X1zpZviOHjqn8LCgoQEhJiVyYyMhIAbPt/OUcgICAAHTt2RF5engNX53m8KhgAgOnTp9c6LLBnz55q+8aMGYMxY8Y0cKt+zTPv2JSkhm+XRjA9jp8T6XQ0gj2RkuCdyy54mjyWaC4d0Znfziw1WwGLWB2C5Ukd7sha2BDDxxERETAajcjIyLB9+JvNZmRlZeHpp58GAERFRSEgIABnzpxBv379AADl5eXIzc31mJ50UV4XDBARkedxV6IitYePJUnCzJkzsXjxYnTp0gURERF46aWXEBoaags49Ho9pk2bhgULFiAsLAwdOnTA0qVLAcANXz7VwWCAiIi8VkMMH8+ePRvFxcWYOnUqCgsL0a9fP6SlpUGn09nKLF26FH5+fnjiiSdw+/ZtREdH44svvkDz5s1dd/EqkmRH7xf8yYQJEzB58mQ88MADDdUmj2I2m2EwGAA0dviWSUkSj7E0gsf4aRsL1xHQyCBUXucXJFyHXhNSf6FfaGMNFa6jW+MgofKDgsX6vkdE5gqVBwD9+N/UX+gXLA/0E65DlHbvV8LHFKV+J1R+x4mI+gv9whcm8W+O396+JlT+quaScB03rTUvUFObkopC4TpKy4uEyldYbgnXYZXLHConyzJk+SaKiooUr2JX9R65LvIFNNYqn+R3y1KKCdlLVWkbOU54dLSoqAhxcXHo0qUL/va3v6m+0iAREXmfqmECNTZftmjRIvz3v/91eb3CwcC2bdtw8eJFPP3009i8eTPCw8MxbNgw/POf/0R5eXlDtJGIiMgnVC09PHLkSGzfvt12G2NDc2redOvWrZGYmIjjx48jKysLnTt3xhNPPIHQ0FA8//zzNSZ0ICKiO1fV3QRqbL5s8eLFaN++PT777DP84Q9/QFhYGF566SW7tMYNQdFNVD/88APS09ORnp4OrVaLBx98EF9//TXuvvtuLF++XK02EhGRh2MwoI558+bh/Pnz+PzzzzFmzBhcu3YNr7zyCjp37oyhQ4di69atqKioUL1e4WCgvLwcW7duxYgRI9ChQwds2bIFM2fOxKVLl7Bu3Trs2rUL//jHP/DXv/5V9cYSERH5gri4OGzatAkXL17Ea6+9hq5du+Lzzz/Hww8/jHbt2mHOnDmq9sILBwMhISGYMmUKOnTogIMHD+Lw4cOYNm2a3azPQYMGISgoSLVGEhGRZ+MEwobRsmVLJCYm4uTJk/jqq6/wyCOP4PLly1i6dCm6deuGwYMH4+OPP1Zcj/A9cMuXL8eYMWPs7rf8taCgIKeSGxERkXeSVeriF7vZ3XecP38en376KTIyMmz72rVrh927d2PPnj2IiorC1q1bqyV4cpRwz8ATTzxRZyBAREREypWXl2PTpk0YPHgwfvOb3+DVV19FRUUFEhMTcfr0aXz//ffYv38/hg0bhsOHD9e6BLMjuAKhgyRJ4/CiQ67gijwDzrBA7PbSMsHyAHC7Quyrw7UysfwHP1wUX+ikyWGxnjDxjAwANGK/c4tgmwDgh0tiC1NdLRNr022L+G1S5ZLYa6QCpcJ1WGWxCVmi5X2B2lkLfdmpU6ewdu1abNiwAdevX4csy4iNjcW0adMwZswYuwyOMTEx2LFjB+6//358+eWXTtfJYICIiBSTZQmyrEJuAhXO4c369euHzMxMyLIMvV6Pp59+GtOmTas322737t1x6NAhp+tlMEBEROQhDhw4gN69e2PatGl49NFH0bixY0vP/+lPf1KUJoDBABERKcZhAnUcOnQIUVFRwsfFxMTY0jA7g8EAEREpptaCQb6+6JAzgYAaPHMWGhEREbkMewaIiEgx+adNjfOQ6zEYICIixSqHCZTfCeBrwwRarVM3GkOSJDRp0gTt27fHoEGDkJiYiPDwcKfbwWECIiIiN5Fl2anNarXixo0bOHnyJFatWoV7770XX3/9tdPtYDBARESKySpuvsRqtTq93b59GydPnsS0adNQVFSEBQsWON0ODhMQEZFivJvA9QICAnDXXXfhrbfewv79+7F3716nz8WeASIiIi/Xo0cPFBYWOn08gwEiIlLMquLmS65fv45bt24pOsetW7fw5ptv4vz5806fg8MEDtMADubZdiqJkAcmHnImGUuFLJYk5pbmpnAd18vFEglduBkoVP7IlRZC5QFAs+uaUPngb78TrkOSxPpPTU4kXDp6tblQ+bxisdnjheVlQuUBoFS6LVRe9DUIABZZLBmSLIt/ZDlzjChJcvQtXVY9VbAsq5N+2NdSGLdu3RoTJ07Eu+++6/Q5nn32WWzYsAEVFc4n0PK8TyAiIiIfUXV3gBrnUYI9A0REpJgMCVYHe0/rO4+v+eqrr/Dkk08qOl4pBgNERKQYhwmcd+7cOZw7d07ROSRJWRDFYICIiMhNdu/e7e4mAGAwQEREKmAKY+cMGDDA3U0A4EUTCK9fv47HHnsMer0eQUFBmDx5Mm7erHsm+sCBAyFJkt02bdo0F7WYiMh3VC06pMZGruc1PQOPPfYYfvjhB6Snp6O8vByTJk3C1KlTsXHjxjqPmzJlCv7617/aHjdu3Lihm0pERORVvCIYOHXqFNLS0nDo0CH06dMHAPDmm2/iwQcfxGuvvYbQ0NBaj23cuDGMRqOrmkpE5JOYwti7ecUwQWZmJoKCgmyBAADExcVBo9EgKyurzmM//PBDtGrVCj169MDcuXPrXemptLQUZrPZbiMiorpxmMC7eUXPgMlkQps2bez2+fn5oUWLFjCZTLUe9+ijj6JDhw4IDQ3FiRMn8P/+3//DmTNn8NFHH9V6THJyMhYtWqRa24mIiDydW4OBOXPm4NVXX62zzKlTp5w+/9SpU23/37NnT4SEhGDw4ME4f/48OnXqVOMxc+fORWJiou2x2WxGWFiY020gIvIFXGfAu7k1GJg1axYmTpxYZ5mOHTvCaDTi8uXLdvsrKipw/fp1ofkA0dHRACoXeKgtGAgICEBAQEC1/RpJ61zOgQbizDrnorkGRNdrB4AyWSzhRokknpvgMn4UKh9wQytYg79geeB6eZv6C/1C6ysthesQda1M9LqBXMFcA9/fEHuN/CiL/77LBHMTOPO6Ff3bcCo3gQtumpMcHvlVvy28tdC7uTUYaN26NVq3bl1vuZiYGBQWFuLIkSOIiooCAHzxxRewWq22D3hHZGdnAwBCQkKcai8REdGdyHO+6tbhrrvuwtChQzFlyhQcPHgQ+/fvx/Tp0zFu3DjbnQQXL15Et27dcPDgQQDA+fPn8fLLL+PIkSPIzc3F9u3bMX78eDzwwAO455573Hk5RER3HE4gdI5Wq3V68/NT7/u8V0wgBCrvCpg+fToGDx4MjUaD0aNHY+XKlbafl5eX48yZM7a7Bfz9/bFr1y6sWLECxcXFCAsLw+jRozF//nx3XQIR0R2LtxY6R0m2QTWyHVbxmmCgRYsWdS4wFB4ebvfEhIWF4csvv3RF04iIiJxitXrGLAmvCQaIiMhzqdXF72vDBJ6CwQARESkmQ4IMZWl0q85DrucVEwiJiIio4bBngIiIFJOhThc/Rwncg8EAEREpxjkD3o3DBERERD6OPQNERKQY1xnwbgwGHCVpKjdHOLFuuSinchNYBXMTSKXCdVQIdjaVSI2E6xCebCz47nKryCBYAVBwWydUvlkj8bwBWsHrvlUh/rb6Y1mZUPnr1mKh8jc0hULlAaBUMJ+BRRZ/3VqtYvkMrE7kPxD+m/WgXCiO4DCBd/OuVxsREdGvrF69GuHh4dDpdIiOjrYtS1+bLVu2oFu3btDpdOjZsyd27txp93NZlpGUlISQkBAEBgYiLi4OZ8+erfFcpaWliIyMhCRJtvw33ojBABERKSar+J+IzZs3IzExEQsWLMDRo0fRq1cvxMfHV8t0W+XAgQN45JFHMHnyZBw7dgwJCQlISEjAN998YyuzZMkSrFy5EikpKcjKykKTJk0QHx+PkpKSauebPXu2LUeON2MwQEREirkrUdHrr7+OKVOmYNKkSbj77ruRkpKCxo0b47333qux/BtvvIGhQ4fihRdewF133YWXX34ZvXv3xqpVqwBU9gqsWLEC8+fPx0MPPYR77rkH69evx6VLl7Bt2za7c/3rX//C559/jtdee82Zp8yjMBggIiKPYzab7bbS0upzQcrKynDkyBHExcXZ9mk0GsTFxSEzM7PG82ZmZtqVB4D4+Hhb+ZycHJhMJrsyBoMB0dHRducsKCjAlClTsGHDBjRu3FjRtXoCBgNERKSYrOIGVCabMxgMti05OblanVevXoXFYkFwcLDd/uDgYJhMphrbaTKZ6ixf9W9dZWRZxsSJEzFt2jT06dOnnmfGO/BuAiIiUkztuwny8/Oh1+tt+wMCApSfXCVvvvkmbty4gblz57q7KaphzwAREXkcvV5vt9UUDLRq1QparRYFBQV2+wsKCmA0Gms8r9ForLN81b91lfniiy+QmZmJgIAA+Pn5oXPnzgCAPn36YMKECU5crfsxGCAiIsVkWb3NUf7+/oiKikJGRoZtn9VqRUZGBmJiYmo8JiYmxq48AKSnp9vKR0REwGg02pUxm83IysqylVm5ciWOHz+O7OxsZGdn225N3Lx5M1555RXHL8CDcJiAiIgUs/60qXEeEYmJiZgwYQL69OmDvn37YsWKFSguLsakSZMAAOPHj0fbtm1tcw5mzJiBAQMGYNmyZRg+fDg2bdqEw4cPY82aNQAASZIwc+ZMLF68GF26dEFERAReeuklhIaGIiEhAQDQvn17uzY0bdoUANCpUye0a9fO+Yt3IwYDRETktcaOHYsrV64gKSkJJpMJkZGRSEtLs00AzMvLg0bzcyd4bGwsNm7ciPnz52PevHno0qULtm3bhh49etjKzJ49G8XFxZg6dSoKCwvRr18/pKWlQacTW2nUm0iyLNIp43vMZjMMBgO02paQGnI5YsGlR7WSv3AVftpAwfLiL/xGGrE6ArTiS//qJH39hX7BILcWKt9SFm9TS3/R5YjFR+hcsxyx2DK7ossR/6i5JlQeAG7KYsfctv4oXEdphVmofLlF7LoBoMJSfcGausiqfM+u5dyyFRbLNRQVFdlN0nNG1XvkU+3mIkCj/MOy1FqCd/6brErbyHHsGSAiIuUEx/vrOg+5HoMBB0nQQHJwvqUsmkjHCVZZLOkQAFRYbjdAS1yvQiP2DcuiEfu2e0sS+5YIAD+Wi/Um6MrEelAAwE9wvm85xF8jtyWx57ZEI/YN+RaKhMoDQJl8S6h8hdWJREWCf0+iSb+c4ej7DZEaGAwQEZFi7ppASOpgMEBERIqJ3hZY13nI9dgPRURE5OPYM0BERIpxmMC7MRggIiLFZFmGGneq82539+AwARERkY9jzwARESmmdtZCci0GA0REpJgMddYLYizgHhwmICIi8nFeEwy88soriI2NRePGjREUFOTQMbIsIykpCSEhIQgMDERcXBzOnj3bsA0lIvJBVcMEamzkel4TDJSVlWHMmDF4+umnHT5myZIlWLlyJVJSUpCVlYUmTZogPj4eJSViS64SEVHdGAx4N6+ZM7Bo0SIAQGpqqkPlZVnGihUrMH/+fDz00EMAgPXr1yM4OBjbtm3DuHHjajyutLQUpaU/r21uNleuUy9JGoGshY4Vsz9E7O5a2YncBBbRO3gtwlVAFszYaNWKX4fGKpZjQbSOMs1NofIAUCaJtclPChCuQyMYu1uduGO7AmLr+lfIYuVF8wwAQIVV7BiLM7kJBHMNNGRGwSoOv98QqeCOfbXl5OTAZDIhLi7Ots9gMCA6OhqZmZm1HpecnAyDwWDbwsLCXNFcIiKvVjmBUI3/yB3u2GDAZDIBAIKDg+32BwcH235Wk7lz56KoqMi25efnN2g7iYjuBBwm8G5uDQbmzJkDSZLq3E6fPu3SNgUEBECv19ttREREdzK3zhmYNWsWJk6cWGeZjh07OnVuo9EIACgoKEBISIhtf0FBASIjI506JxER1YxZC72bW4OB1q1bo3Xr1g1y7oiICBiNRmRkZNg+/M1mM7KysoTuSCAiovrJkGFVYcSfswbcw2vmDOTl5SE7Oxt5eXmwWCzIzs5GdnY2bt78eeZ3t27d8PHHHwMAJEnCzJkzsXjxYmzfvh1ff/01xo8fj9DQUCQkJLjpKoiIiDyP19xamJSUhHXr1tke33vvvQCA3bt3Y+DAgQCAM2fOoKioyFZm9uzZKC4uxtSpU1FYWIh+/fohLS0NOp3OpW0nIrrTcZjAu0ky80XWyWw2w2AwoJFfsMP3/Yreaw84cd+yE3VA8L5lreQvXIVWI3b/vFYrXodGEoth/bVNhco30gQKlQcAnWQQKu8HrjPgcB2C6wyUWYrF67CILURWIbjWBSD+vtCQ6wzIshXlFQUoKipSPEm66j1yVMv/h0aCf/81KbeW4qNrr6rSNnKc1wwTEBERUcPwmmECIiLyXLIsQ42OZnZWuweDASIiUkytBYO46JB7cJiAiIjIx7FnwEEik3+cSWIimnjImURForfvypJ4HVbBdjmTVEajEXvZirapQiPepnKN2IQyrTOJigQnlEnQCtchC2anssjlYuWd+H1XyGVC5UV/35XHiF2HM0QnBEpe9l3NqtI6A2qcg8QxGCAiIsVkqHRrofJTkBO8K/QkIiIi1bFngIiIFOMwgXdjMEBERIrJsjpd/Lyz0D04TEBEROTj2DNARESKcZjAuzEYICIixayySsEAxwncgsMEREREPo49A0REpJj8039qnIdcj8EAEREpJgNOrL1a83nI9ThMQERE5OPYM9AQBPIY2A5xRW4CQSL5GGzHCOYzsEriL0FJFjvGYhVb275CMM8AAGg1YrkGtBp/4To0UiOh8qJr4TtD9DXiTA4A0d+fM/kPRK/Dmb8N0d+HMzlO3Hlu3k3g3RgMEBGRYrKs0pwB3k3gFhwmICIi8nHsGSAiIsU4TODd2DNARETk49gzQEREirFnwLsxGCAiIsXkn8IBNc5DrsdhAiIiIh/HYICIiBSz2voGlG+iVq9ejfDwcOh0OkRHR+PgwYN1lt+yZQu6desGnU6Hnj17YufOnXY/l2UZSUlJCAkJQWBgIOLi4nD27Fnbz3NzczF58mREREQgMDAQnTp1woIFC1BWJrYmhidhMEBERIq5KxjYvHkzEhMTsWDBAhw9ehS9evVCfHw8Ll++XGP5AwcO4JFHHsHkyZNx7NgxJCQkICEhAd98842tzJIlS7By5UqkpKQgKysLTZo0QXx8PEpKSgAAp0+fhtVqxTvvvIOTJ09i+fLlSElJwbx585x/At1MkrnCQ53MZjMMBgP8tK0dXkHMmdUBrbJYROmKFQidiRVFV1mTnFmBUPAYjWB5Z1YH5AqEjvHUFQitVsGVM534+3PF78NRsmxFeUUBioqKoNfrFZ2r6j0yWv80/CSxv4OaVMilyDK/7XDboqOjcd9992HVqlUAAKvVirCwMDz33HOYM2dOtfJjx45FcXExduzYYdt3//33IzIyEikpKZBlGaGhoZg1axb+8pe/AACKiooQHByM1NRUjBs3rsZ2LF26FG+//TYuXLjgzGW7nee8OomIyGtZVfwPqAwyfrmVllYP8srKynDkyBHExcXZ9mk0GsTFxSEzM7PGdmZmZtqVB4D4+Hhb+ZycHJhMJrsyBoMB0dHRtZ4TqAwYWrRo4fgT5mEYDBARkWKyJEOWrCpslZ3VYWFhMBgMti05OblanVevXoXFYkFwcLDd/uDgYJhMphrbaTKZ6ixf9a/IOc+dO4c333wTTz31lAPPlGfymlsLX3nlFXz22WfIzs6Gv78/CgsL6z1m4sSJWLdund2++Ph4pKWlCdcvcruLM12I4omKLMJ1iJKcuMVHFowvnUv4IvhcCQ4TWKwlQuUBQCOJdfu7YphAo/G8P29nft+i3f7O/f0JJipy5vY3wQFZZ4YVHL0OZ34Prpafn283TBAQoHwIoiFcvHgRQ4cOxZgxYzBlyhR3N8dpnvduUYuysjKMGTMGMTExePfddx0+bujQoXj//fdtjz31BUVE5M1klRYdqkp2pNfr650z0KpVK2i1WhQUFNjtLygogNForPEYo9FYZ/mqfwsKChASEmJXJjIy0u64S5cuYdCgQYiNjcWaNWvqvzgP5jXDBIsWLcLzzz+Pnj17Ch0XEBAAo9Fo25o3b95ALSQi8l1qzxlwhL+/P6KiopCRkfFzO6xWZGRkICYmpsZjYmJi7MoDQHp6uq18REQEjEajXRmz2YysrCy7c168eBEDBw5EVFQU3n//fWg0XvNxWiOv6Rlw1p49e9CmTRs0b94cv/vd77B48WK0bNmy1vKlpaV2E1XMZrMrmklERE5ITEzEhAkT0KdPH/Tt2xcrVqxAcXExJk2aBAAYP3482rZta5tzMGPGDAwYMADLli3D8OHDsWnTJhw+fNj2zV6SJMycOROLFy9Gly5dEBERgZdeegmhoaFISEgA8HMg0KFDB7z22mu4cuWKrT219Uh4ujs6GBg6dChGjRqFiIgInD9/HvPmzcOwYcOQmZkJrVZb4zHJyclYtGiRi1tKROTd3LUc8dixY3HlyhUkJSXBZDIhMjISaWlptgmAeXl5dt/aY2NjsXHjRsyfPx/z5s1Dly5dsG3bNvTo0cNWZvbs2SguLsbUqVNRWFiIfv36IS0tDTqdDkBlT8K5c+dw7tw5tGvXzr79Xnq3vlvXGZgzZw5effXVOsucOnUK3bp1sz1OTU3FzJkzHZpA+GsXLlxAp06dsGvXLgwePLjGMjX1DISFhUGrbenwhB6r4H3RACALrzPgigmEkhMHiXaVeeZaBqI4gdAxvjyBUBJ8rTf0BMIKyxVV1xm4J2gCtIJ/BzWxyGU4UbhOlbaR49z6bjFr1ixMnDixzjIdO3ZUrb6OHTuiVatWOHfuXK3BQEBAACcZEhGRT3FrMNC6dWu0bt3aZfX997//xbVr1+xmiBIRkXJWWJ26Hbmm85Drec30x7y8PGRnZyMvLw8WiwXZ2dnIzs7GzZs3bWW6deuGjz/+GABw8+ZNvPDCC/jPf/6D3NxcZGRk4KGHHkLnzp0RHx/vrssgIrojueNuAlKP5w0q1iIpKcluAaF7770XALB7924MHDgQAHDmzBkUFRUBALRaLU6cOIF169ahsLAQoaGhGDJkCF5++WUOAxAREf2C1wQDqampSE1NrbPML+dCBgYG4t///ncDt4qIiAD33U1A6vCaYICIiDyXFRZIUH6Xk1WFc5A4BgMOqryVz9HbdpxIL+wFa4U3DCfyHwjeDCt626YzZMF8CaIpqwHxWyQliwtu2xTNReHU71vwtj9nbi10ybdR0eeq4aZ08ds3/RqDASIiUkz+KTuBGuch12MwQEREilklKySJtxZ6K6+5tZCIiIgaBnsGiIhIscoJhMq/X3ICoXswGCAiIhWoc2uhM5OKSTkOExAREfk49gwQEZFiVtkCNb5fWl2QkZWqYzBARESKcQVC78ZhAiIiIh/HngEiIlJMhkWVVRNl3k3gFgwGiIhIscrFgrjokLfiMAEREZGPY8+AgyqTn0gNV4FoghhXJDYSTIzjVBWC1+0KoolxKo8RS47jTB2AWHIj555bwdehC35/4s9Vw/9tOPP7E36u5AZMVCSa7cuRczI3gVdjMEBERIrJsgWyCl+YZN5a6Bae97WMiIiIXIo9A0REpBgnEHo3BgNERKRY5a2FKgwT8NZCt+AwARERkY9jzwARESkmyyotR+yKO6WoGgYDRESkGOcMeDcOExAREfk49gwQEZFiXGfAuzEYICIixbgCoXdjMFCPqmU7RZbvdG6pT9FjXPAH0wBLltZQiQvqEHPH/P6cqsPzrkP89+GJbQI86bl15n2N7mwMBupx48aNn/6v2DWfjR5FbC18Z/jec+o6rnhu+ftznCc+Vzdu3IDBYFDlXJV3E6gxTMAJhO7AYKAeoaGhyM/PR7NmzSBJDZioyIuZzWaEhYUhPz8fer3e3c3xWnweleNz6BhZlnHjxg2EhoaqeFaLSn0ZnDPgDgwG6qHRaNCuXTt3N8Mr6PV6vgGrgM+jcnwO66dWjwDdGRgMEBGRYpXd+xwm8FYMBoiISDEGA96Niw6RYgEBAViwYAECAgLc3RSvxudROT6HRM6RZN5bQkRETjKbzTAYDGjWuBskSav4fLJswY1bp1FUVMR5Hy7EYQIiIlKMwwTejcMEREREPo49A0REpJhaOQWYm8A9GAwQEZFilTkFmJvAW3GYgFT3yiuvIDY2Fo0bN0ZQUJC7m+MVVq9ejfDwcOh0OkRHR+PgwYPubpLX2bt3L0aOHInQ0FBIkoRt27a5u0lEXoPBAKmurKwMY8aMwdNPP+3upniFzZs3IzExEQsWLMDRo0fRq1cvxMfH4/Lly+5umlcpLi5Gr169sHr1anc3xSfJslW1jVyPtxZSg0lNTcXMmTNRWFjo7qZ4tOjoaNx3331YtWoVAMBqtSIsLAzPPfcc5syZ4+bWeSdJkvDxxx8jISHB3U2541XdWqjzbw9JUv79UpatKCnL462FLsaeASI3Kisrw5EjRxAXF2fbp9FoEBcXh8zMTDe2jIh8CScQErnR1atXYbFYEBwcbLc/ODgYp0+fdlOriMTJKkweVPM8JIY9A+SQOXPmQJKkOjd+eBH5LnfOGRCdgLtlyxZ069YNOp0OPXv2xM6dO391LTKSkpIQEhKCwMBAxMXF4ezZs3Zlrl+/jsceewx6vR5BQUGYPHkybt68Kdx2T8FggBwya9YsnDp1qs6tY8eO7m6m12nVqhW0Wi0KCgrs9hcUFMBoNLqpVUTeQ3QC7oEDB/DII49g8uTJOHbsGBISEpCQkIBvvvnGVmbJkiVYuXIlUlJSkJWVhSZNmiA+Ph4lJSW2Mo899hhOnjyJ9PR07NixA3v37sXUqVMb/HobjEzUQN5//33ZYDC4uxker2/fvvL06dNtjy0Wi9y2bVs5OTnZja3ybgDkjz/+2N3N8AlFRUUyALmRn1H2bxSqeGvkZ5QByEVFRQ7V37dvX/nZZ5+1PbZYLHJoaGitfz8PP/ywPHz4cLt90dHR8lNPPSXLsixbrVbZaDTKS5cutf28sLBQDggIkP/v//5PlmVZ/vbbb2UA8qFDh2xl/vWvf8mSJMkXL1507InzMOwZINXl5eUhOzsbeXl5sFgsyM7ORnZ2tld3oTWkxMRErF27FuvWrcOpU6fw9NNPo7i4GJMmTXJ307zKzZs3ba81AMjJybG9DskVrCpulXcp/HIrLS2tVqMzE3AzMzPtygNAfHy8rXxOTg5MJpNdGYPBgOjoaFuZzMxMBAUFoU+fPrYycXFx0Gg0yMrKcujZ8jScQEiqS0pKwrp162yP7733XgDA7t27MXDgQDe1ynONHTsWV65cQVJSEkwmEyIjI5GWllZtUiHV7fDhwxg0aJDtcWJiIgBgwoQJSE1NdVOryFlhYWF2jxcsWICFCxfa7XNmAq7JZKqxvMlksv28al9dZdq0aWP3cz8/P7Ro0cJWxtswGCDVpaam8s1X0PTp0zF9+nR3N8OrDRw4EDKXTXEbtRYLqjpPfn6+3ToDAQEBqpyfasZggIiIFFP71kK9Xl/vokPOTMA1Go11lq/6t6CgACEhIXZlIiMjbWV+PUGxoqIC169f99qJv5wzQEREXsnf3x9RUVHIyMiw7bNarcjIyEBMTEyNx8TExNiVB4D09HRb+YiICBiNRrsyZrMZWVlZtjIxMTEoLCzEkSNHbGW++OILWK1WREdHq3Z9rsSeASIiUqxyiEaFrIWCQz2JiYmYMGEC+vTpg759+2LFihV2E3DHjx+Ptm3bIjk5GQAwY8YMDBgwAMuWLcPw4cOxadMmHD58GGvWrAFQuZT1zJkzsXjxYnTp0gURERF46aWXEBoaalve+q677sLQoUMxZcoUpKSkoLy8HNOnT8e4ceMQGhqq+DlwBwYDRESkAgsASYXziAUD9U3AzcvLg0bzcyd4bGwsNm7ciPnz52PevHno0qULtm3bhh49etjKzJ49G8XFxZg6dSoKCwvRr18/pKWlQafT2cp8+OGHmD59OgYPHgyNRoPRo0dj5cqVCq/dfZioiIiInFaVqEij0UOSlAcDsizDajUzUZGLsWeAiIgUq7wLQJ1ggFyPwQAREalAnWBAdJiA1MG7CYiIiHwcgwEiN8vNzbVlfqy6j9kRAwcOtB1XtQQvkdvIVvU2cjkGA0QeYteuXXb3Ni9cuBATJ060PZ44caLdcqwfffRRvalaiVxFVvE/cj3OGSDyEC1btkTLli0dLt+iRQuYzeYGbBER+Qr2DBCp6MqVKzAajfjb3/5m23fgwAH4+/tXW/WM6M6ibtZCci32DBCpqHXr1njvvfeQkJCAIUOGoGvXrnjiiSdsi5MQ3blklW4E4DCBOzAYIFLZgw8+iClTpuCxxx5Dnz590KRJE9tSqCJ+na6VmSCJqKEwGCBqAK+99hp69OiBLVu24MiRI0y/Sj6Ak/+8GecMEDWA8+fP49KlS7BarcjNzXV3c4gajL+//09pey2qbUajEf7+/i6+Et/GngEilZWVleHxxx/H2LFj0bVrV/zpT3/C119/jTZt2ri7aUSq0+l0yMnJQVlZmWrn9Pf3t0sKRA2PwQCRyl588UUUFRVh5cqVaNq0KXbu3Iknn3wSO3bscHfTiBqETqfjh7eX4zABkYr27NmDFStWYMOGDdDr9dBoNNiwYQP27duHt99+293NIyKqEXsGiFQ0cOBAlJeX2+0LDw9HUVGRm1pERFQ/9gwQeYjY2FjExsY6XH7YsGHo3r17A7aIiHyFJDN5NJFbVVRU2O44CAgIQFhYmEPHXbx4Ebdv3wYAtG/fnrOvichpDAaIiIh8HIcJiIiIfByDASIiIh/HYICIiMjHMRggIiLycQwGiIiIfByDASIiIh/HYICIiMjHMRggIiLycf8f5hVFtEiNm8IAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(5, 5))\n", + "plot = ax.imshow(I_reshape, cmap=\"magma\", extent=[-r_out_rad, r_out_rad, -r_out_rad, r_out_rad])\n", + "cmap = plt.colorbar(plot)\n", + "cmap.set_label(r'I [Jy $sr^{-1}$]', size = 15)\n", + "ax.set_title(r'$N^{2}$ points, with N = ' + str(N))\n", + "ax.set_xlabel(\"x ['']\")\n", + "ax.set_ylabel(\"y ['']\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "PROD_Frank2DVenv", + "language": "python", + "name": "prod_frank2dvenv" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + }, + "nbTranslate": { + "displayLangs": [ + "*" + ], + "hotkey": "alt-t", + "langInMainMenu": true, + "sourceLang": "en", + "targetLang": "fr", + "useGoogleTranslate": true + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/frank/fourier2d.py b/frank/fourier2d.py index 54d6e6a0..b0c73308 100644 --- a/frank/fourier2d.py +++ b/frank/fourier2d.py @@ -1,5 +1,6 @@ import numpy as np +import time class DiscreteFourierTransform2D(object): def __init__(self, Rmax, N, nu=0): @@ -35,7 +36,8 @@ def __init__(self, Rmax, N, nu=0): def get_collocation_points(self): return np.array([self.Xn, self.Yn]), np.array([self.Un, self.Vn]) - def coefficients(self, u = None, v = None, direction="forward"): + def coefficients(self, u = None, v = None, x = None, y = None, direction="forward"): + #start_time = time.time() if direction == 'forward': ## Normalization is dx*dy since we the DFT to be an approximation ## of the integral (which depends on the area) @@ -51,16 +53,16 @@ def coefficients(self, u = None, v = None, direction="forward"): norm = 1 / (4*self.Xmax*self.Ymax) factor = 2j*np.pi - X, Y = self.Un, self.Vn + X, Y = u, v if u is None: - u = self.Xn - v = self.Yn + X, Y = self.Un, self.Vn + u = self.Xn + v = self.Yn else: raise AttributeError("direction must be one of {}" "".format(['forward', 'backward'])) - H = norm * np.exp(factor*(np.outer(u, X) + np.outer(v, Y))) - + #print("--- %s minutes to calculate 2D-DFT coefficients---" % (time.time()/60 - start_time/60)) return H diff --git a/frank/statistical_models.py b/frank/statistical_models.py index b96218cf..8d4ce21f 100644 --- a/frank/statistical_models.py +++ b/frank/statistical_models.py @@ -25,6 +25,7 @@ from frank.constants import rad_to_arcsec, deg_to_rad from frank.minimizer import LineSearch, MinimizeNewton +import time class VisibilityMapping: r"""Builds the mapping between the visibility and image planes. @@ -206,16 +207,40 @@ def map_visibilities(self, u, v, V, weights, frequencies=None, geometry=None): ks = ki[start:end] ws = wi[start:end] Vs = Vi[start:end] - + X = self._get_mapping_coefficients(qs, ks, us, vs) - - wXT = np.array(X.T * ws, order='C') - - Ms[i] += np.dot(wXT, X) - js[i] += np.dot(wXT, Vs) + X_CT = self._get_mapping_coefficients(qs, ks, us, vs, inverse=True) + wXT = np.matmul(X_CT, np.diag(ws)) # this line is the same as below. + #wXT = np.matmul(np.transpose(np.conjugate(X)), np.diag(ws), dtype = "complex64") + Ms[i] += np.matmul(wXT, X, dtype="complex128") + js[i] += np.matmul(wXT, Vs, dtype="complex128") start = end end = min(Ndata, end + Nstep) + + print("M type", Ms[0].dtype) + print("j type", js[0].dtype) + print("M imag-> ", " max: ", np.max(Ms[0].imag), ", min: ", np.min(Ms[0].imag) , ", mean: ", np.mean(Ms[0].imag) ,", median: ", np.median(Ms[0].imag), ", std: ", np.std(Ms[0].imag)) + print("M real-> ", " max: ", np.max(Ms[0].real), ", min: ", np.min(Ms[0].real) , ", mean: ", np.mean(Ms[0].real) ,", median: ", np.median(Ms[0].real), ", std: ", np.std(Ms[0].real)) + tol = 1e-10 + print("with tol: ", tol, ", M is symmetric?", np.allclose(Ms[0].real, Ms[0].T.real, atol=tol)) + # Ms[0] = np.abs(Ms[0]) + # from scipy.linalg import issymmetric + # print(issymmetric(Ms[0]), "that M is a Symmetric Matrix") + + + + import matplotlib.pyplot as plt + plt.matshow(Ms[0].real, cmap="magma", vmax = np.max(Ms[0].real), vmin = np.mean(Ms[0].real)) + plt.colorbar() + plt.title("M matrix, real part") + plt.show() + #import sys + #sys.exit() + + + #Ms[0] = np.loadtxt(r'.\..\Notebooks\M_N75.txt', dtype = 'c8') + #js[0] = np.loadtxt(r'.\..\Notebooks\j_N75.txt', dtype = 'c8') # Compute likelihood normalization H_0, i.e., the # log-likelihood of a source with I=0. @@ -263,10 +288,10 @@ def check_hash(self, hash, multi_freq=False, geometry=None): self._DFT.Rmax == hash[1].Rmax and self._DFT.size == hash[1].size and self._DFT.order == hash[1].order and - #geometry.inc == hash[2].inc and - #geometry.PA == hash[2].PA and - #geometry.dRA == hash[2].dRA and - #geometry.dDec == hash[2].dDec and + geometry.inc == hash[2].inc and + geometry.PA == hash[2].PA and + geometry.dRA == hash[2].dRA and + geometry.dDec == hash[2].dDec and self._vis_model == hash[3] ) @@ -468,7 +493,6 @@ def interpolate(self, f, r, space='Real'): def _get_mapping_coefficients(self, qs, ks, u, v, geometry=None, inverse=False): """Get :math:`H(q)`, such that :math:`V(q) = H(q) I_\nu`""" - """ if self._vis_model == 'opt_thick': # Optically thick & geometrically thin if geometry is None: @@ -482,7 +506,6 @@ def _get_mapping_coefficients(self, qs, ks, u, v, geometry=None, inverse=False): scale = np.exp(-np.outer(ks*ks, self._H2)) else: raise ValueError("model not supported. Should never occur.") - """ if inverse: scale = np.atleast_1d(1/scale).reshape(1,-1) qs = qs / rad_to_arcsec @@ -490,8 +513,7 @@ def _get_mapping_coefficients(self, qs, ks, u, v, geometry=None, inverse=False): else: direction='forward' - #H = self._DHT.coefficients(qs, direction=direction) * scale - H = self._DFT.coefficients(u, v, direction=direction) #* scale + H = self._DFT.coefficients(u, v, direction=direction)*scale return H @@ -725,15 +747,23 @@ def __init__(self, DHT, M, j, p=None, scale=None, guess=None, self.u, self.v = self._DFT.uv_points self.Ykm = self._DFT.coefficients(direction="backward") self.Ykm_f = self._DFT.coefficients(direction="forward") - m, c , l = -5, 60, 1e4 + + m, c , l = -5, 60, 1e4 + #m, c, l = 0.23651345032212925, 60.28747193555951, 1.000389e+05 + #start_time = time.time() #m, c, l = self.minimizeS() # Finding best parameters to S matrix. - S_real = self.calculate_S(self.u, self.v, l, m, c) + #print("--- %s minutes to minimize S---" % (time.time()/60 - start_time/60)) + S_real = self.calculate_S_real(self.u, self.v, l, m, c) + start_time = time.time() S_real_inv = np.linalg.inv(S_real) + print("--- %s minutes to calculate S_real_inv---" % (time.time()/60 - start_time/60)) self._Sinv = S_real_inv - + print(self._Sinv.dtype, " Sinv dtype") + start_time = time.time() self._fit() + print("--- %s minutes to fit---" % (time.time()/60 - start_time/60)) - def true_squared_exponential_kernel(self, u, v, l, m, c): + def true_squared_exponential_kernel(self, u, v, l, m, c): u1, u2 = np.meshgrid(u, u) v1, v2 = np.meshgrid(v, v) q1 = np.hypot(u1, v1) @@ -748,75 +778,106 @@ def power_spectrum(q, m, c): for i in indexes: q[i] = q_min return np.exp(m*np.log(q) + c) - + p1 = power_spectrum(q1, m, c) p2 = power_spectrum(q2, m, c) - SE_Kernel = np.sqrt(p1 * p2) * np.exp(-0.5*((u1-u2)**2 + (v1-v2)**2)/ l**2) return SE_Kernel - - def calculate_S(self, u, v, l, m, c): + + def calculate_S_real(self, u, v, l, m, c): + start_time = time.time() S_fspace = self.true_squared_exponential_kernel(u, v, l, m, c) - S_real = np.matmul(self.Ykm, np.matmul(S_fspace, self.Ykm_f)) + print("--- %s minutes to calculate S---" % (time.time()/60 - start_time/60)) + start_time = time.time() + S_real = np.matmul(self.Ykm, np.matmul(S_fspace, self.Ykm_f), dtype = "complex64") + print("--- %s minutes to calculate S_real---" % (time.time()/60 - start_time/60)) + + print(S_real.dtype, " S_real dtype") + import matplotlib.pyplot as plt + plt.matshow(S_fspace, cmap="magma", vmin = 0, vmax = 1e-3) + plt.colorbar() + plt.title("S real matrix, real part ") + plt.show() + return S_real - + + def calculate_mu_cholesky(self, Dinv): + print("calculate mu with cholesky") + try: + Dchol = scipy.linalg.cho_factor(Dinv) + mu = scipy.linalg.cho_solve(Dchol, self._j) + + except np.linalg.LinAlgError: + U, s, V = scipy.linalg.svd(Dinv, full_matrices=False) + s1 = np.where(s > 0, 1. / s, 0) + mu = np.dot(V.T, np.multiply(np.dot(U.T, self._j), s1)) + return mu + + def calculate_mu_gc(self, Dinv): + from scipy.sparse.linalg import cg, bicg, bicgstab, gmres + method = "BiConjugate Gradient Method" + #from scipy.sparse import csr_matrix, issparse + #is_sparse = issparse(Dinv) + #print("Is Dinv sparse?: ", is_sparse) + start_time = time.time() + mu, exitCode = bicgstab(Dinv, self._j, atol = 0) + print("--- %s minutes to calculate mu---" % (time.time()/60 - start_time/60)) + print("Succesful ", method, "?: ", exitCode == 0) + print("Is the result correct?: ", np.allclose(np.dot(Dinv, mu), self._j)) + return mu + + def likelihood(self, param, data): + from scipy.special import gamma + m, c, l = param + Wvalues = self._Wvalues + N = np.diag(1/Wvalues) + + alpha = 1.3 + l0 = 1e7 + + # Create an Inverse Gamma distribution function + def inv_gamma_function(l, alpha, beta): + return ((gamma(alpha)*beta)**(-1))*((beta/l)**(alpha + 1))*np.exp(-beta/l) + + S_real = self.calculate_S_real(self.u, self.v, l, m, c) + start_time = time.time() + S_real_inv = np.linalg.inv(S_real) + print("--- %s minutes to calculate S_real_inv ---" % (time.time()/60 - start_time/60)) + Dinv = self._M + S_real_inv + start_time = time.time() + D = np.linalg.inv(Dinv) + print("--- %s minutes to calculate D ---" % (time.time()/60 - start_time/60)) + mu = self.calculate_mu_gc(Dinv) + + start_time = time.time() + logdetS = np.linalg.slogdet(S_real)[1] + logdetD = np.linalg.slogdet(D)[1] + logdetN = np.linalg.slogdet(N)[1] + print("--- %s minutes to calculate determinants ---" % (time.time()/60 - start_time/60)) + factor = np.log(2*np.pi) + + start_time = time.time() + log_likelihood = 2*np.log(np.abs((1/m)*(1/c))) \ + + 2*np.log(inv_gamma_function(l, alpha, l0)) \ + - 0.5*(factor + logdetN) \ + - 0.5*(factor + logdetS) \ + + 0.5*(factor + logdetD) \ + + 0.5*np.dot(np.transpose(self._j), mu) \ + - 0.5*np.dot(np.transpose(np.conjugate(data)), np.dot(np.diag(Wvalues), data)) + print("--- %s minutes to calculate log_likelihood ---" % (time.time()/60 - start_time/60)) + return -log_likelihood + def minimizeS(self): from scipy.optimize import minimize - from scipy.special import gamma V = self._V - - def calculate_D(S): - S_real_inv = np.linalg.inv(S) - Dinv = self._M + S_real_inv - D = np.linalg.inv(Dinv) - return [Dinv, D] - - def calculate_mu(Dinv): - try: - Dchol = scipy.linalg.cho_factor(Dinv) - mu = scipy.linalg.cho_solve(Dchol, self._j) - - except np.linalg.LinAlgError: - U, s, V = scipy.linalg.svd(Dinv, full_matrices=False) - s1 = np.where(s > 0, 1. / s, 0) - mu = np.dot(V.T, np.multiply(np.dot(U.T, self._j), s1)) - return mu - - def likelihood(param, data): - m, c, l = param - Wvalues = self._Wvalues - N = np.diag(1/Wvalues) - - alpha = 1.3 - l0 = 1e7 - - # Create an Inverse Gamma distribution function - def inv_gamma_function(l, alpha, beta): - return ((gamma(alpha)*beta)**(-1))*((beta/l)**(alpha + 1))*np.exp(-beta/l) - - S = self.calculate_S(self.u, self.v, l, m, c) - [Dinv, D] = calculate_D(S) - mu = calculate_mu(Dinv) - logdetS = np.linalg.slogdet(S)[1] - logdetD = np.linalg.slogdet(D)[1] - logdetN = np.linalg.slogdet(N)[1] - factor = np.log(2*np.pi) - - log_likelihood = 2*np.log(np.abs((1/m)*(1/c))) \ - + 2*np.log(inv_gamma_function(l, alpha, l0)) \ - - 0.5*(factor + logdetN) \ - - 0.5*(factor + logdetS) \ - + 0.5*(factor + logdetD) \ - + 0.5*np.dot(np.transpose(self._j), mu) \ - - 0.5*np.dot(np.transpose(np.conjugate(data)), np.dot(np.diag(Wvalues), data)) - return -log_likelihood - - result = minimize(likelihood, x0=np.array([-5, 60, 1e5]), args=(V,), - bounds=[(-6, 6), (1, 70), (1e4, 1e6)], - method="Nelder-Mead", tol=1e-7, - ) + print("minimizing") + start_time = time.time() + result = minimize(self.likelihood, x0=np.array([-5, 60, 1e5]), args=(V,), + method="Nelder-Mead", tol=1e-7, + bounds=[(-7, 7), (1, 80), (1e4, 1e5)]) m, c, l = result.x - print("Best parameters: ", "m: ", m, "c: ", c, "l: ", "{:e}".format(l)) + print("--- %s minutes to minimizing ---" % (time.time()/60 - start_time/60)) + print("Result: ", "m: ", m, "c: ", c, "l: ", "{:e}".format(l)) return [m, c, l] def _fit(self): @@ -828,21 +889,19 @@ def _fit(self): Dinv = self._M + Sinv - try: - self._Dchol = scipy.linalg.cho_factor(Dinv) - self._Dsvd = None - - self._mu = scipy.linalg.cho_solve(self._Dchol, self._j) - - except np.linalg.LinAlgError: - U, s, V = scipy.linalg.svd(Dinv, full_matrices=False) - - s1 = np.where(s > 0, 1. / s, 0) - - self._Dchol = None - self._Dsvd = U, s1, V + """ + #import scipy.linalg as sc + def is_pos_def(x): + return np.all(np.linalg.eigvals(x) > 0) + a = sc.issymmetric(Dinv) # necessary condition to be SPD + b = is_pos_def(Dinv) # necessary condition to be SPD + # there is left one condition necessary to be SPD, which is that xT * A * x > 0 for all x != 0. + print("Is symmetric: ", a) + print("Is positive definite: ", b) + """ - self._mu = np.dot(V.T, np.multiply(np.dot(U.T, self._j), s1)) + #self._mu = self.calculate_mu_cholesky(Dinv) + self._mu = self.calculate_mu_gc(Dinv) # Reset the covariance matrix - we will compute it when needed if self._Nfields > 1: