def METHOD_NAME(sites, ops, hkl, obs, eps=1.e-8):
  result = flex.double()
  sites_eps = list(sites)
  for js,site in enumerate(sites):
    site_eps = list(site)
    for jp in range(3):
      vs = []
      for signed_eps in [eps, -eps]:
        site_eps[jp] = site[jp] + signed_eps
        sites_eps[js] = matrix.col(site_eps)
        sf = exp_i_hx(sites=sites_eps, ops=ops, hkl=hkl)
        tf = least_squares(obs=obs, calc=sf.f())
        vs.append(tf.f())
      result.append((vs[0]-vs[1])/(2*eps))
    sites_eps[js] = sites[js]
  return result