Add sherpa wrapper #52
Draft
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This is based on the convert_xspec_user_model script written by the Science Data Systems group of the Chandra X-ray Center. Issues are - only built on a Linux conda install of CIAO and I have not actually run it (as I don't know of a sensible set of parameters) - should we wrap the utils.cpp::jetinterp routine (it is labelled as being required by ISIS) - how about supporting more-comlex combinations of models (that is supporting models like those used n the Examples/ directory)? The parameter defaults and ranges are based on bhjet.sl.
The code used to create the wrapper routine used by the Sherpa interface is useful for building an XSPEC interface, so separate it out into it's own file. This also changes exactly what is being wrapped: for the Sherpa case we were wrapping the jetmain function, but a review of the bhjet.sl file shows me that we really want to combine jetmain and jetinterp. It would be possible to use the ISIS approach in Sherpa - that is provide a Python model evaluator that calls both jetmain and jetinterp - but it's not obvious if that is sensible for XSPEC. So we wrap a single routine that does both. This makes it harder to tweak (e.g. decide to change the "base" grid the model is evaluated on). Add a "lmodel" interface for the model, to tell XSPEC how to create the model. This is based on the bhjet.sl file. There has been no attempt to build an XSPEC model using this code.
DougBurke
force-pushed
the
add-sherpa-wrapper
branch
from
d6a47f3 to
d3ba3d9
Compare
This doesn't actually fix the problem I am seeing, but is a useful change pypa/pip#9542
Fortunately no FORTRAN code in use here.
|
I've stuck ni sme debugging statements so that the following will create which shows the inputs to |
|
I'm so sorry for getting to this only now! Unfortunately that particular function goes back to the earliest and most horribly documented days of the code (despite my best efforts...). |
Is this correct or should energ/phot be set to NBINS - 1 and egrid to NBINS?
|
For reference, an example "session" in Sherpa - using CIAO 4.14 installed with conda and only installing the sherpa-relevant parts we can now create a model instance [I have a custom sherpa prompt] Let;'s fake a data set: We can now associate the model NOTE I should have been able to say and the model looks like
|
|
Doh, I know why that last line should have returned |
Merged
|
So, with the following file (renamed hack.cpp), which "just" calls If I build this in then run with then the output contains the following, which lists the output of Here's the full output So it's generating NaN's - which suggests I'm doing something wrong calling |
|
This is very very confusing to me still. What happens if you use a wider grid (say, from radio to hard X-rays)? If the issue goes away in that case I would say everything is fine, since the model is not supposed to be used for fitting X-rays only anyway. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
The primary change is to allow the BHJet model to be used from Sherpa. There's also some support for building an XSPEC interface [*], but I have not got the time/energy to do the useful part of getting the build to work (I don't think it would be hard).
There are a few minor tweaks to the README to make it read better on GitHub when parsed as markdown, and the use of gsl-config rather than hard-coding the values when building the BHJet code.
The default parameter values and ranges are based on bhjet.sl, as pointed out by Sera, but all mistakes are my own.
Is there a
Input/ip.datfile that could be included or referred to let users check that they get sensible results.[*] the Sherpa interface is built on Sherpa's XSPEC interface, so the
wrap_xspec.cppfile that I add which should be useful for XSPEC is also useful for Sherpa.Sherpa installation
This is added to the
BHJet/README.mdfile, but installation requires CIAO 4.14, preferably installed with conda, and then you can sayThen you can load it into Sherpa and use it as a model with
Implementation note
So
wrap_xspec.cppprovides an "XSPEC-compatible" interface using the C (rather than FORTRAN or C++) API - see https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixLocal.htmlThis wrapper routine, wonderfully named
xspec_jetinterp, is meant to be close to the ISISbhjet_fitroutine - i.e. it does both the call tojetmainand then the call tojetinterp, rather than having a Python-level routine analogous tobhfit_jet. This is less open to user tweaks, but easier to handle (and I think required if you wanted to use the model in XSPEC; in Sherpa we could use a similar approach to ISIS).Unfortunately as I'm not an ISIS user then I appear to have made a mistake somewhere (see next section). In the S-Lang/ISIS call to
bhjet_fit(), are theloandhiarguments always in Angstroms?Current problems
Importing the Sherpa model "directly" (a "lower-level" interface than users are expected to use) shows that I have not correctly handled the call to both
jetmainand thenjetinterp, but I haven't got the time to track this down just yetFor the really low-level interface, first get the default parameter values as I'm too-lazy to type them in myself: note that we have 28 parameters because the last-one is the normalization.
Now we can access the "C++ code" (essentially this is equivalent to slirp creating a S-Lang functino that calls the C++ code):
This routtine takes three array/list arguments: the parameters, the low-edge, and the high-edge of the bins (it's done this way so I can just re-use the code we have for interfacing to X-SPEC models, but we could change it if need be):
So, am I calling it wrong? The lo/hi edges are assumed to be keV.