# aderrors - Error analysis of Monte Carlo data with Automatic Differentiation
`aderrors` is a fortran implementation of the $`\Gamma`$-method for
analysis of Monte Carlo data. It uses Automatic Differentiation to
perform **exact linear error propagation**. It preforms the computation
of gradients and Hessians of arbitrary functions, which allows a
robust and **exact error propagation even in iterative algorithms**.
- [Features](#features)
- [Examples](#examples)
- [Simple analysis of MC data](#simple-analysis-of-mc-data)
- [A complete example](#a-complete-example)
- [A calculator with uncertainties](#a-calculator-with-uncertainties)
- [Installation](#installation)
- [Requirements](#requirements)
- [Instructions](#instructions)
- [Using the library](#using-the-library)
- [Full documentation](#full-documentation)
- [How to cite](#how-to-cite)
## Features
- **Exact** linear error propagation, even in iterative algorithms
(i.e. error propagation in fit parameters).
- Handles data from **any number of ensembles** (i.e. simulations with
different parameters).
- Support for **replicas** (i.e. several runs with the same simulation
parameters).
- Irrelagular MC measurements are handled transparently.
- Standalone **portable** implementation without any external
dependencies.
- **Fast** computation of autocorrelation functions with the [FFT
package](http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html) (included
in the distribution).
- **Exact** determination of gradients and Hessians of arbitrary
functions.
## Examples
This is just a small collection of examples. Note that the basic data type is
`uwreal`, that is able to handle MC histories and data with
errors. The error is determined by calling the method `uwerr`
on the data type. More examples can be found in the `test` directory
of the distribution. They are explained in the documentation `doc/aderrors.pdf`.
All these examples use the `module simulator` (available in
`test/simulator.f90`) to generate autocorrelated data with
autocorrelation function $`\Gamma(t) = \sum_k \lambda_k
e^{-|t|/\tau_k}`$.
### Simple analysis of MC data
This code performs a simple error analysis of the data with the default
parameters ($`S_{\tau}=4`$ for automatic windowing, no tail added to the
autocorrelation function). The result of the analysis with the
$`\Gamma`$-method is compared with the exact values of the error and
$`\tau_{\rm int}`$ returned by the `module simulator`. This example is
included in the distribution in the file `test/simple.f90`.
```fortran
program simple
use ISO_FORTRAN_ENV, Only : error_unit, output_unit
use numtypes
use aderrors
use simulator
implicit none
integer, parameter :: nd = 20000
type (uwreal) :: x
real (kind=DP) :: data_x(nd), err, ti, tau(5), lam(5)
! Fill arrays data_x(:) with autocorrelated
! data from the module simulator.
tau = (/1.0_DP, 3.0_DP, 4.0_DP, 5.0_DP, 7.34_DP/)
lam = (/1.00_DP, 0.87_DP, 1.23_DP, 0.56_DP, 0.87_DP/)
call gen_series(data_x, err, ti, tau, lam, 0.3_DP)
! Load data_x(:) measurements in variable x. Use
! default settings (Stau=4, texp=0, 1 replica)
x = data_x
! Perform error analysis (optimal window)
call x%uwerr()
! Print results
write(*,'(1A,1I6,1A)')'** Measurements: ', nd, ' ** '
write(*,100)' - Gamma-method: ', x%value(), " +/- ", x%error(), '( tauint: ', &
x%taui(1), " +/- ", x%dtaui(1), ')'
write(*,100)' - Exact: ', 0.3_DP, " +/- ", err, '( tauint: ', ti, " )"
100 FORMAT((2X,1A,1F8.5,1A,1F7.5,5X,1A,1F0.2,1A,1F7.5,1A))
stop
end program simple
```
Running this code gives as output
```
** Measurements: 20000 **
- Gamma-method: 0.29883 +/- 0.04271 ( tauint: 4.13 +/- 0.48344)
- Exact: 0.30000 +/- 0.04074 ( tauint: 3.82 )
```
### A complete example
This example is included in the distribution in the file
`test/complete.f90`. The analysis is performed
for a complicated non-linear function of the two primary
observables. It also uses replica for the MC ensemble labeled 1.
```fortran
program complete
use ISO_FORTRAN_ENV, Only : error_unit, output_unit
use numtypes
use constants
use aderrors
use simulator
implicit none
integer, parameter :: nd = 5000, nrep=4
type (uwreal) :: x, y, z
integer :: iflog, ivrep(nrep)=(/1000,30,3070,900/), i, is, ie
real (kind=DP) :: data_x(nd), data_y(nd/2), err, ti, texp
real (kind=DP) :: tau(4), &
lam_x(4)=(/1.0_DP, 0.70_DP, 0.40_DP, 0.40_DP/), &
lam_y(4)=(/2.3_DP, 0.40_DP, 0.20_DP, 0.90_DP/)
character (len=200) :: flog='history_z.log'
! Fill arrays data_x(:) with autocorrelated
! data from the module simulator. Use nrep replica
tau = (/1.0_DP, 3.0_DP, 12.0_DP, 75.0_DP/)
texp = maxval(tau)
is = 1
do i = 1, nrep
ie = is + ivrep(i) - 1
call gen_series(data_x(is:ie), err, ti, tau, lam_x, 0.3_DP)
is = ie + 1
end do
! Fill data_y(:) with different values of tau also using
! module simulator
forall (i=1:4) tau(i) = real(2*i,kind=DP)
call gen_series(data_y, err, ti, tau, lam_y, 1.3_DP)
! Load data_x(:) measurements in variable x.
! Set replica vector, exponential autocorrelation time
! and ensemble ID.
x = data_x
call x%set_id(1)
call x%set_replica(ivrep)
call x%set_texp(texp)
! Load data_y(:) measurements in variable y
y = data_y
call y%set_id(2)
! Exact, transparent error propagation
z = sin(x)/(cos(y) + 1.0_DP)
! Attach tail in ensemble with ID 1 when signal in the
! normalized auto-correlation function equals its error
call z%set_dsig(1.0_DP,1)
! Set Stau=3 for automatic window in ensemble with ID 2
call z%set_stau(3.0_DP,2)
! Perform error analysis (tails, optimal window,...)
call z%uwerr()
! Print results and output details to flog
write(*,'(1A,1F8.5,1A,1F8.5)')'** Observable z: ', z%value(), " +/- ", z%error()
do i = 1, z%neid()
write(*,'(3X,1A,1I3,3X,1F5.2,"%")',advance="no")'Contribution to error from ensemble ID', &
z%eid(i), 100.0_DP*z%error_src(i)
write(*,'(2X,1A,1F0.4,1A,1F8.4,1A)')'(tau int: ', z%taui(i), " +/- ", z%dtaui(i), ")"
end do
open(newunit=iflog, file=trim(flog))
call z%print_hist(iflog)
close(iflog)
stop
end program complete
```
Running this code gives as output
```
** Observable z: 0.24426 +/- 0.05374
Contribution to error from ensemble ID 1 83.93% (tau int: 5.8333 +/- 2.0772)
Contribution to error from ensemble ID 2 16.07% (tau int: 2.5724 +/- 0.5268)
```
The file `history_z.log` contains the MC histories, normalized
autocorrelation functions and $`\tau_{\rm int}`$. The format is
a simple text file: it can be processed by shell scripts
(see `tools/plot/plot_hist.sh`) to produce the following graphics.
### A calculator with uncertainties
The module also handles simple data with errors (i.e. **not** a MC
history). This is extremely useful since in many occasions we have to
use data from statistically independent sources where only the central
value and the error is available.
With the module `aderrors` these kind of data in handled
transparently, just as if it were another ensemble (see `test/calculator.f90`).
```fortran
program calculator
use aderrors
implicit none
type (uwreal) :: x, y, z, t
x = (/1.223_8, 0.012_8/) ! x = 1.223 +/- 0.012
y = sin(2.0_8*x)
z = 1.0_8 + 2.0_8 * sin(x)*cos(x)
t = z - y
call t%uwerr()
write(*,'(1A,1F18.16,1A,1F18.16)')'Exactly one: ', t%value(), " +/- ", t%error() ! 1.0 +/- 0.0
stop
end program calculator
```
Running this code gives as output
```
Exactly one: 1.0000000000000000 +/- 0.0000000000000000
```
## Installation
### Requirements
The code is strict `fortran 2008` compliant. Any compiler that
supports this standard can be used. The code has been tested with
`gfortran 6.X`, `gfortran 7.X`, `gfortran 8.X`, `intel fortran 17`,
`intel fortran 18`.
Note that `gfortran 4.X` and `gfortran 5.X` do not support
`submodules` (part of the `fortran 2008` standard). This code will not
work whith these versions.
### Instructions
1. Download or clone the repository.
1. Edit the `Makefile` in the `build` directory. Change the compiler
command/options (variables `FC` and `FOPT`).
1. Compile the library with `GNU make`.
1. Optionally build/run the test codes with `make test`. Executabes will
be placed in the `test` directory.
1. If preferred, move the contents of the `include` and `lib`
directories somewhere else.
### Using the library
1. Compile your programs with `-I