code: def caching_call(f, x, x_sofar, 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
docs: Computes f(x) given that f(x_sofar) = x_sofar.
      returns f(x), and new versions of x_sofar and f_sofar.