Migrate gcorr iteration script to library entry point (#33)

Author: arahlin

MANIFEST.in

@@ -1 +1 @@
-recursive-include docs example scripts
+recursive-include docs example

docs/api.rst

@@ -35,9 +35,17 @@ xfaster_exec
   submits the job to a grid system
 
 spec_tools
--------------
+----------
 .. automodule:: xfaster.spec_tools
   :members:
   :undoc-members:
 
-  More generic functions used by methods in xfaster_class
+  More generic functions used by methods in :py:mod:`~xfaster.xfaster_class`
+
+gcorr_tools
+-----------
+.. automodule:: xfaster.gcorr_tools
+  :members:
+  :undoc-members:
+
+  Tools for iteratively computing corrections to the ``g_ell`` factor

docs/gcorr.rst

@@ -0,0 +1,79 @@
+G-Correction Calculation
+========================
+
+Code package for calculating correction factors to g_ell using pre-generated
+simulation maps for each single mask. This procedure assumes you are running
+this code on a cluster-- it is extremely slow to do otherwise.
+
+There is one main command-line utility that calls two functions from the
+:py:mod:`~xfaster.gcorr_tools` library:
+
+1. :py:func:`~xfaster.gcorr_tools.xfaster_gcorr_once` -- This function calls
+   XFaster to run or submit jobs for gcorr runs.
+2. :py:func:`~xfaster.gcorr_tools.process_gcorr` -- A function that computes the
+   gcorr correction from the ensemble of bandpowers, updates the running total,
+   and backs up the necessary files from each iteration.
+3. ``xfaster gcorr`` -- An xfaster command-line utility, that calls function 1
+   and 2 to get a new gcorr.npz file each time. This is the main code you'll
+   run.
+
+There is also a config file with options specific to computing gcorr.
+
+Procedure
+---------
+
+1. Edit the gcorr config file to suit your purposes. Examples are given for
+   signal and null runs in the examples directory. Required fields are:
+
+   * ``null`` -- must be true for null tests and false for signal runs
+   * ``map_tags`` -- comma-separated list of map tags
+   * ``data_subset`` -- the globabble data_subset argument to xfaster_run, but
+     without map tags. So, "full", "chunk*", etc.
+   * ``output_root`` -- the parent directory where your gcorr XFaster runs will
+     be written
+   * ``nsim`` -- the number of simulations to use to compute gcorr
+   * ``[xfaster_opts]`` -- this is where you'll put any options that will be
+     directly input to :py:func:`~xfaster.xfaster_exec.xfaster_run`
+   * ``[submit_opts]`` -- this is where you'll put any options that will be
+     directly input to :py:func:`~xfaster.xfaster_exec.xfaster_submit`, in
+     addition to those in ``[xfaster_opts]``
+
+2. Run ``xfaster gcorr`` once to get the full set of XFaster output files in the
+   output directory.  Since we haven't computed gcorr yet, this will set
+   ``apply_gcorr=False``. Make sure to use as many OMP threads as possible since
+   this is the step where the sims_xcorr file, which benefits the most from
+   extra threads, is computed.  Your command should look like this::
+
+      xfaster gcorr path-to-my-gcorr-config.ini
+
+3. Run ``xfaster gcorr`` until convergence is reached. In practice, you will run
+   the command above and wait for it to finish. If you include the
+   ``--max-iters`` option with a non-zero value, the code will try to determine
+   whether convergence or max_iters has been reached and stop on its own.
+   Otherwise, you can look at the correction-to-the-correction that it both
+   prints and plots (it should converge to 1s for TT, EE, BB), and rerun the
+   same command if it hasn't converged.  In much more detail, here's what the
+   code does:
+
+   1. Call :py:func:`~xfaster.gcorr_tools.xfaster_gcorr_once` for the 0th sim
+      seed while also reloading gcorr (if this is not the first iteration).
+      This does a couple things-- saves the new gcorr in the ``masks_xcorr``
+      file, so later seeds will use the right thing. And recompute the transfer
+      function, which doesn't depend on the ``sim_index``, so is only necessary
+      to do once.
+   2. After the transfer functions are all on disk, submit individual jobs for
+      all the other seeds, just doing the bandpowers step for those.
+   3. Once they're all done, run :py:func:`~xfaster.gcorr_tools.compute_gcal`,
+      and save a correction-to-gcorr as ``gcorr_corr_<tag>_iter<iter>.npz`` in the
+      rundir.
+   4. If not the first iteration, load up the correction-to-gcorr computed for
+      this iteration. Multiply it by the total gcorr, and save that to the
+      output directory as ``gcorr_total_<tag>_iter<iter>.npz``.
+   5. Plot gcorr total and the correction to gcorr total. Save in rundir/plots.
+   6. Clear out rundir bandpowers/transfer functions/logs.
+   7. Exit.
+
+4. After convergence is reached, copy the gcorr_total file for the last
+   iteration from the rundir to the mask directory, labeling it
+   ``mask_map_<tag>_gcorr.npz`` for signal or ``mask_map_<tag>_gcorr_null.npz``
+   for null.

docs/index.rst

@@ -20,6 +20,7 @@ More details of the math can be found in :ref:`Algorithm<Algorithm>`, and a deta
   quickstart
   algorithm
   notebooks/XFaster_Tutorial
+  gcorr
   api
 
 Indices and tables

docs/notebooks/XFaster_Tutorial.ipynb

@@ -420,7 +420,7 @@
    "source": [
     "We have now finished all the set-up. Now it's time to start calculating things. First, we compute the cross spectra of our masks--which will be needed for computing the $K_{\\ell \\ell^{\\prime}}$ mode-coupling matrix--and the $g_\\ell$ mode-counting factor. This is done with [get_mask_weights()](../api.rst#xfaster.xfaster_class.XFaster.get_mask_weights).\n",
     "\n",
-    "For the example, we will not apply an empirical correction to $g_\\ell$, the calibration of which is discussed in Section 2.3.2 of the [XFaster paper](https://arxiv.org/abs/2104.01172).  Otherwise, we would set `apply_gcorr=True`, and the code would look in the masks directory for this correction file for each map. The `reload_gcorr` option is only useful when doing the empirical calibration; it reloads the file, while by-passing the checkpoint tree that is usually performed after `get_mask_weights`."
+    "For the example, we will not apply an empirical correction to $g_\\ell$, the calibration of which is discussed in Section 2.3.2 of the [XFaster paper](https://arxiv.org/abs/2104.01172).  Otherwise, we would set `apply_gcorr=True`, and the code would look in the masks directory for this correction file for each map. The `reload_gcorr` option is only useful when doing the empirical calibration; it reloads the file, while by-passing the checkpoint tree that is usually performed after `get_mask_weights`.  Instructions are [provided](../gcorr.rst) for constructing this calibration using the example dataset."
    ]
   },
   {

scripts/gcorr/gcorr_config_null.ini example/gcorr_config_null.iniscripts/gcorr/gcorr_config_null.ini renamed to example/gcorr_config_null.ini

@@ -3,16 +3,16 @@
 null = true
 map_tags = 95,150
 data_subset = full
-output_root = ../../example/null/gcorr_run
+output_root = null/gcorr_run
 nsim = 100
 
 # Options we can directly pass to XFaster
 [xfaster_opts]
-config = ../../example/config_example.ini
+config = config_example.ini
 # Make null map set using the same set as for the signal run.
-data_root = ../../example/maps_example_half1
+data_root = maps_example_half1
 # Need second data root for null run
-data_root2 = ../../example/maps_example_half2
+data_root2 = maps_example_half2
 signal_type = synfast
 # Noise type needed for nulls
 noise_type = gaussian

scripts/gcorr/gcorr_config_signal.ini example/gcorr_config_signal.iniscripts/gcorr/gcorr_config_signal.ini renamed to example/gcorr_config_signal.ini

@@ -3,13 +3,13 @@
 null = false
 map_tags = 95,150
 data_subset = full
-output_root = ../../example/gcorr_run
+output_root = gcorr_run
 nsim = 100
 
 # Options we can directly pass to XFaster
 [xfaster_opts]
-config = ../../example/config_example.ini
-data_root = ../../example/maps_example
+config = config_example.ini
+data_root = maps_example
 # data_root2 = # change for nulls
 signal_type = synfast
 # noise type ignored for signal gcorr