Skip to content

Commit

Permalink
fix computation of d6 emission
Browse files Browse the repository at this point in the history
  • Loading branch information
xgarrido committed Mar 12, 2020
1 parent 240b570 commit e5e5274
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 19 deletions.
48 changes: 30 additions & 18 deletions pysm/models/dust.py
Original file line number Diff line number Diff line change
Expand Up @@ -185,51 +185,63 @@ def __init__(
unit_mbb_temperature=unit_mbb_temperature,
map_dist=map_dist
)
self.correlation_length = correlation_length
self.correlation_length = correlation_length * u.dimensionless_unscaled

def get_emission(self, freqs):
@u.quantity_input
def get_emission(self, freqs: u.GHz, weights=None) -> u.uK_RJ:
""" Function to calculate the emission of a decorrelated modified black
body model.
"""
freqs = utils.check_freq_input(freqs)
weights = utils.normalize_weights(freqs, weights)
# calculate the decorrelation
(rho_cov_I, rho_mean_I) = get_decorrelation_matrix(
rho_cov_I, rho_mean_I = get_decorrelation_matrix(
self.freq_ref_I, freqs, self.correlation_length
)
(rho_cov_P, rho_mean_P) = get_decorrelation_matrix(
rho_cov_P, rho_mean_P = get_decorrelation_matrix(
self.freq_ref_P, freqs, self.correlation_length
)
nfreqs = freqs.shape[-1]
extra_I = np.dot(rho_cov_I, np.random.randn(nfreqs))
extra_P = np.dot(rho_cov_P, np.random.randn(nfreqs))

decorr = np.zeros((nfreqs, 3))
decorr[:, 0, None] = rho_mean_I + extra_I[:, None]
decorr[:, 1, None] = rho_mean_P + extra_P[:, None]
decorr[:, 2, None] = rho_mean_P + extra_P[:, None]
# apply the decorrelation to the mbb_emission
return decorr[..., None] * super().get_emission(freqs)


@u.quantity_input(freqs=u.GHz, correlation_length=u.dimensionless_unscaled)
def frequency_decorr_model(freqs, correlation_length) -> u.dimensionless_unscaled:
output = np.zeros((3, len(self.I_ref)), dtype=self.I_ref.dtype)
# apply the decorrelation to the mbb_emission for each frequencies before integrating
for i, (freq, weight) in enumerate(zip(freqs, weights)):
temp = decorr[..., None][i] * super().get_emission(freq)
if len(freqs) > 1:
utils.trapz_step_inplace(freqs, weights, i, temp, output)
else:
output = temp
return output << u.uK_RJ


@u.quantity_input
def frequency_decorr_model(
freqs: u.GHz,
correlation_length: u.dimensionless_unscaled
):
""" Function to calculate the frequency decorrelation method of
Vansyngel+17.
"""
log_dep = np.log(freqs[:, None] / freqs[None, :])
return np.exp(-0.5 * (log_dep / correlation_length) ** 2)


@u.quantity_input(
freq_constrained=u.GHz,
freqs_constrained=u.GHz,
correlation_length=u.dimensionless_unscaled,
)
@u.quantity_input
def get_decorrelation_matrix(
freq_constrained, freqs_unconstrained, correlation_length
) -> u.dimensionless_unscaled:
freq_constrained: u.GHz,
freqs_unconstrained: u.GHz,
correlation_length: u.dimensionless_unscaled
):
""" Function to calculate the correlation matrix between observed
frequencies. This model is based on the proposed model for decorrelation
of Vanyngel+17. The proposed frequency covariance matrix in this paper
of Vansyngel+17. The proposed frequency covariance matrix in this paper
is implemented, and a constrained Gaussian realization for the unobserved
frequencies is calculated.
Expand Down Expand Up @@ -269,7 +281,7 @@ def get_decorrelation_matrix(
evals = np.diag(np.sqrt(np.maximum(rho_uu_w, np.zeros_like(rho_uu_w))))
rho_covar = np.dot(rho_uu_v, np.dot(evals, np.transpose(rho_uu_v)))
rho_mean = -np.dot(rho_uu, rho_inv_cu)
return (rho_covar, rho_mean)
return rho_covar, rho_mean


def invert_safe(matrix):
Expand Down
2 changes: 1 addition & 1 deletion pysm/utils/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ def check_freq_input(freqs):
""" Function to check that the input to `Model.get_emission` is a
np.ndarray.
This function will convet input integers or arrays to a single element
This function will convert input integers or arrays to a single element
numpy array.
Parameters
Expand Down

0 comments on commit e5e5274

Please # to comment.