Merge pull request #163 from florisvb/common-kerneldiff

Made `kerneldiff` to address #162

Author: pavelkomarov

README.md

@@ -26,10 +26,10 @@ Python methods for numerical differentiation of noisy data, including multi-obje
 
 PyNumDiff is a Python package that implements various methods for computing numerical derivatives of noisy data, which can be a critical step in developing dynamic models or designing control. There are seven different families of methods implemented in this repository:
 
-1. convolutional smoothing followed by finite difference calculation
-2. polynomial fit methods
-3. basis function fit methods
-4. iterated finite differencing
+1. prefiltering followed by finite difference calculation
+2. iterated finite differencing
+3. polynomial fit methods
+4. basis function fit methods
 5. total variation regularization of a finite difference derivative
 6. generalized Kalman smoothing
 7. local approximation with linear model

docs/source/smooth_finite_difference.rst

@@ -2,4 +2,11 @@ smooth_finite_difference
 ========================
 
 .. automodule:: pynumdiff.smooth_finite_difference
-    :members:
+    :no-members:
+
+.. autofunction:: pynumdiff.smooth_finite_difference.kerneldiff
+.. autofunction:: pynumdiff.smooth_finite_difference.butterdiff
+.. autofunction:: pynumdiff.smooth_finite_difference.meandiff
+.. autofunction:: pynumdiff.smooth_finite_difference.mediandiff
+.. autofunction:: pynumdiff.smooth_finite_difference.gaussiandiff
+.. autofunction:: pynumdiff.smooth_finite_difference.friedrichsdiff

examples/1_basic_tutorial.ipynb

@@ -12,7 +12,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 3,
    "id": "5e5ea3c3-aed3-4b8d-a7c7-1db05fc73f3c",
    "metadata": {},
    "outputs": [],
@@ -31,7 +31,7 @@
     "from pynumdiff.utils.simulate import sine, triangle, pop_dyn, linear_autonomous, pi_cruise_control, lorenz_x\n",
     "from pynumdiff.utils.evaluate import rmse, error_correlation\n",
     "from pynumdiff.finite_difference import finitediff\n",
-    "from pynumdiff.smooth_finite_difference import mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff\n",
+    "from pynumdiff.smooth_finite_difference import kerneldiff, butterdiff\n",
     "from pynumdiff.polynomial_fit import splinediff, polydiff, savgoldiff\n",
     "from pynumdiff.basis_fit import spectraldiff, rbfdiff\n",
     "from pynumdiff.total_variation_regularization import tvrdiff, smooth_acceleration\n",
@@ -49,7 +49,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 4,
    "id": "ef461de3-e615-4c32-8931-51bf60d8e16c",
    "metadata": {},
    "outputs": [
@@ -74,26 +74,23 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 5,
    "id": "1eba51f1-bf31-4f84-b82b-176461319de6",
    "metadata": {},
    "outputs": [],
    "source": [
-    "methods = [(meandiff, 'MeanDiff'),\n",
-    "\t\t\t(mediandiff, 'MedianDiff'),\n",
-    "\t\t\t(gaussiandiff, 'GaussianDiff'),\n",
-    "\t\t\t(friedrichsdiff, 'FriedrichsDiff'),\n",
+    "methods = [(kerneldiff, 'KernelDiff'),\n",
     "\t\t\t(butterdiff, 'ButterDiff'),\n",
     "\t\t\t(splinediff, 'SplineDiff'),\n",
     "\t\t\t(polydiff, 'PolyDiff'),\n",
     "\t\t\t(savgoldiff, 'SavGolDiff'),\n",
-    "\t\t    (spectraldiff, 'SpectralDiff'),\n",
+    "\t\t\t(spectraldiff, 'SpectralDiff'),\n",
     "\t\t\t(rbfdiff, 'RBF'),\n",
     "\t\t\t(finitediff, 'IteratedFD'), # skip first_order, because it's not going to be the best\n",
     "\t\t\t(tvrdiff, 'TVR'),\n",
     "\t\t\t(smooth_acceleration, 'SmoothAccelTVR'), # skip in plotting?\n",
     "\t\t\t(rtsdiff, 'RTS'),\n",
-    "\t\t    (robustdiff, 'RobustDiff')]\n",
+    "\t\t\t(robustdiff, 'RobustDiff')]\n",
     "sims = [(pi_cruise_control, 'Cruise Control'),\n",
     "\t\t(sine, 'Sum of Sines'),\n",
     "\t\t(triangle, 'Triangles'),\n",
@@ -112,7 +109,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 6,
    "id": "275e4552-4479-40e3-b111-41c35d62ee9e",
    "metadata": {},
    "outputs": [],
@@ -121,14 +118,14 @@
     "\t\"\"\"Create and save a dataframe of results over all methods for all sims, given a certain\n",
     "\thyperparameterization. Each one of these takes about 11 minutes on my machine\"\"\"\t\n",
     "\ttry: # short circuit if the file already exists, just read in and display\n",
-    "\t\tres = pandas.read_csv(f\"~/Desktop/results2/{cutoff_frequency}_{dt}_{noise_type}_{noise_scale}_{outliers}_{random_seed}.csv\",\n",
+    "\t\tres = pandas.read_csv(f\"~/Desktop/results/{cutoff_frequency}_{dt}_{noise_type}_{noise_scale}_{outliers}_{random_seed}.csv\",\n",
     "\t\t\t\t\t\t\t index_col=0)\n",
     "\t\tprint(\".\", end='')\n",
     "\texcept FileNotFoundError:\n",
     "\t\ttvgamma = np.exp(-1.6*np.log(cutoff_frequency) -0.71*np.log(dt) - 5.1)\n",
-    "\t\t#res = pandas.DataFrame(index=[x[1] for x in methods], columns=[x[1] for x in sims])\n",
-    "\t\tres = pandas.read_csv(f\"~/Desktop/results/{cutoff_frequency}_{dt}_{noise_type}_{noise_scale}_{outliers}_{random_seed}.csv\",\n",
-    "\t\t\t\t\t\t\t index_col=0)\n",
+    "\t\tres = pandas.DataFrame(index=[x[1] for x in methods], columns=[x[1] for x in sims])\n",
+    "\t\t#res = pandas.read_csv(f\"~/Desktop/results/{cutoff_frequency}_{dt}_{noise_type}_{noise_scale}_{outliers}_{random_seed}.csv\",\n",
+    "\t\t#\t\t\t\t\t index_col=0)\n",
     "\n",
     "\t\tfor sim,sname in tqdm(sims):\n",
     "\t\t\tx, x_truth, dxdt_truth = sim(duration=4, dt=dt, noise_type=noise_type,\n",
@@ -137,7 +134,7 @@
     "\t\t\t#t = np.arange(0,dt*len(x), dt)\n",
     "\t\t\n",
     "\t\t\t#for method,mname in tqdm(methods):\n",
-    "\t\t\tmethod,mname = (robustdiff, 'RobustDiff') # ADDING A METHOD\n",
+    "\t\t\tmethod,mname = (robustdiff, 'KernelDiff') # ADDING A METHOD\n",
     "\t\t\t#best_params, best_rmse = optimize(method, x, dt, dxdt_truth=dxdt_truth, metric='rmse')\n",
     "\t\t\t#best_params, best_ec = optimize(method, x, dt, dxdt_truth=dxdt_truth, metric='error_correlation')\n",
     "\t\t\tbest_params, best_score = optimize(method, x, dt, tvgamma=tvgamma,\n",
@@ -148,7 +145,7 @@
     "\t\t\tec = error_correlation(dxdt_hat, dxdt_truth)\n",
     "\t\t\tres.loc[mname, sname] = f\"RMSE: {rms_dxdt:.5g}<br/>R^2: {ec:.5g}\"\n",
     "\t\n",
-    "\t\tres.to_csv(f\"~/Desktop/results2/{cutoff_frequency}_{dt}_{noise_type}_{noise_scale}_{outliers}_{random_seed}.csv\")\n",
+    "\t\tres.to_csv(f\"~/Desktop/results/{cutoff_frequency}_{dt}_{noise_type}_{noise_scale}_{outliers}_{random_seed}.csv\")\n",
     "\n",
     "\t\tdisplay(HTML(res.style.set_caption(\n",
     "\t\t\tf\"cutoff_frequency={cutoff_frequency}, dt={dt}, noise_type={noise_type}, noise_scale={noise_scale}, outliers={outliers}, random_seed={random_seed}\").set_properties(\n",
@@ -169,36 +166,43 @@
    "execution_count": null,
    "id": "fac1f6c9-f657-492b-8b68-9bdc28298261",
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "  0%|                                                                         | 0/6 [00:00<?, ?it/s]"
+     ]
+    }
+   ],
    "source": [
     "for random_seed in random_seeds:\n",
     "\t# Does dt matter? Vary dt, keeping everything else constant\n",
-    "\t# for dt in dts: # 4\n",
-    "\t# \tmain_loop(dt=dt, cutoff_frequency=3, noise_type='normal', noise_params=[0, 0.1], noise_scale=1,\n",
-    "\t# \t\t\t  outliers=False, random_seed=random_seed)\n",
+    "\tfor dt in dts: # 4\n",
+    "\t\tmain_loop(dt=dt, cutoff_frequency=3, noise_type='normal', noise_params=[0, 0.1], noise_scale=1,\n",
+    "\t\t\t\t  outliers=False, random_seed=random_seed)\n",
     "\n",
-    "\t# # Does noise type matter? Vary noise type, keeping everything else constant\n",
-    "\t# for noise_type,noise_params in noise_types: # 3\n",
-    "\t# \tmain_loop(dt=0.01, cutoff_frequency=3, noise_type=noise_type, noise_params=noise_params, noise_scale=1,\n",
-    "\t# \t\t\t  outliers=False, random_seed=random_seed)\n",
+    "\t# Does noise type matter? Vary noise type, keeping everything else constant\n",
+    "\tfor noise_type,noise_params in noise_types: # 3\n",
+    "\t\tmain_loop(dt=0.01, cutoff_frequency=3, noise_type=noise_type, noise_params=noise_params, noise_scale=1,\n",
+    "\t\t\t\t  outliers=False, random_seed=random_seed)\n",
     "\t\n",
-    "\t# # Does noise scale matter? Vary noise scale, keeping everything else constant\n",
-    "\t# for noise_scale in noise_scales: # 4\n",
-    "\t# \tmain_loop(dt=0.01, cutoff_frequency=3, noise_type='normal', noise_params=[0,0.1], noise_scale=noise_scale,\n",
-    "\t# \t\t\t  outliers=False, random_seed=random_seed)\n",
+    "\t# Does noise scale matter? Vary noise scale, keeping everything else constant\n",
+    "\tfor noise_scale in noise_scales: # 4\n",
+    "\t\tmain_loop(dt=0.01, cutoff_frequency=3, noise_type='normal', noise_params=[0,0.1], noise_scale=noise_scale,\n",
+    "\t\t\t\t  outliers=False, random_seed=random_seed)\n",
     "\t\n",
     "\t# Does the presence of outliers matter? Vary presence of outliers, keeping everything else constant\n",
     "\tfor outliers in [False, True]: # 2 Maybe try this one with more noise types too?\n",
     "\t\tmain_loop(dt=0.01, cutoff_frequency=3, noise_type='normal', noise_params=[0, 0.1], noise_scale=1,\n",
     "\t\t\t\t  outliers=outliers, random_seed=random_seed)\n",
     "\t\n",
-    "\t# #Does the cutoff frequency matter? Vary cutoff, keeping everything else constant\n",
-    "\t# for cutoff_frequency in cutoff_frequencies: # 5\n",
-    "\t# \tmain_loop(dt=0.01, cutoff_frequency=cutoff_frequency, noise_type='normal', noise_params=[0, 0.1], noise_scale=1,\n",
-    "\t# \t\t\t  outliers=False, random_seed=random_seed)\n",
+    "\t# Does the cutoff frequency matter? Vary cutoff, keeping everything else constant\n",
+    "\tfor cutoff_frequency in cutoff_frequencies: # 5\n",
+    "\t\tmain_loop(dt=0.01, cutoff_frequency=cutoff_frequency, noise_type='normal', noise_params=[0, 0.1], noise_scale=1,\n",
+    "\t\t\t\t  outliers=False, random_seed=random_seed)\n",
     "\t\n",
-    "\t# Because they all pass through the central point, the number of unique runs here is 4 + 3 + 2 + 1 + 2 = 12 * 11 minutes = 2 hours\n",
-    "\t# repeat for each random seed = 10 hours?"
+    "\t# Because they all pass through the central point, the number of unique runs here is 4 + 2 + 3 + 1 + 4 = 14"
    ]
   },
   {
@@ -224,9 +228,9 @@
     "\t\"\"\"\n",
     "\tfig, ax = plt.subplots(2, 1, figsize=(30, 12), sharex=True, gridspec_kw={'hspace': 0})\n",
     "\n",
-    "\tmarkers = ['^', 's', 'p', 'h', 'd', '2', '+', 'x', '*',  r'$\\rho$', 'o', r'$\\gamma$', r'$\\tilde{\\gamma}$', '.', r'$M$']\n",
-    "\tfacecolors = ['none' if i not in [5,6,7,13] else None for i in range(len(markers))] # 'none makes transparent'\n",
-    "\tsizes = [90, 70, 90, 90, 80, 120, 100, 70, 120, 90, 70, 90, 170, 70, 120]\n",
+    "\tmarkers = ['p', 'd', '*', '2', '+', 'x', 'o', '.', r'$\\gamma$', r'$\\tilde{\\gamma}$', '^', 's']\n",
+    "\tfacecolors = ['none' if i not in [3,4,5,7] else None for i in range(len(markers))] # 'none makes transparent'\n",
+    "\tsizes = [90, 80, 120, 120, 100, 70, 70, 70, 90, 170, 90, 70]\n",
     "\tcmap = plt.get_cmap('turbo', 6) # Assign a unique color for each simulation\n",
     "\tcolors = [cmap(i) for i in range(6)]; colors[0] = 'purple'; colors[-1] = 'red'\n",
     "\tcolors[2] = [x*0.8 for x in colors[2]]; colors[3] = mcolors.to_rgb('gold'); colors[3] = [x*0.9 for x in colors[3]]\n",

examples/2a_optimizing_parameters_with_dxdt_known.ipynb

@@ -13,7 +13,7 @@
     from .linear_model import lineardiff
 
 from .finite_difference import finitediff, first_order, second_order, fourth_order
-from .smooth_finite_difference import meandiff, mediandiff, gaussiandiff, friedrichsdiff, butterdiff
+from .smooth_finite_difference import kerneldiff, meandiff, mediandiff, gaussiandiff, friedrichsdiff, butterdiff
 from .polynomial_fit import splinediff, polydiff, savgoldiff
 from .basis_fit import spectraldiff, rbfdiff
 from .total_variation_regularization import iterative_velocity

examples/2b_optimizing_parameters_with_dxdt_unknown.ipynb

@@ -65,7 +65,8 @@ def finitediff(x, dt, num_iterations, order):
 
 
 def first_order(x, dt, params=None, options={}, num_iterations=1):
-    """First-order difference method
+    """First-order difference method\n
+    **Deprecated**, prefer :code:`finitediff` with order 1 instead.
 
     :param np.array[float] x: data to differentiate
     :param float dt: step size
@@ -85,11 +86,13 @@ def first_order(x, dt, params=None, options={}, num_iterations=1):
         warn("`params` and `options` parameters will be removed in a future version. Use `num_iterations` instead.", DeprecationWarning)
         num_iterations = params[0] if isinstance(params, list) else params
 
+    warn("`first_order` is deprecated. Call `finitediff` with order 1 instead.", DeprecationWarning)
     return finitediff(x, dt, num_iterations, 1)
 
 
 def second_order(x, dt, num_iterations=1):
-    """Second-order centered difference method, with special endpoint formulas.
+    """Second-order centered difference method, with special endpoint formulas.\n
+    **Deprecated**, prefer :code:`finitediff` with order 2 instead.
 
     :param np.array[float] x: data to differentiate
     :param float dt: step size
@@ -100,11 +103,13 @@ def second_order(x, dt, num_iterations=1):
              - **x_hat** -- original x if :code:`num_iterations=1`, else smoothed x that yielded dxdt_hat
              - **dxdt_hat** -- estimated derivative of x
     """
+    warn("`second_order` is deprecated. Call `finitediff` with order 2 instead.", DeprecationWarning)
     return finitediff(x, dt, num_iterations, 2)
 
 
 def fourth_order(x, dt, num_iterations=1):
-    """Fourth-order centered difference method, with special endpoint formulas.
+    """Fourth-order centered difference method, with special endpoint formulas.\n
+    **Deprecated**, prefer :code:`finitediff` with order 4 instead.
 
     :param np.array[float] x: data to differentiate
     :param float dt: step size
@@ -115,4 +120,5 @@ def fourth_order(x, dt, num_iterations=1):
              - **x_hat** -- original x if :code:`num_iterations=1`, else smoothed x that yielded dxdt_hat
              - **dxdt_hat** -- estimated derivative of x
     """
+    warn("`fourth_order` is deprecated. Call `finitediff` with order 4 instead.", DeprecationWarning)
     return finitediff(x, dt, num_iterations, 4)

examples/4_performance_analysis.ipynb

@@ -170,7 +170,8 @@ def rtsdiff(x, _t, order, qr_ratio, forwardbackward):
 
 
 def constant_velocity(x, dt, params=None, options=None, r=None, q=None, forwardbackward=True):
-    """Run a forward-backward constant velocity RTS Kalman smoother to estimate the derivative.
+    """Run a forward-backward constant velocity RTS Kalman smoother to estimate the derivative.\n
+    **Deprecated**, prefer :code:`rtsdiff` with order 1 instead.
 
     :param np.array[float] x: data series to differentiate
     :param float dt: step size
@@ -195,11 +196,13 @@ def constant_velocity(x, dt, params=None, options=None, r=None, q=None, forwardb
     elif r == None or q == None:
         raise ValueError("`q` and `r` must be given.")
 
+    warn("`constant_velocity` is deprecated. Call `rtsdiff` with order 1 instead.", DeprecationWarning)
     return rtsdiff(x, dt, 1, q/r, forwardbackward)
 
 
 def constant_acceleration(x, dt, params=None, options=None, r=None, q=None, forwardbackward=True):
-    """Run a forward-backward constant acceleration RTS Kalman smoother to estimate the derivative.
+    """Run a forward-backward constant acceleration RTS Kalman smoother to estimate the derivative.\n
+    **Deprecated**, prefer :code:`rtsdiff` with order 2 instead.
 
     :param np.array[float] x: data series to differentiate
     :param float dt: step size
@@ -224,11 +227,13 @@ def constant_acceleration(x, dt, params=None, options=None, r=None, q=None, forw
     elif r == None or q == None:
         raise ValueError("`q` and `r` must be given.")
 
+    warn("`constant_acceleration` is deprecated. Call `rtsdiff` with order 2 instead.", DeprecationWarning)
     return rtsdiff(x, dt, 2, q/r, forwardbackward)
 
 
 def constant_jerk(x, dt, params=None, options=None, r=None, q=None, forwardbackward=True):
-    """Run a forward-backward constant jerk RTS Kalman smoother to estimate the derivative.
+    """Run a forward-backward constant jerk RTS Kalman smoother to estimate the derivative.\n
+    **Deprecated**, prefer :code:`rtsdiff` with order 3 instead.
 
     :param np.array[float] x: data series to differentiate
     :param float dt: step size
@@ -253,6 +258,7 @@ def constant_jerk(x, dt, params=None, options=None, r=None, q=None, forwardbackw
     elif r == None or q == None:
         raise ValueError("`q` and `r` must be given.")
 
+    warn("`constant_jerk` is deprecated. Call `rtsdiff` with order 3 instead.", DeprecationWarning)
     return rtsdiff(x, dt, 3, q/r, forwardbackward)