Weave.c

Go to the documentation of this file.

//
// Copyright (C) 2016, 2017 Karl Wette
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with with program; see the file COPYING. If not, write to the
// Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
// MA 02110-1301 USA
//
 
///
/// \file
/// \ingroup lalpulsar_bin_Weave
///
 
#include "config.h"
 
#include "Weave.h"
#include "SetupData.h"
#include "SearchIteration.h"
#include "ComputeResults.h"
#include "CacheResults.h"
#include "OutputResults.h"
#include "SearchTiming.h"
 
#ifdef LALPULSAR_CUDA_ENABLED
#include <cuda.h>
#include <cuda_runtime_api.h>
#endif
 
#include <lal/LogPrintf.h>
#include <lal/UserInput.h>
#include <lal/Random.h>
 
int main( int argc, char *argv[] )
{
 
  // Set help information
  lalUserVarHelpBrief = "search for gravitational-wave pulsars";
 
  ////////// Parse user input //////////
 
  // Optional arguments to XLALCreateFstatInput(), used to initialise some user input variables
  FstatOptionalArgs Fstat_opt_args = FstatOptionalArgsDefaults;
 
  // Initialise user input variables
  struct uvar_type {
    BOOLEAN validate_sft_files, interpolation, lattice_rand_offset, mean2F_hgrm, segment_info, simulate_search, time_search, cache_all_gc, strict_spindown_bounds;
    CHAR *setup_file, *sft_files, *output_file, *ckpt_output_file;
    LALStringVector *sft_timestamps_files, *sft_noise_sqrtSX, *injections, *Fstat_assume_sqrtSX, *lrs_oLGX;
    REAL8 sft_timebase, semi_max_mismatch, coh_max_mismatch, ckpt_output_period, ckpt_output_exit, lrs_Fstar0sc, nc_2Fth;
    REAL8Range alpha, delta, freq, f1dot, f2dot, f3dot, f4dot;
    REAL8Vector *random_injection;
    UINT4 sky_patch_count, sky_patch_index, freq_partitions, f1dot_partitions, Fstat_run_med_window, Fstat_Dterms, toplist_limit, rand_seed, cache_max_size;
    int lattice, Fstat_method, Fstat_SSB_precision, toplists, extra_statistics, recalc_statistics;
  } uvar_struct = {
    .Fstat_Dterms = Fstat_opt_args.Dterms,
    .Fstat_SSB_precision = Fstat_opt_args.SSBprec,
    .Fstat_method = FMETHOD_RESAMP_BEST,
    .Fstat_run_med_window = Fstat_opt_args.runningMedianWindow,
    .alpha = {0, LAL_TWOPI},
    .delta = {-LAL_PI_2, LAL_PI_2},
    .freq_partitions = 1,
    .f1dot_partitions = 1,
    .interpolation = 1,
    .lattice = TILING_LATTICE_ANSTAR,
    .toplist_limit = 1000,
    .toplists = WEAVE_STATISTIC_MEAN2F,
    .extra_statistics = WEAVE_STATISTIC_NONE,
    .recalc_statistics = WEAVE_STATISTIC_NONE,
    .nc_2Fth = 5.2,
  };
  struct uvar_type *const uvar = &uvar_struct;
 
  // Register user input variables:
  //
  // - General
  //
  XLALRegisterUvarMember(
    setup_file, STRING, 'S', REQUIRED,
    "Setup file generated by lalpulsar_WeaveSetup; the segment list, parameter-space metrics, and other required data. "
  );
  XLALRegisterUvarMember(
    output_file, STRING, 'o', REQUIRED,
    "Output file which stores all quantities computed by lalpulsar_Weave. "
  );
  //
  // - SFT input/generation and signal generation
  //
  lalUserVarHelpOptionSubsection = "SFT input/generation and signal generation";
  XLALRegisterUvarMember(
    sft_files, STRING, 'I', NODEFAULT,
    "Pattern matching the SFT files to be analysed. Possibilities are:\n"
    " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'"
  );
  XLALRegisterUvarMember(
    validate_sft_files, BOOLEAN, 'V', DEVELOPER,
    "Validate the checksums of the SFTs matched by " UVAR_STR( sft_files ) " before loading them into memory. "
  );
  XLALRegisterUvarMember(
    sft_timebase, REAL8, 't', NODEFAULT,
    "Generate SFTs with this timebase instead of loading from files. "
  );
  XLALRegisterUvarMember(
    sft_timestamps_files, STRINGVector, 'T', DEVELOPER,
    "Files containing timestamps for the generated SFTs; if not given, SFTs with contiguous timestamps are generated. "
    "Arguments correspond to the detectors in the setup file given by " UVAR_STR( setup_file )
    "; for example, if the setup file was created with " UVAR_STR( detectors ) " set to 'H1,L1', then an argument of "
    "'t1.txt,t2.txt' to this option will read H1 timestamps from the file 't1.txt', and L1 timestamps from the file 't2.txt'. "
    "The timebase of the generated SFTs is specified by " UVAR_STR( sft_timebase ) ". "
  );
  XLALRegisterUvarMember(
    sft_noise_sqrtSX, STRINGVector, 'p', NODEFAULT,
    "Inject fake Gaussian noise with these amplitude spectral densities [sqrt(Sh)] into the generated SFTs. "
    "Arguments correspond to the detectors in the setup file given by " UVAR_STR( setup_file )
    "; for example, if the setup file was created with " UVAR_STR( detectors ) " set to 'H1,L1', then an argument of "
    "'1.2,3.4' to this option will generate H1 SFTs with a noise sqrt(Sh) of 1.2, and L1 SFTs with a noise sqrt(Sh) of 3.4. "
  );
  XLALRegisterUvarMember(
    injections, STRINGVector, 'J', NODEFAULT,
    "%s", InjectionSourcesHelpString
  );
  XLALRegisterUvarMember(
    random_injection, REAL8Vector, 'H', NODEFAULT,
    "Inject a simulated signal, with phase parameters drawn randomly from the search parameter space, "
    "and with a randomly-generated polarisation angle (psi) and initial phase (phi0). "
    "The amplitude of the injection is specified by either 1 or 2 arguments to " UVAR_STR( random_injection ) ":\n"
    " - h0: generate a random cosine of inclination angle (cosi), then set aPlus=0.5*h0*(1+cosi^2), aCross=h0*cosi;\n"
    " - aPlus,aCross: use the given plus- and cross-polarisation amplitudes."
  );
 
  //
  // - Search parameter space
  //
  lalUserVarHelpOptionSubsection = "Search parameter space";
  XLALRegisterUvarMember(
    alpha, RAJRange, 'a', NODEFAULT,
    "Search parameter space in right ascension. "
    "If not specified, an all-sky search is performed; otherwise " UVAR_STR( delta ) " must also be specified. "
    "Range for a partial-sky search is limited to PI radians. "
  );
  XLALRegisterUvarMember(
    delta, DECJRange, 'd', NODEFAULT,
    "Search parameter space in declination. "
    "If not specified, an all-sky search is performed; otherwise " UVAR_STR( alpha ) " must also be specified. "
  );
  XLALRegisterUvarMember(
    sky_patch_count, UINT4, 'K', DEVELOPER,
    "Divide the entire sky into this number of ~equal-template-count patches. "
    "Requires " UVAR_STR( sky_patch_index ) ", mutually exclusive with " UVAR_STR2AND( alpha, delta ) ". "
  );
  XLALRegisterUvarMember(
    sky_patch_index, UINT4, 'k', DEVELOPER,
    "Search the sky patch given by this index, from zero to one less than " UVAR_STR( sky_patch_count ) ". "
    "Requires " UVAR_STR( sky_patch_count ) ", mutually exclusive with " UVAR_STR2AND( alpha, delta ) ". "
  );
  XLALRegisterUvarMember(
    freq, REAL8Range, 'f', REQUIRED,
    "Search parameter space in frequency, in Hertz. "
  );
  XLALRegisterUvarMember(
    freq_partitions, UINT4, 'F', DEVELOPER,
    "Internally divide the frequency parameter space into this number of ~equal-width partitions. "
  );
  XLALRegisterUvarMember(
    f1dot, REAL8Range, '1', OPTIONAL,
    "Search parameter space in first spindown, in Hertz/second. "
  );
  XLALRegisterUvarMember(
    f1dot_partitions, UINT4, '!', DEVELOPER,
    "Internally divide the first spindown parameter space into this number of ~equal-width partitions. "
  );
  XLALRegisterUvarMember(
    f2dot, REAL8Range, '2', OPTIONAL,
    "Search parameter space in second spindown, in Hertz/second^2. "
  );
  XLALRegisterUvarMember(
    f3dot, REAL8Range, '3', DEVELOPER,
    "Search parameter space in third spindown, in Hertz/second^3. "
  );
  XLALRegisterUvarMember(
    f4dot, REAL8Range, '4', DEVELOPER,
    "Search parameter space in fourth spindown, in Hertz/second^4. "
    "(Just in case a nearby supernova goes off!) "
  );
  XLALRegisterUvarMember(
    strict_spindown_bounds, BOOLEAN, '0', DEVELOPER,
    "Do not add padding to spindown parameter space bounds."
  );
  //
  // - Lattice tiling setup
  //
  lalUserVarHelpOptionSubsection = "Lattice tiling setup";
  XLALRegisterUvarMember(
    semi_max_mismatch, REAL8, 's', REQUIRED,
    "Maximum metric mismatch of the lattice tiling on which semicoherent quantities are computed, e.g. F-statistics averaged over segments. "
  );
  XLALRegisterUvarMember(
    coh_max_mismatch, REAL8, 'c', NODEFAULT,
    "Maximum metric mismatch of the per-segment lattice tilings on which coherent quantities are computed, e.g. coherent F-statistics. "
    "If the search setup contains only 1 segment, then this option must not be specified. "
  );
  XLALRegisterUvarMember(
    interpolation, BOOLEAN, 'i', OPTIONAL,
    "If TRUE, perform interpolation from the semicoherent lattice tiling to the per-segment coherent lattice tilings; "
    UVAR_STR( coh_max_mismatch ) " must also be specified in this case. "
    "If FALSE, turn off interpolation and use the same lattice tiling for both semicoherent and coherent computations; "
    UVAR_STR( coh_max_mismatch ) " must not be specified in this case. "
  );
  XLALRegisterUvarAuxDataMember(
    lattice, UserEnum, &TilingLatticeChoices, 'l', DEVELOPER,
    "Type of lattice used to generate the lattice tilings. "
  );
  XLALRegisterUvarMember(
    lattice_rand_offset, BOOLEAN, 'j', DEVELOPER,
    "Offset the physical parameter-space origin of the lattice tilings by a random fraction of the lattice step size. "
    "This is important when performing mismatch studies to ensure that the mismatch distribution is fully sampled. "
  );
  //
  // - F-statistic computation
  //
  lalUserVarHelpOptionSubsection = "F-statistic computation";
  XLALRegisterUvarAuxDataMember(
    Fstat_method, UserEnum, XLALFstatMethodChoices(), 'm', DEVELOPER,
    "Method used to calculate the F-statistic. "
  );
  XLALRegisterUvarMember(
    Fstat_run_med_window, UINT4, 'w', DEVELOPER,
    "Size of the running median window used to normalise SFTs and compute noise weight. "
  );
  XLALRegisterUvarMember(
    Fstat_assume_sqrtSX, STRINGVector, 'q', DEVELOPER,
    "Assume that the noise in the SFTs have known amplitude spectral densities [sqrt(Sh)], which are given by the arguments to "
    "this option, and normalise the SFTs by these given sqrt(Sh). "
    "Arguments correspond to the detectors in the setup file given by " UVAR_STR( setup_file )
    "; for example, if the setup file was created with " UVAR_STR( detectors ) " set to 'H1,L1', then an argument of "
    "'3.2,4.3' to this option will assume that H1 SFTs contain noise with a sqrt(Sh) of 3.2, and L1 SFTs contain noise with a sqrt(Sh) of 4.3. "
    "If this option is not given, the SFTs are normalised using noise sqrt(Sh) estimated from the SFTs themselves. "
  );
  XLALRegisterUvarMember(
    Fstat_Dterms, UINT4, 0, DEVELOPER,
    "Number of Dirichlet kernel terms to use in computing the F-statistic. May not be available for all F-statistic methods. "
  );
  XLALRegisterUvarAuxDataMember(
    Fstat_SSB_precision, UserEnum, &SSBprecisionChoices, 0, DEVELOPER,
    "Precision in calculating the barycentric transformation. "
  );
  //
  // Various statistics input arguments
  //
  lalUserVarHelpOptionSubsection = "Various statistics input arguments";
  XLALRegisterUvarMember(
    lrs_Fstar0sc, REAL8, 0, OPTIONAL,
    "(Semi-coherent) transition-scale parameter 'Fstar0sc' (=Nseg*Fstar0coh) for B_S/GL.. family of statistics."
  );
  XLALRegisterUvarMember(
    lrs_oLGX, STRINGVector, 0, OPTIONAL,
    "Per-detector line-vs-Gauss prior odds 'oLGX' (Defaults to oLGX=1/Ndet) for B_S/GL.. family of statistics."
  );
  XLALRegisterUvarMember(
    nc_2Fth, REAL8, 0, OPTIONAL,
    "Number count: per-segment 2F threshold value."
  );
  //
  // - Output control
  //
  lalUserVarHelpOptionSubsection = "Output control";
  XLALRegisterUvarMember(
    toplist_limit, UINT4, 'n', OPTIONAL,
    "Maximum number of candidates to return in an output toplist; if 0, all candidates are returned. "
  );
  XLALRegisterUvarAuxDataMember(
    toplists, UserFlag, &WeaveToplistChoices, 'L', OPTIONAL,
    "Sets which combination of toplists to return in the output file given by " UVAR_STR( output_file ) ":\n"
    "%s", WeaveToplistHelpString
  );
  XLALRegisterUvarMember(
    mean2F_hgrm, BOOLEAN, 0, DEVELOPER,
    "Output a histogram of all mean multi-Fstatistics computed by the search. "
  );
  XLALRegisterUvarAuxDataMember(
    extra_statistics, UserFlag, &WeaveStatisticChoices, 'E', OPTIONAL,
    "Sets which ('stage 0') extra statistics to compute and return in the output file given by " UVAR_STR( output_file ) ":\n"
    "%s", WeaveStatisticHelpString
  );
 
  XLALRegisterUvarAuxDataMember(
    recalc_statistics, UserFlag, &WeaveStatisticChoices, 'R', OPTIONAL,
    "Sets which extra *recalc* statistics to compute on final toplist without interpolation. See " UVAR_STR( extra_statistics ) " for statistics descriptions.\n"
  );
 
  XLALRegisterUvarMember(
    segment_info, BOOLEAN, 'M', DEVELOPER,
    "Output various information regarding the segment list, e.g. number of SFTs within each segment. "
  );
  //
  // - Checkpointing
  //
  lalUserVarHelpOptionSubsection = "Checkpointing";
  XLALRegisterUvarMember(
    ckpt_output_file, STRING, 'C', DEVELOPER,
    "File to which to periodically write checkpoints of output results. "
  );
  XLALRegisterUvarMember(
    ckpt_output_period, REAL8, 'z', DEVELOPER,
    "Write checkpoints of output results after this time period (in seconds) has elapsed. "
  );
  XLALRegisterUvarMember(
    ckpt_output_exit, REAL8, 0, DEVELOPER,
    "Write a checkpoint of output results after this fraction of the search has been completed, then exit. "
    "Arguments to this option must be in the range [0,1]. "
    "(This option is only really useful for testing the checkpointing feature.) "
  );
  //
  // - Esoterica
  //
  lalUserVarHelpOptionSubsection = "Esoterica";
  XLALRegisterUvarMember(
    rand_seed, UINT4, 'e', DEVELOPER,
    "Random seed used to initialise random number generators. "
  );
  XLALRegisterUvarMember(
    simulate_search, BOOLEAN, 0, DEVELOPER,
    "Simulate search; perform all search actions apart from computing any results. "
    "If SFT parameters (i.e. " UVAR_STR( sft_files ) " or " UVAR_STR( sft_timebase ) ") are supplied, simulate search with full memory allocation, i.e. with F-statistic input data, cached coherent results, etc. "
    "Otherwise, perform search with minimal memory allocation, i.e. do not allocate memory for any data or results. "
  );
  XLALRegisterUvarMember(
    time_search, BOOLEAN, 0, DEVELOPER,
    "Collect and output detailed timing information from various stages of the search pipeline. "
  );
  XLALRegisterUvarMember(
    cache_max_size, UINT4, 0, DEVELOPER,
    "Limit the size of the internal caches, used to store intermediate results, to this number of items per segment. "
    "If zero, the caches will grow in size to store all items that are still required. "
    "Has no effect when performing a fully-coherent single-segment search, or a non-interpolating search. "
  );
  XLALRegisterUvarMember(
    cache_all_gc, BOOLEAN, 0, DEVELOPER,
    "If TRUE, try to instead remove as many items as possible, provided that they are no longer required. "
    "If FALSE, whenever an item is added to the internal caches, at most one item that may no longer be required is removed. "
    "Has no effect when performing a fully-coherent single-segment search, or a non-interpolating search. "
  );
 
  // Parse user input
  XLAL_CHECK_MAIN( xlalErrno == 0, XLAL_EFUNC, "A call to XLALRegisterUvarMember() failed" );
  BOOLEAN should_exit = 0;
  XLAL_CHECK_MAIN( XLALUserVarReadAllInput( &should_exit, argc, argv, lalPulsarVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
  if ( should_exit ) {
    return EXIT_FAILURE;
  }
 
  // Check user input:
  //
  // - General
  //
 
  //
  // - SFT input/generation and signal generation
  //
  XLALUserVarCheck( &should_exit,
                    uvar->simulate_search || UVAR_SET2( sft_files, sft_timebase ) == 1,
                    "Exactly one of " UVAR_STR2OR( sft_files, sft_timebase ) " must be specified" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( sft_files ) || !UVAR_ALLSET3( sft_timebase, sft_timestamps_files, sft_noise_sqrtSX ),
                    UVAR_STR( sft_files ) " are mutually exclusive with " UVAR_STR3AND( sft_timebase, sft_timestamps_files, sft_noise_sqrtSX ) );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( validate_sft_files ) || UVAR_SET( sft_files ),
                    UVAR_STR( validate_sft_files ) " requires " UVAR_STR( sft_files ) );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( sft_timebase ) || uvar->sft_timebase > 0,
                    UVAR_STR( sft_timebase ) " must be strictly positive" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_ALLSET2( injections, random_injection ),
                    UVAR_STR( injections ) " and " UVAR_STR( random_injection ) " are mutually exclusive" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( random_injection ) || uvar->random_injection->length <= 2,
                    UVAR_STR( random_injection ) " must be passed either 1 or 2 arguments" );
  //
  // - Search parameter space
  //
  XLALUserVarCheck( &should_exit,
                    -LAL_PI_2 <= uvar->delta[0] && uvar->delta[1] <= LAL_PI_2,
                    UVAR_STR( delta ) " must be within range [-PI/2,PI/2]" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET2( sky_patch_count, sky_patch_index ) || !UVAR_ALLSET2( alpha, delta ),
                    UVAR_STR2AND( sky_patch_count, sky_patch_index ) " are mutually exclusive with " UVAR_STR2AND( alpha, delta ) );
  XLALUserVarCheck( &should_exit,
                    UVAR_SET2( sky_patch_count, sky_patch_index ) != 1,
                    UVAR_STR( sky_patch_count ) " requires " UVAR_STR( sky_patch_index ) " and vice versa" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( sky_patch_index ) || uvar->sky_patch_index < uvar->sky_patch_count,
                    UVAR_STR( sky_patch_index ) " must be positive and strictly less than " UVAR_STR( sky_patch_count ) );
  XLALUserVarCheck( &should_exit,
                    uvar->freq_partitions > 0,
                    UVAR_STR( freq_partitions ) " must be strictly positive" );
  XLALUserVarCheck( &should_exit,
                    uvar->f1dot_partitions > 0,
                    UVAR_STR( freq_partitions ) " must be strictly positive" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( f1dot ) || UVAR_SET( freq ),
                    UVAR_STR( freq ) " must be specified if " UVAR_STR( f1dot ) " is specified" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( f2dot ) || UVAR_SET( f1dot ),
                    UVAR_STR( f1dot ) " must be specified if " UVAR_STR( f2dot ) " is specified" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( f3dot ) || UVAR_SET( f2dot ),
                    UVAR_STR( f2dot ) " must be specified if " UVAR_STR( f3dot ) " is specified" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( f4dot ) || UVAR_SET( f3dot ),
                    UVAR_STR( f3dot ) " must be specified if " UVAR_STR( f4dot ) " is specified" );
  //
  // - Lattice tiling setup
  //
  XLALUserVarCheck( &should_exit,
                    uvar->semi_max_mismatch > 0,
                    UVAR_STR( semi_max_mismatch ) " must be strictly positive" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( coh_max_mismatch ) || uvar->coh_max_mismatch > 0,
                    UVAR_STR( coh_max_mismatch ) " must be strictly positive" );
  //
  // - F-statistic computation
  //
  XLALUserVarCheck( &should_exit,
                    uvar->Fstat_SSB_precision < SSBPREC_LAST,
                    UVAR_STR( Fstat_SSB_precision ) " must be in range [0,%u)", SSBPREC_LAST );
  //
  // - Output control
  //
 
  //
  // - Checkpointing
  //
  XLALUserVarCheck( &should_exit,
                    !UVAR_ALLSET2( ckpt_output_period, ckpt_output_exit ),
                    UVAR_STR2AND( ckpt_output_period, ckpt_output_exit ) " are mutually exclusive" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( ckpt_output_period ) || uvar->ckpt_output_period > 0,
                    UVAR_STR( ckpt_output_period ) " must be strictly positive" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_SET( ckpt_output_exit ) || ( 0 <= uvar->ckpt_output_exit && uvar->ckpt_output_exit <= 1 ),
                    UVAR_STR( ckpt_output_exit ) " must be in range [0,1]" );
  //
  // - Esoterica
  //
  XLALUserVarCheck( &should_exit,
                    !UVAR_ALLSET2( time_search, simulate_search ),
                    UVAR_STR2AND( time_search, simulate_search ) " are mutually exclusive" );
  XLALUserVarCheck( &should_exit,
                    !UVAR_ALLSET2( time_search, ckpt_output_file ),
                    UVAR_STR2AND( time_search, ckpt_output_file ) " are mutually exclusive" );
 
  // Exit if required
  if ( should_exit ) {
    return EXIT_FAILURE;
  }
  LogPrintf( LOG_NORMAL, "Parsed user input successfully\n" );
 
  // Allocate random number generator
  RandomParams *rand_par = XLALCreateRandomParams( uvar->rand_seed );
  XLAL_CHECK_MAIN( rand_par != NULL, XLAL_EFUNC );
 
  ////////// Load setup data //////////
 
  // Initialise setup data
  WeaveSetupData XLAL_INIT_DECL( setup );
 
  {
    // Open setup file
    LogPrintf( LOG_NORMAL, "Opening setup file '%s' for reading ...\n", uvar->setup_file );
    FITSFile *file = XLALFITSFileOpenRead( uvar->setup_file );
    XLAL_CHECK_MAIN( file != NULL, XLAL_EFUNC );
 
    // Read setup data
    XLAL_CHECK_MAIN( XLALWeaveSetupDataRead( file, &setup ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Close output file
    XLALFITSFileClose( file );
    LogPrintf( LOG_NORMAL, "Closed setup file '%s'\n", uvar->setup_file );
  }
 
  // Print reference time
  LogPrintf( LOG_NORMAL, "Setup file reference time = %" LAL_GPS_FORMAT "\n", LAL_GPS_PRINT( setup.ref_time ) );
 
  // Number of detectors and segments
  const UINT4 ndetectors = setup.detectors->length;
  const UINT4 nsegments = setup.segments->length;
 
  // Concatenate list of detector into a string
  char *setup_detectors_string = XLALConcatStringVector( setup.detectors, "," );
  XLAL_CHECK_MAIN( setup_detectors_string != NULL, XLAL_EFUNC );
 
  // Compute segment list range
  LIGOTimeGPS segments_start, segments_end;
  XLAL_CHECK_MAIN( XLALSegListRange( setup.segments, &segments_start, &segments_end ) == XLAL_SUCCESS, XLAL_EFUNC );
  LogPrintf( LOG_NORMAL, "Setup file segment list range = [%" LAL_GPS_FORMAT ", %" LAL_GPS_FORMAT "] GPS, segment count = %u\n", LAL_GPS_PRINT( segments_start ), LAL_GPS_PRINT( segments_end ), nsegments );
 
  ////////// Set up calculation of various requested output statistics //////////
 
  WeaveStatisticsParams *statistics_params = XLALCalloc( 1, sizeof( *statistics_params ) );
  XLAL_CHECK_MAIN( statistics_params != NULL, XLAL_ENOMEM );
  statistics_params->detectors = XLALCopyStringVector( setup.detectors );
  XLAL_CHECK_MAIN( statistics_params->detectors != NULL, XLAL_EFUNC );
  statistics_params->nsegments = nsegments;
  statistics_params->ref_time = setup.ref_time;
 
  //
  // Figure out which statistics need to be computed, and when, in order to
  // produce all the requested toplist-statistics in the "main loop" and all remaining
  // extra-statistics in the "completion loop" after the toplist has been computed
  // work out dependency-map for different statistics sets: toplist-ranking, output, total set of dependencies in main/completion loop ...
  XLAL_CHECK_MAIN( XLALWeaveStatisticsParamsSetDependencyMap( statistics_params, uvar->toplists, uvar->extra_statistics, uvar->recalc_statistics ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  // Prepare memory for coherent F-stat argument setups
  statistics_params->coh_input = XLALCalloc( nsegments, sizeof( statistics_params->coh_input[0] ) );
  XLAL_CHECK_MAIN( statistics_params->coh_input != NULL, XLAL_ENOMEM );
 
  // ---------- prepare setup for line-robust statistics if requested ----------
  if ( statistics_params->all_statistics_to_compute & ( WEAVE_STATISTIC_BSGL | WEAVE_STATISTIC_BSGLtL | WEAVE_STATISTIC_BtSGLtL ) ) {
    REAL4 *oLGX_p = NULL;
    REAL4 oLGX[PULSAR_MAX_DETECTORS];
    if ( uvar->lrs_oLGX != NULL ) {
      XLAL_CHECK_MAIN( uvar->lrs_oLGX->length == ndetectors, XLAL_EINVAL, "length(lrs-oLGX) = %d must equal number of detectors (%d)'\n", uvar->lrs_oLGX->length, ndetectors );
      XLAL_CHECK_MAIN( XLALParseLinePriors( &oLGX[0], uvar->lrs_oLGX ) == XLAL_SUCCESS, XLAL_EFUNC );
      oLGX_p = &oLGX[0];
    }
    const BOOLEAN useLogCorrection = 0;
    statistics_params->BSGL_setup = XLALCreateBSGLSetup( ndetectors, uvar->lrs_Fstar0sc, oLGX_p, useLogCorrection, nsegments );
    XLAL_CHECK_MAIN( statistics_params->BSGL_setup != NULL, XLAL_EFUNC );
  }
  // set number-count threshold
  statistics_params->nc_2Fth = uvar->nc_2Fth;
 
  // If asked to return a histogram of mean multi-Fstatistics, ensure they're computed
  if ( uvar->mean2F_hgrm ) {
    statistics_params->mainloop_statistics |= WEAVE_STATISTIC_MEAN2F;
  }
 
  ////////// Set up lattice tilings //////////
 
  // Check seach sky parameter space is consistent with the type of metric.
  LogPrintf( LOG_NORMAL, "Using metric for %s search\n", setup.metric_type );
  const BOOLEAN metric_is_directed = ( strcmp( setup.metric_type, "directed" ) == 0 );
  const BOOLEAN sky_is_a_point = ( uvar->alpha[0] == uvar->alpha[1] && uvar->delta[0] == uvar->delta[1] );
  XLALUserVarCheck( &should_exit,
                    !metric_is_directed || sky_is_a_point,
                    "The 'directed' metric type can only be used with a zero-area sky patch (i.e., --alpha and --delta must specify a single point)" );
 
  // Check interpolation/maximum mismatch options are consistent with the type of search being performed
  if ( nsegments == 1 ) {
    XLALUserVarCheck( &should_exit,
                      !UVAR_SET( coh_max_mismatch ),
                      UVAR_STR( coh_max_mismatch ) " must not be specified if setup file '%s' contains only 1 segment", uvar->setup_file );
    XLALUserVarCheck( &should_exit,
                      !UVAR_SET( interpolation ) || !uvar->interpolation,
                      UVAR_STR( interpolation ) " must either be FALSE or not specified if setup file '%s' contains only 1 segment", uvar->setup_file );
  } else if ( uvar->interpolation ) {
    XLALUserVarCheck( &should_exit,
                      UVAR_SET( coh_max_mismatch ),
                      UVAR_STR( coh_max_mismatch ) " must be specified if " UVAR_STR( interpolation ) " is true" );
  } else {
    XLALUserVarCheck( &should_exit,
                      !UVAR_SET( coh_max_mismatch ),
                      UVAR_STR( coh_max_mismatch ) " must not be set if " UVAR_STR( interpolation ) " is false" );
  }
  if ( should_exit ) {
    return EXIT_FAILURE;
  }
 
  // Decide which mismatches to use, depending of whether this is an interpolating or non-interpolating search
  const BOOLEAN interpolation = ( nsegments > 1 ) ? uvar->interpolation : 0;
  const double semi_max_mismatch = uvar->semi_max_mismatch;
  const double coh_max_mismatch = interpolation ? uvar->coh_max_mismatch : semi_max_mismatch;
  if ( nsegments == 1 ) {
    LogPrintf( LOG_NORMAL, "Performing a fully-coherent single-segment search with maximum (semicoherent) mismatch = %.15g\n", semi_max_mismatch );
  } else if ( !interpolation ) {
    LogPrintf( LOG_NORMAL, "Performing a non-interpolating search with maximum (semicoherent) mismatch = %.15g\n", semi_max_mismatch );
  } else {
    LogPrintf( LOG_NORMAL, "Performing an interpolating search with maximum semicoherent mismatch = %.15g, maximum coherent mismatch = %.15g\n", semi_max_mismatch, coh_max_mismatch );
  }
 
  // Scale metrics to fiducial frequency, given by maximum search frequency
  XLAL_CHECK_MAIN( XLALScaleSuperskyMetricsFiducialFreq( setup.metrics, uvar->freq[1] ) == XLAL_SUCCESS, XLAL_EFUNC );
  LogPrintf( LOG_NORMAL, "Metric fiducial frequency set to maximum search frequency = %.15g Hz\n", uvar->freq[1] );
 
  // Equalise metric frequency spacing, given the specified maximum mismatches
  XLAL_CHECK_MAIN( XLALEqualizeReducedSuperskyMetricsFreqSpacing( setup.metrics, coh_max_mismatch, semi_max_mismatch ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  // Store user input spindown ranges in an array for ease of use
  const double uvarspins[][2] = {
    { uvar->f1dot[0], uvar->f1dot[1] },
    { uvar->f2dot[0], uvar->f2dot[1] },
    { uvar->f3dot[0], uvar->f3dot[1] },
    { uvar->f4dot[0], uvar->f4dot[1] }
  };
 
  // Check that metrics computed in setup file have sufficient spindown dimensions to cover user input
  size_t nmetricspins = 0;
  XLAL_CHECK_MAIN( XLALSuperskyMetricsDimensions( setup.metrics, &nmetricspins ) == XLAL_SUCCESS, XLAL_EFUNC );
  const size_t nuvarspins = XLAL_NUM_ELEM( uvarspins );
  XLAL_CHECK_MAIN( nmetricspins <= nuvarspins, XLAL_EINVAL, "Number of spindowns from metrics (%zu) computed in setup file '%s' must be <= %zu", nmetricspins, uvar->setup_file, nuvarspins );
  const size_t ninputspins = UVAR_SET4( f1dot, f2dot, f3dot, f4dot );
  XLAL_CHECK_MAIN( ninputspins <= nmetricspins, XLAL_EINVAL, "Number of spindowns from user input (%zu) must be <= number of spindowns from metrics (%zu) computed in setup file '%s'", ninputspins, nmetricspins, uvar->setup_file );
 
  // Number of parameter-space dimensions: 2 for sky + 1 for frequency + 'nmetricspins' for spindowns
  const size_t ndim = 2 + 1 + nmetricspins;
 
  // Number of coherent parameter-space tilings
  // - If performing a fully-coherent search (i.e. of a single segment), we only need the semicoherent
  //   tiling; otherwise we need coherent tilings for each segment, plus the semicoherent tiling
  const size_t ncohtiles = ( nsegments > 1 ) ? nsegments : 0;
 
  // Total number of parameter-space tilings; always number of coherent tilings plus the semicoherent tiling
  const size_t ntiles = ncohtiles + 1;
 
  // Index of semicoherent tiling in arrays; always the last element
  const size_t isemi = ntiles - 1;
 
  // Create parameter-space tilings
  LatticeTiling *XLAL_INIT_DECL( tiling, [ntiles] );
  for ( size_t i = 0; i < ntiles; ++i ) {
    tiling[i] = XLALCreateLatticeTiling( ndim );
    XLAL_CHECK_MAIN( tiling[i] != NULL, XLAL_EFUNC );
  }
 
  // Create arrays to store the appropriate parameter-space metrics for each tiling
  gsl_matrix *XLAL_INIT_DECL( rssky_metric, [ntiles] );
  SuperskyTransformData *XLAL_INIT_DECL( rssky_transf, [ntiles] );
  for ( size_t i = 0; i < nsegments; ++i ) {
    rssky_metric[i] = setup.metrics->coh_rssky_metric[i];
    rssky_transf[i] = setup.metrics->coh_rssky_transf[i];
  }
  rssky_metric[isemi] = setup.metrics->semi_rssky_metric;
  rssky_transf[isemi] = setup.metrics->semi_rssky_transf;
 
  // Create arrays to store the range of physical coordinates covered by each tiling
  const PulsarDopplerParams *XLAL_INIT_DECL( min_phys, [ntiles] );
  const PulsarDopplerParams *XLAL_INIT_DECL( max_phys, [ntiles] );
 
  //
  // Set up semicoherent lattice tiling
  //
 
  // Set sky semicoherent parameter-space bounds
  // - Compute area of sky semicoherent parameter space for later output
  double semi_sky_area = 0.0;
  if ( UVAR_SET( sky_patch_count ) ) {
    XLAL_CHECK_MAIN( XLALSetSuperskyEqualAreaSkyBounds( tiling[isemi], rssky_metric[isemi], semi_max_mismatch, uvar->sky_patch_count, uvar->sky_patch_index ) == XLAL_SUCCESS, XLAL_EFUNC );
    LogPrintf( LOG_NORMAL, "Search sky parameter space sky patch = %u of %u\n", uvar->sky_patch_index, uvar->sky_patch_count );
    semi_sky_area = 4.0 * LAL_PI / uvar->sky_patch_count;
  } else {
    XLAL_CHECK_MAIN( XLALSetSuperskyPhysicalSkyBounds( tiling[isemi], rssky_metric[isemi], rssky_transf[isemi], uvar->alpha[0], uvar->alpha[1], uvar->delta[0], uvar->delta[1] ) == XLAL_SUCCESS, XLAL_EFUNC );
    LogPrintf( LOG_NORMAL, "Search sky parameter space right ascension = [%.15g, %.15g] rad\n", uvar->alpha[0], uvar->alpha[1] );
    LogPrintf( LOG_NORMAL, "Search sky parameter space declination = [%.15g, %.15g] rad\n", uvar->delta[0], uvar->delta[1] );
    semi_sky_area = ( uvar->alpha[1] - uvar->alpha[0] ) * ( sin( uvar->delta[1] ) - sin( uvar->delta[0] ) );
  }
 
  // Set frequency/spindown semicoherent parameter-space bounds
  XLAL_CHECK_MAIN( XLALSetSuperskyPhysicalSpinBound( tiling[isemi], rssky_transf[isemi], 0, uvar->freq[0], uvar->freq[1] ) == XLAL_SUCCESS, XLAL_EFUNC );
  LogPrintf( LOG_NORMAL, "Search frequency parameter space = [%.15g, %.15g] Hz\n", uvar->freq[0], uvar->freq[1] );
  for ( size_t s = 1; s <= nmetricspins; ++s ) {
    XLAL_CHECK_MAIN( XLALSetSuperskyPhysicalSpinBound( tiling[isemi], rssky_transf[isemi], s, uvarspins[s - 1][0], uvarspins[s - 1][1] ) == XLAL_SUCCESS, XLAL_EFUNC );
    XLAL_CHECK_MAIN( XLALSetSuperskyPhysicalSpinBoundPadding( tiling[isemi], rssky_transf[isemi], s, !uvar->strict_spindown_bounds ) == XLAL_SUCCESS, XLAL_EFUNC );
    LogPrintf( LOG_NORMAL, "Search %zu-order spindown parameter space = [%.15g, %.15g] Hz/s^%zu %s padding\n", s, uvarspins[s - 1][0], uvarspins[s - 1][1], s, uvar->strict_spindown_bounds ? "without" : "with" );
  }
 
  // Add random offsets to physical origin of semicoherent lattice tiling, if requested
  if ( UVAR_SET( lattice_rand_offset ) ) {
    XLAL_CHECK_MAIN( XLALSetLatticeTilingRandomOriginOffsets( tiling[isemi], rand_par ) == XLAL_SUCCESS, XLAL_EFUNC );
  }
 
  // Set semicoherent parameter-space lattice and metric
  XLAL_CHECK_MAIN( XLALSetTilingLatticeAndMetric( tiling[isemi], uvar->lattice, rssky_metric[isemi], semi_max_mismatch ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  // Print number of (tiled) parameter-space dimensions
  LogPrintf( LOG_NORMAL, "Number of (tiled) parameter-space dimensions = %zu (%zu)\n", ndim, XLALTiledLatticeTilingDimensions( tiling[isemi] ) );
 
  // Register callback to compute range of physical coordinates covered by semicoherent parameter space
  XLAL_CHECK_MAIN( XLALRegisterSuperskyLatticePhysicalRangeCallback( tiling[isemi], rssky_transf[isemi], &min_phys[isemi], &max_phys[isemi] ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  // Register callbacks to compute, for each coherent tiling, range of coherent reduced supersky coordinates which enclose semicoherent parameter space
  // - Arrays are of length 'ntiles' since 'ncohtiles' will be zero for a fully-coherent search
  const gsl_vector *XLAL_INIT_DECL( coh_min_rssky, [ntiles] );
  const gsl_vector *XLAL_INIT_DECL( coh_max_rssky, [ntiles] );
  for ( size_t i = 0; i < ncohtiles; ++i ) {
    XLAL_CHECK_MAIN( XLALRegisterSuperskyLatticeSuperskyRangeCallback( tiling[isemi], rssky_transf[isemi], rssky_transf[i], &coh_min_rssky[i], &coh_max_rssky[i] ) == XLAL_SUCCESS, XLAL_EFUNC );
  }
 
  // Iterate over semicoherent tiling and perform callback actions
  LogPrintf( LOG_NORMAL, "Setting up semicoherent lattice tiling ...\n" );
  XLAL_CHECK_MAIN( XLALPerformLatticeTilingCallbacks( tiling[isemi] ) == XLAL_SUCCESS, XLAL_EFUNC );
  LogPrintf( LOG_NORMAL, "Finished setting up semicoherent lattice tiling\n" );
 
  // Output number of points in semicoherent tiling
  {
    const LatticeTilingStats *stats = XLALLatticeTilingStatistics( tiling[isemi], ndim - 1 );
    XLAL_CHECK_MAIN( stats != NULL, XLAL_EFUNC );
    LogPrintf( LOG_NORMAL, "Number of semicoherent templates = %" LAL_UINT8_FORMAT "\n", stats->total_points );
  }
 
  // Get frequency spacing used by parameter-space tiling
  // - XLALEqualizeReducedSuperskyMetricsFreqSpacing() ensures this is the same for all segments
  const double dfreq = XLALLatticeTilingStepSize( tiling[isemi], ndim - 1 );
 
  //
  // Set up coherent lattice tilings
  //
  LogPrintf( LOG_NORMAL, "Setting up coherent lattice tilings ...\n" );
  for ( size_t i = 0; i < ncohtiles; ++i ) {
 
    // Set coherent parameter-space bounds which enclose semicoherent parameter space
    XLAL_CHECK_MAIN( XLALSetSuperskyRangeBounds( tiling[i], coh_min_rssky[i], coh_max_rssky[i] ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Add random offsets to physical origin of coherent lattice tiling, if requested
    if ( UVAR_SET( lattice_rand_offset ) ) {
      XLAL_CHECK_MAIN( XLALSetLatticeTilingRandomOriginOffsets( tiling[i], rand_par ) == XLAL_SUCCESS, XLAL_EFUNC );
    }
 
    // Ensure that coherent and semicoherent lattice tilings have the same tiled/non-tiled dimensions
    XLAL_CHECK_MAIN( XLALSetTiledLatticeDimensionsFromTiling( tiling[i], tiling[isemi] ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Set coherent parameter-space lattice and metric
    XLAL_CHECK_MAIN( XLALSetTilingLatticeAndMetric( tiling[i], uvar->lattice, rssky_metric[i], coh_max_mismatch ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Register callback to compute range of physical coordinates covered by coherent parameter space
    XLAL_CHECK_MAIN( XLALRegisterSuperskyLatticePhysicalRangeCallback( tiling[i], rssky_transf[i], &min_phys[i], &max_phys[i] ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Iterate over coherent tiling and perform callback actions
    XLAL_CHECK_MAIN( XLALPerformLatticeTilingCallbacks( tiling[i] ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  }
  LogPrintf( LOG_NORMAL, "Finished setting up coherent lattice tilings\n" );
 
  ////////// Load input data //////////
 
  // Load or generate SFTs, unless search is being simulated
  SFTCatalog *sft_catalog = NULL;
  if ( UVAR_SET( sft_files ) ) {
 
    // Load SFT catalog from files given by 'sft_files'
    LogPrintf( LOG_NORMAL, "Loading SFTs matching '%s' into catalog ...\n", uvar->sft_files );
    sft_catalog = XLALSFTdataFind( uvar->sft_files, NULL );
    XLAL_CHECK_MAIN( sft_catalog != NULL, XLAL_EFUNC );
    XLAL_CHECK_MAIN( sft_catalog->length > 0, XLAL_EFUNC );
    LogPrintf( LOG_NORMAL, "Loaded SFT catalog from SFTs matching '%s'\n", uvar->sft_files );
 
    // Validate checksums of SFTs, if requested
    if ( uvar->validate_sft_files ) {
      LogPrintf( LOG_NORMAL, "Validating SFTs ...\n" );
      BOOLEAN crc_check = 0;
      XLAL_CHECK_MAIN( XLALCheckCRCSFTCatalog( &crc_check, sft_catalog ) == XLAL_SUCCESS, XLAL_EFUNC );
      XLAL_CHECK_MAIN( crc_check, XLAL_EFUNC, "Failed to validate checksums of SFTs" );
      LogPrintf( LOG_NORMAL, "Finished validating SFTs\n" );
    }
 
  } else if ( UVAR_SET( sft_timebase ) ) {
 
    // Create timestamps for generated SFTs
    MultiLIGOTimeGPSVector *sft_timestamps = NULL;
    if ( UVAR_SET( sft_timestamps_files ) ) {
 
      // Check that the number of SFT timestamp files is consistent with the number of detectors
      XLAL_CHECK_MAIN( uvar->sft_timestamps_files->length == ndetectors, XLAL_EINVAL, "Number SFT timestamp files (%i) is inconsistent with number of detectors (%i) in setup file '%s'", uvar->sft_timestamps_files->length, ndetectors, uvar->setup_file );
 
      // Load SFT timestamps from files given by 'sft_timestamps_files'
      sft_timestamps = XLALReadMultiTimestampsFiles( uvar->sft_timestamps_files );
      XLAL_CHECK_MAIN( sft_timestamps != NULL, XLAL_EFUNC );
      for ( size_t i = 0; i < ndetectors; ++i ) {
        sft_timestamps->data[i]->deltaT = uvar->sft_timebase;
        LogPrintf( LOG_NORMAL, "Loaded SFT timestamps for detector '%s' from file '%s'\n", setup.detectors->data[i], uvar->sft_timestamps_files->data[i] );
      }
 
    } else {
 
      // Generate identical SFT timestamps for each detector, starting from beginning of segment list, with timebase given by 'sft_timebase'
      sft_timestamps = XLALMakeMultiTimestamps( segments_start, XLALGPSDiff( &segments_end, &segments_start ), uvar->sft_timebase, 0, ndetectors );
      XLAL_CHECK_MAIN( sft_timestamps != NULL, XLAL_EFUNC );
      LogPrintf( LOG_NORMAL, "Generated SFT timestamps for %i detectors, timebase = %.15g sec\n", ndetectors, uvar->sft_timebase );
 
    }
 
    // Generate SFT catalog for detectors 'sft_detectors' and timestamps 'sft_timestamps'
    sft_catalog = XLALMultiAddToFakeSFTCatalog( sft_catalog, setup.detectors, sft_timestamps );
    XLAL_CHECK_MAIN( sft_catalog != NULL, XLAL_EFUNC );
    XLAL_CHECK_MAIN( sft_catalog->length > 0, XLAL_EFUNC );
 
    // Cleanup
    XLALDestroyMultiTimestamps( sft_timestamps );
 
  }
  if ( sft_catalog != NULL ) {
 
    // Check that all SFT catalog detectors were included in setup file
    LALStringVector *sft_catalog_detectors = XLALListIFOsInCatalog( sft_catalog );
    XLAL_CHECK_MAIN( sft_catalog_detectors != NULL, XLAL_EFUNC );
    char *sft_catalog_detectors_string = XLALConcatStringVector( sft_catalog_detectors, "," );
    XLAL_CHECK_MAIN( sft_catalog_detectors_string != NULL, XLAL_EFUNC );
    XLAL_CHECK_MAIN( strcmp( sft_catalog_detectors_string, setup_detectors_string ) == 0, XLAL_EINVAL, "List of detectors '%s' in SFT catalog differs from list of detectors '%s' in setup file '%s'", sft_catalog_detectors_string, setup_detectors_string, uvar->setup_file );
 
    // Log number of SFTs, both in total and for each detector
    MultiSFTCatalogView *sft_catalog_view = XLALGetMultiSFTCatalogView( sft_catalog );
    XLAL_CHECK_MAIN( sft_catalog_view != NULL, XLAL_EINVAL );
    XLAL_CHECK_MAIN( sft_catalog_view->length > 0, XLAL_EFUNC );
    LogPrintf( LOG_NORMAL, "Using %u SFTs in total", sft_catalog->length );
    for ( size_t j = 0; j < sft_catalog_view->length; ++j ) {
      XLAL_CHECK_MAIN( sft_catalog_view->data[j].length > 0, XLAL_EFUNC );
      char *detector_name = XLALGetChannelPrefix( sft_catalog_view->data[j].data[0].header.name );
      XLAL_CHECK_MAIN( detector_name != NULL, XLAL_EFUNC );
      XLAL_CHECK_MAIN( XLALFindStringInVector( detector_name, setup.detectors ) >= 0, XLAL_EFAILED );
      LogPrintfVerbatim( LOG_NORMAL, ", %u SFTs from detector '%s'", sft_catalog_view->data[j].length, detector_name );
      XLALFree( detector_name );
    }
    LogPrintfVerbatim( LOG_NORMAL, "\n" );
 
    // Cleanup
    XLALDestroyStringVector( sft_catalog_detectors );
    XLALFree( sft_catalog_detectors_string );
    XLALDestroyMultiSFTCatalogView( sft_catalog_view );
 
  }
 
  // Decide on search simulation level
  const WeaveSimulationLevel simulation_level = uvar->simulate_search ? ( WEAVE_SIMULATE | ( sft_catalog == NULL ? WEAVE_SIMULATE_MIN_MEM : 0 ) ) : 0;
  if ( simulation_level & WEAVE_SIMULATE ) {
    LogPrintf( LOG_NORMAL, "Simulating search with %s memory usage; no results will be computed\n", simulation_level & WEAVE_SIMULATE_MIN_MEM ? "minimal" : "full" );
  }
 
  // Set up injections
  PulsarParamsVector *injections = NULL;
  if ( UVAR_SET( injections ) ) {
 
    // Parse signal injection string
    injections = XLALPulsarParamsFromUserInput( uvar->injections, &setup.ref_time );
    XLAL_CHECK_MAIN( injections != NULL, XLAL_EFUNC );
 
  } else if ( UVAR_SET( random_injection ) ) {
 
    // Create injection parameters vector
    injections = XLALCreatePulsarParamsVector( 1 );
    XLAL_CHECK_MAIN( injections != NULL, XLAL_EFUNC );
    PulsarParams *const inj = &injections->data[0];
    strcpy( inj->name, "random_injection" );
 
    // Draw random phase parameters from search parameter space
    {
      double XLAL_INIT_DECL( random_point, [ndim] );
      gsl_matrix_view random_point_matrix = gsl_matrix_view_array( random_point, ndim, 1 );
      XLAL_CHECK_MAIN( XLALRandomLatticeTilingPoints( tiling[isemi], 0, rand_par, &random_point_matrix.matrix ) == XLAL_SUCCESS, XLAL_EFUNC );
      gsl_vector_view random_point_vector = gsl_vector_view_array( random_point, ndim );
      XLAL_CHECK_MAIN( XLALConvertSuperskyToPhysicalPoint( &inj->Doppler, &random_point_vector.vector, NULL, rssky_transf[isemi] ) == XLAL_SUCCESS, XLAL_EFUNC );
    }
 
    // Randomly generate and set amplitude parameters
    inj->Amp.psi = LAL_TWOPI * XLALUniformDeviate( rand_par );
    inj->Amp.phi0 = LAL_TWOPI * XLALUniformDeviate( rand_par );
    if ( uvar->random_injection->length == 1 ) {
      const REAL8 h0 = uvar->random_injection->data[0];
      const REAL8 cosi = -1 + 2 * XLALUniformDeviate( rand_par );
      inj->Amp.aPlus = 0.5 * h0 * ( 1.0 + cosi * cosi );
      inj->Amp.aCross = h0 * cosi;
    } else {
      inj->Amp.aPlus = uvar->random_injection->data[0];
      inj->Amp.aCross = uvar->random_injection->data[1];
    }
 
  }
 
  // Set F-statistic optional arguments
  Fstat_opt_args.randSeed = uvar->rand_seed;
  Fstat_opt_args.SSBprec = uvar->Fstat_SSB_precision;
  Fstat_opt_args.Dterms = uvar->Fstat_Dterms;
  Fstat_opt_args.runningMedianWindow = uvar->Fstat_run_med_window;
  Fstat_opt_args.FstatMethod = uvar->Fstat_method;
  Fstat_opt_args.injectSources = injections;
  Fstat_opt_args.prevInput = NULL;
  Fstat_opt_args.collectTiming = uvar->time_search;
 
  // Load input data required for computing coherent results
  const LALStringVector *sft_noise_sqrtSX = UVAR_SET( sft_noise_sqrtSX ) ? uvar->sft_noise_sqrtSX : NULL;
  const LALStringVector *Fstat_assume_sqrtSX = UVAR_SET( Fstat_assume_sqrtSX ) ? uvar->Fstat_assume_sqrtSX : NULL;
  LogPrintf( LOG_NORMAL, "Loading input data for coherent results ...\n" );
  XLAL_INIT_MEM( statistics_params->n2F_det );
  for ( size_t i = 0; i < nsegments; ++i ) {
    statistics_params->coh_input[i] = XLALWeaveCohInputCreate( setup.detectors, simulation_level, sft_catalog, i, &setup.segments->segs[i], min_phys[i], max_phys[i], dfreq, setup.ephemerides, sft_noise_sqrtSX, Fstat_assume_sqrtSX, &Fstat_opt_args, statistics_params, 0 );
    XLAL_CHECK_MAIN( statistics_params->coh_input[i] != NULL, XLAL_EFUNC );
  }
  if ( !( simulation_level & WEAVE_SIMULATE_MIN_MEM ) && ( statistics_params->mainloop_statistics & WEAVE_STATISTIC_COH2F_DET ) ) {
    for ( size_t i = 0; i < ndetectors; ++i ) {
      XLAL_CHECK_MAIN( 0 < statistics_params->n2F_det[i] && statistics_params->n2F_det[i] <= nsegments, XLAL_EFAILED, "Invalid number of per-detector F-statistics (%u) for detector '%s'", statistics_params->n2F_det[i], setup.detectors->data[i] );
    }
  }
 
  LogPrintf( LOG_NORMAL, "Finished loading input data for coherent results\n" );
 
  // Create caches to store intermediate results from coherent parameter-space tilings
  // - If no interpolation, caching is not required so reduce maximum cache size to 1
  WeaveCache *XLAL_INIT_DECL( coh_cache, [nsegments] );
  for ( size_t i = 0; i < nsegments; ++i ) {
    const size_t cache_max_size = interpolation ? uvar->cache_max_size : 1;
    const BOOLEAN cache_all_gc = interpolation ? uvar->cache_all_gc : 0;
    coh_cache[i] = XLALWeaveCacheCreate( tiling[i], interpolation, rssky_transf[i], rssky_transf[isemi], statistics_params->coh_input[i], cache_max_size, cache_all_gc );
    XLAL_CHECK_MAIN( coh_cache[i] != NULL, XLAL_EFUNC );
  }
 
  ////////// Perform search //////////
 
  // Create iterator over the main loop search parameter space
  WeaveSearchIterator *main_loop_itr = XLALWeaveMainLoopSearchIteratorCreate( tiling[isemi], uvar->freq_partitions, uvar->f1dot_partitions );
  XLAL_CHECK_MAIN( main_loop_itr != NULL, XLAL_EFUNC );
 
  // Create storage for cache queries for coherent results in each segment
  WeaveCacheQueries *queries = XLALWeaveCacheQueriesCreate( tiling[isemi], rssky_transf[isemi], dfreq, nsegments, uvar->freq_partitions );
  XLAL_CHECK_MAIN( queries != NULL, XLAL_EFUNC );
 
  // Pointer to final semicoherent results
  WeaveSemiResults *semi_res = NULL;
 
  // Create output results structure
  WeaveOutputResults *out = XLALWeaveOutputResultsCreate( &setup.ref_time, ninputspins, statistics_params, uvar->toplist_limit, uvar->mean2F_hgrm );
  XLAL_CHECK_MAIN( out != NULL, XLAL_EFUNC );
 
  // Create search timing structure
  WeaveSearchTiming *tim = XLALWeaveSearchTimingCreate( uvar->time_search, statistics_params );
  XLAL_CHECK_MAIN( tim != NULL, XLAL_EFUNC );
 
  // Number of times output results have been restored from a checkpoint
  UINT4 ckpt_output_count = 0;
 
  // Try to restore output results from a checkpoint file, if given
  if ( UVAR_SET( ckpt_output_file ) ) {
 
    // Try to open output checkpoint file
    LogPrintf( LOG_NORMAL, "Trying to open output checkpoint file '%s' for reading ...\n", uvar->ckpt_output_file );
    int errnum = 0;
    FITSFile *file = NULL;
    XLAL_TRY( file = XLALFITSFileOpenRead( uvar->ckpt_output_file ), errnum );
    if ( errnum == XLAL_ENOENT ) {
      LogPrintf( LOG_NORMAL, "Output checkpoint file '%s' does not exist; no checkpoint will be loaded\n", uvar->ckpt_output_file );
    } else {
      XLAL_CHECK_MAIN( errnum == 0 && file != NULL, XLAL_EFUNC );
      LogPrintf( LOG_NORMAL, "Output checkpoint file '%s' exists; checkpoint will be loaded\n", uvar->ckpt_output_file );
 
      // Read number of times output results have been restored from a checkpoint
      XLAL_CHECK_MAIN( XLALFITSHeaderReadUINT4( file, "ckptcnt", &ckpt_output_count ) == XLAL_SUCCESS, XLAL_EFUNC );
      XLAL_CHECK_MAIN( ckpt_output_count > 0, XLAL_EIO, "Invalid output checkpoint file '%s'", uvar->ckpt_output_file );
 
      // Read maximum size of toplists
      UINT4 toplist_limit = 0;
      XLAL_CHECK( XLALFITSHeaderReadUINT4( file, "toplimit", &toplist_limit ) == XLAL_SUCCESS, XLAL_EFUNC );
 
      // Read output results
      XLAL_CHECK_MAIN( XLALWeaveOutputResultsReadAppend( file, &out, toplist_limit ) == XLAL_SUCCESS, XLAL_EFUNC, "Invalid output checkpoint file '%s'", uvar->ckpt_output_file );
 
      // Restore state of main loop iterator
      XLAL_CHECK_MAIN( XLALWeaveSearchIteratorRestore( main_loop_itr, file ) == XLAL_SUCCESS, XLAL_EFUNC, "Invalid output checkpoint file '%s'", uvar->ckpt_output_file );
 
      // Close output checkpoint file
      XLALFITSFileClose( file );
      LogPrintf( LOG_NORMAL, "Closed output checkpoint file '%s'\n", uvar->ckpt_output_file );
 
    }
 
  }
 
  // Start timing main search loop
  XLAL_CHECK_MAIN( XLALWeaveSearchTimingStart( tim ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  // Elapsed wall time at which search was last checkpointed
  double wall_ckpt_elapsed = 0;
 
  // Elapsed wall time at which progress was last printed, and interval at which to print progress
  double wall_prog_elapsed = 0;
  double wall_prog_period = 5.0;
 
  // Whether to print predicted remaining time, and previous prediction for total elapsed time
  BOOLEAN wall_prog_remain_print = 0;
  double wall_prog_total_prev = 0;
 
  // Print initial progress
  LogPrintf( LOG_NORMAL, "Starting main loop at %.3g%% complete, peak memory %.1fMB\n", XLALWeaveSearchIteratorProgress( main_loop_itr ), XLALGetPeakHeapUsageMB() );
 
  // Begin main loop
  BOOLEAN search_complete = 0;
  while ( !search_complete ) {
 
    // Switch timing section
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_OTHER, WEAVE_SEARCH_TIMING_ITER ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Get next semicoherent frequency block
    // - Exit main loop if iteration is complete
    // - Expire cache items if requested by iterator
    BOOLEAN expire_cache = 0;
    UINT8 semi_index = 0;
    const gsl_vector *semi_rssky = NULL;
    INT4 semi_left = 0;
    INT4 semi_right = 0;
    UINT4 freq_partition_index = 0;
    XLAL_CHECK( XLALWeaveSearchIteratorNext( main_loop_itr, &search_complete, &expire_cache, &semi_index, &semi_rssky, &semi_left, &semi_right, &freq_partition_index ) == XLAL_SUCCESS, XLAL_EFUNC );
    if ( search_complete ) {
      XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_ITER, WEAVE_SEARCH_TIMING_OTHER ) == XLAL_SUCCESS, XLAL_EFUNC );
      break;
    } else if ( expire_cache ) {
      for ( size_t i = 0; i < nsegments; ++i ) {
        XLAL_CHECK_MAIN( XLALWeaveCacheExpire( coh_cache[i] ) == XLAL_SUCCESS, XLAL_EFUNC );
      }
    }
 
    // Switch timing section
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_ITER, WEAVE_SEARCH_TIMING_QUERY ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Initialise cache queries
    XLAL_CHECK_MAIN( XLALWeaveCacheQueriesInit( queries, semi_index, semi_rssky, semi_left, semi_right, freq_partition_index ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Query for coherent results for each segment
    for ( size_t i = 0; i < nsegments; ++i ) {
      XLAL_CHECK_MAIN( XLALWeaveCacheQuery( coh_cache[i], queries, i ) == XLAL_SUCCESS, XLAL_EFUNC );
    }
 
    // Finalise cache queries
    PulsarDopplerParams XLAL_INIT_DECL( semi_phys );
    UINT4 semi_nfreqs = 0;
    XLAL_CHECK_MAIN( XLALWeaveCacheQueriesFinal( queries, &semi_phys, &semi_nfreqs ) == XLAL_SUCCESS, XLAL_EFUNC );
    if ( semi_nfreqs == 0 ) {
      XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_QUERY, WEAVE_SEARCH_TIMING_OTHER ) == XLAL_SUCCESS, XLAL_EFUNC );
      continue;
    }
 
    // Switch timing section
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_QUERY, WEAVE_SEARCH_TIMING_COH ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Retrieve coherent results from each segment
    const WeaveCohResults *XLAL_INIT_DECL( coh_res, [nsegments] );
    UINT8 XLAL_INIT_DECL( coh_index, [nsegments] );
    UINT4 XLAL_INIT_DECL( coh_offset, [nsegments] );
    for ( size_t i = 0; i < nsegments; ++i ) {
      XLAL_CHECK_MAIN( XLALWeaveCacheRetrieve( coh_cache[i], queries, i, &coh_res[i], &coh_index[i], &coh_offset[i], tim ) == XLAL_SUCCESS, XLAL_EFUNC );
      XLAL_CHECK_MAIN( coh_res[i] != NULL, XLAL_EFUNC );
    }
 
    // Switch timing section
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_COH, WEAVE_SEARCH_TIMING_SEMISEG ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Initialise semicoherent results
    XLAL_CHECK_MAIN( XLALWeaveSemiResultsInit( &semi_res, simulation_level, ndetectors, nsegments, semi_index, &semi_phys, dfreq, semi_nfreqs, statistics_params ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Add coherent results to semicoherent results
    XLAL_CHECK_MAIN( XLALWeaveSemiResultsComputeSegs( semi_res, nsegments, coh_res, coh_index, coh_offset, tim ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Switch timing section
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_SEMISEG, WEAVE_SEARCH_TIMING_SEMI ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Compute all toplist-ranking semicoherent results
    XLAL_CHECK_MAIN( XLALWeaveSemiResultsComputeMain( semi_res, tim ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Switch timing section
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_SEMI, WEAVE_SEARCH_TIMING_OUTPUT ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Add semicoherent results to output
    XLAL_CHECK_MAIN( XLALWeaveOutputResultsAdd( out, semi_res, semi_nfreqs ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Switch timing section
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_OUTPUT, WEAVE_SEARCH_TIMING_OTHER ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Main iterator percentage complete
    const REAL4 prog_per_cent = XLALWeaveSearchIteratorProgress( main_loop_itr );
 
    // Current elapsed wall and CPU times
    double wall_elapsed = 0, cpu_elapsed = 0;
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingElapsed( tim, &wall_elapsed, &cpu_elapsed ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Print iteration progress, if required
    if ( wall_elapsed - wall_prog_elapsed >= wall_prog_period ) {
 
      // Print progress
      LogPrintf( LOG_NORMAL, "%s at %.3g%% complete", simulation_level & WEAVE_SIMULATE ? "Simulation" : "Search", prog_per_cent );
 
      // Print elapsed time
      LogPrintfVerbatim( LOG_NORMAL, ", elapsed %.1f sec", wall_elapsed );
 
      // Print remaining time, if it can be reliably predicted
      const double wall_prog_remain = XLALWeaveSearchIteratorRemainingTime( main_loop_itr, wall_elapsed );
      const double wall_prog_total = wall_elapsed + wall_prog_remain;
      if ( wall_prog_remain_print || fabs( wall_prog_total - wall_prog_total_prev ) <= 0.1 * wall_prog_total_prev ) {
        LogPrintfVerbatim( LOG_NORMAL, ", remaining ~%.1f sec", wall_prog_remain );
        wall_prog_remain_print = 1;   // Always print remaining time once it can be reliably predicted
      } else {
        wall_prog_total_prev = wall_prog_total;
      }
 
      // Print CPU usage
      LogPrintfVerbatim( LOG_NORMAL, ", CPU %.1f%%", 100.0 * cpu_elapsed / wall_elapsed );
 
      // Print memory usage
      LogPrintfVerbatim( LOG_NORMAL, ", peak memory %.1fMB", XLALGetPeakHeapUsageMB() );
#ifdef LALPULSAR_CUDA_ENABLED
      if ( Fstat_opt_args.FstatMethod == FMETHOD_RESAMP_CUDA ) {
        size_t CUDA_free_mem_B = 0;
        size_t CUDA_tot_mem_B = 0;
        XLAL_CHECK_MAIN( cudaMemGetInfo( &CUDA_free_mem_B, &CUDA_tot_mem_B ) == cudaSuccess, XLAL_EERR );
        XLAL_CHECK_MAIN( CUDA_free_mem_B <= CUDA_tot_mem_B, XLAL_EERR );
        const size_t CUDA_used_mem_B = CUDA_tot_mem_B - CUDA_free_mem_B;
        const REAL4 CUDA_used_mem_MB = CUDA_used_mem_B / ( 1024.0 * 1024.0 );
        const REAL4 CUDA_tot_mem_MB = CUDA_tot_mem_B / ( 1024.0 * 1024.0 );
        LogPrintfVerbatim( LOG_NORMAL, ", CUDA memory %.1f/%.1fMB", CUDA_used_mem_MB, CUDA_tot_mem_MB );
      }
#endif
 
      // Print mean maximum size obtained by caches
      {
        REAL4 cache_mean_max_size = 0;
        XLAL_CHECK_MAIN( XLALWeaveGetCacheMeanMaxSize( &cache_mean_max_size, nsegments, coh_cache ) == XLAL_SUCCESS, XLAL_EFUNC );
        LogPrintfVerbatim( LOG_NORMAL, ", cache max ~%0.1f", cache_mean_max_size );
      }
 
      // Finish progress printing
      LogPrintfVerbatim( LOG_NORMAL, "\n" );
 
      // Update elapsed wall time at which progress was last printed, and increase interval at which to print progress
      wall_prog_elapsed = wall_elapsed;
      wall_prog_period = GSL_MIN( 1200, wall_prog_period * 1.5 );
 
    }
 
    // Checkpoint output results, if required
    if ( UVAR_SET( ckpt_output_file ) ) {
 
      // Switch timing section
      XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_OTHER, WEAVE_SEARCH_TIMING_CKPT ) == XLAL_SUCCESS, XLAL_EFUNC );
 
      // Decide whether to checkpoint output results
      const BOOLEAN do_ckpt_output_period = UVAR_SET( ckpt_output_period ) && wall_elapsed - wall_ckpt_elapsed >= uvar->ckpt_output_period;
      const BOOLEAN do_ckpt_output_exit = UVAR_SET( ckpt_output_exit ) && prog_per_cent >= 100.0 * uvar->ckpt_output_exit;
      if ( do_ckpt_output_period || do_ckpt_output_exit ) {
 
        // Open output checkpoint file
        FITSFile *file = XLALFITSFileOpenWrite( uvar->ckpt_output_file );
        XLAL_CHECK_MAIN( file != NULL, XLAL_EFUNC );
        XLAL_CHECK_MAIN( XLALFITSFileWriteVCSInfo( file, lalPulsarVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
        XLAL_CHECK_MAIN( XLALFITSFileWriteUVarCmdLine( file ) == XLAL_SUCCESS, XLAL_EFUNC );
 
        // Write number of times output results have been restored from a checkpoint
        ++ckpt_output_count;
        XLAL_CHECK_MAIN( XLALFITSHeaderWriteUINT4( file, "ckptcnt", ckpt_output_count, "number of checkpoints" ) == XLAL_SUCCESS, XLAL_EFUNC );
 
        // Write output results
        XLAL_CHECK_MAIN( XLALWeaveOutputResultsWrite( file, out ) == XLAL_SUCCESS, XLAL_EFUNC );
 
        // Save state of main loop iterator
        XLAL_CHECK_MAIN( XLALWeaveSearchIteratorSave( main_loop_itr, file ) == XLAL_SUCCESS, XLAL_EFUNC );
 
        // Close output checkpoint file
        XLALFITSFileClose( file );
 
        // Print progress
        LogPrintf( LOG_NORMAL, "Wrote output checkpoint to file '%s' at %.3g%% complete, elapsed %.1f sec\n", uvar->ckpt_output_file, prog_per_cent, wall_elapsed );
 
        // Exit main loop, if checkpointing was triggered by 'do_ckpt_output_exit'
        if ( do_ckpt_output_exit ) {
          LogPrintf( LOG_NORMAL, "Exiting main seach loop after writing output checkpoint\n" );
          XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_CKPT, WEAVE_SEARCH_TIMING_OTHER ) == XLAL_SUCCESS, XLAL_EFUNC );
          break;
        }
 
        // Update elapsed wall time at which search was last checkpointed
        wall_ckpt_elapsed = wall_elapsed;
 
      }
 
      // Switch timing section
      XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_CKPT, WEAVE_SEARCH_TIMING_OTHER ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    }
 
  }   // End of main loop
 
  // Clear all cache items from memory
  for ( size_t i = 0; i < nsegments; ++i ) {
    XLAL_CHECK_MAIN( XLALWeaveCacheClear( coh_cache[i] ) == XLAL_SUCCESS, XLAL_EFUNC );
  }
 
  // Print progress
  double wall_main = 0, cpu_main = 0;
  XLAL_CHECK_MAIN( XLALWeaveSearchTimingElapsed( tim, &wall_main, &cpu_main ) == XLAL_SUCCESS, XLAL_EFUNC );
  LogPrintf( LOG_NORMAL, "Finished main loop at %.3g%% complete, main-loop time %.1f sec, CPU %.1f%%, peak memory %.1fMB\n", XLALWeaveSearchIteratorProgress( main_loop_itr ), wall_main, 100.0 * cpu_main / wall_main, XLALGetPeakHeapUsageMB() );
 
  // Prepare completion-loop calculations:
  // if any 'recalc' (= 'stage 1') statistics have been requested: we'll need 'Demod' Fstatistic setups
  // so if the 'stage 0' calculation used Demod => nothing to do, if it used 'Resamp' then re-compute the setups
  if ( ( uvar->recalc_statistics != WEAVE_STATISTIC_NONE ) ) {
    statistics_params->coh_input_recalc = XLALCalloc( nsegments, sizeof( statistics_params->coh_input_recalc[0] ) );
    XLAL_CHECK_MAIN( statistics_params->coh_input_recalc != NULL, XLAL_ENOMEM );
    FstatOptionalArgs Fstat_opt_args_recalc = Fstat_opt_args;
    Fstat_opt_args_recalc.FstatMethod = FMETHOD_DEMOD_BEST;
    Fstat_opt_args_recalc.prevInput = NULL;
    XLAL_INIT_MEM( statistics_params->n2F_det );
    for ( size_t i = 0; i < nsegments; ++i ) {
      statistics_params->coh_input_recalc[i] = XLALWeaveCohInputCreate( setup.detectors, simulation_level, sft_catalog, i, &setup.segments->segs[i], min_phys[i], max_phys[i], 0, setup.ephemerides, sft_noise_sqrtSX, Fstat_assume_sqrtSX, &Fstat_opt_args_recalc, statistics_params, 1 );
      XLAL_CHECK_MAIN( statistics_params->coh_input_recalc[i] != NULL, XLAL_EFUNC );
    }
    if ( !( simulation_level & WEAVE_SIMULATE_MIN_MEM ) && ( statistics_params->completionloop_statistics[1] & WEAVE_STATISTIC_COH2F_DET ) ) {
      for ( size_t i = 0; i < ndetectors; ++i ) {
        XLAL_CHECK_MAIN( 0 < statistics_params->n2F_det[i] && statistics_params->n2F_det[i] <= nsegments, XLAL_EFAILED, "Invalid number of per-detector F-statistics (%u) for detector '%s'", statistics_params->n2F_det[i], setup.detectors->data[i] );
      }
    }
  }
 
  // Switch timing section
  XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_OTHER, WEAVE_SEARCH_TIMING_CMPL ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  // Completion loop: compute all extra statistics that weren't required in the main loop
  XLAL_CHECK_MAIN( XLALWeaveOutputResultsCompletionLoop( out ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  // Switch timing section
  XLAL_CHECK_MAIN( XLALWeaveSearchTimingSection( tim, WEAVE_SEARCH_TIMING_CMPL, WEAVE_SEARCH_TIMING_OTHER ) == XLAL_SUCCESS, XLAL_EFUNC );
 
  // Stop timing main search loop, and get total wall and CPU times
  double wall_total = 0, cpu_total = 0;
  XLAL_CHECK_MAIN( XLALWeaveSearchTimingStop( tim, &wall_total, &cpu_total ) == XLAL_SUCCESS, XLAL_EFUNC );
  LogPrintf( LOG_NORMAL, "Finished completion-loop, total time %.1f sec, CPU %.1f%%, peak memory %.1fMB\n", wall_total, 100.0 * cpu_total / wall_total, XLALGetPeakHeapUsageMB() );
 
  ////////// Output search results //////////
 
  if ( search_complete ) {
 
    // Open output file
    LogPrintf( LOG_NORMAL, "Opening output file '%s' for writing ...\n", uvar->output_file );
    FITSFile *file = XLALFITSFileOpenWrite( uvar->output_file );
    XLAL_CHECK_MAIN( file != NULL, XLAL_EFUNC );
    XLAL_CHECK_MAIN( XLALFITSFileWriteVCSInfo( file, lalPulsarVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
    XLAL_CHECK_MAIN( XLALFITSFileWriteUVarCmdLine( file ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write number of times output results were restored from a checkpoint
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteUINT4( file, "numckpt", ckpt_output_count, "number of checkpoints" ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write list of detectors
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteStringVector( file, "detect", setup.detectors, "setup detectors" ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write number of segments
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteUINT4( file, "nsegment", nsegments, "number of segments" ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write frequency spacing
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteREAL8( file, "dfreq", dfreq, "frequency spacing" ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write semicoherent parameter-space bounds
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteREAL8( file, "semiparam skyarea [sr]", semi_sky_area, "area of sky parameter space" ) == XLAL_SUCCESS, XLAL_EFUNC );
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteREAL8( file, "semiparam minfreq [Hz]", uvar->freq[0], "minimum frequency range" ) == XLAL_SUCCESS, XLAL_EFUNC );
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteREAL8( file, "semiparam maxfreq [Hz]", uvar->freq[1], "maximum frequency range" ) == XLAL_SUCCESS, XLAL_EFUNC );
    for ( size_t s = 1; s <= ninputspins; ++s ) {
      char keyword[64];
      char comment[64];
      snprintf( keyword, sizeof( keyword ), "semiparam minf%zudot [Hz/s^%zu]", s, s );
      snprintf( comment, sizeof( comment ), "minimum %zu-order spindown range", s );
      XLAL_CHECK_MAIN( XLALFITSHeaderWriteREAL8( file, keyword, uvarspins[s - 1][0], comment ) == XLAL_SUCCESS, XLAL_EFUNC );
      snprintf( keyword, sizeof( keyword ), "semiparam maxf%zudot [Hz/s^%zu]", s, s );
      snprintf( comment, sizeof( comment ), "maximum %zu-order spindown range", s );
      XLAL_CHECK_MAIN( XLALFITSHeaderWriteREAL8( file, keyword, uvarspins[s - 1][1], comment ) == XLAL_SUCCESS, XLAL_EFUNC );
    }
 
    // Write cumulative number of semicoherent templates in each dimension
    for ( size_t i = 0; i < ndim; ++i ) {
      char keyword[64];
      const LatticeTilingStats *semi_stats = XLALLatticeTilingStatistics( tiling[isemi], i );
      XLAL_CHECK_MAIN( semi_stats != NULL, XLAL_EFUNC );
      XLAL_CHECK_MAIN( semi_stats->name != NULL, XLAL_EFUNC );
      XLAL_CHECK_MAIN( semi_stats->total_points > 0, XLAL_EFUNC );
      snprintf( keyword, sizeof( keyword ), "nsemitmpl %s", semi_stats->name );
      XLAL_CHECK_MAIN( XLALFITSHeaderWriteUINT8( file, keyword, semi_stats->total_points, "cumulative number of semicoherent templates" ) == XLAL_SUCCESS, XLAL_EFUNC );
    }
 
    // Write number of computed coherent results, and number of coherent and semicoherent templates
    UINT8 coh_nres = 0, coh_ntmpl = 0, semi_ntmpl = 0;
    XLAL_CHECK_MAIN( XLALWeaveCacheQueriesGetCounts( queries, &coh_nres, &coh_ntmpl, &semi_ntmpl ) == XLAL_SUCCESS, XLAL_EFUNC );
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteUINT8( file, "ncohres", coh_nres, "number of computed coherent results" ) == XLAL_SUCCESS, XLAL_EFUNC );
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteUINT8( file, "ncohtpl", coh_ntmpl, "number of coherent templates" ) == XLAL_SUCCESS, XLAL_EFUNC );
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteUINT8( file, "nsemitpl", semi_ntmpl, "number of semicoherent templates" ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write peak memory usage
    XLAL_CHECK_MAIN( XLALFITSHeaderWriteREAL8( file, "peakmem [MB]", XLALGetPeakHeapUsageMB(), "peak memory usage" ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write timing information
    XLAL_CHECK_MAIN( XLALWeaveSearchTimingWriteInfo( file, tim, queries ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write various information from coherent input data
    XLAL_CHECK_MAIN( XLALWeaveCohInputWriteInfo( file, nsegments, statistics_params->coh_input ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write various information from caches
    XLAL_CHECK_MAIN( XLALWeaveCacheWriteInfo( file, nsegments, coh_cache ) == XLAL_SUCCESS, XLAL_EFUNC );
 
    // Write search results, unless search is being simulated
    if ( simulation_level == 0 ) {
      XLAL_CHECK_MAIN( XLALWeaveOutputResultsWrite( file, out ) == XLAL_SUCCESS, XLAL_EFUNC );
    }
 
    // Write list of injections
    if ( injections != NULL ) {
      XLAL_CHECK_MAIN( XLALFITSWritePulsarParamsVector( file, "injections", injections ) == XLAL_SUCCESS, XLAL_EFUNC );
    }
 
    // Write various segment information from coherent input data
    if ( uvar->segment_info ) {
      XLAL_CHECK_MAIN( XLALWeaveCohInputWriteSegInfo( file, nsegments, statistics_params->coh_input ) == XLAL_SUCCESS, XLAL_EFUNC );
    }
 
    // Close output file
    XLALFITSFileClose( file );
    LogPrintf( LOG_NORMAL, "Closed output file '%s'\n", uvar->output_file );
 
  }
 
  ////////// Cleanup memory and exit //////////
 
  // Cleanup memory from search timing
  XLALWeaveSearchTimingDestroy( tim );
 
  // Cleanup memory from output results
  XLALWeaveOutputResultsDestroy( out );
 
  // Cleanup memory from semicoherent results
  XLALWeaveSemiResultsDestroy( semi_res );
 
  // Cleanup memory from parameter-space iteration
  XLALWeaveSearchIteratorDestroy( main_loop_itr );
 
  // Cleanup memory from computing 'stage 0' coherent results
  XLALWeaveCacheQueriesDestroy( queries );
  for ( size_t i = 0; i < nsegments; ++i ) {
    XLALWeaveCacheDestroy( coh_cache[i] );
  }
 
  // Cleanup memory from loading input data
  XLALDestroySFTCatalog( sft_catalog );
 
  // Cleanup memory from lattice tilings
  for ( size_t i = 0; i < ntiles; ++i ) {
    XLALDestroyLatticeTiling( tiling[i] );
  }
 
  // Cleanup memory from setup data
  XLALWeaveSetupDataClear( &setup );
  XLALFree( setup_detectors_string );
 
  // Cleanup random number generator
  XLALDestroyRandomParams( rand_par );
 
  // Cleanup injection parameters vector
  XLALDestroyPulsarParamsVector( injections );
 
  // Cleanup memory from user input
  XLALDestroyUserVars();
 
  // Check for memory leaks
  LALCheckMemoryLeaks();
 
  if ( search_complete ) {
    LogPrintf( LOG_NORMAL, "Finished successfully!\n" );
  } else {
    LogPrintf( LOG_NORMAL, "Finished but search not completed!\n" );
  }
 
  return EXIT_SUCCESS;
 
}
 
// Local Variables:
// c-file-style: "linux"
// c-basic-offset: 2
// End: