search

desihub / desisim

DESI simulations
BSD 3-Clause "New" or "Revised" License
16 stars 22 forks source link

problem interpolating atmosphere in desisim.specsim.get_simulator #564

Open moustakas opened 2 years ago

moustakas commented 2 years ago 
from desisim.simexp import simulate_spectra
from desisim.templates import ELG
flux, wave, meta, obj = ELG().make_templates(1, seed=1)
sim1 = simulate_spectra(wave, flux)
INFO:iers.py:82:freeze_iers: Freezing IERS table used by astropy time, coordinates.
DEBUG:simexp.py:418:simulate_spectra: loading specsim desi config desi
DEBUG:simexp.py:422:simulate_spectra: creating specsim desi simulator
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-9-f559ab1da86b> in <module>
----> 1 sim1 = simulate_spectra(wave, flux)

~/code/desihub/desisim/py/desisim/simexp.py in simulate_spectra(wave, flux, fibermap, obsconditions, redshift, dwave_out, seed, psfconvolve, specsim_config_file)
    422     log.debug('creating specsim desi simulator')
    423     # desi = specsim.simulator.Simulator(config, num_fibers=nspec)
--> 424     desi = desisim.specsim.get_simulator(config, num_fibers=nspec,
    425         camera_output=psfconvolve)
    426

~/code/desihub/desisim/py/desisim/specsim.py in get_simulator(config, num_fibers, camera_output)
     50         #- New config; create Simulator object
     51         import specsim.simulator
---> 52         qsim = specsim.simulator.Simulator(config, num_fibers,
     53             camera_output=camera_output)
     54

~/code/desihub/specsim/specsim/simulator.py in __init__(self, config, num_fibers, camera_output, verbose)
     72
     73         # Initalize our component models.
---> 74         self.atmosphere = specsim.atmosphere.initialize(config)
     75         self.instrument = specsim.instrument.initialize(config, camera_output)
     76         self.source = specsim.source.initialize(config)

~/code/desihub/specsim/specsim/atmosphere.py in initialize(config)
    709
    710     # Load tabulated data.
--> 711     surface_brightness_dict = config.load_table(
    712         atm_config.sky, 'surface_brightness', as_dict=True)
    713     extinction_coefficient = config.load_table(

~/code/desihub/specsim/specsim/config.py in load_table(self, parent, column_names, interpolate, as_dict)
    486                         kind='linear', copy=False,
    487                         bounds_error=bounds_error, fill_value=fill_value)
--> 488                     interpolated_values = interpolator(self.wavelength.value)
    489                     unit = loaded_columns[column_name].unit
    490                     if unit:

~/opt/anaconda3/envs/desi/lib/python3.8/site-packages/scipy/interpolate/polyint.py in __call__(self, x)
     76         """
     77         x, x_shape = self._prepare_x(x)
---> 78         y = self._evaluate(x)
     79         return self._finish_y(y, x_shape)
     80

~/opt/anaconda3/envs/desi/lib/python3.8/site-packages/scipy/interpolate/interpolate.py in _evaluate(self, x_new)
    682         y_new = self._call(self, x_new)
    683         if not self._extrapolate:
--> 684             below_bounds, above_bounds = self._check_bounds(x_new)
    685             if len(y_new) > 0:
    686                 # Note fill_value must be broadcast up to the proper size

~/opt/anaconda3/envs/desi/lib/python3.8/site-packages/scipy/interpolate/interpolate.py in _check_bounds(self, x_new)
    714                              "range.")
    715         if self.bounds_error and above_bounds.any():
--> 716             raise ValueError("A value in x_new is above the interpolation "
    717                              "range.")
    718

ValueError: A value in x_new is above the interpolation range.
moustakas commented 2 years ago

This ticket is a bit of a rabbit hole and I'm inclined to close it. The crash actually occurs in specsim when the atmospheric extinction is interpolated onto the input wavelength grid-- https://github.com/desihub/specsim/blob/master/specsim/config.py#L488

However, the extinction grid only extends to 9999.9 A while the wavelength grid in the example in this ticket extends to 10000.0 A and so an extrapolation exception is raised. In its defense, specsim does have some options to allow for extrapolated interpolation, but it looks like some of the options in scipy may have changed since the code was written.

Unfortunately, choosing a shorter maximum wavelength results in a different crash that I'm not inspired to pursue:

from desisim.simexp import simulate_spectra
from desisim.templates import ELG
flux, wave, meta, obj = ELG(maxwave=9990).make_templates(1, seed=1)
sim1 = simulate_spectra(wave, flux)
INFO:io.py:971:read_basis_templates: Reading /Users/ioannis/work/desi/spectro/templates/basis_templates/v3.2/elg_templates_v2.2.fits
INFO:iers.py:82:freeze_iers: Freezing IERS table used by astropy time, coordinates.
DEBUG:simexp.py:418:simulate_spectra: loading specsim desi config desi
DEBUG:simexp.py:422:simulate_spectra: creating specsim desi simulator
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/Users/ioannis/code/desihub/desisim/py/desisim/simexp.py", line 423, in simulate_spectra
    desi = desisim.specsim.get_simulator(config, num_fibers=nspec,
  File "/Users/ioannis/code/desihub/desisim/py/desisim/specsim.py", line 52, in get_simulator
    qsim = specsim.simulator.Simulator(config, num_fibers,
  File "/Users/ioannis/code/desihub/specsim/specsim/simulator.py", line 75, in __init__
    self.instrument = specsim.instrument.initialize(config, camera_output)
  File "/Users/ioannis/code/desihub/specsim/specsim/instrument.py", line 622, in initialize
    initialized_cameras.append(specsim.camera.Camera(
  File "/Users/ioannis/code/desihub/specsim/specsim/camera.py", line 145, in __init__
    raise RuntimeError(
RuntimeError: Wavelength grid min does not cover b-camera response.

For the record, this ticket came up as I was trying to understand https://github.com/desihub/desisim/issues/563.