kernel nonsense

def _gaussian_kernel1d(sigma, order, radius):
	'''
	Produces a 1D Gaussian or Gaussian-derivative filter kernel as a numpy array.

	Args:
		sigma (float): The standard deviation of the filter.
		order (int): The derivative-order to use. 0 indicates a Gaussian function, 1 a 1st order derivative, etc.
		radius (int): The kernel produced will be of length (2*radius+1)

	Returns:
		Array of length (2*radius+1) containing the filter kernel.
	'''
	if order < 0:
		raise ValueError('order must non-negative')
	if not (isinstance(radius, int) or radius.is_integer()) or radius <= 0:
		raise ValueError('radius must a positive integer')

	p = np.polynomial.Polynomial([0, 0, -0.5 / (sigma * sigma)])
	x = np.arange(-radius, radius + 1)
	phi_x = np.exp(p(x), dtype=np.double)
	phi_x /= phi_x.sum()

	if order > 0:
		# For Gaussian-derivative filters, the function must be derived one or more times.
		q = np.polynomial.Polynomial([1])
		p_deriv = p.deriv()
		for _ in range(order):
			# f(x) = q(x) * phi(x) = q(x) * exp(p(x))
			# f'(x) = (q'(x) + q(x) * p'(x)) * phi(x)
			q = q.deriv() + q * p_deriv
		phi_x *= q(x)

	return phi_x

def smooth(x, sigma):
	'''
	Smooths a 1D array via a Gaussian filter.
	
	Args:
		x (array of floats): The array to be smoothed.
		sigma (float): The standard deviation of the smoothing filter to use.

	Returns:
		Array of same length as x.
	'''
	
	if len(x) == 0:
		raise ValueError('Cannot smooth an empty array')
		
	# Choose a radius for the filter kernel large enough to include all significant elements. Using
	# a radius of 4 standard deviations (rounded to int) will only truncate tail values that are of
	# the order of 1e-5 or smaller. For very small sigma values, just use a minimal radius.
	kernel_radius = max(1, int(round(4.0 * sigma)))
	filter_kernel = _gaussian_kernel1d(sigma, 0, kernel_radius)
	
	# The filter kernel will be applied by convolution in 'valid' mode, which includes only the
	# parts of the convolution in which the two signals full overlap, i.e. where the shorter signal
	# is entirely contained within the longer signal, producing an output signal of length equal to
	# the difference in length between the two input signals, plus one. So the input signal must be
	# extended by mirroring the ends (to give realistic values for the first and last pixels after
	# smoothing) until len(x_mirrored) - len(w) + 1 = len(x). This requires adding (len(w)-1)/2
	# values to each end of the input. If len(x) < (len(w)-1)/2, then the mirroring will need to be
	# performed over multiple iterations, as the mirrors can only, at most, triple the length of x
	# each time they are applied.
	extended_input_len = len(x) + len(filter_kernel) - 1
	x_mirrored = x
	while len(x_mirrored) < extended_input_len:
		mirror_len = min(len(x_mirrored), (extended_input_len - len(x_mirrored))//2)
		x_mirrored = np.r_[x_mirrored[mirror_len-1 : : -1],
						   x_mirrored,
						   x_mirrored[-1 : -(mirror_len+1) : -1]]
						   
	# Convolve the extended input copy with the filter kernel to apply the filter.
	# Convolving in 'valid' mode clips includes only the parts of the convolution in which the two
	# signals full overlap, i.e. the shorter singal is entirely contained within the longer signal.
	# It produces an output of length equal to the difference in length between the two input
	# signals, plus one. So this relies on the assumption that len(s) - len(w) + 1 >= len(x).
	y = np.convolve(x_mirrored, filter_kernel, mode="valid")
	
	assert(len(y) == len(x))
	
	return y