johannesbuchner / bxa Goto Github PK
View Code? Open in Web Editor NEWBayesian X-ray analysis (nested sampling for Xspec and Sherpa)
Home Page: https://johannesbuchner.github.io/BXA/
License: GNU General Public License v3.0
Bayesian X-ray analysis (nested sampling for Xspec and Sherpa)
Home Page: https://johannesbuchner.github.io/BXA/
License: GNU General Public License v3.0
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?
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?
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.
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.
Running solver.run()
cannot be interrupted in jupyter lab.
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
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.
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?
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)
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.
# Name Version Build Channel
bxa 4.0.5 pyh44b312d_0 conda-forge
xorg-libxau 1.0.9 h7f98852_0 conda-forge
# Name Version Build Channel
ultranest 3.5.7 py310hde88566_0 conda-forge
$ python --version
Python 3.10.9
xspec-modelsonly 12.12.1c h707b37c_0 https://cxc.cfa.harvard.edu/conda/ciao
Sherpa 4.15.0
$ 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
I am in theciao-4.15
environment, 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!
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).
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?
Hi,
there is a Python 3 issue in BXA/bxa/sherpa/background/fitters.py, where the function file() is still used.
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
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)?
At the moment BXA (actually pymultinest) only uses one core. How can I make it use more than one?
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?
PR #51 merged code that shows nice notebooks.
We should
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.
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?
Hi, a few comments/questions on the documentation as part of the review:
cython
, tqdm
, h5py
and ultranest
as a minimal list to get things running.examples/
, but would be helpful if at least one case was treated notebook-style on the documentation pages.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.gal.py
is a typo as refers to another scriptHello,
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!
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,
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.
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!
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?
Here are the issues with getting into a state of being able to run the python example_simplest.py
The main required dependency for the fitting, ultranest
is not installed automatically either by pip or python setup.py install
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
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
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?
I tried to install BXA using this guide.
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.
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.
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.
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)
---------------------------------------------------------------------------
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.
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
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.
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)
Could you add "_{ndata}" to the name of the output file in autobackgroundmodel/init.py, since it is required by autobackgroundmodel/fitbkg.py?
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.
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.
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 ;)
When passing the parameter evidence_tolerance the run does not work
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.
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.
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:
BXA/bxa/sherpa/background/pca.py
Line 359 in 5a6ee5f
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.
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.
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
Hi, some more comments on the example code:
npyinterp
and python 2? Could this be updated to python 3? Please correct me if I am wrong.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.xagnfitter.py
detailed in the docs and fitagn.py
in the docker examples?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)A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.