beta_fitting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "07cd3ac9-5d20-42b2-b0e7-b74d7a173d8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# BASIC LIBRARIES\n",
    "import os, glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm, cauchy\n",
    "\n",
    "# PLOTTING LIBRARIES\n",
    "from corner import corner\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm, LinearSegmentedColormap, Normalize\n",
    "from matplotlib.patches import Ellipse\n",
    "from matplotlib import gridspec\n",
    "\n",
    "# ASTROPY\n",
    "from astropy.io import fits\n",
    "from astropy.wcs import WCS\n",
    "from astropy.nddata import Cutout2D, CCDData\n",
    "from astropy.convolution import Gaussian2DKernel as Gauss\n",
    "from astropy.convolution import convolve\n",
    "\n",
    "# SHERPA 4.14 (installed with CIAO using Conda)\n",
    "from sherpa.models.parameter import Parameter\n",
    "from sherpa_contrib.profiles import *\n",
    "from sherpa.astro.ui import *\n",
    "from sherpa.plot import *\n",
    "from sherpa.astro import xspec\n",
    "\n",
    "# IPYWIDGETS\n",
    "from IPython.display import clear_output\n",
    "from ipywidgets import interact, interactive, interactive_output, Label, Layout, IntSlider, FloatSlider, HBox, VBox, Text, FloatRangeSlider, SelectMultiple, FloatLogSlider\n",
    "import ipywidgets as widgets\n",
    "\n",
    "q1, q3 = 0.15865525, 0.84134475\n",
    "galaxies = glob.glob(\"*.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e2cf46f1-8e88-4311-a204-be2b307bfabf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def load_galaxy(galaxy, select_scale):\n",
    "    reset()\n",
    "    clean()\n",
    "    global path\n",
    "    galaxy = galaxy.replace(\".fits\",\"\")\n",
    "    path = f\"{galaxy}.fits\"\n",
    "    \n",
    "    file = fits.open(path)\n",
    "    data = file[0].data\n",
    "    shape = data.shape\n",
    "    X0 = shape[0] / 2\n",
    "    scale = int(shape[0])\n",
    "\n",
    "    load_image(path)\n",
    "    if select_scale != \"original\":\n",
    "        scale = int(select_scale.split()[0].split(\"x\")[0])\n",
    "        if scale > X0*2:\n",
    "            with output:\n",
    "                print(\"Invalid scale. Original image is smaller.\")\n",
    "        else:\n",
    "            with output:\n",
    "                clear_output()\n",
    "            notice2d(f\"BOX({X0+0.5},{X0+0.5},{scale},{scale})\")\n",
    "\n",
    "    else:\n",
    "        with output_scale:\n",
    "            clear_output()\n",
    "            print(f\"     Original scale: {scale} x {scale}\")\n",
    "            \n",
    "    set_coord(\"image\")\n",
    "    \n",
    "    # for ign in glob.glob(f\"{path}/ignore.reg\"): ignore2d(ign)\n",
    "    \n",
    "    # load_psf(\"psf0\", \"psf.fits\")\n",
    "    # set_psf(\"psf0\")\n",
    "    # psf0.size = [16,16]\n",
    "    # psf0.center = [64,64]\n",
    "    \n",
    "    return X0, shape\n",
    "\n",
    "def choose_model(select_model, X0):\n",
    "    reset()\n",
    "    if select_model == \"Single beta\":\n",
    "        m = set_source((beta2d.b1 + const2d.bkg))\n",
    "    if select_model == \"Single beta + gauss\":\n",
    "        m = set_source((beta2d.b1 + gauss2d.g1 + const2d.bkg))\n",
    "    elif select_model == \"Double beta\":\n",
    "        m = set_source((beta2d.b1 + beta2d.b2 + const2d.bkg))\n",
    "    elif select_model == \"Double beta + gauss\":\n",
    "        m = set_source((beta2d.b1 + beta2d.b2 + gauss2d.g1 + const2d.bkg))\n",
    "    elif select_model == \"Triple beta\":\n",
    "        m = set_source((beta2d.b1 + beta2d.b2 + beta2d.b3 + const2d.bkg))\n",
    "    elif select_model == \"Cuspy beta\":\n",
    "        m = set_source((beta2d.b1 + beta2d.b2 + const2d.bkg))\n",
    "        \n",
    "    set_par(b1.xpos, val=X0, min=0, max=2*X0)\n",
    "    set_par(b1.ypos, val=X0, min=0, max=2*X0)\n",
    "    set_par(b1.ellip, val=0, min=0, max=0.7)\n",
    "    set_par(b1.theta, val=0)\n",
    "    set_par(b1.r0, val=1, min=5e-2, max=1e3)\n",
    "    set_par(b1.ampl, val=50, min=1e-2, max=5e4)\n",
    "    set_par(b1.alpha, val=1, min=0.1, max=10)\n",
    "\n",
    "    if select_model in [\"Double beta\", \"Double beta + gauss\", \"Triple beta\"]:\n",
    "        set_par(b2.xpos, val=X0, min=0, max=2*X0)\n",
    "        set_par(b2.ypos, val=X0, min=0, max=2*X0)\n",
    "        set_par(b2.ellip, val=0, min=0, max=0.7)\n",
    "        set_par(b2.theta, val=0)\n",
    "        set_par(b2.r0, val=20, min=1e-1, max=1e3)\n",
    "        set_par(b2.ampl, val=20, min=1e-2, max=5e4)\n",
    "        set_par(b2.alpha, val=1, min=0.1, max=10)\n",
    "\n",
    "    if select_model in [\"Triple beta\"]:\n",
    "        set_par(b3.xpos, val=X0, min=0, max=2*X0)\n",
    "        set_par(b3.ypos, val=X0, min=0, max=2*X0)\n",
    "        set_par(b3.ellip, val=0, min=0, max=0.7)\n",
    "        set_par(b3.theta, val=0)\n",
    "        set_par(b3.r0, val=100, min=1, max=5e3)\n",
    "        set_par(b3.ampl, val=1, min=1e-2, max=5e4)\n",
    "        set_par(b3.alpha, val=1, min=0.1, max=10)\n",
    "\n",
    "    if select_model in [\"Single beta + gauss\", \"Double beta + gauss\"]:\n",
    "        set_par(g1.xpos, val=X0, min=X0-40, max=X0+40)\n",
    "        set_par(g1.ypos, val=X0, min=X0-40, max=X0+40)\n",
    "        set_par(g1.ellip, val=0, min=0, max=0.7)\n",
    "        set_par(g1.theta, val=0)\n",
    "        set_par(g1.fwhm, val=1, min=1e-3, max=2e1)\n",
    "        set_par(g1.ampl, val=50, min=1e-3, max=1e5)     \n",
    "        \n",
    "    set_par(bkg.c0, val=1e-3, min=1e-6, max=200)\n",
    "\n",
    "    guess(b1)\n",
    "\n",
    "def actualize(b):\n",
    "    if model in [\"Single beta\"]:\n",
    "        x0.value, y0.value, ellip.value, theta.value, r0.value, ampl.value, alpha.value, c0.value = b1.xpos.val, b1.ypos.val, b1.ellip.val, b1.theta.val, b1.r0.val, b1.ampl.val, b1.alpha.val, bkg.c0.val\n",
    "\n",
    "    elif model in [\"Single beta + gauss\"]:\n",
    "        x0.value, y0.value, ellip.value, theta.value, r0.value, ampl.value, alpha.value, c0.value, x0g.value, y0g.value, fwhm.value, ampl_g.value = b1.xpos.val, b1.ypos.val, b1.ellip.val, b1.theta.val, b1.r0.val, b1.ampl.val, b1.alpha.val, bkg.c0.val, g1.xpos.val, g1.ypos.val, g1.fwhm.val, g1.ampl.val\n",
    "\n",
    "    elif model in [\"Double beta\"]:\n",
    "        x0.value, y0.value, ellip.value, theta.value, r0.value, ampl.value, alpha.value, c0.value, x02.value, y02.value, ellip2.value, theta2.value, r02.value, ampl2.value, alpha2.value = b1.xpos.val, b1.ypos.val, b1.ellip.val, b1.theta.val, b1.r0.val, b1.ampl.val, b1.alpha.val, bkg.c0.val, b2.xpos.val, b2.ypos.val, b2.ellip.val, b2.theta.val, b2.r0.val, b2.ampl.val, b2.alpha.val\n",
    "\n",
    "    elif model in [\"Double beta + gauss\"]:\n",
    "        x0.value, y0.value, ellip.value, theta.value, r0.value, ampl.value, alpha.value, c0.value, x02.value, y02.value, ellip2.value, theta2.value, r02.value, ampl2.value, alpha2.value, x0g.value, y0g.value, fwhm.value, ampl_g.value = b1.xpos.val, b1.ypos.val, b1.ellip.val, b1.theta.val, b1.r0.val, b1.ampl.val, b1.alpha.val, bkg.c0.val, b2.xpos.val, b2.ypos.val, b2.ellip.val, b2.theta.val, b2.r0.val, b2.ampl.val, b2.alpha.val, g1.xpos.val, g1.ypos.val, g1.fwhm.val, g1.ampl.val\n",
    "\n",
    "    elif model in [\"Triple beta\"]:\n",
    "        x0.value, y0.value, ellip.value, theta.value, r0.value, ampl.value, alpha.value, c0.value, x02.value, y02.value, ellip2.value, theta2.value, r02.value, ampl2.value, alpha2.value, x03.value, y03.value, ellip3.value, theta3.value, r03.value, ampl3.value, alpha3.value = b1.xpos.val, b1.ypos.val, b1.ellip.val, b1.theta.val, b1.r0.val, b1.ampl.val, b1.alpha.val, bkg.c0.val, b2.xpos.val, b2.ypos.val, b2.ellip.val, b2.theta.val, b2.r0.val, b2.ampl.val, b2.alpha.val, b3.xpos.val, b3.ypos.val, b3.ellip.val, b3.theta.val, b3.r0.val, b3.ampl.val, b3.alpha.val\n",
    "\n",
    "def beta1d_model(x, ampl, rc, alpha, ellip, theta, c0):\n",
    "    f = ampl * (1 + (x/rc)**2)**(-alpha)  \n",
    "    return f + c0\n",
    "\n",
    "def gauss1d_model(x, ampl, fwhm, c0):\n",
    "    f = ampl * np.exp(-4 * np.log(2) * (x)**2 / fwhm**22)\n",
    "    return f + c0\n",
    "\n",
    "def make_plot():\n",
    "    plt.close()\n",
    "    fig = plt.figure(figsize=(5.5,5.5))\n",
    "    \n",
    "    # RADIAL PROFILE - SHERPA\n",
    "    # get_model_prof_prefs()[\"linestyle\"] = \"-\"\n",
    "    # get_data_prof_prefs()[\"xlog\"] = True\n",
    "    # get_data_prof_prefs()[\"ylog\"] = True\n",
    "    # get_delchi_prof_prefs()[\"xlog\"] = True\n",
    "    # prof_fit_delchi(model=b1)\n",
    "\n",
    "    \n",
    "    # RADIAL PROFILE\n",
    "    prof = get_data_prof(model=b1)\n",
    "    rad = (prof.xlo + prof.xhi) / 2\n",
    "    y, yerr = prof.y, prof.yerr\n",
    "    y_fit = get_model_prof(model=b1).y\n",
    "    res = (y-y_fit) / yerr\n",
    "    \n",
    "    # y_fit1 = beta1d_model(rad, b1.ampl.val, b1.r0.val, b1.alpha.val, 0, 0, bkg.c0.val)\n",
    "    # try: y_fit2 = beta1d_model(rad, b2.ampl.val, b2.r0.val, b2.alpha.val, 0, 0, bkg.c0.val)\n",
    "    # except: pass\n",
    "    # try: y_fit3 = beta1d_model(rad, b3.ampl.val, b3.r0.val, b3.alpha.val, 0, 0, bkg.c0.val)\n",
    "    # except: pass\n",
    "        \n",
    "    spec = gridspec.GridSpec(ncols=1, nrows=2, height_ratios=[3, 1], hspace=0.05)\n",
    "    ax = fig.add_subplot(spec[0])\n",
    "    ax_res = fig.add_subplot(spec[1])\n",
    "    \n",
    "    c = 0.492 * 0.492\n",
    "    ax.set_title(\"Profile\")\n",
    "    kw = {\"lw\" : 0, \"elinewidth\" : 2, \"ms\" : 4, \"marker\" : \"o\", \"markerfacecolor\" : \"none\", \"markeredgecolor\" : \"C0\"}\n",
    "    i = yerr / y < 2\n",
    "    ax.errorbar(rad[i], y[i] / c, yerr=yerr[i], zorder=1, **kw)\n",
    "    ax.step(rad, y_fit / c, where=\"mid\", zorder=2, label=\"total\")\n",
    "    \n",
    "    y_fit1 = beta1d_model(rad, b1.ampl.val, b1.r0.val, b1.alpha.val, 0, 0, 0)\n",
    "    ax.step(rad, y_fit1 / c, ls=\"--\", where=\"mid\", zorder=2, label=\"beta2d.b1\")\n",
    "    \n",
    "    # if model in [\"Single beta + gauss\", \"Double beta\", \"Double beta + gauss\", \"Triple beta\"]:\n",
    "    #     y_fit1 = beta1d_model(rad, b1.ampl.val, b1.r0.val, b1.alpha.val, 0, 0, bkg.c0.val)\n",
    "    #     ax.step(rad, y_fit1 / c, ls=\"--\", where=\"mid\", zorder=2, label=\"beta2d.b1\")\n",
    "    if model in [\"Double beta\", \"Double beta + gauss\", \"Triple beta\"]:\n",
    "        y_fit2 = beta1d_model(rad, b2.ampl.val, b2.r0.val, b2.alpha.val, 0, 0, 0)\n",
    "        ax.step(rad, y_fit2 / c, ls=\"--\", where=\"mid\", zorder=2, label=\"beta2d.b2\")\n",
    "    if model in [\"Triple beta\"]:\n",
    "        y_fit3 = beta1d_model(rad, b3.ampl.val, b3.r0.val, b3.alpha.val, 0, 0, 0)\n",
    "        ax.step(rad, y_fit3 / c, ls=\"--\", where=\"mid\", zorder=2, label=\"beta2d.b3\")\n",
    "    if model in [\"Single beta + gauss\", \"Double beta + gauss\"]:\n",
    "        y_fit_g = gauss1d_model(rad, g1.ampl.val, g1.fwhm.val, 0)\n",
    "        ax.step(rad, y_fit_g / c, ls=\"--\", where=\"mid\", zorder=2, label=\"gauss.g1\")\n",
    "    \n",
    "    ax.step(rad, np.ones_like(rad)*bkg.c0.val, ls=\"--\", where=\"mid\", zorder=2, label=\"bkg.c0\")\n",
    "    \n",
    "    ax.legend(loc=\"lower left\")\n",
    "    ax.set_xscale(\"log\")\n",
    "    ax.set_yscale(\"log\")\n",
    "    ax.set_xticklabels([])\n",
    "    ax.set_ylabel(\"Counts per pixel\")\n",
    "    ax.set_ylim(bottom=bkg.c0.val*0.1)\n",
    "    #ax.set_ylabel(\"counts / arcsec$^2$\")\n",
    "\n",
    "    # PROFILE RESIDUALS\n",
    "    ax_res.errorbar(rad[i], res[i], yerr=1, **kw, zorder=1)\n",
    "    ax_res.axhline(0, ls=\"--\", color=\"black\", zorder=2)\n",
    "    ax_res.set_xscale(\"log\")\n",
    "    ax_res.set_xlabel(\"Radius (pixel)\")\n",
    "    ax_res.set_ylabel(\"Residual ($\\sigma$)\")\n",
    "    ax_res.set_xticks([1,10,100])\n",
    "    ax_res.set_xticklabels([1,10,100])\n",
    "    \n",
    "    ax1 = fig.add_axes([1.05, 0.1, 0.8, 0.8])\n",
    "    ax1.set_title(\"Data\")\n",
    "    data = get_data().y.reshape(get_data().shape)\n",
    "    X0 = data.shape[0] / 2\n",
    "\n",
    "    # ELLIPSES AT CORE RADII\n",
    "    if scale == \"original\": SCALE = X0\n",
    "    else: SCALE = float(scale.split()[0].split(\"x\")[0]) / 2\n",
    "\n",
    "    zoom = X0 / SCALE\n",
    "    X1, Y1, R1, E1, A1 = b1.xpos.val / zoom - 0.5, b1.ypos.val / zoom - 0.5, r0.value, ellip.value, theta.value * 180 / np.pi\n",
    "    if model in [\"Double beta\", \"Double beta + gauss\", \"Triple beta\"]:\n",
    "        X2, Y2, R2, E2, A2 = b2.xpos.val / zoom - 0.5, b2.ypos.val / zoom - 0.5, r02.value, ellip2.value, theta2.value * 180 / np.pi\n",
    "    if model in [\"Triple beta\"]:\n",
    "        X3, Y3, R3, E3, A3 = b3.xpos.val / zoom - 0.5, b3.ypos.val / zoom - 0.5, r03.value, ellip3.value, theta3.value * 180 / np.pi\n",
    "    \n",
    "    ticks = np.linspace(0, 2*SCALE // 64 * 64, 5)\n",
    "\n",
    "    ellipse = Ellipse((X1,Y1), R1*2, R1*2 * (1 - E1), angle=A1, ec=\"w\", fill=False, ls=\"--\")\n",
    "    if model in [\"Double beta\", \"Double beta + gauss\", \"Triple beta\"]:\n",
    "        ellipse2 = Ellipse((X2,Y2), R2*2, R2*2 * (1 - E2), angle=A2, ec=\"w\", fill=False, ls=\"--\")\n",
    "    if model in [\"Triple beta\"]:\n",
    "        ellipse3 = Ellipse((X3,Y3), R3*2, R3*2 * (1 - E3), angle=A3, ec=\"w\", fill=False, ls=\"--\")\n",
    "    \n",
    "    size_to_std = {64 : 1, 128 : 1, 192 : 1, 256 : 1.5, 384 : 2, 512 : 3, 640 : 4, 768 : 4, 896 : 5}\n",
    "    smooth_std = size_to_std[SCALE * 2 // 64 * 64]\n",
    "    \n",
    "    data = data[int(X0-SCALE):int(X0+SCALE),int(X0-SCALE):int(X0+SCALE)]\n",
    "    data = convolve(data, boundary = \"extend\", nan_treatment=\"fill\", fill_value = np.amin(data),\n",
    "                    kernel = Gauss(x_stddev=smooth_std, y_stddev=smooth_std))\n",
    "    \n",
    "    for gal in [\"3C449\",\"IC1459\",\"IC4296\",\"IC4765\",\"NGC57\",\"NGC193\",\"NGC1553\",\"NGC4342\",\"NGC4438\",\"NGC4526\",\"NGC4555\",\"NGC4936\"]:\n",
    "        if gal in galaxy: data = np.log10(data+1)\n",
    "    ax1.imshow(np.log10(data+1), origin=\"lower\") #, norm=LogNorm())\n",
    "    ax1.plot(X1, Y1, \"+r\", ms=15)\n",
    "    ax1.set_xlabel(\"x (pixels)\")\n",
    "    ax1.set_ylabel(\"y (pixels)\")\n",
    "    ax1.set_xticks(ticks+np.array([0,0,0,0,-1]))\n",
    "    # ax1.set_xticklabels(ticks)\n",
    "    ax1.set_yticks(ticks)\n",
    "    # ax1.set_yticklabels(ticks)\n",
    "    \n",
    "    ax1.add_patch(ellipse)\n",
    "    ax1.text(X1, Y1 + R1+1, \"r0\", va=\"bottom\", ha=\"center\", color=\"w\")\n",
    "    if model in [\"Double beta\", \"Double beta + gauss\", \"Triple beta\"]:\n",
    "        if b2.r0.val < SCALE * np.sqrt(2):\n",
    "            ax1.add_patch(ellipse2)\n",
    "            ax1.text(X2, Y2 + R2+1, \"r02\", va=\"bottom\", ha=\"center\", color=\"w\")\n",
    "    if model in [\"Triple beta\"]:\n",
    "        if b3.r0.val < SCALE * np.sqrt(2):\n",
    "            ax1.add_patch(ellipse3)\n",
    "            ax1.text(X3, Y3 + R3+1, \"r03\", va=\"bottom\", ha=\"center\", color=\"w\")\n",
    "\n",
    "    # MODEL IMAGE\n",
    "    cmap = plt.cm.viridis\n",
    "    cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "    cmap = LinearSegmentedColormap.from_list(\n",
    "        'Custom cmap', cmaplist, 30)\n",
    "        \n",
    "    ax2 = fig.add_axes([0.10, -0.85, 0.8, 0.8])\n",
    "    ax2.set_title(\"Model\")\n",
    "    data = get_model_image().y.reshape(get_data().shape)\n",
    "    data = data[int(X0-SCALE):int(X0+SCALE),int(X0-SCALE):int(X0+SCALE)]\n",
    "    ax2.imshow(data+1, origin=\"lower\", norm=LogNorm(), cmap=cmap)\n",
    "    ax2.plot(X1, Y1, \"+r\", ms=15)\n",
    "    ax2.set_xlabel(\"x (pixels)\")\n",
    "    ax2.set_ylabel(\"y (pixels)\")\n",
    "    ax2.set_xticks(ticks+np.array([0,0,0,0,-1]))\n",
    "    # ax2.set_xticklabels(ticks)\n",
    "    ax2.set_yticks(ticks)\n",
    "    # ax2.set_yticklabels(ticks)\n",
    "\n",
    "    # RESIDUAL IMAGE\n",
    "    data = get_resid_image().y.reshape(get_data().shape)\n",
    "    data = data[int(X0-SCALE):int(X0+SCALE),int(X0-SCALE):int(X0+SCALE)]\n",
    "    \n",
    "    # data -= np.median(data)\n",
    "    data = convolve(data, boundary = \"extend\", nan_treatment=\"fill\", fill_value = np.amin(data),\n",
    "                    kernel = Gauss(x_stddev=smooth_std, y_stddev=smooth_std))\n",
    "    \n",
    "    array = data.flatten()\n",
    "    x = np.linspace(min(array), max(array))\n",
    "    mean, std = np.median(array), np.std(array) * 4    \n",
    "    data = np.where(data > std, std, data)\n",
    "    data = np.where(data < -std, -std, data)\n",
    "    \n",
    "    ax3 = fig.add_axes([1.05, -0.85, 0.8, 0.8])\n",
    "    ax3.set_title(\"Residual\")\n",
    "    ax3.imshow(data, origin=\"lower\", cmap=\"coolwarm\") #, norm=Normalize(vmin=-std, vmax=std))\n",
    "    ax3.plot(X1, Y1, \"+r\", ms=15)\n",
    "    ax3.set_xlabel(\"x (pixels)\")\n",
    "    ax3.set_ylabel(\"y (pixels)\")\n",
    "    ax3.set_xticks(ticks+np.array([0,0,0,0,-1]))\n",
    "    # ax3.set_xticklabels(ticks)\n",
    "    ax3.set_yticks(ticks)\n",
    "    # ax3.set_yticklabels(ticks)\n",
    "    \n",
    "    # HISTOGRAM OF REDISUAL IMAGE\n",
    "    # ax4 = fig.add_axes([1.61, -0.25, 0.2, 0.14])\n",
    "    # ax4.set_title(f\"{mean:.3f} +/- {std:.3f}\", color=\"white\")\n",
    "    # ax4.hist(array, bins=bins, density=True, ec=\"black\")\n",
    "    # ax4.set_yscale(\"log\")\n",
    "    # ax4.set_ylim(*ax4.get_ylim())\n",
    "    # ax4.plot(x, dist(mean, std).pdf(x))\n",
    "    # ax4.set_yticks([])\n",
    "    # ax4.tick_params(axis='x', colors='white')\n",
    "    # ax4.tick_params(axis='y', colors='white')\n",
    "    # for spine in ax4.spines.values():\n",
    "    #     spine.set_edgecolor('white')\n",
    "    # ax4.set_facecolor('lightgray')\n",
    "    \n",
    "    plt.show()\n",
    "        \n",
    "def save_current_model(b):  \n",
    "    if scale == \"original\": SCALE = scale\n",
    "    else: SCALE = scale.split()[0]\n",
    "    with open(f\"Beta_models/Prefits/{galaxy}.txt\", \"w\") as file:\n",
    "        print(get_model(), file=file)\n",
    "    \n",
    "    with output:\n",
    "        print(\"\\nModel saved to\", f\"{galaxy}.txt\")\n",
    "\n",
    "def save_current_resid(b):\n",
    "    if scale == \"original\": SCALE = scale\n",
    "    else: SCALE = scale.split()[0]\n",
    "    \n",
    "    file = fits.open(path)\n",
    "    wcs = WCS(file[0].header)\n",
    "    \n",
    "    resid = get_resid_image().y.reshape(get_data().shape)\n",
    "    ccd = CCDData(resid, unit=\"adu\", wcs=wcs)\n",
    "    ccd.write(f\"Beta_models/Residuals/{galaxy}_resid.fits\", overwrite=True)\n",
    "    \n",
    "    with output:\n",
    "        print(\"\\nResidual saved to\", f\"{galaxy}_resid.fits\")\n",
    "        \n",
    "def fit_current_model(b):\n",
    "    with output:\n",
    "        clear_output()\n",
    "        \n",
    "        print(\"Making a fit. This might take time.\\n\")\n",
    "        fit()\n",
    "\n",
    "def texg(x, xp, xm, n):\n",
    "    r = -(int(np.log10(xm)) - n)\n",
    "    return \"${0}_{{-{1}}}^{{+{2}}}$\".format(round(x,r), round(xm,r), round(xp,r))\n",
    "\n",
    "def run_MCMC(b):\n",
    "    with output:\n",
    "        covar()\n",
    "\n",
    "        res = get_covar_results()\n",
    "\n",
    "        set_sampler(\"metropolismh\")\n",
    "        set_sampler_opt('defaultprior', True)\n",
    "        \n",
    "        clear_output()\n",
    "        print(f\"Runing MCMC simulation with {length} iterations.\\nFirst {burn} iterations will be burned.\\n\")\n",
    "    \n",
    "    stats, accept, params = get_draws(1, niter=length)\n",
    "\n",
    "    burned = stats[:burn]\n",
    "    stats, accept, params = stats[burn:], accept[burn:], params[:,burn:]\n",
    "\n",
    "    df = pd.DataFrame(data=params.T, columns=res.parnames)\n",
    "    cols = df.columns\n",
    "    size = 12 * (5/len(cols))**0.7\n",
    "\n",
    "    table = df.describe(percentiles=[q1, 0.5, q3]).T[[\"15.9%\", \"50%\", \"84.1%\"]]\n",
    "    table[\"Minus\"] = table[\"50%\"] - table[\"15.9%\"]\n",
    "    table[\"Plus\"] = table[\"84.1%\"] - table[\"50%\"]\n",
    "    table[\"Median\"] = table[\"50%\"]\n",
    "    table = table.round(3)\n",
    "    table[\"Param\"] = table.index\n",
    "    table = table[[\"Param\", \"Median\", \"Minus\", \"Plus\"]]\n",
    "    \n",
    "    df.to_csv(f\"{galaxy}_chain.csv\", index=False)\n",
    "        \n",
    "    with output:\n",
    "        # CORNER PLOT\n",
    "        fig = plt.figure(figsize=(8,8))\n",
    "        \n",
    "        fig = corner(df, color=\"C0\", fig=fig, #quantiles=[0.16, 0.5, 0.84],\n",
    "                     labels=cols, #label_kwargs={\"fontsize\" : size-2}, \n",
    "                     labelpad=0.09,\n",
    "                     show_titles=True) #, title_kwargs={\"fontsize\": size-2})\n",
    "\n",
    "        # QUANTILES\n",
    "        axes = np.array(fig.axes).reshape((len(cols), len(cols)))\n",
    "        for i, col in enumerate(df.columns):\n",
    "            ax = axes[i, i]\n",
    "            ax.set_title(texg(table.Median.loc[col], table.Plus.loc[col], table.Minus.loc[col], 2), size=size)\n",
    "            ax.axvline(np.quantile(df[col], q1), color=\"C2\", ls=\"--\", lw=2.5)\n",
    "            ax.axvline(np.quantile(df[col], 0.5), color=\"C1\", ls=\"--\", lw=2.5)\n",
    "            ax.axvline(np.quantile(df[col], q3), color=\"C2\", ls=\"--\", lw=2.5)\n",
    "\n",
    "        dof = get_stat_info()[0].dof\n",
    "        ax1 = fig.add_axes([0.33, 0.871, 0.30, 0.1])\n",
    "        ax1.set_title(\"Trace\")\n",
    "        ax1.plot(stats/dof, c=\"C2\")\n",
    "        ax1.plot(stats[:burn]/dof, c=\"C1\", label=\"burned\")\n",
    "        #ax1.set_xlabel(\"Iteration\")\n",
    "        ax1.set_ylabel(statistic, labelpad=2)\n",
    "        #ax1.set_yticklabels([])\n",
    "        ax1.legend(frameon=False)\n",
    "        \n",
    "        axt = fig.add_axes([0.64, 0.97, 0.34, 0.001])\n",
    "        axt.table(cellText=table.values, colLabels=table.columns, fontsize=size+1,\n",
    "                  loc='bottom', cellLoc=\"right\", colLoc=\"center\").scale(1, 2.5)\n",
    "        \n",
    "            \n",
    "        plt.show()\n",
    "        fig.savefig(f\"{galaxy}_chain.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a003ab8a-9a3e-49e5-8593-b798d0f0cabd",
   "metadata": {},
   "outputs": [],
   "source": [
    "layout = Layout(width='350px')\n",
    "Float = lambda x,xm,xp,s: FloatSlider(value=x, min=xm, max=xp, step=round((xp-xm)/50,7), layout=layout, continuous_update=False, description=s)\n",
    "Log = lambda x,xm,xp,s: FloatLogSlider(value=x, base=10, min=xm, max=xp, step=round((xp-xm)/50,7), layout=layout, continuous_update=False, description=s)\n",
    "\n",
    "def button(string): \n",
    "    return widgets.Button(description=string,disabled=False,\n",
    "                          value=False,\n",
    "                          button_style='',\n",
    "                          tooltip=string,\n",
    "                          icon='check')\n",
    "\n",
    "def box(string, value=False, width='65px'): \n",
    "    return widgets.Checkbox(value=value, \n",
    "                            description=string, \n",
    "                            disabled=False, \n",
    "                            indent=False, \n",
    "                            layout=Layout(width=width))\n",
    "\n",
    "select_model = widgets.Dropdown(\n",
    "    options=['Single beta', 'Single beta + gauss', 'Double beta', 'Double beta + gauss', 'Triple beta'],\n",
    "    value='Single beta', \n",
    "    description='Model:', \n",
    "    disabled=False)\n",
    "\n",
    "select_galaxy = widgets.Dropdown(\n",
    "    options=galaxies,\n",
    "    value=galaxies[0], #\"NGC4649_512\",\n",
    "    description=\"Galaxy:\",\n",
    "    disabled=False)\n",
    "\n",
    "select_scale = widgets.Dropdown(\n",
    "    options=[\"original\", \"128 x 128 pixels\", \"256 x 256 pixels\", \"384 x 384 pixels\", \"512 x 512 pixels\", \"640 x 640 pixels\", \"768 x 768 pixels\"],\n",
    "    value=\"original\",\n",
    "    description=\"Scale:\",\n",
    "    disabled=False)\n",
    "\n",
    "select_method = widgets.Dropdown(\n",
    "    options=[\"levmar\", \"neldermead\", \"moncar\"],\n",
    "    value=\"levmar\",\n",
    "    description=\"Method:\",\n",
    "    disabled=False)\n",
    "\n",
    "select_statistic = widgets.Dropdown(\n",
    "    options=[\"cstat\", \"cash\", \"chi2\", \"chi2xspecvar\"],\n",
    "    value=\"cstat\",\n",
    "    description=\"Statistic:\",\n",
    "    disabled=False)\n",
    "\n",
    "button_layout = Layout(width='150px')\n",
    "\n",
    "save_button = widgets.Button(\n",
    "    description=\"Save to txt\",\n",
    "    value=False,\n",
    "    layout=button_layout)\n",
    "\n",
    "save_resid_button = widgets.Button(\n",
    "    description=\"Save residual\",\n",
    "    value=False,\n",
    "    layout=button_layout)\n",
    "\n",
    "length_text = widgets.Text(\n",
    "    value='1000',\n",
    "    placeholder='100000',\n",
    "    description='Chain length:',\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "burn_text = widgets.Text(\n",
    "    value='100',\n",
    "    placeholder='10000',\n",
    "    description='Chain burn:',\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "beta1 = widgets.HTML(value=\"<h3 style='text-align:center;margin:auto'>Beta2d.b1</h3>\", description='')\n",
    "beta2 = widgets.HTML(value=\"<h3 style='text-align:center;margin:auto'>Beta2d.b2</h3>\", description='')\n",
    "beta3 = widgets.HTML(value=\"<h3 style='text-align:center;margin:auto'>Beta2d.b3</h3>\", description='')\n",
    "gauss1 = widgets.HTML(value=\"<h3 style='text-align:center;margin:auto'>Gauss2d.g1</h3>\", description='')\n",
    "bkg1 = widgets.HTML(value=\"<h3 style='text-align:center;margin:auto'>Const2d.bkg</h3>\", description='')\n",
    "\n",
    "def h_blank(h=16):\n",
    "    blank = HBox([])\n",
    "    blank.layout = Layout(height=f\"{h}px\")\n",
    "    return blank\n",
    "\n",
    "def v_blank(h=50):\n",
    "    blank = VBox([])\n",
    "    blank.layout = Layout(width=f\"{h}px\")\n",
    "    return blank\n",
    "\n",
    "output = widgets.Output()\n",
    "\n",
    "output2 = widgets.Output()\n",
    "\n",
    "output3 = widgets.Output()\n",
    "\n",
    "output_scale = widgets.Output()\n",
    "\n",
    "output.layout = Layout(width='700px')\n",
    "\n",
    "output_scale.layout = Layout(width='265px')\n",
    "\n",
    "fit_button = widgets.Button(\n",
    "    description=\"Fit\",\n",
    "    value=False,\n",
    "    button_style='success',\n",
    "    layout=button_layout)\n",
    "\n",
    "MCMC_button = widgets.Button(\n",
    "    description=\"Run MCMC\",\n",
    "    value=False,\n",
    "    layout=button_layout)\n",
    "\n",
    "# ON-CLICK BUTTONS\n",
    "save_button.on_click(save_current_model)\n",
    "\n",
    "save_resid_button.on_click(save_current_resid)\n",
    "\n",
    "fit_button.on_click(fit_current_model)\n",
    "fit_button.on_click(actualize)\n",
    "\n",
    "MCMC_button.on_click(run_MCMC)\n",
    "\n",
    "def scale_to_original(change):\n",
    "    select_scale.value=\"original\"\n",
    "    \n",
    "def model_to_single(change):\n",
    "    select_model.value=\"Single beta\"\n",
    "\n",
    "def method_to_levmar(change):\n",
    "    select_method.value=\"levmar\"\n",
    "    \n",
    "def statistic_to_cstat(change):\n",
    "    select_statistic.value=\"cstat\"\n",
    "    \n",
    "select_galaxy.observe(scale_to_original)\n",
    "\n",
    "select_galaxy.observe(model_to_single)\n",
    "\n",
    "select_galaxy.observe(method_to_levmar)\n",
    "\n",
    "select_galaxy.observe(statistic_to_cstat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ceb775d6-5203-4c22-b2b6-76240616f66b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d6bb20b93f8042afbbe903af27d04472",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(VBox(layout=Layout(width='0px')), VBox(children=(VBox(children=(HBox(children=(Dropdown(descrip…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def main(select_galaxy, select_model, select_scale):\n",
    "    with output:\n",
    "        clear_output()\n",
    "    \n",
    "    X0, shape = load_galaxy(select_galaxy, select_scale)\n",
    "    \n",
    "    global galaxy\n",
    "    galaxy = select_galaxy\n",
    "    \n",
    "    global scale\n",
    "    scale = select_scale\n",
    "    \n",
    "    global model\n",
    "    model = select_model\n",
    "\n",
    "    choose_model(select_model, X0)\n",
    "\n",
    "    # CREATE SLIDERS & BUTTONS\n",
    "    global x0, y0, ellip, theta, r0, ampl, alpha, c0\n",
    "    x0 = Float(b1.xpos.val, b1.xpos.min, b1.xpos.max, \"X0\")\n",
    "    y0 = Float(b1.ypos.val, b1.ypos.min, b1.ypos.max, \"Y0\")\n",
    "    x0_f = box(\"freeze\"); y0_f = box(\"freeze\")\n",
    "\n",
    "    ellip = Float(b1.ellip.val, b1.ellip.min, b1.ellip.max, \"ellip\")\n",
    "    theta = Float(b1.theta.val, b1.theta.min, b1.theta.max, \"theta\")\n",
    "    ellip_f = box(\"freeze\", value=True); theta_f = box(\"freeze\", value=True); \n",
    "\n",
    "    r0 = Log(np.log10(b1.r0.val), np.log10(b1.r0.min), np.log10(b1.r0.max), \"r0\")\n",
    "    ampl = Log(b1.ampl.val, np.log10(b1.ampl.min), np.log10(b1.ampl.max), \"ampl\")\n",
    "    alpha = Log(b1.alpha.val, np.log10(b1.alpha.min), np.log10(b1.alpha.max), \"alpha\")\n",
    "    r0_f = box(\"freeze\"); ampl_f = box(\"freeze\"); alpha_f = box(\"freeze\")\n",
    "\n",
    "    if select_model in [\"Double beta\", \"Double beta + gauss\", \"Triple beta\"]:\n",
    "        global x02, y02, ellip2, theta2, r02, ampl2, alpha2\n",
    "        x02 = Float(b2.xpos.val, b2.xpos.min, b2.xpos.max, \"X0_2\")\n",
    "        y02 = Float(b2.ypos.val, b2.ypos.min, b2.ypos.max, \"Y0_2\")\n",
    "        x02_f = box(\"freeze\"); y02_f = box(\"freeze\")\n",
    "        x02_l = box(\"link\", value=True, width='45px'); y02_l = box(\"link\", value=True, width='45px')\n",
    "\n",
    "        ellip2 = Float(b2.ellip.val, b2.ellip.min, b2.ellip.max, \"ellip_2\")\n",
    "        theta2 = Float(b2.theta.val, b2.theta.min, b2.theta.max, \"theta_2\")\n",
    "        ellip2_f = box(\"freeze\", value=True); ellip2_l = box(\"link\", value=True, width='45px')\n",
    "        theta2_f = box(\"freeze\", value=True); theta2_l = box(\"link\", value=True, width='45px')\n",
    "\n",
    "        r02 = Log(b2.r0.val, np.log10(b2.r0.min), np.log10(b2.r0.max), \"r0_2\")\n",
    "        ampl2 = Log(b2.ampl.val, np.log10(b2.ampl.min), np.log10(b2.ampl.max), \"ampl_2\")\n",
    "        alpha2 = Log(b2.alpha.val, np.log10(b2.alpha.min), np.log10(b2.alpha.max), \"alpha_2\")\n",
    "        r02_f = box(\"freeze\"); ampl2_f = box(\"freeze\"); alpha2_f = box(\"freeze\"); alpha2_l = box(\"link\", value=True, width='45px')\n",
    "\n",
    "    if select_model in [\"Triple beta\"]:\n",
    "        global x03, y03, ellip3, theta3, r03, ampl3, alpha3\n",
    "        x03 = Float(b3.xpos.val, b3.xpos.min, b3.xpos.max, \"X0_3\")\n",
    "        y03 = Float(b3.ypos.val, b3.ypos.min, b3.ypos.max, \"Y0_3\")\n",
    "        x03_f = box(\"freeze\"); y03_f = box(\"freeze\")\n",
    "        x03_l = box(\"link\", value=True, width='45px'); y03_l = box(\"link\", value=True, width='45px')\n",
    "\n",
    "        ellip3 = Float(b3.ellip.val, b3.ellip.min, b3.ellip.max, \"ellip_3\")\n",
    "        theta3 = Float(b3.theta.val, b3.theta.min, b3.theta.max, \"theta_3\")\n",
    "        ellip3_f = box(\"freeze\", value=True); ellip3_l = box(\"link\", value=True, width='45px')\n",
    "        theta3_f = box(\"freeze\", value=True); theta3_l = box(\"link\", value=True, width='45px')\n",
    "        \n",
    "        r03 = Log(b3.r0.val, np.log10(b3.r0.min), np.log10(b3.r0.max), \"r0_3\")\n",
    "        ampl3 = Log(b3.ampl.val, np.log10(b3.ampl.min), np.log10(b3.ampl.max), \"ampl_3\")\n",
    "        alpha3 = Log(b3.alpha.val, np.log10(b3.alpha.min), np.log10(b3.alpha.max), \"alpha_3\")\n",
    "        r03_f = box(\"freeze\"); ampl3_f = box(\"freeze\"); alpha3_f = box(\"freeze\"); alpha3_l = box(\"link\", value=True, width='45px')\n",
    "    \n",
    "    if select_model in [\"Single beta + gauss\", \"Double beta + gauss\"]:\n",
    "        global x0g, y0g, fwhm, ampl_g\n",
    "        x0g = Float(g1.xpos.val, g1.xpos.min, g1.xpos.max, \"X0_g\")\n",
    "        y0g = Float(g1.ypos.val, g1.ypos.min, g1.ypos.max, \"Y0_g\")\n",
    "        x0g_f = box(\"freeze\"); y0g_f = box(\"freeze\")\n",
    "        x0g_l = box(\"link\", value=True, width='45px'); y0g_l = box(\"link\", value=True, width='45px')\n",
    "\n",
    "        fwhm = Log(g1.fwhm.val, np.log10(g1.fwhm.min), np.log10(g1.fwhm.max), \"fwhm\")\n",
    "        ampl_g = Log(g1.ampl.val, np.log10(g1.ampl.min), np.log10(g1.ampl.max), \"ampl_g\")\n",
    "        fwhm_f = box(\"freeze\"); ampl_g_f = box(\"freeze\")\n",
    "\n",
    "    c0 = Log(bkg.c0.val, np.log10(bkg.c0.min), np.log10(bkg.c0.max), \"c0\")\n",
    "    c0_f = box(\"freeze\")\n",
    "    \n",
    "\n",
    "    # CHOOSING PARAMS AND SLIDERS\n",
    "    params = {\"length_text\" : length_text, \"burn_text\" : burn_text, \"select_statistic\" : select_statistic, \"select_method\" : select_method,\n",
    "              \"x0\" : x0, \"x0_f\" : x0_f, \"y0\" : y0, \"y0_f\" : y0_f, \"c0\" : c0, \"c0_f\" : c0_f,\n",
    "              \"ellip\" : ellip, \"ellip_f\" : ellip_f, \"theta\" : theta, \"theta_f\" : theta_f,\n",
    "              \"r0\" : r0, \"r0_f\" : r0_f, \"ampl\" : ampl, \"ampl_f\" : ampl_f, \"alpha\" : alpha, \"alpha_f\" : alpha_f}\n",
    "    \n",
    "    buttons = VBox([h_blank(16),\n",
    "                    HBox([select_method, length_text]),\n",
    "                    HBox([select_statistic, burn_text]),\n",
    "                    h_blank(16),\n",
    "                    HBox([fit_button, MCMC_button, save_button, save_resid_button])]) #, plot_button])])\n",
    "    \n",
    "    single_beta = VBox([beta1,\n",
    "                        HBox([x0, x0_f]), \n",
    "                        HBox([y0, y0_f]),\n",
    "                        HBox([ellip, ellip_f]),\n",
    "                        HBox([theta, theta_f]),\n",
    "                        HBox([r0, r0_f]),\n",
    "                        HBox([ampl, ampl_f]), \n",
    "                        HBox([alpha, alpha_f])])\n",
    "    \n",
    "    background = VBox([bkg1,\n",
    "                       HBox([c0, c0_f])])\n",
    "                       \n",
    "    \n",
    "    if select_model in [\"Double beta\", \"Double beta + gauss\", \"Triple beta\"]:\n",
    "        double_beta = VBox([beta2,\n",
    "                            HBox([x02, x02_f, x02_l]), \n",
    "                            HBox([y02, y02_f, y02_l]),\n",
    "                            HBox([ellip2, ellip2_f, ellip2_l]),\n",
    "                            HBox([theta2, theta2_f, theta2_l]),\n",
    "                            HBox([r02, r02_f]), \n",
    "                            HBox([ampl2, ampl2_f]), \n",
    "                            HBox([alpha2, alpha2_f, alpha2_l])])\n",
    "    \n",
    "    if select_model in [\"Triple beta\"]:\n",
    "        triple_beta = VBox([beta3,\n",
    "                            HBox([x03, x03_f, x03_l]), \n",
    "                            HBox([y03, y03_f, y03_l]),\n",
    "                            HBox([ellip3, ellip3_f, ellip3_l]),\n",
    "                            HBox([theta3, theta3_f, theta3_l]),\n",
    "                            HBox([r03, r03_f]), \n",
    "                            HBox([ampl3, ampl3_f]), \n",
    "                            HBox([alpha3, alpha3_f, alpha3_l])])\n",
    "\n",
    "    if select_model in [\"Single beta + gauss\", \"Double beta + gauss\"]:\n",
    "        gauss = VBox([gauss1,\n",
    "                      HBox([x0g, x0g_f, x0g_l]), \n",
    "                      HBox([y0g, y0g_f, y0g_l]),\n",
    "                      HBox([fwhm, fwhm_f]),\n",
    "                      HBox([ampl_g, ampl_g_f])])\n",
    "    \n",
    "    if select_model in [\"Single beta\"]:\n",
    "        sliders = HBox([single_beta, background])\n",
    "        \n",
    "    elif select_model in [\"Single beta + gauss\"]:\n",
    "        sliders = HBox([single_beta, gauss, background])\n",
    "        \n",
    "        params = {**params, \"x0g\" : x0g, \"x0g_f\" : x0g_f, \"x0g_l\" : x0g_l, \"y0g\" : y0g, \"y0g_f\" : y0g_f,  \"y0g_l\" : y0g_l,\n",
    "                  \"fwhm\" : fwhm, \"fwhm_f\" : fwhm_f, \"ampl_g\" : ampl_g, \"ampl_g_f\" : ampl_g_f}\n",
    "            \n",
    "    elif select_model in [\"Double beta\"]:\n",
    "        sliders = HBox([single_beta, double_beta, background])\n",
    "        \n",
    "        params = {**params, \"x02\" : x02, \"x02_f\" : x02_f, \"x02_l\" : x02_l, \"y02\" : y02, \"y02_f\" : y02_f,  \"y02_l\" : y02_l,\n",
    "                  \"ellip2\" : ellip2, \"ellip2_f\" : ellip2_f, \"ellip2_l\" : ellip2_l, \"theta2\" : theta2, \"theta2_f\" : theta2_f, \"theta2_l\" : theta2_l,\n",
    "                  \"r02\" : r02, \"r02_f\" : r02_f, \"ampl2\" : ampl2, \"ampl2_f\" : ampl2_f, \"alpha2\" : alpha2, \"alpha2_f\" : alpha2_f, \"alpha2_l\" : alpha2_l}\n",
    "\n",
    "    elif select_model in [\"Double beta + gauss\"]:\n",
    "        sliders = HBox([single_beta, double_beta, gauss, background])\n",
    "        \n",
    "        params = {**params, \"x02\" : x02, \"x02_f\" : x02_f, \"x02_l\" : x02_l, \"y02\" : y02, \"y02_f\" : y02_f,  \"y02_l\" : y02_l,\n",
    "                  \"ellip2\" : ellip2, \"ellip2_f\" : ellip2_f, \"ellip2_l\" : ellip2_l, \"theta2\" : theta2, \"theta2_f\" : theta2_f, \"theta2_l\" : theta2_l,\n",
    "                  \"r02\" : r02, \"r02_f\" : r02_f, \"ampl2\" : ampl2, \"ampl2_f\" : ampl2_f, \"alpha2\" : alpha2, \"alpha2_f\" : alpha2_f, \"alpha2_l\" : alpha2_l,\n",
    "                  \"x0g\" : x0g, \"x0g_f\" : x0g_f, \"x0g_l\" : x0g_l, \"y0g\" : y0g, \"y0g_f\" : y0g_f,  \"y0g_l\" : y0g_l,\n",
    "                  \"fwhm\" : fwhm, \"fwhm_f\" : fwhm_f, \"ampl_g\" : ampl_g, \"ampl_g_f\" : ampl_g_f}\n",
    "        \n",
    "    elif select_model in [\"Triple beta\"]:\n",
    "        sliders = HBox([single_beta, double_beta, triple_beta, background])\n",
    "\n",
    "        params = {**params, \"x02\" : x02, \"x02_f\" : x02_f, \"x02_l\" : x02_l, \"y02\" : y02, \"y02_f\" : y02_f,  \"y02_l\" : y02_l,\n",
    "                  \"ellip2\" : ellip2, \"ellip2_f\" : ellip2_f, \"ellip2_l\" : ellip2_l, \"theta2\" : theta2, \"theta2_f\" : theta2_f, \"theta2_l\" : theta2_l,\n",
    "                  \"r02\" : r02, \"r02_f\" : r02_f, \"ampl2\" : ampl2, \"ampl2_f\" : ampl2_f, \"alpha2\" : alpha2, \"alpha2_f\" : alpha2_f, \"alpha2_l\" : alpha2_l,\n",
    "                  \"x03\" : x03, \"x03_f\" : x03_f, \"x03_l\" : x03_l, \"y03\" : y03, \"y03_f\" : y03_f,  \"y03_l\" : y03_l,\n",
    "                  \"ellip3\" : ellip3, \"ellip3_f\" : ellip3_f, \"ellip3_l\" : ellip3_l, \"theta3\" : theta3, \"theta3_f\" : theta3_f, \"theta3_l\" : theta3_l,\n",
    "                  \"r03\" : r03, \"r03_f\" : r03_f, \"ampl3\" : ampl3, \"ampl3_f\" : ampl3_f, \"alpha3\" : alpha3, \"alpha3_f\" : alpha3_f, \"alpha3_l\" : alpha3_l}\n",
    "    \n",
    "    out = interactive_output(change_params, params)\n",
    "            \n",
    "    display(VBox([sliders, VBox([HBox([VBox([buttons, output]), out]), output2])]))\n",
    "\n",
    "def change_params(select_statistic=\"cstat\", select_method=\"levmar\", length_text=1000, burn_text=100,\n",
    "                  x0=0, x0_f=False, y0=0, y0_f=False, c0=0, c0_f=False,\n",
    "                  ellip=0, ellip_f=False, theta=0, theta_f=False,\n",
    "                  r0=0, r0_f=False, ampl=0, ampl_f=False, alpha=0, alpha_f=False,\n",
    "                  x02=0, x02_f=False, x02_l=False, y02=0, y02_f=False, y02_l=False,\n",
    "                  ellip2=0, ellip2_f=False, ellip2_l=False, theta2=0, theta2_f=False, theta2_l=False,\n",
    "                  r02=0, r02_f=False, ampl2=0, ampl2_f=False, alpha2=0, alpha2_f=False, alpha2_l=False,\n",
    "                  x03=0, x03_f=False, x03_l=False, y03=0, y03_f=False, y03_l=False,\n",
    "                  ellip3=0, ellip3_f=False, ellip3_l=False, theta3=0, theta3_f=False, theta3_l=False,\n",
    "                  r03=0, r03_f=False, ampl3=0, ampl3_f=False, alpha3=0, alpha3_f=False, alpha3_l=False,\n",
    "                  x0g=0, x0g_f=False, x0g_l=False, y0g=0, y0g_f=False, y0g_l=False,\n",
    "                  fwhm=0, fwhm_f=False, ampl_g=0, ampl_g_f=False):\n",
    "         \n",
    "    set_stat(select_statistic)\n",
    "    set_method(select_method)\n",
    "    \n",
    "    global statistic\n",
    "    statistic = select_statistic\n",
    "        \n",
    "    global length, burn\n",
    "    length, burn = int(length_text), int(burn_text)\n",
    "\n",
    "    # SINGLE BETA\n",
    "    if x0_f: freeze(b1.xpos)\n",
    "    else: \n",
    "        thaw(b1.xpos)\n",
    "        set_par(b1.xpos, val=x0)\n",
    "    \n",
    "    if y0_f: freeze(b1.ypos)\n",
    "    else: \n",
    "        thaw(b1.ypos)\n",
    "        set_par(b1.ypos, val=y0)\n",
    "    \n",
    "    if ellip_f: freeze(b1.ellip)\n",
    "    else: \n",
    "        thaw(b1.ellip)\n",
    "        set_par(b1.ellip, val=ellip)\n",
    "        \n",
    "    if theta_f: freeze(b1.theta)\n",
    "    else: \n",
    "        thaw(b1.theta)\n",
    "        set_par(b1.theta, val=theta)\n",
    "    \n",
    "    if r0_f: freeze(b1.r0)\n",
    "    else: \n",
    "        thaw(b1.r0)\n",
    "        set_par(b1.r0, val=r0)\n",
    "    \n",
    "    if ampl_f: freeze(b1.ampl)\n",
    "    else: \n",
    "        thaw(b1.ampl)\n",
    "        set_par(b1.ampl, val=ampl)\n",
    "    \n",
    "    if alpha_f: freeze(b1.alpha)\n",
    "    else: \n",
    "        thaw(b1.alpha)\n",
    "        set_par(b1.alpha, val=alpha)\n",
    "\n",
    "    # DOUBLE BETA\n",
    "    try:\n",
    "        if r02_f: freeze(b2.r0)\n",
    "        else:\n",
    "            thaw(b2.r0)\n",
    "            set_par(b2.r0, val=r02)\n",
    "\n",
    "        if ampl2_f: freeze(b2.ampl)\n",
    "        else:\n",
    "            thaw(b2.ampl)\n",
    "            set_par(b2.ampl, val=ampl2)\n",
    "\n",
    "        if alpha2_l: \n",
    "            b2.alpha = b1.alpha\n",
    "            #alpha2.value = b1.alpha.val\n",
    "        else:\n",
    "            unlink(b2.alpha)\n",
    "            if alpha2_f: freeze(b2.alpha)\n",
    "            else:\n",
    "                thaw(b2.alpha)\n",
    "                set_par(b2.alpha, val=alpha2)\n",
    "\n",
    "        if ellip2_l: \n",
    "            b2.ellip = b1.ellip\n",
    "            #ellip2.value = b1.ellip.val\n",
    "        else:\n",
    "            unlink(b2.ellip)\n",
    "            if ellip2_f: freeze(b2.ellip)\n",
    "            else: \n",
    "                thaw(b2.ellip)\n",
    "                set_par(b2.ellip, val=ellip2)\n",
    "\n",
    "        if theta2_l: \n",
    "            b2.theta = b1.theta\n",
    "            #theta2.value = b1.theta.val\n",
    "        else:\n",
    "            unlink(b2.theta)\n",
    "            if theta2_f: freeze(b2.theta)\n",
    "            else: \n",
    "                thaw(b2.theta)\n",
    "                set_par(b2.theta, val=theta2)\n",
    "\n",
    "        if x02_l: \n",
    "            b2.xpos = b1.xpos\n",
    "            #x02.vlaue = b1.xpos.val\n",
    "        else:\n",
    "            unlink(b2.xpos)\n",
    "            if x02_f: freeze(b2.xpos)\n",
    "            else: \n",
    "                thaw(b2.xpos)\n",
    "                set_par(b2.xpos, val=x02)\n",
    "\n",
    "        if y02_l: \n",
    "            b2.ypos = b1.ypos\n",
    "            #y02.value = b1.ypos.val\n",
    "        else:\n",
    "            unlink(b2.ypos)\n",
    "            if y02_f: freeze(b2.ypos)\n",
    "            else: \n",
    "                thaw(b2.ypos)\n",
    "                set_par(b2.ypos, val=y02)\n",
    "    except: pass\n",
    "\n",
    "    # TRIPLE BETA\n",
    "    try:\n",
    "        if r03_f: freeze(b3.r0)\n",
    "        else:\n",
    "            thaw(b3.r0)\n",
    "            set_par(b3.r0, val=r03)\n",
    "\n",
    "        if ampl3_f: freeze(b3.ampl)\n",
    "        else:\n",
    "            thaw(b3.ampl)\n",
    "            set_par(b3.ampl, val=ampl3)\n",
    "\n",
    "        if alpha3_l: \n",
    "            b3.alpha = b1.alpha\n",
    "            #alpha3.value = b1.alpha.val\n",
    "        else:\n",
    "            unlink(b3.alpha)\n",
    "            if alpha3_f: freeze(b3.alpha)\n",
    "            else:\n",
    "                thaw(b3.alpha)\n",
    "                set_par(b3.alpha, val=alpha3)\n",
    "            \n",
    "        if ellip3_l: \n",
    "            b3.ellip = b1.ellip\n",
    "            #ellip3.value = b1.ellip.val\n",
    "        else:\n",
    "            unlink(b3.ellip)\n",
    "            if ellip3_f: freeze(b3.ellip)\n",
    "            else: \n",
    "                thaw(b3.ellip)\n",
    "                set_par(b3.ellip, val=ellip3)\n",
    "\n",
    "        if theta3_l: \n",
    "            b3.theta = b1.theta\n",
    "            #theta3.value = b1.theta.val\n",
    "        else:\n",
    "            unlink(b3.theta)\n",
    "            if theta3_f: freeze(b3.theta)\n",
    "            else: \n",
    "                thaw(b3.theta)\n",
    "                set_par(b3.theta, val=theta3)\n",
    "            \n",
    "        if x03_l: \n",
    "            b3.xpos = b1.xpos\n",
    "            #x03.value = b1.xpos.val\n",
    "        else:\n",
    "            unlink(b3.xpos)\n",
    "            if x03_f: freeze(b3.xpos)\n",
    "            else: \n",
    "                thaw(b3.xpos)\n",
    "                set_par(b3.xpos, val=x03)\n",
    "\n",
    "        if y03_l: \n",
    "            b3.ypos = b1.ypos\n",
    "            #y03.value = b1.ypos.val\n",
    "        else:\n",
    "            unlink(b3.ypos)\n",
    "            if y03_f: freeze(b3.ypos)\n",
    "            else: \n",
    "                thaw(b3.ypos)\n",
    "                set_par(b3.ypos, val=y03)\n",
    "    except: pass\n",
    "\n",
    "    # GAUSS\n",
    "    try:\n",
    "        if fwhm_f: freeze(g1.fwhm)\n",
    "        else:\n",
    "            thaw(g1.fwhm)\n",
    "            set_par(g1.fwhm, val=fwhm)\n",
    "\n",
    "        if ampl_g_f: freeze(g1.ampl)\n",
    "        else:\n",
    "            thaw(g1.ampl)\n",
    "            set_par(g1.ampl, val=ampl_g)\n",
    "\n",
    "        if x0g_l: \n",
    "            g1.xpos = b1.xpos\n",
    "            #x0g.value = b1.xpos.val\n",
    "        else:\n",
    "            unlink(g1.xpos)\n",
    "            if x0g_f: freeze(g1.xpos)\n",
    "            else: \n",
    "                thaw(g1.xpos)\n",
    "                set_par(g1.xpos, val=x0g)\n",
    "\n",
    "        if y0g_l: \n",
    "            g1.ypos = b1.ypos\n",
    "            #y0g.value = b1.ypos.val\n",
    "        else:\n",
    "            unlink(g1.ypos)\n",
    "            if y0g_f: freeze(g1.ypos)\n",
    "            else:\n",
    "                thaw(g1.ypos)\n",
    "                set_par(g1.ypos, val=y0g)\n",
    "    except: pass\n",
    "\n",
    "    # BACKGROUND\n",
    "    if c0_f: freeze(bkg.c0)\n",
    "    else: \n",
    "        thaw(bkg.c0)\n",
    "        set_par(bkg.c0, val=c0)\n",
    "\n",
    "    make_plot()\n",
    "\n",
    "\n",
    "out2 = interactive_output(main, {\"select_galaxy\" : select_galaxy, \n",
    "                                 \"select_model\" : select_model, \n",
    "                                 \"select_scale\" : select_scale})\n",
    "\n",
    "tvl = HBox([v_blank(0),\n",
    "            VBox([VBox([HBox([select_galaxy, output_scale]),\n",
    "            HBox([select_model, select_scale])]), out2])])\n",
    "\n",
    "# tvl2 = widgets.Text(description=\"nichts\")\n",
    "\n",
    "# tab = widgets.Tab([tvl, tvl2])\n",
    "# tab.set_title(0, 'Beta modelling')\n",
    "# tab.set_title(1, 'Cavity detection')\n",
    "\n",
    "# display(tab)\n",
    "display(tvl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6423b3a7-4f44-4b63-a76f-875589e4e404",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get_model().pars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "18241280-f4a5-4d9b-b3f6-058fbbe7215b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fec8cbe-fd4c-4f22-8fdd-82ac4e434a55",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "818ee2f2-a150-4bbd-bb75-90bc319c3303",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "991d5774-86c4-42b4-8031-573666727fbb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45d7a6b3-872d-4792-9136-4befb7c19934",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81f209f8-e52b-4f8a-9b39-4e39d83fc8e1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1b984aa-46bc-4723-9789-3deb367b76ef",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}