text: import numpy as np
  from numpy.testing import (assert_allclose, assert_equal,
                             assert_almost_equal, assert_array_equal,
                             assert_array_almost_equal)
  
  from scipy.ndimage import convolve1d
  
  from scipy.signal import savgol_coeffs, savgol_filter
  from scipy.signal._savitzky_golay import _polyder
  
  
  def check_polyder(p, m, expected):
      dp = _polyder(p, m)
      assert_array_equal(dp, expected)
  
  
  def test_polyder():
      cases = [
          ([5], 0, [5]),
          ([5], 1, [0]),
          ([3, 2, 1], 0, [3, 2, 1]),
          ([3, 2, 1], 1, [6, 2]),
          ([3, 2, 1], 2, [6]),
          ([3, 2, 1], 3, [0]),
          ([[3, 2, 1], [5, 6, 7]], 0, [[3, 2, 1], [5, 6, 7]]),
          ([[3, 2, 1], [5, 6, 7]], 1, [[6, 2], [10, 6]]),
          ([[3, 2, 1], [5, 6, 7]], 2, [[6], [10]]),
          ([[3, 2, 1], [5, 6, 7]], 3, [[0], [0]]),
      ]
      for p, m, expected in cases:
          check_polyder(np.array(p).T, m, np.array(expected).T)
  
  
  #--------------------------------------------------------------------
  # savgol_coeffs tests
  #--------------------------------------------------------------------
  
  def alt_sg_coeffs(window_length, polyorder, pos):
      """This is an alternative implementation of the SG coefficients.
  
      It uses numpy.polyfit and numpy.polyval. The results should be
      equivalent to those of savgol_coeffs(), but this implementation
      is slower.
  
      window_length should be odd.
  
      """
      if pos is None:
          pos = window_length // 2
      t = np.arange(window_length)
      unit = (t == pos).astype(int)
      h = np.polyval(np.polyfit(t, unit, polyorder), t)
      return h
  
  
  def test_sg_coeffs_trivial():
      # Test a trivial case of savgol_coeffs: polyorder = window_length - 1
      h = savgol_coeffs(1, 0)
      assert_allclose(h, [1])
  
      h = savgol_coeffs(3, 2)
      assert_allclose(h, [0, 1, 0], atol=1e-10)
  
      h = savgol_coeffs(5, 4)
      assert_allclose(h, [0, 0, 1, 0, 0], atol=1e-10)
  
      h = savgol_coeffs(5, 4, pos=1)
      assert_allclose(h, [0, 0, 0, 1, 0], atol=1e-10)
  
      h = savgol_coeffs(5, 4, pos=1, use='dot')
      assert_allclose(h, [0, 1, 0, 0, 0], atol=1e-10)
  
  
  def compare_coeffs_to_alt(window_length, order):
      # For the given window_length and order, compare the results
      # of savgol_coeffs and alt_sg_coeffs for pos from 0 to window_length - 1.
      # Also include pos=None.
      for pos in [None] + list(range(window_length)):
          h1 = savgol_coeffs(window_length, order, pos=pos, use='dot')
          h2 = alt_sg_coeffs(window_length, order, pos=pos)
          assert_allclose(h1, h2, atol=1e-10,
                          err_msg=("window_length = %d, order = %d, pos = %s" %
                                   (window_length, order, pos)))
  
  
  def test_sg_coeffs_compare():
      # Compare savgol_coeffs() to alt_sg_coeffs().
      for window_length in range(1, 8, 2):
          for order in range(window_length):
              compare_coeffs_to_alt(window_length, order)
  
  
  def test_sg_coeffs_exact():
      polyorder = 4
      window_length = 9
      halflen = window_length // 2
  
      x = np.linspace(0, 21, 43)
      delta = x[1] - x[0]
  
      # The data is a cubic polynomial.  We'll use an order 4
      # SG filter, so the filtered values should equal the input data
      # (except within half window_length of the edges).
      y = 0.5 * x ** 3 - x
      h = savgol_coeffs(window_length, polyorder)
      y0 = convolve1d(y, h)
      assert_allclose(y0[halflen:-halflen], y[halflen:-halflen])
  
      # Check the same input, but use deriv=1.  dy is the exact result.
      dy = 1.5 * x ** 2 - 1
      h = savgol_coeffs(window_length, polyorder, deriv=1, delta=delta)
      y1 = convolve1d(y, h)
      assert_allclose(y1[halflen:-halflen], dy[halflen:-halflen])
  
      # Check the same input, but use deriv=2. d2y is the exact result.
      d2y = 3.0 * x
      h = savgol_coeffs(window_length, polyorder, deriv=2, delta=delta)
      y2 = convolve1d(y, h)
      assert_allclose(y2[halflen:-halflen], d2y[halflen:-halflen])
  
  
  def test_sg_coeffs_deriv():
      # The data in `x` is a sampled parabola, so using savgol_coeffs with an
      # order 2 or higher polynomial should give exact results.
      i = np.array([-2.0, 0.0, 2.0, 4.0, 6.0])
      x = i ** 2 / 4
      dx = i / 2
      d2x = np.full_like(i, 0.5)
      for pos in range(x.size):
          coeffs0 = savgol_coeffs(5, 3, pos=pos, delta=2.0, use='dot')
          assert_allclose(coeffs0.dot(x), x[pos], atol=1e-10)
          coeffs1 = savgol_coeffs(5, 3, pos=pos, delta=2.0, use='dot', deriv=1)
          assert_allclose(coeffs1.dot(x), dx[pos], atol=1e-10)
          coeffs2 = savgol_coeffs(5, 3, pos=pos, delta=2.0, use='dot', deriv=2)
          assert_allclose(coeffs2.dot(x), d2x[pos], atol=1e-10)
  
  
  def test_sg_coeffs_deriv_gt_polyorder():
      """
      If deriv > polyorder, the coefficients should be all 0.
      This is a regression test for a bug where, e.g.,
          savgol_coeffs(5, polyorder=1, deriv=2)
      raised an error.
      """
      coeffs = savgol_coeffs(5, polyorder=1, deriv=2)
      assert_array_equal(coeffs, np.zeros(5))
      coeffs = savgol_coeffs(7, polyorder=4, deriv=6)
      assert_array_equal(coeffs, np.zeros(7))
  
  
  def test_sg_coeffs_large():
      # Test that for large values of window_length and polyorder the array of
      # coefficients returned is symmetric. The aim is to ensure that
      # no potential numeric overflow occurs.
      coeffs0 = savgol_coeffs(31, 9)
      assert_array_almost_equal(coeffs0, coeffs0[::-1])
      coeffs1 = savgol_coeffs(31, 9, deriv=1)
      assert_array_almost_equal(coeffs1, -coeffs1[::-1])
  
  
  #--------------------------------------------------------------------
  # savgol_filter tests
  #--------------------------------------------------------------------
  
  
  def test_sg_filter_trivial():
      """ Test some trivial edge cases for savgol_filter()."""
      x = np.array([1.0])
      y = savgol_filter(x, 1, 0)
      assert_equal(y, [1.0])
  
      # Input is a single value. With a window length of 3 and polyorder 1,
      # the value in y is from the straight-line fit of (-1,0), (0,3) and
      # (1, 0) at 0. This is just the average of the three values, hence 1.0.
      x = np.array([3.0])
      y = savgol_filter(x, 3, 1, mode='constant')
      assert_almost_equal(y, [1.0], decimal=15)
  
      x = np.array([3.0])
      y = savgol_filter(x, 3, 1, mode='nearest')
      assert_almost_equal(y, [3.0], decimal=15)
  
      x = np.array([1.0] * 3)
      y = savgol_filter(x, 3, 1, mode='wrap')
      assert_almost_equal(y, [1.0, 1.0, 1.0], decimal=15)
  
  
  def test_sg_filter_basic():
      # Some basic test cases for savgol_filter().
      x = np.array([1.0, 2.0, 1.0])
      y = savgol_filter(x, 3, 1, mode='constant')
      assert_allclose(y, [1.0, 4.0 / 3, 1.0])
  
      y = savgol_filter(x, 3, 1, mode='mirror')
      assert_allclose(y, [5.0 / 3, 4.0 / 3, 5.0 / 3])
  
      y = savgol_filter(x, 3, 1, mode='wrap')
      assert_allclose(y, [4.0 / 3, 4.0 / 3, 4.0 / 3])
  
  
  def test_sg_filter_2d():
      x = np.array([[1.0, 2.0, 1.0],
                    [2.0, 4.0, 2.0]])
      expected = np.array([[1.0, 4.0 / 3, 1.0],
                           [2.0, 8.0 / 3, 2.0]])
      y = savgol_filter(x, 3, 1, mode='constant')
      assert_allclose(y, expected)
  
      y = savgol_filter(x.T, 3, 1, mode='constant', axis=0)
      assert_allclose(y, expected.T)
  
  
  def test_sg_filter_interp_edges():
      # Another test with low degree polynomial data, for which we can easily
      # give the exact results. In this test, we use mode='interp', so
      # savgol_filter should match the exact solution for the entire data set,
      # including the edges.
      t = np.linspace(-5, 5, 21)
      delta = t[1] - t[0]
      # Polynomial test data.
      x = np.array([t,
                    3 * t ** 2,
                    t ** 3 - t])
      dx = np.array([np.ones_like(t),
                     6 * t,
                     3 * t ** 2 - 1.0])
      d2x = np.array([np.zeros_like(t),
                      np.full_like(t, 6),
                      6 * t])
  
      window_length = 7
  
      y = savgol_filter(x, window_length, 3, axis=-1, mode='interp')
      assert_allclose(y, x, atol=1e-12)
  
      y1 = savgol_filter(x, window_length, 3, axis=-1, mode='interp',
                         deriv=1, delta=delta)
      assert_allclose(y1, dx, atol=1e-12)
  
      y2 = savgol_filter(x, window_length, 3, axis=-1, mode='interp',
                         deriv=2, delta=delta)
      assert_allclose(y2, d2x, atol=1e-12)
  
      # Transpose everything, and test again with axis=0.
  
      x = x.T
      dx = dx.T
      d2x = d2x.T
  
      y = savgol_filter(x, window_length, 3, axis=0, mode='interp')
      assert_allclose(y, x, atol=1e-12)
  
      y1 = savgol_filter(x, window_length, 3, axis=0, mode='interp',
                         deriv=1, delta=delta)
      assert_allclose(y1, dx, atol=1e-12)
  
      y2 = savgol_filter(x, window_length, 3, axis=0, mode='interp',
                         deriv=2, delta=delta)
      assert_allclose(y2, d2x, atol=1e-12)
  
  
  def test_sg_filter_interp_edges_3d():
      # Test mode='interp' with a 3-D array.
      t = np.linspace(-5, 5, 21)
      delta = t[1] - t[0]
      x1 = np.array([t, -t])
      x2 = np.array([t ** 2, 3 * t ** 2 + 5])
      x3 = np.array([t ** 3, 2 * t ** 3 + t ** 2 - 0.5 * t])
      dx1 = np.array([np.ones_like(t), -np.ones_like(t)])
      dx2 = np.array([2 * t, 6 * t])
      dx3 = np.array([3 * t ** 2, 6 * t ** 2 + 2 * t - 0.5])
  
      # z has shape (3, 2, 21)
      z = np.array([x1, x2, x3])
      dz = np.array([dx1, dx2, dx3])
  
      y = savgol_filter(z, 7, 3, axis=-1, mode='interp', delta=delta)
      assert_allclose(y, z, atol=1e-10)
  
      dy = savgol_filter(z, 7, 3, axis=-1, mode='interp', deriv=1, delta=delta)
      assert_allclose(dy, dz, atol=1e-10)
  
      # z has shape (3, 21, 2)
      z = np.array([x1.T, x2.T, x3.T])
      dz = np.array([dx1.T, dx2.T, dx3.T])
  
      y = savgol_filter(z, 7, 3, axis=1, mode='interp', delta=delta)
      assert_allclose(y, z, atol=1e-10)
  
      dy = savgol_filter(z, 7, 3, axis=1, mode='interp', deriv=1, delta=delta)
      assert_allclose(dy, dz, atol=1e-10)
  
      # z has shape (21, 3, 2)
      z = z.swapaxes(0, 1).copy()
      dz = dz.swapaxes(0, 1).copy()
  
      y = savgol_filter(z, 7, 3, axis=0, mode='interp', delta=delta)
      assert_allclose(y, z, atol=1e-10)
  
      dy = savgol_filter(z, 7, 3, axis=0, mode='interp', deriv=1, delta=delta)
      assert_allclose(dy, dz, atol=1e-10)
  