cosmodesi / pyrecon

package for BAO reconstruction
BSD 3-Clause "New" or "Revised" License
9 stars 8 forks source link

Small issue with the DistanceToRedshift utility #14

Open epaillas opened 4 months ago

epaillas commented 4 months ago

If you provide an integer to the DistanceToRedshift utility function, it will always spit out an integer redshift:

In [1]: from cosmoprimo.fiducial import DESI
In [2]: from pyrecon.utils import DistanceToRedshift

In [3]: cosmo = DESI()

In [4]: distance = cosmo.comoving_radial_distance

In [5]: d2r = DistanceToRedshift(distance)

In [6]: d2r(1000)
Out[6]: array(0)

In [7]: d2r(1000.0)
Out[7]: array(0.36653534)

which I don't think is desired. This is probably caused by the astype call function, which I guess you introduced to distinguish between arrays and scalars:

def __call__(self, distance):
    """Return (interpolated) redshift at distance ``distance`` (scalar or array)."""
    distance = np.asarray(distance)
    return self.interp(distance).astype(distance.dtype, copy=False)

Maybe adding a distance = float(distance) somewhere internally would take care of this?

adematti commented 4 months ago

Hi! This astype call was to keep the same input type (float32 vs float64). I agree we need to enforce float!

def __call__(self, distance):
    """Return (interpolated) redshift at distance ``distance`` (scalar or array)."""
    distance = np.asarray(distance)
    if not np.issubdtype(distance.dtype, np.floating): distance = distance.astype(np.float64)
    return self.interp(distance).astype(distance.dtype, copy=False)

Do you want to push this?