johannesbuchner / bxa Goto Github PK

In the documentation several important features are listed but without explanation. It would be helpful if these were explained quantitatively with examples / comparisons:

when systematically analysing a large data-set

How does this differ from what is available in XSPEC/Sherpa?

when comparing multiple models

This is explained well for those initiated with Bayes factor model comparison, but a link to that part of the docs would be useful to make the point

when analysing low counts data-set

This claim is made here and here as well as the normal use of proper Poisson likelihoods. What specific features makes BXA shine in comparison? Can links be provided as well as comparative examples?

when you don’t want to babysit your fits

What is meant by this statement? As stated above, can examples be linked/ shown that validate this?
I'm assuming it is because there is a stopping criteria in nested sampling, but is this meant to inform the user that they do not need to access the correctness of their posterior?

when you don’t want to test MCMC chains for their convergence

I suppose this could be answered with the previous question.

In general, it could be important to fully separate what is a feature or strength of BXA vs what is inherited from the package it relies on / compliments.

If I load the chain built with bxa into xspec, it is not used to compute the errors or the flux uncertainties.
I tried with the chain produced by xspec and it works as expected.

It must be something in the format that makes the difference, but finding exactly what is challenging.

Hi Johannes,

I am trying to use BXA to compare two models on an XMM-Newton data set and I have issues of convergence, probably due to some parameter degeneracy, even if models are not very complex:

constant * TBabs * pow (This one converged in an hour or so)
constant * TBabs * pcfabs * pow

This last one has been running for more than one day and it does not look to be converging.
The warning is:
/pyenv/versions/3.8.2/lib/python3.8/site-packages/ultranest/integrator.py:1632: UserWarning: Sampling from region seems inefficient (0/40 accepted in iteration 2500). To improve efficiency, modify the transformation so that the current live points (stored for you in BXA-pcfabs/extra/sampling-stuck-it%d.csv) are ellipsoidal, or use a stepsampler, or set frac_remain to a lower number (e.g., 0.5) to terminate earlier.
warnings.warn(warning_message)

I have already put frac_reami=0.5 in my parameters, but I do not understand well this thing of ellipsoids...

Do you have some suggestions to increase efficiency ?

You find attached some files from the run

debug.log
sampling-stuck-it%d.csv

The WIP notebook is in
https://gitlab.astro.unige.ch/xmm-newton-workflows/0862410301
which uses the call to BXA in
https://gitlab.astro.unige.ch/ferrigno/pysas/-/blob/master/pyxmmsas/__init__.py#L420

• BXA version: 4.1.1
• UltraNest version: 3.5.7
• Python version: 3.8.17
• Xspec or Sherpa and version: 12.13.0c (Xspec),
• Operating System: MacOS Ventura 13.5

Description

I tried to install BXA using this guide.

What I Did

During initial tests the following error appeared:

python -c 'import xspec'
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/xspec/__init__.py", line 71, in <module>
    from mxspec import callModelFunction
  File "/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/__init__.py", line 13, in <module>
    from .mmodel import callModelFunction
  File "/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/mmodel.py", line 4, in <module>
    from . import _pymXspec
ImportError: dlopen(/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/_pymXspec.so, 0x0002): tried: '/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/_pymXspec.so' (no such file), '/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/_pymXspec.so' (mach-o file, but is an incompatible architecture (have 'arm64', need 'x86_64')), '/System/Volumes/Preboot/Cryptexes/OS/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/_pymXspec.so' (no such file), '/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/_pymXspec.so' (mach-o file, but is an incompatible architecture (have 'arm64', need 'x86_64')), '/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/lib_pymXspec.dylib' (no such file), '/Users/mike/Soft/heasoft-6.31.1/Xspec/aarch64-apple-darwin22.5.0/lib/python/mxspec/lib_pymXspec.dylib' (mach-o file, but is an incompatible architecture (have 'arm64', need 'x86_64')), '/System/Volumes/Preboot/Cryptexes/OS/Users/mike/Soft/heasoft-6.31.1/Xspec/aarch64-apple-darwin22.5.0/lib/python/mxspec/lib_pymXspec.dylib' (no such file), '/Users/mike/Soft/heasoft-6.31.1/Xspec/aarch64-apple-darwin22.5.0/lib/python/mxspec/lib_pymXspec.dylib' (mach-o file, but is an incompatible architecture (have 'arm64', need 'x86_64'))

python -c 'import bxa.xspec as bxa'
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/Users/mike/opt/anaconda3/envs/xspecbxa/lib/python3.8/site-packages/bxa/xspec/__init__.py", line 12, in <module>
    from . import qq
  File "/Users/mike/opt/anaconda3/envs/xspecbxa/lib/python3.8/site-packages/bxa/xspec/qq.py", line 11, in <module>
    from xspec import Plot
  File "/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/xspec/__init__.py", line 71, in <module>
    from mxspec import callModelFunction
  File "/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/__init__.py", line 13, in <module>
    from .mmodel import callModelFunction
  File "/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/mmodel.py", line 4, in <module>
    from . import _pymXspec
ImportError: dlopen(/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/_pymXspec.so, 0x0002): tried: '/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/_pymXspec.so' (no such file), '/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/_pymXspec.so' (mach-o file, but is an incompatible architecture (have 'arm64', need 'x86_64')), '/System/Volumes/Preboot/Cryptexes/OS/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/_pymXspec.so' (no such file), '/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/python/mxspec/_pymXspec.so' (mach-o file, but is an incompatible architecture (have 'arm64', need 'x86_64')), '/Users/mike/Soft/heasoft-6.31.1/aarch64-apple-darwin22.5.0/lib/lib_pymXspec.dylib' (no such file), '/Users/mike/Soft/heasoft-6.31.1/Xspec/aarch64-apple-darwin22.5.0/lib/python/mxspec/lib_pymXspec.dylib' (mach-o file, but is an incompatible architecture (have 'arm64', need 'x86_64')), '/System/Volumes/Preboot/Cryptexes/OS/Users/mike/Soft/heasoft-6.31.1/Xspec/aarch64-apple-darwin22.5.0/lib/python/mxspec/lib_pymXspec.dylib' (no such file), '/Users/mike/Soft/heasoft-6.31.1/Xspec/aarch64-apple-darwin22.5.0/lib/python/mxspec/lib_pymXspec.dylib' (mach-o file, but is an incompatible architecture (have 'arm64', need 'x86_64'))

I can run Xspec and PyXspec without any problems.

• BXA version: 2.9
• Python version: 3.8.12
• Xspec version: 12.12.1
• Operating System: macOS Monterey v12.6

Description

I am fitting 9 spectra simultaneously with PyXspec (one XMM/PN, one XMM/MOS1, one XMM/MOS2 for three different epochs). Since the spectra are all low signal to noise (source is <~50% of total counts), I am fitting a simple model but have a cross-calibration constant for each dataset that are free to vary relative to the first one (i.e. 8 free parameters in addition to the free parameters of the main model). For these cross-calibration constants I am using a custom log-Gaussian prior (code attached) in the BXA fit.

After fitting, the posteriors for a lot of the constants and some of the model parameters are all a single value which then gives an error when trying to produce the corner plot with corner after the fit since the parameter posteriors have no dynamic range.

Are there some suggested practices to overcome this? E.g., would increasing the number of live points help for the fit? This issue does not occur when I only fit 3 spectra simultaneously (i.e. with two cross-calibration constants), so I am wondering if the custom log-Gaussian prior may be too restrictive for fitting with a lot of spectra?

Code

Here is the code I am using for the custom log-Gaussian prior

def create_loggaussian_prior_for(model, par, mean, std):
	"""
	This combines the Gaussian and loguniform BXA priors to give
	a Gaussian distribution of the log-parameter.

	NOTE: mean and std are also in log units
	"""
	import scipy.stats
	pval, pdelta, pmin, pbottom, ptop, pmax = par.values

	if pmin == 0:
		raise Exception('You forgot to set reasonable parameter limits on %s' % par.name)
	low = np.log10(pmin)
	spread = np.log10(pmax) - np.log10(pmin)
	hi = np.log10(pmax)
	if spread > 10:
		print('   note: this parameter spans *many* dex. Double-check the limits are reasonable.')

	print('  log-gaussian prior for %s of %f +- %f' % (par.name, mean, std))
	rv = scipy.stats.norm(mean, std)
	def loggauss_transform(x): 
		return max(low, min(hi, rv.ppf(x * spread + low)))
	def loggauss_after_transform(x): return 10**x
	return dict(model=model, index=par._Parameter__index, name='log(%s)' % par.name, 
		transform=loggauss_transform, aftertransform=loggauss_after_transform)

After using the auto_background function in Sherpa, the data and fit plots use counts in the y axis. Since the default behavior of Sherpa is using count rates, this change after calculating PCA models can be quite annoying. It is something easy to fix by the user, but it can take some time to realize what happen.

As far as I understand this is caused by the following line in the fit method of the PCAFitter class:

It would be helpful if a note is included in the documentation so the users are aware of that change when using PCA models.

Ideally, the best would be to recover the initial behavior of Sherpa after running a PCA fit. I can give it a try, if you think is something worth it.

I use BXA in pyxspec. In some cases, there are nan values in lowmodel_parts. So the assert error is reported in line 85 of gof.py.
I find some values of kplusones equal 0 (line 78, 79).

Hi,

there is a Python 3 issue in BXA/bxa/sherpa/background/fitters.py, where the function file() is still used.

Hi Johannes,

Thank you for the wonderful package! This is probably a more conceptual question than a github "issue":

I was wondering if you have any suggestions for dealing with degeneracy in identical parameters in multi-component models. E.g., fitting with a model comprised of two apec components, and aiming to penalize when temperatures in both models are very similar (e.g., it can't find two temperatures, so the best-fit just splits the amplitude between the two components with similar temperature). I was trying to implement some trick, e.g., using the create_custom_prior_for(), but I couldn't come up with anything viable. Right now, I just let both parameters run wild and if the posteriors for both are really consistent and there are no other interesting features in any of the posteriors, I just conclude that the second component is unnecessary (which is probably ok, but not the best thing to do). I think without such a penalization, model evidence LogZ might naturally prefer the two-component model. Am I correct?

Any suggestion would be appreciated.
Arash

P.S. My apologies if github issue is not a proper venue for such a question.

I am running BXA v2.10 on macOS Mojave v10.14.6, and am trying to create a new PCA model for the Athena/X-IFU from the publicly available background file (available from http://x-ifu.irap.omp.eu/resources/for-the-community) with the command:
python autobackgroundmodel/__init__.py bkg1.pha bkg2.pha bkg3.pha

First, the spectrum did not have a COUNTS column, so I created one with astropy based on the pre-existing RATE column and EXPOSURE header keywords. Since the spectrum appears to contain many counts, I wondered if it would be sufficient by itself to create a PCA model, instead of requiring multiple individual observations stacked together.

I then ran the following command:
python3 /astrosoftware/BXA/autobackgroundmodel/__init__.py Total_pointsources_XIFU_CC_BASELINECONF_2018_10_10_COUNTS.pha

But I get the error message:

fetching counts from 1 files ...
  reading "Total_pointsources_XIFU_CC_BASELINECONF_2018_10_10_COUNTS.pha"
combining ... 10
Traceback (most recent call last):
  File "/astrosoftware/BXA/autobackgroundmodel/__init__.py", line 131, in <module>
    assert len(telescopes) == 1, "trying to combine multiple telescopes!"
AssertionError: trying to combine multiple telescopes!

Is it possible this error is arising because I am using a much older version of BXA than the current one? If not, is there a way I can generate the PCA background for the X-IFU?

I just installed bxa (pip install on Mac OS High Sierra), but I get an error when I'm trying to import bxa.sherpa to my ipython:

$ ipython
Python 3.7.2 (default, Dec 29 2018, 00:00:04) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.2.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: from __future__ import print_function                                                                           

In [2]: import bxa.sherpa as bxa                                                                                        
Traceback (most recent call last):

  File "/Users/aneta/anaconda3/envs/test_py3.7/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 3267, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)

  File "<ipython-input-2-19cf673f1dd0>", line 1, in <module>
    import bxa.sherpa as bxa

  File "/Users/aneta/anaconda3/envs/test_py3.7/lib/python3.7/site-packages/bxa/sherpa/__init__.py", line 63
    print 'Exception in log_likelihood function: ', e
                                                 ^
SyntaxError: Missing parentheses in call to 'print'. Did you mean print('Exception in log_likelihood function: ', e)?

• BXA version: 4.0.6
• UltraNest version: 3.5.5
• Python version:3.8.10
• Xspec or Sherpa and version: XSPEC version: 12.12.0
• Operating System: Linux

Description

When passing the parameter evidence_tolerance the run does not work

What I Did

line 194 of solver.py put
run_kwargs['dlogz'] = run_kwargs.pop('dlogz', evidence_tolerance)
instead of
run_kwargs['evidence_tolerance'] = run_kwargs.pop('dlogz', evidence_tolerance)

It fixes easily the issue.

Is it possible there is a discrepancy between the version of the sinning.py module installed via pip and the version available in the source code on GitHub?

Specifically, at line 106, the pip-installed version uses:
pxlo = xlo[(nstats[:,1] * n).astype(numpy.int)]

while the source code uses:
pxlo = xlo[(nstats[:,1] * n).astype(int)]

I wonder if this correction has not yet been incorporated into the version distributed through pip? I installed BXA using pip 24.0.

• BXA version:

# Name                    Version                   Build  Channel
bxa                       4.0.5              pyh44b312d_0    conda-forge
xorg-libxau               1.0.9                h7f98852_0    conda-forge

• UltraNest version:

# Name                    Version                   Build  Channel
ultranest                 3.5.7           py310hde88566_0    conda-forge

• Python version:

$ python --version
   Python 3.10.9

• Xspec or Sherpa and version:

xspec-modelsonly          12.12.1c             h707b37c_0    https://cxc.cfa.harvard.edu/conda/ciao
Sherpa 4.15.0

• Operating System:

$ uname -a
Linux DESKTOP-10UCA6U 5.15.90.1-microsoft-standard-WSL2 #1 SMP Fri Jan 27 02:56:13 UTC 2023 x86_64 x86_64 x86_64 GNU/Linux

Description

I am in theciao-4.15environment, using conda install -c conda-forge XPA shows that the installation has been successful, when I use python -c "import xspec" the following error occurs

$ python -c "import xspec"
terminate called after throwing an instance of 'XSModelFunction::NoSuchComponent'
Aborted

Can someone help me out with this problem? Thanks!

Hi,

Is there an option to save the model file similar to the one in Xspec? It becomes difficult to track the exact models/configurations while doing many fits. I currently solve this issue by simply adding Xset.save(filename, 'm') to the programs, but it would be good to have that option within BXA if possible.

Thanks!

Hi, a few comments/questions on the documentation as part of the review:

• Installation: Add the full list of dependencies to the install steps. I also needed cython, tqdm, h5py and ultranest as a minimal list to get things running.
• For both the sherpa/Xspec cases, it would be nice to see a full standalone demo including some example data. I know this is done in the scripts in examples/ , but would be helpful if at least one case was treated notebook-style on the documentation pages.
• Similarly for the xagnfitter.py script - you have a case with some example data in the Travis CI, it could be nice to also have this in the docs as a walkthrough for users.
• The default ultranest output is cool, but could be distracting/annoying when running in certain settings. Could you show in the docs an easy interface in BXA to switch it off?
• In the "convenience features" documentation
  • Seems the example call of gal.py is a typo as refers to another script
  • Would be nice to have some quick code showing exactly how to import chains in Xspec or run a script in parallel

Since Ultranest is online now, will the next update of BXA also using it instead of Multinest?
Under the hood both packages are doing the same, I assume (or are there some performance improvements as well?).
However, with Ultranest the installation of BXA should be much easier, right?

PS: nice package naming ;)

• BXA version: 4.0.5
• UltraNest version: 3.3.3
• Python version: 3.8.12
• Xspec version: 12.12.1
• Operating System: macOS Monterey 12.6

Description

When using the BXA-generated chain files in Xspec after running a BXA fit, Xspec cannot estimate uncertainties from the chain with e.g., "flux 2. 10. err 1000". The issue seems to be that the BXA chain file column headings are incompatible with the model used.

What I Did

The model fit was:

Model constant<1>*TBabs<2>(zTBabs<3>*cabs<4>*zpowerlw<5> + zgauss<6>) Source No.: 1   Active/On
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   constant   factor              1.00000      frozen
   2    2   TBabs      nH         10^22    4.00000E-02  frozen
   3    3   zTBabs     nH         10^22    37.1630      +/-  0.0          
   4    3   zTBabs     Redshift            5.00000E-02  frozen
   5    4   cabs       nH         10^22    37.1630      = p3
   6    5   zpowerlw   PhoIndex            1.88462      +/-  0.0          
   7    5   zpowerlw   Redshift            5.00000E-02  frozen
   8    5   zpowerlw   norm                1.37403E-03  +/-  0.0          
   9    6   zgauss     LineE      keV      6.41373      +/-  0.0          
  10    6   zgauss     Sigma      keV      8.41334E-02  +/-  0.0          
  11    6   zgauss     Redshift            5.00000E-02  frozen
  12    6   zgauss     norm                1.58220E-05  +/-  0.0          

And the free parameter numbers in the BXA fit were: 3, 6, 8, 9, 10, 12. The corresponding chain.fits file generated by BXA had column headings (note the nH prior was added out of order of the other parameters):

PhoIndex__6, log(norm)__20, log(nH)__27, LineE__45, log(Sigma)__58, log(norm)__72, FIT_STATISTIC

The problem was fixed by manually changing the BXA chain file column headings, column orders and parameter units to match an equivalent Xspec-generated MCMC chain file:

Columns = nH__3, PhoIndex__6, norm__8, LineE__9, Sigma__10, norm__12, FIT_STATISTIC
Units = 10^22, , , keV, keV, , C-Statistic

Could you add "_{ndata}" to the name of the output file in autobackgroundmodel/init.py, since it is required by autobackgroundmodel/fitbkg.py?

• BXA version: 4.1.1
• UltraNest version: 3.5.7
• Python version: 3.11.3
• Xspec or Sherpa and version: xspec 12.13.0c
• Operating System: Mac OS Ventura 13.2.1 (also on Ubuntu 22)

Description

Running solver.run() cannot be interrupted in jupyter lab.

What I Did

Running this in a jupyter lab cell:

solver = bxa.BXASolver(transformations=transformations, outputfiles_basename=outputfiles_basename)
results = solver.run(resume=True)

Trying to interrupt the process/kernel does not work, the iterations go on and on. The only remedy is to restart the kernel which leads to losing the session.

I am not sure if this is BXA or ultranest issue, there should be a possibility to make the iterations interruptible somehow.

Thanks,
Ivan

More than an issue, it is a request for clarification: should the Y-axis label for the convolved spectrum in the online documentation be counts/keV instead of counts/s/cm^2?

Is there a link to the CI tests that can be included in README showing which software platforms BXA currently runs on (eg. linux/macos) as well the python versions?

These could be useful for building examples that work out-of-the-box as well as informing the users of pitfalls in installing the software.

Hi Johannes,

Just checking to see if there's a plan to make BXA work with the latest Astropy. It looks like things are getting wonky due to a change to the astropy.fits.io API.

Thanks,

-Brian

The gen.xspec script functions properly and the examples all run correctly.

It would be nice if the large amount of output from the sampler was removed from the docs / README as it is difficult to find the relevant information without a lot of scrolling. Is the output needed for any purpose?

The documentation refers to "prior transformations" throughout the docs. This is a reference to transforming the nested sampling focus on transforming all priors from a unit cube.

This could be confusing as users who do not realize that this is specifically for nested sampling and not a form of MCMC that is provided in both XSPEC (https://asd.gsfc.nasa.gov/XSPECwiki/Markov_Chain_Monte_Carlo) and Sherpa (https://sherpa.readthedocs.io/en/latest/mcmc/index.html). It might be useful at the introduction (since it is not in the name) to clarify this? The connection to nested sampling is clearly there, but a note or something similar could be more instructive.

Perhaps is would useful to define what a transformation is here:
https://github.com/JohannesBuchner/BXA/blob/master/bxa/xspec/solver.py#L104

There are references to pymultinest throughout the docs.

What are all the available priors and can a user add custom priors?

In the examples, a Gaussian prior must be specified with limits. Is it not possible to have an unbounded Gaussian prior?

It seems there is a wealth of functionality with priors in the examples that is perhaps not documented and a naive user could miss these capabilities.

I am running BXA 4.0.5, UltraNest 3.4.4, Python 3.8.13, Sherpa 4.14.0, on MacOS High Sierra 10.13.6.

I just installed BXA over the weekend and was trying out the xagnfitter.py script, but ran into an error early in the process. After setting up the environment as well as the input files (filenames.txt, the galNH and redshift files), and linking the model directory, I ran the script as directed using 'python xagnfitter.py'. The code crashed somewhere around the 'except' statement on line 262 of the script. Below I've included a partial log of the terminal output and the traceback.

Increasing parameters again...
3 parameters, aic=297.37
4 parameters, aic=297.45
5 parameters, aic=296.76
Final choice: 2 parameters, aic=295.52

Adding Gaussian#1
largest remaining discrepancy at 8.660keV[99], need 447 counts
placing gaussian at 8.66keV, with power 0.01026600881157122
Traceback (most recent call last):
  File "xagnfitter.py", line 259, in <module>
    bkg_model = CachedModel(get_bkg_model(id))
  File "<string>", line 1, in get_bkg_model
  File "/anaconda3/envs/ciao-4.14/lib/python3.8/site-packages/sherpa/astro/ui/utils.py", line 9374, in get_bkg_model
    raise ModelErr('nobkg', bkg_id, id)
sherpa.utils.err.ModelErr: background model 1 for data set 1 has not been set

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "xagnfitter.py", line 265, in <module>
    bkg_model = auto_background(id)
  File "/anaconda3/envs/ciao-4.14/lib/python3.8/site-packages/bxa/sherpa/background/pca.py", line 416, in auto_background
    bkgmodel.fit()
  File "/anaconda3/envs/ciao-4.14/lib/python3.8/site-packages/bxa/sherpa/background/pca.py", line 385, in fit
    g.Sigma.max = x[-1] - x[0]
  File "/anaconda3/envs/ciao-4.14/lib/python3.8/site-packages/sherpa/utils/__init__.py", line 171, in __setattr__
    object.__setattr__(self, name, val)
  File "/anaconda3/envs/ciao-4.14/lib/python3.8/site-packages/sherpa/models/parameter.py", line 230, in _set_limit
    raise ParameterErr('edge', self.fullname,
sherpa.utils.err.ParameterErr: parameter g_1_99.Sigma has a hard maximum of 20

I spoke to Peter Boorman about this, and he ran into the same issue described above with a separate set of data and bxa installation when diagnosing this issue. He discovered a workaround in the pca.py file (/anaconda3/envs/ciao-4.14/lib/python3.8/site-packages/bxa/sherpa/background/pca.py). Changing in pca.py the line:

g.Sigma.max = x[-1] - x[0]

to:

if x[-1] - x[0] >= 20.:
	g.Sigma.max = 20.
else:
	g.Sigma.max = x[-1] - x[0]

seems to work for both of us and bypasses what appears to be a hard limit for the sigma maximum for the Gaussian model in Sherpa. After implementing the workaround, the xagnfitter.py script executes properly and produces the desired plots.

I've tried to follow the docs and the API in order to visualise the different components in my model and did not succeed. I have 3 model parameters. Here is what I tried:

m = xspec.Model("constant*(bknpower + gauss)")
Then I set the priors on some parameters and solved and get a nice corner plot etc. And trying to visualise the two models I'm interested in with this:

data = solver.posterior_predictions_convolved(component_names=['ignore','bknpower','gauss'],nsamples=100)

And I can only see on black curve passing through the data points, no separtion of the components. Adding at the end plt.gca().legend() only showed one legend entry for the data and nothing for the models.

I even tried to add plot_args=[{'color': 'red'},{'color:'green'},{'color':'blue'}] in the call to solver.posterior_predictions_convolved() with no success.

I suppose I'm doing something wrong.

In the same way of insufficient info in the docs, how can I get the actual total model (or per component) with the best-fit parameters? Currently it is plotted with plt.figure() but if I want to use the predicted model in subplots then it's not very convenient. I can see the dataset data['models'] with shape (100,1,1306) but I cannot make any use of it with the available info.

Can you please add some examples in the documentations?

Thanks in advance,

It seems like it is impossible to have the posterior plots when the joint likelihood of the source+background and the background are used. BXAsolver() performs fine, but the binning expects integral counts in the plot background subtraction. This fails each time because it is likely that one or more source channels have counts lower than the corresponding background channel. Or the source may have 0 counts in a channel.

Could the background having 0 counts lead to the mess? Should a tiny bit of rebinning of the original data (to have at least 1 count per bin, like XSPEC recommends) be allowed to take this into account? Ignoring artificial features or smearing of information, what other problems may this lead to?

When using more than one data-group in XSPEC, which happens when you are fitting more than one spectrum, and using the same model for different data-groups, one will end up in having model parameters with the same name and model-index for different data-groups.
So far this isn't a problem. However, the BXA function "store_chain" creates a FITS-Table to be used for the XPSEC task chain, where it names each column with parameter name and its model-index:
columns = [pyfits.Column(name='%s__%d' % (t['name'], t['index']), (Line 33 of __init__.py)
Hence, when letting the same parameter be free for more than one data-group (e.g. Instrumental Background Normalization), "pyfits" will crash because it can't create a FITS-Table with identical columns. I guess, one simple workaround would be to create the column names beforehand and clean for identical names. I just removed this part from the BXA code, since I don't really know how to make use of those chain-FITS-Files.

• BXA version: 4.0.5
• UltraNest version: 3.4.4
• Python version: 3.9
• Xspec or Sherpa and version: Xspec 12.11.1
• Operating System: Debian GNU/Linux 11 (bullseye)

Description

While attempting to parallelize BXA with mpi, h5py file is created but locked. After following recommendation in previous (closed) bxa issue thread here, and attempting to reinstall all dependencies, problem persists, but with new error.

I read many forums regarding the errors, and they have recommended reinstalling dependencies, it seems as though the h5py file is corrupted while being created.

What I Did

Old error:

Traceback (most recent call last):
  File "/home/jpbreuer/Scripts/bxa_test.py", line 373, in <module>
    results = solver.run(resume=True)
  File "/home/jpbreuer/.local/lib/python3.9/site-packages/bxa/xspec/solver.py", line 188, in run
    self.results = solve(
  File "/home/jpbreuer/.local/lib/python3.9/site-packages/ultranest/solvecompat.py", line 55, in pymultinest_solve_compat
    sampler = ReactiveNestedSampler(
  File "/home/jpbreuer/.local/lib/python3.9/site-packages/ultranest/integrator.py", line 1077, in __init__
    self.pointstore = HDF5PointStore(storage_filename, storage_num_cols, mode='a' if resume else 'w')
  File "/home/jpbreuer/.local/lib/python3.9/site-packages/ultranest/store.py", line 187, in __init__
    self.fileobj = h5py.File(filepath, **h5_file_args)
  File "/home/jpbreuer/.local/lib/python3.9/site-packages/h5py/_hl/files.py", line 507, in __init__
    fid = make_fid(name, mode, userblock_size, fapl, fcpl, swmr=swmr)
  File "/home/jpbreuer/.local/lib/python3.9/site-packages/h5py/_hl/files.py", line 232, in make_fid
    fid = h5f.open(name, h5f.ACC_RDWR, fapl=fapl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5f.pyx", line 106, in h5py.h5f.open
BlockingIOError: [Errno 11] Unable to open file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable')

Updated error:

Traceback (most recent call last):
  File "/home/jpbreuer/Scripts/bxa_test.py", line 128, in <module>
    results = solver.run(resume=True)
  File "/usr/local/lib/python3.9/dist-packages/bxa/xspec/solver.py", line 188, in run
    self.results = solve(
  File "/usr/local/lib/python3.9/dist-packages/ultranest/solvecompat.py", line 55, in pymultinest_solve_compat
    sampler = ReactiveNestedSampler(
  File "/usr/local/lib/python3.9/dist-packages/ultranest/integrator.py", line 1077, in __init__
    self.pointstore = HDF5PointStore(storage_filename, storage_num_cols, mode='a' if resume else 'w')
  File "/usr/local/lib/python3.9/dist-packages/ultranest/store.py", line 187, in __init__
    self.fileobj = h5py.File(filepath, **h5_file_args)
  File "/usr/lib/python3/dist-packages/h5py/_debian_h5py_serial/_hl/files.py", line 387, in __init__
    fid = make_fid(name, mode, userblock_size,
  File "/usr/lib/python3/dist-packages/h5py/_debian_h5py_serial/_hl/files.py", line 187, in make_fid
    fid = h5f.create(name, h5f.ACC_EXCL, fapl=fapl, fcpl=fcpl)
  File "h5py/_debian_h5py_serial/_objects.pyx", line 54, in h5py._debian_h5py_serial._objects.with_phil.wrapper
  File "h5py/_debian_h5py_serial/_objects.pyx", line 55, in h5py._debian_h5py_serial._objects.with_phil.wrapper
  File "h5py/_debian_h5py_serial/h5f.pyx", line 108, in h5py._debian_h5py_serial.h5f.create
OSError: Unable to create file (unable to open file: name = 'bxatest/results/points.hdf5', errno = 17, error message = 'File exists', flags = 15, o_flags = c2)
Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/h5py/_debian_h5py_serial/_hl/files.py", line 185, in make_fid
    fid = h5f.open(name, h5f.ACC_RDWR, fapl=fapl)
  File "h5py/_debian_h5py_serial/_objects.pyx", line 54, in h5py._debian_h5py_serial._objects.with_phil.wrapper
  File "h5py/_debian_h5py_serial/_objects.pyx", line 55, in h5py._debian_h5py_serial._objects.with_phil.wrapper
  File "h5py/_debian_h5py_serial/h5f.pyx", line 88, in h5py._debian_h5py_serial.h5f.open
OSError: Unable to open file (truncated file: eof = 96, sblock->base_addr = 0, stored_eof = 2048)

Hello,

I have encountered the following error while running the BXA script using mpiexec -np 2 python3 script.py.

Traceback (most recent call last):
  File "BXA_borus02D_FF_noSz_allsets.py", line 73, in <module>
    solver = bxa.BXASolver(transformations=priors, outputfiles_basename=outputfiles_basename)
  File "/home/dl/.local/lib/python3.8/site-packages/bxa/xspec/solver.py", line 132, in __init__
    os.mkdir(outputfiles_basename)
FileExistsError: [Errno 17] File exists: 'BXA_borus02D_FF_noSz_allsets/'

The script is an unmodified version of what I am able to run without MPI. I have tried deleting the directory and starting over again but the issue persists. Is there a fix for this?

Following are the version numbers of various tools.
BXA: 4.0.5
mpiexec: 4.0.3
OS: Ubuntu 20.04

Please let me know if you need any other information.

Thank you in advance!

• BXA version: 4.1.1
• UltraNest version: 3.5.7
• Python version: 3.9
• Xspec or Sherpa and version: ciao 4.15
• Operating System: MacOS 12.6

Description

As a way to understand how the NS works and to visualize the evolution of likelihood as a function of iterations, I'm looking at the logL column in the weighted_post.txt.

I'm using unbinned X-ray spectral data and the set_stat("wstat") in Sherpa.
The background is extracted from the same observation.
The BXA file weighted_post.txt looks this :

In this case, the BXA maximized the logL value (the last rows with higher value are better) whereas in the sherpa or Xspec implementation, we minize Cstat or Wstat.
How is this logL value is calculated ? Can it be related to Wstat value computed by sherpa ?

One the thing I'm trying is comparing the Wstat values from sherpa's fit to the solutions found by BXA.

When I use chi statistic in XSPEC, if the reduced chisq is approximately 1 then I will consider the model is appropriate.
I am wondering if there is such a convenient criterion in BXA? If there is no such criterion, so how to quantitively justify whether the model fits the data?

At the moment BXA (actually pymultinest) only uses one core. How can I make it use more than one?

• BXA version: 4.0.5
• UltraNest version: 3.3.3
• Python version: 3.10.8
• Xspec or Sherpa and version: Xspec 12.12.1
• Operating System: Arch Linux, Kernel 6.0.2-arch1-1

Description

I have a complicated model set up with several named models, 49 data groups with different spectra, and many linked and free parameters spread out between several data groups. I have no problems defining the priors, but I begin having problems when trying to run the solver to start fitting the free parameters.

At first, the problem could be reproduced before running the solver in the set_parameters(transformations=solver.transformations,values=values) line, but the problem persists afterwards to the solver instance.

I have tried renaming the models, and simplifying to a simpler use case with only a few free parameters in a few data groups. I expect that the issue comes from the priors. Despite the priors being able to be defined with no problems, the solver does not seem to like that the free parameters to fit are spread out over several different data groups. I think the error comes from line 96 in solver.py, and I have tried fiddling with it but can not seem to resolve it.

What I Did

cxb = Model("constant(apec + phabs(apec + apec + powerlaw))constant + constant*constant(phabs*apec)constant + constant*constant(phabs*apec)constant + constant*constant(phabs*apec)constant")

# .... Freeze everything except a few params ....

prior_gh_kT = bxa.create_uniform_prior_for(cxb, AllModels(1).apec_4.kT)
prior_gh_norm = bxa.create_loguniform_prior_for(cxb, AllModels(1).apec_4.norm)
prior_pow_norm = bxa.create_loguniform_prior_for(cxb, AllModels(1).powerlaw.norm)
prior_clst_kT = bxa.create_uniform_prior_for(cxb, AllModels(14).apec_11.kT)
prior_clst_norm = bxa.create_loguniform_prior_for(cxb, AllModels(14).apec_11.norm)
prior_grp_kT = bxa.create_uniform_prior_for(cxb, AllModels(26).apec_16.kT)
prior_grp_norm = bxa.create_loguniform_prior_for(cxb, AllModels(26).apec_16.norm)
prior_fil_kT = bxa.create_uniform_prior_for(cxb, AllModels(38).apec_21.kT)
prior_fil_norm = bxa.create_loguniform_prior_for(cxb, AllModels(38).apec_21.norm)

priors = [prior_gh_kT, prior_gh_norm, prior_pow_norm, prior_clst_kT, prior_clst_norm, prior_grp_kT, prior_grp_norm, prior_fil_kT, prior_fil_norm]

solver = bxa.BXASolver(transformations=priors, outputfiles_basename="bxa_all_cxbclustgroupfil")
# values = solver.prior_function(np.random.uniform(size=len(solver.paramnames)))
# set_parameters(transformations=solver.transformations,values=values)
results = solver.run(resume=True)

Traceback

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In [33], line 1
----> 1 results = solver.run(resume=True)

File ~/.local/lib/python3.10/site-packages/bxa/xspec/solver.py:193, in BXASolver.run(self, evidence_tolerance, n_live_points, wrapped_params, **kwargs)
    190 Lepsilon = kwargs.pop('Lepsilon', 0.1)
    192 with XSilence():
--> 193         self.results = solve(
    194                 self.log_likelihood, self.prior_function, n_dims,
    195                 paramnames=self.paramnames,
    196                 outputfiles_basename=self.outputfiles_basename,
    197                 resume=resume, Lepsilon=Lepsilon,
    198                 n_live_points=n_live_points, evidence_tolerance=evidence_tolerance,
    199                 seed=-1, max_iter=0, wrapped_params=wrapped_params, **kwargs
    200         )
    201         logls = [self.results['weighted_samples']['logl'][
    202                 numpy.where(self.results['weighted_samples']['points'] == sample)[0][0]]
    203                 for sample in self.results['samples']]
    204         self.posterior = self.results['samples']

File ~/.local/lib/python3.10/site-packages/ultranest/solvecompat.py:47, in pymultinest_solve_compat(LogLikelihood, Prior, n_dims, paramnames, outputfiles_basename, resume, n_live_points, evidence_tolerance, seed, max_iter, wrapped_params, verbose, speed, **kwargs)
     44 if not verbose:
     45     outputkwargs = dict(viz_callback=False, show_status=False)
---> 47 sampler = ReactiveNestedSampler(
     48     paramnames, LogLikelihood, transform=Prior,
     49     log_dir=outputfiles_basename, resume='resume' if resume else 'overwrite',
     50     wrapped_params=wrapped_params, draw_multiple=False, vectorized=False,
     51     **outputkwargs)
     53 if speed == "safe":
     54     pass

File ~/.local/lib/python3.10/site-packages/ultranest/integrator.py:1095, in ReactiveNestedSampler.__init__(self, param_names, loglike, transform, derived_param_names, wrapped_params, resume, run_num, log_dir, num_test_samples, draw_multiple, num_bootstraps, vectorized, ndraw_min, ndraw_max, storage_backend, warmstart_max_tau)
   1093 self.ndraw_max = ndraw_max
   1094 self.build_tregion = transform is not None
-> 1095 if not self._check_likelihood_function(transform, loglike, num_test_samples):
   1096     assert self.log_to_disk
   1097     if resume_similar and self.log_to_disk:

File ~/.local/lib/python3.10/site-packages/ultranest/integrator.py:1157, in ReactiveNestedSampler._check_likelihood_function(self, transform, loglike, num_test_samples)
   1153 p = transform(u) if transform is not None else u
   1154 assert np.shape(p) == (num_test_samples, self.num_params), (
   1155     "Error in transform function: returned shape is %s, expected %s" % (
   1156         np.shape(p), (num_test_samples, self.num_params)))
-> 1157 logl = loglike(p)
   1158 assert np.logical_and(u > 0, u < 1).all(), (
   1159     "Error in transform function: u was modified!")
   1160 assert np.shape(logl) == (num_test_samples,), (
   1161     "Error in loglikelihood function: returned shape is %s, expected %s" % (np.shape(logl), (num_test_samples,)))

File ~/.local/lib/python3.10/site-packages/ultranest/utils.py:131, in vectorize.<locals>.vectorized(args)
    129 def vectorized(args):
    130     """Vectorized version of function."""
--> 131     return np.asarray([function(arg) for arg in args])

File ~/.local/lib/python3.10/site-packages/ultranest/utils.py:131, in <listcomp>(.0)
    129 def vectorized(args):
    130     """Vectorized version of function."""
--> 131     return np.asarray([function(arg) for arg in args])

File ~/.local/lib/python3.10/site-packages/bxa/xspec/solver.py:157, in BXASolver.log_likelihood(self, params)
    156 def log_likelihood(self, params):
--> 157         set_parameters(transformations=self.transformations, values=params)
    158         like = -0.5 * Fit.statistic
    159         # print("like = %.1f" % like)

File ~/.local/lib/python3.10/site-packages/bxa/xspec/solver.py:96, in set_parameters(transformations, values)
     94         pars += [t['model'], {t['index']:v}]
     95 # print(pars)
---> 96 AllModels.setPars(*pars)

File /opt/heasoft-6.30/x86_64-pc-linux-gnu-libc2.35/lib/python/xspec/model.py:1173, in ModelManager.setPars(self, *args)
   1168 if (iArg == len(args) or
   1169       isinstance(args[iArg], Model)):
   1170     if len(rawValArgs):
   1171         # These args belong to the previous model.
   1172         # It's time to collate them.
-> 1173         indAndVals = _collateSetPars(mod, *rawValArgs)
   1174         nPars = len(indAndVals[0])
   1175         allModNames += nPars*[mod.name]

File /opt/heasoft-6.30/x86_64-pc-linux-gnu-libc2.35/lib/python/xspec/model.py:1287, in _collateSetPars(mod, *argsTuple)
   1271 """Organize the user's various input args into containers compatible with
   1272       the C++ setPars and setGlobalPars functions.
   1273    
   (...)
   1284    
   1285 """
   1286 # iPar is the 1-based GLOBAL index.
-> 1287 startIdx = mod.startParIndex
   1288 iPar=startIdx
   1289 # First gather the indices and vals in a dictionary to prevent
   1290 #  multiple appearances of a single par index in the final output
   1291 #  lists.  When multiple entries are found, the later value will step
   1292 #  on the earlier one and the user will be given a warning.
   1293 # Following any dictionary entry, the iPar count will be set to the
   1294 #   last iPar key in the dictionary PLUS ONE. 

File /opt/heasoft-6.30/x86_64-pc-linux-gnu-libc2.35/lib/python/xspec/model.py:467, in Model._getStartParIndex(self)
    466 def _getStartParIndex(self):
--> 467     return _pyXspec.getModelTuple(self._handle)[4]

Exception: Error: Python Model object reference no longer corresponds to
          an actual XSPEC model.

PR #51 merged code that shows nice notebooks.

We should

• add these on the website under "Bayesian workflow" using nb-sphinx, so that they are run when the documentation is built
• link from the BXA/xspec manual, so that users can switch to it
• add them to the CI, so that the workflow is run and is ensured to work

I just released a new version (3.4.0) @alexkolo @ruizca @anetasie @bwgref @cdeil. Please test it if you can.

Hi, some more comments on the example code:

• It seems that some of the sherpa examples rely on npyinterp and python 2? Could this be updated to python 3? Please correct me if I am wrong.
• The example code in docker/ also seems a bit outdated (e.g. doesn't use ultranest and python 3). OK if you want to keep this for reference, but then an up-to-date docker example should be added too.
• What is the difference between xagnfitter.py detailed in the docs and fitagn.py in the docker examples?
• When running gal.py I get a connection timeout to swift.ac.uk (looks likely to be a temporary problem on their side, so just to notify you)
• Can you fix the latest Travis CI build? Seems that the examples run fine and it should work

• BXA version: 4.1.1
• UltraNest version: 3.6.2
• Python version:3.11.5
• Xspec or Sherpa and version: xspec 12.13.1
• Operating System: Linux CentOS 7

Description

Hi,

I am trying to fit a single Athena xifu spectrum, much like the example at https://github.com/JohannesBuchner/BXA/blob/master/examples/xspec/example_simplest.py.

This works brilliantly for the case where I am fitting the whole channel range, but as soon as I ignore some channels/energies I get a segmentation fault e.g.

[  +0.008827] python[254852]: segfault at 793c000 ip 00002b84c0e8c87c sp 00007ffe5926bed0 error 4
[  +0.000006]  in libXSStat.so[2b84c0e6a000+42000]

The line that triggers it is

spectrum.ignore("**"); spectrum.notice("0.5-3.0")

and from the output this gives 29931 channels (1,29931) ignored in spectrum.

The segfault seems to occur after all of the usual xspec- like output, and the program then prints out the number of channels ignored and then instantly crashes. I don't know if I am doing something wrong, if there is a problem with my install, or some other issue.

This happens both when I run with mpiexec in parallel, or when I run with just one task executing the python file.

Fred

Is this open source?
If yes, could you please add a license?

Here are the issues with getting into a state of being able to run the python example_simplest.py

1. The main required dependency for the fitting, ultranest is not installed automatically either by pip or python setup.py install

2. The setup currently does not run on MacOS.

All of the following is run in a clean virtual environment on Ubuntu 20.04.2 with python 3.8 (no conda)

After completing the install for CIAO with the following steps:

bash ciao_install
pip install sherpa # (installs numpy)
python setup.py install # inside BXA folder cloned from GitHub

the following steps need to be run to avoid errors:

pip install ultranest

From this point I try to run python example_simplest.py

Result:

WARNING: imaging routines will not be available,
failed to import sherpa.image.ds9_backend due to
'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH'
WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available
read RMF file athenapp_ir_b4c_wfi_withfilter_fov40.0arcmin_avg.rsp
Traceback (most recent call last):
  File "example_simplest.py", line 20, in <module>
    set_model(xspowerlaw.mypow)
NameError: name 'xspowerlaw' is not defined

It would be useful to modify the following:

• add all required dependencies to the install/setup script

• include examples that work with Sherpa alone as it is the easiest entry point into testing the code.

• as BXA is a plugin to the XSPEC/Sherpa analysis software, warning to the user upon installation that the environment does not detect the required essential dependencies to run upon install.

• perhaps an end-to-end example in the documentation (including data) would be helpful.

While it is not BXA should not be responsible for upkeep of XSPEC/Sherpa's installation procedures and capabilities, its focus as a plugin to these software means that it needs to perhaps guide the user thru the installation of its core functionality.

In regards to the software not installing on MacOS, this is again an issue of the core packages, but should be reflected in the caveats / platforms listed in the README and documentation. A few of these may be addressed with #21

I just wanted make users aware that the code crashes if the string outputfiles_basename is too long. I think it can't be more than 100 characters.
As a workaround, one can create all the files at a temporary short path and move them afterwards to the desired path, for instance with shutil.move(Current, Target).
Note, that some BXA functions will try to access files at the temporary short path. Hence, one should complete the BXA analysis before moving/archiving the files.

Hi,
Thank you for developing the bxa software.

About bxa.sherpa.background.pca.auto_background() feature, there are 2400 channels in spectra extracted from EPIC/MOS detectors. However, there are only xmm_4096.json and xmm_1024.json.

How could the pca based background model be used for MOS spectra?

-Shifu

Currently the background model examples depend on both Sherpa and XSPEC.

Is there a way to separate these into analysis framework independent examples or is the functionality tied to having both frameworks installed?