numerical gradient for arbitrary functions that have codomain as the set of all real numbers

gradient.py

#! /usr/bin/env python3

from copy import deepcopy


def median(l):
    l = sorted(l)

    if len(l) % 2 == 0:
        return (l[len(l)//2] + l[len(l)//2 - 1])/2
    else:
        return l[len(l)//2]


def gradient(function, input, epsilon=0.00001):
    grad = [0.0] * len(input)

    for ix in range(len(input)):
        icopy = deepcopy(input)
        icopy[ix] += epsilon

        fx_plus_h = function(icopy)

        icopy = deepcopy(input)
        icopy[ix] -= epsilon

        fx_minus_h = function(icopy)

        grad[ix] = (fx_plus_h - fx_minus_h) / (2 * epsilon)

    return grad


if __name__ == '__main__':
    l = [11, 5, 23, 4, 56, 32]
    g = gradient(median, l)

    print(g)