import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

def objective_function(t):
    X, Y = t
    return 1-1/((1+X**2+Y**2)*(1+np.sin(2*X)**2+1.5*np.sin(2*Y)**2))
    # return 1-1/((1+X**2+Y**2))

optimization_path = []

def objective_function_with_memory(t):
    x, y = t
    z = objective_function(t)
    optimization_path.append((x,y,z))
    return z

t0 = np.array([np.pi/4, np.pi/4])
# t0 = np.array([np.pi/2, np.pi/4])

result = minimize(objective_function_with_memory, 
                  t0, method='COBYLA', options={'maxiter': 25} )

print(result)
print(len(optimization_path))

xpath, ypath, zpath = zip(*optimization_path)


x_vals = np.linspace(-3*np.pi/4, 3*np.pi/4, 200)
y_vals = np.linspace(-3*np.pi/4, 3*np.pi/4, 200)
X, Y = np.meshgrid(x_vals, y_vals)
Z = objective_function([X,Y])

fig = plt.figure(figsize=(13, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.5) #type: ignore
ax.scatter(result.x[0], result.x[1], 
           result.fun, color='red', marker='x')

ax.plot(xpath, ypath, zpath, color='purple')
ax.scatter(xpath, ypath, zpath, color='purple', marker='o')
plt.show()

 message: Maximum number of function evaluations has been exceeded.
 success: False
  status: 2
     fun: 1.1292277039109955e-05
       x: [-1.397e-03  4.678e-04]
    nfev: 25
   maxcv: 0.0
25

steps = list(range(len(zpath)))
plt.plot(steps, zpath, marker='o', linestyle='-')
# plt.plot(steps[20:], zpath[20:], marker='o', linestyle='-')
plt.grid(True)
plt.show()