{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting Gaussians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# standard imports\n",
    "import matplotlib.pyplot as plt\n",
    "import bettermoments as bm\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load up the data\n",
    "path = '../../../gofish/docs/user/TWHya_CS_32.fits'\n",
    "data, velax, bunits = bm.load_cube(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply spectral smoothing\n",
    "smoothed_data = bm.smooth_data(data=data,\n",
    "                               smooth=0,\n",
    "                               polyorder=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# estimate the RMS of the cube\n",
    "rms = bm.estimate_RMS(data=data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load up a user-defined mask\n",
    "user_mask = bm.get_user_mask(data=data,\n",
    "                             user_mask_path=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a threshold-based mask\n",
    "threshold_mask = bm.get_threshold_mask(data=data,\n",
    "                                       clip=15.0,\n",
    "                                       smooth_threshold_mask=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a channel-based mask.\n",
    "channel_mask = bm.get_channel_mask(data=data,\n",
    "                                   firstchannel=0,\n",
    "                                   lastchannel=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine the three masks\n",
    "mask = bm.get_combined_mask(user_mask=user_mask,\n",
    "                            threshold_mask=threshold_mask,\n",
    "                            channel_mask=channel_mask,\n",
    "                            combine='and')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7ff4ed182d50>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATAAAAD8CAYAAADwpviIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAW2ElEQVR4nO3de7BdZX3G8e9zbrmRkIRADEk0odx0bBV6CmFwHCSKiIzQGXRgbEkpncy02uKlo9B2xunlD+10vM100IygsVWURiwUURoBh9qp0XA1JGCO4ZaYCyEXyP2cvX/9Y78HtiEne5+z1tl7v4fnM/POWetd79nrzUr48V7XVkRgZpajrnZXwMxsrBzAzCxbDmBmli0HMDPLlgOYmWXLAczMstXT7gqYmdWT9AzwMlABhiKif6SyDmBm1oneFRE7GxVyF9LMsqVOXonfp0kxmWntrobZhPYyu3dGxMlFPuO975oWL+6qNCz30OOHnwAO1WWtiIgV9WUkPQ3sBgL46tHX63V0F3Iy0zhfS9tdDbMJ7cex6tmin7FzV4U19y5oWK533q8PHW9MK3lHRGyRdAqwWtKTEfHgsQq6C2lmJQgqUW2YmvqkiC3p5w7g+8B5I5V1ADOzwgKoEg1TI5KmSZo+fAxcAqwbqXxHdyHNLB9VmmthNTAX+L4kqMWnb0fEj0Yq7ABmZoUFwWCTXcTjfk7EJuBtzZZ3ADOzwgKoNNFFLJsDmJmVopkxrrI5gJlZYQFU2rCm1AHMzEpRyhD+KDmAmVlhQXgMzMzyFAGDbdiV6ABmZiUQFdTyuzZciS/pVkk7JK2ry5stabWkjennrJQvSV+WNCDpcUnn1v3OslR+o6Rl4/PHMbN2CKAajVPZmtlK9A3g0qPybgTui4gzgPvSOcD7gDNSWg7cDLWAB3wGOJ/avqbPDAc9M5sYKqkVdrxUtoYBLO0C33VU9hXAynS8EriyLv+bUfMzYKakecB7gdURsSsidgOreW1QNLNM1Raytj6AjXUMbG5EbE3H26jtXwKYDzxfV25zyhsp/zUkLafWemMyU8dYPTNrpQAGo/Xvhig8iB8RIam03m16edkKgBma3blvWzSzVwSi0oaX24z1jttT15D0c0fK3wIsrCu3IOWNlG9mE0Q11DCVbawB7C5geCZxGXBnXf61aTZyCbA3dTXvBS6RNCsN3l+S8sxsAujYMTBJtwEXAXMkbaY2m/hZ4HZJ1wPPAh9Kxe8BLgMGgAPAdQARsUvSPwK/SOX+ISKOnhgws2yJSieOgUXENSNces3L6qP2DSEfGeFzbgVuHVXtzCwLtTeydmAAMzNrJEIcie6W39cBzMxKUW3DViIHMDMrrDaI7y6kmWWpQwfxzcwa8SC+mWWtMg4LVRtxADOzwgIxGK0PJw5gZlaYB/HNLFuB3IU0s3x5EN/MshSBl1GYWZ5qg/jeSmRmmfIgvpllKRifFxY24gBmZqVwC8zMslT7XkgHMDPLUnu+mdsBzMwKq32tmmchzSxDEXIX0szyVeZCVkndwFpgS0RcPlI5BzAzK6z2PrBSx8BuADYAM45XqPVtPjObgGpvZG2UmvokaQHwfuBrjcq6BWZmhdWWUTTVApsjaW3d+YqIWHFUmS8CnwKmN/owBzAzK2wUeyF3RkT/SBclXQ7siIiHJF3U6MMcwMysFCW9TudC4AOSLgMmAzMk/XtE/NGxCnsMzMwKq71ORw1T48+JmyJiQUQsAq4G7h8peIFbYGZWEm/mNrMs1d5GUW6HLiJ+AvzkeGUcwMyssNpWIq/EN7MseSuRmWWs5JX4TSkUMiV9XNITktZJuk3SZEmLJa2RNCDpu5L6UtlJ6XwgXV9Uyp/AzNqurFnI0RpzAJM0H/groD8i3gp0U5v2/BzwhYg4HdgNXJ9+5Xpgd8r/QipnZhNENboaprIV/cQeYIqkHmAqsBW4GFiVrq8ErkzHV6Rz0vWlklrf5jSz0g2/E79RKtuYA1hEbAH+BXiOWuDaCzwE7ImIoVRsMzA/Hc8Hnk+/O5TKnzTW+5tZ5whgKLoaprIV6ULOotaqWgycCkwDLi1aIUnLJa2VtHaQw0U/zsxaJLcu5LuBpyPihYgYBO6gto9pZupSAiwAtqTjLcBCgHT9RODFoz80IlZERH9E9PcyqUD1zKxlmug+dlQXklrXcYmkqWksaymwHngAuCqVWQbcmY7vSuek6/dHRBS4v5l1iOEXGjZKZRvzOrCIWCNpFfAwMAQ8AqwAfgB8R9I/pbxb0q/cAvybpAFgF7UZSzObILLbCxkRnwE+c1T2JuC8Y5Q9BHywyP3MrDON4oWGpfJKfDMrLBBDVW8lMrNMtWMrkQOYmRUX7kKaWaY8BmZmWXMAM7MsBaLiQXwzy5UH8c0sS+FBfDPLWTiAmVmexmezdiMOYGZWCrfAzCxLEVCpOoCZWaY8C2lmWQrchTSzbHkQ38wy1o73KzuAmVkp3IU0syzVZiG9F9LMMuUupJlly11IM8tSIAcwM8tXO77k1QHMzIoLiBK2EkmaDDwITKIWn1alr288JgcwMytFSV3Iw8DFEbFPUi/wU0k/jIifHauwA5iZlaKMWciICGBfOu1NacRPbv3CDTObcIb3QjZKwBxJa+vS8qM/S1K3pEeBHcDqiFgz0n3dAjOz4gJorgu5MyL6j/tRERXg7ZJmAt+X9NaIWHessm6BmVkpIhqn0X1e7AEeAC4dqYwDmJmVQES1cWr4KdLJqeWFpCnAe4AnRyrvLqSZlaOchWDzgJWSuqk1sG6PiLtHKuwAZmbFRTnLKCLiceCcZss7gJlZOdqwFL/QGJikmZJWSXpS0gZJF0iaLWm1pI3p56xUVpK+LGlA0uOSzi3nj2BmnUFNpHIVHcT/EvCjiDgbeBuwAbgRuC8izgDuS+cA7wPOSGk5cHPBe5tZJ6k2kUo25gAm6UTgncAtABFxJE17XgGsTMVWAlem4yuAb0bNz4CZkuaN9f5m1kGG14E1SiUr0gJbDLwAfF3SI5K+JmkaMDcitqYy24C56Xg+8Hzd729Oeb9F0vLhVbqDHC5QPTNrpbLXgTWjSADrAc4Fbo6Ic4D9vNpdBF7Z1zSqakfEiojoj4j+XiYVqJ6ZtVQ0kUpWJIBtBjbX7VNaRS2gbR/uGqafO9L1LcDCut9fkPLMbCLIqQsZEduA5yWdlbKWAuuBu4BlKW8ZcGc6vgu4Ns1GLgH21nU1LXcSPW9aSPeZv4N6vDrn9UjROJWt6L+0vwS+JakP2ARcR1o9K+l64FngQ6nsPcBlwABwIJW13HV10336IrZdfAqzrtrClJ4DbP/2HzD39vVU9uxtd+2sVUJQwgsNR6tQAIuIR4Fj7SxfeoyyAXykyP2sw0j0vHE+u/pP5sDF+7hu3mNsPXIiT100l+4jb+Hke59maOu2dtfSWiW3haz2+tU1bRpHLvl9fnP5AkIweKSH/9l9OgDvPG2AwT/czbPXnkbX296MevvaXFtricwG8e31SuLIkrN54Zw+QjBt6yCLvt7Fk3ecxeFqD3/9hv/m7958D30XvshvLppF1+KFHhd7PWhDAPO/KhuT7oNDHJ7VR/dBMemRTVR27+bUB/u4/8gSfvejmzl/8jMcHuxBU+DAGbPpm3MC3Y8NUN2/v91Vt/HQ/AsNS+UAZqPSPfNEmP8G4pebWByn0bt9L0O7dwMQg0eY+5Wf853/fQ9fPXs6k2aK/fOD3Wf2Mu9/DlI9eKjNtbfxNB6zjI04gNmoVPa+RNfgENX9+9H/PcbQUddjaAgeeYIZj0D33FN4Zvnp9L4UaP0molppS52tRTyIb51OfX10nTKnqbJx4CCznqrSty84sPStdM85aZxrZ+3UjnVgDmA2KnH4MNWt2+meeWLD2cXqvn30vlzh5Td28fKCbjT9hBbV0tqiDSvx3YW0UemaPJl9l72N6IIZj71AZeBp1NdHHDny2t26EUzduBOYw9Rn9jD07Oa21NlaYJxmGRtxALNRqR46xIw1r75UZMefX8Ce/iNM3djHm24doLJ9x2+V174DHJ7ZzbSDh8FjYBObA5jlYGjLb1BPD9v+4jy0dBfvP/VpBt/axYO95/DGH81Fg1W04dfohGnEoUPM+sV2Yv/BdlfbxpnG4YWFjTiA2Zh0TZ9OZRLs3zuVNd1voq9niKGzD7DxlKlMe7abUyedycG5k6lM6mLGD37p9V+vB56FtByot4+D559OtQe0s4+dL0xnz/4pvHn+Nuad/gLVXqAaTN55hH2ndjU9a2n5amYGshPfRmGvQzF4hClrN3HCyWdSmSIGu3o5MNjFM91VBEzZHvQ8t4PKi7tZsOUNVDb/pt1VtlbwSnzLRWXni5z0X+vpe9fZHDi5C1V7mLxnBic8vY/u7c8xtGMnVCsMPfNcu6tqreJBfMtJZc9epv7nz5nW3U1UA6oVAl6zOt9eH7yVyPITUds+ZK9v4VlIM8uZW2Bmli0HMDPLVTvGwLwOzMyy5RaYmZXDXUgzy5JnIc0sa26BmVmOhBeymlnOPAtpZlkq6W0UkhZKekDSeklPSLrheOXdAjOzcpQziD8EfDIiHpY0HXhI0uqIWH+swm6BmVkpymiBRcTWiHg4Hb8MbADmj1TeLTAzK0dzY2BzJK2tO18RESuOVVDSIuAcYM1IH+YAZmbFNf+tRDsjor9RIUknAN8DPhYRL41UzgHMzEpR1jIKSb3Ugte3IuKO45UtPAYmqVvSI5LuTueLJa2RNCDpu5L6Uv6kdD6Qri8qem8z6yDRRGpAkoBbgA0R8flG5csYxL+B2kDbsM8BX4iI04HdwPUp/3pgd8r/QipnZhOEqo1TEy4E/hi4WNKjKV02UuFCAUzSAuD9wNfSuYCLgVWpyErgynR8RTonXV+ayptZ7pppfTU3C/nTiFBE/F5EvD2le0YqX7QF9kXgU7y6AuQkYE9EDL9jeDOvToHOB55PlRwC9qbyv0XScklrJa0d5HDB6plZK6jJVLYxBzBJlwM7IuKhEutDRKyIiP6I6O9lUpkfbWbjqYQW2GgVmYW8EPhA6p9OBmYAXwJmSupJrawFwJZUfguwENgsqQc4EXixwP3NrINk9UbWiLgpIhZExCLgauD+iPgw8ABwVSq2DLgzHd+VzknX74+INvyRzWxctKEFNh5biT4NfELSALUxrltS/i3ASSn/E8CN43BvM2uHKG0WclRKWcgaET8BfpKONwHnHaPMIeCDZdzPzDqQ3wdmZrnyCw3NLF8OYGaWK7fAzCxPQVkvNBwVBzAzK8xf6mFmeXMAM7NcqQ3r0h3AzKy4cVpp34gDmJmVwmNgZpat8dgq1IgDmJmVwy0wM8tSk9/7WDYHMDMrhwOYmeXIC1nNLGuqeh2YmeXI68DMLGdeRmFm+XILzMxy5UF8M8tTAN7MbWa58hiYmWXJ68DMLF8R7kKaWb7cAjOzfDmAmVmu3AIzszwFUPEYmJllqh0tsK7W39LMJqThmcjjpSZIulXSDknrGpV1ADOzUigapyZ9A7i0mYIOYGZWXDSZmvmoiAeBXc2UHXMAk7RQ0gOS1kt6QtINKX+2pNWSNqafs1K+JH1Z0oCkxyWdO9Z7m1lnEaBKNEzAHElr69LyIvctMog/BHwyIh6WNB14SNJq4E+A+yLis5JuBG4EPg28DzgjpfOBm9NPM5sAmvxm7p0R0V/WPcfcAouIrRHxcDp+GdgAzAeuAFamYiuBK9PxFcA3o+ZnwExJ88Z6fzPrICV2IUejlDEwSYuAc4A1wNyI2JoubQPmpuP5wPN1v7Y55R39WcuHm5eDHC6jemY27pqYgRyHvZKFA5ikE4DvAR+LiJfqr0XEqONuRKyIiP6I6O9lUtHqmVmLlDULKek24P+AsyRtlnT9SGULLWSV1EsteH0rIu5I2dslzYuIramLuCPlbwEW1v36gpRnZhNBSS2siLim2bJFZiEF3AJsiIjP1126C1iWjpcBd9blX5tmI5cAe+u6mmaWs2h6FrJURVpgFwJ/DPxS0qMp72+AzwK3p2bfs8CH0rV7gMuAAeAAcF2Be5tZp8lpM3dE/JTa8o9jWXqM8gF8ZKz3M7PO1uQyilJ5M7eZlcMBzMyyFIC/1MPMciTCXUgzy1i19U0wBzAzK85dSDPLmbuQZpYvBzAzy5O/2NbMcuVvJTKznHkMzMzy5QBmZlkKoOoAZmZZ8iC+meXMAczMshRAxVuJzCxLAeEAZma5chfSzLLkWUgzy5pbYGaWLQcwM8tSBFQqLb+tA5iZlcMtMDPLlgOYmeUpPAtpZpkKCC9kNbNseSuRmWUpwl+rZmYZ8yC+meUq3AIzszz5hYZmlitv5jazXAUQbdhK1NXqG0q6VNJTkgYk3djq+5vZOIj0QsNGqQmjiREtDWCSuoF/Bd4HvAW4RtJbWlkHMxsfUY2GqZHRxohWt8DOAwYiYlNEHAG+A1zR4jqY2XgopwU2qhjR6jGw+cDzdeebgfPrC0haDixPp4d/HKvWtahuYzUH2NnuShyH61dMp9cPitfxTUUr8DK77/1xrJrTRNHJktbWna+IiBV15w1jRL2OG8RPf5gVAJLWRkR/m6t0XJ1eR9evmE6vH3RGHSPi0nbct9VdyC3AwrrzBSnPzAxGGSNaHcB+AZwhabGkPuBq4K4W18HMOteoYkRLu5ARMSTpo8C9QDdwa0Q8cZxfWXGca52i0+vo+hXT6fWDPOrYlNHGCEUblv+bmZWh5QtZzczK4gBmZtnq2ADWCVuOJC2U9ICk9ZKekHRDyp8tabWkjennrJQvSV9OdX5c0rktqme3pEck3Z3OF0tak+rx3TQYiqRJ6XwgXV/UovrNlLRK0pOSNki6oJOeoaSPp7/fdZJukzS5nc9Q0q2SdkhaV5c36uclaVkqv1HSsrLr2REiouMStcG7XwOnAX3AY8Bb2lCPecC56Xg68Ctq2xv+Gbgx5d8IfC4dXwb8EBCwBFjTonp+Avg2cHc6vx24Oh1/BfjzdPwXwFfS8dXAd1tUv5XAn6XjPmBmpzxDagsnnwam1D27P2nnMwTeCZwLrKvLG9XzAmYDm9LPWel4Viv+vluZ2l6BEf4CLwDurTu/CbipA+p1J/Ae4ClgXsqbBzyVjr8KXFNX/pVy41inBcB9wMXA3ekf8k6g5+hnSW1m54J03JPKaZzrd2IKEDoqvyOeIa+u/J6dnsndwHvb/QyBRUcFsFE9L+Aa4Kt1+b9VbqKkTu1CHms7wfw21QWA1FU4B1gDzI2IrenSNmBuOm5Hvb8IfAoY3mh2ErAnIoaOUYdX6peu703lx9Ni4AXg66mb+zVJ0+iQZxgRW4B/AZ4DtlJ7Jg/RWc8QRv+8Ou6/ofHQqQGso0g6Afge8LGIeKn+WtT+99aWtSiSLgd2RMRD7bh/k3qodYdujohzgP3UukCvaPMznEVts/Bi4FRgGtCWbTHNaufz6jSdGsA6ZsuRpF5qwetbEXFHyt4uaV66Pg/YkfJbXe8LgQ9Ieobarv2LgS8BMyUNL1Kur8Mr9UvXTwReHMf6Qe3//JsjYk06X0UtoHXKM3w38HREvBARg8Ad1J5rJz1DGP3z6pj/hsZTpwawjthyJEnALcCGiPh83aW7gOFZnWXUxsaG869NM0NLgL11zf7SRcRNEbEgIhZRe0b3R8SHgQeAq0ao33C9r0rlx/X/5BGxDXhe0lkpaymwng55htS6jkskTU1/38P165hneIz7NvO87gUukTQrtTIvSXkTS7sH4UZK1GZXfkVtNvJv21SHd1Brqj8OPJrSZdTGPO4DNgI/Bman8qL2MrZfA78E+ltY14t4dRbyNODnwADwH8CklD85nQ+k66e1qG5vB9am5/if1GbFOuYZAn8PPAmsA/4NmNTOZwjcRm08bpBaC/b6sTwv4E9TPQeA61r1b7GVyVuJzCxbndqFNDNryAHMzLLlAGZm2XIAM7NsOYCZWbYcwMwsWw5gZpat/wdWHGrNX/gulgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "im = ax.imshow(np.sum(mask, axis=0), origin='lower')\n",
    "plt.colorbar(im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mask the smoothed data\n",
    "masked_data = smoothed_data * mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "# identify the pixels to fit with the MCMC\n",
    "indices = bm.get_finite_pixels(masked_data, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10/10 [00:43<00:00,  4.35s/it]\n",
      "100%|██████████| 11/11 [00:47<00:00,  4.28s/it]\n",
      "100%|██████████| 10/10 [00:45<00:00,  4.57s/it]\n",
      "100%|██████████| 11/11 [00:48<00:00,  4.39s/it]\n",
      "100%|██████████| 10/10 [00:44<00:00,  4.49s/it]\n",
      "100%|██████████| 11/11 [00:47<00:00,  4.32s/it]\n",
      "100%|██████████| 10/10 [00:42<00:00,  4.25s/it]\n",
      "100%|██████████| 11/11 [00:45<00:00,  4.15s/it]\n"
     ]
    }
   ],
   "source": [
    "# collapse the data\n",
    "\n",
    "fits = bm.collapse_analytical(velax, data, rms, 'gaussian',\n",
    "                              indices=indices, chunks=8,\n",
    "                              mcmc='emcee',\n",
    "                              nwalkers=32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}