Skip to content

Commit

Permalink
Address review comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@tskisner
tskisner committed Aug 16, 2024
1 parent 5ff3aee commit 57a9e46
Show file tree
Hide file tree
Showing 2 changed files with 113 additions and 16 deletions.
103 changes: 102 additions & 1 deletion src/toast/hwp_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,19 @@


def hwpss_samples(n_samp, comm):
"""Helper function to distribute samples."""
"""Helper function to distribute samples.
This distributes slices of samples uniformly among the processes
of an MPI communicator.
Args:
n_samp (int): The number of samples
comm (MPI.Comm): The communicator.
Returns:
(slice): The local slice on this process.
"""
if comm is None:
slc = slice(0, n_samp, 1)
return slc
Expand All @@ -29,6 +41,23 @@ def hwpss_samples(n_samp, comm):


def hwpss_sincos_buffer(angles, flags, n_harmonics, comm=None):
"""Precompute sin / cos terms.
Compute the harmonic terms involving sin / cos of the HWP angle.
This is done once in parallel and used by each process for its
local detectors.
Args:
angles (array): The HWP angle at each time stamp
flags (array): The flags indicating bad angle samples.
n_harmonics (ing): The number of harmonics in the model.
comm (MPI.Comm): The optional communicator to parallelize
the calculation.
Returns:
(array): The full buffer of sin/cos factors on all processes.
"""
slc = hwpss_samples(len(angles), comm)
ang = np.copy(angles[slc])
good = flags[slc] == 0
Expand All @@ -45,6 +74,46 @@ def hwpss_sincos_buffer(angles, flags, n_harmonics, comm=None):


def hwpss_compute_coeff_covariance(times, flags, sincos, comm=None):
"""Build covariance of HWPSS model coefficients.
The HWPSS model for this function is the one used by the Maxipol
and EBEX experiments. See for example Joy Didier's thesis, equation
8.17: https://academiccommons.columbia.edu/doi/10.7916/D8MW2HCG
The model with N harmonics and HWP angle H can be written as:
h(t) = Sum(i=1...N) { [ C_i0 + C_i1 * t ] cos(i * H(t)) +
[ C_i2 + C_i3 * t ] sin(i * H(t)) ] }
Writing this in terms of the vector of coefficients and the matrix of
of sin / cos factors:
h(t) = M x C
h(t) =
[ cos(H(t_0)) t_0 cos(H(t_0)) sin(H(t_0)) t_0 sin(H(t_0)) cos(2H(t_0)) ...]
[ cos(H(t_1)) t_1 cos(H(t_1)) sin(H(t_1)) t_1 sin(H(t_1)) cos(2H(t_1)) ...]
[ ... ] X [ C_10 C_11 C_12 C_13 C_20 C_21 C_22 C_23 C_30 ... ]^T
The least squares solution for the coefficients is then
C = (M^T M)^-1 M^T h(t)
We then assume that h(t) is just the input data and that it is dominated by the
HWPSS. This function computes the covariance matrix and factors it for later
use.
Args:
times (array): The **relative** timestamps of the samples from the start
of the observation.
flags (array): The flags indicating bad angle samples
sincos (array): The pre-computed sin / cos terms.
comm (MPI.Comm): The optional communicator to parallelize
the calculation.
Returns:
(tuple): The LU factorization and pivots.
"""
slc = hwpss_samples(len(times), comm)
my_sincos = sincos[slc, :]
my_times = times[slc]
Expand Down Expand Up @@ -141,6 +210,25 @@ def hwpss_compute_coeff_covariance(times, flags, sincos, comm=None):


def hwpss_compute_coeff(detdata, flags, times, sincos, cov_lu, cov_piv):
"""Compute the HWPSS model coefficients.
See docstring for `hwpss_compute_coeff_covariance`. This function computes
the expression M^T h(t) and then uses the LU factorization of the covariance
to solve for the coefficients.
Args:
detdata (array): The detector data for one detector.
flags (array): The detector flags.
times (array): The **relative** timestamps of the samples from the start
of the observation.
sincos (array): The pre-computed sin / cos terms.
cov_lu (array): The covariance LU factorization.
cov_piv (array): The covariance pivots.
Returns:
(array): The coefficients.
"""
n_samp = len(times)
n_harmonics = sincos.shape[1] // 2
good = flags == 0
Expand All @@ -164,6 +252,19 @@ def hwpss_compute_coeff(detdata, flags, times, sincos, cov_lu, cov_piv):


def hwpss_build_model(times, flags, sincos, coeff):
"""Construct the HWPSS template from coefficients.
Args:
times (array): The **relative** timestamps of the samples from the start
of the observation.
flags (array): The flags indicating bad angle samples
sincos (array): The pre-computed sin / cos terms.
coeff (array): The model coefficents for this detector.
Returns:
(array): The template.
"""
n_samp = len(times)
n_harmonics = sincos.shape[1] // 2
good = flags == 0
Expand Down
26 changes: 11 additions & 15 deletions src/toast/ops/noise_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
from ..noise import Noise
from ..noise_sim import AnalyticNoise
from ..observation import default_values as defaults
from ..mpi import flatten
from ..timing import Timer, function_timer
from ..traits import Float, Int, Quantity, Unicode, trait_docs
from ..utils import Environment, Logger
Expand Down Expand Up @@ -518,9 +519,6 @@ def _fit_log_psd(self, freqs, data, guess=None):
bad = np.logical_not(good)
n_bad = np.count_nonzero(bad)
if n_bad > 0:



msg = "Some PSDs have negative values. Consider changing "
msg += "noise estimation parameters."
log.warning(msg)
Expand Down Expand Up @@ -635,7 +633,7 @@ class FlagNoiseFit(Operator):
sigma_rms = Float(
None,
allow_none=True,
help="In addition to flagging based on estimated model, also apply overall TOD cut"
help="In addition to flagging based on estimated model, also apply overall TOD cut",
)

sigma_NET = Float(5.0, help="Flag detectors with NET values outside this range")
Expand Down Expand Up @@ -717,16 +715,16 @@ def _exec(self, data, detectors=None, **kwargs):
else:
proc_vals = obs.comm_col.gather(local_net, root=0)
if obs.comm_col_rank == 0:
all_net = np.array([val for plist in proc_vals for val in plist])
all_net = np.array(flatten(proc_vals))
proc_vals = obs.comm_col.gather(local_fknee, root=0)
if obs.comm_col_rank == 0:
all_fknee = np.array([val for plist in proc_vals for val in plist])
all_fknee = np.array(flatten(proc_vals))
proc_vals = obs.comm_col.gather(local_rms, root=0)
if obs.comm_col_rank == 0:
all_rms = np.array([val for plist in proc_vals for val in plist])
all_rms = np.array(flatten(proc_vals))
proc_vals = obs.comm_col.gather(local_names, root=0)
if obs.comm_col_rank == 0:
all_names = [val for plist in proc_vals for val in plist]
all_names = flatten(proc_vals)

# Iteratively cut
all_flags = None
Expand Down Expand Up @@ -761,8 +759,8 @@ def _exec(self, data, detectors=None, **kwargs):
# Already cut
continue
if (
np.absolute(fknee - fknee_mean) >
fknee_std * self.sigma_fknee
np.absolute(fknee - fknee_mean)
> fknee_std * self.sigma_fknee
):
msg = f"obs {obs.name}, det {name} has f_knee "
msg += f"{fknee} that is > {self.sigma_fknee} "
Expand All @@ -777,10 +775,7 @@ def _exec(self, data, detectors=None, **kwargs):
if not all_good[idet]:
# Already cut
continue
if (
np.absolute(rms - rms_mean) >
rms_std * self.sigma_rms
):
if np.absolute(rms - rms_mean) > rms_std * self.sigma_rms:
msg = f"obs {obs.name}, det {name} has TOD RMS "
msg += f"{rms} that is > {self.sigma_rms} "
msg += f"x {rms_std} from {rms_mean}"
Expand All @@ -791,7 +786,8 @@ def _exec(self, data, detectors=None, **kwargs):
log.debug(msg)
flag_pass += 1
all_flags = {
x: self.outlier_flag_mask for i, x in enumerate(all_names)
x: self.outlier_flag_mask
for i, x in enumerate(all_names)
if not all_good[i]
}
msg = f"obs {obs.name}: flagged {len(all_flags)} / {len(all_names)}"
Expand Down

0 comments on commit 57a9e46

Please sign in to comment.