Background improvements and fixes (#355)

* Forcing max logl bkg factor to be the integration limits.

* Better avoidance of zero background upper limit in module generator.

Searching now non-zero background channels from both higher and lower energies. Fixes #350

* Syntax error fix in CustomBackground import

* Example modules updated.

* Raising errors if the given background support is not reasonable.

Fixes #353

* Likelihood checking simplified in the test runs.

* Allowing Nones in truth values.

Fixes #351

* Changelog and version update

* Updated `main.py` to include an `else` clause for the background support, and corrected X-PSI version

* Update `module_generator.py` to include an `else` clause for the background support

* bkg support `else` clause removed from `main.py`

* bkg support `else` clause removed from `module_generator.py`

* Changelog date changed.

---------

Co-authored-by: Devarshi Choudhury <[email protected]>

CHANGELOG.rst

@@ -34,6 +34,42 @@ and this project adheres to
 .. ^^^^^^^^^^^
 
 
+[v2.1.2] - 2024-02-05
+~~~~~~~~~~~~~~~~~~~~~
+
+Summary
+^^^^^^^
+
+* Updates and fixes were done to background marginalisation, post-processing, and module generator routines.
+
+Added
+^^^^^
+
+* Added error messages in the background marginalisation if providing a background support that does not fulfill the documented requirements (T.S.).
+
+Fixed
+^^^^^
+
+* Fixed the sometimes incorrect value of the factor (``B_for_integrand``) that should ensure that the exponent in the likelihood integrand is at some point unity within the integration domain in ``xpsi/likelihoods/default_background_marginalisation.pyx``. This was not always working when the background bounds were based both on the user-given bounds and on the default bounds leading to numerical problems in the integration and unnecessarily bad likelihood values in some example cases. Now ``B_for_integrand`` is forced to be within the integration limits (T.S.).
+
+* Fixed module imports in ``xpsi/module_generator.py`` (D.C., T.S.).
+
+* Fixed the background support upper limit zero replacements to work even when all the highest energy channels have zero background in ``xpsi/module_generator.py`` (T.S., S.V.).
+
+Changed
+^^^^^^^
+
+* Changed ``xpsi/PostProcessing/_metadata.py`` so that ``None`` can be given as truth value for a parameter, which is not wanted to be shown in a corner plot (T.S., Y.K.).
+
+Attribution
+^^^^^^^^^^^
+
+Tuomo Salmi (T.S.),
+Devarshi Choudhury (D.C.),
+Serena Vinciguerra (S.V.),
+Yves Kini (Y.K.)
+
+
 [v2.1.1] - 2023-11-10
 ~~~~~~~~~~~~~~~~~~~~~
 

docs/source/conf.py

@@ -80,7 +80,7 @@
 # The short X.Y version.
 version = '2.1'
 # The full version, including alpha/beta/rc tags.
-release = '2.1.1'
+release = '2.1.2'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

examples/examples_modeling_tutorial/TestRun_BB.py

@@ -490,8 +490,6 @@ def transform(self, p, **kwargs):
 
 # let's require that checks pass before starting to sample
 true_logl = -116504.074
-
-#print(likelihood(p)) #Need to print this if not using force_update in the following line.
 likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True)
 
 if __name__ == '__main__': # sample from the posterior

examples/examples_modeling_tutorial/TestRun_Num.py

@@ -627,22 +627,8 @@ def transform(self, p, **kwargs):
                   'verbose': True}
 
 # let's require that checks pass before starting to sample
-try:
-	true_logl = -68147.0113542
-	#print(likelihood(p))#Need to print this if not using force_update in the following line.
-	likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True)
-except:
-	print("Likelihood check did not pass. Checking if wrong atmosphere model installed.")
-	true_logl = -116504.074
-	#print(likelihood(p))#Need to print this if not using force_update in the following line.
-	try:
-		likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True)
-		print("Seems that blacbkody atmosphere extension was used instead of numerical.")
-		print("Please re-install X-PSI using numerical atmosphere extension if want to use this test run.")
-	except:
-		print("Seems that neither of the likelihood checks passed, so something must be wrong.")
-		exit()
-	exit()
+true_logl = -68147.0113542
+likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True)
 
 if __name__ == '__main__': # sample from the posterior
     xpsi.Sample.nested(likelihood, prior,**runtime_params)

examples/examples_modeling_tutorial/TestRun_NumBeam.py

@@ -594,22 +594,8 @@ def transform(self, p, **kwargs):
                   'verbose': True}
 
 # let's require that checks pass before starting to sample
-try:
-	true_logl = -6.84930649e+04 
-	#print(likelihood(p))#Need to print this if not using force_update in the following line.
-	likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True)
-except:
-	print("Likelihood check did not pass. Checking if blackbody atmosphere model installed.")
-	true_logl = -116504.074
-	#print(likelihood(p))#Need to print this if not using force_update in the following line.
-	try:
-		likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True)
-		print("Seems that blacbkody atmosphere extension was used instead of numerical.")
-		print("Please re-install X-PSI using numerical atmosphere extension including beaming parameters if want to use this test run.")
-	except:
-		print("Seems that neither of the likelihood checks passed, which can mean that the basic numerical atmosphere extension was used instead of the one with additional beaming parameters (required by this test run). Please re-install X-PSI with the latter if that is the case. Otherwise, something else must be wrong.")
-		exit()
-	exit()
+true_logl = -6.84930649e+04
+likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True)
 
 if __name__ == '__main__': # sample from the posterior
     xpsi.Sample.nested(likelihood, prior,**runtime_params)

examples/examples_module_generator/_auto_modules/CustomPrior.py

@@ -429,4 +429,4 @@ def transform(self, p, **kwargs):
         # compactness ratio M/R_eq
         p += [gravradius(ref['mass']) / ref['radius']]
 
-        return p
+        return p

examples/examples_module_generator/_auto_modules/main.py

@@ -1,4 +1,4 @@
-""" Main module for NICER PSR J0030+0451 <- X-PSI 2.1.0 CST+PDT"""
+""" Main module for NICER PSR J0030+0451 <- X-PSI 2.1.2 CST+PDT"""
 import os
 import argparse
 import re
@@ -855,7 +855,6 @@ def str_to_bool(x):
 from xpsi import HotRegions
 
 print('Rank reporting: %d' % xpsi._rank)
-
 if __name__ == '__main__':
     from CustomInstrument import CustomInstrument
     from CustomSignal import CustomSignal
@@ -1018,9 +1017,12 @@ def parse_value(value):
 
     for i in range(support.shape[0]):
         if support[i,1] == 0.0:
-            for j in range(i, support.shape[0]):
-                if support[j,1] > 0.0:
-                    support[i,1] = support[j,1]
+            for j in range(1, support.shape[0]):
+                if i+j < support.shape[0] and support[i+j,1] > 0.0:
+                    support[i,1] = support[i+j,1]
+                    break
+                elif i-j >= 0 and support[i-j,1] > 0.0:
+                    support[i,1] = support[i-j,1]
                     break
 
     support *= (XTI.data.exposure_time / args.XTI_background_exposure_time) * float(eval(args.XTI_background_scaling_factor)) # exposure ratio * scaling

setup.py

@@ -261,7 +261,7 @@ def EXTENSION(modname):
 
     setup(
         name = 'xpsi',
-        version = '2.1.1',
+        version = '2.1.2',
         author = 'The X-PSI Core Team',
         author_email = '[email protected]',
         url = 'https://github.com/xpsi-group/xpsi',

xpsi/PostProcessing/_metadata.py

@@ -55,7 +55,7 @@ def __init__(self, ID, names, bounds=None, labels=None,
         if kde_settings is not None:
             self.kde_settings = kde_settings
 
-        if truths is not None and None not in truths.values():
+        if truths is not None:
             self.truths = truths
         else:
             self._truths = dict(list(zip(self.names, [None] * len(self.names))))
@@ -210,7 +210,6 @@ def truths(self, obj):
             if len(obj) != len(self.names):
                 raise TypeError
             for key in obj:
-                float(obj[key])
                 if key not in self.names:
                     raise TypeError
         except TypeError:

xpsi/__init__.py

@@ -4,7 +4,7 @@
 
 """
 
-__version__ = "2.1.1"
+__version__ = "2.1.2"
 __author__ = "The X-PSI Core Team"
 
 try:

xpsi/likelihoods/default_background_marginalisation.pyx

@@ -298,6 +298,15 @@ def eval_marginal_likelihood(double exposure_time,
         double SCALE = exposure_time / n
         double _val
 
+    if(support.shape[0]!=counts.shape[0]):
+        raise TypeError('The number of energy channels in the background support does not match to that of the data.')
+
+    for i in range(<size_t>support.shape[0]):
+        if(support[i,1] == 0):
+            raise TypeError('Background upper limit cannot be set to 0.')
+        if(support[i,1] > 0 and support[i,1] - support[i,0] < 0):
+            raise TypeError('Background upper limit must be higher than the lower limit.')
+
     if background is not None:
         _background = background
 
@@ -531,6 +540,12 @@ def eval_marginal_likelihood(double exposure_time,
                     lower = support[i,0]
                     B_for_integrand = support[i,1]
 
+            #Ensuring that the maximum log-likelihood background for the integrand is within the integration limits
+            if B_for_integrand < lower:
+                B_for_integrand = lower
+            elif B_for_integrand > upper:
+                B_for_integrand = upper
+
             f.params = &a
 
             a.A = 0.0

xpsi/module_generator.py

@@ -1460,6 +1460,9 @@ def __call__(self, parser, namespace, values, option_string=None):
     from {0} import {1}_CustomBackground
         '''.format(args.custom_background_module,
                    _instrument)
+        )
+        module += (
+        '''
 else:
     from .{0} import {1}_CustomBackground
         '''.format(args.custom_background_module,
@@ -1650,9 +1653,12 @@ def parse_value(value):
 
     for i in range(support.shape[0]):
         if support[i,1] == 0.0:
-            for j in range(i, support.shape[0]):
-                if support[j,1] > 0.0:
-                    support[i,1] = support[j,1]
+            for j in range(1, support.shape[0]):
+                if i+j < support.shape[0] and support[i+j,1] > 0.0:
+                    support[i,1] = support[i+j,1]
+                    break
+                elif i-j >= 0 and support[i-j,1] > 0.0:
+                    support[i,1] = support[i-j,1]
                     break
 
     support *= ({0}.data.exposure_time / args.{0}_background_exposure_time) * float(eval(args.{0}_background_scaling_factor)) # exposure ratio * scaling