How to use the numdifftools.nd_algopy.Hessian function in numdifftools

def _run_hamiltonian(verbose=True):
    c = classicalHamiltonian()
    if verbose:
        print(c.potential(array([-0.5, 0.5])))
        print(c.potential(array([-0.5, 0.0])))
        print(c.potential(array([0.0, 0.0])))

    xopt = optimize.fmin(c.potential, c.initialposition(), xtol=1e-10)

    hessian = nd.Hessian(c.potential)

    H = hessian(xopt)
    true_H = np.array([[5.23748385e-12, -2.61873829e-12],
                       [-2.61873829e-12, 5.23748385e-12]])
    error_estimate = np.NAN
    if verbose:
        print(xopt)
        print('H', H)
        print('H-true_H', np.abs(H - true_H))
        # print('error_estimate', info.error_estimate)

        eigenvalues = linalg.eigvals(H)
        normal_modes = c.normal_modes(eigenvalues)

        print('eigenvalues', eigenvalues)
        print('normal_modes', normal_modes)

def test_hessian_cosIx_yI_at_I0_0I(self):
        # cos(x-y), at (0,0)

        def fun(xy):
            return np.cos(xy[0] - xy[1])
        htrue = [[-1., 1.], [1., -1.]]
        methods = ['forward', ]  # 'reverse']

        for method in methods:
            Hfun2 = nd.Hessian(fun, method=method)
            h2 = Hfun2([0, 0])
            # print(method, (h2-np.array(htrue)))
            assert_array_almost_equal(h2, htrue)

#
    # Scientific
    f = benchmark1.F(N)
    t = time.time()
    hessian = scientific.Hessian(f)
    preproc_time = time.time() - t
    t = time.time()
    H = hessian(3 * np.ones(N))
    run_time = time.time() - t
    results_hessian[method['scientific']] = run_time, np.linalg.norm(
        (H - ref_H).ravel()) / np.linalg.norm((ref_H).ravel()), preproc_time

    # algopy forward utpm variant
    f = benchmark1.F(N)
    t = time.time()
    hessian = algopy.Hessian(f, method='forward')
    preproc_time = time.time() - t
    t = time.time()
    H = hessian(3 * np.ones(N))
    run_time = time.time() - t
    results_hessian[method['algopy_forward']] = run_time, np.linalg.norm(
        (H - ref_H).ravel()) / np.linalg.norm((ref_H).ravel()), preproc_time

    # numdifftools
    f = benchmark1.F(N)
    t = time.time()
    hessian = numdifftools.Hessian(f, **options)
    preproc_time = time.time() - t
    t = time.time()
    H = hessian(3 * np.ones(N))
    run_time = time.time() - t
    results_hessian[method['numdifftools']] = run_time, np.linalg.norm(

# HESSIAN COMPUTATION
# -------------------
print 'starting hessian computation '
results_hessian_list = []
hessian_N_list = gradient_N_list
# hessian_N_list = [1, 2, 4, 8, 16, 32, 64]
# hessian_N_list = [2]

for N in hessian_N_list:
    print 'N=', N
    results_hessian = np.zeros((4, 3))

    f = benchmark1.F(N)
    t = time.time()
    hessian = algopy.Hessian(f,  method='forward')
    preproc_time = time.time() - t
    t = time.time()
    ref_H = hessian(3 * np.ones(N))
    run_time = time.time() - t
    results_hessian[method['algopy_reverse']] = run_time, 0.0, preproc_time

#
    # Scientific
    f = benchmark1.F(N)
    t = time.time()
    hessian = scientific.Hessian(f)
    preproc_time = time.time() - t
    t = time.time()
    H = hessian(3 * np.ones(N))
    run_time = time.time() - t
    results_hessian[method['scientific']] = run_time, np.linalg.norm(