Computes f(x) given that f(x_sofar) = x_sofar.
    returns f(x), and new versions of x_sofar and f_sofar.

def caching_call(f, x, x_sofar, f_sofar):
    """
    Computes f(x) given that f(x_sofar) = x_sofar.
    returns f(x), and new versions of x_sofar and f_sofar.
    """
    lenx = x.shape[0]

    nr, rf, rt, nu, xu, ui = remove_duplicates(x)

    unique_indices = ui[:nu]
    x_unique = xu[:nu]
    repeat_from = rf[:nr]
    repeat_to = rt[:nr]

    # Check which observations have already been made.
    if x_sofar is not None:

        f_unique, new_indices, N_new_indices = check_repeats(
            x_unique, x_sofar, f_sofar)

        # If there are any new input points, draw values over them.
        if N_new_indices > 0:

            x_new = x_unique[new_indices[:N_new_indices]]
            f_new = f(x_new)

            f_unique[new_indices[:N_new_indices]] = f_new

            # Record the new values
            x_sofar = vstack((x_sofar, x_new))
            f_sofar = hstack((f_sofar, f_new))
        else:
            f = f_unique

    # If no observations have been made, don't check.
    else:
        f_unique = f(x_unique)
        x_sofar = x_unique
        f_sofar = f_unique

    f = empty(lenx)
    f[unique_indices] = f_unique
    f[repeat_to] = f[repeat_from]

    return f, x_sofar, f_sofar