{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic pyerrors example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import pyerrors, as well as autograd wrapped numpy and matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pyerrors as pe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('./base_style.mplstyle')\n",
    "plt.rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use numpy to generate some fake data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_sample1 = np.random.normal(2.0, 0.5, 1000)\n",
    "test_sample2 = np.random.normal(1.0, 0.1, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this we can construct `Obs`, which are the basic object of `pyerrors`. For each sample we give to the obs, we also have to specify an ensemble/replica name. In this example we assume that both datasets originate from the same gauge field ensemble labeled 'ens1'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs1 = pe.Obs([test_sample1], ['ensemble1'])\n",
    "obs2 = pe.Obs([test_sample2], ['ensemble1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now combine these two observables into a third one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs3 = np.log(obs1 ** 2 / obs2 ** 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pyerrors` overloads all basic math operations, the user can work with these `Obs` as if they were real numbers. The proper resampling is performed in the background via automatic differentiation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are now interested in the error of `obs3`, we can use the `gamma_method` to compute it and then print the object to the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.367(20)\n"
     ]
    }
   ],
   "source": [
    "obs3.gamma_method()\n",
    "print(obs3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the method `details` we can take a look at the integrated autocorrelation time estimated by the automatic windowing procedure as well as the detailed content of the `Obs` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 1.36706932e+00 +/- 2.04253682e-02 +/- 1.02126841e-03 (1.494%)\n",
      " t_int\t 5.01998002e-01 +/- 4.47213595e-02 S = 2.00\n",
      "1000 samples in 1 ensemble:\n",
      "  · Ensemble 'ensemble1' : 1000 configurations (from 1 to 1000)\n"
     ]
    }
   ],
   "source": [
    "obs3.details()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected the random data from numpy exhibits no autocorrelation ($\\tau_\\text{int}\\,\\approx0.5$). It can still be interesting to have a look at the window size dependence of the integrated autocorrelation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 15456 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This figure shows the windowsize dependence of the integrated autocorrelation time. The red vertical line signalizes the window chosen by the automatic windowing procedure with $S=2.0$.\n",
    "Choosing a larger windowsize would not significantly alter $\\tau_\\text{int}$, so everything seems to be correct here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlated data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now generate fake data with a given covariance matrix and integrated autocorrelation times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov = np.array([[0.5, -0.2], [-0.2, 0.3]])  # Covariance matrix\n",
    "tau = [4, 8]  # Autocorrelation times\n",
    "c_obs1, c_obs2 = pe.misc.gen_correlated_data([2.8, 2.1], cov, 'ens1', tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and once again combine the two `Obs` to a new one with arbitrary math operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 3.27194697e-01 +/- 1.53249111e+00 +/- 2.49471479e-01 (468.373%)\n",
      " t_int\t 4.75187177e+00 +/- 1.33949719e+00 S = 2.00\n",
      "1000 samples in 1 ensemble:\n",
      "  · Ensemble 'ens1' : 1000 configurations (from 1 to 1000)\n"
     ]
    }
   ],
   "source": [
    "c_obs3 = np.sin(c_obs1 / c_obs2 - 1)\n",
    "c_obs3.gamma_method()\n",
    "c_obs3.details()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time we see a significant autocorrelation so it is worth to have a closer look at the normalized autocorrelation function (rho) and the integrated autocorrelation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22852 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 17424 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c_obs3.plot_rho()\n",
    "c_obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now redo the error analysis and alter the value of S or attach a tail to the autocorrelation function to take into account long range autocorrelations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 3.27194697e-01 +/- 1.67779862e+00 +/- 2.08884244e-01 (512.783%)\n",
      " t_int\t 5.69571763e+00 +/- 2.09295390e+00 tau_exp = 20.00,  N_sigma = 1\n",
      "1000 samples in 1 ensemble:\n",
      "  · Ensemble 'ens1' : 1000 configurations (from 1 to 1000)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23544 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c_obs3.gamma_method(tau_exp=20)\n",
    "c_obs3.details()\n",
    "c_obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
