{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 6. Exploring model output: orbits and weights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook we will look at the output of a Schwarzschild model in more detail. We will access the orbit library and the orbital weights for the model that we ran previously.\n", "\n", "As a side note, in this tutorial we will also see the that our current description of the line-of-sight velocity distribution - i.e. the Gauss Hermite expansion - has some undesirable features. This is the motivation for one of our upcoming code developments: using better descriptions of the LOSVD.\n", "\n", "Let's re-create the model that we previously ran" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import dynamite as dyn\n", "\n", "fname = 'NGC6278_config.yaml'\n", "c = dyn.config_reader.Configuration(fname)\n", "\n", "parset = c.parspace.get_parset()\n", "\n", "model = dyn.model.Model(config=c, parset=parset)\n", "orblib0 = model.get_orblib()\n", "weight_solver0 = model.get_weights(orblib0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To read the orbit library for this model, we can do the following," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "orblib0.read_losvd_histograms()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This creates the object``Orblib.LegacyOrbitLibrary.losvd_histograms`` with the following properties" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f'Orbits are stored in a {type(orblib0.losvd_histograms)}')\n", "print(f'That contains this object: {orblib0.losvd_histograms[0]}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The histogram is stored in a list, so to access it, we need to call the first list element. This object has two attributes, ``x``, and ``y``, which are the velocity array, and and array of LOSVDs of the orbit library respectively. These arrays have the following shapes," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f'velocity array has shape {orblib0.losvd_histograms[0].x.shape}')\n", "print(f'LOSVD has shape {orblib0.losvd_histograms[0].y.shape}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Where do these numbers come from? They are set by values in the configuration file." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The size of the first dimension of the LOSVD - 360 - is the number of orbits in our library. This comes from three values which we also specified by values in the configuration file, i.e." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tmp = c.settings.orblib_settings\n", "print('nE = ', tmp['nE'])\n", "print('nI2 = ', tmp['nI2'])\n", "print('nI3 = ', tmp['nI3'])\n", "print('ndithering = ', tmp['dithering'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "E, I2 and I3 are orbital integrals of motion. A grid over these values is used to specify initial conditions for the orbit library. Dithering is a the number of additional orbits which are calculated to regularise the solution. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n_orbit_bundles = tmp['nE'] * tmp['nI2'] * tmp['nI3']\n", "n_orbits_per_bundle = tmp['dithering']**3\n", "size_of_box_orblib = n_orbits_per_bundle * n_orbit_bundles\n", "# tube orbits can be reflected, so we have twice as many tube orbits as box orbits\n", "size_of_tube_orblib = 2 * n_orbit_bundles * n_orbits_per_bundle \n", "size_of_total_orblib = size_of_box_orblib + size_of_tube_orblib\n", "print('Total orbit library size =', size_of_total_orblib)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The size of the second dimension of the LOSVD corresponds to the size of the velocity array.\n", "\n", "The third dimension size of the LOSVD - 152 - is the number of spatial apertures. It should be equal to the number of spatial apertures for which we have provided kinematic data, i.e. the number of rows of the kinematic data file, i.e." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "len(c.system.cmp_list[2].kinematic_data[0].data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So - in summary - the LOSVD histogram has shape\n", "\n", "$$\n", "(n_\\mathrm{orbits}, n_\\mathrm{velocityBins}, n_\\mathrm{apertures})\n", "$$\n", "\n", "i.e. we have one LOSVD per orbit and per spatial aperture. Let's write a simple plotting routine to look at some examples," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "def plot_losvds(losvd_histogram, orb_idx, aperture_idx_list):\n", " v = losvd_histogram[0].x\n", " losvd = losvd_histogram[0].y[orb_idx, :, :]\n", " plt.plot(v, np.sum(losvd, 1), label='total')\n", " for aperture_idx in aperture_idx_list:\n", " plt.plot(v,\n", " losvd[:, aperture_idx],\n", " '--',\n", " label=f'aperture {aperture_idx}')\n", " plt.gca().set_title(f'LOSVD of orbit {orb_idx}')\n", " plt.gca().set_xlabel('v [km/s]')\n", " plt.gca().set_yscale('log')\n", " plt.gca().legend()\n", " plt.tight_layout()\n", " return" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "orb_idx = 25\n", "aperture_idx_list = [0, 2, 20, 30]\n", "plot_losvds(orblib0.losvd_histograms, orb_idx, aperture_idx_list)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some apertures, this orbit hardly contributes at all. One of these orbits looks interesting - the one in aperture 20. Let's plot this one on its own," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: plot the LOSVD of orbit 19 in apertures 0, 7, 9\n", "\n", "**Exercise**: which orbit has the narrowest LOSVD in aperture 0? Which orbit has the broadest LOSVD?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the LOSVDs that we plotted above (orbit 25 in aperture 20) looks interesting. Let's focus in on it," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "orb_idx, aperture_idx = 25, 20\n", "v_arr = orblib0.losvd_histograms[0].x\n", "losvd = orblib0.losvd_histograms[0].y[orb_idx, :, aperture_idx]\n", "\n", "plt.plot(v_arr, losvd)\n", "plt.gca().set_xlabel('v [km/s]')\n", "title = f'LOSVD of orbit {orb_idx} in aperture {aperture_idx}'\n", "_ = plt.gca().set_title(title)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's bimodal!\n", "\n", "**Exercise**: can you find any other bimodal LOSVDs in the orbit library?\n", "\n", "In order to compare the orbits we have calculated to the observations, we must transform the kinematics of the orbit LOSVDs into their Gauss Hermite representations. Let's see how well Gauss Hermite expansion can do at reproducing the bimodal LOSVD above. The function to do the transformation from orbits to observed kinematics is held in the kinematcs object itself, so let's extract this," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kinematics = c.system.cmp_list[2].kinematic_data[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and now extract the observed $v$ and sigma in the aperture of interest," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extracting observed v, sigma in a given aperture\n", "\n", "# get the row of the table kinematics.dat which corresponds to the aperture of interest\n", "row_idx = np.where(kinematics.data['vbin_id']==aperture_idx)[0][0]\n", "# extract v and sigma from that for of the table\n", "v, sigma = kinematics.data[row_idx]['v', 'sigma']\n", "# print result\n", "print(f'In aperture {aperture_idx}:')\n", "print(f' v = {v} km/s')\n", "print(f' sigma = {sigma} km/s')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's fit a Gauss Hermite distribution to the orbit's LOSVD. When we solved for the orbital weights earlier, this was all done \"under-the-hood\". Now we can see it in action," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# create a histogram object of the losvds for all apertures and orbits first\n", "losvd = orblib0.losvd_histograms[0].y[:, :, :]\n", "\n", "## velocity histograms where vel_hist.y has shape (n_orbits, n_vbins, n_regions)\n", "velhist = dyn.kinematics.Histogram(\n", " xedg=orblib0.losvd_histograms[0].xedg,\n", " y=losvd)\n", "\n", "v = kinematics.data[:]['v']\n", "sigma = kinematics.data[:]['sigma']\n", "\n", "# get a gauss hermite expansion to this LOSVD\n", "gh_coefficients = kinematics.get_gh_expansion_coefficients(\n", " v_mu=v,\n", " v_sig=sigma,\n", " vel_hist=velhist,\n", " max_order=4)\n", "\n", "# remove unused empty dimension\n", "tmp = np.squeeze(gh_coefficients)\n", "\n", "# print best-fit gauss hermite coefficients\n", "print(f'GH coefficients for orbit {orb_idx} and aperture {aperture_idx}:')\n", "i = 1\n", "for hi in tmp[orb_idx,aperture_idx,:]:\n", " print(f' h_{i} = {hi:.4f}')\n", " i+=1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's look at the LOSVD given by these Gauss Hermite coefficients," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# evaluate the losvd\n", "gh_losvd = kinematics.evaluate_losvd(velhist.x,\n", " np.array(v),\n", " np.array(sigma),\n", " gh_coefficients)\n", "# remove unused empty dimension\n", "gh_losvd = np.squeeze(gh_losvd)\n", "\n", "# normalise the two LOSVDs so that we can plot them on the same axis\n", "losvd /= np.sum(losvd) \n", "gh_losvd /= np.sum(gh_losvd)\n", "\n", "# plot the LOSVDs\n", "plt.plot(v_arr, losvd[orb_idx,:,aperture_idx])\n", "plt.plot(v_arr, gh_losvd[orb_idx,aperture_idx,:], label='GH(4) expansion')\n", "plt.gca().legend()\n", "plt.gca().set_xlabel('v [km/s]')\n", "title = f'LOSVD of orbit {orb_idx} in aperture {aperture_idx}'\n", "_ = plt.gca().set_title(title)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Gauss Hermite parameterisation of the LOSVD does an OK job at reproducing the bimodality of the LOSVD, although the peaks are somewhat underpredicted. The wings of the profile are significantly negative however. This isn't good! Negative LOSVDs are unphysical. The fact that Gauss Hermite expansions become negative could bias our modelling. This is one of our key motivations for wanting to replace Gauss Hermites with different parametrisations of the LOSVD." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: fit higher order GH expansions to the above LOSVD. How high do you need to go to improve the match and eliminate negative wings?\n", "\n", "**Exercise**: fit a GH expansion to the LOSVD of orbit 14 in aperture 19. Compare this to the histogrammed LOSVD in a plot." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Looking at the orbital weights\n", "\n", "Now let's look at the best fitting orbital weights," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model.weights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is a row for every orbit, a weight for every orbit." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What does the LOSVD look like for this set of weights? To produce this, we need to weigh the orbit library LOSVDs by the orbital weights." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract the orbital weights\n", "w = model.weights\n", "\n", "# extract the model LOSVD\n", "orbit_losvds = orblib0.losvd_histograms[0].y\n", "\n", "# dot product of model losvd and weights\n", "model_losvd = np.dot(orbit_losvds.T, w).T\n", "\n", "# normalise the model losvds so that later we can plot them on the same y-axis as the previous plots\n", "model_losvd /= np.sum(model_losvd, 0)\n", "\n", "print('model_losvd has shape = ', model_losvd.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This array is the model LOSVD. The first dimension corresponds to the velocity array, the second is the number of spatial apertures.\n", "\n", "**Exercise**: plot the model LOSVD in apertures 0, 20 and 50." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's compare the model LOSVD to to observed Gauss Hermite LOSVD in aperture 20," ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "aperture_idx = 20\n", "\n", "# extract v, sigma and GH corefficients of the observed LOSVD in out aperture\n", "v, sigma, h3, h4 = kinematics.data[aperture_idx]['v', 'sigma', 'h3', 'h4']\n", "# set h0, h1 and h2 to their defaults\n", "h0, h1, h2 = 1., 0., 0.\n", "# combine all gh coefficients to one array\n", "gh_coefficients = np.array([[[h0, h1, h2, h3, h4]]])\n", "\n", "# evaluate the LOSVD of this \n", "gh_losvd_obs = kinematics.evaluate_losvd(velhist.x,\n", " np.array([v]),\n", " np.array([sigma]),\n", " gh_coefficients)\n", "# remove un-used dimensions\n", "gh_losvd_obs = np.squeeze(gh_losvd_obs)\n", "# normalise\n", "gh_losvd_obs /= np.sum(gh_losvd_obs)\n", "\n", "# plot the \n", "plt.plot(v_arr, np.squeeze(gh_losvd_obs), label='Observed LOSVD')\n", "plt.plot(v_arr, model_losvd[:, aperture_idx], label='Model')\n", "plt.gca().legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model and observed LOSVD have similar widths. But the sharp, spiky peak in the model LOSVDs is not seen in the obseved LOSVD.\n", "\n", "The model LOSVD we just plotted is generated directly from the histogrammed LOSVDs of the orbits. We could instead look at the model LOSVD coming from the Gauss Hermite representation of the orbit library." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "aperture_idx = 20\n", "\n", "# extract the losvd of all orbits.\n", "tmp = orblib0.losvd_histograms[0].y[:,:,]\n", "\n", "# store them in a histogram object\n", "vel_hist = dyn.kinematics.Histogram(\n", " xedg=orblib0.losvd_histograms[0].xedg,\n", " y=tmp,\n", " normalise=False)\n", "\n", "\n", "# get the GH expansion coefficients for all of these orbits\n", "gh_coefficients = kinematics.get_gh_expansion_coefficients(\n", " v_mu=kinematics.data['v'][aperture_idx],\n", " v_sig=kinematics.data['sigma'][aperture_idx],\n", " vel_hist=vel_hist,\n", " max_order=4)\n", "\n", "\n", "# evaluate the losvd given these GH coefficients\n", "gh_model_losvd = kinematics.evaluate_losvd(\n", " velhist.x,\n", " kinematics.data['v'][[aperture_idx]],\n", " kinematics.data['sigma'][[aperture_idx]],\n", " gh_coefficients)\n", "\n", "\n", "# remove un-used dimensionsy