Starfish-develop / Starfish

Tools for Flexible Spectroscopic Inference
https://starfish.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
69 stars 22 forks source link

Creating interpolator fails. #146

Closed cosmocaos closed 2 years ago

cosmocaos commented 2 years ago

After making a HDF5 grid, and loading successfully the flux for combinations of parameters in the grid. I tried to interpolate a spectrum:

myHDF5 = HDF5Interface(filename=inputfile) flux = myHDF5.load_flux(np.array([6200, 4.5, 0.0])) myInterpolator = Interpolator(myHDF5)

Which gives the error:

keeping grid as is /opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Starfish/grid_tools/interpolators.py:22: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray. parameter_list = np.asarray(parameter_list)

ValueError Traceback (most recent call last) ~/GDrive/RESEARCH/GCK/code/galex111.py in ----> 1 myInterpolator = Interpolator(myHDF5)

/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Starfish/grid_tools/interpolators.py in init(self, interface, wl_range, cache_max, cache_dump) 88 self._determine_chunk_log() 89 ---> 90 self._setup_index_interpolators() 91 self.cache = OrderedDict([]) 92 self.cache_max = cache_max

/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Starfish/grid_tools/interpolators.py in _setup_index_interpolators(self) 138 # Given a parameter value, produce fractional index between two points 139 # Store the interpolators as a list --> 140 self.index_interpolator = IndexInterpolator(self.interface.points) 141 142 lenF = self.interface.ind[1] - self.interface.ind[0]

/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Starfish/grid_tools/interpolators.py in init(self, parameter_list) 22 parameter_list = np.asarray(parameter_list) 23 self.npars = parameter_list.shape[-1] ---> 24 self.parameter_list = np.unique(parameter_list) 25 self.index_interpolator = interp1d( 26 self.parameter_list, np.arange(len(self.parameter_list)), kind="linear"

<__array_function__ internals> in unique(*args, **kwargs) ~/Library/Python/3.7/lib/python/site-packages/numpy/lib/arraysetops.py in unique(ar, return_index, return_inverse, return_counts, axis) 270 ar = np.asanyarray(ar) 271 if axis is None: --> 272 ret = _unique1d(ar, return_index, return_inverse, return_counts) 273 return _unpack_tuple(ret) 274 ~/Library/Python/3.7/lib/python/site-packages/numpy/lib/arraysetops.py in _unique1d(ar, return_index, return_inverse, return_counts) 331 aux = ar[perm] 332 else: --> 333 ar.sort() 334 aux = ar 335 mask = np.empty(aux.shape, dtype=np.bool_) ValueError: operands could not be broadcast together with shapes (12,) (73,)
mileslucas commented 2 years ago

Hi @cosmocaos thanks for reaching out!

Can you tell me a little more about the spectral model grid you're using? At first glance from this error, I want to make sure that all of the spectra are evaluated on the same wavelength grid- if they aren't the interpolator and spectral emulator will fail, but there's a workaround.

cosmocaos commented 2 years ago

Hi @mileslucas, thanks for replying,

I created the spectral grid with the code below:

from Starfish.grid_tools import PHOENIXGridInterfaceNoAlpha as PHOENIX from Starfish.grid_tools import HDF5Creator, Instrument from Starfish.grid_tools import Interpolator

infile = '/Volumes/cosmocaos/Research/libraries/raw/PHOENIX/' mygrid = PHOENIX(path=infile, air=False, wl_range=[880,23000])
instrumento = Instrument(name='miinstrumento', FWHM=1000, wl_range=(880,23000), oversampling=4.0) outfile = computer + '/GDrive/RESEARCH/GCK/data/spectral_libraries/PHOENIX_lowres_Starfish3.hdf5' creator = HDF5Creator(mygrid,instrument=instrumento,filename=outfile) creator.process_grid()

The PHOENIX spectral library I'm using do not have alpha enhancement. With this format:

Captura de Pantalla 2022-06-05 a la(s) 7 26 05 p  m
mileslucas commented 2 years ago

Okay you're using the PHOENIX grid- that means there's some bug for me to track down! It might take me a day or two but I'll let you know when I find it and start working on a fix. Thanks for your patience in advance

cosmocaos commented 2 years ago

Thanks!

mileslucas commented 2 years ago

Found the bug! I had test code that was different than how the interpolator is actually used, so I thought everything was working when it actually wasn't. I've tested by hand that you can interpolate spectra after creating an interpolator from an interface-

In [1]: from Starfish.grid_tools import HDF5Interface, Interpolator

In [2]: grid = HDF5Interface("F_SPEX_grid.hdf5")

In [3]: itp = Interpolator(grid)
keeping grid as is

In [4]: itp([6000, 4.1, 0.05])
Out[4]: 
array([2069306.05347025, 2881213.91452267, 3505401.2814753 , ...,
        628004.53663218, 1251478.9462981 , 2062678.80636186])
cosmocaos commented 2 years ago

Thanks @mileslucas, its working now!