<p align="center"> <img src="docs/logo.jpg" alt="" width="512"/> </p>

mwydyn (Welsh: worm; pronounced: muy-din; IPA: [ˈmʊi̯dɪn])

About

mwydyn is a fully-automated multiple velocity component hyperfine line-fitting code.

Motivation: The code was written to perform automatic spectral-decomposition of N₂H⁺ (1-0) data cubes using a hyperfine model. We were finding that in much of our data, the 'isolated' component of N₂H⁺ (1-0) that is often used for Gaussian decomposition was heavily blended, and so we had to model the full hyperfine structure for effective decomposition. This code largely follows the hyperfine modeling implemented by GILDAS/CLASS, and is subject to the same assumptions and caveats. While, by default we allow up to 3 separate N₂H⁺ (1-0) velocity components per spectrum, in principle this is extensible to an arbitrary number, as well as other molecules or transitions with hyperfine structure.

Credit

mwydyn was developed by:

• Michael Anderson
• Andrew Rigby

With contributions from:

• Nicolas Peretto
• Gwenllian Williams
• Sarah Ragan

If you use mwydyn please consider citing Rigby et al. (2024) and Anderson et al. (in prep).

Performance statistics

• 57 spectra per minute per core (Apple M3 Pro, 2023)
• 33 spectra per minute per core (Apple M2, 2022)
• 10 spectra per minute per core (Intel Xeon E5-2640 v4 2.40GHz, 2016)

Dependencies

mwydyn was written and tested in python version 3.9.18. It requires the following packages, and has been tested on the version in parentheses.

• astropy (5.1.1)
• lmfit (1.3.2)
• matplotlib (3.8.2)
• numpy (1.26.4)
• packaging (23.2)
• tqdm (4.66.4)

For more up-to-date versions of python, (e.g. 3.11.x), some of the parallel processing implementation has changed, so for the time being we recommend running mwydyn in a 3.9 virtual environment if you wish to use that feature.

Installation

While we hope to eventually convert mwydyn into a pip-installable package, it is simple to use as a script. The recommended procedure is to install directly from this github page:

$ cd /your/project/directory
$ git clone https://github.com/mphanderson/mwydyn.git

Usage

The code can be run either from the terminal, or from within an ipyhton environment (recommended). You can specify either a config file, a data file, or both as command line arguments:

$ python mwydyn.py myconfig.ini mydata.fits

We recommend running the code from inside an ipython environment in order to use some of the interactive plotting features:

In [1]: run mwydyn.py myconfig.ini mydata.fits

This will run the code using the named config file and input data filename. If there is an input filename in the config, then it will be ignored and the command line argument will be used instead.

In [1]: run mwydyn.py myconfig.ini

This will run the code with the named config file, and use whatever input data filename is listed in the config. If there isn’t one in there, it will complain and quit.

In [1]: run mwydyn.py mydata.fits

This will run the code with the default config.ini file in the main directory and use the command line input data filename.

How it works

mwydyn runs with two basic steps:

1. Initial fitting
2. Spatial refinement

In the initial fitting step, the code will start off by cycling through each of the spectra in your data cube that satisfy the signal-to-noise ratio criterion snrlim, and attempts to fit up to N_max (defined in the config file) components per spectrum using the lmfit module. It devises some initial guesses for the fit parameters based on the location of the brightest channel and the spread in velocity of channels which are deemed to be detections. The code then determines which combination of components is the optimum, preferring fewer components unless as significant improvement in the BIC value (our goodness of fit metric) is obtained with more components.

Until this point, each spectrum has been treated independently. The spatial refinement step allows some level of communication between neighbouring pixels, which will not be independent if your pixels are sampling your telescope beam sensibly. The spectra are cycled through again up to a number of times specified by the CleanIter argument, and neighbouring spectra are inspected for better fits. If a local spectrum has a better fit, then we now attempt to refit the current spectrum using the fit parameters of the better fit as the initial guesses. If a better fit is obtained, we keep the new fit, or else keep the old one. This cycling allows the best fit parameters to percolate through the data cube. If later iterations of spatial refinement are making no difference, this step is terminated early to save time.

Demonstration:

mwydyn is packaged with some small test N₂H⁺ (1-0) data cubes from ALMA (Anderson et al. in prep.) and NOEMA (Rigby et al. 2024) to demonstrate the functionality. We include config_demo_alma.ini and config_demo_noema.ini to demonstrate how mwydyn works on these cubes. Let's try fitting the ALMA cube:

In [1]: run mwydyn.py config_demo_alma.ini

 Input data not specified in command line, using configuration file entry.

 ======================= Running mwydyn ===========================
 File: SDC326_n2h+_TP+7m+12m.subcube.fits
 Config file: config_demo_alma.ini
 Maximum no. components: 3
 Minimum S/N: 10.0
 Critical BIC value: 20.0
 Initial FWHM guess: 0.5
 Initial tau guess: 0.2
 Fitting method: leastsq
 Fit constraints: True
 Parallel processing on 4 processors
 Cleaning iterations: 10
 Cleaning radius: 2 pixels

 Fitting each spectrum : 100%|██████████████████████████████████████████| 256/256 [01:09<00:00,  3.70it/s]
 Spatial refinement    :  40%|█████████████████▍                           | 4/10 [00:21<00:33,  5.65s/it]
 No new results on refinement iteration 4: stopping.

 Data products saved to: /your/project/directory/mwydyn/data_products/

 * Click on a pixel to inspect the spectrum and fit

 Finished in 1.69 minutes
 ==================================================================

This has run mwydyn over the data cube, and stored the products in the data_products/ directory. If running in iPython, some interactive figures pop up at completion. The first is the summary figure:

<p align="center"> <img src="docs/summary.jpg" alt="" width="776"/> </p>

Which gives some quick-look images comparing, for example, the integrated intensity of the input cube (top left) with the integrated intensity of the resulting model cube (top middle left), and the number of fitted components per spectrum (top right).

In this figure, you can also click on individual pixels in the various maps to inspect the fits to the spectra. For example, the model for the spectrum at pixel coordinates (8, 7) is:

<p align="center"> <img src="docs/spectrum.jpg" alt="" width="512"/> </p>

Which has been fitted with 2 components, one broad (FWHM = 3.7 km/s) and one narrow (1.0 km/s).

Secondly, a figure showing the results of the spatial refinement step appears:

<p align="center"> <img src="docs/spatial.jpg" alt="" width="776"/> </p>

This figure shows the integrated intensity of the model cube after each step of spatial refinement. The top row shows the integrated intensity of the model at each step compared to the original data on the left. The second row shows the difference between the model and data for the relevant stage, and the fourth row shows the residual of the model refined at each iteration compared to the previous one. The third row shows the corresponding maps of BIC (lower values = better fit). Two figures in the left column show a map of the number of refits used per spectrum and the iteration at which those pixels were last refined: in this case, no spectrum was refit more than once. These allow the user to check how the refits percolate across the map, and give an idea of whether you might need to increase the value of CleanIter. Although we ran with CleanIter = 10 in this run, no further refinement was made after the third iteration, and so the code finished early.

Configuring mwydyn

The main code is mwydyn.py, which should be called with either a configuration file as an argument and, optionally, a filename for the cube you wish to fit. Most of the behaviour is controlled by a configuration file. By default, it is configured for fitting N₂H⁺ (1-0) spectra which should run very smoothly. For other species, some minimal modification will be required.

The variables that should be defined to run and customise the code are stored in a configuration (.ini) file which, by default, will be the supplied file called config.ini. The configuration file is organised in five blocks: [paths], [input], [processing], [fitting], and [output], and here we summarise what each of the variables do. The default values are listed in square brackets.

Paths

• path - The path to main directory where mwydyn.py is located. A value of default will attempt to do this automatically. [default]
• input_dir - The location of the input data directory relative to path/. [data_in/]
• prod_dir - The location to store the output of mwydyn relative to path/. [data_products/]
• model_dir - The location of the model files for the molecule relative to path/. [models/]

Input

• input_fn - The name of the .fits file on which mwydyn will be run, or set to user to input your filename as a command line argument. [