Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

enh: optimizing the reading of eigenstates from WFSX file #388

Merged
merged 3 commits into from
Nov 15, 2021

Conversation

pfebrer
Copy link
Contributor

@pfebrer pfebrer commented Oct 20, 2021

Addresses #147 by making python keep track of the fortran file unit.

Now if one wants to iterate over all eigenstates the whole file is only read once, in contrast to the previous implementation, which read the file from the beginning at each iteration. This quick code snippet shows the difference in performance for iterating over all eigenstates for the file spin_unpolarized.bands.WFSX, which contains 95 k points:

import sisl

wfsx = sisl.get_sile("spin_unpolarized.bands.WFSX")
geom = sisl.geom.diamond()

def iterate(which):
    it = {"new": wfsx.yield_eigenstate, "old": wfsx.old_yield_eigenstate}[which]
    
    for state in it(parent=geom):
        pass
    
%timeit iterate("new")
%timeit iterate("old")
1.87 ms ± 90.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
35.7 ms ± 97.1 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

I have some doubts:

  • I think the spin and k indices are not written on the WFSX file as the old implementation assumes (spin-loop[k-loop]). I believe it's the opposite (k-loop[spin-loop]). You can test this by trying to read a polarized WFSX file, which will fail on the second eigenstate. I attach one so that you can test it. If this is true, I need to make some changes to the old implementation.
  • You already mentioned in Siesta read wfs and wfsx files #147 that the reading is unstable. I think I stumble upon this when reading the polarized file attached here. The binary check raises an error at ik = 68. Why is this? 😅

Files (too big to attach): https://drive.google.com/drive/folders/1lgqckfkpq0mZo-3ucBxtXvz5fMKlXCQv?usp=sharing

Copy link
Owner

@zerothi zerothi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll have another closer look at your comments later, this week is havock here ;)

Great work!

sisl/io/siesta/binaries.py Outdated Show resolved Hide resolved
sisl/io/siesta/binaries.py Outdated Show resolved Hide resolved
sisl/io/siesta/binaries.py Outdated Show resolved Hide resolved
sisl/io/siesta/binaries.py Outdated Show resolved Hide resolved
sisl/io/siesta/binaries.py Outdated Show resolved Hide resolved
@zerothi
Copy link
Owner

zerothi commented Oct 20, 2021

Thanks for your comments. Indeed this shows some inconsistencies in the entire SiestaBin suite!!! :)

I agree this might be fitting in the binary sile for siesta, for coherence! Thanks for noticing.

Regarding the read_eigenstate, this is probably my mistake!!! ;)
The idea should be that all methods have a similar interface, so it should definitely be k and not ik. So this needs changing. I just recall I did this quickly because somebody wanted to use the wavefunctions. So I think it is ok to break API to make it stable.

@pfebrer
Copy link
Contributor Author

pfebrer commented Oct 20, 2021

Regarding the read_eigenstate, this is probably my mistake!!! ;)

No, no. It wasn't there. I added it and then removed it. But couldn't it be useful in some cases to be able to get an eigenstate by its index? It would use the old subroutines and therefore avoid going back and forth between fortran and python.

@zerothi
Copy link
Owner

zerothi commented Oct 20, 2021

Regarding the read_eigenstate, this is probably my mistake!!! ;)

No, no. It wasn't there. I added it and then removed it. But couldn't it be useful in some cases to be able to get an eigenstate by its index? It would use the old subroutines and therefore avoid going back and forth between fortran and python.

Since your method is so fast I think that having two methods for the same thing is not ideal. In any case the end-user would still require some method of knowing the k-point. I can't really see a really good use case that can't be solved by explicit k specification? I think we should stick with one...

@pfebrer
Copy link
Contributor Author

pfebrer commented Oct 20, 2021

Ok! Is this last change fine or was _bin_check used somewhere else?

@pfebrer
Copy link
Contributor Author

pfebrer commented Oct 21, 2021

I added the possibility to read the information (not values) of all k points at once. It is useful for debugging but could be also useful for other things (you can get a BandStructure / BrillouinZone object easily for example, or find the best match to your desired k value).

Since all information is read in one call this is significantly faster than iterating over the eigenstates:

import sisl

wfsx = sisl.get_sile("spin_unpolarized.bands.WFSX")
geom = sisl.geom.diamond()

def all_fortran():
    wfsx.read_file_info(parent=geom)

def python_loop():
    for state in wfsx.yield_eigenstate(parent=geom):
        state.info
        
%timeit all_fortran()
%timeit python_loop()
335 µs ± 9.52 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
1.99 ms ± 112 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Here is an example of what is returned by read_file_info:

{'k': array([[0.     , 0.     , 0.     ],
        [0.     , 0.0714 , 0.0714 ],
        [0.     , 0.1428 , 0.1428 ],
        [0.     , 0.2142 , 0.2142 ],
        [0.     , 0.2856 , 0.2856 ],
        [0.     , 0.357  , 0.357  ],
        [0.     , 0.4284 , 0.4284 ],
        [0.     , 0.4998 , 0.4998 ],
        [0.     , 0.5712 , 0.5712 ],
        [0.     , 0.6426 , 0.6426 ],
        [0.     , 0.714  , 0.714  ],
        [0.     , 0.7854 , 0.7854 ],
        [0.     , 0.8568 , 0.8568 ],
        [0.     , 0.9282 , 0.9282 ],
        [0.     , 0.9996 , 0.9996 ],
        [0.     , 1.071  , 1.071  ],
        [0.     , 1.1424 , 1.1424 ],
        [0.     , 1.2138 , 1.2138 ],
        [0.     , 1.2852 , 1.2852 ],
        [0.     , 1.3566 , 1.3566 ],
        [0.     , 1.428  , 1.428  ],
        [0.     , 1.4994 , 1.4994 ],
        [0.     , 1.5708 , 1.5708 ],
        [0.     , 1.6422 , 1.6422 ],
        [0.     , 1.7136 , 1.7136 ],
        [0.     , 1.785  , 1.785  ],
        [0.08925, 1.785  , 1.87425],
        [0.1785 , 1.785  , 1.9635 ],
        [0.26775, 1.785  , 2.05275],
        [0.357  , 1.785  , 2.142  ],
        [0.44625, 1.785  , 2.23125],
        [0.5355 , 1.785  , 2.3205 ],
        [0.62475, 1.785  , 2.40975],
        [0.714  , 1.785  , 2.499  ],
        [0.80325, 1.785  , 2.58825],
        [0.8925 , 1.785  , 2.6775 ],
        [0.952  , 1.785  , 2.618  ],
        [1.0115 , 1.785  , 2.5585 ],
        [1.071  , 1.785  , 2.499  ],
        [1.1305 , 1.785  , 2.4395 ],
        [1.19   , 1.785  , 2.38   ],
        [1.2495 , 1.785  , 2.3205 ],
        [1.309  , 1.785  , 2.261  ],
        [1.3685 , 1.785  , 2.2015 ],
        [1.428  , 1.785  , 2.142  ],
        [1.4875 , 1.785  , 2.0825 ],
        [1.547  , 1.785  , 2.023  ],
        [1.6065 , 1.785  , 1.9635 ],
        [1.666  , 1.785  , 1.904  ],
        [1.7255 , 1.785  , 1.8445 ],
        [1.785  , 1.785  , 1.785  ],
        [1.69575, 1.69575, 1.69575],
        [1.6065 , 1.6065 , 1.6065 ],
        [1.51725, 1.51725, 1.51725],
        [1.428  , 1.428  , 1.428  ],
        [1.33875, 1.33875, 1.33875],
        [1.2495 , 1.2495 , 1.2495 ],
        [1.16025, 1.16025, 1.16025],
        [1.071  , 1.071  , 1.071  ],
        [0.98175, 0.98175, 0.98175],
        [0.8925 , 0.8925 , 0.8925 ],
        [0.80325, 0.80325, 0.80325],
        [0.714  , 0.714  , 0.714  ],
        [0.62475, 0.62475, 0.62475],
        [0.5355 , 0.5355 , 0.5355 ],
        [0.44625, 0.44625, 0.44625],
        [0.357  , 0.357  , 0.357  ],
        [0.26775, 0.26775, 0.26775],
        [0.1785 , 0.1785 , 0.1785 ],
        [0.08925, 0.08925, 0.08925],
        [0.     , 0.     , 0.     ],
        [0.1071 , 0.1071 , 0.1071 ],
        [0.2142 , 0.2142 , 0.2142 ],
        [0.3213 , 0.3213 , 0.3213 ],
        [0.4284 , 0.4284 , 0.4284 ],
        [0.5355 , 0.5355 , 0.5355 ],
        [0.6426 , 0.6426 , 0.6426 ],
        [0.7497 , 0.7497 , 0.7497 ],
        [0.8568 , 0.8568 , 0.8568 ],
        [0.9639 , 0.9639 , 0.9639 ],
        [1.071  , 1.071  , 1.071  ],
        [1.1781 , 1.1781 , 1.1781 ],
        [1.2852 , 1.2852 , 1.2852 ],
        [1.3923 , 1.3923 , 1.3923 ],
        [1.4994 , 1.4994 , 1.4994 ],
        [1.6065 , 1.6065 , 1.6065 ],
        [1.7136 , 1.7136 , 1.7136 ],
        [1.8207 , 1.8207 , 1.8207 ],
        [1.9278 , 1.9278 , 1.9278 ],
        [2.0349 , 2.0349 , 2.0349 ],
        [2.142  , 2.142  , 2.142  ],
        [2.2491 , 2.2491 , 2.2491 ],
        [2.3562 , 2.3562 , 2.3562 ],
        [2.4633 , 2.4633 , 2.4633 ],
        [2.5704 , 2.5704 , 2.5704 ],
        [2.6775 , 2.6775 , 2.6775 ]]),
 'weight': array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]),
 'nwf': array([[8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
         8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
         8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
         8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
         8, 8, 8, 8, 8, 8, 8, 8]], dtype=int32)}

On the 505 MB file that I attached, reading the info takes 140 ms. So I'd say this is very usable in practice even for huge WFSX files.

Copy link
Owner

@zerothi zerothi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll have another closer look next week on this. Looks interesting.

sisl/io/siesta/_src/wfsx_read.f90 Outdated Show resolved Hide resolved
sisl/io/siesta/_src/wfsx_read.f90 Outdated Show resolved Hide resolved
@pfebrer
Copy link
Contributor Author

pfebrer commented Oct 21, 2021

Ok, all this is because I want to allow using WFSX files for plotting wavefunctions, (fat)bands and PDOS.

For the wavefunctions and (fat)bands I will just add entry points in the respective plots. For the PDOS, do you think it would be useful to provide a method outside the plotting suite as well? E.g. : wfsxSileSiesta.read_PDOS().

Or maybe in a more abstract way provide apply dispatchers that would work similarly as the BrillouinZone applies? (I like this idea)

@pfebrer
Copy link
Contributor Author

pfebrer commented Oct 21, 2021

In general, an eigenstate iterator may be interesting(?)

Then the code could be generalized and the only thing to change between reading from different sources would be the generator.

EigenstateIterator(provider=wfsx.yield_eigenstate(), wrap=calc_PDOS)

And if you wanted to do the same with the hamiltonian/BZ you just change the provider:

EigenstateIterator(provider=bz.apply.iter.eigenstate(), wrap=calc_PDOS)

In the future, for any other code that has a file with eigenstate outputs, it would be trivial to implement the calculation of PDOS.

In the raw case it seems pretty stupid but you would have iterators ready to calculate the PDOS, the bands, etc... (e.g. parallelization could be built-in) I don't know exactly how it should work or if it would really be useful, but just dropping the idea here 😅

@pfebrer
Copy link
Contributor Author

pfebrer commented Oct 28, 2021

I have added the possibility to read the basis. Although there's not much information in it, it might be useful sometimes as a last resort, who knows (also it probably can be done much more efficiently).

Also, I have reorganized things so that they are (hopefully) easier to read and understand. All functions that need external handling of the file unit are called _read_next_* in python and read_wfsx_next_* in fortran.

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 8, 2021

Can I generate a BandStructure object from a list of k points? I want to be able to plot bands and fatbands just from the WFSX file, but I am missing some points 😅

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 9, 2021

Also, one question for PDOS. If I want to generate the PDOS from the WFSX file, I need to read the overlap matrix as well, right?

@zerothi
Copy link
Owner

zerothi commented Nov 9, 2021

Can I generate a BandStructure object from a list of k points? I want to be able to plot bands and fatbands just from the WFSX file, but I am missing some points sweat_smile

Yes, I can't see why not. Except that the bandstructure object requires a parent with cell information.

Also, one question for PDOS. If I want to generate the PDOS from the WFSX file, I need to read the overlap matrix as well, right?

Yes, you need the overlap matrix as well.

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 9, 2021

Yes, I can't see why not. Except that the bandstructure object requires a parent with cell information.

So just pass the list of points to the point argument and set division to 1? It would be nice to have something to detect corners for this case, wouldn't it? 😄

@zerothi
Copy link
Owner

zerothi commented Nov 9, 2021

Yes, I can't see why not. Except that the bandstructure object requires a parent with cell information.

So just pass the list of points to the point argument and set division to 1? It would be nice to have something to detect corners for this case, wouldn't it? smile

Hmm, I see. It may not be so viable after all... Is there a particular reason for not just using the BrillouinZone object directly? I.e. what methods from BandStructure do you need?

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 9, 2021

I think I just needed BandStructure.lineark to plot the bands proportionally. I mean that dx in the bands plot is a constant amount.

And BandStructure.lineartick to get the ticks, but this is probably not needed since there are no ticks in the WFSX file anyway.

@zerothi
Copy link
Owner

zerothi commented Nov 9, 2021

I think I just needed BandStructure.lineark to plot the bands proportionally. I mean that dx in the bands plot is a constant amount.

And BandStructure.lineartick to get the ticks, but this is probably not needed since there are no ticks in the WFSX file anyway.

Latest commit has a linspace_bz which converts BZ object k-points into linear distances (accepts any BZ-object).

@pfebrer pfebrer force-pushed the optimize_wfsx branch 3 times, most recently from 062c642 to 854079f Compare November 10, 2021 17:38
@zerothi
Copy link
Owner

zerothi commented Nov 10, 2021

Could we perhaps limit this PR to the reading of the eigenstates, you are adding too many things that ought to be in a separate MR ;)

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 11, 2021

Hmm ok, but I wanted to make sure that the API is useful for creating the plots, which was my initial motivation.

Maybe when I finish with the plots I can just keep the changes in the wfsx sile if you want.

@zerothi
Copy link
Owner

zerothi commented Nov 11, 2021

Hmm ok, but I wanted to make sure that the API is useful for creating the plots, which was my initial motivation.

Maybe when I finish with the plots I can just keep the changes in the wfsx sile if you want.

It is fine for now, but generally it would be better to create a new branch off of this one and test it there. Smaller merge requests that are self-contained is much better and easier to manage.

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 11, 2021

This is ready for review! After the review, all that's left is to add tests and remove the old reading function.

I squashed everything into two commits so that you can check it easier: one for WFSX reading and the other to incorporate it into the plots.

Generally one should limit the return types as dictionaries
as it imposes naming conventions on certain variables
that may not be optimal for the end user. In these cases
namedtuple returns are better since users can do whatever they
wish with them.

Removed lots of not-used variables in the fortran code.

Changed names of many of the _bin_* methods to explicitly
mention the fortran interface. This clarifies belonging.

Also ran pep8 which fixed lots of indentation problems.
@zerothi
Copy link
Owner

zerothi commented Nov 15, 2021

I have pushed a new commit which you can check, also added a test for this.

@zerothi
Copy link
Owner

zerothi commented Nov 15, 2021

Ok, all this is because I want to allow using WFSX files for plotting wavefunctions, (fat)bands and PDOS.

For the wavefunctions and (fat)bands I will just add entry points in the respective plots. For the PDOS, do you think it would be useful to provide a method outside the plotting suite as well? E.g. : wfsxSileSiesta.read_PDOS().

Or maybe in a more abstract way provide apply dispatchers that would work similarly as the BrillouinZone applies? (I like this idea)

No, I don't think we should plaster the file with this, it can be done with f = wfsxSileSiesta(parent=Hamiltonian) ; for state in f.yield_eigenstate(): state.PDOS(...)

I clearly made a mistake in the Hamiltonian object to add DOS and PDOS methods there. Better to have these decoupled.

@zerothi zerothi merged commit 1f10e1a into zerothi:main Nov 15, 2021
@zerothi
Copy link
Owner

zerothi commented Nov 15, 2021

Note that now each method does not have the parent argument. Instead the creation of the object has this interface:

wfsxSileSiesta(..., parent=?, geometry=?, sc=?)

where the sc is automatically fetched from parent if supplied. This makes calls more strict and easier. :)

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 15, 2021

Great! I'll submit another PR in a moment then because in the plots I was using yield_eigenstate(parent=...) :)

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 15, 2021

Also, could you add spin-polarized and noncolinear spin versions of the .bands.WFSX file if you have time? I would like to add tests for all the spin possibilities in plots

@zerothi
Copy link
Owner

zerothi commented Nov 15, 2021

Great! I'll submit another PR in a moment then because in the plots I was using yield_eigenstate(parent=...) :)

Great, I have a helper script pep8.sh that I think you should run before MR's, quite often it catches indentation problems in some of the plotting regions.

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 15, 2021

Great, I have a helper script pep8.sh that I think you should run before MR's, quite often it catches indentation problems in some of the plotting regions.

Yes, I'll do that

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 16, 2021

Generally one should limit the return types as dictionaries
as it imposes naming conventions on certain variables
that may not be optimal for the end user. In these cases
namedtuple returns are better since users can do whatever they
wish with them.

I remember why I returned dictionaries. Dictionaries allow returning additional items in the future without breaking backwards compatibility. If you return a tuple, you are basically stuck with the choice you made in the first place. What if there's an extra size that you think could be useful to return in the future? If you make the tuple larger, working code will get errors when updating the sisl version, which can be annoying.

@zerothi
Copy link
Owner

zerothi commented Nov 16, 2021

I remember why I returned dictionaries. Dictionaries allow returning additional items in the future without breaking backwards compatibility. If you return a tuple, you are basically stuck with the choice you made in the first place. What if there's an extra size that you think could be useful to return in the future? If you make the tuple larger, working code will get errors when updating the sisl version, which can be annoying.

Yes, but dictionaries is really annoying to work with when dealing with simple data. I agree that the flexibility is good. But I am thinking about the end users who are needing to write code. And using dictionaries is in my opinion a needless hurdle. It obfuscates code and complicates things.

for out in ...iter:
    no_u = out["no_u"]
    data = out["data"]

could be done with

for no_u, data in ...iter:
    ...

it is clear and meaningful. Also when you see other big, succesfull code projects (numpy, scipy). They only ever return dictionaries when you know you have variable data associated. And here we must have some kind of consistency.

If we later wishes to return something else, then this could be solved by adding flags to the sile itself. wfsxSile.VERSION or something like that.
Secondly, how often do you think this code will change? I don't think it will actually change a lot during the next several years ;)

@zerothi
Copy link
Owner

zerothi commented Nov 16, 2021

Actually, this is why I would prefer namedtuple returns, it has a direct _todict method allowing you to use it directly as a dict.

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 16, 2021

Secondly, how often do you think this code will change? I don't think it will actually change a lot during the next several years ;)

No, but it may change 10 years from now and then break codes whose authors may not be maintaining anymore.

Actually, this is why I would prefer namedtuple returns, it has a direct _todict method allowing you to use it directly as a dict.

I think this could be the best option.

@zerothi
Copy link
Owner

zerothi commented Nov 16, 2021

Secondly, how often do you think this code will change? I don't think it will actually change a lot during the next several years ;)

No, but it may change 10 years from now and then break codes whose authors may not be maintaining anymore.

Actually, this is why I would prefer namedtuple returns, it has a direct _todict method allowing you to use it directly as a dict.

I think this could be the best option.

You can see that that is what I have done in this PR

@pfebrer
Copy link
Contributor Author

pfebrer commented Nov 16, 2021

Aah ok yes, I was looking at read_info, not read_sizes 👍

@pfebrer pfebrer deleted the optimize_wfsx branch December 8, 2021 10:22
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants