Merge pull request #26 from Zstone19/mica

Add MICA2 module

README.md

@@ -10,6 +10,8 @@
 [![python-ver](https://img.shields.io/pypi/pyversions/pypetal)](https://img.shields.io/pypi/pyversions/pypetal)
 [![pypi-downloads](https://static.pepy.tech/badge/pypetal)](https://pepy.tech/project/pypetal)
 
+<a href="https://ascl.net/2401.004"><img src="https://img.shields.io/badge/ascl-2401.004-blue.svg?colorB=262255" alt="ascl:2401.004" /></a>
+
 
 pyPETaL is a time-series data analysis pipeline for AGN reverberation mapping (RM) data. It combines multiple different popular softwares using for AGN RM analysis, including PyCCF, PyZDCF, JAVELIN, and PyROA. This pipeline also implements outlier rejection using Damped Random Walk Gaussian proces fitting, and detrending through the LinMix algorithm. pyPetal implements a weighting scheme (for all modules) in order to mitigate aliasing in peaks of time lag distributions between light curves.
 
@@ -40,6 +42,9 @@ PLIKE is an optional algorithm that is used in pyPetal. There is a script availa
 ```
 
 
+
+### Detrending
+
 pyPetal offers the option to detrend the input light curves via the [LinMix](https://github.com/jmeyers314/linmix.git) algorithm. This package is not offered in the base version of pyPetal, but can be installed with pip:
 ```
     pip install "linmix @ git+https://github.com/jmeyers314/linmix.git"
@@ -48,8 +53,22 @@ pyPetal offers the option to detrend the input light curves via the [LinMix](htt
 Or with ``poetry``:
 ```
     poetry install --with extra
+```
+
 
+### MICA2
+
+pyPetal also offers [MICA2](https://github.com/LiyrAstroph/MICA2) as an optional module. This package is more complex to install than the others - to find out more, read the [MICA2 installation guide](https://mica2.readthedocs.io/en/latest/getting_started.html).
+
+``poetry`` can help to install some easily-installable python dependencies for MICA2:
 ```
+    poetry install --with extra
+```
+
+The functionality and inputs of the MICA2 module are identical to the original software, so the best way to learn how this module functions is to read the MICA2 documentation.
+
+__NOTE__: Seeing as this is an optional module, pyPetal will still compile and run if MICA2's dependencies are not installed. Each pyPetal run will assume ``run_mica2=False``.
+
 
 
 ### pyPetal and JAVELIN
@@ -82,3 +101,14 @@ __NOTE:__ The user may need to install ``NumPy`` before installing pyPetal-jav t
 ```
     pip install numpy
 ```
+
+
+
+## Citing pyPetal
+
+To cite the pyPetal code itself, use the [ASCL reference](https://ascl.net/2401.004):
+```
+    Stone Z., 2024, version xxxx, Astrophysics Source Code Library, record ascl:2401.004
+```
+
+Cit the paper pyPetal was used in: [Shen et al. (2023)](https://ui.adsabs.harvard.edu/abs/2023arXiv230501014S/abstract) 

docs/index.rst

@@ -8,7 +8,7 @@ pyPETaL: A Pipeline for Estimating AGN Time Lags
 =================================================
 **pyPetal** is a pipeline made in python to obtain time lags from multi-band AGN time-series data. Normally, AGN geometric and kinematic analyses (e.g. reverberation mapping) utilize a variety of different tools to obtain time lags between two (or more) light curves. This package combines four popular algorithms for estimating time lags (pyCCF, pyZDCF, PyROA, and JAVELIN), and uses the popular Damped Random Walk algorithm to model input light curves for outlier rejection.
 
-Currently, pyPetal has combined the functionality of pyCCF, pyZDCF, PyROA, and JAVELIN to produce cross-correlation functions, discrete correlation functions, and mean time lags. This is only made to run on Linux-based operating systems, though this may be improved in the future.
+Currently, pyPetal has combined the functionality of pyCCF, pyZDCF, PyROA, MICA2m and JAVELIN to produce cross-correlation functions, discrete correlation functions, and mean time lags. This is only made to run on Linux-based operating systems, though this may be improved in the future.
 
 
 .. toctree::
@@ -29,6 +29,7 @@ Currently, pyPetal has combined the functionality of pyCCF, pyZDCF, PyROA, and J
    notebooks/pyzdcf
    notebooks/plike
    pyroa_toc
+   mica2_toc
    javelin_toc
    notebooks/weighting
    notebooks/all_together
@@ -57,6 +58,17 @@ Currently, pyPetal has combined the functionality of pyCCF, pyZDCF, PyROA, and J
 
 
 
+
+Citing pyPetal
+--------------
+
+To cite the pyPetal code itself, use the ASCL reference: `</https://ascl.net/2401.004>`__
+
+To cite the paper pyPetal was used in: `<https://ui.adsabs.harvard.edu/abs/2023arXiv230501014S>`__
+
+
+
+
 References
 ----------
 
@@ -66,4 +78,5 @@ pyPetal makes use of a multitude of packages:
 * **pyZDCF**: `<https://doi.org/10.5281/zenodo.7253034>`__,
 * **PLIKE**: `<http://arxiv.org/abs/1302.1508>`__
 * **PyROA**: `<https://ui.adsabs.harvard.edu/abs/2021MNRAS.508.5449D>`__
+* **MICA2**: `<https://ui.adsabs.harvard.edu/abs/2016ApJ...831..206L>`__
 * **JAVELIN**: `<http://adsabs.harvard.edu/abs/2013ApJ...765..106Z>`__, `<http://adsabs.harvard.edu/abs/2011ApJ...735...80Z>`__, `<http://adsabs.harvard.edu/abs/2016ApJ...819..122Z>`__

docs/installation.rst

@@ -26,6 +26,21 @@ PLIKE is an optional algorithm that is used in pyPetal. There is a script availa
 
 
 
+MICA2
+------
+
+pyPetal also offers `MICA2 <https://github.com/LiyrAstroph/MICA2>`_ as an optional module. This package is more complex to install than the others - to find out more, read the `MICA2 installation guide <https://mica2.readthedocs.io/en/latest/getting_started.html>`_.
+
+``poetry`` can help to install some easily-installable python dependencies for MICA2:
+::
+    poetry install --with extra
+
+
+The functionality and inputs of the MICA2 module are identical to the original software, so the best way to learn how this module functions is to read the MICA2 documentation.
+
+.. note:: Seeing as this is an optional module, pyPetal will still compile and run if MICA2's dependencies are not installed. Each pyPetal run will assume ``run_mica2=False``.
+
+
 
 
 pyPetal and JAVELIN

docs/mica2_toc.rst

@@ -0,0 +1,34 @@
+=======
+MICA2
+=======
+
+MICA2 is a more complex tool used to correlate two light curves in a non-parametric fashion to find the transfer function between the two, 
+as a sum of either Gaussians or tophat functions. MICA2 is much more flexible in its functionality, though this comes with much more complex
+customization from the user in terms of parameters chosen and types of analysis possible. These tutorials provide instructions
+on in-depth usage of MICA2 and understanding of its arguments and outputs.
+
+While these tutorials do a sufficient job of introducing the user to the basic functionality of MICA2, the user is encouraged to read the
+MICA2 `documentation <https://mica2.readthedocs.io/en/latest/index.html>`_ and/or `repository <https://github.com/LiyrAstroph/MICA2>`_ for 
+more in-depth information on the software. In addition, the user can read the associated `paper <https://ui.adsabs.harvard.edu/abs/2016ApJ...831..206L/abstract>`_
+for any additional details.
+
+.. warning::
+    MICA2 is a software that can normally be run on multiple cores by using ``mpiexec``, such as with::
+
+        mpiexec -n 4 python mica2_script.py
+
+
+    pyPetal extends this functionality, so that the MICA2 module can utilize ``mpiexec``. However, if this is done, all other module will
+    only be run on one core (as they are optimized for python multi-threading instead of MPI). If the user wishes to only use the MICA2 module,
+    it may be best to run using ``mpiexec``. If the user wishes to use the other modules, it may be best to run the MICA2 module separately 
+    from the other modules.
+
+
+.. toctree::
+    :maxdepth: 1
+    :caption: MICA2 Tutorials
+
+    notebooks/mica2/basic
+    notebooks/mica2/numcomp
+    notebooks/mica2/noorder
+    notebooks/mica2/together

docs/pl_args.rst

@@ -288,6 +288,128 @@ Module: pyROA (``run_pyroa``)
       - :python:`None`
 
 
+
+Module: MICA2 (``run_mica2``)
+-------------------------------
+
+MICA2 arguments:
+
+.. list-table::
+    :widths: 16 50 16 16
+    :header-rows: 1
+
+    * - Argument
+      - Description
+      - Type
+      - Default
+    * - ``type_tf``
+      - The type of transfer function to use in MICA2. Can be either "gaussian" or "tophat".
+      - :python:`str`
+      - :python:`"gaussian"`
+    * - ``max_num_saves``
+      - The number of saves to use in CDNest for the sampling.
+      - :python:`int`
+      - 2000
+    * - ``flag_uniform_var_params``
+      - Whether each dataset has the same variability parameters.
+      - :python:`bool`
+      - :python:`False`
+    * - ``flag_uniform_transfuns``
+      - Whether each dataset has the same line parameters.
+      - :python:`bool`
+      - :python:`False`
+    * - ``flag_trend``
+      - Whether or not to include a trend in the transfer function. 0 for constant trend, 1 for linear trend, 2 for conic trend.
+      - :python:`int`
+      - 0
+    * - ``flag_lag_posivity``
+      - Whether or not to force Gaussians to be located at positive lags.
+      - :python:`bool`
+      - :python:`False`
+    * - ``flag_negative_resp``
+      - Whether or not to turn on negative response.
+      - :python:`bool`
+      - :python:`False`
+    * - ``number_component``
+      - The lower and upper limits to use for the number of Gaussians. If only one value is given, it will be assumed for both limits.
+      - list of :python:`int`
+      - [1, 1]
+    * - ``width_limit``
+      - The lower and upper limits for the width of the Gaussians. If only one value is given, it will be assumed for both limits. If :python:`None`, these will be determined by MICA2.
+      - list of :python:`float`
+      - :python:`None`
+    * - ``flag_con_sys_err``
+      - Whether or not to include a constant systematic error for the continuum light curve.
+      - :python:`bool`
+      - :python:`False`
+    * - ``flag_line_sys_err``
+      - Whether or not to include a systematic error for the line light curve(s).
+      - :python:`bool`
+      - :python:`False`
+    * - ``type_lag_prior``
+      - The type of prior to use for the lags. See the MICA2 documentation.
+      - :python:`int`
+      - 0
+    * - ``lag_prior``
+      - The parameters to use for the lag prior. See the MICA2 documentation.
+      - list of :python:`float`
+      - :python:`None`
+    * - ``together``
+      - Whether or not to fit the time lags of all light curves together, with the same transfer function.
+      - :python:`bool`
+      - :python:`False`
+    * - ``no_order``
+      - Whether or not the order of the time lags makes a difference (i.e., whether or not the time lags are allowed to be equal). If order doesn't matter, MICA2 will be able to fit all light curves in one run. If set to :python:`True`, a separate MICA2 run will be performed for each line.
+      - :python:`bool`
+      - :python:`True`
+
+
+
+
+CDNest arguments:
+
+.. list-table::
+    :widths: 16 50 16 16
+    :header-rows: 1
+
+    * - Argument
+      - Description
+      - Type
+      - Default
+    * - ``num_particles``
+      - See the MICA2 or CDNest documentation.
+      - :python:`int`
+      - 1
+    * - ``thread_steps_factor``
+      -
+      - :python:`int`
+      - 1
+    * - ``new_level_interval_factor``
+      -
+      - :python:`int`
+      - 1
+    * - ``save_interval_factor``
+      -
+      - :python:`int`
+      - 1
+    * - ``lam``
+      -
+      - :python:`float`
+      - 10
+    * - ``beta``
+      -
+      - :python:`float`
+      - 100
+    * - ``ptol``
+      -
+      - :python:`float`
+      - 0.1
+    * - ``max_num_levels``
+      -
+      - :python:`int`
+      - 0
+
+
 Module: JAVELIN (``pypetal_jav.run_pipeline``)
 ----------------------------------------------
 

docs/pl_files.rst

@@ -12,13 +12,15 @@ In general, running all of the modules with three lines (named "cont", "line1",
     │   └── detrend.pdf
     ├── line1/
     │   ├── drw_rej/
+    │   ├── mica2/
     │   ├── pyccf/
     │   ├── pyzdcf/
     │   ├── javelin/
     │   ├── weights/
     │   └── detrend.pdf
     ├── line2/
     │   ├── drw_rej/
+    │   ├── mica2/
     │   ├── pyccf/
     │   ├── pyzdcf/
     │   ├── javelin/
@@ -28,6 +30,7 @@ In general, running all of the modules with three lines (named "cont", "line1",
     ├── pyroa/
     ├── pyroa_lcs/
     ├── light_curves/
+    ├── mica2_weights_res.pdf
     ├── pyccf_weights_res.pdf
     ├── pyroa_weights_res.pdf
     └── javelin_weights_res.pdf
@@ -262,6 +265,69 @@ If :python:`together=False`, the columns will be the same as for :python:`togeth
 
 
 
+Module: MICA2
+-------------
+
+Like the PyROA module, the output of this module depends on the ``together`` parameter.
+
+If :python:`together=True`, the output directory for all lines will be ``output_directory/mica2/``. If :python:`together=False`, each line will have it's MICA2 results in its own subdirectory, labeled ``mica2/``.
+
+Each PyROA directory (whether ``together`` is :python:`True` or :python:`False`) will have both a ``data/`` and ``param/`` directory, which were used by MICA2 to store the CDNest sampling and output information. To learn more about this data, see the MICA2 and CDNest documentation.
+
+In general, the names of the files will depend on the number of gaussians/tophats used in the analysis. There will be a file for every gaussian used, indicated by a number indexed at 1.
+
+Only a few files will be of note in the ``data/`` directory, which are two figures:
+
+.. list-table::
+    :widths: 50 50
+    :header-rows: 1
+
+    * - Filename
+      - Description
+    * - ``cdnest_{ngauss}.pdf``
+      - A figure showing the post-processing analysis of the diffusive nested sampling process.
+    * - ``fig_{ngauss}.pdf``
+      - A figure showing the quality of the MICA2 fits, including the center/centroid histograms, the transfer function, and the fits to the light curves.
+
+
+Additionally, pyPetal will save the following files in the ``mica2/`` directory:
+
+.. list-table::
+    :widths: 30 30 10 30
+    :header-rows: 1
+
+    * - Filename
+      - Description
+      - Format
+      - Columns
+    * - ``cont_recon.dat``
+      - The reconstructed continuum light curve.
+      - CSV
+      - time, value, uncertainty
+    * - ``{line_name}_recon.dat``
+      - The reconstructed line light curve.
+      - CSV
+      - time, value, uncertainty
+    * - ``{line_name}_centers_{ngauss}.dat``
+      - The output samples for the centers of the gaussians for a given line and gaussian/tophat.
+      - CSV
+      - value
+    * - ``{line_name}_centroids_{ngauss}.dat``
+      - The output samples for the centroids of the gaussians for a given line and gaussian/tophat.
+      - CSV
+      - value
+    * - ``{line_name}_transfunc.dat``
+      - The transfer function for a given line and gaussian/tophat.
+      - CSV
+      - tau, transfer_function, lower_uncertainty, upper_uncertainty
+
+If :python:`together=True`, the only difference will be that the transfer function file will be named ``transfunc.dat``.
+
+If :python:`together=False` and :python:`no_order=False`, the ``data/`` and ``param/`` directories will be located in ``output_directory/mica2/`` and the individual sample files will be located in ``output_directory/{line_name}/mica2/``.
+
+
+
+
 Module: JAVELIN
 ---------------
 

docs/pl_output.rst

@@ -212,6 +212,12 @@ The ``MyFit`` object will have the following attributes:
       - list of :python:`float`
 
 
+Module: MICA2
+--------------
+
+The MICA2 module's output will be a list of ``pymica.gmodel`` objects output from MICA2 itself, in the order of the line names input into pyPetal.
+
+
 
 Module: JAVELIN
 ---------------