Postprocessing bugfixes and updates - tag v2.0.3. (#312)

thjsal · yveskini · sguillot · web-flow · commit 3e8cbd5f3e6c · 2023-07-11T19:18:17.000+02:00
* Fixing the bug with cornerplot parameter limits for X-PSI versions above 2.0 Fixes #301 * Option for setting CI precision for corner plots. * CHANGELOG update * Postprocess fix for many photospheres (#313) * Photosphere setter function added in Star.py for postprocessing. * Extra line removed. * Error and special case handling for precision settings. * Precisions given as list of ints or Nones, giving warning instead of error if not using it correctly, and fixing a bug in implementation. * Postprocessing figure adjustments (#315) * Minor ticks back to corner plots by default. * Install getdist and nestcheck from source suggested for latest updates. * Fix to combining multiple runs. (#316) * Ready for version update Co-authored-by: Yves Kini <kiniyves@gmail.com> Co-authored-by: Sebastien Guillot <sebguillot60@gmail.com> * Use NestCheck both with and without hyphens in file names. (#317) * Use NestCheck both with and without hyphens in file names. * Error handling improved and simplified. * Date update in Changelog. --------- Co-authored-by: Yves Kini <kiniyves@gmail.com> Co-authored-by: Sebastien Guillot <sebguillot60@gmail.com>
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -34,6 +34,37 @@ and this project adheres to
 .. ^^^^^^^^^^^
 
 
+[v2.0.3] - 2023-07-11
+~~~~~~~~~~~~~~~~~~~~~
+
+Summary
+^^^^^^^
+
+* Updates and bug fixes in post-processing.
+
+Fixed
+^^^^^
+
+* Fixed a bug when defining ``param_plot_lims`` in ``xpsi/PostProcessing/_corner.py`` caused by ``tight_gap_fraction`` being only defined in the customized GetDist version that is not used anymore. That parameter is now defined in X-PSI instead (T.S., Y.K., S.G.).
+
+* Fixed a bug when combining multiple runs in ``xpsi/PostProcessing/_runs.py``, which caused the combination sometimes fail since PolyChord (instead of MultiNest) default was used for the initial live point likelihoods in dead-birth files. This bug appeared after switching to use a non-customized version of NestCheck (after X-PSI version 2.0.0). Now the newest NestCheck version allows to change this value, and this change is now done within X-PSI. If trying to use an older NestCheck version, an error is raised (T.S., Y.K.).
+
+* Fixed the hyphens in the file names in ``xpsi/PostProcessing/_backends.py`` when reading MultiNest output files with the newest NestCheck version from GitHub, although trying still to read the filenames also with the older syntax to allow older NestCheck versions for other things than combining runs (T.S.).
+
+Added
+^^^^^
+* Added a keyword argument in ``xpsi/PostProcessing/_corner.py`` to allow user to define the decimal precisions for all the credible intervals printed in the figures (T.S.).
+
+* Added a photosphere setter in ``xpsi/Star.py`` which should allow producing residual and signal plots for models with multiple photosphere objects as explained in ``https://github.com/xpsi-group/xpsi/issues/304`` (Y.K, T.S.).
+
+* Added minor ticks back to corner plots in ``xpsi/PostProcessing/_corner.py``. Previously, these ticks were produced by the customized older GetDist version (T.S.).
+
+Attribution
+^^^^^^^^^^^
+
+Tuomo Salmi, Yves Kini, Sebastien Guillot
+
+
 [v2.0.2] - 2023-06-09
 ~~~~~~~~~~~~~~~~~~~~~
 
@@ -55,6 +86,7 @@ Fixed
 
 Attribution
 ^^^^^^^^^^^
+
 Tuomo Salmi
 
 
@@ -77,6 +109,7 @@ Added
 
 Attribution
 ^^^^^^^^^^^
+
 Tuomo Salmi,
 Evert Rol,
 Martin Heemskerk
diff --git a/docs/source/Post-processing.ipynb b/docs/source/Post-processing.ipynb
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -80,7 +80,7 @@
 # The short X.Y version.
 version = '2.0'
 # The full version, including alpha/beta/rc tags.
-release = '2.0.2'
+release = '2.0.3'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/docs/source/install.rst b/docs/source/install.rst
@@ -115,6 +115,13 @@ post-processing:
 
     conda install -c conda-forge getdist h5py nestcheck fgivenx
 
+.. note::
+
+    However, to get the most updated versions of getdist and nestcheck (which may be needed by
+    some of the X-PSI post-processing features), they should be installed from the source
+    (https://github.com/cmbant/getdist and https://github.com/ejhigson/nestcheck)
+    by cloning the repositories and running ``python setup.py install`` in them.
+
 In addition, some optional miscellaneous packages are:
 
 #. `jupyter <https://jupyter-notebook.readthedocs.io/en/stable/>`_ if you want to run X-PSI in a notebook.
diff --git a/setup.py b/setup.py
@@ -289,7 +289,7 @@ def EXTENSION(modname):
 
     setup(
         name = 'xpsi',
-        version = '2.0.2',
+        version = '2.0.3',
         author = 'The X-PSI Core Team',
         author_email = 'A.L.Watts@uva.nl',
         url = 'https://github.com/xpsi-group/xpsi',
diff --git a/xpsi/PostProcessing/_backends.py b/xpsi/PostProcessing/_backends.py
@@ -103,12 +103,20 @@ def __init__(self, root, base_dir, use_nestcheck, transform=None,
             except KeyError:
                 print('Root %r sampling implementation not specified... '
                       'assuming MultiNest for nestcheck...')
-                self._nc_bcknd = process_multinest_run(root,
+                try:
+                    self._nc_bcknd = process_multinest_run(root,
                                                        base_dir=base_dir)
+                except FileNotFoundError:
+                    self._nc_bcknd = process_multinest_run(root+"-",
+                                                   base_dir=base_dir)
             else:
                 if kwargs['implementation'] == 'multinest':
-                    self._nc_bcknd = process_multinest_run(root,
+                    try:
+                        self._nc_bcknd = process_multinest_run(root,
                                                            base_dir=base_dir)
+                    except FileNotFoundError:
+                        self._nc_bcknd = process_multinest_run(root+"-",
+                                                       base_dir=base_dir)
                 elif kwargs['implementation'] == 'polychord':
                     self._nc_bcknd = process_polychord_run(root,
                                                            base_dir=base_dir)
diff --git a/xpsi/PostProcessing/_corner.py b/xpsi/PostProcessing/_corner.py
@@ -148,7 +148,13 @@ def plot(self,
             distributions. If ``bootstrap and not separate_plots`` then
             the density distribution linewidth is set to zero if not
             explicitly specified with kwarg ``lw_1d``.
-            In addition, keyword arguments for avoiding unnecessary re-drawing of prior samples (``force_draw``, ``prior_samples_fnames`` and ``priors_identical``).
+            In addition, keyword arguments for avoiding unnecessary re-drawing of prior samples
+            (``force_draw``, ``prior_samples_fnames`` and ``priors_identical``).
+            Param``precisions`` (a list of integers or Nones) can be used to define the decimal number
+            precision for each credible interval plotted. In case of 2 parameters, one can do e.g.
+            precisions=[2,None] to use 2 digit decimal precision for the first parameter and use the
+            default automatic precision for the second.
+
 
         """
         self.set_subset(IDs, combine, combine_all,
@@ -445,6 +451,7 @@ def _plot_triangle(self,
         """
         #self.val_cred = []
         self.credible_interval_1d_all_show=credible_interval_1d_all_show
+        tight_gap_fraction = 0.13 # space between ticks and the edge
 
         if credible_interval_1d_all_show:
             KL_divergence=False
@@ -546,6 +553,15 @@ def _plot_triangle(self,
         else:
             plotter.legend.set_frame_on(legend_frameon)
 
+        params = self.params
+        for j in range(len(params.names)):
+                for i in range(j,len(params.names)):
+                    ax = plotter.subplots[i,j]
+                    ax.xaxis.set_minor_locator(AutoMinorLocator())
+                for i in range(j):
+                    ax = plotter.subplots[j,i]
+                    ax.yaxis.set_minor_locator(AutoMinorLocator())
+
         # add custom parameter plotting limits and updated autolocation
         with fragile(verbose(param_plot_lims,
                              'Applying bespoke parameter viewing intervals',
@@ -575,7 +591,7 @@ def _plot_triangle(self,
                 axis = plotter.subplots[-1,j].xaxis
                 xmin, xmax = axis.get_view_interval()
                 width = xmax - xmin
-                gap_wanted = width * plotter.settings.tight_gap_fraction
+                gap_wanted = width * tight_gap_fraction
                 tick = [x for x in axis.get_major_ticks() if xmin <= x.get_loc() <= xmax]
 
                 if tick[0].get_loc() - xmin < gap_wanted:
@@ -588,7 +604,7 @@ def _plot_triangle(self,
                     axis = plotter.subplots[j,0].yaxis
                     xmin, xmax = axis.get_view_interval()
                     width = xmax - xmin
-                    gap_wanted = width * plotter.settings.tight_gap_fraction
+                    gap_wanted = width * tight_gap_fraction
                     tick = [x for x in axis.get_major_ticks() if xmin <= x.get_loc() <= xmax]
 
                     if tick[0].get_loc() - xmin < gap_wanted:
@@ -640,6 +656,19 @@ def _plot_triangle(self,
 
         self.credible_intervals=OrderedDict()
 
+        if "precisions" in kwargs:
+            precisions = kwargs.get('precisions')
+            if (not isinstance(precisions, list)) or (not all((isinstance(element, int) or element==None) for element in precisions)):
+                print("Warning: Precisions need to be given as a list of integers or Nones. " +
+                "Using the automatic default precisions instead.")
+                precisions = [None]*plotter.subplots.shape[0]
+            elif len(precisions) != plotter.subplots.shape[0]:
+                print("Warning: Precisions list has wrong number of dimensions. " +
+                "Using the automatic default precisions instead.")
+                precisions = [None]*plotter.subplots.shape[0]
+        else:
+            precisions = [None]*plotter.subplots.shape[0]
+
         if credible_interval_1d_all_show and self.all_same(self.get_attr("parent_ID")):
             for r in range(len(self.subset_to_plot)):
 
@@ -656,7 +685,8 @@ def _plot_triangle(self,
                                             annotate_xy=annotate_xy,
                                             sixtyeight=sixtyeight,
                                             ninety=ninety,
-                                            compute_all_intervals=compute_all_intervals)
+                                            compute_all_intervals=compute_all_intervals,
+                                            precisions=precisions)
 
                 self.credible_intervals[id]=self.val_cred
         else:
@@ -673,7 +703,8 @@ def _plot_triangle(self,
                                             annotate_xy=annotate_xy,
                                             sixtyeight=sixtyeight,
                                             ninety=ninety,
-                                            compute_all_intervals=compute_all_intervals)
+                                            compute_all_intervals=compute_all_intervals,
+                                            precisions=precisions)
 
                 self.credible_intervals[id]=self.val_cred
 
@@ -846,7 +877,7 @@ def KL(ns_run, logw):
                   'Added 1D marginal credible intervals')
     def _add_credible_interval(self, plotter, posterior, bootstrap, n_simulate,
                                annotate, annotate_xy, sixtyeight,
-                               ninety, compute_all_intervals):
+                               ninety, compute_all_intervals, precisions=None):
         """
         Estimate 1-:math:`\sigma` credible interval in one-dimension on a
         combined run, or if such a run does not exist, on the run with
@@ -863,6 +894,10 @@ def _add_credible_interval(self, plotter, posterior, bootstrap, n_simulate,
             ``False`` plots 95% credible interval about the median -- i.e.,
             symmetric quantiles about the median.
         """
+
+        if precisions is None:
+            precisions = [None]*plotter.subplots.shape[0]
+
         diag = [plotter.subplots[i,i] for i in range(plotter.subplots.shape[0])]
 
         run = self.run #posterior.subset_to_plot[0]
@@ -882,7 +917,7 @@ def estimator(*args, **kwargs):
 
         quantiles = [0.159, 0.5, 0.841] if sixtyeight else ([0.05,0.5,0.95] if ninety else [0.025, 0.5, 0.975])
 
-        def format_CI(name, cred, summary, additional=2, sscript=False):
+        def format_CI(name, cred, summary, additional=2, sscript=False, precision=None):
 
             if len(cred.shape) > 1:
                 _qs = (cred[1,1],
@@ -893,9 +928,12 @@ def format_CI(name, cred, summary, additional=2, sscript=False):
                 _qs = (cred[1],
                        cred[1] - cred[0],
                        cred[2] - cred[1])
-
             _p = max(_precision(_qs[0]), _precision(_qs[1]), _precision(_qs[2]))
-            _f = '%.' + str(_p + additional) + 'f'
+
+            if precision is None:
+                _f = '%.' + str(_p + additional) + 'f'
+            else:
+                _f = '%.' + str(precision) + 'f'
 
             if name: name += ' '
 
@@ -955,7 +993,8 @@ def calculate_intervals(quantiles):
                                       cred,
                                       68 if sixtyeight else (90 if ninety else 95),
                                       additional=1,
-                                      sscript=True)
+                                      sscript=True,
+                                      precision=precisions[i])
 
                     title = ax.get_title()
                     if title:
@@ -1004,6 +1043,7 @@ def calculate_intervals(quantiles):
 
 
                 if compute_all_intervals:
+
                     yield format_CI(self.params.names[i],
                                     cred,
                                     68 if sixtyeight else (90 if ninety else 95))
@@ -1061,7 +1101,8 @@ def calculate_intervals(quantiles):
                                       cred,
                                       68 if sixtyeight else (90 if ninety else 95),
                                       additional=1,
-                                      sscript=True)
+                                      sscript=True,
+                                      precision=precisions[i])
 
                     title = ax.get_title()
                     if title:
diff --git a/xpsi/PostProcessing/_runs.py b/xpsi/PostProcessing/_runs.py
@@ -185,12 +185,24 @@ def _combine(self, combine_all, overwrite):
             run['output']['base_dir'] = base_dir
             run['output']['file_root'] = file_root
 
-            write_run_output(run,
-                             write_dead = True,
-                             write_stats = True,
-                             posteriors = True,
-                             stats_means_errs = True,
-                             n_simulate = 1000)
+            try:
+                #use MultiNest initial likelihood (logl_init) instead of the default PolyChord
+                write_run_output(run,
+                                 write_dead = True,
+                                 write_stats = True,
+                                 posteriors = True,
+                                 stats_means_errs = True,
+                                 n_simulate = 1000,
+                                 logl_init = -0.179769313486231571E+309)
+            except TypeError as e:
+                if str(e) == "Unexpected **kwargs: {'logl_init': -1.7976931348623157e+308}":
+                    raise TypeError("The used nestcheck version does not support combining "
+                    "MultiNest runs in X-PSI. To use this feature, nestcheck version newer "
+                    "than in this commit: "
+                    "https://github.com/ejhigson/nestcheck/commit/513ef962ef7b0d66377686f9fe0a9e354dad48b3 "
+                    "should be used (see installation instructions for installing it from github).")
+                else:
+                    raise
 
         kwargs = {'kde_settings': self.kde_settings,
                   'ID': 'combined',
diff --git a/xpsi/Star.py b/xpsi/Star.py
@@ -62,6 +62,10 @@ def photospheres(self):
         """ Get the list of photosphere objects. """
         return self._photospheres
 
+    @photospheres.setter
+    def photospheres(self, obj):
+        self._photospheres = [obj]
+
     def activate_fast_mode(self, activate):
         try:
             for photosphere in self._photospheres:
diff --git a/xpsi/__init__.py b/xpsi/__init__.py
@@ -4,7 +4,7 @@
 
 """
 
-__version__ = "2.0.2"
+__version__ = "2.0.3"
 __author__ = "The X-PSI Core Team"
 
 try:
