Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add redshift effect #16

Merged
merged 12 commits into from
Jul 17, 2024
Merged
11 changes: 0 additions & 11 deletions src/tdastro/effects/effect_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,17 +9,6 @@ class EffectModel(ParameterizedNode):
def __init__(self, **kwargs):
super().__init__(**kwargs)

def required_parameters(self):
"""Returns a list of the parameters of a PhysicalModel
that this effect needs to access.

Returns
-------
parameters : `list` of `str`
A list of every required parameter the effect needs.
"""
return []

def apply(self, flux_density, wavelengths=None, physical_model=None, **kwargs):
"""Apply the effect to observations (flux_density values)

Expand Down
86 changes: 86 additions & 0 deletions src/tdastro/effects/redshift.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,86 @@
from tdastro.effects.effect_model import EffectModel


class Redshift(EffectModel):
"""A redshift effect model.

This contains a "pre-effect" method, which is used to calculate the emitted wavelengths/times
needed to give us the observed wavelengths and times given the redshift. Times are calculated
with respect to the t0 of the given model.

Notes
-----
Conversions used are as follows:
- rest_frame_wavelength = observation_frame_wavelength / (1 + redshift)
- rest_frame_time = (observation_frame_time - t0) / (1 + redshift) + t0
- observation_frame_flux = rest_frame_flux / (1 + redshift)
"""

def __init__(self, redshift=None, t0=None, **kwargs):
"""Create a Redshift effect model.

Parameters
----------
redshift : `float`
The redshift.
t0 : `float`
The reference epoch; e.g. the time of the maximum light of a supernova or the epoch of zero phase
for a periodic variable star.
**kwargs : `dict`, optional
Any additional keyword arguments.
"""
super().__init__(**kwargs)
self.add_parameter("redshift", redshift, required=True, **kwargs)
self.add_parameter("t0", t0, required=True, **kwargs)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I hadn't thought about doing things this way, but this might be preferable.

The way the code is currently structured these are not meant to be included as parameters for the effect model, but rather to pull them from the physical model (via required_parameters). So for applying an effect, we would do something like: physical_model.redshift.

But this approach could be cleaner. We could get rid of required_parameters altogether and just pass in the PhysicalModel the way you do in the test scripts.

I need to think through the pros and cons a bit.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe that redshift and t0 are physical model's parameters, not effect's parameters


def __str__(self) -> str:
"""Return a string representation of the Redshift effect model."""
return f"RedshiftEffect(redshift={self.redshift})"

def pre_effect(self, observer_frame_times, observer_frame_wavelengths, **kwargs):
"""Calculate the rest-frame times and wavelengths needed to give us the observer-frame times
and wavelengths (given the redshift).

Parameters
----------
observer_frame_times : numpy.ndarray
The times at which the observation is made.
observer_frame_wavelengths : numpy.ndarray
The wavelengths at which the observation is made.
**kwargs : `dict`, optional
Any additional keyword arguments.

Returns
-------
tuple of (numpy.ndarray, numpy.ndarray)
OliviaLynn marked this conversation as resolved.
Show resolved Hide resolved
The rest-frame times and wavelengths needed to generate the rest-frame flux densities,
which will later be redshifted back to observer-frame flux densities at the observer-frame
times and wavelengths.
"""
observed_times_rel_to_t0 = observer_frame_times - self.t0
rest_frame_times_rel_to_t0 = observed_times_rel_to_t0 / (1 + self.redshift)
rest_frame_times = rest_frame_times_rel_to_t0 + self.t0
rest_frame_wavelengths = observer_frame_wavelengths / (1 + self.redshift)
return (rest_frame_times, rest_frame_wavelengths)

def apply(self, flux_density, wavelengths, physical_model=None, **kwargs):
"""Apply the effect to observations (flux_density values).

Parameters
----------
flux_density : `numpy.ndarray`
A length T X N matrix of flux density values.
wavelengths : `numpy.ndarray`, optional
A length N array of wavelengths.
physical_model : `PhysicalModel`
A PhysicalModel from which the effect may query parameters such as redshift, position, or
distance.
**kwargs : `dict`, optional
Any additional keyword arguments.

Returns
-------
flux_density : `numpy.ndarray`
The redshifted results.
"""
return flux_density / (1 + self.redshift)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See the overall comment in the __init__() function about overall approach.

In the current set up, this should be physical_model.redshift since we want to use a specific redshift for each object.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For this and my comments throughout: Olivia and I chatted about it today. If we make effects behave like other nodes in the graph, they will have consistent sampling behavior. In that case we will want to have the effect take in the PhysicalModel's parameters.

The pros of this approach is that 1) it keeps everything in the same graph formulation, 2) keeps the same sampling counter throughout the entire graph, 3) doesn't require a specialized interface for effects (the required_parameters list), and 4) allows for stateful effects.

The cons are: 1) as @hombit) notes there is a chance this could be misspecified during model creation, 2) it requires creating a unique effect for each physical model.

Let's go with the current approach for now. If this causes too much confusion, we can change back to stateless effects.

11 changes: 5 additions & 6 deletions src/tdastro/sources/physical_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,12 +79,6 @@ def add_effect(self, effect, allow_dups=True, **kwargs):
if effect_type == type(prev):
raise ValueError("Added the effect type to a model {effect_type} more than once.")

required: list = effect.required_parameters()
for parameter in required:
# Raise an AttributeError if the parameter is missing or set to None.
if getattr(self, parameter) is None:
raise AttributeError(f"Parameter {parameter} unset for model {type(self).__name__}")

self.effects.append(effect)

def _evaluate(self, times, wavelengths, **kwargs):
Expand Down Expand Up @@ -129,6 +123,11 @@ def evaluate(self, times, wavelengths, resample_parameters=False, **kwargs):
if resample_parameters:
self.sample_parameters(kwargs)

# Pre-effects are adjustments done to times and/or wavelengths, before flux density computation.
for effect in self.effects:
if hasattr(effect, "pre_effect"):
times, wavelengths = effect.pre_effect(times, wavelengths, **kwargs)

# Compute the flux density for both the current object and add in anything
# behind it, such as a host galaxy.
flux_density = self._evaluate(times, wavelengths, **kwargs)
Expand Down
71 changes: 71 additions & 0 deletions tests/tdastro/effects/test_redshift.py
OliviaLynn marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
import numpy as np
from tdastro.effects.redshift import Redshift
from tdastro.sources.step_source import StepSource


def get_no_effect_and_redshifted_values(times, wavelengths, t0, t1, brightness, redshift) -> tuple:
"""Get the values for a source with no effects and a redshifted source."""
model_no_effects = StepSource(brightness=brightness, t0=t0, t1=t1)
model_redshift = StepSource(brightness=brightness, t0=t0, t1=t1, redshift=redshift)
model_redshift.add_effect(Redshift(redshift=model_redshift, t0=model_redshift))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See the overall comment in the __init__() function about overall approach.

If we stick with the current approach, we could just make this: model_redshift.add_effect(Redshift())

However I'm thinking your approach might be cleaner.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with Jeremy, we shouldn't pass redshifts and t0s twice to both classes, even implicitly. We might see bad bugs with these approach when parameters are random number generators...


values_no_effects = model_no_effects.evaluate(times, wavelengths)
values_redshift = model_redshift.evaluate(times, wavelengths)

# Check shape of output is as expected
assert values_no_effects.shape == (len(times), len(wavelengths))
assert values_redshift.shape == (len(times), len(wavelengths))

return values_no_effects, values_redshift


def test_redshift() -> None:
"""Test that we can create a Redshift object and it gives us values as expected."""
times = np.array([1, 2, 3, 5, 10])
wavelengths = np.array([100.0, 200.0, 300.0])
t0 = 1.0
t1 = 2.0
brightness = 15.0
redshift = 1.0

# Get the values for a redshifted step source, and a step source with no effects for comparison
(values_no_effects, values_redshift) = get_no_effect_and_redshifted_values(
times, wavelengths, t0, t1, brightness, redshift
)

# Check that the step source activates within the correct time range:
# For values_no_effects, the activated values are in the range [t0, t1]
for i, time in enumerate(times):
if t0 <= time and time <= t1:
assert np.all(values_no_effects[i] == brightness)
else:
assert np.all(values_no_effects[i] == 0.0)

# With redshift = 1.0, the activated values are *observed* in the range [(t0-t0)*(1+redshift)+t0,
# (t1-t0*(1+redshift+t0] (the first term reduces to t0). Also, the values are scaled by a factor
# of (1+redshift).
for i, time in enumerate(times):
if t0 <= time and time <= (t1 - t0) * (1 + redshift) + t0:
assert np.all(values_redshift[i] == brightness / (1 + redshift))
else:
assert np.all(values_redshift[i] == 0.0)


def test_other_redshift_values() -> None:
"""Test that we can create a Redshift object with various other redshift values."""
times = np.linspace(0, 100, 1000)
wavelengths = np.array([100.0, 200.0, 300.0])
t0 = 10.0
t1 = 30.0
brightness = 50.0

for redshift in [0.0, 0.5, 2.0, 3.0, 30.0]:
model_redshift = StepSource(brightness=brightness, t0=t0, t1=t1, redshift=redshift)
model_redshift.add_effect(Redshift(redshift=model_redshift, t0=model_redshift))
values_redshift = model_redshift.evaluate(times, wavelengths)

for i, time in enumerate(times):
if t0 <= time and time <= (t1 - t0) * (1 + redshift) + t0:
assert np.all(values_redshift[i] == brightness / (1 + redshift))
else:
assert np.all(values_redshift[i] == 0.0)
Loading