sersic_fit_test.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "source": [
        "%load_ext autoreload\n",
        "%autoreload 2"
      ],
      "outputs": [],
      "execution_count": 1,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.signal import convolve2d\n",
        "from scipy.optimize import minimize\n",
        "from scipy.special import gamma\n",
        "import scipy.stats as st\n",
        "from sklearn.metrics import mean_squared_error\n",
        "from gzbuilder_analysis.parsing import sanitize_param_dict\n",
        "import gzbuilder_analysis.rendering as rendering\n",
        "from gzbuilder_analysis.rendering.sersic import sersic2d, oversampled_sersic_component, _b\n",
        "from gzbuilder_analysis.fitting.jupyter import UpdatableDisplay\n",
        "import lib.galaxy_utilities as gu\n",
        "from tqdm import tqdm\n",
        "from numba import jit, prange, vectorize"
      ],
      "outputs": [],
      "execution_count": 2,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "sid_list = np.loadtxt('lib/subject-id-list.csv', dtype='u8')\n",
        "simard_df = pd.read_csv('lib/simard-catalog.csv', index_col=0)\n",
        "best_models = pd.read_pickle('lib/best_individual.pickle')"
      ],
      "outputs": [],
      "execution_count": 3,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "PARAMS = np.array(('roll', 'rEff', 'axRatio', 'i0'))\n",
        "def comp_to_p(d):\n",
        "  mu = np.array(d['mu'])\n",
        "  return np.concatenate((mu, [d[k] for k in PARAMS]))\n",
        "\n",
        "def comp_from_p(p):\n",
        "  return dict(mu=p[:2], **{k: v for k, v in zip(PARAMS, p[2:])})"
      ],
      "outputs": [],
      "execution_count": 4,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "s0 = 'float32('+('float32,'*9)+')'\n",
        "s1 = 'float64('+('float64,'*9)+')'\n",
        "@vectorize([s0, s1], target='parallel')\n",
        "def _sersic_func(x, y, mux, muy, roll, rEff, axRatio, i0, n):\n",
        "  out = 0.0\n",
        "  ds = [-0.4, -0.2, 0.0, 0.2, 0.4]\n",
        "  for dx in ds:\n",
        "    for dy in ds:\n",
        "      out += sersic2d(\n",
        "        x + dx, y + dy,\n",
        "        mu=(mux, muy), roll=roll, rEff=rEff,\n",
        "        axRatio=axRatio, i0=i0, c=2, \n",
        "        n=n\n",
        "      ) / 25\n",
        "  return out\n"
      ],
      "outputs": [],
      "execution_count": 5,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def bulge_disk_model(p, cx, cy):\n",
        "    disk = _sersic_func(cx, cy, *p[:6], 1)\n",
        "    bulge = _sersic_func(cx, cy, *p[6:12], 4)\n",
        "#     disk = sersic2d(cx, cy, mu=p[:2], roll=p[2], rEff=p[3],\n",
        "#                     axRatio=p[4], c=2, i0=p[5], n=1.0)\n",
        "#     bulge = sersic2d(cx, cy, mu=p[6:8], roll=p[8], rEff=p[9],\n",
        "#                      axRatio=p[10], c=2, i0=p[11], n=4.0)\n",
        "    return disk + bulge"
      ],
      "outputs": [],
      "execution_count": 25,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def f(p, target, cx, cy, psf=np.ones((1,1)), pixel_mask=np.ones(1),\n",
        "      loss=mean_squared_error):\n",
        "  im = bulge_disk_model(p, cx, cy)\n",
        "  im_psf = convolve2d(\n",
        "    im, psf, mode='same', boundary='symm'\n",
        "  )\n",
        "  return loss(target * pixel_mask, im_psf * pixel_mask)"
      ],
      "outputs": [],
      "execution_count": 26,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "display(_m)\n",
        "display(p0)"
      ],
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "{'disk': {'mu': array([32.53491211, 32.5390625 ]),\n",
              "  'roll': -0.025635408521677425,\n",
              "  'rEff': 25.90694972290857,\n",
              "  'axRatio': 0.6071298668687596,\n",
              "  'i0': 0.12,\n",
              "  'c': 2,\n",
              "  'n': 1},\n",
              " 'bulge': {'mu': array([35.05834961, 33.0703125 ]),\n",
              "  'roll': 0.08314123188844122,\n",
              "  'rEff': 5.8334124206994575,\n",
              "  'axRatio': 0.48698760942986274,\n",
              "  'i0': 0.71,\n",
              "  'c': 2,\n",
              "  'n': 0.5},\n",
              " 'bar': None,\n",
              " 'spiral': array([], dtype=float64)}"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "array([ 3.25349121e+01,  3.25390625e+01, -2.56354085e-02,  2.59069497e+01,\n",
              "        6.07129867e-01,  1.20000000e-01,  3.50583496e+01,  3.30703125e+01,\n",
              "        8.31412319e-02,  5.83341242e+00,  4.86987609e-01,  7.10000000e-01])"
            ]
          },
          "metadata": {}
        }
      ],
      "execution_count": 29,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.imshow(oversampled_sersic_component(m['disk'], image_size))"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 32,
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x12c6b6fd0>"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": [
              "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD7CAYAAACscuKmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAde0lEQVR4nO2df8xkVXnHv8/M++4uu4DLot1sWOLSSCT8IYvZIERiFIqh1ug/xmhNQxsS/rENpiYKbdLUpE3qPyp/NCbEX/xh/VHUSohR6YoxTZqVpaACK4IWZenCYgsiLLv7vjNP/5i77zznzH3ue+beuXdm9nw/yZv3zr333HPmzpy53+c5zzmPqCoIIWc/vXk3gBDSDezshGQCOzshmcDOTkgmsLMTkgns7IRkQqPOLiI3isjjIvKkiNw2q0YRQmaP1B1nF5E+gF8AuAHAUQAPAPigqj42u+YRQmbFSoOyVwF4UlV/BQAi8lUA7wXgdvYtslW3YUeDKgkhVZzEKzitp6TsWJPOfhGAp83rowDeUlVgG3bgLXJ9gyoJIVUc0oPusSadPQkRuQXALQCwDdvbro4Q4tDEQfcMgIvN673FvgBVvVNVD6jqgVVsbVAdIaQJTTr7AwAuFZFLRGQLgA8AuGc2zSKEzJraMl5V10XkLwF8D0AfwBdU9dGZtYwQMlMa2eyq+h0A35lRWwghLcIIOkIygZ2dkExgZyckE9jZCckEdnZCMoGdnZBMYGcnJBPY2QnJBHZ2QjKh9VlvxCCl04xJFUxiMjP4ZCckE9jZCckEyvhpoAzvnqb3nGbABnyyE5IJ7OyEZEK+Mp6SPA/qfM5nqfTnk52QTGBnJyQT2NkJyYSz02ZfRntcMv/d1eG8WzAm9fuzZLZ95t8wQvKBnZ2QTFheGb9IUj13CT4Lmt7DeZgB3ndwQeX9pndYRL4gIsdF5BGzb5eI3CciTxT/L2i3mYSQpqT8nH4JwI3RvtsAHFTVSwEcLF4TQhaYTWW8qv5IRPZFu98L4O3F9l0Afgjg4zNr1bwl+pxlufQWyERpCR3OWOpO+5m1KfsX1Jtf91u9W1WPFdvPAtg9o/YQQlqi8SNMVRWA+xMlIreIyGERObyGU02rI4TUpK43/jkR2aOqx0RkD4Dj3omqeieAOwHgfNk1/lGYh1TvSJ7PRYYv0ohAgkRueo8amwGp92secr8leV/3G3IPgJuK7ZsAfHs2zSGEtEXK0NtXAPwngDeKyFERuRnAPwG4QUSeAPBHxWtCyAKT4o3/oHPo+lo1diHfZyhpZyrJu5LaXZgRVTJ6lu/TkdF1Ppda0r/qvbQl8VuS9wtk6BFC2oSdnZBMYGcnJBOWeCJMs9+pxrZ4U7t0xna1dD2U2Q/r01kOF1nbus59bmjnJ9v2Xtu6suWnvOd8shOSCezshGTC4sv4ecj1aeusUUdj2d1brN9p990Mp5e0WuetTSv9K6S2951ZOHlf9h2qaOJifWMIIa3Bzk5IJiyGjF9kqZ5w7WRJXkd615H7XUj8VHne75fvr/AkB+84sZ4k6V/Hy2+kd9X3LEnizyMaz8AnOyGZwM5OSCbMT8bXkO6tyfXE67pyvUo2z1Di1/LgtzUppkZbkgNvrHS3ZkBD6Z/s5U+V+wkSf5E8+HyyE5IJ7OyEZEL3Mn5K+Z4k3Rt61ivlsSeva5Sprqfh++x6KazYyZ4gV6VKktryRrq70r/KS+/d/4oyQT3evYzfY4L0nqkHv6Gk55OdkExgZyckE9jZCcmExYigi2hsp09rm9cZOjNlkm3xOja3d+1Fst8B32YP7Exrf0bnWx+AuZZr50f3Jcm2rxjG8+z5JFse8IfrEifc1IrAm9KG55OdkExgZyckE+Ym45Oj4RpGwCVJ98RhtOBaVe23bQ5k/PSS3m1/vD9l8kvTOfR1IuCcsLUJ2R3I4ATpL5EMt5I2ZRgvlu2OxK81XNdVBN6Uw3IpSSIuFpH7ReQxEXlURG4t9jNHOyFLRIqMXwfwUVW9HMDVAD4sIpeDOdoJWSpSMsIcA3Cs2P69iBwBcBHq5GiXhp72tiLgonOS5Lon1eNrO+dNtNlrjzhmQFWdKfubUiXpA7mspftlQsaXS3/1riWRbFV7n2ydRirDuZcT9ZhrV5hHrgffK1IVgZewIm6ypHeYykEnIvsAXAngEJijnZClIrmzi8i5AL4B4COq+pI9VpWjPcjPrszPTsi8SPLGi8gqRh39y6r6zWJ3Uo72ID97b1e5DpllgAzgy+DUCSpeIIzd3+uX74+v59Uft8X12pdva6Xp4MtVl9TJHxsNUP914M0u319Zfuh4wyu8/K7cHw7K9/dC2SzBpZ3vYyzb7XnmmP38K4NypgzEiU3gaRNVpnjjBcDnARxR1U+ZQ8zRTsgSkfJkfyuAPwPwMxF5uNj3NxjlZP96ka/91wDe304TCSGzIMUb/x/wcwDUy9FOCOmcxVyDrulEFsdOdofU4rYkRL1JvyICz0Zj2WMVZbTv2POp5YP3hvLzZhlBNzF0ZqoxxwKbdWANY3/oLbDzrZ1cUV4GxjY3Q2+Bb8Oeo7HPw5l8M7TXir4nCUN01ucwETWYEnU3w4k0jI0nJBPY2QnJhI5lvJTL91lPZJlWusf1O8Nq4g2dxRLOkeSuVI/KBHLR5EEP9se30dZj346VeoGkRxqOOowj4GTgSXdnGG0Qlbdyd1A+3CY9ZxgOCL8DdugvkPfmXgzi8rZtZn9vWHoKEH0eQ8fEsN/FqM1TT6ThGnSEkBTY2QnJhMVclmpa6Z46kcVKosCzHZX3vN7Wy14hyQO57pTRlXgEoPyY9SZ7+0fXLo+u86R77IxOwc4vkdj7a49Z6T4o97LLeiRpB+XSOzAP1m00XBQB50l8+/6tSRC2Hmqlu10iK7hR9qSKqLsEST8q73jqU+fGl82Hr1ph2z9ECDmbYGcnJBPmJ+NTl6VKkO61JrJ4ATJxPVaGr1gvveNlj87ToLwv4wOJHsj9cVuG5pzhSvieA+luJX2//Jw63ngJvNzhaYF0tyrayPWePWc98sbb89bN+zf7xbyvCTPANmjorBuwDhd7O9RKd+vNn8g1b8yNaSV9XL83eWaG8MlOSCawsxOSCd3L+BL5PinDp/wNmpjbPeUc9KrYdk+6rzjnIJbkRtKvmv2rYZnharlEt/sDSb8atjkoY6W7I+MnvPGerHdWgurFQTGD8mO99fGFe2u+jO+tmfe5ZqS7lfdr9rOMvPGBXB83xtYSvMWhb8fY0YDADBvEtov9npUH37iSHnDnw4fnOJ75uP5ZrC5LCDk7YGcnJBPY2QnJhIWMoAtIHW6zePa4N6klHlLxJrI4dvrEMNpqv3Q7GDrbEtvsxubeMj5vsKXcTh/ENvuqs21s+dB+R4h3O+0wmrXZo2Es+9ra1tZO76+hdD8A9E+ZaDJ7+/v2szT7qz5/cywYUrOnxGV6zvKI1n6PI9jsR5gweaZyQM2ZMJO8hl0CfLITkgns7IRkwvwSOzZdIik177k40j9xPnw4qaV8zrlGZkCKdB9sDds8NHLdHhtsGZ+zvtWeE8n4reNtW8ZKejWf9rAfL5GEcgLpbuV5WH/PSPS+SQ/QO12+fyVKIWDvc98OEa45Q0rRZxasxBVeuLR4LIAnMtScIVh6KpbR1t6wc/jLLzUxpJwYXZfExvfc71d8shOSCezshGTCYnjj66wOG5SPI+i8+ezOdtVEGG/5qRVnG9FEFhsZt8Vuh2220n19m5htc47df07Y5IE975yxJBxsHW8PtxrZuBLLeEdGDoxn20S59U5FMtq8Xnl1vN0/aefjj8+P1bWbDNGaYVVK18jgnpedpmp13J5zzJhosdTXILzQm3xVsawVHFI88zXY9MkuIttE5Mci8pMiP/sniv2XiMghEXlSRL4mIls2uxYhZH6kyPhTAK5T1SsA7Adwo4hcDeCTAD6tqm8A8AKAm9trJiGkKSkZYRTAy8XL1eJPAVwH4E+L/XcB+HsAn519E2swrQe+MmGDN0+8fBuIJqyYY8OKiSzWgz6wnnUr3bfDbIeSbn2Hke47zOSRc8du8u3njF3j27ca9zmALSvlk71Pro2/Iq+eGjfy5IlQyA1eGZ83NKMRKzaox87biE0y3Xxb7Dz1gf+ZqTVR7IQdI9Ul0tR2kpCtpzLJRpCAImF/TI157279s0oSISL9Is/bcQD3AfglgBdV9cw35CiAi9JaSAiZB0mdXVUHqrofwF4AVwG4LLWCMD/7yZrNJIQ0ZSpvvKq+KCL3A7gGwE4RWSme7nsBPOOUMfnZL5zdejtVgTQppC6LVSNXursibBDU4s9Ht3HvnryPZfzg3LFeXD1/LNcv3PnyxvbF5724sf367f8XlH/t6sso4/jp8za2f31i18b2b166IDjvhZUd47YFvlrjzQ5keFiPja0frlvTx+w390iiZbmsrFc7gtB3JHU8AuOtSBs0Mlod1ua0C+aW28QUtkyVpm+fFG/860RkZ7F9DoAbABwBcD+A9xWnMT87IQtOypN9D4C7RKSP0Y/D11X1XhF5DMBXReQfADwE4PMttpMQ0pAUb/xPAVxZsv9XGNnvhJAlYDEi6OZNbKOlTNJJzHtu13rTwM6PzrOvrflnbdZVEw1nJrgAgJjouHPMENueHS9tbL/p/LFb5dodvwjKX7W13Hl6+PTY/v7R6tgvO4wWsVtbH9vmvzORdsPTxv4+Vb6eHhDPu7c50WG20/wkKX6WCbzzguEtv3gXxBGk00bUMTaekExgZyckEyjjgeZzidvEiyxLZOjMMOlFEWRbZfxV6Jvhom0yjrRbNeNlKxNLOdtQN6cxgUnjtTg65m3XIVXSn6XwyU5IJrCzE5IJ+cp4O3Eg/snz5j1721W5ym1O8mA7jgCz215GFTNPPFrWaf3k2Bt+4sQ41O7ZLeMIuId7eze2fxdNiD+85XfjeozE/+3auPxTJy7c2P6fl18TlD9xcuy111PjtvRNm8W+lzi5irOKbbht56ZH5avmqpftX2TTrSX4ZCckE9jZCcmE5ZXxQSK7eI0ju5SPnafsTHaIJZ0n0W2VNtd4tMRTKMNNgMh6+X4gTKBgA1HswrXDYFmnOKjELGuFsaR+fjCW4b9/dSzvf7MtnMiyzcxn7xmNfGowbsArZj67NRUAYPDyOMqn94qR8cESVePze5EZYpNG2EkxQd53a+pEplPw2h7y5olXfeZVJtocaX1ZKkLI2QE7OyGZsJgyPpBRNgeYWW7KRljEUqvvHPOWJ1VfEqpdyshKQutNX48863aJo4H1ppu8Z/FqrtaqcDMeWJMkLN5bNzLemAGDE+P9r241Mnw1ukC82uwZ7GmmDjkVPidWTlq5Xi7dV16158e53sz2aZsfrjyne5zfPTg2MCu6BmZA1QhKyqhLtJRVcJ4TOJ+QN70r+GQnJBPY2QnJBHZ2QjJhbja7tXdkGBugCb9BQQRcnGrDTnhw6hzYcZx4Prtjm9u1zSomVfTc+dTlTSxKxTsmmhIMPcX50e1wncnOYtez0yBve5TFxvETePXHiR1tewL7O9jW0v1AaKevmFztPVvmtNkfJXwM7XT7+SVsx2VS7HfAH5ZLHa5rmtiR+dkJIWWwsxOSCYs59Gbxhj765Xm3R+cFM1HMfpsRxh9SsVi57oqmSMbbV8GvaSAJo0u4mU9MMkUjlQdhQhcMjCweBnLdS6wYReB588adW9mLzAgJogPNeTZvuzOkNjrPSnxzLUe699bCmTRiXsv6sHQbntQHws/GMQMmItjUk/62vBnGTZXqFd/HJvDJTkgmsLMTkgndy/gzssauFDoRwWb1YvnvUejNjw5a73yQBcSWt1K5IlOHldS2fiPdZS3StHHy8TPNqpLxNtLLRqrZiTRGkvdPhzI8SBoZJJY05W0yynh125QFdZ155gDQcyas9JzJP71Yxlu5HpxXLr0lUcZj3ZxX4Y1Hije+KurO8YxrRQReige+6eQXS/KTvUju+JCI3Fu8Zn52QpaIaWT8rRilfToD87MTskQkyXgR2QvgTwD8I4C/lpGLull+9iAopkJDepLemSADRBLTSnqrVY1018iz78p6OynG7LZJDWKC5ZLMaICVvQAAE+QSeOBNwgUv7/vomDkvmAPvJFaI06OnyHhnua1Rm+12uYll5XU8nz/woAfrATjSfT2UxIHpt24lvSPjY298cMx8NwJ5H5fZ3AOfTJXXP+W6G/X756Q+2T8D4GMYz4G6EMzPTshSkZLF9d0Ajqvqg3UqYH52QhaDFBn/VgDvEZF3AdgG4HwAd6DF/OyBp91bzL9qnrGd9x4cMpLOSvdItltZ70p6u3RSRcy0Gukuaq8blrHn2bnpCOS5mRvfj2Pbh6XnedJ9Qran5KuvWN3VXVE38NIPS88B4uCX8jkMwT1bjz6XYbnXPvhsPY97dJ4r3ScCcRKCZ6ri31OCZ5KkexqbPtlV9XZV3auq+wB8AMAPVPVDYH52QpaKJkE1H8fIWfckRjY887MTssBMFVSjqj8E8MNim/nZCVkiOo6g07GdIXYYLbJLnOg6N7KuwhZKst/FFzhqjgVWbWCLRRFzfafNA2tXR8OFNmrOHNM1u9/a4lGbPTvd3ktnbv0EVctsb+yPijjRZIFtPnDuS9UxZ1LKxEQWr0wwwcXY5XH93jBa4tBbsp3uMMtIOQ/GxhOSCezshGTC/OazW0kUy2gnui5J0o9ONNdKkPTxslZ2VKdnpVr5fPqJoTcr660k9uQ9ELwHtctf2fdm5P3EkKSdd2+j6xJzkqtzbOK9bRSoGhJylm+uWuLJk+vDRHmcIt2d/RNlEiPjppbuVUtRB+d597zZPHc+2QnJBHZ2QjJhMZaliuWJ56lPkfQxKZNnomg2seFhwUq1ztz6eL8XATiwkj7y4Bu7QnqO9LZSP454s/PrU6R7VD4hfi6kMrJrc+keR9C5EZFD51oVMlxTJHlFBF3S3PS4nVXz1r3ywbWaSXc9U76iCj7ZCckEdnZCMqFbGa9juSFVEy88T32CpAciGet57SuSRCR57QMVGXvWvUCWcpNidMgcs9tB28ykkFiee9I/qKTinlcdO0Nq4EedXOcpyRQ9qR5f25Pr3jnRscqlpLw2z1O6J8InOyGZwM5OSCbML9ebjZ9uKunjIs5PmOu1j4NyrNfeuXDgsY/bbyW+bbPNyR5NKA8lmV2edUovO+DnykuR6nVJDCTxzk8KMKnypqdI4tSEDamx7U2XkrIkSPdpZXsMn+yEZAI7OyGZwM5OSCYsRARdlS0S2POeXVNjIk1QPDU/fIItD8TDYsb+rlqvObDt7dChLZ5Y3rH5UwmSWTadZ11neClluK5qWefgUjWG0YLr+mWWwU638MlOSCawsxOSCd3L+LJlqSpPTxiiS51IY0lNLOnhReOhQm4HQ3fxOYlyf+MC8f2zkYLlZkAqM10gqamkTZ2UEpRJG+5LKVMrAi64QDvRcHXq4ZOdkExgZyckExZzWSq3iONNj+V9ite+RgRecKnUOfQWu0RSlWc9oQEykZLF1mPb0mLUnEcNSep7thPtkIblu/CsTxaZ8j41XJYqNYvrUwB+j9HXaF1VD4jILgBfA7APwFMA3q+qLzRqDSGkNaaR8e9Q1f2qeqB4fRuAg6p6KYCDxWtCyILSRMa/F8Dbi+27MMoU8/FaV6qSJwkSv9WgHK/OqmYNKuade1R4983FxvWnXTXNJmmTVBluSQnkaSrPJ643O896WKR9L3sqqd8EBfB9EXlQRG4p9u1W1WPF9rMAds+sVYSQmZP6ZL9WVZ8RkT8AcJ+I/NweVFUVx2NU/DjcAgDbsL1RYwkh9Unq7Kr6TPH/uIh8C6OEjs+JyB5VPSYiewAcd8qO87PLrhpu2kTp7RZP8OCnSqVEb35Qf6J2SgrkqSJhBGChqPF+a8Xp15HRSy7XPTb9KorIDhE578w2gHcCeATAPRjlZQeYn52QhSflyb4bwLcKR9MKgH9R1e+KyAMAvi4iNwP4NYD3t9dMQkhTNu3sRR72K0r2/y+A69toFCFk9izEfPZa1IjAC4snDtd5dVbRom3vUme4r2Namxuf3IBmdvEy2OVVMDaekExgZyckE+Ywn72QQrOUmnXkdeXlNpdryctfpzKc3e+u2vi6eUyEscxwWaWAGUviWS7/NBe5nmAi8clOSCawsxOSCXOcz14hO9ryJjeMxgsvlSb7KuV+eMGp21BR6Xi7LRldhznI25nK88mLt3ft0vqYEYYQkgA7OyGZsJhBNZ5c6VreT9Tfntz3q6zxnuccvDFLWpXhfqXd1xnU38575pOdkExgZyckE9jZCcmExbTZPbq25Sfqr2HL1bDzwyoXaOhsGZm3/V1FS7a5B5/shGQCOzshmbBcMt4jVQ7NY553UxnZ0AxYehZZhnt0LM9TyfybREg+sLMTkglnh4xPZZHlvscyytizlQWV56nwyU5IJrCzE5IJecn4VKaVa4sk+0k9llyip5D0ZBeRnSJyt4j8XESOiMg1IrJLRO4TkSeK/xe03VhCSH1SZfwdAL6rqpdhlDDiCJifnZClYlMZLyKvAfA2AH8OAKp6GsBpEZldfvZlZxYSkKZAMzKQ4U1JebJfAuB5AF8UkYdE5HNFgkfmZydkiUjp7CsA3gzgs6p6JYBXEEl2HeX1cfOzi8hhETm8hlNN20sIqUlKZz8K4KiqHipe341R53+uyMuOzfKzq+oBVT2wiq2zaDMhpAabdnZVfRbA0yLyxmLX9QAeA/OzzxbVbv+Wvf1tv5+zkNRx9r8C8GUR2QLgVwD+AqMfCuZnJ2RJSOrsqvowgAMlh5ifnZAlgRF0uULpmx2MjSckE9jZCckEdnZCMoGdnZBMYGcnJBPY2QnJBHZ2QjKBnZ2QTGBnJyQT2NkJyQR2dkIygZ2dkExgZyckE9jZCckEdnZCMoGdnZBMYGcnJBPY2QnJBNEOlycSkecxWnf+t51VOslrWT/rP4vrf72qvq7sQKedHQBE5LCqli1eyfpZP+tvEcp4QjKBnZ2QTJhHZ79zDnWyftafff2d2+yEkPlAGU9IJnTa2UXkRhF5XESeFJHbNi/RuL4viMhxEXnE7NslIveJyBPF/wtarP9iEblfRB4TkUdF5NYu2yAi20TkxyLyk6L+TxT7LxGRQ8Xn8LUih19riEhfRB4SkXu7rl9EnhKRn4nIwyJyuNjX5Xdgp4jcLSI/F5EjInJNl/VbOuvsItIH8M8A/hjA5QA+KCKXt1ztlwDcGO27DcBBVb0UwEFEueZnzDqAj6rq5QCuBvDh4j131YZTAK5T1SsA7Adwo4hcDeCTAD6tqm8A8AKAm1uq/wy3AjhiXndd/ztUdb8Z8uryO3AHgO+q6mUArsDoPnRZ/xhV7eQPwDUAvmde3w7g9g7q3QfgEfP6cQB7iu09AB7v8B58G8AN82gDgO0A/gvAWzAK6lgp+1xaqHcvRl/o6wDcC0A6rv8pAK+N9nVy/wG8BsB/o/CNzfs72KWMvwjA0+b10WJf1+xW1WPF9rMAdndRqYjsA3AlgENdtqGQ0A8DOA7gPgC/BPCiqq4Xp7T9OXwGwMcADIvXF3ZcvwL4vog8KCK3FPu6uv+XAHgewBcLM+ZzIrKjw/oDsnbQ6eintfXhCBE5F8A3AHxEVV/qsg2qOlDV/Rg9Ya8CcFlbdcWIyLsBHFfVB7uqs4RrVfXNGJmPHxaRt9mDLd//FQBvBvBZVb0So1DxQLJ39R0Euu3szwC42LzeW+zrmudEZA8AFP+Pt1mZiKxi1NG/rKrfnEcbAEBVXwRwP0ayeaeInEnX3ebn8FYA7xGRpwB8FSMpf0eH9UNVnyn+HwfwLYx+8Lq6/0cBHFXVQ8XruzHq/J1//kC3nf0BAJcWntgtAD4A4J4O6z/DPQBuKrZvwsiObgUREQCfB3BEVT/VdRtE5HUisrPYPgcjf8ERjDr9+9quX1VvV9W9qroPo8/7B6r6oa7qF5EdInLemW0A7wTwCDq6/6r6LICnReSNxa7rATzWVf1lDersD8C7APwCI7vxbzuo7ysAjgFYw+hX9maMbMaDAJ4A8O8AdrVY/7UYSbSfAni4+HtXV20A8CYADxX1PwLg74r9fwjgxwCeBPCvALZ28Fm8HcC9XdZf1POT4u/RM9+5jr8D+wEcLj6DfwNwQZf12z9G0BGSCVk76AjJCXZ2QjKBnZ2QTGBnJyQT2NkJyQR2dkIygZ2dkExgZyckE/4fATo+hxEOxMIAAAAASUVORK5CYII=\n"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "execution_count": 32,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "i = best_models.index[0]\n",
        "m = best_models.loc[i]['Model']\n",
        "p0 = np.concatenate((comp_to_p(m['disk']), comp_to_p(m['bulge'])))\n",
        "_m = {'disk': dict(m['disk']), 'bulge': dict(m['bulge']), 'bar': None, 'spiral': np.array([])}\n",
        "_m['bulge']['n'] = 4\n",
        "image_size = gu.get_diff_data(i)['width']\n",
        "psf = gu.get_psf(i)\n",
        "cx, cy = np.mgrid[0:image_size, 0:image_size]\n",
        "\n",
        "im = bulge_disk_model(p0, cy, cx)\n",
        "\n",
        "im2 = rendering.calculate_model(\n",
        "   _m, psf=None, image_size=image_size\n",
        ")\n",
        "f, ax = plt.subplots(ncols=3, figsize=(20, 5))\n",
        "ax[0].imshow(im)\n",
        "ax[1].imshow(im2)\n",
        "c = ax[2].imshow(im - im2)\n",
        "plt.colorbar(c, ax=ax)"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 44,
          "data": {
            "text/plain": [
              "<matplotlib.colorbar.Colorbar at 0x12f4565f8>"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1440x360 with 4 Axes>"
            ],
            "image/png": [
              "iVBORw0KGgoAAAANSUhEUgAAA/kAAAE1CAYAAAC88CDqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dfbBkd3kn9u8zIyEQbwILywKxEWsoCOUY4ZpgY1ybNQJb8bpW7MZ22eUiSqKUdmtNCldc5YCp2uxbqrxxYq+34nUyMRj9wRowhqBiiUHIbHldsWUG82KBzCKzEKQIZC2wwrwIaebJH7cF97Zmbr93n+7+fKpOzT19us957u07z8zTv9/zO9XdAQAAALbfiU0HAAAAACyHIh8AAAB2hCIfAAAAdoQiHwAAAHaEIh8AAAB2hCIfAAAAdoQiHwAAgJ1VVW+oqvuq6o4lne93q+pLVfWuCxz/51X1l8u41jwU+QAAAOyyNya5bonn+6Ukrzzfgao6leQpS7zWzBT5AAAA7Kzu/v0kXzj8WFV952hE/oNV9W+q6nkznO+2JF8ef7yqTubgA4CfXzTmRVy0yYsDAADABpxO8ne7+5NV9b1J/kWSly54zlcluaW7762qhQOclyIfAACAvVFVT0jy/Ul++1Axfsno2N9O8o/O87J7uvuHjznn05P8eJK/vtRg56DIBwAAYJ+cSPKl7r5m/EB3vz3J2+c45wuTPDvJXaMPDi6tqru6+9kLRToHPfkAAADsje5+IMm/q6ofT5I68IIFz/mvuvs7uvvq7r46yVc3UeAninwAAAB2WFX9VpI/TPLcqrq7qm5M8tNJbqyqjyT5WJLrZzjfv0ny20muHZ3vgtP4N6G6e9MxAAAAAEtgJB8AAAB2hCIfAAAAdoTV9QEAAGCCH/7Bx/e//8LZmV/3wY8++J7uvm4FIZ2XIh8AAAAmuP8LZ3P7e66a+XUXX/nnl68gnAtS5AMAAMBEnbN9btNBTLRQT35VXVdVn6iqu6rqNcsKCmBI5DpgH8h1AMfrJOfSM2/rNvdIflWdTPJrSV6e5O4kH6iqW7r74xd6zWPqkn5sHj/vJYEd9PV8Jd/oB2vTcVyIXAcsy5fzxfu7+2mbjuN85sl1Fz3p0r742y9bV4jAlvj6n9872Fy3DOcy/JH8RabrvyjJXd39qSSpqjcnuT7JBf8xeGwen++taxe4JLBrbu/bNh3CJHIdsBTv67d9ZtMxHGPmXHfxt1+W7/zl/3ZN4QHb4mPX/+Mh57qFdDpne/0j87NaZLr+M5J89tD+3aPHjqiqm6rqTFWdeSgPLnA5gI2Q64B9MHOuO/vAV9cWHMBQbMN0/YV68qfR3ae7+1R3n7o4l6z6cgAbIdcB++Bwrjv5pEs3HQ4A57HIdP17kjzz0P5Vo8cAdolcB+wDuQ5ggk5ydgMj87NaZCT/A0meU1XPqqrHJPnJJLcsJyyAwZDrgH0g1wFMYRum6889kt/dD1fVq5K8J8nJJG/o7o8tLTKAAZDrgH0g1wFM1slWLLy3yHT9dPe7k7x7SbEADJJcB+wDuQ5gsuHfQG/BIh8AAAD2Qae3oidfkQ8AAACTdHJ2+DW+Ih8AAAAm6ZiuDwAAADuicja16SAmUuQDAADABJ3knOn6AAAAsBuM5AMAAMAO6CjyAQAAYGeca0U+AAAAbD0j+QAAALAjOpWzObHpMCZS5AMAAMAUTNcHAACAHbAt0/WHP9cAAAAAmIqRfAAAAJiocraHP06uyAcAAIAJOsm5LZgMr8gHAACAKWxDT74iHwAAACbo3o7p+sOPEAAAAAbgXGrmbZKqemZVvb+qPl5VH6uqVy8So5F8AAAAmODgFnorGSd/OMnPdfefVNUTk3ywqm7t7o/PczJFPgAAAEy0mun63X1vkntHX3+5qu5M8owkinwAAABYhQVW17+8qs4c2j/d3afP98SqujrJC5PcPs+FEkU+AAAATOVsz7W6/v3dfWrSk6rqCUl+J8nPdvcD81woUeQDAADARJ1aVU9+quriHBT4b+ruty9yLkU+AAAATOHcCnryq6qSvD7Jnd39y4ueT5EPAAAAE6xwdf2XJHllkj+tqg+PHvuF7n73PCdT5AMAAMAEnZq3J//483b/QZKlnViRDwAAAFOYc3X9tVLkAwAAwATdydkV9OQvmyIfAAAAJqqcW96s+pWZ+DFEVb2hqu6rqjsOPfbUqrq1qj45+vMpqw0TYLXkOmBfyHcAu22auQZvTHLd2GOvSXJbdz8nyW2jfYBt9sbIdcB+eGPkO4CZdQ6m68+6rdvEK3b37yf5wtjD1ye5efT1zUleseS4ANZKrgP2hXwHML+zOTHztm7z9uRf0d33jr7+XJIrlhQPwJDIdcC+kO8AJuhUzq3gFnrLtvDCe93dVdUXOl5VNyW5KUkem0sXvRzARsh1wL44Lt8dznUXP+3Ja40LYAg2MTI/q3kj/HxVXZkkoz/vu9ATu/t0d5/q7lMX55I5LwewEXIdsC+myneHc93JJ/lAE9gvneRcn5h5W7d5r3hLkhtGX9+Q5J3LCQdgUOQ6YF/IdwATVc7Osa3bNLfQ+60kf5jkuVV1d1XdmOQXk7y8qj6Z5GWjfYCtJdcB+0K+A5jPtozkT+zJ7+6fusCha5ccC8DGyHXAvpDvAOa3iZH5WS288B4AAADsuu7ayMj8rBT5AAAAMIWzinwAAADYfp3knOn6AAAAsAvKSD4AAADsgoPV9Y3kAwAAwE44O/ku9BunyAcAAIAJOmUkHwAAAHbFOSP5AAAAsP26k7NbMJI//I8hAAAAgKkYyQcAAIApbENPvpF8AAAAmOBg4b0TM2+TVNUbquq+qrpjGXEq8gEAAGAKZ1Mzb1N4Y5LrlhWj6foAAAAwQWc10/W7+/er6uplnU+RDwAAABPVVNPvN02RDwAAAFM4N930+3GXV9WZQ/unu/v0kkJ6FEU+AAAATNCdnJ1vuv793X1q2fFciCIfAAAAprAN0/WHHyEAAABs2MEt9GbfJqmq30ryh0meW1V3V9WNi8RpJB8AAACmMGdP/rG6+6eWeT5FPgAAAEywqlvoLZsiHwAAAKawDT35inwAAACYZMoe+01T5AMAAMAEndX05C+bIh8AAACmYCQfAAAAdsC2LLw3/FUDAAAAgKkYyQcAAIApbMNIviIfAAAAJuhYXR8AAAB2xjasrj+xJ7+qnllV76+qj1fVx6rq1aPHn1pVt1bVJ0d/PmX14QKshlwH7AO5DmABfTBdf9Zt3aZZeO/hJD/X3c9P8n1Jfqaqnp/kNUlu6+7nJLlttA+wreQ6YB/IdQBzemR1/a0v8rv73u7+k9HXX05yZ5JnJLk+yc2jp92c5BWrChJg1eQ6YB/IdQCL2YYif6ae/Kq6OskLk9ye5Iruvnd06HNJrrjAa25KclOSPDaXzhsnwNrIdcA+WDTXXfy0J68+SIAB2ZaF96aZrp8kqaonJPmdJD/b3Q8cPtbdnYPZC4/S3ae7+1R3n7o4lywULMCqyXXAPlhGrjv5JB9oAvunu2be1m2qkfyqujgH/xC8qbvfPnr481V1ZXffW1VXJrlvVUGyQjX8T6Lm0uf9vwkcS67bYXIdfJNcBzC/XVldv5K8Psmd3f3Lhw7dkuSG0dc3JHnn8sMDWA+5DtgHch3A/HpLVtefZiT/JUlemeRPq+rDo8d+IckvJnlrVd2Y5DNJfmI1IQKshVwH7AO5DmABm5h+P6uJRX53/0FywTkJ1y43HIDNkOuAfSDXASxiOxbem2l1fbbArvadzmqRn4MeVxg+ue6AXAcAa7UTI/kAAACw7zoxkg8AAAA7obdjIpwiHwAAAKawDbfQU+QPnb7T9Zv0M9+Gj+9g28h16yfXAcBOOrHpAAAAAGDoOgcL7826TaOqrquqT1TVXVX1mkXiNJIPAAAAE63mFnpVdTLJryV5eZK7k3ygqm7p7o/Pcz4j+QAAADCF7tm3KbwoyV3d/anu/kaSNye5ft4YjeQPwbb2otYaPyPqc+u71iT6WGE+ct1kch0ADNq00+/HXF5VZw7tn+7u04f2n5Hks4f2707yvfNcKFHkAwAAwEQHI/NzFfn3d/epZcdzIYp8AAAAmMIqevKT3JPkmYf2rxo9Nhc9+QAAADCFFfXkfyDJc6rqWVX1mCQ/meSWeWM0kr8JQ+pLXWev6SIWjXOdfa7Hvb96WNknct3s5DoAGLQ5p+tPOGc/XFWvSvKeJCeTvKG7Pzbv+RT5AAAAMEFn+vvez3zu7ncnefcyzqXIBwAAgClsw1w1RT4AAABMMv/q+mulyF+HdfalrrHvtE6s7/vqcwt+ZjbLz2WVPa3jvwv6Vtklct3C5DoAGLgt+CdNkQ8AAABTMJIPAAAAO2IbJqcp8ldhi6asrnQa6hKnjS4a50xTYCfFvcwprqa0ss3kukdOPv1z5brRueU6ALZPx0g+AAAA7IZOsgVF/vpWLgIAAABWykg+AAAATGEbOs4U+cuwyr7UJd8maqZ+z0WvPVNv6cmju4veRmqsl/S473vmW1aN/1z0rbIv5Lrzk+tmJ9fB4CzSZ1x19O/w+LlmPQ6DtgW/rop8AAAAmKgsvAcAAAA7w0g+AAAA7IB2C73dtc57Q89o5nssH9eLOuO5apk/l5PHn6sn9W+eO+b7mqGH9eDpE66lb5VdJded/1Ry3XnPvxC5DpZu1YXI4fNP6rkfN+n5J8aPj7/+mGOrNstP9bi4z3ecLbEFb5wiHwAAAKYy3EGQRyjyAQAAYBpbMJI/8b5BVfXYqvrjqvpIVX2sqv7h6PFnVdXtVXVXVb2lqh6z+nABVkOuA/aBXAewoJ5jW7NpRvIfTPLS7v7Lqro4yR9U1f+d5L9P8ivd/eaq+t+T3Jjk11cY636Y8X7Nq+xLndh3emLsXLM+fxbnxnpLx4+P92/WubHDh44f18OaTOxj3WjfKqsk162TXHd+ch2rJ9ftsEX65I871/nOd+7c4ePH55ujz01OnjzuXEmdnBD7oVgm9bmvcmL1rLXbFgwAM0kn2YKF9yb+T6QP/OVo9+LR1klemuRto8dvTvKKlUQIsAZyHbAP5DqAxXTPvq3bVMMNVXWyqj6c5L4ktyb58yRf6u6HR0+5O8kzLvDam6rqTFWdeSgPLiNmgJWQ64B9sKxcd/aBr64nYIAh2YLp+lMV+d19truvSXJVkhcled60F+ju0919qrtPXZxL5gwTYPXkOmAfLCvXnXzSpSuLEWCwumbf1mym1fW7+0tV9f4kL05yWVVdNPrU96ok96wiwL0wY2/qUs994nDv1tgv4KQ+1LHjdXL8+WP7h4+PHxvvtx3vBR3v9Tw7vn/26NPH9uvQfo//SMavNeQ+U/eSXgu5bkXkOrluWnLdWsh1+2XWnvyzZ4/miJMnj+aIb3zjogsee9S5Hj55ZP+JT/jakf2Hx2I7Od7/f3yoR0z6Lidlk1l+SouWbjLbdqoteOOmWV3/aVV12ejrxyV5eZI7k7w/yY+NnnZDkneuKkiAVZPrgH0g1wEsYJ6p+gNdXf/KJDdX1ckcfCjw1u5+V1V9PMmbq+qfJPlQktevME6AVZPrgH0g1wHMbTPT72c1scjv7o8meeF5Hv9UDvq4ALaeXAfsA7kOYLtU1Y8n+QdJ/uMkL+ruM5NeM1NPPps3872ixx13f+hZ+1IvGvv1GetTrUvGFh+75DHf/PLcpY89cujcpRcffe1YH+qJB472btVXv35kv79+dP9R90w91M9ZY/ehnti3Ombme0kDM5PrRsflOmBJxu9zP96jP770xfi963ssiZz98rfy2bkHx9cmObp78ReOHj9xzVeOjfXsWGzLHDfd5BjseBZ9VA5fVyAsZv1v1B1J/naS/2PaFyjyAQAAYBprLvK7+87kPAsHH0ORDwAAANOYr8i/vKoOT7M/3d2nlxPQoynyAQAAYJLOvAvv3d/dpy50sKrel+Q7znPodd09891OFPlM7dH3lh7rlxrvWz3Ul5okD1/+xG9+fd9/+sQjxz70un9xZP+7/vnfO7L/jN872st10UMPH73WN8aavzJ2fIbpLcB+k+uAbXW47368537ya2e82MlvXeuS+4/mpr/yj/6fI/tP/6OjufDzXzu6/7XHHl2v5OsPH82zZ8+NrQdwaL2AR68tcHT/3ApXQp/UYz9uUg/+ca/Xrz8ctYI3o7tftszzKfIBAABgGlvwicv4WrsAAADAAFTV36qqu5O8OMm/qqr3THqNkfwtM34Lo4VvM7WIc+O3YhmLbWya6ckvf+vWT0/70NHPl6595Y1H9q/6Dw8cfe0X/vLotb/+4NH9h49eazyWw/eFedSxGW8L5TZSsHpy3YhcB8xp/JZ5sx6/5JKjt+E8Ofb8x1/xrfz0917y+0eOveH7v//I/m/+lbcf2f87d7/4yP4nvnTFkf2LTxy99oMPXbhkmbUtYZkWvbJb6G2nVUzXP053vyPJO2Z5jSIfAAAAprHBD5ampcgHAACASTpbMeVCkQ8AAADTUOQzlT7Ud1QrXgtxvLf00OXq3NH+p5w4GkufPXtk/1ETVR48+voe6x093Ft68ReP9qFeNKF39FHnGotlPLY86nxj39txeobnrtv49wXbRK6T66Yl18F5LfMWbOM99uO3mjs51hd/8sTR5z989lu5863/39Hbf3/tG0dviffDd/7okf2vPHT01qNfG+u5H7/WeN/9pPUEDjsx4ftcpWXfYo9hWHdP/jwU+QAAADANRT4AAADsCEU+AAAAbL9q0/VZg4n3kh7vuRzvgz30+h479Ki+1fFrj/dM1ti1z46//qFvvfbYM2dyX+mEPtSZ7g89oS915ntFD7nPFbaUXHf+58t1sN9m6S6ftRN9vA/+3DFrZXzx64879lxf+NqlR2MZ75M/N75Oy6S8fOHvZp0995PM/DM/5vVbUFfujwH9jl2IIh8AAACmsQWfuKx4eWMAAABgXYzkAwAAwBT05O+qSf2ZC517Ql/pxJcv0Lc6fr/m8bko4+cav5f0o/pUx+7nvIBH9Z2Om9RLOkPv6KD6Ut0rmk2S6w7IdTOfe2ZyHSRZ733SJ/Wujx+vR/Xof2v/4bNHc/h4z/zDZ4+/z/2jYnlUj/7xse0qmXGgtuCNUeQDAADAJFbXBwAAgB2iyAcAAIAdocjfE7vUt3qcCf1Rj+prXadV9qEucK3Zz70FWYP9JdcdnEquW5xcB1MZUuf5cfemHz823jN/Ymx+83HnOp+h9uAfd1/7aZ4/yeHzyZrDsQ3T9d1CDwAAAHaEkXwAAACYxhaM5CvyAQAAYBKr6++xdfatjpvQxzpLj+ZCPa1rtnDv6dGTLe9cE6+1BVkCLkSuWzu5DnbLeNactcd7nY7ri5/UMz9rT/1Qe/DHzRrlpPd7nMw5UFvwxijyAQAAYBpbUORPvfBeVZ2sqg9V1btG+8+qqtur6q6qektVPWZ1YQKsh1wH7AO5DmB2lYPp+rNu6zbL6vqvTnLnof1/muRXuvvZSb6Y5MZlBgawIXIdsA/kOoB59Bzbmk1V5FfVVUn+RpLfGO1XkpcmedvoKTcnecUqAtwJ3Ue3lV7r3PHbTKfqrdk29TOa2Tp/F5iZXLcguU6u++a15bohk+tWqw5ty7bhuoEJFn1/Kqv9/WEJ5hjFX3Qkv6p+qar+rKo+WlXvqKrLJr1m2pH8f5bk55M88r+Cb0vype5+eLR/d5JnzBwxwLDIdcA+kOsA5rX+kfxbk3xXd393kn+b5LWTXjCxyK+qH01yX3d/cJ6IquqmqjpTVWceyoPznAJg5eQ6YB8sM9edfeCrS44OYAusucjv7vce+hD2j5JcNek106yu/5Ikf7OqfiTJY5M8KcmvJrmsqi4aXfCqJPdcIKjTSU4nyZPqqWYWAUMl1wH7YGm57nHPfrpcB7Be/02St0x60sSR/O5+bXdf1d1XJ/nJJL/X3T+d5P1Jfmz0tBuSvHP+WPfMJnsVJ/Vorqtfc9mG/H3pS90Kct0KyHXLN+TvS67bCnLdZo33XI8P9unJ3m2TBnetubAd5uzJv/yRWVCj7aYj56x6X1XdcZ7t+kPPeV2Sh5O8aVKM04zkX8j/kOTNVfVPknwoyesXOBfAUMl1wD6Q6wCmMd8nMPd396kLnrL7Zce9uKr+qyQ/muTa7smfps9U5Hf3v07yr0dffyrJi2Z5PcA2kOuAfSDXAcxoA9Msquq6HCyW+p9191SLoSwyks+yHPdhTG1wsta2TWPdFFNTYTpy3XaT62Aqx/1NmfS3aPy4KfvDM0smlDV306K3xJvD/5bkkiS3HtzxNH/U3X/3uBco8gEAAGAa615mqPvZs75GkQ8AAABT2MBI/swU+QAAADANRT4Lm9QDuck+1n2hDxVWT67bPLkO1m48sy3zb6H+/gPL/jkc93pZdA9syf0NFfkAAAAwQWU7PjBT5AMAAMA0jOQDAADAbrDwHqs3aw/lvva16jWF7SbXTUeug0GbtQd/1r/Rh5+/aBZcZi/7svviZ/k+Jx2fFNtx78Ge/kvDFvxTq8gHAACAaSjyAQAAYAe06foAAACwOxT5DI5+TWAfyHXAAC07M81yvll7z2fpZV91X/ws51ukx35W/qXZT9swkn9i0wEAAAAAy2EkHwAAAKaxBSP5inwAAACYwjZM11fkAwDAjll1//8y++IX6dmf9H3O+nPYgvqNTepsxS+JIh8AAACmocgHAACA7VcxXR8AAAB2hyIfAAAYmkV72Wepc/TNs0uqh/8bqMgHAACASSy8BwAAALtDTz4AAADsCkU+AACwy8bvc79oDbTs88EyGckHAACAXbEFRf6JTQcAAAAAg9cHI/mzbouoqn9cVR+tqg9X1Xur6umTXqPIBwAAgGn0HNtifqm7v7u7r0nyriR/f9ILTNcHAADmtuzZy1swG5o9VVl/T353P3Bo9/GZ4q+IIh8AAACm0XNV+ZdX1ZlD+6e7+/S0L66q/ynJf5nkPyT5wUnPn6rIr6pPJ/lykrNJHu7uU1X11CRvSXJ1kk8n+Ynu/uK0gQIMjVwH7AO5DmDt7u/uUxc6WFXvS/Id5zn0uu5+Z3e/Lsnrquq1SV6V5H887mKz9OT/YHdfcyi41yS5rbufk+S20T7AtpPrgH0g1wHMYRUL73X3y7r7u86zvXPsqW9K8l9MOt8iC+9dn+Tm0dc3J3nFAucCGCq5DtgHch3AJPMsurf46vrPObR7fZI/m/SaaYv8TvLeqvpgVd00euyK7r539PXnklxxgaBuqqozVXXmoTw45eUANkKuA/bBUnLd2Qe+uo5YAQalzs2+LegXq+qOqvpokh9K8upJL5h24b0f6O57qurbk9xaVUc+Pejurjr/RITRggKnk+RJ9VSLZQJDJtcB+2Apue5xz366XAfsn/Wvrj9xev64qUbyu/ue0Z/3JXlHkhcl+XxVXZkkoz/vm/XiAEMi1wH7QK4DmN8qevKXbWKRX1WPr6onPvJ1DqYI3JHkliQ3jJ52Q5LxRQEAtoZcB+wDuQ5gAZ2DW+jNuq3ZNNP1r0jyjqp65Pn/srt/t6o+kOStVXVjks8k+YnVhQmwcnIdsA/kOoAFbGJkflYTi/zu/lSSF5zn8X+f5NpVBAWwbnIdsA/kOoAF7UKRDwAAAPuusiMj+QAAALD3NtRjPytFPgAAAEzBSD4AAADsCkU+AAAA7AYj+QAAALALOsm54Vf5inwAAACYxvBrfEU+AAAATGMbpuuf2HQAAAAAwHIYyQcAAIBp9PCH8hX5AAAAMIVtmK6vyAcAAIBJOhbeAwAAgF1QScp0fQAAANgR5zYdwGSKfAAAAJiCkXwAAADYBXryAQAAYFe0W+gBAADArnALPQAAANgVRvIBAABgB3RSW7C6/olNBwAAAABboXv2bQmq6ueqqqvq8knPNZIPAAAA09jAbP2qemaSH0ry/07zfCP5AAAAMIXqnnlbgl9J8vOZ8iMGI/kAAAAwjfmK9sur6syh/dPdfXqaF1bV9Unu6e6PVNVUF1PkAwAAwCSdZL6F9+7v7lMXOlhV70vyHec59Lokv5CDqfpTU+QDAADABJWlTb8/ortfdt7rVf0nSZ6V5JFR/KuS/ElVvai7P3eh8ynyAQAAYGC6+0+TfPsj+1X16SSnuvv+416nyAcAAIBprGAkf9kU+QAAADCNDRb53X31NM+b6hZ6VXVZVb2tqv6squ6sqhdX1VOr6taq+uToz6csFDHAhsl1wD6Q6wDm9MjCe7NuazZVkZ/kV5P8bnc/L8kLktyZ5DVJbuvu5yS5bbQPsM3kOmAfyHUAc6rumbd1m1jkV9WTk/y1JK9Pku7+Rnd/Kcn1SW4ePe3mJK9YVZAAqybXAftArgNYUPfs25pNM5L/rCR/keQ3q+pDVfUbVfX4JFd0972j53wuyRWrChJgDeQ6YB/IdQBzm6PAH2iRf1GS70ny6939wiRfydgUru7uHHQoPEpV3VRVZ6rqzEN5cNF4AVZFrgP2wdJy3dkHvrryYAEGpbMzRf7dSe7u7ttH+2/LwT8On6+qK5Nk9Od953txd5/u7lPdferiXLKMmAFWQa4D9sHSct3JJ126loABBmUXFt7r7s8l+WxVPXf00LVJPp7kliQ3jB67Ick7VxIhwBrIdcA+kOsAFrMNC+9dNOXz/rskb6qqxyT5VJL/OgcfELy1qm5M8pkkP7GaEAHWRq4D9oFcBzCvDRTts5qqyO/uDyc5dZ5D1y43HIDNkeuAfSDXAcypk5zbkSIfAAAA9ttmFtKblSIfAAAApqHIBwAAgB2hyAcAAIAdoCcfAAAAdkUnvYEb38/oxKYDAAAAAJbDSD4AAABMQ08+AAAA7AA9+QAAALBDjOQDAADAjlDkAwAAwC5oRT4AAADshE5ybvi30FPkAwAAwDSM5AMAAMCO2IIi/8SmAwAAAIDh64Nb6M26LaCq/kFV3VNVHx5tPzLpNUbyAQAAYJJOujfSk/8r3f2/TPtkRT4AAABMY8GR+XUwXR8AAACm0T37trhXVdVHq+oNVfWUSU9W5AMAAMAk3Qe30Jt1Sy6vqjOHtpsOn7aq3ldVd5xnuz7Jryf5ziTXJLk3yf86KUzT9QEAAGAa8ymddekAAAZ0SURBVI3M39/dpy58yn7ZNCepqv8zybsmPU+RDwAAAFPoc+tdeK+qruzue0e7fyvJHZNeo8gHAACAiZbWYz+L/7mqrjm4eD6d5O9MeoEiHwAAAAaou18562sU+QAAADBJZytuoVe9xukGVfUXST6T5PIk96/twtMbalyJ2OYx1LiS4ca2ibj+o+5+2pqvuVJy3ULENruhxpWIbdxO5bstyHXJcGMbalzJcGMbalyJ2MbtVK477Mknvq2/7zHXzfy69z74Lz943MJ7y7bWkfxH3uyqOrPOb3JaQ40rEds8hhpXMtzYhhrXtpHr5ie22Q01rkRsu27ouS4ZbmxDjSsZbmxDjSsR2z7pJL0FI/mm6wMAAMAk3Umvd3X9eSjyAQAAYApG8i/s9IauO8lQ40rENo+hxpUMN7ahxrWthvrzHGpcidjmMdS4ErHtiyH/LIca21DjSoYb21DjSsS2N76cL77nfefeevkcL13rughrXXgPAAAAWJ0Tmw4AAAAAWI61FvlVdV1VfaKq7qqq16zz2ueJ5Q1VdV9V3XHosadW1a1V9cnRn0/ZQFzPrKr3V9XHq+pjVfXqAcX22Kr646r6yCi2fzh6/FlVdfvofX1LVT1m3bGN4jhZVR+qqncNLK5PV9WfVtWHq+rM6LGNv5+jOC6rqrdV1Z9V1Z1V9eKhxLbN5LqpYxtkvht6rhvFMrh8J9ftH7lu6tjkuvljHFyuG8UxyHwn1/GItRX5VXUyya8l+c+TPD/JT1XV89d1/fN4Y5Lxmxy+Jslt3f2cJLeN9tft4SQ/193PT/J9SX5m9HMaQmwPJnlpd78gyTVJrquq70vyT5P8Snc/O8kXk9y4gdiS5NVJ7jy0P5S4kuQHu/uaQ7cwGcL7mSS/muR3u/t5SV6Qg5/fUGLbSnLdTIaa74ae65Lh5ju5bk/IdTOR6+Y31FyXDDPfyXUc6O61bElenOQ9h/Zfm+S167r+BWK6Oskdh/Y/keTK0ddXJvnEJuMbxfHOJC8fWmxJLk3yJ0m+NwcLSVx0vvd5jfFclYPE9dIk70pSQ4hrdO1PJ7l87LGNv59Jnpzk32W0NseQYtvmTa5bKM7B5buh5brRtQeZ7+S6/drkuoXilOumi2mQuW507cHlO7nOdnhb53T9ZyT57KH9u0ePDckV3X3v6OvPJblik8FU1dVJXpjk9gwkttG0qQ8nuS/JrUn+PMmXuvvh0VM29b7+syQ/n+SRG1d+20DiSpJO8t6q+mBV3TR6bAjv57OS/EWS3xxNhfuNqnr8QGLbZnLdHIaW7wac65Lh5ju5br/IdXOQ62Yy1FyXDDPfyXV8k4X3LqC7Owd/gTeiqp6Q5HeS/Gx3P3D42CZj6+6z3X1NDj5dfVGS520ijsOq6keT3NfdH9x0LBfwA939PTmY0vgzVfXXDh/c4Pt5UZLvSfLr3f3CJF/J2BSuTf89YPWG8B4PMd8NMdclg893ch2DNYT3WK6b3sBzXTLMfCfX8U3rLPLvSfLMQ/tXjR4bks9X1ZVJMvrzvk0EUVUX5+AfgTd199uHFNsjuvtLSd6fg6lSl1XVRaNDm3hfX5Lkb1bVp5O8OQfTun51AHElSbr7ntGf9yV5Rw7+ER3C+3l3kru7+/bR/tty8I/DEGLbZnLdDIae7waW65IB5zu5bu/IdTOQ62Y22FyXDDbfyXV80zqL/A8kec5oVczHJPnJJLes8frTuCXJDaOvb8hBz9RaVVUleX2SO7v7lwcW29Oq6rLR14/LQT/ZnTn4R+HHNhVbd7+2u6/q7qtz8Hv1e93905uOK0mq6vFV9cRHvk7yQ0nuyADez+7+XJLPVtVzRw9dm+TjQ4hty8l1UxpqvhtqrkuGm+/kur0k101JrpvdUHNdMtx8J9dxxDoXAEjyI0n+bQ76fV63zmufJ5bfSnJvkody8MnXjTno9bktySeTvC/JUzcQ1w/kYBrNR5N8eLT9yEBi++4kHxrFdkeSvz96/K8m+eMkdyX57SSXbPB9/etJ3jWUuEYxfGS0feyR3/shvJ+jOK5Jcmb0nv5fSZ4ylNi2eZPrpo5tkPluG3LdKJ7B5Du5bj83uW7q2OS6xeIcTK47FMMg851cZ3tkq9EvBAAAALDlLLwHAAAAO0KRDwAAADtCkQ8AAAA7QpEPAAAAO0KRDwAAADtCkQ8AAAA7QpEPAAAAO0KRDwAAADvi/wfYGHc2SfWstgAAAABJRU5ErkJggg==\n"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "execution_count": 44,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "BOUNDS = [\n",
        "  (0, 512), (0, 512),  # mux, muy\n",
        "  (-np.inf, np.inf), (-np.inf, np.inf), # roll, rEff\n",
        "  (0, np.inf),  # axRatio\n",
        "  (0, np.inf), # i0\n",
        "] * 2"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def get_fitted_bd_model(subject_id, progress=True):\n",
        "  psf = gu.get_psf(subject_id)\n",
        "  diff_data = gu.get_diff_data(subject_id)\n",
        "  pixel_mask = 1 - np.array(diff_data['mask'])[::-1]\n",
        "  galaxy_data = np.array(diff_data['imageData'])[::-1]\n",
        "\n",
        "  m = best_models.loc[subject_id]['Model']\n",
        "  if m['disk'] is None or m['bulge'] is None:\n",
        "    return None\n",
        "  p0 = comp_to_p(m['disk']) + comp_to_p(m['bulge'])\n",
        "  image_size = galaxy_data.shape[0]\n",
        "  cx, cy = np.mgrid[0:image_size, 0:image_size]\n",
        "  if progress:\n",
        "    with tqdm(desc='Fitting model', leave=False) as pbar:\n",
        "      def update_bar(p):\n",
        "        pbar.update(1)\n",
        "      res = minimize(\n",
        "        f, p0, callback=update_bar,\n",
        "        args=(galaxy_data, cx, cy, psf, pixel_mask),\n",
        "        bounds=BOUNDS\n",
        "      )\n",
        "  else:\n",
        "    res = minimize(\n",
        "      f, p0, callback=update_bar,\n",
        "      args=(galaxy_data, cx, cy, psf, pixel_mask)\n",
        "    )\n",
        "  fit_disk = sanitize_param_dict(\n",
        "    comp_from_p(res['x'][:6].tolist())\n",
        "  )\n",
        "  fit_bulge = sanitize_param_dict(\n",
        "    comp_from_p(res['x'][6:].tolist())\n",
        "  )\n",
        "  return fit_disk, fit_bulge, res"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "out = pd.Series([]).rename('fit_results')"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "d = UpdatableDisplay('')\n",
        "for i, subject_id in enumerate(sid_list):\n",
        "  d('{} / {} : {}'.format(i, len(sid_list), subject_id))\n",
        "  if subject_id in out.index:\n",
        "    continue\n",
        "  try:\n",
        "    out[subject_id] = get_fitted_bd_model(subject_id)\n",
        "  except KeyboardInterrupt:\n",
        "    break"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Calculate bulge to total fractions:"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def get_flux(**comp):\n",
        "  re = comp['rEff'] / 3\n",
        "  Se = comp['i0'] / 2\n",
        "  n = comp['n']\n",
        "  k = _b(n)\n",
        "  q = comp['axRatio']\n",
        "  return 2 * np.pi * re**2 * Se * np.exp(k) * n * k**(-2*n) * gamma(2*n) / q\n",
        "\n",
        "def get_bt(disk, bulge):\n",
        "  b = get_flux(n=4, **bulge)\n",
        "  t = b + get_flux(n=1, **disk)\n",
        "  return b / t"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "bt = out.dropna().apply(lambda r: get_bt(r[0], r[1]))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "st.pearsonr(bt, simard_df['__B_T_r'].loc[bt.index])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(6, 6), dpi=75)\n",
        "plt.plot(simard_df['__B_T_r'].loc[bt.index], bt, '.')\n",
        "plt.gca().add_artist(plt.Line2D((0, 1), (0, 1), color='k', alpha=0.2))\n",
        "plt.axis('equal')\n",
        "plt.xlim(0, .9)\n",
        "plt.ylim(0, .9)\n",
        "plt.tight_layout()"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(10, 8), dpi=100)\n",
        "subject_id = np.random.choice(out.dropna().index)\n",
        "psf = gu.get_psf(subject_id)\n",
        "diff_data = gu.get_diff_data(subject_id)\n",
        "pixel_mask = 1 - np.array(diff_data['mask'])[::-1]\n",
        "galaxy_data = np.array(diff_data['imageData'])[::-1]\n",
        "r = comp_to_p(out.loc[subject_id][0]) + comp_to_p(out.loc[subject_id][1])\n",
        "image_size = galaxy_data.shape[0]\n",
        "cx, cy = np.mgrid[0:image_size, 0:image_size]\n",
        "im_psf = convolve2d(\n",
        "  bulge_disk_model(r, cx, cy), psf, mode='same', boundary='symm'\n",
        ")\n",
        "plt.subplot(231)\n",
        "plt.imshow(galaxy_data, vmin=0, vmax=1)\n",
        "plt.subplot(232)\n",
        "s = mean_squared_error(galaxy_data * pixel_mask, im_psf * pixel_mask)\n",
        "plt.title('{}: {:.4f}'.format(subject_id, s))\n",
        "plt.imshow(im_psf, vmin=0, vmax=1)\n",
        "plt.subplot(233)\n",
        "full_model = rendering.calculate_model(\n",
        "    best_models.loc[subject_id]['Model'],\n",
        "    psf=gu.get_psf(subject_id),\n",
        "    image_size=gu.get_diff_data(subject_id)['width'],\n",
        ")\n",
        "s2 = mean_squared_error(galaxy_data * pixel_mask, full_model * pixel_mask / 0.8)\n",
        "plt.title(s2)\n",
        "plt.imshow(full_model, vmin=0, vmax=1)\n",
        "plt.subplot(235)\n",
        "d = galaxy_data * pixel_mask - im_psf * pixel_mask\n",
        "plt.imshow(d / 0.8, cmap='RdGy', vmin=-np.abs(d).max(), vmax=np.abs(d).max())\n",
        "plt.subplot(236)\n",
        "d2 = rendering.compare_to_galaxy(\n",
        "    full_model, galaxy_data, psf=None, pixel_mask=pixel_mask, stretch=False\n",
        ")\n",
        "plt.imshow(d2, cmap='RdGy', vmin=-np.abs(d).max(), vmax=np.abs(d).max())"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "out.loc[subject_id][2]"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "best_models.loc[20902040]"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "simard_df.loc[subject_id].keys()"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "simard_df.loc[subject_id]['__B_T_r']"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "sanitize_param_dict(fit_disk)['axRatio'], 1 / np.arctan(simard_df.loc[subject_id]['i'])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "sanitize_param_dict(fit_bulge)['axRatio'], 1 - simard_df.loc[subject_id]['e']"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    }
  ],
  "metadata": {
    "kernel_info": {
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.6.7",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "kernelspec": {
      "name": "python3",
      "language": "python",
      "display_name": "Python 3"
    },
    "nteract": {
      "version": "0.14.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}