Add xspec convolution api #750
Conversation
DougBurke
force-pushed
the
add-xspec-convolution-api
branch
from
7bb02b0
to
052ab5f
Compare
DougBurke
force-pushed
the
add-xspec-convolution-api
branch
from
052ab5f
to
ddf7a18
Compare
DougBurke
force-pushed
the
add-xspec-convolution-api
branch
from
ddf7a18
to
41e8e97
Compare
|
I've rebased this to check the tests still work after the regrid changes in #747 |
|
Will rebase this soon because of #732 (a change to a comment line, so no user-visible changes) |
DougBurke
force-pushed
the
add-xspec-convolution-api
branch
from
41e8e97
to
27f6c50
Compare
There was a problem hiding this comment.
I like where this is going and I think it's an important addition to Sherpa supporting all the XSPEC models out there. I can't say that I understand all the code, but the tests look pretty extensive.
From an UI perspective, I wonder if we need to do more for the way that models are defined.
@DougBurke says
That is, you say xscflux.cfl(mdlexpr) to create a model instance called cfl which is applied to the model stored in mdlexpr (which could be a single component or multiple components). You can not use a model instance (e.g. cfl) directly - that is cfl * xspowerlaw.pl will not work.
I agree the "*" should not work, but how about clf(xspowerlaw)? I guess that won't work, because "(..)" calls the model for evaluation, but it's what I naively expected to do. On the other hand, that's how the OO layer allows it (https://github.com/sherpa/sherpa/pull/750/files#diff-ac7835014238b3cefc8770664f5af6e0R493). Is there a way we can make cfl(mdlexpr) (i.e. with a cfl that's been used before e.g. for some other dataset) work, in addition to xscflux.cfl(mdlexpr)? (Or maybe I'm mistaken and it already works?). If not, can we use some other binary operator to express convulution models, e.g. cfl @ mdlexpr or clf | mldexpr (although the later one looks like a pipe and thus the ordering would be confusing (i.e. should it be mldexpr | clf).
I do think that this needs to be explained in the rst docs (I know that this PR is not at the stage yet, but I'm just saying). Most of the text could just be taken from the description of he issue.
|
What I meant was you can not say |
|
@DougBurke Fantastic. I agree that that's exactly how it should be. I wonder (not as part if this PR) if we can treat the pile-up model the same way. jdpileup is essentially a convolution model, it's the only pile-up model that's implemented and it currently has all this extra UI associated with just that one model (set_pileup, there is an issue somewhere whose number I don' remember now that you can't unset a pile-up model, etc.). |
DougBurke
force-pushed
the
add-xspec-convolution-api
branch
from
06bb77a
to
aa8162d
Compare
|
I've rebased to pick up the python 3.8 changes. It has also pointed out that the calc_energy_flux-related test had to be changed (thanks to #756) so I've added a commit to fix that. |
|
@hamogu - I've added more documentation. I don't really know what else to add here. edited to add although saying that, I'm not sure how if the "ui" documentation needs more work, as this has been mainly focussed on the lower level parts of the API. Some of the models do contain hopefully-useful docstrings: e.g. |
|
@DougBurke I build it and tested. The API is good. Naming a convolution model is consistent with the Sherpa general treatment of the models, so it should be clear to the users. One question that I have is whether in case of I know in your example the |
|
@anetasie - the cflux parameter is the log of it, so and also I can't remember what the default limits that cflux uses (its emin/max parameters, or whatever they're called) EDITED TO ADD just checked, and the default Emin/Emax are 0.5 to 10. It's a bit hard to compare your results as you need to a) make sure that pl.norm = 1 and is frozen (this is mentioned in the XSPEC help, where the link https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node280.html is valid for 12.11.0, but it's not particularly clear: if you don't have the norm frozen then it is degenerate with the lg10Flux value, and you want it set to 1 otherwise the flux will be off by 1/norm b) it's not clear from above what happened in the lines in-between the fit and the calc_energy_flux call, so hard to know what has differed, but you'd want to |
|
@hamogu - do you still have comments on this, or are the doc updates sufficient? |
|
Fine with me! |
|
This will need rebasing now that XSPEC 12.11.0 support (#772) has been merged. |
There was a problem hiding this comment.
Comment should explain why one does (or does not) do something, as opposed to mimic what the code does:
def calc(self, pars, rhs, *args, **kwargs):
"""Evaluate the convolved model.
Note that this method is not cached by
sherpa.models.modelCacher1d (may change in the future).
I still don't know why you decided not to cache these functions.
|
@dtnguyen2 FYI I'm going to be rebasing this patch so that it will work now we've got support for XSPEC 12.11.0. It foesn't change any of the logic, but the commits will get changed. |
DougBurke
force-pushed
the
add-xspec-convolution-api
branch
from
979983a
to
a93b6eb
Compare
The Python "standard" for docstrings says that the first line should be a short synopsis of the function.
A number of reasons:
|
|
I agree with @DougBurke : The docstring of a function is meant to users on how to use that particular function. For that, they don't have to know why something was done, just how. Of course, inline comments in the code that are only seen by other developers by all means should explain why something was done as @dtnguyen2 suggests. |
This commit takes code that was used in the SDS "sherpa contribute"
routines and moves them into Sherpa itself. It adds Python-level
support for the XSPEC convolution models. In the SDS version I had
explicit load_xs<convolution model> routines, following on from
the existing load_psf/table style routines, but I have decided that
this adds an extra layer of complexity that may not be needed. So
this commit emphasizes "direct" creation of these convolution
models. The load_xs<> style routines could be added on top of this
commit if needed.
The XSPEC cflux convolution model (that calculates the energy
flux of a model) can be used like:
from sherpa.astro import ui
ui.load_pha('3c273.pi')
ui.subtract()
ui.notice(0.5, 7)
ui.set_source(ui.xsphabs.gal * ui.xscflux.cfl(ui.powlaw1d.pl))
pl.ampl.frozen = True
ui.fit()
ui.conf()
You can also use them with regridded models, as shown below
convsource = ui.xszashift.zsh(ui.box1d.box + ui.const1d.bgnd)
zsh.redshift = 0.8
rgrid = np.arange(0.1, 23, 0.01)
box.xlo = 10
box.xhi = 14
shiftsource = convsource.regrid(rgrid[:-1], rgrid[1:])
ui.set_source(shiftsource)
Documentation relies on the XSPEC documentation, but some additions
have been made (e.g. the need for a fixed normalisation/amplitude
parameter for the flux/lumin models).
There was a section that had
if XSPEC 12.9.1
...
if XSPEC > 12.9.1
...
endif
endif
I have moved them so they do not overlap.
DougBurke
force-pushed
the
add-xspec-convolution-api
branch
from
a93b6eb
to
f1c1627
Compare
Change some links to match the 12.11.0 version (changed since 12.10.1). I have asked team XSPEC if we can have stable URLs for these (they used to exist but disappeared recently). Fixed some URLs: - http to https - add gsfc to host name (the old version no-longer works)
|
@dtnguyen2 - I've finished rebasing / fixing things for the XSPEC 12.11.0 changes which were recently in #772 |
|
@Marie-Terrell is in the process of building xspec 12.11.0 which is needed to test this PR, I will resume the review this PR once the build is complete. |
|
@dtnguyen2 - just to note that this can be run against XSPEC 12.10.1. There's only one convolution model that would be excluded by using 12.10.1 (and thanks to the code that @olaurino wrote ages ago, it will still work woith 12.10.1, it'll just error out at runttime if you try to use this model. |
|
@DougBurke, the first example that you wrote, doesn't work with 12.10.1n |
|
@dtnguyen2 - it should do. What was the error you got? |
This just removes the use of SherpaTest and changes things to use pytest or basic asserts, rather than the SherpaTest functionality. The tests are the same.
|
|
I can't replicate this: The difference between patch m and n here is not relevant. |
The XSPEC team have restored the additive/multiplicative/convolution links so we can refer to them using stable URLS. This also makes some "old" updates just to add documentation for some 12.9.1 and later models (not a significant change).
|
|
It looks to me like a build or environment issue. If you run the following it will point you to the init.py file You can then see if it contains the string Convolution - e.g. |
|
I had to merge a bunch of files to get it to work. How in the world did the sherpa command cause Box1D is defined as: |
|
That's just a typo - it should have been box.xlow.
You shouldn't need to be merging files just to check out a PR. This one has
seen several rebase runs, so you do need to check for these changing, but
if you are regularly needing rebases then there's a problem in your setup.
Doug
…On Thu, Jul 23, 2020 at 1:20 PM dtnguyen2 ***@***.***> wrote:
I had to merge a bunch of files to get it to work. How in the world did
the sherpa command box.xlo = 10 work for you?
Traceback (most recent call last):
File "tmp.py", line 28, in <module>
box.xlo = 10
File "/export/sherpa/sherpa/models/model.py", line 199, in __setattr__
NoNewAttributesAfterInit.__setattr__(self, name, val)
File "/export/sherpa/sherpa/utils/__init__.py", line 146, in __setattr__
(type(self).__name__, name))
AttributeError: 'Box1D' object has no attribute 'xlo'
cause Box1D is defined as:
def __init__(self, name='box1d'):
self.xlow = Parameter(name, 'xlow', 0)
self.xhi = Parameter(name, 'xhi', 0)
self.ampl = Parameter(name, 'ampl', 1, -1, 1)
ArithmeticModel.__init__(self, name, (self.xlow, self.xhi, self.ampl))
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#750 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/AABW27TGRZJ6RBCSEBLTEN3R5BWMBANCNFSM4LATUOQA>
.
|
|
I rebased one of my PRs and it messed up my setup, will fix it soon. |
|
Yay - thanks @dtnguyen2 and @wmclaugh |
Summary
Add support for XSPEC convolution-style models - this link is valid for XSPEC 12.10.1 documentation, it may need changing once the next version is released - such as
cfluxandreflect.Fixes #68
Notes
Since these models change the results of other models they are applied to models rather than combined with the +, *, -, or / operators. That is, you say
xscflux.cfl(mdlexpr)to create a model instance calledcflwhich is applied to the model stored inmdlexpr(which could be a single component or multiple components). You can not use a model instance (e.g.cfl) directly - that iscfl * xspowerlaw.plwill not work.Examples
The following were run with
Previous versions used version 12.10.1
Example: cflux
The
cfluxmodel is a "hack" that lets a user find out the integrated flux for a model (over the emin to emax range of thecfluxcomponent) via the modelslg10Fluxparameter. It can be applied to multiple components but here I show it calculating the unabsorbed flux for a powerlaw fit. Note that you need to freeze the amplitude of the model expression it convolves (otherwise the parameters are degenerate and the optimiser will have a fit if you try to fit).In the following I get the unabsorbed flux for the 0.5-7 keV range from
cfluxand then from Sherpa'scalc_energy_fluxcommand. I haven't tried any of the other "sherpa flux" functionality, which might be interesting to look at how the errors compare, but that is outside of this PR, as I just want to show thatcfluxis behaving sensibly here.which can be compared to the Sherpa fit without the component:
We can see that
nHandgammaare essentially the same, and the reduced statistic is the same (0.567), which you would hope ifcfluxis working correctly. We can compare the flux tooSo they agree to the third decimal place (after rounding).
Example: zashift
The
zashiftmodel redshifts an "additive" style model. Here we use it to redshift a flat background plus box (10-14 keV) from a redshift of 1 to the observed band. I chose these limits since the RMF cuts off at 11 keV, so we can see the need for the "regrid" support (and check how this works).We can see what these look like (manually using matplotlib to create two plots since
ui.plotdoes not handle PHA plots well, which is an issue somewhere but I'm not going to look for it as it's irrelevant to this PR):You can see that the box model has been redshifted from 10-14 keV to 5-7 keV (the background is a constant so it appears unshifted), except that the RMF cut off at 11 keV means that there's no signal above 11/(1+z)=5.5 keV. This is not an issue with this PR, but is my original driver for the regrid code.
Example: zshift + regrid
Can we use the regrid support to extend the array over which the convolution is performed? To avoid known issues with the current regrid support when there's partial overlap of the grid (which are being worked on) I am going to make sure the "regrid" grid covers the RMF grid (which is 0.1 to 11 keV with a spacing of 0.01 keV).
I think this example - namely how
regridis used to create a new model - has a certain amount of logic behind it, but I am not convinced it is the best API for users of the UI layer, but this is related to theregridsupport and is not part of this PR (i.e. I don't want bike-shedding on this part of Sherpa to slow down this PR :-)After this we can see that the source is evaluated "correctly", in that the box is defined over the 5-7 keV range, and the background is still there at higher energies (this is not meant to be a realistic model, which is why the
plot_modeloutput looks strange):What happens when the models are used incorrectly
That is, how useful are the error messages? I am going to try and compare to existing behavior where possible.
Do not supply an instance name
We see the same error as with the existing models:
and
Do not include a model instance
The following is incorrect, since
fluxmodelhas no model expression to work with, but there's no error when we callset_source. Instead, we see the error when the model is evaluated (e.g.plot_source,calc_stat,fit, ...), and unfortunately they are different, and neither are particularly useful.It looks like more-complicated incorrect expressions such as
return the "
numpy.ndarray" type error even forui.plot_source.Details
Sherpa contains "C level" bindings to the XSPEC convolution models, but we have not provided "python level" access to these from Sherpa. I have maintained such an interface in the SDS contributed code, but now that the regrid support is in (needed to use a number of XSPEC convolution models), it makes sense to move the code into Sherpa.
Originally I had added a number of "load" routines (based on
load_psf), one for each convolution model, to create an instance of the model, but in this PR I have decided not to provide this interface, as it isn't needed. There's nothing in this PR that means we couldn't add in such routines (or a single load routine that takes the convolution-model name as an argument, for instance), if we think it's useful (but I'd rather not as I don't think it helps users and just gives us more things to document/test and increases the API).Each model is an instance of the
XSConvolutionKernelclass, and when applied to other models it creates aXSConvolutionModelinstance. I based this on how the PSF convolution kernel is set up. It works, but is an area of the design that needs checking (ie that it makes sense). The code changes are quite small however, which suggests that it fits into Sherpa's "design".