Incorporate filters to enhance LD plots appearance

[broncio123]

For the future, this can be done by defining our own filters, OR Use other filters implemented by other Python libraries for image-processing.

[broncio123]

So far, the use of gradient-based filters for extraction of manifold curves from LD maps has been explored. See Issues #17 #18 #20 #25

However, the above seems ineffective to provide neat curves as graphical outputs highlighting the distinction of stable and unstable manifolds. Instead, the use of a Laplacian-based filter - which relies on second derivatives, seem to be more effective to this end.

@vkrajnak ... below, I provide example code-lines implementing this, alongside their graphical output. This was found thanks to @VikJGG who kindly shared with me his implementation in MatLab.

After computing LD values

LD_forward = compute_lagrangian_descriptor(grid_parameters, vector_field, tau, p_value)
LD_backward = compute_lagrangian_descriptor(grid_parameters, vector_field, -tau, p_value)

We used pcolor instead of contourf (Need to verify if this makes a difference)

import matplotlib.pyplot as plt
from scipy.ndimage import laplace

fig, (ax1,ax2) = plt.subplots(1,2,figsize=(8,4),dpi=200,sharex=True,sharey=True)

ax1_min, ax1_max, N1 = slice_parameters[0]
ax2_min, ax2_max, N2 = slice_parameters[1]
points_ax1 = np.linspace(ax1_min, ax1_max, N1)
points_ax2 = np.linspace(ax2_min, ax2_max, N2)

X, Y = np.meshgrid(points_ax1, points_ax2)

LD = LD_forward + LD_backward
ax1.contourf(points_ax1, points_ax2,LD,cmap='bone',levels=200)

tol = 0.01
#compute laplacian
LD = LD_forward
scalar = laplace(LD)
scalar = scalar/scalar.max()
scalar = scalar**2
LDm = np.ma.masked_where(scalar < tol, scalar )
ax2.pcolor(X, Y, LDm, cmap='bwr')

# compute laplacian
LD = LD_backward
scalar = laplace(LD)
scalar = scalar/scalar.max()
scalar = scalar**2
LDm = np.ma.masked_where(scalar < tol, scalar )
ax2.pcolor(X, Y, -LDm, cmap='bwr')

#ax1.set_xlim(-4,4)
#ax1.set_ylim(-1,1)
ax1.set_xlim(-5.5, 5.5)
ax1.set_ylim(-1.25, 1.25)

fig.suptitle("$ PES(k = k_c), \Delta H_0 = "+str(dH)+" :  LD_p^{total}, p = 1/2, \\tau$ = "+str(tau), fontsize=14, y=1.04)
ax1.set_xlabel("$x$")
ax2.set_xlabel("$x$")
ax1.set_ylabel("$p_x$")

fig.tight_layout()
plt.show()

[broncio123]

@vkrajnak Just remembered I posted some code where I show how to use the Laplacian to highlight the location of invariant manifolds. Maybe keep this in mind when working on the implementation of filters.

[vkrajnak]

A proof of concept can be found on the max_gradient branch