Note: This page was generated from a Jupyter notebook which can be found at docs/tutorial_notebooks/2_quickstart.ipynb in your DYNAMITE directory. —-

2. Quickstart: Running a single Schwarzschild model

By the end of the notebook, you will have run a Schwarzschild model. This will involve, 1. understanding the configuration file 2. executing commands to create and run a Schwarzschild model 3. plotting some output for this model

Setup

You should be in the directory docs/tutorial_notebooks. If this is the first time you have run this tutorial, this directory should have the following structure,

tutorial_notebooks
  |-
  | |- NGC6278_input
  |     |- mge.ecvs
  |     |- gauss_hermite_kins.ecvs
  |     |- bins.dat
  |     |- aperture.dat
  |
  |- NGC6278_config.yaml
  |- *.ipynb

We have provided example input data files and a configution file for the galaxy NGC6278. The input data files are:

• mge.ecvs - the Multi Gaussian Expansion (MGE) describing the stellar surface density

• gauss_hermite_kins.ecvs - the kinematics extracted for this galaxy

• aperture.dat - information about the spatial apertures/bins using for kinematic extraction

• bins.dat - information about the spatial apertures/bins using for kinematic extraction

The MGE and kinematics files must be in the form of Astropy ECSV files. This is a convenient and flexible file format, so we use it wherever possible for storing our input and output data.

We have provided test data for one galaxy. A very basic set of instructions for generating your own input data are as follows,

• fit and MGE to a photometric image e.g. using mge

• Voronoi bin your IFU datacube e.g. using e.g. vorbin

• extract kinematics from the binned datacube e.g. using e.g. PPXF

• make the apertures and binfile files. An example script which has created these files is the file generate_input_califa.py in this directory

In the future, we will provide more examples of test data, and more detailed instructions for preparing your own data.

To get started, let’s import DYNAMITE and print the version and installation path,

import dynamite as dyn

print('DYNAMITE')
print('    version', dyn.__version__)
print('    installed at ', dyn.__path__)

DYNAMITE
    version 2.0.0
    installed at  ['/Users/leahskusa/Desktop/uni/dynamics_assistant/dynamite/dynamite_test/DYN_TEST/lib/python3.9/site-packages/dynamite-2.0.0-py3.9.egg/dynamite']

Reading the configuration file

The configuration file controls all of the settings that you may wish to vary when running your Schwarschild models, e.g.

• specifying the components of the gravitational potential

• specifying the potential parameters, or parameter ranges, etc

• specify what type of kinematic data you are providing, e.g.

  • discrete vs continuous,

  • Gauss-Hermite vs Histograms

• the location of the input and output files

• the number of models you want to run

This list of options is incomplete - for a more detailed description of the configuration file, see the documentation.

The configuration file for this tutorial is

NGC6278_config.yaml

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 files are pairs of keys and values

key : value

which can be organised into hierarchical levels separated by tabs

main_key:
    sub_key1 : value1
    sub_key2 : value2

Comments begin with a #. Values can be any type of variable e.g. integers, floats, strings, booleans etc.

To read in the congfiguration file we can use the following command, creating a configuration object which here we call c,

fname = 'NGC6278_config.yaml'
c = dyn.config_reader.Configuration(fname, reset_logging=True)

[INFO] 12:45:28 - dynamite.config_reader.Configuration - Config file NGC6278_config.yaml read.
[INFO] 12:45:28 - dynamite.config_reader.Configuration - io_settings...
[INFO] 12:45:28 - dynamite.config_reader.Configuration - Output directory tree: NGC6278_output/.
[INFO] 12:45:28 - dynamite.config_reader.Configuration - system_attributes...
[INFO] 12:45:28 - dynamite.config_reader.Configuration - model_components...
[INFO] 12:45:29 - dynamite.config_reader.Configuration - system_parameters...
[INFO] 12:45:29 - dynamite.config_reader.Configuration - orblib_settings...
[INFO] 12:45:29 - dynamite.config_reader.Configuration - weight_solver_settings...
[INFO] 12:45:29 - dynamite.config_reader.Configuration - Will attempt to recover partially run models.
[INFO] 12:45:29 - dynamite.config_reader.Configuration - parameter_space_settings...
[INFO] 12:45:29 - dynamite.config_reader.Configuration - multiprocessing_settings...
[INFO] 12:45:29 - dynamite.config_reader.Configuration - ... using 4 CPUs for orbit integration.
[INFO] 12:45:29 - dynamite.config_reader.Configuration - ... using 4 CPUs for weight solving.
[INFO] 12:45:29 - dynamite.config_reader.Configuration - legacy_settings...
[INFO] 12:45:29 - dynamite.config_reader.Configuration - System assembled
[INFO] 12:45:29 - dynamite.config_reader.Configuration - Configuration validated
[INFO] 12:45:29 - dynamite.config_reader.Configuration - Instantiated parameter space
[INFO] 12:45:29 - dynamite.model.AllModels - Previous models have been found: Reading NGC6278_output/all_models.ecsv into AllModels.table
[INFO] 12:45:29 - dynamite.config_reader.Configuration - Instantiated AllModels object

On making this object, some output is printed telling us whether any previous models have been found. Assuming that you’re running this tutorial for the first time, then no models will be found and an empty table is created at AllModels.table. This table holds holds information about all the models which have been run so far.

The configuration object c is structured in a similar way to the 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,

print(type(c.system))
print(type(c.settings))

<class 'dynamite.physical_system.System'>
<class 'dynamite.config_reader.Settings'>

The physical system is comprised of components, which are stored in a list c.system.cmp_list

print(f'cmp_list is a {type(c.system.cmp_list)}')
print(f'cmp_list has length {len(c.system.cmp_list)}')

cmp_list is a <class 'list'>
cmp_list has length 3

Let’s print information about the components,

for i in range(3):

    print(f'Information about component {i}:')

    # extract component i from the component list
    component = c.system.cmp_list[i]

    # print the name
    print(f'   name =  {component.name}')

    # print a list of the names of the parameters of this component
    parameters = component.parameters
    parameter_names = [par0.name for par0 in parameters]
    string = '   has parameters : '
    for name in parameter_names:
        string += f'{name} , '
    print(string)

    # print the type of this component
    print(f'   type =  {type(component)}')

    # does it contribute to the potential?
    print(f'   contributes to the potential? -  {component.contributes_to_potential}')

Information about component 0:
   name =  bh
   has parameters : m-bh , a-bh ,
   type =  <class 'dynamite.physical_system.Plummer'>
   contributes to the potential? -  True
Information about component 1:
   name =  dh
   has parameters : c-dh , f-dh ,
   type =  <class 'dynamite.physical_system.NFW'>
   contributes to the potential? -  True
Information about component 2:
   name =  stars
   has parameters : q-stars , p-stars , u-stars ,
   type =  <class 'dynamite.physical_system.TriaxialVisibleComponent'>
   contributes to the potential? -  True

Each component has a name, some parameters, and a type. Currently, dynamite only supports this exact combination of component types, i.e.

• one Plummer component representing the black hole

• one NFW component representing the dark halo (as of v1.0.0, DYNAMITE additionally supports Hernquist, TriaxialCoredLogPotential, and GeneralisedNFW dark halos, see documentation)

• one TriaxialVisibleComponent representing the stellar body of the galaxy

In the future, we want to support other components, and more flexible combinations of components.

For the stars - i.e. component 2 - we must provide some input data files. The location of these files is specified in the configuration file, at

settings -> io_settings -> input_directory

which takes the value,

c.settings.io_settings['input_directory']

'NGC6278_input/dynamite_input/'

The names of the following files are also specified in the configuration file, in the locations

system_components -> stars -> mge_pot
system_components -> stars -> mge_lum
system_components -> stars -> kinematics -> <kinematics name> --> datafile
system_components -> stars -> kinematics -> <kinematics name> --> aperturefile
system_components -> stars -> kinematics -> <kinematics name> --> binfile

which take values of the appropriate filenames. You are free to give these files whatever name you like, as long as it is specified in the configuration file.

Let’s have a look at the MGE,

c.system.cmp_list[2].mge_pot, c.system.cmp_list[2].mge_lum

(MGE({'name': None, 'datafile': 'mge.ecsv', 'input_directory': 'NGC6278_input/dynamite_input/', 'data': <Table length=6>
    I      sigma      q    PA_twist
 float64  float64  float64 float64
 -------- -------- ------- --------
 26819.14  0.49416 0.89541      0.0
  2456.39  2.04299 0.79093      0.0
    456.8  2.44313  0.9999      0.0
   645.49   6.5305 0.55097      0.0
    14.73 17.41488  0.9999      0.0
   123.85 21.84711 0.55097      0.0, 'logger': <Logger dynamite.mges.MGE (DEBUG)>}),
 MGE({'name': None, 'datafile': 'mge.ecsv', 'input_directory': 'NGC6278_input/dynamite_input/', 'data': <Table length=6>
    I      sigma      q    PA_twist
 float64  float64  float64 float64
 -------- -------- ------- --------
 26819.14  0.49416 0.89541      0.0
  2456.39  2.04299 0.79093      0.0
    456.8  2.44313  0.9999      0.0
   645.49   6.5305 0.55097      0.0
    14.73 17.41488  0.9999      0.0
   123.85 21.84711 0.55097      0.0, 'logger': <Logger dynamite.mges.MGE (DEBUG)>}))

The kinematics are stored here,

type(c.system.cmp_list[2].kinematic_data)

list

Note that this object has type list. This is because a single component can have multiple different sets of kinematics. In this example, the first (and only) entry in the list is

type(c.system.cmp_list[2].kinematic_data[0])

dynamite.kinematics.GaussHermite

We see that this kinematics object has type GaussHermite. This has also been specified in the configuration file, under

system_components -> stars -> kinematics --> kinset1 --> type

In addition to GaussHermite, DYNAMITE also supports BayesLOSVD kinematics (see documentation).

The kinemtic data itself can be accessed as follows,

c.system.cmp_list[2].kinematic_data[0].data

Table length=152

vbin_id	v	dv	sigma	dsigma	h3	dh3	h4	dh4
int64	float64	float64	float64	float64	float64	float64	float64	float64
1	-44.7518	2.0988	211.7829	1.9654	0.0	0.005	0.0	0.005
2	-69.2798	2.3608	188.883	2.1766	0.0	0.005	0.0	0.005
3	-70.8007	3.9338	177.095	3.7142	0.0	0.005	0.0	0.005
4	-47.9191	3.4617	184.1865	3.0684	0.0	0.005	0.0	0.005
5	-30.5932	2.386	208.2472	2.2207	0.0	0.005	0.0	0.005
6	-76.5504	5.5494	162.3144	4.572	0.0	0.005	0.0	0.005
7	-49.8285	4.4211	174.0705	3.6461	0.0	0.005	0.0	0.005
8	-35.7514	3.5335	180.28	3.3672	0.0	0.005	0.0	0.005
9	-75.4144	3.5041	139.9438	3.9235	0.0	0.005	0.0	0.005
10	-43.617	4.1029	175.4291	3.7975	0.0	0.005	0.0	0.005
11	-80.8035	4.942	150.4998	4.1311	0.0	0.005	0.0	0.005
12	-0.6747	3.1619	193.1385	2.6654	0.0	0.005	0.0	0.005
13	-27.9752	5.4621	173.8336	4.8414	0.0	0.005	0.0	0.005
14	-9.841	2.2988	216.2178	2.4602	0.0	0.005	0.0	0.005
15	6.9234	2.4226	211.2376	2.5402	0.0	0.005	0.0	0.005
...	...	...	...	...	...	...	...	...
137	88.5377	6.215	125.0646	5.1352	0.0	0.005	0.0	0.005
138	83.5245	3.5661	117.6799	3.7398	0.0	0.005	0.0	0.005
139	102.0751	4.1575	109.2731	4.4547	0.0	0.005	0.0	0.005
140	97.3282	4.9919	128.3343	4.4056	0.0	0.005	0.0	0.005
141	112.2471	4.6657	119.3019	5.0147	0.0	0.005	0.0	0.005
142	136.8294	4.5934	107.5522	4.9049	0.0	0.005	0.0	0.005
143	130.3479	8.6243	119.4998	4.3771	0.0	0.005	0.0	0.005
144	142.0869	13.3664	95.5229	6.7464	0.0	0.005	0.0	0.005
145	111.6049	4.8763	120.454	4.4527	0.0	0.005	0.0	0.005
146	120.9348	3.9434	107.9628	4.4038	0.0	0.005	0.0	0.005
147	123.2916	5.2241	110.2003	5.0809	0.0	0.005	0.0	0.005
148	146.468	5.5686	102.5112	4.2288	0.0	0.005	0.0	0.005
149	185.5138	6.2494	78.0394	4.6495	0.0	0.005	0.0	0.005
150	153.1489	4.2502	89.5659	4.0578	0.0	0.005	0.0	0.005
151	114.6925	5.827	100.1592	6.2472	0.0	0.005	0.0	0.005
152	197.743	3.9056	65.0177	5.4738	0.0	0.005	0.0	0.005

Creating a Schwarzschild model

Our next step will be to create a model i.e. a dyn.model.Model object. To help understand what the Model object is, let’s read the internal documentation (i.e. the docstring) for the class:

class Model(object):
    """A DYNAMITE model.

    The model can be run by running the methods (i) get_orblib, (ii) get_weights
    and (iii) (in the future) do_orbit_colouring. Running each of these methods
    will return the appropriate object, e.g. model.get_orblib() --> returns an
    OrbitLibrary object model.get_weights(...) --> returns a WeightSolver object

    Parameters
    ----------
    config : a ``dyn.config_reader.Configuration`` object
    parset : row of an Astropy Table
        contains the values of the potential parameters for this model
    directory : str
        The model directory name (without path). If None or not specified,
        the all_models_file will be searched for the directory name. If the
        all_models file does not exist, the model directory will be set to
        ``orblib_000_000/ml{ml}``.

    Returns
    -------
    Nothing returned. Attributes holding outputs are are added to the
    object when methods are run.

    """
    def __init__(self, config=None, parset=None, directory=None):

Looking at the init signature for this class, we see that a Model requires 3 input arguments. We’ve already met the first, the configuration object c.

The remaining input parameter we need to provide is 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 a parameter set specified by these values as follows,

parset = c.parspace.get_parset()
print(parset)

  m-bh    a-bh c-dh f-dh q-stars p-stars u-stars  ml
-------- ----- ---- ---- ------- ------- ------- ---
100000.0 0.001  8.0 10.0    0.54    0.99  0.9999 5.0

Note that, compared to values specified in the configuration file, parameters which have been specified as logarithmic, i.e. those configured as

parameters -> XXX -> logarithmic : True

have been exponentiated in this table. More details can be found in the tutorial/documentation parameter_space.ipynb.

With the first 2 input arguments at hand, we can now create our model object,

model = dyn.model.Model(config=c, parset=parset)

As this is our first model and we did not specify an explicit model directory, a standard directory name has been assigned to the model. Note that this path has not been created in the file system so far.

Having created the model object, we can now run it. First, let’s create directories for the output,

model.setup_directories()

This should have created a directory, where the name has been specified in the config file,

c.settings.io_settings['output_directory']

'NGC6278_output/'

inside of which you should find

NGC6278_output/models/

inside of which a unique directory for this particular model has been created, called

model.get_model_directory()

'NGC6278_output/models/orblib_000_000/ml5.00/'

The directory name is constructed from the standard name for new, directly created models and their ml value.

The next step is to calculate the orbit library. This step will take a few minutes,

model.get_orblib()

<dynamite.orblib.LegacyOrbitLibrary at 0x11c5b9460>

Having calculated an orbit library, we now need to find out which orbits are useful for reproducing the observations. This is an Non-Negative Least Squares (NNLS) optimization problem, which can be solved as follows

model.get_weights()

[INFO] 12:51:56 - dynamite.weight_solvers.LegacyWeightSolver - Using WeightSolver: LegacyWeightSolver
[INFO] 12:51:56 - dynamite.weight_solvers.LegacyWeightSolver - NNLS solution read from existing output

<dynamite.weight_solvers.LegacyWeightSolver at 0x11c5b9580>

Congratulations! You have run your first Schwarzschild model using DYNAMITE. The chi-squared of this model is,

model.chi2

36462.868022611096

Is that good? I don’t know! To find out, we’ll have to run more models and compare their chi-squared values. For information about this, see the tutorial running_a_grid_of_models.ipynb.

Plot the Models

Now let’s look at some output for the model that we have just run

plotter = dyn.plotter.Plotter(config=c)

figure = plotter.plot_kinematic_maps(model)

[INFO] 12:52:16 - dynamite.plotter.Plotter - Plotting kinematic maps for kin_set no 0: califa
/Users/leahskusa/Desktop/uni/dynamics_assistant/dynamite/dynamite_test/DYN_TEST/lib/python3.9/site-packages/plotbin-3.1.3-py3.9.egg/plotbin/sauron_colormap.py:105: UserWarning: Trying to register the cmap 'sauron' which already exists.
/Users/leahskusa/Desktop/uni/dynamics_assistant/dynamite/dynamite_test/DYN_TEST/lib/python3.9/site-packages/plotbin-3.1.3-py3.9.egg/plotbin/sauron_colormap.py:106: UserWarning: Trying to register the cmap 'sauron_r' which already exists.

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, \(V\), \(\sigma\), \(h_3\) and \(h_4\). We can see the following features in the fit,

• the model and data surface densities are very similar

• the sense of rotation of the \(v\) map is reproduced well, though the amplitude is lower than observed

• the \(\sigma\), \(h_3\) and \(h_4\) maps are less well reproduced

While the fit is certainly not perfect, it’s reassuring to see that some features are reproduced well. To improve the fit, we will have to explore parameter space more fully. See the tutorial on “running a grid of models” for more details.

Exercise

Change one of the potential parameters, then run another Schwarzschild model. You could do this by creating a new configuration file, or manually changing one of the parameters in this notebook. For example, to change the parameter ml - which is the mass-to-light ratio of the stars - to six, then you can uncomment the next line,

# parset['f'] = 6.

Change any parameter you would like, then re-run a model and save the output plots.