Mani Chandra wrote: > We're dealing with a 5000 x 5000 site system, and the time to > interpolate the densities and currents using > kwant.plotter.interpolate_density and > kwant.plotter.interpolate_current are by far the rate limiting steps > in the calculation:
You must be applying Kwant’s interpolator to the biggest system anyone tried ever! Why are you interpolating such a huge system? What are you trying to achieve? What is the width of the interpolation “bumps” that you are using and what is the size of your output array? (Or in other words: how are you calling the interpolator?) > Is there any way to speed up the interpolation routines? I see that > the present implementation in plotter.py involves a for loop (in > python) over each site/hopping. Are you using streamplot or just interpolating currents/densities? Kwant’s streamplot routine piggy backs on matplotlib’s streamplot and in a typical use case (e.g. the attached script) matplotlib’s streamplot will take half of the time that Kwant’s interpolation routine takes. So even if the interpolator produced its results instantly, current plotting would become only three times faster in such applications. The interpolator could be optimized further in several ways - as always it’s a compromise between development effort and speed. But the application that it was created for (i.e. figures for publications) tends to involve systems of size 100^2 or so and then its speed is quite satisfactory (see attached example). Creating figures from much larger systems only makes sense when the interpolation is smoothing things out a lot (i.e. the “bumps” are much larger than a hopping), otherwise there would be too much detail in the figure. In this regime the current implementation will be inefficient because its vectorization (see further below) will be over small subarrays. > Perhaps a vectorized implementation may help? Actually, it is already vectorized - otherwise it would be much slower! The main loop inside _interpolate_field() goes indeed sequentially over all hoppings/densities, but for each such hopping/density the output field in its vicinity is updated using vectorized numpy operations. A straightforward brute-force way to speed up the interpolation would be to split up the work of that loop over multiple cores - each would have its own output array, and in the end they would be all summed together. That should be quite easy to do [1]. A complementary approach would be to rewrite the inner loop in a fast compiled language. This would help especially in the case where the interpolation bump does not cover many points of the output array and as such the existing vectorization is not very good (the involved subarrays are small). Yet another approach would be to create a special implementation for the common case of a rectangular lattice that is commensurate with the output array of the interpolator. All three approaches could be combined. To go further, someone who is really motivated could even implement the interpolation algorithm on GPUs... If you or someone else would like to work on this I would be happy to guide you. Cheers Christoph
import math
from cmath import exp
import numpy
from matplotlib import pyplot
import kwant
from kwant.digest import gauss
def hopping(sitei, sitej, phi, salt):
xi, yi = sitei.pos
xj, yj = sitej.pos
return -exp(-0.5j * phi * (xi - xj) * (yi + yj))
def onsite(site, phi, salt):
return 0.05 * gauss(repr(site), salt) + 4
def make_system(L=50):
def central_region(pos):
x, y = pos
return -L < x < L and \
abs(y) < L - 37.5 * math.exp(-x**2 / 12**2)
lat = kwant.lattice.square(norbs=1)
sys = kwant.Builder()
sys[lat.shape(central_region, (0, 0))] = onsite
sys[lat.neighbors()] = hopping
sym = kwant.TranslationalSymmetry((-1, 0))
lead = kwant.Builder(sym)
lead[(lat(0, y) for y in range(-L + 1, L))] = 4
lead[lat.neighbors()] = hopping
sys.attach_lead(lead)
sys.attach_lead(lead.reversed())
return sys.finalized()
sys = make_system()
energy = 0.15
kwant.plot(sys)
args = [1/40.0, ""]
psi = kwant.wave_function(sys, energy, args)(0)
J = kwant.operator.Current(sys).bind(args)
current = sum(J(p) for p in psi)
kwant.plotter.current(sys, current)
[1] https://medium.com/analytics-vidhya/using-numpy-efficiently-between-processes-1bee17dcb01
smime.p7s
Description: S/MIME cryptographic signature
