From 1ff24717eb34cc777ee77ec04b53a1427e55d4b4 Mon Sep 17 00:00:00 2001 From: Tuomo Salmi Date: Wed, 11 Dec 2024 16:14:24 +0200 Subject: [PATCH] Example scripts fixed (#493) * TestRun_BB.py matched better with the new version of the corresponding modeling tutorial. * Likelihood values and intro text updated in other examples. --- .../examples_modeling_tutorial/TestRun_BB.py | 220 ++++++++---------- .../examples_modeling_tutorial/TestRun_Num.py | 4 +- .../TestRun_NumBeam.py | 5 +- 3 files changed, 101 insertions(+), 128 deletions(-) diff --git a/examples/examples_modeling_tutorial/TestRun_BB.py b/examples/examples_modeling_tutorial/TestRun_BB.py index aa5b3be8..0164c3e9 100644 --- a/examples/examples_modeling_tutorial/TestRun_BB.py +++ b/examples/examples_modeling_tutorial/TestRun_BB.py @@ -2,11 +2,6 @@ Test script to check that X-PSI installation is working (with the default blackbody atmosphere extension). The run is succesful if it passes the likelihood check (in ~ 1 minute) and sampling finishes without errors (in ~ 5 min). The model and synthetic data are based on those shown in the Modelling tutorial notebook. - -Prequisities: -Before running the script, add the NICER instrument files to the model_data subdirectory: -nicer_v1.01_arf.txt, nicer_v1.01_rmf_energymap.txt, and nicer_v1.01_rmf_matrix.txt -(found from https://doi.org/10.5281/zenodo.7094144). ''' @@ -28,6 +23,7 @@ phases=np.linspace(0.0, 1.0, 33), first=0, last=180, exposure_time=984307.6661) + data = xpsi.Data(**settings) class CustomInstrument(xpsi.Instrument): @@ -47,61 +43,88 @@ def __call__(self, signal, *args): return self._folded_signal - @classmethod - def from_response_files(cls, ARF, RMF, channel_edges, max_input, - min_input=0, - ): - """ Constructor which converts response files into :class:`numpy.ndarray`s. - :param str ARF: Path to ARF which is compatible with - :func:`numpy.loadtxt`. - :param str RMF: Path to RMF which is compatible with - :func:`numpy.loadtxt`. - :param str channel_edges: Path to edges which is compatible with - :func:`numpy.loadtxt`. - """ - link ="https://doi.org/10.5281/zenodo.7094144" - if min_input != 0: - min_input = int(min_input) +NICER = CustomInstrument.from_ogip_fits(RMF_path='nicer_20170601v003.rmf', + ARF_path='nicer_20170601v005.arf', + datafolder='../../examples/examples_modeling_tutorial/model_data') - max_input = int(max_input) - try: - ARF = np.loadtxt(ARF, dtype=np.double, skiprows=3) - except: - print("ERROR: You miss the following file: {}\nThe file is found from here: {}".format(ARF, link)) - exit() - try: - RMF = np.loadtxt(RMF, dtype=np.double) - except: - print("ERROR: You miss the following file: {}\nThe file is found from here: {}".format(RMF, link)) - exit() +from xpsi.likelihoods.default_background_marginalisation import eval_marginal_likelihood +from xpsi.likelihoods.default_background_marginalisation import precomputation - try: - channel_edges = np.loadtxt(channel_edges, dtype=np.double, skiprows=3)[:,1:] - except: - print("ERROR: You miss the following file: {}\nThe file is found from here: {}".format(channel_edges, link)) - exit() +class CustomSignal(xpsi.Signal): + """ A custom calculation of the logarithm of the likelihood. + + We extend the :class:`xpsi.Signal.Signal` class to make it callable. + + We overwrite the body of the __call__ method. The docstring for the + abstract method is copied. + + """ + + def __init__(self, workspace_intervals = 1000, epsabs = 0, epsrel = 1.0e-8, + epsilon = 1.0e-3, sigmas = 10.0, support = None, *args, **kwargs): + """ Perform precomputation. """ - matrix = np.ascontiguousarray(RMF[min_input:max_input,20:201].T, dtype=np.double) + super(CustomSignal, self).__init__(*args, **kwargs) - edges = np.zeros(ARF[min_input:max_input,3].shape[0]+1, dtype=np.double) + try: + self._precomp = precomputation(self._data.counts.astype(np.int32)) + except AttributeError: + print('No data... can synthesise data but cannot evaluate a ' + 'likelihood function.') + else: + self._workspace_intervals = workspace_intervals + self._epsabs = epsabs + self._epsrel = epsrel + self._epsilon = epsilon + self._sigmas = sigmas - edges[0] = ARF[min_input,1]; edges[1:] = ARF[min_input:max_input,2] + if support is not None: + self._support = support + else: + self._support = -1.0 * np.ones((self._data.counts.shape[0],2)) + self._support[:,0] = 0.0 - for i in range(matrix.shape[0]): - matrix[i,:] *= ARF[min_input:max_input,3] + @property + def support(self): + return self._support - channels = np.arange(20, 201) + @support.setter + def support(self, obj): + self._support = obj - return cls(matrix, edges, channels, channel_edges[20:202,-2]) + def __call__(self, *args, **kwargs): + self.loglikelihood, self.expected_counts, self.background_signal, self.background_signal_given_support = \ + eval_marginal_likelihood(self._data.exposure_time, + self._data.phases, + self._data.counts, + self._signals, + self._phases, + self._shifts, + self._precomp, + self._support, + self._workspace_intervals, + self._epsabs, + self._epsrel, + self._epsilon, + self._sigmas, + kwargs.get('llzero')) -NICER = CustomInstrument.from_response_files(ARF = 'model_data/nicer_v1.01_arf.txt', - RMF = 'model_data/nicer_v1.01_rmf_matrix.txt', - channel_edges = 'model_data/nicer_v1.01_rmf_energymap.txt', - max_input = 500, - min_input = 0) +NICER.trim_response(min_channel=20, max_channel=200, threshold=2.) +signal = CustomSignal(data = data, + instrument = NICER, + background = None, + interstellar = None, + workspace_intervals = 1000, + cache = True, + epsrel = 1.0e-8, + epsilon = 1.0e-3, + sigmas = 10.0, + min_channel = 20, + max_channel = 200, + support = None) bounds = dict(distance = (0.1, 1.0), # (Earth) distance mass = (1.0, 3.0), # mass @@ -127,7 +150,6 @@ def from_response_files(cls, ARF, RMF, channel_edges, max_input, num_leaves=100, num_rays=200, atm_ext="BB", - image_order_limit=3, prefix='p') class derive(xpsi.Derive): @@ -164,7 +186,6 @@ def __call__(self, boundto, caller = None): num_rays=200, do_fast=False, atm_ext="BB", - image_order_limit=3, is_antiphased=True, prefix='s') @@ -231,78 +252,6 @@ def global_variables(self): print("Bolometric pulse for 2nd spot:") print(repr(np.sum(photosphere.signal[1][0], axis=0))) -from xpsi.likelihoods.default_background_marginalisation import eval_marginal_likelihood -from xpsi.likelihoods.default_background_marginalisation import precomputation - -class CustomSignal(xpsi.Signal): - """ A custom calculation of the logarithm of the likelihood. - - We extend the :class:`xpsi.Signal.Signal` class to make it callable. - - We overwrite the body of the __call__ method. The docstring for the - abstract method is copied. - - """ - - def __init__(self, workspace_intervals = 1000, epsabs = 0, epsrel = 1.0e-8, - epsilon = 1.0e-3, sigmas = 10.0, support = None, *args, **kwargs): - """ Perform precomputation. """ - - super(CustomSignal, self).__init__(*args, **kwargs) - - try: - self._precomp = precomputation(self._data.counts.astype(np.int32)) - except AttributeError: - print('No data... can synthesise data but cannot evaluate a ' - 'likelihood function.') - else: - self._workspace_intervals = workspace_intervals - self._epsabs = epsabs - self._epsrel = epsrel - self._epsilon = epsilon - self._sigmas = sigmas - - if support is not None: - self._support = support - else: - self._support = -1.0 * np.ones((self._data.counts.shape[0],2)) - self._support[:,0] = 0.0 - - @property - def support(self): - return self._support - - @support.setter - def support(self, obj): - self._support = obj - - def __call__(self, *args, **kwargs): - self.loglikelihood, self.expected_counts, self.background_signal, self.background_signal_given_support = \ - eval_marginal_likelihood(self._data.exposure_time, - self._data.phases, - self._data.counts, - self._signals, - self._phases, - self._shifts, - self._precomp, - self._support, - self._workspace_intervals, - self._epsabs, - self._epsrel, - self._epsilon, - self._sigmas, - kwargs.get('llzero')) - -signal = CustomSignal(data = data, - instrument = NICER, - background = None, - interstellar = None, - workspace_intervals = 1000, - cache = True, - epsrel = 1.0e-8, - epsilon = 1.0e-3, - sigmas = 10.0, - support = None) from scipy.stats import truncnorm class CustomPrior(xpsi.Prior): @@ -461,7 +410,7 @@ def transform(self, p, **kwargs): num_energies=128, threads=1, prior=prior, - externally_updated=True) + externally_updated=False) wrapped_params = [0]*len(likelihood) wrapped_params[likelihood.index('p__phase_shift')] = 1 @@ -488,8 +437,31 @@ def transform(self, p, **kwargs): 'max_iter': 100,# manual termination condition for short test 'verbose': True} +xpsi.set_phase_interpolant('Akima') +p = [1.4, + 12.5, + 0.2, + math.cos(1.25), + 0.0, + 1.0, + 0.075, + 6.2, + 0.025, + math.pi - 1.0, + 0.2] + +for h in hot.objects: + h.set_phases(num_leaves = 100) + +likelihood.threads = 3 +likelihood.reinitialise() +likelihood.clear_cache() + +# inform source code that parameter objects updated when inverse sampling +likelihood.externally_updated = True + # let's require that checks pass before starting to sample -true_logl = -116504.074 +true_logl = -26947.474150 likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True) if __name__ == '__main__': # sample from the posterior diff --git a/examples/examples_modeling_tutorial/TestRun_Num.py b/examples/examples_modeling_tutorial/TestRun_Num.py index c4a8eb70..336f4bc9 100644 --- a/examples/examples_modeling_tutorial/TestRun_Num.py +++ b/examples/examples_modeling_tutorial/TestRun_Num.py @@ -1,7 +1,7 @@ ''' Test script to check that X-PSI installation is working (with the numerical atmosphere extension). The run is succesful if it passes the likelihood check (in ~ 1 minute) and sampling finishes without errors (in ~ 5 min). -The model and synthetic data are based on those shown in the Modelling tutorial notebook. +The model and instrument response are not exactly the same that were used to produce the synthetic data, but this is just an example. Prequisities: Before running the script, add the NICER instrument files to the model_data subdirectory: @@ -620,7 +620,7 @@ def transform(self, p, **kwargs): 'verbose': True} # let's require that checks pass before starting to sample -true_logl = -68147.0113542 +true_logl = -6.6981424550e+04 likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True) if __name__ == '__main__': # sample from the posterior diff --git a/examples/examples_modeling_tutorial/TestRun_NumBeam.py b/examples/examples_modeling_tutorial/TestRun_NumBeam.py index db27714c..3d5add59 100644 --- a/examples/examples_modeling_tutorial/TestRun_NumBeam.py +++ b/examples/examples_modeling_tutorial/TestRun_NumBeam.py @@ -1,7 +1,8 @@ ''' Test script to check that X-PSI installation is working (with the numerical atmosphere extension with free beaming options). The run is succesful if it passes the likelihood check (in ~ 1 minute) and sampling finishes without errors (in ~ 5 min). -The model and synthetic data are based on those shown in the Modelling tutorial notebook, although adding here new parameters on the atmospheric beaming. +The model and instrument response are not exactly the same that were used to produce the synthetic data, but this is just an example. +The model here includes additional free parameters on the atmospheric beaming. Prequisities: Before running the script, add the NICER instrument files to the model_data subdirectory: @@ -589,7 +590,7 @@ def transform(self, p, **kwargs): 'verbose': True} # let's require that checks pass before starting to sample -true_logl = -6.84930649e+04 +true_logl = -6.7210627096e+04 likelihood.check(None, [true_logl], 1.0e-6,physical_points=[p],force_update=True) if __name__ == '__main__': # sample from the posterior