issue when setting multiple bands in bandflux / bandmag

[MickaelRigault]

Hello guys,

given the documentation of sncosmo, I should be able to input a list of bands when I want to get a lightcurve (bandflux or bandmag) for a given model.
I fails here: sncosmo/models.py:151, in _bandflux
because:

File ~/miniforge3/lib/python3.9/site-packages/sncosmo/models.py:151, in _bandflux(model, band, time_or_phase, zp, zpsys)
    148 bandflux = np.zeros(time_or_phase.shape, dtype=float)
    150 # Loop over unique bands.
--> 151 for b in set(band):
    152     mask = band == b
    153     b = get_bandpass(b)

TypeError: unhashable type: 'numpy.ndarray'

Here is a minimal example:

import sncosmo # I have version 2.8.0
import numpy as np # I have version 1.22.4
model = sncosmo.Model("salt2")
times = np.linspace(-20, 50, 100)

This works:

model.bandflux("ztfg", times)

This does not for broadcast reasons, but it makes sense

model.bandflux(["ztfg", "ztfr"], times)

This should work but does not:

model.bandflux(["ztfg", "ztfr"], times[:,None])

[rbiswas4]

Hi Mickael,

We can take a look at the documentation and see if it can be made a little clearer, but the issue here is that the length of times must match the length of the filters list. Without this, it is unclear which filter to use at which of time.

Hence:

model.bandflux(["ztfg", "ztfr"], times[:2])
>>> array([0.        , 0.07692057])

works.

If I wanted to do what I think you want (fluxes at times in both ztfg and ztfr), I would probably do

 b1 = np.full_like(times, fill_value='ztfg', dtype='<U4')
b2 = np.full_like(times, fill_value='ztfr', dtype='<U4')
b = np.hstack((b1, b2))
t = np.hstack((times, times))
fluxes = model.bandflux(b, t)

[MickaelRigault]

Hi @rbiswas4,

Thanks for the quick reply. This is indeed what I want, and what you provide almost works.
It does when adding:
fluxes.reshape(2, len(times))

It is however not very clean and satisfying.

I would like to do something like this:

from scipy import stats
xx = np.linspace(-10,10, 100)
means = [0,1,2]
pdfs = stats.norm.pdf(xx[:,None], loc=means, scale=1).T
# here pdfs is a 3*1000 array, 1 pdf for the full xx for each "loc"

Does this make sense ?
I find it more numpy/scipy python-spirit

If this could be implemented that would be nice, if not, the solution you made (+ the reshape) works.

[rbiswas4]

Hi @MickaelRigault, (reply updated with url for sncosmo-eps mentioned before)

Thanks for your suggestions. While I agree that the kind of API you show is closer to what we expect in numpy, the API here is meant to support getting band fluxes and magnitudes of supernovae observed in real surveys. And surveys do not observe in all bands at the same time. So, a numpy array of shape (times, num_bands) is not appropriate for this kind of data (unless we then apply masks). The expectation here is that this function is most often called during simulations and light curve fitting and must reflect the real observations there. Additionally, since this API has been stable for a while and is used by other codes, I am not sure it is advisable to change this (say, for example, using a numpy array and masks) at this point.

That said, I have needed to do exactly what you were trying many times before. So, I think having functionality that achieves this is beneficial. There are ways to achieve what you are suggesting, without changing the API behavior in use. One possibility is to overload the function, so that when the length of times exceeds the number of elements in the list of bandpasses, this behavior is used. I feel that would be complicated. A second possibility, that I would personally prefer is having a separate function like multibandflux, and multibandmag provide this. What do you feel about that? There is now a process of discussion for such feature requests that have to be followed.

BTW, I also saw your message about a different way to do what I suggested and actually prefer your resize method to what I had done.

[benjaminrose]

@MickaelRigault, is this now a request to clarify the best practice in the documentation?

[MickaelRigault]

that could be nice indeed if there is no aim to change the API.
Examples would be enough I beleive.