Skip to content

Commit

Permalink
Fix bug with unit conversion in bandpasses which were not fully fixed…
Browse files Browse the repository at this point in the history
… in 0.8.1
  • Loading branch information
MetinSa committed Jan 11, 2023
1 parent 861e51f commit 788beeb
Show file tree
Hide file tree
Showing 9 changed files with 59 additions and 38 deletions.
Binary file modified docs/img/bandpass_integrated.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified docs/img/center_freq.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
[tool.poetry]
name = "zodipy"
homepage = "https://github.com/Cosmoglobe/zodipy"
version = "0.8.2"
version = "0.8.3"
description = "Software for simulating zodiacal emission"
authors = ["Metin San <[email protected]>"]
readme = "README.md"
Expand Down
2 changes: 0 additions & 2 deletions tests/test_get_emission.py
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,6 @@ def test_get_emission_pix(
observer = data.draw(obs(model, time))
frequency = data.draw(freq(model))
pix = data.draw(pixels(nside))

emission = model.get_emission_pix(
frequency,
pixels=pix,
Expand Down Expand Up @@ -344,7 +343,6 @@ def test_bandpass_integration(
data: DataObject,
) -> None:
"""Property test for bandpass integrations."""

theta, phi = angles
model.extrapolate = True
observer = data.draw(obs(model, time))
Expand Down
28 changes: 14 additions & 14 deletions tests/test_validators.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
from hypothesis.strategies import DataObject, data

from zodipy._types import FrequencyOrWavelength
from zodipy._validators import validate_and_normalize_weights, validate_frequencies
from zodipy._validators import get_validated_and_normalized_weights, get_validated_freq
from zodipy.zodipy import Zodipy

from ._strategies import model, random_freq, weights
Expand All @@ -18,25 +18,25 @@
@given(model(extrapolate=False))
def test_validate_frequencies(model: Zodipy) -> None:
with pytest.raises(TypeError):
validate_frequencies(
get_validated_freq(
freq=BANDPASS_FREQUENCIES.value,
model=model.ipd_model,
extrapolate=model.extrapolate,
)
with pytest.raises(TypeError):
validate_frequencies(
get_validated_freq(
freq=25,
model=model.ipd_model,
extrapolate=model.extrapolate,
)
with pytest.raises(u.UnitsError):
validate_frequencies(
get_validated_freq(
freq=BANDPASS_FREQUENCIES.value * u.g,
model=model.ipd_model,
extrapolate=model.extrapolate,
)
with pytest.raises(u.UnitsError):
validate_frequencies(
get_validated_freq(
freq=25 * u.g,
model=model.ipd_model,
extrapolate=model.extrapolate,
Expand All @@ -46,57 +46,57 @@ def test_validate_frequencies(model: Zodipy) -> None:
@given(random_freq(bandpass=True), data())
def test_validate_weights(freq: FrequencyOrWavelength, data: DataObject) -> None:
bp_weights = data.draw(weights(freq))
bp_weights = validate_and_normalize_weights(
bp_weights = get_validated_and_normalized_weights(
weights=bp_weights,
freq=freq,
)
assert np.trapz(bp_weights, freq.value) == pytest.approx(1.0)

with pytest.raises(ValueError):
validate_and_normalize_weights(
get_validated_and_normalized_weights(
weights=None,
freq=BANDPASS_FREQUENCIES,
)


def test_validate_weights_numbers() -> None:
validate_and_normalize_weights(
get_validated_and_normalized_weights(
weights=None,
freq=BANDPASS_FREQUENCIES[0],
)
with pytest.raises(ValueError):
validate_and_normalize_weights(
get_validated_and_normalized_weights(
weights=np.array([1, 2, 3]),
freq=BANDPASS_FREQUENCIES,
)
with pytest.raises(ValueError):
validate_and_normalize_weights(
get_validated_and_normalized_weights(
weights=BANDPASS_WEIGHTS[:10],
freq=BANDPASS_FREQUENCIES[0],
)
with pytest.raises(ValueError):
validate_and_normalize_weights(
get_validated_and_normalized_weights(
weights=None,
freq=BANDPASS_FREQUENCIES,
)


def test_validate_weights_strictly_increasing() -> None:
with pytest.raises(ValueError):
validate_and_normalize_weights(
get_validated_and_normalized_weights(
weights=BANDPASS_WEIGHTS,
freq=np.flip(BANDPASS_FREQUENCIES),
)


def test_validate_weights_shape() -> None:
weights = validate_and_normalize_weights(
weights = get_validated_and_normalized_weights(
weights=BANDPASS_WEIGHTS,
freq=BANDPASS_FREQUENCIES,
)
assert weights.shape == BANDPASS_WEIGHTS.shape

weights = validate_and_normalize_weights(
weights = get_validated_and_normalized_weights(
weights=None,
freq=BANDPASS_FREQUENCIES[0],
)
Expand Down
9 changes: 4 additions & 5 deletions zodipy/_bandpass.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
from zodipy._ipd_model import InterplanetaryDustModel
from zodipy._source_funcs import get_blackbody_emission
from zodipy._types import FrequencyOrWavelength
from zodipy._validators import validate_and_normalize_weights, validate_frequencies
from zodipy._validators import get_validated_and_normalized_weights, get_validated_freq


@dataclass
Expand Down Expand Up @@ -46,8 +46,8 @@ def validate_and_get_bandpass(
extrapolate: bool,
) -> Bandpass:
"""Validate user inputted bandpass and return a Bandpass object."""
validate_frequencies(freq, model, extrapolate)
normalized_weights = validate_and_normalize_weights(weights, freq)
freq = get_validated_freq(freq, model, extrapolate)
normalized_weights = get_validated_and_normalized_weights(weights, freq)

return Bandpass(freq, normalized_weights)

Expand All @@ -60,10 +60,9 @@ def get_bandpass_interpolation_table(
) -> npt.NDArray[np.float64]:
"""Pre-compute the bandpass integrated blackbody emission for a grid of temperatures."""
# Prepare bandpass to be integrated in power units and in frequency convention.

if not bandpass.frequencies.unit.is_equivalent(u.Hz):
bandpass.switch_convention()
else:
bandpass.frequencies = bandpass.frequencies.to(u.Hz)

integrals = np.zeros(n_points)
temp_grid = np.linspace(min_temp, max_temp, n_points)
Expand Down
20 changes: 16 additions & 4 deletions zodipy/_interpolate_source.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,13 @@ def get_source_parameters_kelsall_comp(
if not bandpass.frequencies.unit.is_equivalent(model.spectrum.unit):
bandpass.switch_convention()

interpolator = partial(interp1d, x=model.spectrum.value, fill_value="extrapolate")
spectrum = (
model.spectrum.to_value(u.Hz)
if model.spectrum.unit.is_equivalent(u.Hz)
else model.spectrum.to_value(u.micron)
)

interpolator = partial(interp1d, x=spectrum, fill_value="extrapolate")

source_parameters: dict[ComponentLabel | str, dict[str, Any]] = {}
for comp_label in model.comps.keys():
Expand Down Expand Up @@ -91,10 +97,16 @@ def get_source_parameters_rmm(
if not bandpass.frequencies.unit.is_equivalent(model.spectrum.unit):
bandpass.switch_convention()

spectrum = (
model.spectrum.to_value(u.Hz)
if model.spectrum.unit.is_equivalent(u.Hz)
else model.spectrum.to_value(u.micron)
)

source_parameters: dict[ComponentLabel | str, dict[str, Any]] = {}
calibration = interp1d(
x=model.spectrum.value, y=model.calibration, fill_value="extrapolate"
)(bandpass.frequencies.value)
calibration = interp1d(x=spectrum, y=model.calibration, fill_value="extrapolate")(
bandpass.frequencies.value
)
calibration = u.Quantity(calibration, u.MJy / u.AU).to_value(u.Jy / u.cm)

if bandpass.frequencies.size > 1:
Expand Down
26 changes: 19 additions & 7 deletions zodipy/_validators.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,22 @@
from ._types import FrequencyOrWavelength, Pixels, SkyAngles


@u.quantity_input(equivalencies=u.spectral())
def validate_frequencies(
def get_validated_freq(
freq: FrequencyOrWavelength, model: InterplanetaryDustModel, extrapolate: bool
) -> None:
) -> FrequencyOrWavelength:
"""Validate user inputted frequency."""
if not isinstance(freq, u.Quantity):
raise TypeError("Frequency must be an astropy Quantity.")

if freq.unit.is_equivalent(u.Hz):
freq = freq.to(u.Hz)
elif freq.unit.is_equivalent(u.micron):
freq = freq.to(u.micron)
else:
raise u.UnitsError("Frequency must be in units compatible with Hz or micron.")

if extrapolate:
return
return freq

freq_in_spectrum_units = freq.to(model.spectrum.unit, equivalencies=u.spectral())
lower_freq_range = model.spectrum.min()
Expand All @@ -37,8 +46,10 @@ def validate_frequencies(
f" {upper_freq_range}] range."
)

return freq

def validate_and_normalize_weights(

def get_validated_and_normalized_weights(
weights: Union[Sequence[float], npt.NDArray[np.floating], None],
freq: FrequencyOrWavelength,
) -> npt.NDArray[np.float64]:
Expand All @@ -47,6 +58,7 @@ def validate_and_normalize_weights(
raise ValueError(
"Bandpass weights must be specified if more than one frequency is given."
)

if weights is not None:
if freq.size != len(weights):
raise ValueError("Number of frequencies and weights must be the same.")
Expand All @@ -64,7 +76,7 @@ def validate_and_normalize_weights(


@u.quantity_input(theta=[u.deg, u.rad], phi=[u.deg, u.rad])
def validate_ang(
def get_validated_ang(
theta: SkyAngles, phi: SkyAngles, lonlat: bool
) -> Tuple[SkyAngles, SkyAngles]:
"""Validate user inputted sky angles."""
Expand All @@ -79,7 +91,7 @@ def validate_ang(
return theta, phi


def validate_pixels(pixels: Pixels, nside: int) -> Pixels:
def get_validated_pix(pixels: Pixels, nside: int) -> Pixels:
"""Validate user inputted pixels."""
if (np.max(pixels) > hp.nside2npix(nside)) or (np.min(pixels) < 0):
raise ValueError("invalid pixel number given nside")
Expand Down
10 changes: 5 additions & 5 deletions zodipy/zodipy.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
from zodipy._sky_coords import get_obs_and_earth_positions
from zodipy._types import FrequencyOrWavelength, Pixels, SkyAngles
from zodipy._unit_vectors import get_unit_vectors_from_ang, get_unit_vectors_from_pixels
from zodipy._validators import validate_ang, validate_pixels
from zodipy._validators import get_validated_ang, get_validated_pix
from zodipy.model_registry import model_registry

PLATFORM = platform.system().lower()
Expand Down Expand Up @@ -171,7 +171,7 @@ def get_emission_ang(
emission: Simulated zodiacal emission in units of 'MJy/sr'.
"""
theta, phi = validate_ang(theta=theta, phi=phi, lonlat=lonlat)
theta, phi = get_validated_ang(theta=theta, phi=phi, lonlat=lonlat)

unique_angles, indicies = np.unique(
np.asarray([theta, phi]), return_inverse=True, axis=1
Expand Down Expand Up @@ -233,7 +233,7 @@ def get_emission_pix(
emission: Simulated zodiacal emission in units of 'MJy/sr'.
"""
pixels = validate_pixels(pixels=pixels, nside=nside)
pixels = get_validated_pix(pixels=pixels, nside=nside)

unique_pixels, indicies = np.unique(pixels, return_inverse=True)
unit_vectors = get_unit_vectors_from_pixels(
Expand Down Expand Up @@ -304,7 +304,7 @@ def get_binned_emission_ang(
emission: Simulated zodiacal emission in units of 'MJy/sr'.
"""
theta, phi = validate_ang(theta=theta, phi=phi, lonlat=lonlat)
theta, phi = get_validated_ang(theta=theta, phi=phi, lonlat=lonlat)

unique_angles, counts = np.unique(
np.asarray([theta, phi]), return_counts=True, axis=1
Expand Down Expand Up @@ -371,7 +371,7 @@ def get_binned_emission_pix(
emission: Simulated zodiacal emission in units of 'MJy/sr'.
"""
pixels = validate_pixels(pixels=pixels, nside=nside)
pixels = get_validated_pix(pixels=pixels, nside=nside)

unique_pixels, counts = np.unique(pixels, return_counts=True)
unit_vectors = get_unit_vectors_from_pixels(
Expand Down

0 comments on commit 788beeb

Please sign in to comment.