Note: This tutorial was generated from an IPython notebook that can be downloaded here.

Radex

from spectralradex import radex
from multiprocessing import Pool
import numpy as np
import time

The simplest use case for SpectralRadex is to be a simple python wrapper for RADEX. This allows large grids of RADEX models or complex parameter inference procedures to be run in an environment suited to those tasks.

If one wishes to run radex, we simply need a dictionary of the parameters RADEX expects. An example can be obtained using the get_default_parameters() function like so

params = radex.get_default_parameters()
print("{")
for key,value in params.items():
    print(f"\t{key} : {value}")
print("}")

{
    molfile : co.dat
    tkin : 30.0
    tbg : 2.73
    cdmol : 10000000000000.0
    h2 : 100000.0
    h : 0.0
    e- : 0.0
    p-h2 : 0.0
    o-h2 : 0.0
    h+ : 0.0
    linewidth : 1.0
    fmin : 0.0
    fmax : 1000.0
    geometry : 1
}

Note that each possible collisional partner has a separated entry for its density. You can check the collisional partners in your datafile with get_collisional_partners()

radex.get_collisional_partners("co.dat")

{'Number': 2, 'Partners': ['p-h2', 'o-h2']}

You only need to provide densities for the partners you wish to include in the calculation but you must include at least one of the partners. Two species cases are:

  • RADEX will sum the o-h2 and p-h2 density values to produce a h2 value.
  • A small number of datafiles have p-H2 collsions only and you may wish to place your total h2 density in that entry to approximate the o-h2 collisions

Once your parameter dictionary is set up, we pass that to the run() function.

output = radex.run(params)
output.head()
E_UP (K) freq WAVEL (um) T_ex tau T_R (K) POP UP POP LOW FLUX (K*km/s) FLUX (erg/cm2/s) Qup Qlow
0 5.53 115.271202 2600.757633 31.666252 0.000223 0.006275 0.246666 0.097917 0.006680 1.317591e-10 1 0
1 16.60 230.538000 1300.403656 29.262261 0.000735 0.017551 0.281677 0.246666 0.018683 2.947981e-09 2 1
2 33.19 345.795990 866.963374 26.640080 0.001112 0.021294 0.211510 0.281677 0.022667 1.207049e-08 3 2
3 55.32 461.040768 650.251515 24.363876 0.001022 0.015261 0.109663 0.211510 0.016246 2.050309e-08 4 3
4 82.97 576.267931 520.231028 22.798547 0.000605 0.007078 0.039845 0.109663 0.007535 1.856956e-08 5 4

Parameter Grids

It is more likely that one will want to run the code over many combinations of input parameters. This can be achieved via the run_grid() function. This function also takes a parameter dictionary of the same format as run(). However, variables which are too be varied over the grid should be supplied as iterables.

Furthermore, to keep things simple, the desired RADEXtakes iterables for the three variables (density, temperature and column density) as well as fixed values for the other RADEX parameters. It then produces the RADEX output for all combinations of the three iterables.

We’ll use an example grid which can be acquired using the get_example_grid_parameters() function.

parameters=radex.get_example_grid_parameters()
parameters

{'molfile': 'co.dat',
 'tkin': array([ 10. ,  82.5, 155. , 227.5, 300. ]),
 'tbg': 2.73,
 'cdmol': array([1.e+14, 1.e+15, 1.e+16, 1.e+17, 1.e+18]),
 'h2': array([   10000.        ,    56234.13251903,   316227.76601684,
         1778279.41003892, 10000000.        ]),
 'h': 0.0,
 'e-': 0.0,
 'p-h2': 0.0,
 'o-h2': 0.0,
 'h+': 0.0,
 'linewidth': 1.0,
 'fmin': 0.0,
 'fmax': 800.0,
 'geometry': 1}

tic = time.perf_counter()

grid_df = radex.run_grid(parameters,target_value="T_R (K)")
toc = time.perf_counter()
print(f"run_grid took {toc-tic:0.4f} seconds without a pool")

run_grid took 2.8573 seconds without a pool

grid_df.iloc[:,0:6].head()
tkin cdmol h2 (1)-(0)[115.2712018 GHz] (2)-(1)[230.538 GHz] (3)-(2)[345.7959899 GHz]
0 10.0 1.000000e+14 10000.0 0.114622 0.108152 0.022018
1 10.0 1.000000e+15 10000.0 1.048925 0.958338 0.215099
2 10.0 1.000000e+16 10000.0 5.189712 4.045272 1.567682
3 10.0 1.000000e+17 10000.0 6.561081 5.156221 3.411413
4 10.0 1.000000e+18 10000.0 6.639451 5.259944 3.822848

Parallelization

In order to be as flexible as possible, SpectralRadex has no built in multiprocessing. However, the run_grid() function does take the optional parameter pool which should be an object with map(), join(), and close() methods that allow functions to be evaluated in parallel. For example, the python standard multiprocessing.pool obect or Schwimmbad’s MPIPool.

If such an object is supplied, the grid will be evaluated in parallel. Note the time in the example below compared to the grid above.

tic = time.perf_counter()
pool=Pool(8)
grid_df = radex.run_grid(parameters,target_value="T_R (K)",pool=pool)
toc = time.perf_counter()
print(f"run_grid took {toc-tic:0.4f} seconds with a pool of 8 workers")
grid_df.iloc[:,0:6].head()

run_grid took 0.7338 seconds with a pool of 8 workers
tkin cdmol h2 (1)-(0)[115.2712018 GHz] (2)-(1)[230.538 GHz] (3)-(2)[345.7959899 GHz]
0 10.0 1.000000e+14 10000.0 0.114622 0.108152 0.022018
1 10.0 1.000000e+15 10000.0 1.048925 0.958338 0.215099
2 10.0 1.000000e+16 10000.0 5.189712 4.045272 1.567682
3 10.0 1.000000e+17 10000.0 6.561081 5.156221 3.411413
4 10.0 1.000000e+18 10000.0 6.639451 5.259944 3.822848