Changes to plot multimodal posteriors when the option is enabled. To …

[yveskini]

…enable the multi_mode option turn multi_mode=True when loading the runs

[lmauviard]

Good and works for me. I suggest to keep the default xpsi loaded file so the user can still access it without reloading everything

[thjsal]

"defaul" should be "default" here. But I agree @lmauviard that keeping the default X-PSI file could be handy if want to switch to showing the normal corner plots and not wanting to re-load the runs. But if having that is complicated or breaks things, I am fine for leaving it out.

[thjsal]

Otherwise it looks good to me and seems to pretty much work for in my tests... I just have 1 out of my 4 runs that sometimes complains that KeyError: 'No run with ID matching request.' . But the same run works if I just only load that one, so maybe I do something wrong when loading many runs (it used to work before though).

One easy thing that could be added is a docstring to the load_runs function that would explain the new input arguments.

[thjsal]

Well, I actually I think there is a bug in this code. If I first load a run with 3 modes for one model and then load another run, with 4 modes, for another model, I get the above-mentioned KeyError when trying to only plot the run (and the model) that has 3 modes. I suspect it for some reason thinks that even that run should have 4 modes. But I am not sure yet what in the code is making that to happen. Here is the full error message:

Executing posterior density estimation...
Curating set of runs for posterior plotting...
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[29], line 3
      1 NC_flag_space = NC_flag
----> 3 _ = pp.plot(
      4      params=['radius', 'mass'], 
      5      IDs=OrderedDict([(root_filename_c3, ['run 1a']),]), 
      6      prior_density=False,#True,
      7      KL_divergence=True,
      8      ndraws=5e4,
      9      combine=True, combine_all=True, only_combined=False, overwrite_combined=True,
     10      param_plot_lims={'mass': (1.0,1.61)},#,'radius': (10.0,14.0)},#,'compactness': (0.085,0.16)},
     11      bootstrap_estimators=False,
     12      bootstrap_density=False,
     13      n_simulate=200,
     14      crosshairs=False,
     15      ext=extension,
     16      maxdots=3000,
     17      credible_interval_1d=True,
     18      annotate_credible_interval=True,
     19      credible_interval_1d_all_show=True,
     20      compute_all_intervals=False,
     21      sixtyeight=True,
     22      axis_tick_x_rotation=45.0,
     23      num_plot_contours=3,
     24      subplot_size=4.0,
     25      legend_corner_coords=(0.565,0.9),
     26      legend_frameon=False,
     27      scale_attrs=OrderedDict([('legend_fontsize', 1.9),
     28                               ('axes_labelsize', 1.9),
     29                               ('axes_fontsize', 'axes_labelsize'),
     30                              ]
     31                             ),
     32      colormap='Reds',
     33      shaded=True,#False,
     34      shade_root_index=-1,
     35      rasterized_shade=True,
     36      no_ylabel=True,
     37      no_ytick=True,
     38      lw=2.0,
     39      lw_1d=2.0,
     40      filled=False,
     41      normalize=True,
     42      veneer=True,
     43      tqdm_kwargs={'disable': False},
     44      lengthen=2.0,
     45      embolden=1.0,
     46      nx=500,
     47      scale_ymax=1.1,
     48      write=False,
     49      directory=dirname,
     50      ci_gap=0.13,
     51      precisions=[2,2],#[2,3,2],
     52      line_colors=['darkred', 'navy', 'darkorange', 'green','black','deeppink','dimgrey'],
     53      priors_identical=False,
     54      force_draw=[False,False,False,False],
     55      prior_samples_fnames = ["prior_samples_NicerxXMM_PDTU_prmcosi2_logg_R1014.npy",],
     56 )
     58 #dpi2 = 3000.0/(_.settings.subplot_size_inch*2)
     59 #print("dpi2=",dpi2)
     60 #plt.savefig(dirname+"/"+root_filename_c4+'_spacetime_MR_main.png',dpi=dpi2,bbox_inches='tight') 

File ~/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-2.2.7-py3.12-linux-x86_64.egg/xpsi/PostProcessing/_global_imports.py:92, in fix_random_seed(func, instance, args, kwargs)
     90 state = _np.random.get_state()
     91 _np.random.seed(random_seed)
---> 92 output = func(*args, **kwargs)
     93 _np.random.set_state(state)
     94 return output

File ~/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-2.2.7-py3.12-linux-x86_64.egg/xpsi/utils.py:32, in make_verbose.<locals>.decorator(func, instance, args, kwargs)
     30         final = msg  # catch last yield if generator
     31 else:
---> 32     final = func(*args, **kwargs)
     33 if _verbose and not deactivate_verbosity:
     34     if exit_msg and isinstance(exit_msg, _six.string_types):

File ~/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-2.2.7-py3.12-linux-x86_64.egg/xpsi/PostProcessing/_corner.py:163, in CornerPlotter.plot(self, params, IDs, combine, combine_all, only_combined, force_combine, overwrite_combined, bootstrap_estimators, bootstrap_density, separate_plots, write, root_filename, directory, ext, dpi, maxdots, **kwargs)
     41 @fix_random_seed
     42 @make_verbose('Executing posterior density estimation',
     43               'Posterior density estimation complete')
   (...)
     60          maxdots=2000,
     61          **kwargs):
     62     """ Generate posterior density plots.
     63 
     64     Invokes :mod:`getdist` and :mod:`nestcheck` for nested sampling runs.
   (...)
    161 
    162     """
--> 163     self.set_subset(IDs, combine, combine_all,
    164                     force_combine, only_combined,
    165                     overwrite_combined)
    166     self.set_params(params)
    168     if bootstrap_density and not separate_plots:

File ~/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-2.2.7-py3.12-linux-x86_64.egg/xpsi/utils.py:32, in make_verbose.<locals>.decorator(func, instance, args, kwargs)
     30         final = msg  # catch last yield if generator
     31 else:
---> 32     final = func(*args, **kwargs)
     33 if _verbose and not deactivate_verbosity:
     34     if exit_msg and isinstance(exit_msg, _six.string_types):

File ~/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-2.2.7-py3.12-linux-x86_64.egg/xpsi/PostProcessing/_postprocessor.py:185, in PostProcessor.set_subset(self, IDs, combine, combine_all, force_combine, only_combined, overwrite)
    179     print('Warning: Only the first eight positional runs will be '
    180           'plotted individually, with IDs to prevent the plot from overcrowding '
    181           + ', '.join(ids[:8]))
    183     ids = ids[:8]
--> 185 posterior.set_subset(ids,
    186                      combine = combine,
    187                      combine_all = combine_all,
    188                      force_combine = force_combine,
    189                      only_combined = only_combined,
    190                      only_principal = False,
    191                      overwrite = overwrite)

File ~/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-2.2.7-py3.12-linux-x86_64.egg/xpsi/PostProcessing/_runs.py:206, in Runs.set_subset(self, IDs, combine, combine_all, force_combine, only_combined, only_principal, overwrite)
    204                 for mode in range(self.mode_len):
    205                     IDs.append(root_ID + f"{self.mode_label} {mode+1}")
--> 206     self._subset = [self[ID] for ID in IDs]
    208 if combine and force_combine: # create new run object
    209     self._combine(combine_all, overwrite)

File ~/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-2.2.7-py3.12-linux-x86_64.egg/xpsi/PostProcessing/_runs.py:311, in Runs.__getitem__(self, ID)
    308     return self._runs[ID]
    310 # if get this far there was lookup error
--> 311 raise KeyError('No run with ID matching request.')

KeyError: 'No run with ID matching request.'

[thjsal]

Could it it be this line that is the problem? Is this saving the number of modes in some model specific variable, or to something that gets overwritten every time a new model is loaded?

[thjsal]

I guess the mode_len should be run-specific, since even different runs can have different number of modes.

[thjsal]

So I think there are a few issues. First of all, the number of modes (and the other mode info too) is saved to a class variable that gets overwritten every time runs to new models are loaded. Secondly, this does not seem to work (at least in a general case) either if loading many runs for the same model, since it assumes that all the runs have the same number of modes.

I am not sure if there is an easy way to fix it, but maybe we can just raise errors/warnings if somebody tries to use this with multiple models/runs, or at least warn about this in the docstring. EDIT: I have done that now.

[lmauviard]

A solution for the 1st issue would be to make multi-modes variables (like mode_len) not class variables. This can be achieved by replacing return cls(runs, likelihood, ID, **kwargs) by something like
Runs = cls(runs, likelihood, ID, **kwargs)
Runs.mode_len =len(modes)
return Runs

For the second issue, this can be fixed by storing the mode_len and other such values in form of a dictionary, with the keys being the run ID and the values the mode_len ?

If you agree I can do it quickly

[thjsal]

Yes, those sounds like good ideas to me. Please go ahead if you like. I think you should then also change mode_label to not be a class variable either (so that different models can have different labels).

And I guess we still need to assume that if we set multi_mode=True, all the loaded runs should have been run using the MultiNest set to mmodal=True? EDIT: Or if not, we can make that an array of booleans instead of one boolean.

[thjsal]

Instead of fixing any of the issues I mentioned above, I just updated the branch by adding errors or warnings for situations where the code does not seem to work. You can remove them if you manage (or want) to solve the issues. If not, I am ok for merging this.

[lmauviard]

This should fix the issues you mentioned. Now, the parameter split_modes in cornerplot tells the plotter to plot all the modes if you give it a runID with multi-mode.
Assuming you give ['run1','run2_mode 2'] : if split_modes is True, it will plot all the run1 modes and only mode 2 of run2. If it is False, it will plot the full run1 and only mode 2 of run2 !

I will re-check this soon.

[thjsal]

Looks good to me and everything I tested seems to work now!