

# Plotting of the Retrieval Results

To plot retrieval results, run the following command from a terminal window:
The routine uses the same YAML configuration file as the retrieval itself. The file was copied to the output directory. Calling the plotting routine must include:
- **res_dir /dir_of_output_folder/**: is a required client argument that points the plotting routine to the folder where the retrieval outputs were saved. (output_folder variable in the configuration file)

```
python /dir_to_pyRetLife_folder/scripts/main_plotting.py --res_dir /dir_of_output_folder/
```

In the following, we outline the main workflow for plotting the retrieval results followed in the `main_plotting.py` file. The same order should be followed in case you want to create your own routine, which uses the retrieval plotting routine, please import the following dictionary
```
from pyretlife.retrieval_plotting.run_plotting import retrieval_plotting_object
```

The retrieval plotting routine is organized in a class structure, which inherits from the retrieval_object class.

## Initialization

In order to initialize an instance of the retrieval plotting class the following parameters need to be specified:

- **results_directory [str]:** This parameter is a requirement and indicates the directory of the folder that contains the output files from the retrieval run.

When initializing of the retrieval_plotting_object, the following steps are performed:
    - The YAML file in the results directory is read in exactly the same way as when running a retrieval.
    - Unit conversions from the input units to the retrieval units are performed.
    - Assign the known values like when running the retrieval.

The posterior distribution is loaded from the `post_equal_weights.dat` file into a pandas dataframe. The dataframe can be accessed via the ‘posteriors’ class variable. Each column contains the posterior for a different retrieved parameter. The column name is the same as the one used in the YAML file.
A new folder `/dir_of_output_folder/Plots/` is generated where all the files generated by the plotting routine are saved.

```
plotting = retrieval_plotting_object(results_directory = '/dir_of_output_folder/')

```

After initializing the `retrieval_plotting object`, we need to calculate the P–T profiles and the Spectra corresponding to the posterior distribution. There are two functions to perform these tasks. Both the P-T profiles and the Spectra are calculated on multiple CPUs to reduce the calculation time required.

- **`calculate_posterior_pt_profile`**: Calculates the P–T profiles corresponding to the posterior distribution. The calculated P–T profiles are saved in the `/dir_of_output_folder/Plots/Ret_PT_Skip_#.pkl`. The following parameters can be specified:
    - **skip [int] (Default: 1)**: To reduce computation time, points in the posterior distribution can be skipped when calculating the P–T profiles. The skip value is specified in the name of the file where the calculated P–T profiles are saved. The default skip value of 1 means that P–T profiles for all points in the posterior distribution are calculated. A skip value of 100 would mean that the P–T profile for every 100th point in the posterior is calculated.
    - **n_processes [int] (Default: 50)**: The number of CPUs to use for the P-T profile calculation. Please ensure that n_processes is not larger than the number of CPUs that your machine has.
    - **reevaluate_PT [bool] (Default: False)**: In general, the plotting routine checks if a `/dir_of_output_folder/Plots/Ret_PT_Skip_#.pkl` exists (i.e. P–T profiles have been calculated previously) and loads the calculated profiles. The P–T profiles can be force-recalculated by setting the reevaluate_PT flag to True.
    - **layers [int] (Default: 500)**: Number of layers to evaluate the P–T profiles on. To ensure smooth P–T profile plots, this value should be significantly higher than the number of layers specified in the configuration file.
    - **p_surf [float] (Default: 4)**: The high-pressure limit in (log10([bar])) up to which the P–T profiles are calculated. The value should be significantly larger than the retrieved surface or cloud-top pressure.

- **`calculate_posterior_pt_spectrum`**: Calculates the spectra corresponding to the posterior distribution. The calculated spectra are saved in the `/dir_of_output_folder/Plots/Ret_Spec_Skip_#.pkl`. The following parameters can be specified:
    - **skip [int] (Default: 1)**: To reduce computation time, points in the posterior distribution can be skipped when calculating the spectra. The same applies as for the skip variable in the P–T profile calculation routine above.
    - **n_processes [int] (Default: 50)**: The number of CPUs to use for the P-T profile calculation. Please ensure, that n_processes is not larger than the number of CPUs that your machine has.
    - **reevaluate_PT [bool] (Default: False)**: In general, the plotting routine checks if a `/dir_of_output_folder/Plots/Ret_Spec_Skip_#.pkl` exists (i.e. spectra have been calculated previously) and loads the calculated spectra. The spectra can be force-recalculated by setting the reevaluate_spectra flag to True.

```
plotting.calculate_posterior_pt_profile(skip=1,n_processes=50,reevaluate_PT=False,layers=500,p_surf=4)

plotting.calculate_posterior_spectrum(self,skip=1,n_processes=50,reevaluate_spectra=False)

```
Optionally, the spectra and P–T profiles corresponding to the provided true values can also be calculated. For this to work, true values for all retrieved parameters have to be provided in the configuration file. If not all are provided, the calculation is skipped.

- **`calculate_true_pt_profile`**: If all true values are provided, the corresponding P-T profile is calculated. If one or more are missing, the calculation is skipped.

- **`calculate_true_spectrum`**: If all true values are provided, the corresponding spectrum is calculated. If one or more are missing, the calculation is skipped.


```
plotting.calculate_true_pt_profile()
plotting.calculate_true_spectrum()
```


There are functions to derive parameters from the posteriors. It is important to call these functions only after the spectra and the P–T profiles corresponding to the posterior distributions have been calculated/loaded, since the functions rely on these quantities. The calculated values are appended to the posteriors variable of the retrieval_plotting_object.
- **`deduce_bond_albedo`**: Calculates equilibrium temperature and Bond albedo from the retrieved planet radius and the posterior spectra. The calculated values are saved in the `/dir_of_output_folder/Plots/Ret_Bond_Albedo.pkl`. It requires:

    - **stellar_luminosity [float]**: Luminosity of the host star in solar luminosities (e.g. for the Sun the value would be 1.0). Has to be greater than 0.
    - **error_stellar_luminosity [float]**: Uncertainty on stellar_luminosity (e.g. a value of 0.05 means that the error is 0.05 solar luminosities). Has to be greater than 0.
    - **planet_star_separation [float]**: semimajor axis of the planet’s orbit around its host star in astronomical units (e.g. for the Earth the value would be 1.0).
    - **error_planet_star_separation [float]**: Uncertainty on planet_star_separation (e.g. a value of 0.05 means that the error is 0.05 astronomical units). Has to be greater than 0.
    - **reevaluate_bond_albedo [float] (Default: None)**: In general, the plotting routine checks if a `/dir_of_output_folder/Plots/Ret_Bond_Albedo.pkl` exists (i.e. equilibrium temperatures and Bond albedos have been calculated previously). The values can be force-recalculated by setting the reevaluate_bond_albedo flag to True.
    - **true_equilibrium_temperature [float] (Default: None)**: Allows to pass a true value for the equilibrium temperature to the retrieval plotting routine.
    - **true_bond_albedo [float] (Default: None)**: Allows to pass a true value for the Bond albedo to the retrieval plotting routine.

- **`deduce_gravity`**: routine that calculates the surface gravity from the retrieved planetary radius and mass. It requires:
    - **true_gravity [float] (Default: None)**: Allows to pass a true value for the surface gravity to the retrieval plotting routine.
- **`deduce_surface_temperature`**: routine that calculates the surface temperature from the retrieved P–T profiles and the retrieved surface pressure. It requires:
    - **true_surface_temperature [float] (Default: None)**: Allows to pass a true value for the surface temperature to the retrieval plotting routine.

```
plotting.deduce_bond_albedo(stellar_luminosity,
                            error_stellar_luminosity,
                            planet_star_separation,
                            error_planet_star_separation,
                            reevaluate_bond_albedo=False,
                            true_equilibrium_temperature = None,
                            true_bond_albedo = None)

plotting.deduce_gravity(true_gravity = None)

plotting.deduce_surface_temperature(true_surface_temperature = None)
```

There are the following routines to visualize the retrieved posterior distributions:

- **Posteriors**: plots histograms/corner plot of the retrieved parameters. It requires:
    - **save [bool] (Default: False)**: If ‘True’, the created plot is saved in the ‘/dir_of_output_folder/Plots/’ folder. Otherwise, the figure (figure and the axes) are returned by the function and can be modified by the user. 
    - **plot_corner [bool] (Default: False)**: It ‘True’ a corner plot is created. If ‘False’, the posterior of each parameter is plotted in a 1d histogram separately.
    - **log_pressures [bool] (Default: True)**: If ‘True’, the posteriors of all parameters with pressure units are plotted in logspace. 
    - **log_mass [bool] (Default: True)**: If ‘True’, the posteriors of all parameters with mass units are plotted in logspace.
    - **log_abundances [bool] (Default: True)**: If ‘True’, the posteriors of all abundance parameters are plotted in logspace.
    - **log_particle_radius [bool] (Default: True)**: If ‘True’, the posteriors of the cloud particle radii are plotted in logspace.
    - **plot_pt [bool] (Default: True)**: If ‘True’, the posteriors of the P–T profile parameters are plotted.
    - **plot_physparam [bool] (Default: True)**: If ‘True’, the posteriors of the physical parameters are plotted.
    plot_clouds [bool] (Default: True): If ‘True’, the posteriors of the parameters describing the atmospheric clouds are plotted.
    - **plot_chemcomp [bool] (Default: True)**: If ‘True’, the posteriors of the chemical composition parameters are plotted.
    - **plot_scatt [bool] (Default: True)**: If ‘True’, the posteriors of the scattering parameters are plotted.
    - **plot_moon [bool] (Default: True)**: If ‘True’, the posteriors of the moon parameters are plotted.
    - **plot_secondary_parameters [bool] (Default: True)**: If ‘True’, the posteriors of the derived parameters are plotted.
    - **bins [int] (Default: 20)**: Specifies the number of bins for the 1d histograms.
    - **quantiles1d [list] (Default: [0.16, 0.5, 0.84])**: Specifies the quantiles of the parameter posteriors, which are plotted as dotted vertical lines in the plot.
    - **color [str] (Default: 'k')**: colours in which the posteriors will be plotted. Any colour valid in pyplot.matplotlib will work also RGB colours.
    - **add_table [bool] (Default: False)**: If True‘, a table summarizing the retrieved values of all parameters is added to the top right corner of the corner plot.
    - **color_truth [str] (Default: 'C3')**: colours in which the provided true values will be plotted. Any colour valid in pyplot.matplotlib will work also RGB colours.
    - **ULU_lim [list] (Default: [-0.15,0.75])**: 
      ```{warning}
  Documentation work in progress. 
  ```
    - **parameter_units='input' [str, dict]**: Units for the parameters that are to be used in the plots. If ‘input’, the units specified in the configuration file are used. Otherwise, the standard retrieval units are used. Alternatively, a dictionary can be passed to specify units for individual parameters. The key specifies the name of the parameter as defined in the configuration file. The corresponding value should be an astropy.units unit.
    - **custom_unit_titles [dict] (Default: {})**: Allows the user to provide custom titles for the units of specific parameters. The key specifies the name of the parameter as defined in the configuration file. The corresponding value should be a string (TeX format) describing the unit of the parameter. If no unit title is provided, it will be generated from the astropy unit automatically.
    - **custom_parameter_titles [dict] (Default: {})**: Allows the user to provide custom titles for specific parameters. The key specifies the name of the parameter as defined in the configuration file. The corresponding value should be a string (TeX format) describing the name of the parameter. If no unit title is provided, the name from the configuration file will be used

```
Posteriors(save=False,
           plot_corner=True,

           log_pressures=True,
           log_mass=True,
           log_abundances=True,
           log_particle_radii=True,

           plot_pt=True,
           plot_physparam=True,
           plot_clouds=True,
           plot_chemcomp=True,
           plot_scatt=True,
           plot_moon=False,
           plot_secondary_parameters=True,

           parameter_units='input',
           custom_unit_titles={},
           custom_parameter_titles={}

           bins=20,
           quantiles1d=[0.16, 0.5, 0.84],
           color='k',
           add_table=False,
           color_truth='C3',
           ULU_lim=[-0.15,0.75])
```
       
There are the following routines to visualize the spectra corresponding to the retrieved posterior distributions of the parameters:
- **plot_retrieved_flux**: plots spectra corresponding to the retrieved parameter posterior distributions. It requires:
    - **ax [matplotlib.pyplot.axes] (Default: None)**: The subplot object.
    - **color [str] (Default: 'C2')**: colours in which the spectra corresponding to the retrieved posterior distributions will be plotted. Any colour valid in matplotlib.pyplol will work also RGB colours.
    - **case_identifier [string] (Default: None)**: Option to print a text specifying the retrieval run in the plot.
    - **plot_noise [bool] (Default: False)**: specifies if the noise on the spectrum that was used as input for the retrieval should be shown in the spectrum plot.
    - **plot_true_spectrum [bool] (Default: False)**: Specifies whether the input spectrum used in the retrieval should be added to the spectrum plot.
    - **plot_datapoints [bool] (Default: False)**: If ‘True’ data point with error bars for the input spectrum used in the retrieval will be added to the plot.
    - **quantiles [list] (Default: [0.05,0.15,0.25,0.35,0.65,0.75,0.85,0.95])**: A list that specifies the quantiles of the spectra corresponding to the retrieved parameter posteriors that are shown in the plot. The number of quantiles specified must always be even, and each value must lie between 0 and 1. A symmetric choice of quantiles is recommended (e.g. [0.05,...,0.95]).
    - **quantiles_title [list] (Default: None)**: A list specifying the titles for the plotted quantiles as strings. These titles will be used in the legend of the plot. The list must have half as many elements as the quantiles list.                                   
    - **figsize [tuple] (Default: (12,2))**: size of the generated matplotlib.pyplot figure. 
    - **legend_loc [str] (Default: 'best')**: specifies where the legend is to be placed within the plot. the same strings as for the matplotlib.pyplot.legend function are allowed.
    - **x_lim [list] (Default=None)**: Specifies the x-range plotted in the figure. 
    - **y_lim [list] (Default=None)**: Specifies the y-range plotted in the figure.              
    - **noise_title [str] (Default: 'Observation Noise')**: the title for the noise, which will be printed in the legend of the flux plot.  
    - **plot_instruments_separately [bool] (Default: False)**: If True, a separate plot for each instrument is created. Additionally, the flux of spectra that are generated by the retrieval are automatically rebinned to the input spectrum.
    - **plot_residual [bool] (Default: False)**: If True, the residual for the retrieved flux relative to the input spectrum is plotted. Otherwise, the absolute spectrum is plotted. 
    - **plot_log_wavelength [bool] (Default: False)**: If True, the x-axis of the plot (wavelengths) will be plotted in log-scale.
    - **plot_log_flux [bool] (Default: False)**: If True, the y-axis of the plot (planet flux) will be plotted in log-scale.
    - **plot_unit_wavelength [astropy.unit] (Default = None)**: Manually set the unit of the x-axis in the plot. If not provided, standard retrieval units are used.
    - **plot_unit_flux [astropy.unit] (Default = None)**: Manually set the unit of the y-axis in the plot. If not provided, standard retrieval units are used.
    - **plot_retrieved_median [bool] (Default: False)**: If True, the median of the retrieved spectra is plotted in addition to the specified quantiles.

```
plot_retrieved_flux(ax = None,
                color='C2',
                case_identifier = None,
                plot_noise = False,
                plot_true_spectrum = False,
                plot_datapoints = False,
                quantiles=[0.05,0.15,0.25,0.35,0.65,0.75,0.85,0.95],
                quantiles_title = None,        
                figsize=(12,2),
                legend_loc = 'best',
                x_lim = None,
                y_lim = None,                    
                noise_title = 'Observation Noise',  
                plot_instruments_separately=False,     
                plot_residual = False,
                plot_log_wavelength=False,
                plot_log_flux=False,
                plot_unit_wavelength=None,
                plot_unit_flux=None,
                plot_retrieved_median=False)
```

There are the following routines to visualize the P–T profiles corresponding to the retrieved posterior distributions of the parameters:
- **plot_retrieved_flux**: plots spectra corresponding to the retrieved parameter posterior distributions. It requires:
    - **save [bool] (Default: False)**: If ‘True’, the created plot is saved in the `/dir_of_output_folder/Plots/` folder. Otherwise, the figure (figure and the axes)
    - **x_lim [list] (Default: [0,1000])**: Specifies the x-range plotted in the figure.
    - **y_lim [list] (Default: [1e-6,1e4])**: Specifies the y-range plotted in the figure.
    - **quantiles [list] (Default: [0.05,0.15,0.25,0.35,0.65,0.75,0.85,0.95])**:  A list that specifies the quantiles of the P–T profiles corresponding to the retrieved parameter posteriors that are shown in the plot. The number of quantiles specified must be even, and each value must lie between 0 and 1. A symmetric choice of quantiles is recommended (e.g. [0.05,...,0.95]).
    - **quantiles_title [list] (Default: None)**: A list specifying the titles for the plotted quantiles as strings. These titles will be used in the legend of the plot. The list must have half as many elements as the quantiles list. 
    - **inlay_loc [str] (Default: 'upper right')**: specifies where the small inlay plot for the surface pressure should be located within the plot. The same strings as for the matplotlib.pyplot.legend function are allowed.
    - **bins_inlay [int] (Default: 20)**: Specifies the number of bins used for the 2d histogram in the surface conditions inlay plot.
    - **x_lim_inlay [list] (Default: None)**: Specifies the x-range plotted in the inlay plot. If none is provided, it will be automatically chosen.
    - **y_lim_inlay [list] (Default: None)**: Specifies the y-range plotted in the inlay plot. If none is provided, it will be automatically chosen.
    - **figure [matplotlib.pyplot.figure] (Default: None)**: The figure object.
    - **ax [matplotlib.pyplot.axes] (Default: None)**: The subplot object.
    - **color [string] (Default: 'C2')**: colours in which the P–T profiles corresponding to the retrieved posterior distributions will be plotted. Any colour valid in matplotlib.pyplot will work also RGB colours.
    - **case_identifier [string] (Default: '')**: Option to print a text specifying the retrieval run in the plot.
    - **legend_n_col [int] (Default: 2)**: Number of columns in the legend of the plot.
    - **legend_loc [str] (Default: 'best')**: Specifies where the figure legend should be located within the plot. The same strings as for the matplotlib.pyplot.legend function are allowed.
    - **figsize [tuple] (Default: (6.4, 4.8))**: Option to change the size of the generated matplotlib.pyplot figure.
    - **h_cover [float] (Default: 0.45)**: Percentage of the figure width that the small inlay plot covers. Values between 0 and 1 (1 stands for 100 percent).
    - **true_cloud_top [list] (Default: [None,None])**: Option to pass the true position of the cloud top to the retrieval plotting routine. The first list element is the pressure, the second the temperature at which the cloud top is located.
    - **plot_truth [bool] (Default: False)**: Option to plot the true pt profile in the plot.
    - **plot_residual [bool] (Default: False)**: If True, the residual for the retrieved flux relative to the input spectrum is plotted. Otherwise, the absolute spectrum is plotted. 
    - **plot_clouds [bool] (Default: False)**: Option to plot the retrieved position of the atmospheric clouds in the P–T profile in the plots.
    - **plot_unit_temperature [astropy.unit] (Default = None)**: Manually set the unit of the x-axis in the plot. If not provided, standard retrieval units are used.
    - **plot_unit_pressure [astropy.unit] (Default = None)**: Manually set the unit of the y-axis in the plot. If not provided, standard retrieval units are used.

```
plot_retrieved_pt_profile(save=False,
                   x_lim =[0,1000],
                   y_lim = [1e-6,1e4],      
                   quantiles=[0.05,0.15,0.25,0.35,0.65,0.75,0.85,0.95],
                   quantiles_title = None,
                   inlay_loc='upper right',
                   bins_inlay = 20,
                   x_lim_inlay =None,
                   y_lim_inlay = None,
                   figure = None, 
                   ax = None,
                   color='C2',
                   case_identifier = '',
                   legend_n_col = 2,
                   legend_loc = 'best',
                   n_processes=50,
                   figsize=(6.4, 4.8),
                   h_cover=0.45,
                   reevaluate_PT = False,
                   true_cloud_top=[None,None],

                   plot_truth = False,
                   plot_residual = False,
                   plot_clouds = False,
                   plot_unit_temperature=None,
                   plot_unit_pressure=None,)

```