import numpy as np
from scipy.differentiate import derivative, hessian, jacobian


x0 = np.pi / 4
scalar = derivative(np.sin, x0)

print(f"derivative: {scalar.df:.8f}")
print(f"expected: {np.cos(x0):.8f}")
print(f"status: {int(scalar.status)}")
print(f"error_estimate: {scalar.error:.2e}")


def vector_function(x):
    return np.array([
        x[0]**2 + x[1],
        x[0] * x[1],
    ])


point = np.array([2.0, 3.0])
expected_jacobian = np.array([
    [4.0, 1.0],
    [3.0, 2.0],
])
jac = jacobian(vector_function, point)

print("jacobian:")
for index, row in enumerate(jac.df):
    print(f"row {index}: {row}")
print("jacobian success:", bool(jac.success.all()))
print(f"jacobian max_abs_diff: {np.max(np.abs(jac.df - expected_jacobian)):.2e}")


def quadratic(x):
    return x[0]**2 + 3 * x[0] * x[1] + 2 * x[1]**2


expected_hessian = np.array([
    [2.0, 3.0],
    [3.0, 4.0],
])
hes = hessian(quadratic, point)

print("hessian:")
for index, row in enumerate(hes.ddf):
    print(f"row {index}: {row}")
print("hessian success:", bool(hes.success.all()))
print(f"hessian max_abs_diff: {np.max(np.abs(hes.ddf - expected_hessian)):.2e}")