{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 2. Quickstart: Running a single Schwarzschild model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By the end of the notebook, you will have run a Schwarzschild model. This will involve:\n", "1. understanding the configuration file\n", "2. executing commands to create and run a Schwarzschild model\n", "3. plotting some output for this model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get started, let's import DYNAMITE and print the version and installation path:" ] }, { "cell_type": "code", "execution_count": null, "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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Input files" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should be in the directory ```docs/tutorial_notebooks```. To prepare the input data files, you should first run the CALIFA section of the previous tutorial \"Data Preparation for Gauss Hermite kinematics\" (``1_data_prep_for_gauss_hermite.ipynb``). \n", "\n", "If this is the first time you run this tutorial, this directory should have the following structure:\n", "\n", "```\n", "tutorial_notebooks\n", " |- ...\n", " |- NGC6278_input\n", " |- dynamite_input\n", " |- aperture.dat\n", " |- bins.dat\n", " |- gauss_hermite_kins.ecvs\n", " |- kinmaps.pdf\n", " |- mge.ecsv\n", " |- pafit.pdf\n", " |- NGC6278.V1200.rscube_INDOUSv2_SN20_stellar_kin.fits\n", " |- NGC6278_config.yaml\n", " |- NGC6278_config_single.yaml\n", " |- *.ipynb\n", " |- ...\n", "```\n", "\n", "The input data files needed to run DYNAMITE are:\n", "\n", "- ``mge.ecvs`` - the Multi Gaussian Expansion (MGE) describing the stellar surface density of the galaxy\n", "- ``gauss_hermite_kins.ecsv`` - the kinematics extracted for this galaxy\n", "- ``aperture.dat`` and ``bins.dat`` - information about the spatial apertures/bins used for kinematic extraction\n", "\n", "The MGE and kinematics files must be in the form of [Astropy ECSV files](https://docs.astropy.org/en/stable/api/astropy.io.ascii.Ecsv.html). This is a convenient and flexible file format, so we use it whenever possible for storing our input and output data.\n", "\n", "We have provided example input data file for a specific galaxy (NGC 6278), and in the directory `NGC6278_input` you can also see some diagnostic plots. A very basic set of instructions to generate your own input data is as follows:\n", "\n", "- fit a MGE to a photometric image, e.g. using [mge](http://www-astro.physics.ox.ac.uk/~mxc/software/#mge)\n", "- create a Voronoi binning of your IFU datacube, e.g. using [e.g. vorbin](http://www-astro.physics.ox.ac.uk/~mxc/software/#binning)\n", "- extract kinematics from the binned datacube, e.g. using [e.g. PPXF](http://www-astro.physics.ox.ac.uk/~mxc/software/#ppxf)\n", "- make the `aperture.dat` and `bins.dat` files; an example script to do this is ``dynamite/data_prep/generate_kin_input.py``, which is used by the previous tutorial \"Data Preparation for Gauss Hermite kinematics\" (``1_data_prep_for_gauss_hermite.ipynb``)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Configuration file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The configuration file controls all the settings that you may wish to vary when running your Schwarschild models. For example, in the configuration file you can specify:\n", "\n", "- the components of the gravitational potential\n", "- the parameters describing the potential, the parameter ranges, etc.\n", "- what type of kinematic data you are providing, e.g.\n", " - discrete vs continuous,\n", " - Gauss-Hermite vs Histograms\n", "- the location of the input and output files\n", "- the number of models you want to run\n", "\n", "This list of options is incomplete - for a more detailed description of the configuration file, see the documentation.\n", "\n", "The configuration file for this tutorial is\n", "```\n", "NGC6278_config_single.yaml\n", "```\n", "Open this file in a text editor, alongside this notebook, to see how it is organised. The file is in `yaml` format. The basic structure of a yaml file are pairs of keys and values\n", "```\n", "key : value\n", "```\n", "which can be organised into hierarchical levels separated by tabs\n", "```\n", "main_key:\n", " sub_key1 : value1\n", " sub_key2 : value2\n", "```\n", "Comments begin with a ``#``. Values can be any type of variable (integers, floats, strings, booleans, etc.)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To read the congfiguration file we can use the following command, which creates a configuration object (called `c` here):" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "fname = 'NGC6278_config_single.yaml'\n", "c = dyn.config_reader.Configuration(fname, reset_logging=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When making this object, some output is printed, telling you whether any previous models have been found. If you are running this tutorial for the first time, no models will be found and an empty table will be created at ``AllModels.table``. This table holds information about all the models which have been run so far.\n", "\n", "Remark: if you run this tutorial again and wish to start from empty output, you can either manually delete the output directory tree (`NGC6278_output`) or add the parameter `reset_existing_output=True` to the configuration object constructor: `c = dyn.config_reader.Configuration(fname, reset_logging=True, reset_existing_output=True)`.\n", "\n", "The configuration object ``c`` is structured in a similar way to the configuration file itself. For example, the configuration file is split into two sections. The top section defines aspects the physical system we wish to model - e.g. the globular cluster, galaxy or galaxy cluster - while the second section contains all other settings we need for running a model - e.g. settings about the orbit integration and input/output options. The two sections are stored in the ``system`` and ``settings`` attributes of the configuration object, respectively." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(type(c.system))\n", "print(type(c.settings))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The physical system is comprised of components, which are stored in a list ``c.system.cmp_list``." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f'cmp_list is a {type(c.system.cmp_list)}')\n", "print(f'cmp_list has length {len(c.system.cmp_list)}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's print some information about the components:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i in range(len(c.system.cmp_list)):\n", " \n", " print(f'Information about component {i}:')\n", "\n", " # Extract component i from the component list\n", " component = c.system.cmp_list[i]\n", "\n", " # Print the name of the component\n", " print(f' name = {component.name}')\n", "\n", " # Print a list of the names of the parameters of this component\n", " parameters = component.parameters\n", " parameter_names = [par0.name for par0 in parameters]\n", " string = ' has parameters : '\n", " for name in parameter_names:\n", " string += f'{name} , '\n", " print(string)\n", "\n", " # Print the type of this component\n", " print(f' type = {type(component)}')\n", "\n", " # Does it contribute to the potential?\n", " print(f' contributes to the potential? - {component.contributes_to_potential}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each component has a name, some parameters, and a type. **Currently, DYNAMITE only supports this exact combination of component types**:\n", "\n", "- one ``Plummer`` component representing the black hole\n", "- one component representing the dark halo (DYNAMITE supports ``NFW`` and a number of additional dark halos, see API documentation)\n", "- one ``TriaxialVisibleComponent`` representing the stellar body of the galaxy\n", "\n", "In the future, we plan to support other components, and more flexible combinations of components.\n", "\n", "For the stars - i.e. component 2 above - we must provide some input data files. The location of these files is specified in the configuration file, at\n", "```\n", "settings -> io_settings -> input_directory\n", "```\n", "and takes the value:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c.settings.io_settings['input_directory']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The names of the necessary input files are also specified in the configuration file, at the following locations:\n", "```\n", "system_components -> stars -> mge_pot\n", "system_components -> stars -> mge_lum\n", "system_components -> stars -> kinematics -> --> datafile\n", "system_components -> stars -> kinematics -> --> aperturefile\n", "system_components -> stars -> kinematics -> --> binfile\n", "```\n", "which take as values the relevant filenames. You are free to give these files whatever name you like, as long as it is specified in the configuration file.\n", "\n", "Let's have a look at the MGE:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c.system.cmp_list[2].mge_pot, c.system.cmp_list[2].mge_lum" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The kinematic data are stored here:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "type(c.system.cmp_list[2].kinematic_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that this object has type ``list``. This is because a single component can have multiple different sets of kinematic data. In this example, the first (and only) entry in the list is" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "type(c.system.cmp_list[2].kinematic_data[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that this kinematics object has type ``GaussHermite``. This was also specified in the configuration file, under\n", "```\n", "system_components -> stars -> kinematics --> kinset1 --> type\n", "```\n", "\n", "In addition to ``GaussHermite``, DYNAMITE also supports ``BayesLOSVD`` kinematics (see documentation).\n", "\n", "The kinemtic data itself can be accessed as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c.system.cmp_list[2].kinematic_data[0].data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating a Schwarzschild model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next step will be to create a model, i.e. a `dyn.model.Model` object for a particular `parset`. This is a particular set of values for each of the parameters of the model. In the configuration file, every parameter has a been given a `value`. We can extract and inspect a parameter set specified by these values as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "parset = c.parspace.get_parset()\n", "print(parset)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that parameters which have been specified as logarithmic, i.e. those configured as \n", "```\n", "parameters -> XXX -> logarithmic : True\n", "```\n", "are exponentiated in this table, while their value is logarithmic in the configuration file. \n", "\n", "More details can be found in the documentation and in the tutorial \"Parameter Space\" (`5_parameter_space.ipynb`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "DYNAMITE uses the `ModelIterator` object to run models. In this case, only one model will be run because the configuration file fixes all parameters to their respective values (`fixed: True`). In the following tutorial \"Model Iterations and Plots\" (`3_model_iterations_and_plots.ipynb`) we will see how to run a suite of models. To better understand what the model iterator does, let's read the internal documentation (i.e. the docstring) for the class:\n", "```\n", "class ModelIterator(object):\n", " \"\"\"Iterator for models\n", "\n", " Creating this ``ModelIterator`` object will (i) generate parameters sets,\n", " (ii) run models for those parameters, (iii) check stopping criteria, and\n", " iterate this procedure till a stopping criterion is met. This is implemented\n", " by creating a ``ModelInnerIterator`` object whose ``run_iteration`` method\n", " is called a number of times.\n", "\n", " Parameters\n", " ----------\n", " config : a ``dyn.config_reader.Configuration`` object\n", " model_kwargs : dict\n", " other kewyord argument required for this model\n", " do_dummy_run : Bool\n", " whether this is a dummy run - if so, dummy_chi2_funciton is executed\n", " instead of the model (for testing!)\n", " dummy_chi2_function : function\n", " a function of model parameters to be executed instead of the real model\n", " plots : bool\n", " whether or not to make plots\n", "\n", " \"\"\"\n", " def __init__(self,\n", " config=None,\n", " model_kwargs={},\n", " do_dummy_run=None,\n", " dummy_chi2_function=None,\n", " plots=True):\n", "```\n", "\n", "All we need to do is pass the configuration object `c` when we instantiate the `ModelIterator`. The following step will take a few minutes (we will not need the `ModelIterator` object itself, so we discard the return value):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "_ = dyn.model_iterator.ModelIterator(config=c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A lot happened here: when run for the first time, the particular `parset` is new, so DYNAMITE creates a new entry in AllModels.table and assigns a model directory that will hold the results. The next step is to calculate the orbit library and then DYNAMITE finds out which orbits are useful for reproducing the observations. This is an Non-Negative Least Squares (NNLS) optimization problem resulting in a certain chi-squared value. Let's first have a look at the all models table:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c.all_models.table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This table is stored in the `all_models_file` in the `output_directory` as defined in the configuration file:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c.settings.io_settings['output_directory'], c.settings.io_settings['all_models_file']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the output directory we find the `AllModels.table` and a `models/` subfolder (in our case, `NGC6278_output/all_models.ecsv` and `NGC6278_output/models/`, respectively), which hold the model parameters, status, and results in the so-called model directories as stated in the `AllModels.table`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look into how to access information about specific models in the `AllModels.table` (see the Model API documentation for more information). We start with instantiating a `dyn.model.Model` object from the table (here we use the only model in the table, which has index 0):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model = c.all_models.get_model_from_row(0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model.parset # Verify the model's parameter set" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model's directory name is constructed by considering the standard name for new models and their `ml` value. Here, we are at the first iteration `000` of the `ModelIterator` and at the first model `000` inside this iteration:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model.directory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Congratulations! You have run your first Schwarzschild model using DYNAMITE! \n", "\n", "The chi-squared of this model is:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "orblib = model.get_orblib()\n", "model.get_weights(orblib) # Read weights if already calculated. If not, calculate weights.\n", "model.chi2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Is that good? At this point, we don't know! To find out, we will have to run more models and compare their chi-squared values. For information on how to do this, see the following tutorial \"Model Iterations and Plots\" (`3_model_iterations_and_plots.ipynb`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot the Models\n", "\n", "Now let's look at some output for the model that we have just run. You might have noticed that the `ModelIterator` also created a directory for plots:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c.settings.io_settings['plot_directory']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This directory contains some diagnostic plots. In a run with multiple models, these plots will be updated as the calculations progress and reflect the latest state in model parameters, the chi-squared values, and the best-fit kinematic maps." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we want to plot kinematic maps for our specific model. After instantiating a `dyn.Plotter` object, we call the appropriate plotting method on our model. Please refer to the Plotting API documentation for more information." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotter = dyn.plotter.Plotter(config=c)\n", "figure = plotter.plot_kinematic_maps(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The top row shows the data, the middle row shows the model, and the bottom row shows the residuals. The columns, from left to right, are the stellar surface density, the mean velocity $V$, the velocity dispersion $\\sigma$, and the $h_3$ and $h_4$ moments. In the figure, we can see that:\n", "\n", "- the model and data surface densities are very similar\n", "- the sense of rotation of the $V$ map is reproduced well, even though the amplitude is lower than observed\n", "- the $\\sigma$, $h_3$ and $h_4$ maps are less well reproduced\n", "\n", "While the fit is certainly not perfect, it is reassuring to see that some features are already reproduced well. To improve the fit, we will have to explore the parameter space more fully. See the following tutorial \"Model Iterations and Plots\" (`3_model_iterations_and_plots.ipynb`) for more details." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise\n", "\n", "Manually adjust a model parameter, run this new model adding it to the all models table, inspect the new table, and plot the new model's kinematic maps:\n", "\n", "1. Adjust a model parameter, e.g. change the black hole mass from 5.0 to 5.2 (in the configuration file: `system_components -> bh -> parameters -> m -> value`).\n", "2. Read the new configuration file and look at the existing all models table." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c = dyn.config_reader.Configuration(fname, reset_logging=True, reset_existing_output=False)\n", "c.all_models.table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Inspect the new parameter values given in the configuration file." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c.parspace.get_parset()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4. Run the new model defined by this new parameter set and inspect the new all models table. Did the chi-squared improve?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "_ = dyn.model_iterator.ModelIterator(config=c)\n", "c.all_models.table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "5. Plot the new model's kinematic maps (the new model is in the last row of the table)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model = c.all_models.get_model_from_row(-1)\n", "plotter = dyn.plotter.Plotter(config=c)\n", "figure = plotter.plot_kinematic_maps(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Experiment some more with adding new models to the all models table..." ] } ], "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": 4 }