{ "cells": [ { "cell_type": "markdown", "id": "unexpected-patrol", "metadata": {}, "source": [ "# 3. Model Iterations and Plots\n", "\n", "By the end of the notebook, you will have run a grid of Schwarzschild models. This will involve the following steps:\n", "1. run a small grid of orbit-based models with DYNAMITE\n", "2. understand the structure of the output\n", "3. plot the output\n", "4. look at the iteration process while running the grid of models\n", "\n", "You should run this from the directory ``docs/tutorial_notebooks``.\n", "\n", "## Setup\n", "\n", "We will run DYNAMITE on CALIFA data of NGC 6278. To prepare the input data files, you should first run the tutorial \"Data Preparation for Gauss Hermite kinematics\" (``1_data_prep_for_gauss_hermite.ipynb``). The relevant files you need for this tutorial are:\n", "\n", "```\n", "| tutorial_notebooks\n", "| ├── NGC6278_input \n", "| │ ├── dynamite_input \n", "| │ │ ├── gauss_hermite_kins.ecsv\n", "| │ │ ├── aperture.dat\n", "| │ │ ├── bins.dat \n", "| │ │ ├── mge.ecsv\n", "| │ │ └── ...\n", "| │ └── ...\n", "| │ └── ...\n", "| └── NGC6278_config.yaml\n", "| └── *.ipynb\n", "|\n", "```\n", "\n", "## Read the configuration file " ] }, { "cell_type": "code", "execution_count": null, "id": "front-macintosh", "metadata": {}, "outputs": [], "source": [ "import dynamite as dyn\n", "\n", "print('DYNAMITE')\n", "print(' version', dyn.__version__)\n", "#print(' installed at ', dyn.__path__) # Uncomment to print the complete DYNAMITE installation path\n", "\n", "fname = 'NGC6278_config.yaml'\n", "c = dyn.config_reader.Configuration(fname, reset_logging=True)" ] }, { "cell_type": "markdown", "id": "packed-television", "metadata": {}, "source": [ "All the options in the configuration file are held in the object `c`. For example, let's look at the `io_settings`. Output from this tutorial will be saved in the `output_directory`." ] }, { "cell_type": "code", "execution_count": null, "id": "aac2bc0e", "metadata": {}, "outputs": [], "source": [ "# delete previous output if available\n", "c.remove_existing_orblibs()\n", "c.remove_existing_all_models_file(wipe_other_files=False)\n", "c.backup_config_file(keep=3, delete_other=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "aggregate-freight", "metadata": {}, "outputs": [], "source": [ "c.settings.io_settings" ] }, { "cell_type": "markdown", "id": "induced-occasion", "metadata": {}, "source": [ "In fact, by creating the configuration object `c`, we have also created the `output_directory` and copied a version of the configuration file there," ] }, { "cell_type": "code", "execution_count": null, "id": "complimentary-framing", "metadata": {}, "outputs": [], "source": [ "ls NGC6278_output" ] }, { "cell_type": "markdown", "id": "spread-prediction", "metadata": {}, "source": [ "## Run the models\n", "\n", "Making the `ModelIterator` object will start running a grid of orbit-based models. This next step will take about 5-10 minutes using 4 cpus " ] }, { "cell_type": "code", "execution_count": null, "id": "canadian-isolation", "metadata": { "scrolled": true }, "outputs": [], "source": [ "import time\n", "\n", "t = time.perf_counter()\n", "\n", "smi = dyn.model_iterator.ModelIterator(config=c)\n", "\n", "delta_t = time.perf_counter()-t\n", "print(f'Computation time: {delta_t} seconds = {delta_t/60} minutes')" ] }, { "cell_type": "markdown", "id": "forced-policy", "metadata": {}, "source": [ "The following files have been created in the models directory," ] }, { "cell_type": "code", "execution_count": null, "id": "innovative-introduction", "metadata": {}, "outputs": [], "source": [ "ls NGC6278_output/models" ] }, { "cell_type": "markdown", "id": "thick-holly", "metadata": {}, "source": [ "Each directory holds a different orbit library \n", "\n", " orblib_XXX_YYY\n", "\n", "where `XXX` labels the iteration when it was created, and `YYY` labels the position within that iteration. Looking inside one of these directories, we see the following files:" ] }, { "cell_type": "code", "execution_count": null, "id": "egyptian-serbia", "metadata": {}, "outputs": [], "source": [ "ls NGC6278_output/models/orblib_000_000" ] }, { "cell_type": "markdown", "id": "partial-group", "metadata": {}, "source": [ "which are:\n", "\n", "- `cmd_*`: bash scripts for running Fortran programs\n", "- `datfil/`: directory holding the orbit library for the reference potential\n", "- `infil/`: input files for running Fortran programs\n", "- `ml*/`: directories containing output orbital weights (and other results) for different values of `ml`\n", "\n", "Each `ml*` directory hold outputs for a re-scaled version of the same potential, where the value of `ml` is a mass scaling applied to a reference potential. The reference potential uses the the first value of `ml` encountered in the parameter search.\n", "\n", "Some plots are automatically created," ] }, { "cell_type": "code", "execution_count": null, "id": "intimate-oakland", "metadata": {}, "outputs": [], "source": [ "ls NGC6278_output/plots" ] }, { "cell_type": "markdown", "id": "f4204047", "metadata": {}, "source": [ "and they represent the following quantities:\n", "\n", "1. `kinchi2_progress_plot.png`: chi2 values vs model ID,\n", "2. `kinchi2_plot.png` : model parameters vs chi2 values,\n", "3. `kinematic_map_califa.png`: the kinematic maps for the current minimum-chi2 model.\n", "\n", "As an alternative to opening the plot files, we can use the model iterator's method `get_plots()`. It returns these plots as a tuple of objects," ] }, { "cell_type": "code", "execution_count": null, "id": "09e66bbb", "metadata": {}, "outputs": [], "source": [ "plots = smi.get_plots()\n", "len(plots)" ] }, { "cell_type": "code", "execution_count": null, "id": "55b7e938", "metadata": {}, "outputs": [], "source": [ "# The chi2 vs. model ID plot (kinchi2_progress_plot)\n", "print(type(plots[0]))\n", "plots[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "bf513b59", "metadata": {}, "outputs": [], "source": [ "# The model parameters vs chi2 values plot (kinchi2_plot)\n", "# If more than 2 parameters were left free, this would be a traingle plot of chi2 values.\n", "print(type(plots[1]))\n", "plots[1]" ] }, { "cell_type": "markdown", "id": "d56df462", "metadata": {}, "source": [ "The last element of the plots tuple is a list of kinematic maps for the current minimum-chi2 model with one entry for each kinematics dataset. Each list member is itself a tuple of type (`matplotlib.figure.Figure`, `str`) with the string indicating the name of the kinematics." ] }, { "cell_type": "code", "execution_count": null, "id": "799f0580", "metadata": {}, "outputs": [], "source": [ "# The kinematic maps (kinematic_map_califa)\n", "print(type(plots[2]))\n", "print(plots[2]) # If we had more than one kinematics, the list would be longer\n", "print(plots[2][0]) # The figure and the associated kinematics dataset\n", "print('First (and only) kinematics: ' + plots[2][0][1]) # The name of the kinematics\n", "plots[2][0][0] # The plot" ] }, { "cell_type": "markdown", "id": "collaborative-america", "metadata": {}, "source": [ "A summary of all the models run so far is saved in the file `NGC6278_output/all_models.ecsv`. This is an Astropy ECSV file. A table holding this data is stored in `c`, " ] }, { "cell_type": "code", "execution_count": null, "id": "tender-google", "metadata": {}, "outputs": [], "source": [ "c.all_models.table" ] }, { "cell_type": "markdown", "id": "pursuant-battlefield", "metadata": {}, "source": [ "At this stage, you could:\n", " \n", "- run more models, perhaps first adjusting settings in the configuration file, as explained in `2_quickstart.ipynb`,\n", " - increasing the `n_max_mods` and/or `n_max_iter`\n", " - adjust parameter bounds and/or which parameters are kept free\n", "- plot other visualisations\n", "\n", "## Plotting \n", "\n", "DYNAMITE provides other plotting methods in the `Plotter`:" ] }, { "cell_type": "code", "execution_count": null, "id": "preliminary-composition", "metadata": {}, "outputs": [], "source": [ "plotter = dyn.plotter.Plotter(config=c)" ] }, { "cell_type": "markdown", "id": "guided-squad", "metadata": {}, "source": [ "In all the functions that require to use the values of the $\\chi^2$ to make the plots, the user can choose which $\\chi^2$ to use, by specifying the value of the parameter ``which_chi2``. The recommended value to use is ``which_chi2='kinchi2'``. \n", "\n", "The plots produced by the functions introduced below are all saved in the plot directory specified in the directory ``plots`` within the output directory specified in the configuration file. After running all the cells in this notebook, you will find your plots in the directory ``NGC6278_output/plots``.\n", "\n", "The ``mass_plot`` function generates a cumulative mass plot, showing the enclosed mass profiles for the mass-follows-light component (red), for the dark matter (blue), and for the sum of the two (black). The solid lines correspond to the best-fit model, the shaded areas represent 1 sigma uncertainties. You can specify the radial extent of the plot and the type of file you want to be saved with the figure (e.g., ``'.png'``, ``'.pdf'``, ... if ``figtype=None``, the default is used and a ``'.png'`` figure is created)." ] }, { "cell_type": "code", "execution_count": null, "id": "photographic-debut", "metadata": {}, "outputs": [], "source": [ "fig1 = plotter.mass_plot(which_chi2='kinchi2', Rmax_arcs=50, figtype=None)" ] }, { "cell_type": "markdown", "id": "protective-abraham", "metadata": {}, "source": [ "The ``orbit_plot`` function generates a plot showing the stellar orbit distribution, described as probability density of orbits; circularity ($\\lambda_z$) is represented here as a function of the distance from the galactic centre r (in arcsec). You can specify the type of file you want to be saved with the figure (e.g., ``'.png'``, ``'.pdf'``, ...). In this case, ``Rmax_arcs`` represents the upper radial limit for orbit selection, in arcsec, meaning that only orbits extending up to ``Rmax_arcs`` are plotted." ] }, { "cell_type": "code", "execution_count": null, "id": "affecting-debut", "metadata": {}, "outputs": [], "source": [ "fig2 = plotter.orbit_plot(Rmax_arcs=50)" ] }, { "cell_type": "markdown", "id": "recorded-school", "metadata": {}, "source": [ "The ``beta_plot`` function generates two plots, showing the intrinsic and projected anisotropy profiles." ] }, { "cell_type": "code", "execution_count": null, "id": "statutory-assumption", "metadata": {}, "outputs": [], "source": [ "fig3, fig4 = plotter.beta_plot(which_chi2='kinchi2', Rmax_arcs=50)" ] }, { "cell_type": "markdown", "id": "painful-program", "metadata": {}, "source": [ "The ``qpu_plot`` function creates a plot showing the intrinsic flattenings $q$ and $p$, with the blue and black lines respectively, as a function of the distance from the galactic centre (in arcsec). The value of $T = (1-p^2)/(1-q^2)$ is also shown (red line)." ] }, { "cell_type": "code", "execution_count": null, "id": "internal-ridge", "metadata": {}, "outputs": [], "source": [ "fig5 = plotter.qpu_plot(which_chi2='kinchi2', Rmax_arcs=50,figtype =None)" ] }, { "cell_type": "markdown", "id": "0b8188ca", "metadata": {}, "source": [ "## The Iteration Process\n", "\n", "We will now have a look at the iteration process by examining the $\\chi^2$-value.\n", "\n", "In the table above (`c.all_models.table`), the model parameters are given for each step of the iteration. Most parameters are fixed and therefore constant, but the mass-to-light ratio `ml` and the dark matter parameter `f-dh` are free. We will examine how these parameters are changed with each iteration.\n", "\n", "\n", "First, we get the upper and lower limits of the parameters from the configuration file (for more details on this step have a look at the notebook `parameter_space.ipynb`):" ] }, { "cell_type": "code", "execution_count": null, "id": "1e451334", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# extract the lo/hi limits for the two free parameters f and ml\n", "\n", "f = c.parspace.get_parameter_from_name('f-dh')\n", "f_lims_raw = [f.par_generator_settings['lo'], f.par_generator_settings['hi']]\n", "f_lims = [f.get_par_value_from_raw_value(lim0) for lim0 in f_lims_raw]\n", "\n", "ml = c.parspace.get_parameter_from_name('ml')\n", "ml_lims_raw = [ml.par_generator_settings['lo'], ml.par_generator_settings['hi']]\n", "ml_lims = [ml.get_par_value_from_raw_value(lim0) for lim0 in ml_lims_raw]" ] }, { "cell_type": "markdown", "id": "b928876c", "metadata": {}, "source": [ "From the table, we can get a list of all iterations that were made:" ] }, { "cell_type": "code", "execution_count": null, "id": "943f0585", "metadata": {}, "outputs": [], "source": [ "# get list of iterations\n", "iterations = np.unique(c.all_models.table['which_iter'])" ] }, { "cell_type": "markdown", "id": "385abfa9", "metadata": {}, "source": [ "Now plot the $\\chi^2$-value of the models for every iteration:" ] }, { "cell_type": "code", "execution_count": null, "id": "39630de4", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# plot chi2 of models vs iterations\n", "for iter0 in iterations:\n", " table = c.all_models.table\n", " table = table[table['which_iter']==iter0]\n", " plt.figure(figsize=(6,5))\n", " plt.scatter(table['f-dh'],\n", " table['ml'],\n", " c=table['chi2'],\n", " cmap=plt.cm.viridis_r,\n", " s=200)\n", " \n", " \n", " cbar = plt.colorbar()\n", " cbar.set_label('$\\chi^2$')\n", " plt.gca().set_title(f'iteration {iter0}')\n", " \n", " plt.gca().set_xlim(*f_lims)\n", " plt.gca().set_ylim(*ml_lims)\n", " plt.gca().set_xlabel(f.LaTeX)\n", " plt.gca().set_ylabel(ml.LaTeX)\n", " plt.gca().set_xscale('log')\n", " plt.tight_layout()\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "100b652a", "metadata": {}, "source": [ "We can see how the parameter values get closer together for each iteration. The next plot gives an overview over all iterations:" ] }, { "cell_type": "code", "execution_count": null, "id": "8f7c9287", "metadata": {}, "outputs": [], "source": [ "# plot the models: f-dh vs ml altogether\n", "plt.scatter(c.all_models.table['f-dh'],\n", " c.all_models.table['ml'],\n", " c=c.all_models.table['chi2'] - np.min(c.all_models.table['chi2']),\n", " cmap=plt.cm.viridis_r,\n", " s=200)\n", "cbar = plt.colorbar()\n", "cbar.set_label('$\\chi^2$')\n", "plt.gca().set_title(f'all iterations')\n", "plt.gca().set_xlim(*f_lims)\n", "plt.gca().set_ylim(*ml_lims)\n", "plt.gca().set_xlabel(f.LaTeX)\n", "plt.gca().set_ylabel(ml.LaTeX)\n", "plt.gca().set_xscale('log')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "316f9225", "metadata": {}, "source": [ "We can see how with every iteration, the $\\chi^2$-value becomes smaller, as does the scatter between the values.\n", "\n", "The plot `kinchi2_plot` already shown above is another representation of the iteration process. The different models are plotted on a grid in the parameter space, where the best fit model is marked with a cross. Above we stored this plot in `plots[1]`," ] }, { "cell_type": "code", "execution_count": null, "id": "3d777ae4", "metadata": {}, "outputs": [], "source": [ "plots[1]" ] } ], "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.9.13" } }, "nbformat": 4, "nbformat_minor": 5 }