###############
Advanced Topics
###############

==============================
Consolidating fit dictionaries
==============================

To reduce the amount of disk I/O, by default the individual fit dictionaries will be gathered up (along with some useful
meta information) into an overall fit dictionary. If you're working with the full samplers for huge cubes, this file
may be too large to handle on a machine. In which case, you may wish to set ``consolidate_fit_dict = false`` in the
``fitting_parameters`` section of the config file.

==================================
Initial fitting using downsampling
==================================

For complex datasets, doing an initial fit to get a first guess number of components
can speed up the fitting by a significant amount (15% or more). For this, we offer
a "downsampling" fitter, which will rebin the data in the x/y axis by a factor (default 10)
and fit these, before passing those to the full resolution fitting. This can also help to
improve the fits in areas of low signal-to-noise, but generally keeps the fit quality the same.
Note that by doing this, the only thing passed to the fitter is the number of components. It will
still then fit additional components and remove them independently per-spaxel.

.. code-block:: python

  fit_dict_fitting_filename = os.path.join(fit_dir, downsample_fit_dir, fit_dict_filename)
  n_comp_fitting_filename = os.path.join(fit_dir, downsample_fit_dir, n_comp_filename)
  likelihood_fitting_filename = os.path.join(fit_dir, downsample_fit_dir, likelihood_filename)

  hf_fitter.downsample_fitter(fit_dict_filename=fit_dict_fitting_filename,
                              n_comp_filename=n_comp_fitting_filename,
                              likelihood_filename=likelihood_fitting_filename,
                              downsample_factor=10,
                              )

=====================
Different fit methods
=====================

The default behaviour of mcfine is to use MCMC to explore every potential fit after the
initial guesses for parameters. This is robust, but slow. By setting ``fit_method = "leastsq"``
in the ``fit_parameters`` section of the config, you can instead just make use of the least-squares
guesses for the parameters, and then mcfine will only run an MCMC once, after a final number of
components and initial guesses for the parameters for each component has been calculated. This
leads to a speed up of around a factor 4 with a negligible change in the fit quality, so is useful
for large datasets.

===================
Different fit types
===================

We expect that most users will use LTE modelling, there also
exists options for RADEX and pure Gaussian fitting. These
are accessed by the ``fit_type`` parameter:

.. code-block:: toml

  [fitting_params]
  fit_type = "lte"  # Can be "lte", "radex", or "pure_gauss"

===============
Initial guesses
===============

A significant amount of time is spent in producing the initial
guesses via ``lmfit``. This is because by default we use basinhopping,
which works but is slow. The majority of this is because getting initial
guesses for the velocities is critical, but difficult. We also offer
a method called ``iterative``, which finds peaks in the spectrum using
zero-crossings in the derivatives (via specutils). For the total number
of components, it will iteratively subtract off a model for the flux
(from previous fits in this iterative step), find peaks in the
model-subtracted flux (i.e. any unaccounted for emission), and then fit
all the components simultaneously once again (to avoid biases that can arise
by fitting one component at a time) in a small window around the various peaks.
This generally produces identical results in testing to the default method, but
runs significantly (about a factor 5 for a 4-component fit) faster.

If using ``iterative``, the strong recommendation for the minimization algorith
is ``powell``.

======================
``emcee`` optimization
======================

By default, we use a fixed number of walkers and steps in the emcee run. This
may not be optimal, so we also offer a couple of adaptive methods.

The first is you can adapt the number of walkers to the number of parameters
being fitted. A typical rule of thumb is to use **at least** twice the number
of walkers as parameters. This can be specified in the config file by using
a string for ``n_walkers`` under the ``[mcmc]`` section: if you use something like

.. code-block:: toml

   [mcmc]
   n_walkers = "3*n_params"

then the number of walkers will be adapted to the number of components being fit.
This can effectively speed up the fit and reduce sampler filesizes, without affecting
the final fit quality

The second is using an adaptive number of steps, checking against an autocorrelation criterion.
By default, we use a fixed number of steps, using a quarter of these as "burn-in", before a full
run. Then, when calculating parameters, we will discard the first half of this chain as parameters
may still be moving around. However, it may be more optimal to use an adaptive number of steps,
and tune this to some multiple of the autocorrelation length. More details of this can be found
in `the emcee docs <https://emcee.readthedocs.io/en/stable/tutorials/autocorr/>`_. You
can switch to this mode like so:

.. code-block:: toml

   [mcmc]
   emcee_run_method = "adaptive"

this will use another few parameters:

.. code-block:: toml

   [mcmc]
   convergence_factor = 10
   tau_change = 0.01
   thin = 0.5
   burn_in = 2
   max_steps = 100000

The first two control when fitting will stop. We will step out of the MCMC run once the chain
is ``convergence_factor`` times the autocorrelation length, and the change in that factor varies
by less than a factor of ``tau_change``. To ensure the MCMC doesn't run forever, the maximum number
of steps is ``max_steps``.

After this, for both methods we then use ``thin`` and ``burn_in`` to control how the sampler is pared down for
parameter estimation. We use a factor of ``burn_in`` time the autocorrelation length to define the burn-in, and then
thin by a factor of ``thin`` times the autocorrelation length to thin out the chains. For more details on
these, see the `emcee docs <https://emcee.readthedocs.io/en/stable/tutorials/monitor/>`_.

============
Space saving
============

In some cases it may be worthwhile to try and save as much space as possible
for the output fit dictionaries. This is especially true for huge datasets
where the total space requirements can ramp up into the TB. By default,
the ``fit_dict`` files will contain both the full ``emcee`` sampler and
information for a covariance matrix. If your corner plots look well-behaved
(i.e. mostly Gaussian), then you can save a **significant** amount of space
by setting ``keep_sampler = false`` in ``[fitting_params]``.

=======================
Configurable parameters
=======================

To see the configurable parameters in McFine, there is an
in-built convenience function:

.. code-block:: python

    from mcfine.utils import print_config_params
    print_config_params()

This will list all the parameters, as well as their type
and default values. You can also see the default values
McFine will use in ``mcfine/toml/`` in the GitHub repository,
which may also be useful for those who aren't too familiar
with toml.

==================
Exploring samplers
==================

Although we expose convenience functions for exploring
the fits (see :doc:`here <tutorials/exploring_cube_fits>`), you
can also directly access the `emcee` sampler object:

.. code-block:: python

    from mcfine.utils import load_fit_dict

    fit_dict = load_fit_dict(file_name)
    sampler = fit_dict["sampler"]

from there, you can mess around with this as you'd like.

=============================
Exploring covariance matrices
=============================

If you set ``keep_covariance`` to True in the config, then
along with parameter dictionaries a covariance dictionary
will also be saved out. You can load this in just like
a fit dictionary, and for each pixel fit (in i, j notation),
there will be a median array and covariance matrix. From this,
you can generate samples like:

.. code-block:: python

   import numpy as np

   samples = np.random.multivariate_normal(
                    par_dict[f"{i}_{j}"]["cov_med"],
                    par_dict[f"{i}_{j}"]["cov_matrix"],
                    size=10000,
             )

===================
Adding another line
===================

It is possible to add other lines to McFine relatively
simply. The majority of the info just needs to be put
into ``line_info.py``. For the LTE case, these are
``v_lines`` and ``strength_lines``. For a single-peak
line, this is just 0 and 1. For RT there's
``transition_lines`` and ``freq_lines``. These should
descend from the RADEX naming scheme. Once you've added
those, include your new line in ``ALLOWED_LINES`` in
``fitting.py``, and edit the config file to use this
new line.

============================
Limiting number of processes
============================

McFine is highly multi-processed, but so are a number of
packages that McFine relies on which can cause issues,
especially on larger machines. To limit the number of threads
packages such as numpy will use, before you call your code you can
put (in the shell):

.. code-block:: shell

    setenv MKL_NUM_THREADS 1
    setenv NUMEXPR_NUM_THREADS 1
    setenv OMP_NUM_THREADS 1
    setenv OPENBLAS_NUM_THREADS 1
    setenv GOTO_NUM_THREADS 1

or the equivalent EXPORT call, depending on the shell you use.

For macOS >= 10.13, you may also need to set

.. code-block:: shell

    OBJC_DISABLE_INITIALIZE_FORK_SAFETY=YES