Fitting custom spectrum with energy range below 1 keV with thermal (f_vth) model.

[abhilash-sw]

Describe the bug

I am trying to fit a custom spectrum which has spectral data below 1 keV (though not usable). I am trying to fit the spectra with above 1 keV using spec.energy_fitting_range attribute.

I am getting "ValueError: All input energy values must be within the range 1.0002920302956426--200.15819869050395 keV."

Subsequently I tried to limit my input spectrum above 1 keV. I am faced with IndexError (full error attached below.)

To Reproduce

spec_data_dict = {'count_channel_bins':channel_bins, 'counts': counts, 'count_error':counts_err, 'effective_exposure':exposure, 'srm': rsp_mat}

spec = Fitter(spec_data_dict)

spec.loglikelihood="poisson"
spec.model = "C*f_vth"

fiter=[1.83,15.0]

spec.energy_fitting_range = fiter
spec.params["C_spectrum1"] = {"Status":"frozen"}

spec.params["T1_spectrum1"] = {"Value":10.0, "Bounds":(5.0, 25.0)} #Plasma temperature in megakelvin
spec.params["EM1_spectrum1"] = {"Value":1e-1, "Bounds":(0.01, 1.0)} #Emission measure in units of 1e46 cm^-3.

print(spec.params)

minimised_params = spec.fit()

Screenshots

---

IndexError Traceback (most recent call last)
Cell In[41], line 1
----> 1 minimised_params = spec.fit()

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/sunxspex_fitting/fitter.py:2135, in Fitter.fit(self, **kwargs)
2127 kwargs.pop("_abs_hess_step", None)
2129 # this has been replaced by a _run_minimiser_core so that this can be swapped out easily down the line
2130 # soltn = minimize(self._fit_stat_minimize,
2131 # free_params_list,
2132 # args=stat_args,
2133 # bounds=free_bounds,
2134 # **kwargs)
-> 2135 soltn = self._run_minimiser_core(minimise_func=self._fit_stat_minimize,
2136 free_parameter_list=free_params_list,
2137 statistic_args=stat_args,
2138 free_param_bounds=free_bounds,
2139 **kwargs)
2141 self._minimize_solution = soltn
2143 std_err = self._calc_minimize_error(stat_args, step=_hess_step, _abs_step=_abs_hess_step)

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/sunxspex_fitting/fitter.py:2097, in Fitter._run_minimiser_core(self, minimise_func, free_parameter_list, statistic_args, free_param_bounds, **kwargs)
2063 def _run_minimiser_core(self, minimise_func, free_parameter_list, statistic_args, free_param_bounds, **kwargs):
2064 """ Allows user (or us) to define their own (different) minimiser easily.
2065
2066 This should return the same type of output as Scipy's minimize at the minute. This just passes
(...)
2095 the parameter results in the order of free_parameter_list.
2096 """
-> 2097 return minimize(minimise_func, free_parameter_list, args=statistic_args, bounds=free_param_bounds, **kwargs)

File ~/anaconda3/lib/python3.10/site-packages/scipy/optimize/_minimize.py:684, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
681 bounds = standardize_bounds(bounds, x0, meth)
683 if meth == 'nelder-mead':
--> 684 res = _minimize_neldermead(fun, x0, args, callback, bounds=bounds,
685 **options)
686 elif meth == 'powell':
687 res = _minimize_powell(fun, x0, args, callback, bounds, **options)

File ~/anaconda3/lib/python3.10/site-packages/scipy/optimize/_optimize.py:845, in _minimize_neldermead(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, bounds, **unknown_options)
843 try:
844 for k in range(N + 1):
--> 845 fsim[k] = func(sim[k])
846 except _MaxFuncCallError:
847 pass

File ~/anaconda3/lib/python3.10/site-packages/scipy/optimize/_optimize.py:569, in _wrap_scalar_function_maxfun_validation..function_wrapper(x, *wrapper_args)
567 ncalls[0] += 1
568 # A copy of x is sent to the user function (gh13740)
--> 569 fx = function(np.copy(x), *(wrapper_args + args))
570 # Ideally, we'd like to a have a true scalar returned from f(x). For
571 # backwards-compatibility, also allow np.array([1.3]),
572 # np.array([[1.3]]) etc.
573 if not np.isscalar(fx):

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/sunxspex_fitting/fitter.py:1923, in Fitter._fit_stat_minimize(self, *args, **kwargs)
1909 def _fit_stat_minimize(self, *args, **kwargs):
1910 """ Return the chosen fit statistic defined to minimise for the best fit.
1911
1912 I.e., returns -2ln(L).
(...)
1921 Float.
1922 """
-> 1923 return self._fit_stat(*args, maximize_or_minimize="minimize", **kwargs)

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/sunxspex_fitting/fitter.py:1661, in Fitter._fit_stat(self, free_params_list, photon_channels, count_channel_mids, srm, livetime, ph_e_binning, observed_counts, observed_count_errors, tied_or_frozen_params_list, param_name_list_order, maximize_or_minimize, **kwargs)
1606 """ Calculate the fit statistic from given parameter values.
1607
1608 Calculates the model count spectrum and calculates the fit statistic
(...)
1657 Float, the combined ln(L) or -2ln(L) for all spectra loaded and model.
1658 """
1660 # make sure only the free parameters are getting varied so put them first
-> 1661 mu = self._pseudo_model(free_params_list,
1662 tied_or_frozen_params_list,
1663 param_name_list_order,
1664 photon_channels=photon_channels,
1665 photon_channel_widths=ph_e_binning,
1666 count_channel_mids=count_channel_mids,
1667 total_responses=srm,
1668 **kwargs)
1670 ll = 0
1671 for m, o, l, err in zip(mu, observed_counts, livetime, observed_count_errors):
1672
1673 # calculate the count rate model for each spectrum

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/sunxspex_fitting/fitter.py:1431, in Fitter._pseudo_model(self, free_params_list, tied_or_frozen_params_list, param_name_list_order, **other_inputs)
1428 # change the tied parameter's value to the value of the param it is tied to
1429 dictionary = self._tie_params(dictionary)
-> 1431 return self._counts_model(**dictionary, **other_inputs)

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/sunxspex_fitting/fitter.py:1367, in Fitter._counts_model(self, **kwargs)
1364 sep_params = dict(zip(self._orig_params, ordered_kwarg_values))
1366 # calculate the [photon s^-1 cm^-2]
-> 1367 m = self._model(**sep_params, energies=kwargs["photon_channels"][s]) * kwargs["photon_channel_widths"][s] # np.diff(kwargs["photon_channels"][s]).flatten() # remove energy bin dependence
1369 # fold the photon model through the SRM to create the count rate model, [photon s^-1 cm^-2] * [count photon^-1 cm^2] = [count s^-1]
1370 cts_model = make_model(energies=kwargs["photon_channels"][s],
1371 photon_model=m,
1372 parameters=None,
1373 srm=kwargs["total_responses"][s])

File :2, in C_f_vthT1EM1C(T1, EM1, C, energies)

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/sunxspex_fitting/photon_models_for_fitting.py:50, in f_vth(temperature, emission_measure46, energies)
48 temperature = temperature * 1e6 << u.K
49 emission_measure = emission_measure46 * 1e46 << u.cm**(-3)
---> 50 return thermal_emission(energies, temperature, emission_measure).value

File ~/anaconda3/lib/python3.10/site-packages/astropy/units/decorators.py:313, in QuantityInput.call..wrapper(*func_args, **func_kwargs)
311 # Call the original function with any equivalencies in force.
312 with add_enabled_equivalencies(self.equivalencies):
--> 313 return_ = wrapped_function(*func_args, **func_kwargs)
315 # Return
316 ra = wrapped_signature.return_annotation

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/thermal.py:221, in thermal_emission(energy_edges, temperature, emission_measure, abundance_type, relative_abundances, observer_distance)
219 # Calculate fluxes.
220 continuum_flux = _continuum_emission(energy_edges_keV, temperature_K, abundances)
--> 221 line_flux = _line_emission(energy_edges_keV, temperature_K, abundances)
222 flux = ((continuum_flux + line_flux) * emission_measure /
223 (4 * np.pi * observer_distance**2))
224 if temperature.isscalar and emission_measure.isscalar:

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/thermal.py:438, in _line_emission(energy_edges_keV, temperature_K, abundances)
431 ##### Weight Line Emission So Peak Energies Maintained Within Input Energy Binning #####
432 # Split emission of each line between nearest neighboring spectral bins in
433 # proportion such that the line centroids appear at the correct energy
434 # when averaged over neighboring bins.
435 # This has the effect of appearing to double the number of lines as regards
436 # the dimensionality of the line_intensities array.
437 line_peaks_keV = LINE_GRID["line peaks keV"][energy_roi_indices]
--> 438 split_line_intensities, line_spectrum_bins = _weight_emission_bins_to_line_centroid(
439 line_peaks_keV, energy_edges_keV, line_intensities)
441 #### Calculate Flux #####
442 # Use binned_statistic to determine which spectral bins contain
443 # components of line emission and sum over those line components
444 # to get the total emission is each spectral bin.
445 flux = stats.binned_statistic(line_spectrum_bins, split_line_intensities,
446 "sum", n_energy_bins, (0, n_energy_bins-1)).statistic

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/thermal.py:626, in _weight_emission_bins_to_line_centroid(line_peaks_keV, energy_edges_keV, line_intensities)
624 new_iline = copy.deepcopy(iline)
625 if len(neg_deviation_indices) > 0:
--> 626 neg_line_intensities, neg_neighbor_intensities, neg_neighbor_iline = _weight_emission_bins(
627 line_deviations_keV, neg_deviation_indices,
628 energy_center_diffs, line_intensities, iline, negative_deviations=True)
629 # Combine new line and neighboring bin intensities and indices into common arrays.
630 new_line_intensities[:, neg_deviation_indices] = neg_line_intensities

File ~/anaconda3/lib/python3.10/site-packages/sunkit_spex-0.1.dev156+gcfd41b9-py3.10.egg/sunkit_spex/thermal.py:677, in _weight_emission_bins(line_deviations_keV, deviation_indices, energy_center_diffs, line_intensities, iline, negative_deviations)
672 b = 1
674 # Calculate difference between line energy and the spectrum bin center as a
675 # fraction of the distance between the spectrum bin center and the
676 # center of the nearest neighboring bin.
--> 677 wghts = np.absolute(line_deviations_keV[deviation_indices]) / energy_center_diffs[iline[deviation_indices+a]]
678 # Tile/replicate wghts through the other dimension of line_intensities.
679 wghts = np.tile(wghts, tuple([line_intensities.shape[0]] + [1] * wghts.ndim))

IndexError: index 58 is out of bounds for axis 0 with size 58

System Details

No response

Installation method

No response

[ianan]

The lookup files that the thermal model uses were originally created for ospex in sswidl using CHIANTI and cover 1-250 keV - they are the v71 files here https://hesperia.gsfc.nasa.gov/ssw/packages/xray/dbase/chianti/

So you would have to use other files, that directory contains v901 ones that go down to either 0.1 or 0.7 keV (01_12 or 07_12 files), though there are still other issues at the moment - see issue #99 for more info.

[abhilash-sw]

Thank you for the suggestion.

I tried to use v901 files (chianti_lines_01_12_unity_v901_t41.sav and chianti_cont_01_250_unity_v901.geny).
But the result doesn't seem to be accurate (attaching generated spectrum for 10MK). Following is the code I used to generate the spectrum.

from sunkit_spex import thermal
thermal.CONTINUUM_GRID=thermal.setup_continuum_parameters('idl/chianti_cont_01_250_unity_v901.geny')
thermal.LINE_GRID=thermal.setup_line_parameters('idl/chianti_lines_01_12_unity_v901_t41.sav')
engs=np.arange(0.4,21,0.01)
f10=thermal.thermal_emission(engs << u.keV,10*1e6<<u.K,1e49<< u.cm**(-3))

The continuum file is in geny format with different keywords. I had to modify load_chianti_continuum function in io.

def load_chianti_continuum():
    """
    Read X-ray continuum emission info from an IDL sav file produced by CHIANTI

    Returns
    -------
    continuum_intensities: `xarray.DataArray`
        Continuum intensity as a function of element, temperature and energy/wavelength
        and associated metadata and coordinates.

    Notes
    -----
    By default, this function uses the file located at
    https://hesperia.gsfc.nasa.gov/ssw/packages/xray/dbase/chianti/chianti_cont_1_250_v71.sav
    To use a different file call this function in the following way:
    >>> from sunpy.data import manager # doctest: +SKIP
    >>> with manager.override_file("chianti_continuum", uri=filename): # doctest: +SKIP
    ...    line_info = load_chianti_lines_light() # doctest: +SKIP

    where filename is the location of the file to be read.
    """
    # Define units
    intensity_unit = u.ph * u.cm**3 / (u.s * u.keV * u.sr)
    temperature_unit = u.K
    wave_unit = u.AA
    # Read file
    contfile = manager.get("chianti_continuum")
    contents = scipy.io.readsav(contfile)
    # Concatenate low and high wavelength intensity arrays.
    if contfile.suffix=='.geny':
        totcont_lo_key = 'p1'
        totcont_key = 'p2'
        zindex_key = 'p0'
        ctemp_key = 'p4'
        edge_str_key = 'p3'
        chianti_doc_key = 'p5'
    else:
        totcont_lo_key = 'totcont_lo'
        totcont_key = 'totcont'
        zindex_key = 'zindex'
        ctemp_key = 'ctemp'
        edge_str_key = 'edge_str'
        chianti_doc_key = 'chianti_doc'

    intensities = np.concatenate((contents[totcont_lo_key], contents[totcont_key]), axis=-1)
    # Integrate over sphere surface of radius equal to observer distance
    # to get intensity at source. This means that physical intensities can
    # be calculated by dividing by 4 * pi * R**2 where R is the observer distance.
    intensities *= 4 * np.pi
    intensity_unit *= u.sr
    # Put file data into intuitive structure and return data.
    continuum_intensities = xarray.DataArray(
        intensities,
        dims=["element_index", "temperature", "wavelength"],
        coords={"element_index": contents[zindex_key],
                "temperature": contents[ctemp_key],
                "wavelength": _clean_array_dims(contents[edge_str_key]["WVL"])},
        attrs={"units": {"data": intensity_unit,
                         "temperature": temperature_unit,
                         "wavelength": wave_unit},
               "file": contfile,
               "wavelength_edges": _clean_array_dims(contents[edge_str_key]["WVLEDGE"]) * wave_unit,
               "chianti_doc": _clean_chianti_doc(contents[chianti_doc_key])
               })
    return continuum_intensities

[ianan]

I've generated new files that go down to 0.1 keV and can be found here: https://gla-my.sharepoint.com/personal/iain_hannah_glasgow_ac_uk/_layouts/15/onedrive.aspx?id=%2Fpersonal%2Fiain%5Fhannah%5Fglasgow%5Fac%5Fuk%2FDocuments%2Fchxdb&ga=1

I tested them in sunkit-spex and seem to be fine (though takes a while to load) https://github.com/ianan/fvth_stuff/blob/main/better_chxdb/check_01files.ipynb