{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-25T22:56:28.405373Z",
     "start_time": "2020-05-25T22:56:27.866309Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sys import path\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['lines.markersize'] = 6\n",
    "mpl.rcParams['scatter.marker'] = '.'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-25T22:56:28.671566Z",
     "start_time": "2020-05-25T22:56:28.447539Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "path.append('../../')\n",
    "import derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-25T22:58:54.742421Z",
     "start_time": "2020-05-25T22:58:54.735182Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_example(diff_method, t, data_f, res_f, sigmas, y_label=None):\n",
    "    '''Utility function for concise plotting of examples.'''\n",
    "    fig, axes = plt.subplots(1, len(sigmas), figsize=[len(sigmas)*4, 3])\n",
    "    \n",
    "    # Compute the derivative\n",
    "    res = diff_method.d(np.vstack([data_f(t, s) for s in sigmas]), t, axis=1)\n",
    "    for i, s in enumerate(sigmas):\n",
    "        axes[i].plot(t, res_f(t))\n",
    "        axes[i].plot(t, res[i])\n",
    "        axes[i].set_title(r\"Noise: $\\sigma$={}\".format(s))\n",
    "        axes[i].set_ylim([-1.25,1.3])\n",
    "    if y_label:\n",
    "        axes[0].set_ylabel(y_label, fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two ways to interact with the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-25T22:56:42.981130Z",
     "start_time": "2020-05-25T22:56:42.977150Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "t = np.linspace(0, 2, 50)\n",
    "x = np.sin(2*np.pi*t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first way is to do a specific import of the desired Derivative object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-25T22:56:43.865597Z",
     "start_time": "2020-05-25T22:56:43.667400Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from derivative import FiniteDifference\n",
    "\n",
    "fig,ax = plt.subplots(1, figsize=[5,3])\n",
    "kind = FiniteDifference(k=1)\n",
    "ax.plot(t, kind.d(x,t));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second way is top use the functional interface and pass the kind of derivative as an argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-25T22:56:46.128870Z",
     "start_time": "2020-05-25T22:56:45.934615Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Use the functional interface and pass the kind as an argument\n",
    "from derivative import dxdt\n",
    "\n",
    "fig,ax = plt.subplots(1, figsize=[5,3])\n",
    "ax.plot(t, dxdt(x, t, \"finite_difference\", k=1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Smooth Derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first example is a sine function with Gaussian noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-25T22:56:55.127034Z",
     "start_time": "2020-05-25T22:56:54.538403Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def noisy_sin(t, sigma):\n",
    "    '''Sine with gaussian noise.'''\n",
    "    np.random.seed(17)\n",
    "    return np.sin(t) + np.random.normal(loc=0, scale=sigma, size=t.shape)\n",
    "\n",
    "sigmas = [0, 0.01, 0.1]\n",
    "fig, ax = plt.subplots(1, len(sigmas), figsize=[len(sigmas)*4, 3])\n",
    "\n",
    "t = np.linspace(0, 2*np.pi, 50, endpoint=False)\n",
    "for axs, s in zip(ax, sigmas): \n",
    "    axs.scatter(t, noisy_sin(t, s))\n",
    "    axs.set_title(r\"Noise: $\\sigma$={}\".format(s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finite differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T23:08:12.474754Z",
     "start_time": "2020-03-25T23:08:11.970077Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fd = derivative.FiniteDifference(3)\n",
    "plot_example(fd, t, noisy_sin, np.cos, sigmas, 'k: 3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Savitzky-Golay filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T23:08:15.637803Z",
     "start_time": "2020-03-25T23:08:14.156468Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sg = derivative.SavitzkyGolay(left=.5, right=.5, order=2, periodic=False)\n",
    "plot_example(sg, t, noisy_sin, np.cos, sigmas, 'Periodic: False, window: 1')\n",
    "\n",
    "sg = derivative.SavitzkyGolay(left=.5, right=.5, order=2, periodic=True)\n",
    "plot_example(sg, t, noisy_sin, np.cos, sigmas, 'Periodic: True, window: 1')\n",
    "\n",
    "sg = derivative.SavitzkyGolay(left=2, right=2, order=3, periodic=False)\n",
    "plot_example(sg, t, noisy_sin, np.cos, sigmas, 'Periodic: False, window: 4')\n",
    "\n",
    "sg = derivative.SavitzkyGolay(left=2, right=2, order=3, periodic=True)\n",
    "plot_example(sg, t, noisy_sin, np.cos, sigmas, 'Periodic: True, window: 4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Splines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-25T22:57:29.126819Z",
     "start_time": "2020-05-25T22:57:28.186814Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "spl = derivative.Spline(.5)\n",
    "plot_example(spl, t, noisy_sin, np.cos, sigmas, 's: 0.5, periodic: False')\n",
    "spl = derivative.Spline(.5, periodic=True)\n",
    "plot_example(spl, t, noisy_sin, np.cos, sigmas, 's: 0.5, periodic: True')\n",
    "spl = derivative.Spline(1, periodic=True)\n",
    "plot_example(spl, t, noisy_sin, np.cos, sigmas, 's: 1, periodic: True')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral method - Fourier basis\n",
    "Add your own filter!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T23:08:22.706408Z",
     "start_time": "2020-03-25T23:08:21.998725Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "no_filter =  derivative.Spectral()\n",
    "yes_filter = derivative.Spectral(filter=np.vectorize(lambda k: 1 if abs(k) < 3 else 0))\n",
    "\n",
    "plot_example(no_filter, t, noisy_sin, np.cos, sigmas, 'No filter')\n",
    "plot_example(yes_filter, t, noisy_sin, np.cos, sigmas, 'Low-pass filter')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral method - Chebyshev basis\n",
    "\n",
    "Now let's do with the Chebyshev basis, which requires cosine-spaced points on [a, b] rather than equispaced points on [a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_cos = np.cos(np.pi * np.arange(50) / 49) * np.pi + np.pi # choose a = 0, b = 2*pi\n",
    "no_filter = derivative.Spectral(basis='chebyshev')\n",
    "yes_filter = derivative.Spectral(basis='chebyshev', filter=np.vectorize(lambda k: 1 if abs(k) < 6 else 0))\n",
    "\n",
    "plot_example(no_filter, t_cos, noisy_sin, np.cos, sigmas, 'No filter')\n",
    "plot_example(yes_filter, t_cos, noisy_sin, np.cos, sigmas, 'Low-pass filter')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trend-filtered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tvd =  derivative.TrendFiltered(alpha=1e-3, order=0, max_iter=int(1e6))\n",
    "plot_example(tvd, t, noisy_sin, np.cos, sigmas, 'order: 0')\n",
    "\n",
    "tvd =  derivative.TrendFiltered(alpha=1e-3, order=1, max_iter=int(1e6))\n",
    "plot_example(tvd, t, noisy_sin, np.cos, sigmas, 'order: 1')\n",
    "\n",
    "tvd =  derivative.TrendFiltered(alpha=1e-3, order=2, max_iter=int(1e6))\n",
    "plot_example(tvd, t, noisy_sin, np.cos, sigmas, 'order: 2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kalman smoothing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "kal =  derivative.Kalman(alpha=0.01)\n",
    "plot_example(kal, t, noisy_sin, np.cos, sigmas, 'alpha: 0.01')\n",
    "\n",
    "kal =  derivative.Kalman(alpha=0.5)\n",
    "plot_example(kal, t, noisy_sin, np.cos, sigmas, 'alpha: 0.1')\n",
    "\n",
    "kal =  derivative.Kalman(alpha=1)\n",
    "plot_example(kal, t, noisy_sin, np.cos, sigmas, 'alpha: 1.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jump Derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second example is the absolute value function with Gaussian noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T23:08:44.039665Z",
     "start_time": "2020-03-25T23:08:44.032588Z"
    }
   },
   "outputs": [],
   "source": [
    "def noisy_abs(t, sigma):\n",
    "    '''Abs with gaussian noise.'''\n",
    "    np.random.seed(17)\n",
    "    return np.abs(t) + np.random.normal(loc=0, scale=sigma, size=x.shape)\n",
    "\n",
    "d_abs = lambda t: t/abs(t)\n",
    "\n",
    "sigmas = [0, 0.01, 0.1]\n",
    "fig, ax = plt.subplots(1, len(sigmas), figsize=[len(sigmas)*4, 3])\n",
    "\n",
    "t = np.linspace(-1, 1, 50)\n",
    "for axs, s in zip(ax, sigmas): \n",
    "    axs.scatter(t, noisy_abs(t, s))\n",
    "    axs.set_title(r\"Noise: $\\sigma$={}\".format(s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finite differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T23:08:46.763052Z",
     "start_time": "2020-03-25T23:08:46.289060Z"
    }
   },
   "outputs": [],
   "source": [
    "fd = derivative.FiniteDifference(k=3)\n",
    "plot_example(fd, t, noisy_abs, d_abs, sigmas, 'k: 3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Savitzky-Galoy filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T23:08:49.122429Z",
     "start_time": "2020-03-25T23:08:47.948831Z"
    }
   },
   "outputs": [],
   "source": [
    "sg = derivative.SavitzkyGolay(left=.125, right=.125, order=2, periodic=True, T=2)\n",
    "plot_example(sg, t, noisy_abs, d_abs, sigmas, 'window size: .25')\n",
    "\n",
    "sg = derivative.SavitzkyGolay(left=.25, right=.25, order=3, periodic=True, T=2)\n",
    "plot_example(sg, t, noisy_abs, d_abs, sigmas, 'window size: .5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Splines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T21:30:51.279530Z",
     "start_time": "2020-03-25T21:30:50.704947Z"
    }
   },
   "outputs": [],
   "source": [
    "spl = derivative.Spline(.1, periodic=False)\n",
    "plot_example(spl, t, noisy_abs, d_abs, sigmas, 's: 0.1, periodic: False')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral method - Fourier basis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T21:30:52.281590Z",
     "start_time": "2020-03-25T21:30:51.283796Z"
    }
   },
   "outputs": [],
   "source": [
    "no_filter =  derivative.Spectral()\n",
    "yes_filter = derivative.Spectral(filter=np.vectorize(lambda k: 1 if abs(k) < 6 else 0))\n",
    "\n",
    "plot_example(no_filter, t, noisy_abs, d_abs, sigmas, 'No filter')\n",
    "plot_example(yes_filter, t, noisy_abs, d_abs, sigmas, 'Low-pass filter')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral method - Chebyshev basis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_cos = np.cos(np.pi * np.arange(50)/49)\n",
    "no_filter = derivative.Spectral(basis='chebyshev')\n",
    "yes_filter = derivative.Spectral(basis='chebyshev', filter=np.vectorize(lambda k: 1 if abs(k) < 15 else 0))\n",
    "\n",
    "plot_example(no_filter, t_cos, noisy_abs, d_abs, sigmas, 'No filter')\n",
    "plot_example(yes_filter, t_cos, noisy_abs, d_abs, sigmas, 'Low-pass filter')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trend-filtered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-25T21:30:52.966158Z",
     "start_time": "2020-03-25T21:30:52.284317Z"
    }
   },
   "outputs": [],
   "source": [
    "tvd =  derivative.TrendFiltered(alpha=1e-3, order=0, max_iter=int(1e5))\n",
    "plot_example(tvd, t, noisy_abs, d_abs, sigmas, 'order: 0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kalman smoothing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "kal =  derivative.Kalman(alpha=0.01)\n",
    "plot_example(kal, t, noisy_abs, d_abs, sigmas, 'alpha: 0.01')\n",
    "\n",
    "kal =  derivative.Kalman(alpha=0.1)\n",
    "plot_example(kal, t, noisy_abs, d_abs, sigmas, 'alpha: 0.1')\n",
    "\n",
    "kal =  derivative.Kalman(alpha=1.)\n",
    "plot_example(kal, t, noisy_abs, d_abs, sigmas, 'alpha: 1.')"
   ]
  }
 ],
 "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.13.1"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "319.809px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}