search

AlexandrosAntoniadis / ODUSSEAS

MIT License
6 stars 1 forks source link

Interpolation error found for one spectrum that I tried #11

Closed sousasag closed 1 month ago

sousasag commented 2 months ago

The error is below. I patch this on my instalation allowing extrapolation in the interp1d function:

f2 = interp1d(self.wavelength, self.flux, kind="linear", fill_value='extrapolate')

f2 = interp1d(wave_rv, self.flux, kind="linear", fill_value='extrapolate')

 /home/sousasag/Programas/anaconda3/envs/saitama/lib/python3.11/site-packages/odusseas/utils.py:104 in __post_init__                                                 │
│                                                                                                                                                                     │
│   101 │   │   if round(dw, 4) != 0.010:                                                                                                                             │
│   102 │   │   │   f2 = interp1d(self.wavelength, self.flux, kind="linear")                                                                                          │
│   103 │   │   │   self.wavelength = np.arange(self.wavelength[0], self.wavelength[-1], cdelt1)                                                                      │
│ ❱ 104 │   │   │   self.flux = f2(self.wavelength)                                                                                                                   │
│   105 │   │                                                                                                                                                         │
│   106 │   │   rv = find_rv(self.wavelength, self.flux)                                                                                                              │
│   107 │   │   wave_rv = self.wavelength / (1 + rv / 3.0e5) if abs(rv) > 0 else self.wavelength                                                                      │
│                                                                                                                                                                     │
│ ╭────────────────────────────────────────────────────────────── locals ───────────────────────────────────────────────────────────────╮                             │
│ │ cdelt1 = 0.01                                                                                                                       │                             │
│ │     dw = 0.005                                                                                                                      │                             │
│ │     f2 = <scipy.interpolate._interpolate.interp1d object at 0x783107741360>                                                         │                             │
│ │    hdr = SIMPLE  =                    T / conforms to FITS standard                                                                 │                             │
│ │          BITPIX  =                  -64 / array data type                                                                           │                             │
│ │          NAXIS   =                    1 / number of array dimensions                                                                │                             │
│ │          NAXIS1  =               824989                                                                                             │                             │
│ │          CDELT1  =                0.005                                                                                             │                             │
│ │          CRVAL1  =    3771.560370205989                                                                                             │                             │
│ │          CTYPE1  = 'Angstrom'                                                                                                       │                             │
│ │          MERGED  = 'Merged with python Astropy combine_espresso_spectra.py'                                                         │                             │
│ │          SPIKES  = 'UVES Spikes removed with sigclop aka RemoveSpikes, sigma=3'                                                     │                             │
│ │          FILE000 = 'ADP.2022-11-04T11:55:39.934.fits'                                                                               │                             │
│ │          FILE001 = 'ADP.2022-12-05T10:14:29.458.fits'                                                                               │                             │
│ │          FILE002 = 'ADP.2022-12-05T11:56:17.792.fits'                                                                               │                             │
│ │          FILE003 = 'ADP.2022-11-04T11:18:13.623.fits'                                                                               │                             │
│ │          FILE004 = 'ADP.2022-11-04T12:31:16.427.fits'                                                                               │                             │
│ │          FILE005 = 'ADP.2022-11-04T11:29:07.002.fits'                                                                               │                             │
│ │          FILE006 = 'ADP.2022-11-04T11:50:15.727.fits'                                                                               │                             │
│ │          FILE007 = 'ADP.2022-11-04T11:44:56.979.fits'                                                                               │                             │
│ │          FILE008 = 'ADP.2022-11-04T11:39:04.070.fits'                                                                               │                             │
│ │          FILE009 = 'ADP.2022-11-04T11:25:46.995.fits'                                                                               │                             │
│ │          FILE010 = 'ADP.2022-12-05T12:01:03.894.fits'                                                                               │                             │
│ │          FILE011 = 'ADP.2022-11-04T12:26:46.784.fits'                                                                               │                             │
│ │          FILE012 = 'ADP.2022-11-04T12:06:53.519.fits'                                                                               │                             │
│ │          FILE013 = 'ADP.2022-11-04T12:14:18.430.fits'                                                                               │                             │
│ │          FILE014 = 'ADP.2022-12-05T12:04:05.087.fits'                                                                               │                             │
│ │      N = 824989                                                                                                                     │                             │
│ │   self = Spectrum(fname='/home/sousasag/Data/ESPRESSO/Melissa_espresso_parameters/combined_spectra/aresMo'+30, resolution='115000') │                             │
│ │     w0 = 3771.560370205989                                                                                                          │                             │
│ ╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯                             │
│                                                                                                                                                                     │
│ /home/sousasag/Programas/anaconda3/envs/saitama/lib/python3.11/site-packages/scipy/interpolate/_polyint.py:81 in __call__                                           │
│                                                                                                                                                                     │
│    78 │   │                                                                                                                                                         │
│    79 │   │   """                                                                                                                                                   │
│    80 │   │   x, x_shape = self._prepare_x(x)                                                                                                                       │
│ ❱  81 │   │   y = self._evaluate(x)                                                                                                                                 │
│    82 │   │   return self._finish_y(y, x_shape)                                                                                                                     │
│    83 │                                                                                                                                                             │
│    84 │   def _evaluate(self, x):                                                                                                                                   │
│                                                                                                                                                                     │
│ ╭───────────────────────────────────── locals ──────────────────────────────────────╮                                                                               │
│ │    self = <scipy.interpolate._interpolate.interp1d object at 0x783107741360>      │                                                                               │
│ │       x = array([3771.56037021, 3771.57037021, 3771.58037021, ..., 7896.4803703 , │                                                                               │
│ │           │      7896.4903703 , 7896.5003703 ])                                   │                                                                               │
│ │ x_shape = (412495,)                                                               │                                                                               │
│ ╰───────────────────────────────────────────────────────────────────────────────────╯                                                                               │
│                                                                                                                                                                     │
│ /home/sousasag/Programas/anaconda3/envs/saitama/lib/python3.11/site-packages/scipy/interpolate/_interpolate.py:533 in _evaluate                                     │
│                                                                                                                                                                     │
│    530 │   │   x_new = asarray(x_new)                                                                                                                               │
│    531 │   │   y_new = self._call(self, x_new)                                                                                                                      │
│    532 │   │   if not self._extrapolate:                                                                                                                            │
│ ❱  533 │   │   │   below_bounds, above_bounds = self._check_bounds(x_new)                                                                                           │
│    534 │   │   │   if len(y_new) > 0:                                                                                                                               │
│    535 │   │   │   │   # Note fill_value must be broadcast up to the proper size                                                                                    │
│    536 │   │   │   │   # and flattened to work here                                                                                                                 │
│                                                                                                                                                                     │
│ ╭──────────────────────────────────── locals ─────────────────────────────────────╮                                                                                 │
│ │  self = <scipy.interpolate._interpolate.interp1d object at 0x783107741360>      │                                                                                 │
│ │ x_new = array([3771.56037021, 3771.57037021, 3771.58037021, ..., 7896.4803703 , │                                                                                 │
│ │         │      7896.4903703 , 7896.5003703 ])                                   │                                                                                 │
│ │ y_new = array([[nan],                                                           │                                                                                 │
│ │         │      [nan],                                                           │                                                                                 │
│ │         │      [nan],                                                           │                                                                                 │
│ │         │      ...,                                                             │                                                                                 │
│ │         │      [nan],                                                           │                                                                                 │
│ │         │      [nan],                                                           │                                                                                 │
│ │         │      [nan]])                                                          │                                                                                 │
│ ╰─────────────────────────────────────────────────────────────────────────────────╯                                                                                 │
│                                                                                                                                                                     │
│ /home/sousasag/Programas/anaconda3/envs/saitama/lib/python3.11/site-packages/scipy/interpolate/_interpolate.py:566 in _check_bounds                                 │
│                                                                                                                                                                     │
│    563 │   │   │   │   │   │   │    f"the interpolation range's minimum value ({self.x[0]}).")                                                                      │
│    564 │   │   if self.bounds_error and above_bounds.any():                                                                                                         │
│    565 │   │   │   above_bounds_value = x_new[np.argmax(above_bounds)]                                                                                              │
│ ❱  566 │   │   │   raise ValueError(f"A value ({above_bounds_value}) in x_new is above "                                                                            │
│    567 │   │   │   │   │   │   │    f"the interpolation range's maximum value ({self.x[-1]}).")                                                                     │
│    568 │   │                                                                                                                                                        │
│    569 │   │   # !! Should we emit a warning if some values are out of bounds?                                                                                      │
│                                                                                                                                                                     │
│ ╭─────────────────────────────────────────── locals ───────────────────────────────────────────╮                                                                    │
│ │       above_bounds = array([False, False, False, ..., False, False,  True])                  │                                                                    │
│ │ above_bounds_value = 7896.5003702960275                                                      │                                                                    │
│ │       below_bounds = array([False, False, False, ..., False, False, False])                  │                                                                    │
│ │               self = <scipy.interpolate._interpolate.interp1d object at 0x783107741360>      │                                                                    │
│ │              x_new = array([3771.56037021, 3771.57037021, 3771.58037021, ..., 7896.4803703 , │                                                                    │
│ │                      │      7896.4903703 , 7896.5003703 ])                                   │                                                                    │
│ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯                                                                    │
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯
ValueError: A value (7896.5003702960275) in x_new is above the interpolation range's maximum value (7896.500370205989).
DanielAndreasen commented 2 months ago

Hi Sérgio,

1) So you fixed it with adding fill_value='extrapolate'? Do you think this is a suitable idea in general? 2) How was the installation process?

sousasag commented 2 months ago

Hi,

  1. I think the extrapolation is only done in very few cases, I tried several spectra and I only had 1, maybe 2, that raised this error. On top of that, when we need extrapolation is just for the borders of the spectrum, which I am not even sure they are used for the pseudo ews, but even in case they are, it should be just for very few Angstroms at most. It is not the perfect solution, but would say more than reasonable.
  2. The installation was very smooth and easy. No problem at all so far :)
DanielAndreasen commented 2 months ago
  1. I'm fine with implementing this. Let's see what Alex thinks.
  2. Glad to hear. We wanted someone else to try the new installation as well :)
AlexandrosAntoniadis commented 1 month ago

Hi,

Yes, let's implement the automatic solution proposed by Sergio to this potential issue regarding the spectral borders, for which in the previous version of the code the user had to adjust manually a variable.

Thank you both for your feedback.

DanielAndreasen commented 1 month ago

The pull request is ready.