{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "24c23400fad3c24d",
   "metadata": {},
   "source": [
    "# Estimate noise model parameters\n",
    "\n",
    "Objects from `thztools`: `noisefit`, `NoiseModel`, `NoiseModel.noise_amp`, `NoiseModel.noise_sim`,`scaleshift`, `timebase`, `wave`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "initial_id",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import thztools as thz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbebea0ddbcd093e",
   "metadata": {},
   "source": [
    "## Set the simulation parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c483f358851818c",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 256  # Number of samples\n",
    "m = 50  # Number of simulated waveforms\n",
    "dt = 0.05  # Sampling time [ps]\n",
    "\n",
    "sigma_alpha = 1e-4  # Additive noise amplitude [signal units]\n",
    "sigma_beta = 1e-2  # Multiplicative noise amplitude [dimensionless]\n",
    "sigma_tau = 1e-3  # Time base noise amplitude [ps]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1af9ac8d12a51a16",
   "metadata": {},
   "source": [
    "## Simulate an input waveform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60b9671fbdafa099",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = thz.timebase(n, dt=dt)\n",
    "mu = thz.wave(n, dt=dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24989e9f70b948f4",
   "metadata": {},
   "source": [
    "## Simulate measurements\n",
    "\n",
    "Define amplitude drift parameters `a` and delay drift parameters `eta`, each with length `m - 1`. Use the `numpy.repeat` function to generate an array of `m` copies of `mu`, then use `scaleshift` to rescale and shift all but the first of them by the elements `a` and `eta`, respectively. Use the transpose operation to orient the array so that the waveforms are arranged columnwise.\n",
    "\n",
    "Next, create an instance of the `NoiseModel` class and use the `noise_sim` method to add simulated noise to each waveform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa74c4d3e9cc730e",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(0)\n",
    "a = 1.0 + 1e-2 * rng.standard_normal(m - 1)\n",
    "eta = 1e-3 * rng.standard_normal(m - 1)\n",
    "\n",
    "z = thz.scaleshift(\n",
    "    np.repeat(np.atleast_2d(mu), m, axis=0), \n",
    "    dt=dt,\n",
    "    a=np.insert(a, 0, 1.0),\n",
    "    eta=np.insert(eta, 0, 0.0)\n",
    ").T\n",
    "\n",
    "noise_model = thz.NoiseModel(\n",
    "    sigma_alpha=sigma_alpha,\n",
    "    sigma_beta=sigma_beta,\n",
    "    sigma_tau=sigma_tau,\n",
    "    dt=dt\n",
    ")\n",
    "x = z + noise_model.noise_sim(z, axis=0, seed=12345)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef21995635b43d02",
   "metadata": {},
   "source": [
    "## Fit a noise model to the simulated measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cb9de4ebbc5d67c",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_res = thz.noisefit(\n",
    "    x,\n",
    "    sigma_alpha0=sigma_alpha,\n",
    "    sigma_beta0=sigma_beta,\n",
    "    sigma_tau0=sigma_tau,\n",
    "    dt=dt\n",
    ")\n",
    "print(f\"{noise_res.noise_model.sigma_alpha=}\")\n",
    "print(f\"{noise_res.noise_model.sigma_beta=}\")\n",
    "print(f\"{noise_res.noise_model.sigma_tau=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8fd1c067ecfd900",
   "metadata": {},
   "source": [
    "## Compare the fitted noise model with an empirical noise estimate\n",
    "\n",
    "Use the fitted values `noise_res.a` and `noise_res.eta` with `scaleshift` to correct for drift in `x`. Determine the standard deviation at each time point of the resulting `x_corrected` array. Compare this to the fitted noise model using the `noise_amp` method of `noise_res.noise_model` with `noise_res.mu`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52f86a563d2b496c",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_corrected = thz.scaleshift(\n",
    "    x, a=1 / noise_res.a, eta=-noise_res.eta, axis=0\n",
    ")\n",
    "plt.plot(t, np.std(x_corrected, axis=1), \"-\", label=\"Data\")\n",
    "plt.plot(t, noise_res.noise_model.noise_amp(noise_res.mu), \"--\", label=\"Fit\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"t (ps)\")\n",
    "plt.ylabel(r\"$\\sigma(t)$\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}