TRIQS / triqs

a Toolbox for Research on Interacting Quantum Systems
https://triqs.github.io
GNU General Public License v3.0
142 stars 72 forks source link

improve `.values()` on meshes #958

Open the-hampel opened 3 months ago

the-hampel commented 3 months ago

After adding the recent feature to obtain a np.array of all mesh values via .values() (81abab72d62846c4558c2af169ef4f9a5b68d5ef and 1ba26952d4ed8af054a2c29c7ec91ee1874475ab) it would be helpful to improve the functionality for MeshImFreq further. Right now calling .values() on a MeshImFreq gives back an array of mesh points:

mesh = MeshImFreq(beta=40, statistic='Fermion', n_iw=2)
print(mesh.values())
mesh_cmplx = 
> [MatsubaraFreq(n: -2, beta: 40.0, statistic: Fermion) MatsubaraFreq(n: -1, beta: 40.0, statistic: Fermion) MatsubaraFreq(n: 0, beta: 40.0, statistic: Fermion) MatsubaraFreq(n: 1, beta: 40.0, statistic: Fermion)]

It would be helpful to return back directly a list of complex numbers. This can right now be achieved only via:

np.vectorize(lambda x: x.imag)(mesh.values())
> array([-0.235619, -0.07854 ,  0.07854 ,  0.235619])

Calling .imag directly on the array does not work, as numpy does not know how to use the dispatcher on a triqs mesh point object. I tried to find a solution for this on the python level but it seems no so straight forward to teach numpy to do this. If someone knows this would be the easiest! Other options: (a) create a class that inherits from np.ndarray that is returned when .values() is called. Then this new class has a member attribute called .imag that will just do the np.vectorize call. (b) implement a real C++ function that returns the values as an array (right now this only lives in the meshes_desc file) and add an optional keyword that switches between returning meshpoints or cmplx numbers, i.e. something like this:

print(mesh.values(return_mesh_points=True))
> [MatsubaraFreq(n: -2, beta: 40.0, statistic: Fermion) MatsubaraFreq(n: -1, beta: 40.0, statistic: Fermion) MatsubaraFreq(n: 0, beta: 40.0, statistic: Fermion) MatsubaraFreq(n: 1, beta: 40.0, statistic: Fermion)]
print(mesh.values(return_mesh_points=False))
> array([-0.235619, -0.07854 ,  0.07854 ,  0.235619])

where the default can be to return mesh points.

krivenko commented 3 months ago

Hey Alexander,

Here are a couple remarks that you might find useful (or not).

This can right now be achieved only via:

np.vectorize(lambda x: x.imag)(mesh.values())
array([-0.235619, -0.07854 ,  0.07854 ,  0.235619])

This can be shortened to np.vectorize(np.imag)(mesh.values()).

Also, one does not really need the .values() method to get the imaginary parts.

map(np.imag, mesh) # -> an iterator returning float
list(map(np.imag, mesh)) # -> a list of float
np.fromiter(map(np.imag, mesh), float) # -> an np.ndarray of float

Calling .imag directly on the array does not work, as numpy does not know how to use the dispatcher on a triqs mesh point object. I tried to find a solution for this on the python level but it seems no so straight forward to teach numpy to do this.

How about simply defining a free function?

def imag(mesh: MeshImFreq):
    return np.fromiter(map(np.imag, mesh), float)