Importing Bugs # 13901 - Seg.Fault with group/ignore about sherpa

[olaurino]

@anetasie could you please include information on how to reproduce this issue? Thanks!

[anetasie]

The data and the script are at /data/lenin2/Bugs/Sherpa/CrashWhenIgnore/bug.py

It seems that this one is related to the issue #38

[olaurino]

Here is a minimal script that reproduces the issue:

from sherpa.astro.ui import *
infile    = "493.dat"
rmffile   = "rmf.fits"
arffile   = "out_1062.arf" 
set_source(xswabs.absrb * xsapec.thrm )
set_method("moncar")
load_data( "493.dat", dstype=DataPHA)
load_arf( arffile)
load_rmf( rmffile)
group_counts(5)
ignore(5.0, None)
absrb.nH = 0.0314
thrm.kT = 0.9
thrm.kT.max = 15
thrm.redshift = 0.009
freeze(absrb.nH)
fit()

Files and scripts to reproduce the issue are here:
https://gist.github.com/olaurino/bf55c41d6c10c7782775

Note that one needs xspec in order to reproduce the original report.

[olaurino]

I can reproduce the issue on OSX 10.9.

The issue seems to be triggered when xspec models are involved, and in particular only when the xsapec model is present in the model expression. For instance, the following code does not segfault:

from sherpa.astro.ui import *
infile    = "493.dat"
rmffile   = "rmf.fits"
arffile   = "out_1062.arf" 
set_source(xswabs.absrb)
set_method("moncar")
load_data( "493.dat", dstype=DataPHA)
load_arf( arffile)
load_rmf( rmffile)
group_counts(5)
ignore(5.0, None)
absrb.nH = 0.0314
fit()

[olaurino]

Since the issue was reported by Jamie on 10.9, I tried the new Mavericks build for CIAO and issue seems to be gone. However, I get a gracefully-handled overflow Warning, and I am wondering if that's the origin of the issue:

$ sherpa bug.py
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.8 Sherpa version 1 Tuesday, June 16, 2015

/proj/xena/ciaot_install/osxm.150616/ciao-4.8/lib/python2.7/site-packages/sherpa/models/model.py:469: RuntimeWarning: overflow encountered in multiply
  val = self.op(lhs, rhs)
Dataset               = 1
Method                = moncar
Statistic             = chi2gehrels
Initial fit statistic = 4.33057
Final fit statistic   = 2.47986e+20 at function evaluation 4144
Data points           = 2
Degrees of freedom    = 0
Reduced statistic     = nan
Change in statistic   = 2.47986e+20
   thrm.kT        3.32009     
   thrm.norm      7.52533e-07 
sherpa-1> 

In any case I am going to build Sherpa from the git sources with xspec support on a 10.9 system to see if the issue shows up or not.

[olaurino]

I currently have issues with building standalone Sherpa with the new OSX Mavericks build of XSPEC. However, as reported in the previous comment, the daily CIAOD build for Mavericks does not segfault.

On the other hand, I don't think the result is correct, and there may still be an underlying memory issue.

@anetasie I need more information in order to make sure that the tests meaningfully pass, whether or not a segfault is triggered. It may be that on Linux there is no segfault but the results do not make any sense, as in the example above on 10.9, where there is no segfault but the results are clearly non-sensical.

Is there something I can check after the fit has been performed? Anything from the expected parameter values to the statistic. I need an answer that make scientific sense and that I can use as an "oracle" of what the result of the test should be. Only this way I can make sure the tests correctly pass.

In the meantime, I'll keep working on the XSPEC build and fix the issues I currently have on Mavericks.

[anetasie]

Grouping and filtering in this test leaves 2 data channels for the fit. This should not be a problem, as there are about 100 individual channels for the model evaluation and the final parameter values seem reasonable. I plot_fit() to visually inspect the results. However, the returned reduced statistics in your example above shows nan - as there are 0 degrees-of-freedom in this case, so the division by 0 in calculating reduce stat. Perhaps this is the source of the problem.

[DougBurke]

I believe there is a problem with the XSpec apec model when called with a non-contiguous grid; this may be related to #42. I have a (not-yet-converted-onto-github) RFE about improving this part of our interface to Xspec models (the RFE is because I don't agree with how it's currently done). Here's an example of the problem (note that it may be tricky to trigger, so I show all the steps):

% ipython --pylab
Python 2.7.10 |Continuum Analytics, Inc.| (default, May 28 2015, 17:02:03) 
Type "copyright", "credits" or "license" for more information.

IPython 3.2.0 -- An enhanced Interactive Python.
Anaconda is brought to you by Continuum Analytics.
Please check out: http://continuum.io/thanks and https://anaconda.org
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.
from sherpa.astro import xspec
froUsing matplotlib backend: Qt4Agg

In [1]: from sherpa.astro import xspec

In [2]: xspec.set_xschatter(25)

In [3]: from sherpa.astro import ui

In [4]: ui.create_model_component('xsapec', 'mdl')

In [5]: print(mdl)
xsapec.mdl
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   mdl.kT       thawed            1        0.008           64        keV
   mdl.Abundanc frozen            1            0            5           
   mdl.redshift frozen            0       -0.999           10           
   mdl.norm     thawed            1            0        1e+24           

In [6]: egrid = np.arange(0.1, 10, 0.01)

In [7]: y1 = mdl(egrid)
SUMDEM : itype = 4

SUMDEM : itype = 4
SUMDEM : TraceAbund = 1

Thermal broadening : F
Velocity broadening : 0
Interpolating between 0.8617 and 1.08482 for 1
with coefficients 0.380143 and 0.619857

Reading data for temperature 30 0.8617 0.861734

Reading data for temperature 30 0.8617 0.861734
Read 29478 element lines and 30  continuum points for 30 unique elements.

Reading data for temperature 31 1.08482 1.08486

Reading data for temperature 31 1.08482 1.08486
Read 26326 element lines and 30  continuum points for 30 unique elements.


In [8]: y2 = mdl(egrid)
SUMDEM : itype = 4

SUMDEM : itype = 4
SUMDEM : TraceAbund = 1

Thermal broadening : F
Velocity broadening : 0
Interpolating between 0.8617 and 1.08482 for 1
with coefficients 0.380143 and 0.619857


In [9]: (y1 == y2).all()
Out[9]: True

In [10]: xlo = egrid[:-1]

In [11]: xhi = egrid[1:]

In [12]: y2 = mdl(xlo, xhi)
SUMDEM : itype = 4

SUMDEM : itype = 4
SUMDEM : TraceAbund = 1

Thermal broadening : F
Velocity broadening : 0
Interpolating between 0.8617 and 1.08482 for 1
with coefficients 0.380143 and 0.619857


In [13]: (y1[:-1] == y2).all()
Out[13]: True

In [14]: mdl([0.1,0.2,0.3,0.4,0.5])
SUMDEM : itype = 4

SUMDEM : itype = 4
SUMDEM : TraceAbund = 1

Thermal broadening : F
Velocity broadening : 0
Interpolating between 0.8617 and 1.08482 for 1
with coefficients 0.380143 and 0.619857

Out[14]: array([ 2.01070237,  0.31376797,  0.21693133,  0.12136491,  0.        ], dtype=float32)

In [15]: mdl([0.1,0.2,0.3,0.4,0.5],[0.2,0.3,0.4,0.5,0.6])
SUMDEM : itype = 4

SUMDEM : itype = 4
SUMDEM : TraceAbund = 1

Thermal broadening : F
Velocity broadening : 0
Interpolating between 0.8617 and 1.08482 for 1
with coefficients 0.380143 and 0.619857

Out[15]: array([ 2.01070237,  0.31376797,  0.21693133,  0.12136491,  0.08586448], dtype=float32)

In [16]: mdl([0.1,0.2,0.4,0.5],[0.2,0.3,0.5,0.6])
SUMDEM : itype = 4

SUMDEM : itype = 4
SUMDEM : TraceAbund = 1

Thermal broadening : F
Velocity broadening : 0
Interpolating between 0.8617 and 1.08482 for 1
with coefficients 0.380143 and 0.619857

SUMDEM : itype = 4

SUMDEM : itype = 4
SUMDEM : TraceAbund = 1

Thermal broadening : F
Velocity broadening : 0
Interpolating between 0.8617 and 1.08482 for 1
with coefficients 0.380143 and 0.619857

Segmentation fault (core dumped)

This is with

% python -c 'import sherpa; import sherpa.astro.xspec; print(sherpa.__version__); print(sherpa.astro.xspec.get_xsversion())'
4.7+475.g0ee48e5.dirty
12.8.2q

EDIT: this is on a linux machine

[DougBurke]

So, apec (and a lot of other models) can crash if called with only two bins, as shown below. I think this is what causes the crash in my previous report, since the way the current Xspec interface code deals with non-continuous grids is to run the model several times, for the different grids (needs checking, as this is going of memory). I don't know if this is related to the issues @olaurino is seeing above.

In [1]: from sherpa.astro import xspec

In [2]: apec = xspec.XSapec()

In [3]: apec([0.1,0.2])
Segmentation fault (core dumped)

EDIT: I've added this as a separate issue (#62)
EDIT: investigating the code above in #55 (comment), I do not think that the crash above is related to #62, as the model is being evaluated on a large grid (~900 bins)

[DougBurke]

I have also checked the ARF and RMF to make sure that their grids do not contain any gaps (which would trigger #62, as discussed in #56 (comment) ) but they do not seem to:

CIAO 4.7 Sherpa version 1 Thursday, December 4, 2014

sherpa-1> x = np.arange(1, 1025, 1)
sherpa-2> y = np.ones(x.size)
sherpa-3> load_arrays(1, x, y, DataPHA)
sherpa-4> load_rmf('rmf.fits')
sherpa-5> set_source(powlaw1d.mdl1)
sherpa-6> m = get_model()
sherpa-7> xlo = m.xlo
sherpa-8> xhi = m.xhi
sherpa-9> (xlo < xhi).all()
          True
sherpa-10> (xlo[1:] == xhi[:-1]).all()
           True
sherpa-11> load_arf('out_1062.arf')
sherpa-12> m2 = get_model()
sherpa-13> alo = m2.xlo
sherpa-14> ahi = m2.xhi
sherpa-15> (alo == xlo).all()
           True
sherpa-16> (ahi == xhi).all()
           True

[DougBurke]

This should be fixed by #101 (and going to XSPEC 12.9.0b or later) but I'll leave this open until we can get confirmation that it's no longer a problem on OS-X 10.9

[olaurino]

I can confirm that the issue does not show up for the current CIAOT, i.e. with code from the release/4.8.0 branch built in CIAO on November 23rd, 2015.
Shall I close this issue?

[olaurino]

This issue can probably be closed. Right?

[DougBurke]

I am happy for this to be closed.