2.4. Inversion

Because inversion is often an ill-defined operation, getting a good inversion contains sometimes a large fraction of inspiration. In order to get a satisfactory fit, one usually needs to tune the number of nodes (for photospheric atmospheres), the weights of the Stokes parameters, the initial configuration, etc. So don’t expect to get a good inversion on your first run of Hazel. But if you keep insisting, you will get something useful out of the code.

2.4.1. Configuration file

A typical configuration file reads something like in the following. Remember that you have to adapt it to the specific details of your spectrum.

# Hazel configuration File

[Working mode]
Output file = output.h5
Number of cycles = 1
Save all cycles = False


[Spectral regions]
    [[Region 1]]
    Name = spec1
    Wavelength = 10826, 10833, 150
    Topology = ph1 -> ch1
    LOS = 0.0, 0.0, 90.0
    Boundary condition = 1.0, 0.0, 0.0, 0.0       # I/Ic(mu=1), Q/Ic(mu=1), U/Ic(mu=1), V/Ic(mu=1)
    Wavelength file = 'observations/10830.wavelength'
    Wavelength weight file = 'observations/10830.weights'
    Observations file = 'observations/10830_stokes.h5'
    Straylight file = None
    Mask file = None
    Weights Stokes I = 1.0, 1.0, 1.0, 1.0
    Weights Stokes Q = 0.0, 1.0, 1.0, 1.0
    Weights Stokes U = 0.0, 1.0, 1.0, 1.0
    Weights Stokes V = 1.0, 1.0, 1.0, 1.0

[Atmospheres]

    [[Photosphere 1]]
    Name = ph1
    Reference atmospheric model = 'photospheres/model_photosphere.1d'
    Spectral region = spec1
    Wavelength = 10826, 10833
    Spectral lines = 300,

        [[[Ranges]]]
        T      = -3000.0, 3000.0
        vmic   = 0.0, 3.0
        v      = -10.0, 10.0
        Bx     = -1000.0, 1000.0
        By     = -1000.0, 1000.0
        Bz     = -1000.0, 1000.0
        ff     = 0.0, 1.0
        vmac   = 0.0, 5.0

        [[[Nodes]]]
        T      = 3, 3, 5, 5
        vmic   = 1, 1, 1, 1
        v      = 1, 1, 1, 1
        Bx     = 1, 1, 1, 1
        By     = 1, 1, 1, 1
        Bz     = 1, 1, 1, 1
        ff     = 0, 0, 0, 0
        vmac   = 0, 0, 0, 0

        [[Regularization]]
        T      = None
        vmic   = None
        v      = None
        Bx     = None
        By     = None
        Bz     = None

    [[Chromosphere 1]]
    Name = ch1                                              # Name of the atmosphere component
    Spectral region = spec1                                 # Spectral region to be used for synthesis
    Height = 3.0                                            # Height of the slab
    Line = 10830                                            # 10830, 5876
    Wavelength = 10826, 10833                         # Wavelength range used for synthesis
    Reference atmospheric model = 'chromospheres/model_chromosphere.1d'    # File with model parameters

        [[[Ranges]]]
        Bx     = -500, 500
        By     = -500, 500
        Bz     = -500, 500
        tau    = 0.1, 2.0
        v      = -10.0, 10.0
        deltav = 3.0, 12.0
        beta   = 1.0, 2.0
        a      = 0.0, 1.0
        ff     = 0.0, 1.0


        [[[Nodes]]]
        Bx     = 0, 0, 1, 1
        By     = 0, 0, 1, 1
        Bz     = 0, 0, 1, 1
        tau    = 0, 0, 0, 0
        v      = 0, 0, 0, 0
        deltav = 0, 0, 0, 0
        beta   = 0, 0, 0, 0
        a      = 0, 0, 0, 0
        ff     = 0, 0, 0, 0

This configuration file defines a spectral region in the near infrared that contains a photosphere and a chromosphere, one after the other. The physical conditions of the photosphere are read from the photospheres/model_photosphere.1d file. Note that all files in this example are 1D so that the input atmosphere is then shared among all observed pixels. You might as well use 3D files (as described in input) with the same number of pixels as the input, and then a different initial atmosphere can be used for every pixel.

In inversion mode you need to define the observations and the weighting scheme, which is done through the Wavelength file, Wavelength weight file and Observations file, all of them explained in input.

2.4.2. Single inversions

2.4.2.1. Without randomization

If the files with the observations are 1D, you can simply carry out the inversion by invoking, assuming that the configuration file is saved on conf.ini:

mod = hazel.Model('configurations/conf_single.ini', working_mode='inversion')
mod.read_observation()
mod.open_output()
mod.invert()
mod.write_output(
mod.close_output()

We refer to Output files for a description of the output.

2.4.2.2. With randomization

Carrying out robust inversions is a difficult task. Randomization is a technique appropriate to capture the sensitivity of the result to the initial conditions. In the ideal case, one should obtain the same solution irrespectively of the initial values of the parameters. However, this will never be the case. The best one can do is to get the sensitivity of the results to the initial conditions. To this end, Hazel v2.0 provides a method to randomize the initial conditions. Just pass the maximum number of randomizations when instantiating the Model and then you can carry out inversions by looping over randomizations and saving each result to the output file. We prefer the user to deal with randomization because of the added flexibility.

mod = hazel.Model('configurations/conf_single.ini', working_mode='inversion', randomization=2)
mod.read_observation()
mod.open_output()
for loop in range(2):
    mod.invert(randomize=True)
    mod.write_output(randomization=loop)
mod.close_output()

2.4.3. Multiple inversions

For 3D cases, we recommend the user to use iterators. Otherwise, one needs to manually read the observations at every pixel.

2.4.3.1. Serial inversions

In serial mode, one can do the inversion of many pixels using 3D atmospheres by using the following strategy.

2.4.3.1.1. Without randomization

iterator = hazel.Iterator(use_mpi=False)
mod = hazel.Model('configurations/conf_nonmpi_inv1d.ini', working_mode='inversion', verbose=2)
iterator.use_model(model=mod)
iterator.run_all_pixels(start=0)

The first thing to do is to instantiate an Iterator which will be used to iterate over all pixels in the file with the observations. The Iterator class admits a single keyword use_mpi, which is set to False by default. If True, it will use the Message Passing Interface machinery, for which you need to have the mpi4py package installed. If False, it will serially iterate over all pixels. Then we instantiate the Model as usual and tell the Iterator to use this model. Finally, we run the calculation by calling the run_all_pixels method of the Iterator. run_all_pixels admit two parameters: start, set by default to 0 that marks the starting point of the inversions and end, set by default to None (equivalent to using the number of pixels in the file with the observations), that mark the last pixel to invert.

2.4.3.1.2. With randomization

The serial mode can also be used with randomization. Because of the complexity, the iterator takes fully care of randomizations and will simplify your life. It can be used simply by passing the randomization keyword during the Model instantiation.

iterator = hazel.Iterator(use_mpi=False)
mod = hazel.Model('configurations/conf_nonmpi_inv1d.ini', working_mode='inversion', verbose=2, randomization=2)
iterator.use_model(model=mod)
iterator.run_all_pixels()

2.4.3.2. Parallel inversions

Carrying out inversions of maps in serial mode is probably not the way to go. The desired strategy is to use computers with several nodes, which can work in parallel. The specific paralellization in Hazel v2.0 is such that practically a linear scaling with the number of nodes is sure. It uses a parent-worker architecture, in which a parent node reads the input and broadcasts it to the nodes, which are then in charge of carrying out the synthesis/inversions. When finished, the results are sent back to the parent, who saves the results and sends a new observation to the available worker.

2.4.3.2.1. Without randomization

Running an MPI work needs to be done by using mpiexec, so that you need to define the following lines in a script:

iterator = hazel.Iterator(use_mpi=True)
mod = hazel.Model('configurations/conf_mpi_invh5.ini', working_mode='inversion', rank=iterator.get_rank())
iterator.use_model(model=mod)
iterator.run_all_pixels()

The only difference with respect to the serial mode is that you need to pass the rank of each worker during the instantiation of the Model because the parent and the worker will need to do different things with the model. Then just run it with mpiexec, defining the number of nodes used in this work:

mpiexec -n n_nodes python inversion.py

Sometimes you need to let OpenMPI or MPICH use all available cores (this is important if your CPUs are using hyperthreading). To this end, write a file (for example mf) with the following content

localhost max_slots=16

changing the number to the appropriate number of available cores you have. Then run the code with:

mpiexec -machinefile mf -n n_nodes python inversion.py

2.4.3.2.2. With randomization

Randomization can also be used in parallel. For that, just follow the same approach as before and pass the randomization keyword to the Model:

iterator = hazel.Iterator(use_mpi=True)
mod = hazel.Model('configurations/conf_mpi_invh5.ini', working_mode='inversion', rank=iterator.get_rank(), randomization=2)
iterator.use_model(model=mod)
iterator.run_all_pixels()

2.4.4. Analysis of the results

We provide a simple way to analyze the results of the inversions. The following code reads the output file and plots the results of the inversion for a given pixel.

f = h5py.File('output.h5', 'r')
wl_syn, stokes_syn = hazel.analysis.synth_model('conf_single.ini', f, active=['ph1', 'ch1', 'ch2', 'te1'], verbose=3)

where we pass the configuration file, the output file struct and the active components to be computed. The output of this function is the synthetic spectrum and the synthetic Stokes profiles.