README_long.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7b3b2bb6-42f4-4db3-a8cb-8ba256c872f9",
   "metadata": {},
   "source": [
    "# TRADES\n",
    "  \n",
    "**`TRADES` v2.19.0 by Luca Borsato - 2016-2023**  \n",
    "\n",
    "Most of the information can be found in the papers:  \n",
    "\n",
    "- [Borsato et al. (2014)][Borsato2014]  \n",
    "- [Malavolta et al. (2017)][Malavolta2017]  \n",
    "- [Borsato et al. (2019)][Borsato2019]  \n",
    "\n",
    "and at the webpage [TRADES@ESPG][TRADESESPG].  \n",
    "Feel free to use or modify the code, but please cite the three above papers.  \n",
    "Comments are welcome!  \n",
    "\n",
    "**WARNING:**  \n",
    "**only tested on a Unix/Linux machine (i.e., Centos Rocks 5.3, Ubuntu > 12.04 and derivatives)**  \n",
    "**Please use `gfortran` version greater than `4.8.1`. Previous versions could fail to compile.**  \n",
    "**`openMP` required to compile parallel versions of `TRADES`**  \n",
    "**Required `python` packages:**  \n",
    "\n",
    "- `numpy` (`conda`),  \n",
    "- `matplotlib` (`conda`),  \n",
    "- `h5py` (`conda`),  \n",
    "- `yaml` (`conda`),  \n",
    "- `scipy` (`conda`),  \n",
    "- `pandas` (`conda`),  \n",
    "- `scikit-learn` (`conda-forge`),  \n",
    "- `statsmodels` (`conda-forge`),  \n",
    "- [`emcee`][emcee] (`conda-forge`) with `tqdm` (`conda-forge`),  \n",
    "- `pygtc` (`conda-forge`),  \n",
    "<!-- - `corner`,   -->\n",
    "- `PyDE` (it is preferred to install `pytransit` with `pip` or `github` that store an updated version of `PyDE`; be aware that `pytransit` needs `astropy` and `arviz` and sometimes you have to re-install the `semantic_version` package...),  \n",
    "- `ultranest` only for a testing script (currently not working).  \n",
    "\n",
    "---\n",
    "\n",
    "## TRADES user guide\n",
    "\n",
    "1. [Introduction](#intro)  \n",
    "2. [Install and Compile](#install_compile)  \n",
    "3. [How to run TRADES](#run_trades)  \n",
    "4. [Files needed by TRADES](#in_files)  \n",
    "5. [Output files by TRADES](#out_files)  \n",
    "6. [Python Library](#pylib)  \n",
    "  \n",
    "**For the impatient search for: COMMON WAY TO COMPILE**  \n",
    "\n",
    "### Introduction {#intro}\n",
    "\n",
    "We have developed a computer program (in Fortran 90, openMP, and MPI) for determining \n",
    "the possible physical and dynamical configurations of extra-solar planetary \n",
    "systems from observational data, known as `TRADES`, which stands for TRAnsits and\n",
    "Dynamics of Exoplanetary Systems.  \n",
    "The program `TRADES` models the dynamics of multiple planet systems and\n",
    "reproduces the observed transit times ($T_{0}$, or mid-transit times),\n",
    "full transit durations ($T_{14}$ or $T_{41}$),\n",
    " and radial velocities (RVs).\n",
    "These $T_{0}\\textrm{s}$ and RVs are computed during the integration of the planetary orbits.  \n",
    "<!-- Hereafter, I will use the reduced chi-squared ($\\chi^{2}_\\textrm{r} = \\chi^{2}/\\textrm{dof}$)\n",
    "as the fitness parameter. -->\n",
    "<!-- , but it would be also possible to use a $\\chi^{2}_\\textrm{w,dof}$ weighted/scaled -->\n",
    "<!-- by the number of each type of data set: transits and radial velocities.   -->\n",
    "We have developed `TRADES` from scratch because we want to avoid black-box programs,\n",
    "it would be easier to parallelise it with openMP, and include additional algorithms.  \n",
    "To solve the inverse problem, `TRADES` can be run in different modes:  \n",
    "\n",
    "- **integration:**  \n",
    "  it runs a simple integration of the orbits of the planetary system calculating the $T_0$s and the RVs.  \n",
    "\n",
    "- **grid search:**  \n",
    "  TRADES samples the orbital elements of one perturbing body \n",
    "  (all the parameters are allowed in the grid, but remember that the number of simulations\n",
    "  will increase hugely).\n",
    "  The grid parameters can be evenly sampled on a fixed grid by setting\n",
    "  the number of steps, or the step size, or by a number of points\n",
    "  chosen randomly within the parameter bounds.\n",
    "  For any given set of values, the orbits are integrated, and the residuals\n",
    "  between the observed and computed $T_0$s and RVs are computed.\n",
    "  For each combination of the parameters, the `LM` algorithm can be called and\n",
    "  the best case is the one with the lowest residuals\n",
    "  (lowest reduced chi-squared = $\\chi^{2}/\\textrm{dof} = \\chi^{2}_\\textrm{r}$).  \n",
    "\n",
    "- **Levenberg-Marquardt (`LM`, `lmdif` from [`MINPACK`][MINPACK]) algorithm:**  \n",
    "  After an initial guess on the orbital parameters of the perturber,\n",
    "  which could be provided by the previously described grid approach,\n",
    "  the LM algorithm exploits the Levenberg-Marquardt minimization method\n",
    "  to find the solution with the lowest residuals.\n",
    "  The LM algorithm requires the analytic derivative of the model with respect to\n",
    "  the parameters to be fitted.\n",
    "  Since the $T_0$s are determined by an iterative method and the\n",
    "  radial velocities are computed using the numerical integrator,\n",
    "  we cannot express these as analytic functions of fitting parameters.\n",
    "  We have adopted the method described in [Moré et al. (1980)][Moré1980]\n",
    "  to compute the Jacobian matrix, which is determined by a forward-difference\n",
    "  approximation.\n",
    "  The `epsfcn` parameter, which is the parameter that determines the first\n",
    "  Jacobian matrix, is automatically selected in a logarithmic range from\n",
    "  the machine precision up to $10^{−6}$ ; the best value is the one\n",
    "  that returns the lower $\\chi^2_\\textrm{r}$.\n",
    "  This method has the advantage to be scale invariant, but it assumes that\n",
    "  each parameter is varied by the same `epsfcn` value (e.g., a variation of\n",
    "  10% of the period has a different effect than a variation of the same\n",
    "  percentage of the argument of pericenter).  \n",
    "  \n",
    "- **genetic algorithm (`GA` or `PIK`, we used the implementation named [`PIKAIA`][PIKAIA],\n",
    "  Charbonneau 1995):**\n",
    "  the `GA` mode searches for the best orbit by performing a genetic optimization\n",
    "  (e.g. Holland 1975; Goldberg 1989), where the fitness parameter was\n",
    "  the inverse of the $\\chi^2_\\textrm{r}$, now it is set to the log-likelihood ($\\log{\\mathcal{L}}$).\n",
    "  This algorithm is inspired by natural selection which is the biological process\n",
    "  of evolution.\n",
    "  Each generation is a new population of *offspring* orbital configurations,\n",
    "  that are the result of *parent* pairs of orbital configurations that are\n",
    "  ranked following the fitness parameter. \n",
    "  A drawback of the `GA` is the slowness of the algorithm, when compared to other\n",
    "  optimizers.\n",
    "  However, the `GA` should converge to a global solution (if it exists)\n",
    "  after the appropriate number of iterations.  \n",
    "  <!-- At the moment if the `GA` after the predefine iteration has a fitness\n",
    "  (or $\\chi^2_\\textrm{r}$) higher than 1000,\n",
    "  it will continue (it stops if it reaches the iteration step with fitness lower than 1000 or\n",
    "   if it reaches fitness <= 1.).   -->\n",
    "  \n",
    "- **particle swarm optimization ([`PSO`][PSO], Tada 2007):**  \n",
    "  the `PSO` is another optimization algorithm that searches for the global solution\n",
    "  of the problem; this approach is inspired by the social behavior of\n",
    "  bird flock and fish school (e.g., Kennedy & Eberhart 1995; Eberhart 2007).\n",
    "  The fitness parameter used is the same as the `GA`,\n",
    "  <!-- the inverse of the `$\\chi^2_\\textrm{r}$`. -->\n",
    "  the log-likelihood ($\\log{\\mathcal{L}}$).\n",
    "  For each *particle*, the next step (or iteration) in the space of the\n",
    "  fitted parameters is mainly given by the combination of three terms:\n",
    "  random walk, best *particle* position (combination of parameters),\n",
    "  and best *global* position (best orbital configuration of the all particles\n",
    "  and all iterations).  \n",
    "\n",
    "- **emcee ([`emcee`][emcee]):**  \n",
    "  see [`emcee`][emcee] docs and papers.  \n",
    "\n",
    "- **PyDE ([`PyDE`][PyDE]):**  \n",
    "  see [`PyDE`][PyDE] docs and papers.  \n",
    "  \n",
    "In each mode, `TRADES` compares observed transit times ($T_{0,\\textrm{obs}}$) and \n",
    "radial velocities ($\\textrm{RV}_{\\textrm{obs}}$) with the simulated ones ($T_{0,\\textrm{sim}}$ and $\\textrm{RV}_\\textrm{sim}$).\n",
    "From version 1.1.2 of `TRADES`, it is possible to use different set of RV, with\n",
    "different RV offset (the so-called _gamma_ point);\n",
    "`TRADES` will compute a $\\gamma$ for each RV data set.  \n",
    "The *grid* search is a good approach in case that we want to explore a limited\n",
    "subset of the parameter space or if we want to analyze the behavior of the system\n",
    "by varying some parameters, for example to test the effects of a growing mass\n",
    "for the perturbing planet.  \n",
    "`GA` and `PSO` are good methods to be used in case of a wider space of parameters.\n",
    "The orbital solution determined with the `GA` or the `PSO` method is eventually\n",
    "refined with the `LM` mode.  \n",
    "`PC` is well described in the paper by Handle et al. (2015), and it uses a Bayesian approach\n",
    "with the nested sampling. It works also on parameter bounds, but it would be better to\n",
    "limit the boundaries, and not used them as wide as those for `GA` and `PSO`.  \n",
    "For each mode, but `PC`, `TRADES` can perform a **bootstrap** analysis to calculate the interval\n",
    "of confidence of the best-fit parameter set. We generate a set of $T_0$s and RVs\n",
    "from the fitted parameters, and we add a Gaussian noise having the calculated value\n",
    "(of $T_0$s and RVs) as the mean and the corresponding measurement error as variance,\n",
    "scaled by the $\\sqrt{\\chi^{2}_\\textrm{reduced}}$.\n",
    "We fit each new set of observables with the `LM` with default options.\n",
    "We iterate the whole process thousands of times to analyze the distribution\n",
    "for each fitted parameter.  \n",
    "\n",
    "For the mathematical and computational description see  [Borsato et al. (2014)][Borsato2014].  \n",
    "\n",
    "---\n",
    "\n",
    "### Install and Compile {#install_compile}\n",
    "\n",
    "1. `TRADES` source is available at [github.com/lucaborsato/trades](ttps://github.com/lucaborsato/trades).  \n",
    "   Download the .zip file or clone the repository (see github help).  \n",
    "2. Extract the .zip file in your drive or enter the cloned repository.\n",
    "  It should contain a `README.md` file, `bin/`, `src/`, and `pytrades/` folders.  \n",
    "  The `src/` folder should countains the following f90 source files:  \n",
    "\n",
    "    - **Module source files:** see `Makefile`\n",
    "    - **simple integration+`LM`+bootstrap main:** `trades_int_lm_bootstrap.f90`\n",
    "    - **simple grid main:** `trades_grid.f90`  \n",
    "    - **simple `PIK` and `PSO` main:** `trades_pik_pso.f90`  \n",
    "    - **Makefile:** `Makefile`  \n",
    "    - **python script to create example simulation folder:** `createSimFile.py`  \n",
    "  \n",
    "  Look at section 6. for the `pytrades/` folder (install and execution).  \n",
    "3. Edit the `Makefile` with your Fortran 90 - MPI compiler, and with the needed compiling options.  \n",
    "    - flag `CC` for the compiler to use. From the implementation of `PC` the `mpif90` must be used;  \n",
    "    - flag `CFLAGS` for the compiler options;  \n",
    "      from the implementation of `PC` the `-ccp` preprocessor option must be used   \n",
    "      add options or uncomment following rows for other debugging options  \n",
    "    - flag `COPT` for the compiler optimization  \n",
    "      from the implementation of `PC` the `-ccp` preprocessor option must be used   \n",
    "    - flag `CFLAGS2` for the opemMP version  \n",
    "    - flag `TARGET_SER_X` is relative path and executable name for the serial program   \n",
    "    - flag `TARGET_OMP_X` is relative path and executable name for the `openMP` parallel program  \n",
    "4. To compile:  \n",
    "  in the `src/` folder type `make` and with `Tab` the user can see all the possible compile options. Check the `Makefile` for more information.  \n",
    "  \n",
    "  **Remember to type `make clean` to remove `*.o` and `*.mod` files before re-compiling `TRADES`.**  \n",
    "  **Remember to type `make cleanall` to remove `*.o`, `*.mod`, and all the executable files before re-compiling `TRADES`.**  \n",
    "  **To compile in parallel mode the `openMP` libraries must be properly installed (as suggested by your Linux distribution).**  \n",
    "  Please, change properly your `Fortran` compiler in the `Makefile`.  \n",
    "  \n",
    "  **COMMON WAY TO COMPILE**\n",
    "  ```bash\n",
    "  cd /path/to/trades/src\n",
    "  make cleanall\n",
    "  make full_parallel_release\n",
    "  ```\n",
    "  In `/path/to/trades/pytrades` there will be the `*.so`  and `pytrades.py` files to be loaded with python.  \n",
    "\n",
    "---\n",
    "\n",
    "### How to run TRADES {#run_trades}\n",
    "\n",
    "Main way to run `TRADES` with `python+PyDE+emcee`:  \n",
    "\n",
    "```bash\n",
    "  cd /path/to/exoplanet/system/path\n",
    "  python /path/to/trades/pytrades/trades_emcee.py --input configuration.yml\n",
    "```\n",
    "\n",
    "After it finished you can obtain the analysis and human-readable output with:  \n",
    "\n",
    "```bash\n",
    "  cd /path/to/exoplanet/system/path/\n",
    "  python /path/to/trades/pytrades/trades_emcee_analysis.py --input configuration.yml\n",
    "```\n",
    "\n",
    "An example of a `configuration.yml` file in [configuration.yml](trades_example/base_2p/configuration.yml).  \n",
    "The folder `/path/to/exoplanet/system/path/` must have the files described later in _Files needed by TRADES_  \n",
    "\n",
    "Other different ways to launch TRADES:\n",
    "\n",
    "- export the path of the executables in your `~/.bahsrc` or `~/.profile`:  \n",
    "  > export PATH=$PATH:/path/to/trades/executables  \n",
    "- it is possible to execute trades from the `bin/` folder by typing:  \n",
    "  > `./trades_s_xxx`  \n",
    "  > `./trades_o_xxx`  \n",
    "\n",
    "**WARNING:**  \n",
    "Before running TRADES in parallel with `open-MP` (`trades_o_xxx`) remember to set the number of cpus (`Ncpu`) to use in the `arg.in` file or by exporting:  \n",
    "`OMP_NUM_THREADS=Ncpu`  \n",
    "`export OMP_NUM_THREADS`  \n",
    "and type it in a terminal; in short way:  \n",
    "`export OMP_NUM_THREADS=Ncpu`  \n",
    "\n",
    "If TRADES has been launched without any arguments, it will search for the needed files in the current folder:  \n",
    "e.g.:  \n",
    "\n",
    "```bash\n",
    "`cd /home/user/Simulation/`  \n",
    "`trades_s_xxx`  \n",
    "```\n",
    "\n",
    "it is equal to type:  \n",
    "\n",
    "```bash\n",
    "`cd /home/user/Simulation/`  \n",
    "`trades_s_xxx .`  \n",
    "```\n",
    "\n",
    "or:  \n",
    "\n",
    "```bash\n",
    "`cd /home/user/`  \n",
    "`trades_s_xxx /home/user/Simulation/`  \n",
    "```\n",
    "\n",
    "In any of these three cases, `TRADES` will write the output files in the folder `/home/user/Simulation/`  \n",
    "  \n",
    "Provide orbital elements (and boundaries) in the bodies file.  \n",
    "Algorithm selection in `arg.in`:  \n",
    "`progtype`:\n",
    "\n",
    "1. grid  \n",
    "2. integration + `LM` + bootstrap  \n",
    "3. `PIKAIA` (+ `LM` + bootstrap)  \n",
    "4. `PSO` (+ `LM` + bootstrap)  \n",
    "\n",
    "---\n",
    "\n",
    "### Files needed by TRADES {#in_files}\n",
    "\n",
    "In the `src/` folder you can find the `createSimFile.py` python script that allows\n",
    "to create all the files needed by `TRADES`. The files are based on Kepler-9 system,\n",
    "with the original data ($T_0$ and RV) from the discovery paper by Holman et al. (2010).  \n",
    "  \n",
    "List of the files with explanation:  \n",
    "`arg.in bodies.lst star.dat b.dat c.dat`  \n",
    "`lm.opt pikaia.opt pso.opt`  \n",
    "`obsRV.dat NB2_observations.dat NB3_observations.dat`  \n",
    "\n",
    "1. `arg.in`: file with program arguments, needed for the integration, fitting type, output files.  \n",
    "  Example file [`arg.in`](trades_example/Kepler-9_example/arg.in).  \n",
    "\n",
    "2. `bodies.lst`: file with list of the files with the parameters for each body.  \n",
    "  The first column is always the file name of the body, followed by `0` or `1` for each parameter. `0` means do not fit it, `1` means fit it.  \n",
    "  The first row is always the star file with the Mass and Radius fitting parameter type.  \n",
    "  From the second row, each line is the body file name followed by the parameters to fit in this order:  \n",
    "  `mass radius period eccentricity argument_of_pericenter mean_anomaly inclination longitude_of_node`.  \n",
    "  Remember that the number of the lines of this file has to match the `NB` parameter in the `arg.in` file.  \n",
    "  Example file [`bodies.lst`](trades_example/Kepler-9_example/bodies.lst)  \n",
    "\n",
    "    - `star.dat`: Mass (and sigma) and Radius (and sigma) of the star in Solar units. First row the Mass, second the Radius.  \n",
    "    In the code will be identified with the id == 1.  \n",
    "    Example file [`star.dat`](trades_example/Kepler-9_example/star.dat)  \n",
    "    - `b.dat`: file with parameters of the planet in the second row of the bodies.lst, that is in the code will be identified with the id == 2.  \n",
    "    Each row is a different parameter (3 columns + flag for `grid`), in the order:  \n",
    "    `mass min max [M_Jup]`  \n",
    "    `radius min max [R_Jup]`  \n",
    "    `period min max [days]`  \n",
    "    `semi-major_axis min max [au]`  \n",
    "    `eccentricity min max`  \n",
    "    `argument_of_pericenter min max [deg]`  \n",
    "    `mean_anomaly min max [deg]`  \n",
    "    `time_of_pericenter_passage min max [JD]`  \n",
    "    `inclination min max [deg]`  \n",
    "    `longitude_of_node min max [deg]`  \n",
    "    When not specified the `min max` values are set by default (i.e. eccentricity < 1, angles between 0 and 360 deg, inclination between 0 and 180 deg, radius < 5 `R_Jup`, and mass < 1 `M_Sun` in `M_Jup`.  \n",
    "    In the `grid` mode the mean of the columns is: `min max step type`. The `step` is used as true step, or step number, and so on accordingly to the fourth column `type`.  \n",
    "    Keywords for column 4 are: `ss` (means step size), `rn` (random number), and `sn` (step number).  \n",
    "    Alternatives:  \n",
    "    `semi-major axis` in au instead of `period` (set `period` greater than `9000000.`, while setting semi-major axis equal to `999.` will let you use the period);  \n",
    "    `time of passage at pericentre` is used if `mean anomaly` is greater than `999.`, while set `time of passage at pericentre` greater than `9.e8` to use `mean anomaly`.\n",
    "    Example file [`b.dat`](trades_example/Kepler-9_example/b.dat)  \n",
    "    - `c.dat`: same as file `b.dat`, but with different parameter values for the body in the third row in `bodies.lst`, that is in the code will be identified with the id == 3.  \n",
    "    Example file [`c.dat`](trades_example/Kepler-9_example/c.dat)  \n",
    "\n",
    "3. `lm.opt`: parameter options for the Levenberg-Marquardt algorithm; they are based on the original manual. Keep it as it for standard analysis.  \n",
    "  Example file [`lm.opt`](trades_example/Kepler-9_example/lm.opt)  \n",
    "\n",
    "4. `pikaia.opt`: parameter options for the `GA` algorithm. The most important parameters to tune are  \n",
    "  (1pik) the number of individuals (row 1, ctrl(1)), that is the number of set of parameters for each generation;  \n",
    "  (2pik) the number of generation (row 2, ctrl(2)), that is the number of iteration that `GA` has to perform, the last iteration returns the best set of parameters;  \n",
    "  (3pik) the seed (row 13), that is a integer number that defines the seed for the random generator, if you keep the same value it repeat the same analysis;  \n",
    "  (4pik) the wrtAll (row 14) is a parameter that defines if you want to that the `GA` writes all the individuals for each generation, set it to `1` to write, `0` not write;  \n",
    "  (5pik) the nGlobal (row 15) is the number of global search to perform with the `GA`, each search returns a solution, and the seed of each analysis is different (seed + i, i=1..nGlobal).  \n",
    "  Example file [`pikaia.opt`](trades_example/Kepler-9_example/pikaia.opt)  \n",
    "\n",
    "5. `pso.opt`: parameter options for the `PSO` algorithm. The rows 1, 2, 4, 5, and 6 are the same parameters explained for the `pikaia.opt` file.  \n",
    "  In particular:  \n",
    "  (1pso) row 1 and (2pik) 2 are exactly the same as in `pikaia.opt`;  \n",
    "  (3pso) row 3 is an integer that specifies when write a summary during the `PSO` analysis;  \n",
    "  (4pso) row 4 is the same as (4pik) row 14 in `pikaia.opt`;  \n",
    "  (5pso) row 5 is the same as (5pik) row 15 in `pikaia.opt`;  \n",
    "  (6pso) row 6 is the same as (3pik) row 13 in `pikaia.opt`.  \n",
    "  Example file [`pso.opt`](trades_example/Kepler-9_example/pso.opt)  \n",
    "\n",
    "6. `obsRV.dat`: list of radial velocities (RVs) data.  \n",
    "  Columns description:  \n",
    "  (1) RV observation time (JD, or time in same units of the integration time);  \n",
    "  (2) observed RVs in meter per seconds;  \n",
    "  (3) observed RV uncertainties in meter per seconds;  \n",
    "  (4) ID of the RV dataset, so if you have only one dataset set all column to `1`, else increas value untill the number of different datasets, i.e., `1`, `2` for 2 datasets (2 different facilities, or one facility before and after upgrade).  \n",
    "  Example file [`obsRV.dat`](trades_example/Kepler-9_example/obsRV.dat)  \n",
    "\n",
    "7. `NB2_observations.dat`: list of transit times (`$T_0$s`) observed for planet in the second row of `bodies.lst`, i.e., `b.dat` is planet 2.  \n",
    "  Columns description:  \n",
    "  (1) transit epoch (N), an integer number that identifies the transit w.r.t. a reference transit time `$T_\\textrm{ref}$` and refined by a linear ephemeris of kind: `$T_N = T_\\textrm{ref} + P_\\textrm{ref} \\times N$`;  \n",
    "  (2) the transit time (`$T_0$`) in JD or the same time unit of the integration/epoch/start time;  \n",
    "  (3) the uncertainty on the `$T_0$`.\n",
    "  Example file [`NB2_observations.dat`](trades_example/Kepler-9_example/NB2_observations.dat)  \n",
    "\n",
    "8. `NB3_observations.dat`: same kind of file `NB2_observations.dat`, but for the planet in the third row of `bodies.lst`, i.e., `c.dat` is planet 3.  \n",
    "  Example file [`NB3_observations.dat`](trades_example/Kepler-9_example/NB3_observations.dat)  \n",
    "  \n",
    "9. `derived_boundaries.dat`: this file a special file.  \n",
    "  If you have some derived parameters (or other values) that can reduce your parameter space you have to create this file in your simulation folder.  \n",
    "  if the argument `secondary_parameters` in `arg.in` file is set to `0`, this file will not be used.  \n",
    "  If `secondary_parameters = 1`, but the file does not exist it will not be used (no derived parameters will be checked).  \n",
    "  The file should have a line for each parameter, the name in the first column (please keep it short), the min and the max value in the 2nd and 3rd column.\n",
    "  Keep last line empty so the code can determine the end of file.  \n",
    "  `derived_boundaries.dat` example:  \n",
    "  `# name phase_min phase_max`  \n",
    "  `ph2 34. 180.`  \n",
    "  `ph3 180. 270.`  \n",
    "  In order to use the derived parameters you have to modify by your own the `derived_parameters_mod.f90`.\n",
    "  In fitness_module.f90 check which subroutines/fuctions will be called to 'check_derived' or 'fix_derived' parameters.  \n",
    "\n",
    "10. `priors.in`: place here priors for Bayesian analysis with `emcee`.  \n",
    "  Follow the instructions in the file [`priors.in`](trades_example/base_2p/priors.in).  \n",
    "\n",
    "<!-- 11. `de.yml`: file `yaml` used for `PyDE` configuration. See an example in the file [`de.yml`](trades_example/base_2p/de.yml).   -->\n",
    "\n",
    "---\n",
    "\n",
    "### Output files by TRADES {#out_files}\n",
    "\n",
    "Each algorithm will write different files, and depends on the flag used in the `arg.in` and in the `*.opt` files.  \n",
    "Each file _should_ have an self-explaning header.  \n",
    "\n",
    "1. **integration:** depends only on `arg.in` file  \n",
    "  `#ID_#LM_rotorbit.dat`, `#ID_#LM_constants.dat`, `#ID_#LM_NB#_elements.dat`, `#ID_#LM_NB#_tra.dat`, `#ID_#LM_gls_output.dat`  \n",
    "  `#ID_#LM_rotorbit.dat`, if `wrtorb = 1`, where `#ID` is the simulation ID, `#LM` is the Levenberg-Marquardt flag (`lmon = 0` or `1`). Columns: 1 Time in JD; 2 Light-Time Travel Effect in days (LTE\\_d); 3:3+NB\\*6 {X,Y,Z,VX,VY,VZ} for each body (NB=number of bodies); last column is the radial velocity (RV) of the star due to the planets in m/s.  \n",
    "  `#ID_#LM_constants.dat`, if `wrtcon = 1`, naming convenction as previous. Columns: 1 Time in JD; 2 momentum; 3 delta between initial and current momentum; 4 Total Energy; 5 delta between initial and current Total Energy.  \n",
    "  `#ID_#LM_gls_output.dat`, it is the output of the General Lomb-Scargle [(GLS)][GLS] applied to the residuals: $\\textrm{res}=\\textrm{RV}_\\textrm{obs}-\\textrm{RV}_\\textrm{sim}$.  \n",
    "  `#ID_#LM_NB#_elements.dat`, if `wrtel = 1`, naming convenction as previous, plus the body id NB#, starting from 2 to the number of bodies used. Columns: 1 Time in JD; 2 Period in days, 3 semi-major axis in astronomical unit (au), 3 eccentricity, 4 inclination in degrees, 4 mean anomaly in degrees, 5 argument of the pericenter in degrees, 6 longitude of the node in degrees, 7 true anomaly in degrees, 8 difference between time of refence (epoch) and time of the passage of pericenter tau in days  \n",
    "  `#ID_#LM_NB#_tra.dat`, if `idtra > 0`, naming convenction as previous. Columns: 1 trasit time, 2 LTE, 3 firt contact time, 4 second contact time, 5 third contact time, 6 fourth contact time, 7:7+NB*6 state vector {X,Y,Z,VX,VY,VZ} for each body (NB=number of bodies).  \n",
    "  \n",
    "---\n",
    "\n",
    "### Python Library {#pylib}\n",
    "\n",
    "`TRADES` can now be used in python scripts, thanks to `f2py`.  \n",
    "In the `src/` folder there is the `pytrades_lib.f90` that is the library that interface `TRADES` to python.  \n",
    "**Compile:** see `Makefile`  \n",
    "  \n",
    "It creates the `pytrades_lib.so` and it copies it into the folder `../pytrades/`  \n",
    "Check if in the `src/` folder there is an hidden file `.f2py_f2cmap` that contains one row: `dict(real=dict(sp='float', dp='double'))`  \n",
    "**Run scripts:**  \n",
    "<!--In folder `pytrades/`: `pytrades_lib.so`, `constants.py`, `ancillary.py`, `trades_pso2emcee.py`.  \n",
    "The script to run is `trades_pso2emcee.py`, that shows how to call the `pytrades_lib` and how to combine `PSO` algorithm with `emcee`.  \n",
    "Run as `python trades_pso2emcee.py -h` for instruction on the all the command line arguments to provide.  -->\n",
    "In folder `pytrades/` there are different scripts, look for the script to run beginning with `trades_`, that shows how to call the `pytrades_lib` and how to combine with `emcee`.  \n",
    "Run as `python trades_XXX.py -h` for instruction on the all the command line arguments to provide.  \n",
    "The `ancillary.py` file has different functions to write and read output files from `PSO/PyDE` and `emcee` and other stuff helpful to manage the output of the simulations.  \n",
    "The `full_emcee_analysis.py` script is able to create all the summary reports and\n",
    "images of a `TRADES+EMCEE` analysis; it works also on a running simulation.  \n",
    "Check also other files in the `pytrades` folder to understand how to use the python library.\n",
    "\n",
    "I am working on some scripts and notebook to use some possibility of `TRADES`, such as an implementation that could be used within [`PyORBIT`][PyORBIT], see folder [python_examples](trades_example/python_examples)  \n",
    "\n",
    "**TO BE CONTINUED**  \n",
    "\n",
    "---\n",
    "\n",
    "### Changes/Log\n",
    "\n",
    "**sorry, I will not be able to report all the small changes...**  \n",
    "\n",
    "#### `TRADES 2.19.0`  \n",
    "\n",
    "Better management of ode subroutines and transit times detection and computation.  \n",
    "Main python script files:  \n",
    "`trades_emcee.py` that run [`PyDE`][PyDE] and [`emcee`][emcee]  \n",
    "`trades_emcee_analysis.py` to analyse the [`emcee`][emcee] output  \n",
    "Added the [](trades_example/python_examples/) folder with `jupyter notebook` to show some python ways to use `TRADES`.  \n",
    "\n",
    "#### `TRADES 2.18.0`  \n",
    "\n",
    "`TRADES` now fit the combination of `ecccentricity` ($e$) and `argument of pericenter` ($\\omega$) as\n",
    "$(\\sqrt{e}\\cos\\omega, \\sqrt{e}\\sin\\omega)$ (also in `PSO` and `GA`).  \n",
    "Added the script that uses [`PyDE`][PyDE] + [`emcee`][emcee].  \n",
    "The default fitness parameter for the `PSO` is now the $\\log{\\mathcal{L}}$.\n",
    "Added python modules to plot the evolution of the `PSO` and `PyDE`.  \n",
    "Converted example script `trades_test_2022-03-04.py` to `trades_python_test_read_data.py` and created `trades_python_test.py`.  \n",
    "\n",
    "#### `TRADES 2.17.1`  \n",
    "\n",
    "`pytrades_lib.f90` now has `kelements_to_data` subroutine that returns by default also the duration (T4 - T1, equal to zero if passed dur_check=0)\n",
    "and the keplerian orbital elements at each transit time.  \n",
    "The `pytrades/trades_test_2022-03-04.py` shows how to use it and it is based on the Kepler-9 example in `trades_example/Kepler-9_example`.  \n",
    "\n",
    "\n",
    "#### `TRADES 2.17.0`  \n",
    "\n",
    "New management of how it determines the transiting body to check during the integration of the orbits.  \n",
    "Removed the parameters (also in the `arg.in`) for weight the chi square in the fit.  \n",
    "Now fitting RV offset (called gamma) and jitter parameters for each RV data-set,\n",
    "and it can fit a trend on RV of specified order.  \n",
    "Some better management of plots and confidence intervals python scripts.  \n",
    "Updates in the `Makefile`.  \n",
    "Testing `ultranest` version.  \n",
    "\n",
    "#### `TRADES 2.16.0`  \n",
    "\n",
    "Moved unused python files in two folders (do not used them, just there as control).  \n",
    "Updated `plot_oc/rv.py` files for better x/y limits and visualisation.  \n",
    "Updated `confidence_intervals.py` with better HDI computation (based on [PyAstronomy][PyAstronomy]), mode determination.  \n",
    "`PyPolyChord` not tested anymore, moved all files in old/control folders.  \n",
    "Added `priors.in` file for Bayesian analysis with [emcee][emcee].  \n",
    "\n",
    "#### `TRADES 2.15.2`  \n",
    "\n",
    "Fixed bug in Makefile that does create a segfault running TRADES executables.  \n",
    "\n",
    "#### `TRADES 2.15.1`  \n",
    "\n",
    "Changes:  \n",
    "\n",
    "- moved `pytrades` from `python 2` to `python 3`;\n",
    "- removed unused variables;\n",
    "- added a subroutine/function to get only the reduced chi square (`$\\chi^{2}_{r}$`).  \n",
    "\n",
    "#### `TRADES 2.15.0`  \n",
    "\n",
    "Changed how TRADES stores and manages the observed and simulated data (Transit Times and Radial Velocities)\n",
    "and it uses new data type.  \n",
    "TRADES can now fit the transit duration in minutes defined as the difference between \n",
    "the last (T4) and the first (t1) contact time.  \n",
    "Improved the T0 and the T4/T1 determination, and now it uses the integrator in the Bisection-Newton-Raphson,\n",
    "in this way it avoids numericals issues with analytical f and g function when \n",
    "assuming constants keplerian orbital elements during transit.  \n",
    "\n",
    "#### `TRADES 2.14.1`  \n",
    "\n",
    "Updated computation of the initial state vector from the Keplerian orbital elements,\n",
    "and implemented the subroutine adapted from [`MERCURY`][MERCURY] to compute the \n",
    "Keplerian orbital elements from the state vector.  \n",
    "\n",
    "#### `TRADES 2.14`  \n",
    "\n",
    "Old `PolyChord` version has been removed, now the 1.10 version can be used.\n",
    "Follow the `PolyChord` compile and install instruction, with python support.\n",
    "It will be available as python module and used in the script ``.  \n",
    "  \n",
    "Checking the conversion of the state vector to orbital elements during integration.  \n",
    "\n",
    "#### `TRADES 2.13.6`\n",
    "\n",
    "Check few variable allocations when no data available.  \n",
    "Ready to move PolyChord to v1.12 with python support.  \n",
    "\n",
    "#### `TRADES 2.13.5`\n",
    "\n",
    "Now TRADES uses the `pool` in `emcee` instead of `threads` to avoid memory overloading\n",
    "in loops.  \n",
    "\n",
    "#### `TRADES 2.13.4`\n",
    "\n",
    "It is possible to save some parameter samples of the `emcee` analysis to be overplotted\n",
    "on the O-C and RV diagram.  \n",
    "Updated how it computes the orbits when the write step is smaller than the integration step.  \n",
    "\n",
    "#### `TRADES 2.13.3`\n",
    "\n",
    "Nothing special, code-styling and better management of planet file names.  \n",
    "\n",
    "#### `TRADES 2.13.2`\n",
    "\n",
    "In `arg.in` file added keyword: `ncpu`.  \n",
    "Set it to define the number of cpus to use in parallel with `TRADES`.  \n",
    "The flag `-c` passed to the python scripts override/overwrite it.  \n",
    "Updated the initialisation of the particles in `PSO` and `PIK`.   \n",
    "\n",
    "#### `TRADES 2.12.1`\n",
    "\n",
    "Python scripts: `trades_pso2emcee.py` and `trades_pso2emcee_sqrte.py`.  \n",
    "Counterparts of `trades_emcee.py` and `trades_emcee_sqrte.py`, but it runs or reads `PSO` before `emcee`.  \n",
    "Testing.  \n",
    "\n",
    "#### `TRADES 2.11.1`\n",
    "\n",
    "Python script `trades_emcee_sqrte.py`, that allow to fit `$\\sqrt{e}\\cos\\omega,\\sqrt{e}\\sin\\omega$` (in the python module, not in the fortran subroutines).  \n",
    "Added script `plot_rv.py` that plot RV model and RV-phased curves.  \n",
    "Now the degrees of freedom (dof) take into account the free parameters (nfre, determined within the program, but not fitted). In particular they are the number of `$\\gamma$` RV offset (`nfree = n_RVset`).  \n",
    "`dof = ndat - nfit - nfree`\n",
    "\n",
    "#### `TRADES 2.10.1`\n",
    "\n",
    "More updates on python scripts.  \n",
    "\n",
    "#### `TRADES 2.10.0`\n",
    "\n",
    "From the python library (`pytrades_lib`) it is possible to initialise `TRADES` without read files and given a set of orbital parameters, times of radial velocities, and transit epochs (integer numbers that identify each transit w.r.t. a transit of reference) it returns the simulated radial velocities and transit times.  \n",
    "To use this feature I added in the folder `pytrades` a test file `trades_test_20161116.py`. Be aware that the radial velocities and transit times have been read from 2 files, so modify it accordingly to your file/data.  \n",
    "Added more files in `pytrades` folder and added files to plot and summarise simulation with `emcee`.  \n",
    "\n",
    "#### `TRADES 2.9.1`\n",
    "\n",
    "Small bugfixes in check of the physical parameters.  \n",
    "Added optional argument for debug in the output subroutine in `ode_run`.   \n",
    "\n",
    "#### `TRADES 2.9.0`\n",
    "\n",
    "Introduced flag in `arg.in` file: `oc_fit`.  \n",
    "`oc_fit = F` is false, it means that the difference between the observed transit times (`T0_obs`) and the simulated (`T0_sim`) is used to compute the `$\\chi^2$`;  \n",
    "`oc_fit = T` is true, it means that `$\\chi^2$` is computed from the amplitude of the TTV, based on the O-C for the observed (`OC_obs`) and simulated data (`OC_sim`):  \n",
    "  `OC_obs = T0_obs - T0_lin,obs`, where `T0_lin,obs` is the linear ephemeris of the observed data and  \n",
    "  `OC_sim = T0_sim - T0_lin,sim`, where `T0_lin,sim` is the linear ephemeris of the simulated transit times.  \n",
    "\n",
    "#### `TRADES 2.8.0`\n",
    "\n",
    "Python library of `TRADES` has been tested and now is public.  \n",
    "  \n",
    "In `bodies.lst` file user can specify if a planet should transit or not.  \n",
    "At the end of each planet row (after 0/1 fit-flag and before the #), user adds a T (or leaves blank) if planet is allowed to transit, otherwise adds an F.  \n",
    "Example:  \n",
    "in `bodies.lst` row `b.dat 1 0 0 0 0 0 0 0 T` is equal to `b.dat 1 0 0 0 0 0 0 0`: planet b could transit.   \n",
    "in `bodies.lst` row `b.dat 0 0 0 1 0 0 0 0 F`: planet b could not transit (in case it discards the solution).  \n",
    "  \n",
    "Converted all logical values in a consistent way.  \n",
    "\n",
    "#### `TRADES 2.7.0`\n",
    "\n",
    "Fitting parameters examples:  \n",
    "in `bodies.lst` row `b.dat 1 0 0 0 0 0 0 0`:\n",
    "`mass_planet / mass_star`  \n",
    "in `bodies.lst` row `b.dat 0 0 0 1 0 0 0 0`:\n",
    "`eccentricity`  \n",
    "in `bodies.lst` row `b.dat 1 0 0 1 1 0 0 0`:\n",
    "`mass_planet / mass_star`, `e*cos(w)`, `e*sin(w)`  \n",
    "in `bodies.lst` row `b.dat 0 0 0 0 1 1 1 0`:\n",
    "`argument_of_pericenter`, `mean longitude (lambda)`, `inclination`  \n",
    "in `bodies.lst` row `b.dat 1 0 0 1 1 1 1 1`:\n",
    "`mass_planet / mass_star`, `e*cos(w)`, `e*sin(w)`, `lambda`, `i*cos(lN)`, `i*sin(lN)`  \n",
    "\n",
    "\n",
    "#### `TRADES 2.6.0`\n",
    "\n",
    "The old main `trades.f90` has been replaced (still in the sources, it compiles)\n",
    "by small easier-to-understand mains.\n",
    "As for the version `2.5.1` the user has to change the `arg.in` accordingly to the purpose.\n",
    "See [](How to run TRADES).  \n",
    "During `RV` fit, `TRADES` will run a `GLS` periodogram [(Zechmeister and Kurster, 2009)][GLS] and look for peaks close (`$\\Delta P = \\pm 0.5\\$ d`) to the periods of the planets of the simulated system.\n",
    "In case of a positive signal, it means the period has been induced (bad `RV` fit) and the fitness (`$\\chi^{2}_\\textrm{r}$`) will be set to max value (bad fit).  \n",
    "  \n",
    "---  \n",
    "\n",
    "[emcee]: https://emcee.readthedocs.io/en/stable/  \n",
    "[TRADES@github]: https://github.com/lucaborsato/trades  \n",
    "[Borsato2014]: https://ui.adsabs.harvard.edu/abs/2014A%26A...571A..38B/abstract  \n",
    "[Malavolta2017]: https://ui.adsabs.harvard.edu/abs/2017AJ....153..224M/abstract  \n",
    "[Borsato2019]: https://ui.adsabs.harvard.edu/abs/2019MNRAS.484.3233B/abstract  \n",
    "[TRADESESPG]: http://groups.dfa.unipd.it/ESPG/trades.html  \n",
    "[MINPACK]: http://www.netlib.org/minpack/  \n",
    "[Moré1980]: http://www.mcs.anl.gov/~more/ANL8074a.pdf  \n",
    "[PIKAIA]: http://www.hao.ucar.edu/modeling/pikaia/pikaia.php  \n",
    "[PSO]: http://www.nda.ac.jp/cc/users/tada/  \n",
    "[GLS]: http://adsabs.harvard.edu/abs/2009A%26A...496..577Z  \n",
    "[MERCURY]: http://www.arm.ac.uk/~jec/  \n",
    "[PyAstronomy]: https://www.hs.uni-hamburg.de/DE/Ins/Per/Czesla/PyA/PyA/index.html#  \n",
    "[PyDE]: https://github.com/hpparvi/PyDE  \n",
    "[PyORBIT]: https://github.com/LucaMalavolta/PyORBIT  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fe2d037-805e-4b1f-81c6-e04cfd11e377",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}