Error during grid.show_dashboard

[itroitskaya]

Hi Gully, I am a PhD student working with @benjaminpope the The University of Queensland. I am trying to use interactive fitting with PHOENIX dashboard and get this error during grid.show_dashboard. Could you please have a look?

from gollum.phoenix import PHOENIXGrid
from specutils import Spectrum1D
import pandas as pd
import astropy.units as u
import numpy as np
from readmultispec import readmultispec

from IPython.display import HTML
from IPython.display import Image

data = readmultispec('./98Tau_2023-02-10_04h48m13s_cb.spec.fits')
data_path = '~/Dev/HiResFITS/PHOENIX-ACES-AGSS-COND-2011'

bdss_spectrum = Spectrum1D(spectral_axis=data['wavelen'][0]*u.micron,
                           flux=data['flux'][0]*u.ct)

wl_lo, wl_hi = (bdss_spectrum.wavelength.value.min(),
                bdss_spectrum.wavelength.value.max())

grid = PHOENIXGrid(wl_lo=wl_lo, wl_hi=wl_hi, path=data_path)

grid.show_dashboard(data=bdss_spectrum)

ERROR:tornado.application:Uncaught exception GET /autoload.js?bokeh-autoload-element=p1002&bokeh-absolute-url=http://localhost:54194&resources=none (::1)
HTTPServerRequest(protocol='http', host='localhost:54194', method='GET', uri='/autoload.js?bokeh-autoload-element=p1002&bokeh-absolute-url=http://localhost:54194&resources=none', version='HTTP/1.1', remote_ip='::1')
Traceback (most recent call last):
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/tornado/web.py", line 1786, in _execute
    result = await result
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/bokeh/server/views/autoload_js_handler.py", line 62, in get
    session = await self.get_session()
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/bokeh/server/views/session_handler.py", line 145, in get_session
    session = await self.application_context.create_session_if_needed(session_id, self.request, token)
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/bokeh/server/contexts.py", line 242, in create_session_if_needed
    self._application.initialize_document(doc)
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/bokeh/application/application.py", line 192, in initialize_document
    h.modify_document(doc)
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/bokeh/application/handlers/function.py", line 142, in modify_document
    self._func(doc)
  File "/Users/uqitroit/Dev/gollum/src/gollum/phoenix.py", line 302, in create_interact_ui
    scalar_norm = np.percentile(self[0].flux.value, 95)
  File "<__array_function__ internals>", line 200, in percentile
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/numpy/lib/function_base.py", line 4205, in percentile
    return _quantile_unchecked(
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/numpy/lib/function_base.py", line 4473, in _quantile_unchecked
    return _ureduce(a,
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3752, in _ureduce
    r = func(a, **kwargs)
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/numpy/lib/function_base.py", line 4639, in _quantile_ureduce_func
    result = _quantile(arr,
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/numpy/lib/function_base.py", line 4745, in _quantile
    take(arr, indices=-1, axis=DATA_AXIS)
  File "<__array_function__ internals>", line 200, in take
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 190, in take
    return _wrapfunc(a, 'take', indices, axis=axis, out=out, mode=mode)
  File "/Users/uqitroit/anaconda3/envs/gollum_dev/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 57, in _wrapfunc
    return bound(*args, **kwds)
IndexError: cannot do a non-empty take from an empty axes.

[Sujay-Shankar]

Hi there! I believe the error you've run into may be stemming from NaN values somewhere in the flux array. As a quick test, in the gollum source code, try going to phoenix.py line 302, and try changing np.percentile to np.nanpercentile, and see if that fixes the issue.

[itroitskaya]

Thanks. I've checked the data, and it turned out that I've mixed up the wavelength units and PHOENIXSpectrum returned empty flux and wavelength arrays with these wrong limits. It may be beneficial in phoenix.py:192 to check not only grid_points, but wavelengths and fluxes as well. However, I've encountered additional problems. First, there was an "Data should overlap the models, double check your wavelength limits." assert fail. The limits were wl_lo = 3850.128, wl_hi = 3925.986 and new_lo: 3850.122160382083, new_hi: 3925.9863831754383 As these were taken directly from the data and Spectrum1D object, they look like a rounding error. What's the best way to deal with it? Next, I've commented out the assertion, and got the following error

File "/Users/uqitroit/Dev/gollum/src/gollum/phoenix.py", line 343, in create_interact_ui
    legend_label=data.meta["header"]["OBJECT"],
KeyError: 'header'

Looks like it is trying to access the header of the Spectrum1D object, which does not have one. What can I do next?

[Sujay-Shankar]

With regard to the first assertion fail, the PHOENIX dashboard requires the following ordering of wavelength limits: wl_lo < new_lo < new_hi < wl_hi, meaning that your current limits just barely throw the error. The assert message may not have communicated this requirement clearly enough. For the header error, currently the dashboard assumes that the data spectrum has an object name stored in its metadata. In hindsight, this is a little brittle.

I've pushed a patch that makes the error message clearer, and no longer breaks when there is no object name in metadata. You will have to either shrink your object's wavelength range or expand the grid's wavelength range, to satisfy the above inequality, and things should work then.

(Edit) Sorry I forgot to address this one initially! We are aware of the inconsistencies with inputting units other than Angstroms into the PHOENIXSpectrum/PHOENIXGrid constructors, and are working on those changes. Until the update for this has been implemented, we ask that you input all values as Angstroms.

If you run into anything else, please let us know!

[itroitskaya]

Thanks for the update. Regarding the assert, I am a bit confused as to how the limits are being set and how should I fix my code for it to work. So currently I have the following code. Spectrum is generated from data, and grid is generated from spectrum. I don't manually set the limits anywhere.

bdss_spectrum = Spectrum1D(spectral_axis=data['wavelen'][0]*u.angstrom, flux=data['flux'][0]*u.ct)
wl_lo, wl_hi = (bdss_spectrum.wavelength.value.min(), bdss_spectrum.wavelength.value.max())
grid = PHOENIXGrid(wl_lo=wl_lo, wl_hi=wl_hi, path=data_path)

print(f'data: {data["wavelen"].min()} - {data["wavelen"].max()}')
print(f'spectrum: {bdss_spectrum.wavelength.value.min()} - {bdss_spectrum.wavelength.value.max()}')
print(f'grid: {grid.wavelength.value.min()} - {grid.wavelength.value.max()}')

data: 3850.122160382083 - 3925.9863831754383
spectrum: 3850.122160382083 - 3925.9863831754383
grid: 3850.128 - 3925.986

if I do grid = PHOENIXGrid(wl_lo=wl_lo-1, wl_hi=wl_hi+1, path=phoenix_path), then the dashboard is displayed, but I feel that's maybe not the best solution.

[Sujay-Shankar]

For those limits, the intuition was that since we have rotational broadening functionality, if we had the PHOENIXGrid have the exact same wavelength axis as the data spectrum, you would see edge effects creeping in from the sides as the convolution did its thing, which would not be desirable to users. To prevent this, we enforce the inequalities such the wavelength axis of the data must lie within that of the grid. Providing an ample buffer on either side will hide edge effects such that they do not visually propagate into the region in which the data spectrum lies. In light of that, the last solution you presented is actually probably the way to go, although we recommend a buffer more like 50 Angstroms on either side if I remember correctly. Something like buffer = 50 grid = PHOENIXGrid(wl_lo=wl_lo-buffer, wl_hi=wl_hi+buffer, path=phoenix_path) should work, and is the intended use when dealing with data.

[gully]

Welcome @itroitskaya ! 👋 Thank you for your interest in using gollum! We promise the dashboard is really great and intuitive once you getting working so thank you for your patience in working around some of the friction points and hardcode. @SujayShankarUT's guidance is correct and should get you up-and-running. We value your feedback, and so if you're okay if it we'd love to have your feedback on how we can improve the dashboard. We will use this bug report to make some more usability improvements, so please do report your experiences, either positive or negative, with the code.

Thanks again!