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.

Sign up for GitHub

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

Jump to bottom

Use native mapmaker #352

Merged
merged 7 commits into from
Jun 2, 2020
Merged

Use native mapmaker #352

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
a5c7fde
initial support for the native mapmaker in toast_ground_sim.py
keskitalo May 29, 2020
a1bc81e
smooth some rough edges
keskitalo May 29, 2020
cf0f37e
run source formatter
keskitalo May 29, 2020
7ede2b8
Handle an unusual case where the user wants to subtract zeroth order …
keskitalo Jun 1, 2020
d04410e
use the native mapmaker in the satellite pipeline
keskitalo Jun 1, 2020
3ee0c07
Write out the mapmaker timing profile, if the user enabled function t…
keskitalo Jun 1, 2020
b7df3a6
Fix merging errors
keskitalo Jun 1, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
103 changes: 73 additions & 30 deletions pipelines/toast_ground_sim.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ def parse_arguments(comm):
pipeline_tools.add_noise_args(parser)
pipeline_tools.add_gainscrambler_args(parser)
pipeline_tools.add_madam_args(parser)
pipeline_tools.add_mapmaker_args(parser)
pipeline_tools.add_sky_map_args(parser)
pipeline_tools.add_pysm_args(parser)
pipeline_tools.add_sss_args(parser)
Expand All @@ -81,6 +82,22 @@ def parse_arguments(comm):
"--outdir", required=False, default="out", help="Output directory"
)

parser.add_argument(
"--madam",
required=False,
action="store_true",
help="Use libmadam for map-making",
dest="use_madam",
)
parser.add_argument(
"--no-madam",
required=False,
action="store_false",
help="Do not use libmadam for map-making [default]",
dest="use_madam",
)
parser.set_defaults(use_madam=False)

parser.add_argument(
"--focalplane",
required=False,
Expand Down Expand Up @@ -387,9 +404,9 @@ def main():

args, comm = parse_arguments(comm)

# Initialize madam parameters

madampars = pipeline_tools.setup_madam(args)
if args.use_madam:
# Initialize madam parameters
madampars = pipeline_tools.setup_madam(args)

# Load and broadcast the schedule file

Expand Down Expand Up @@ -505,30 +522,7 @@ def main():

# Bin and destripe maps

pipeline_tools.apply_madam(
args,
comm,
data,
madampars,
outpath,
detweights,
totalname_freq,
freq=freq,
time_comms=time_comms,
telescope_data=telescope_data,
first_call=(mc == firstmc),
)

if args.apply_polyfilter or args.apply_groundfilter:

# Filter signal

pipeline_tools.apply_polyfilter(args, comm, data, totalname_freq)

pipeline_tools.apply_groundfilter(args, comm, data, totalname_freq)

# Bin filtered maps

if args.use_madam:
pipeline_tools.apply_madam(
args,
comm,
Expand All @@ -540,10 +534,59 @@ def main():
freq=freq,
time_comms=time_comms,
telescope_data=telescope_data,
first_call=False,
extra_prefix="filtered",
bin_only=True,
first_call=(mc == firstmc),
)
else:
pipeline_tools.apply_mapmaker(
args,
comm,
data,
outpath,
totalname_freq,
time_comms=time_comms,
telescope_data=telescope_data,
first_call=(mc == firstmc),
)

if args.apply_polyfilter or args.apply_groundfilter:

# Filter signal

pipeline_tools.apply_polyfilter(args, comm, data, totalname_freq)

pipeline_tools.apply_groundfilter(args, comm, data, totalname_freq)

# Bin filtered maps

if args.use_madam:
pipeline_tools.apply_madam(
args,
comm,
data,
madampars,
outpath,
detweights,
totalname_freq,
freq=freq,
time_comms=time_comms,
telescope_data=telescope_data,
first_call=False,
extra_prefix="filtered",
bin_only=True,
)
else:
pipeline_tools.apply_mapmaker(
args,
comm,
data,
outpath,
totalname_freq,
time_comms=time_comms,
telescope_data=telescope_data,
first_call=False,
extra_prefix="filtered",
bin_only=True,
)

gt.stop_all()
if mpiworld is not None:
Expand Down
160 changes: 62 additions & 98 deletions pipelines/toast_satellite_sim.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,7 @@ def parse_arguments(comm, procs):
)

pipeline_tools.add_madam_args(parser)
pipeline_tools.add_mapmaker_args(parser)
pipeline_tools.add_binner_args(parser)

parser.add_argument(
Expand Down Expand Up @@ -377,132 +378,95 @@ def main():
if comm.world_rank == 0:
tmr.report_clear("Dumping data distribution")

# Mapmaking.

if not args.use_madam:
# in debug mode, print out data distribution information
if args.debug:
handle = None
if comm.world_rank == 0:
log.info("Not using Madam, will only make a binned map")

npp, zmap = pipeline_tools.init_binner(args, comm, data, detweights)
handle = open(os.path.join(args.outdir, "distdata.txt"), "w")
data.info(handle)
if comm.world_rank == 0:
handle.close()
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear("Initialize binned map-making")

# Loop over Monte Carlos

firstmc = args.MC_start
nmc = args.MC_count

for mc in range(firstmc, firstmc + nmc):
mctmr = Timer()
mctmr.start()

outpath = os.path.join(args.outdir, "mc_{:03d}".format(mc))

pipeline_tools.simulate_noise(
args, comm, data, mc, "tot_signal", overwrite=True
)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Simulate noise {:04d}".format(mc))

# add sky signal
pipeline_tools.add_signal(args, comm, data, "tot_signal", signalname)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Add sky signal {:04d}".format(mc))

if gain is not None:
op_apply_gain = OpApplyGain(gain, name="tot_signal")
op_apply_gain.exec(data)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Apply gains {:04d}".format(mc))

if mc == firstmc:
# For the first realization, optionally export the
# timestream data. If we had observation intervals defined,
# we could pass "use_interval=True" to the export operators,
# which would ensure breaks in the exported data at
# acceptable places.
pipeline_tools.output_tidas(args, comm, data, "tot_signal")
pipeline_tools.output_spt3g(args, comm, data, "tot_signal")
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Write TOD snapshot {:04d}".format(mc))

pipeline_tools.apply_binner(
args, comm, data, npp, zmap, detweights, outpath, "tot_signal"
)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Apply binner {:04d}".format(mc))
tmr.report_clear("Dumping data distribution")

if comm.world_rank == 0:
mctmr.report_clear(" Map-making {:04d}".format(mc))
else:
# Mapmaking.

if args.use_madam:
# Initialize madam parameters

madampars = pipeline_tools.setup_madam(args)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear("Initialize madam map-making")

# Loop over Monte Carlos
# Loop over Monte Carlos

firstmc = args.MC_start
nmc = args.MC_count
firstmc = args.MC_start
nmc = args.MC_count

for mc in range(firstmc, firstmc + nmc):
mctmr = Timer()
mctmr.start()
for mc in range(firstmc, firstmc + nmc):
mctmr = Timer()
mctmr.start()

# create output directory for this realization
outpath = os.path.join(args.outdir, "mc_{:03d}".format(mc))
# create output directory for this realization
outpath = os.path.join(args.outdir, "mc_{:03d}".format(mc))

pipeline_tools.simulate_noise(
args, comm, data, mc, "tot_signal", overwrite=True
)
pipeline_tools.simulate_noise(
args, comm, data, mc, "tot_signal", overwrite=True
)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Simulate noise {:04d}".format(mc))

# add sky signal
pipeline_tools.add_signal(args, comm, data, "tot_signal", signalname)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Add sky signal {:04d}".format(mc))

if gain is not None:
op_apply_gain = OpApplyGain(gain, name="tot_signal")
op_apply_gain.exec(data)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Simulate noise {:04d}".format(mc))

# add sky signal
pipeline_tools.add_signal(args, comm, data, "tot_signal", signalname)
tmr.report_clear(" Apply gains {:04d}".format(mc))

if mc == firstmc:
# For the first realization, optionally export the
# timestream data. If we had observation intervals defined,
# we could pass "use_interval=True" to the export operators,
# which would ensure breaks in the exported data at
# acceptable places.
pipeline_tools.output_tidas(args, comm, data, "tot_signal")
pipeline_tools.output_spt3g(args, comm, data, "tot_signal")
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Add sky signal {:04d}".format(mc))

if gain is not None:
op_apply_gain = OpApplyGain(gain, name="tot_signal")
op_apply_gain.exec(data)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Apply gains {:04d}".format(mc))
tmr.report_clear(" Write TOD snapshot {:04d}".format(mc))

if args.use_madam:
pipeline_tools.apply_madam(
args, comm, data, madampars, outpath, detweights, "tot_signal"
)
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Apply madam {:04d}".format(mc))
else:
pipeline_tools.apply_mapmaker(
args, comm, data, outpath, "tot_signal"
)

if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
mctmr.report_clear(" Map-making {:04d}".format(mc))
if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
tmr.report_clear(" Map-making {:04d}".format(mc))

if comm.comm_world is not None:
comm.comm_world.barrier()
if comm.world_rank == 0:
mctmr.report_clear(" Monte Carlo loop {:04d}".format(mc))

gt.stop_all()
if comm.comm_world is not None:
Expand Down
14 changes: 7 additions & 7 deletions src/libtoast/src/toast_atm.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@
// All rights reserved. Use of this source code is governed by
// a BSD-style license that can be found in the LICENSE file.

//#if !defined(NO_ATM_CHECKS)
//# define NO_ATM_CHECKS
//#endif // if !defined(NO_ATM_CHECKS)
// #if !defined(NO_ATM_CHECKS)
// # define NO_ATM_CHECKS
// #endif // if !defined(NO_ATM_CHECKS)

#include <toast/sys_utils.hpp>
#include <toast/sys_environment.hpp>
Expand Down Expand Up @@ -1660,10 +1660,10 @@ cholmod_sparse * toast::atm_sim::build_sparse_covariance(long ind_start,
if (fabs(colcoord[2] - rowcoord[2]) > rcorr) continue;

double val = cov_eval(colcoord, rowcoord);
if (icol == irow) {
// Regularize the matrix by promoting the diagonal
val *= 1.01;
}
if (icol == irow) {
// Regularize the matrix by promoting the diagonal
val *= 1.01;
}

// If the covariance exceeds the threshold, add it to the
// sparse matrix
Expand Down
Loading