Ensure that initialize_lebedev is called --> Clean up spaghetti in cubature #2736

tallakahath · 2022-10-05T02:08:05Z

Previously, initialize_lebedev was never called, and in fact was getting optimized out of the module completely upon compilation. When lebedev_mappping_[] is then accessed across multiple OpenMP threads, the std::map is empty, and a deadlock can happen where two threads try to access-write (since [key] fills if key is not found), and the slightly slower thread ends up in a Bad State where it thinks there is a value but ends up infinitely looping on the lookup (the program will hang on []).

This only happens once every several thousand runs, and only when running with a high degree of parallelism in a system with many atoms. I cannot induce it in captivity, but I have observed it in the wild.

Anyway, [] accesses on std::map aren't thread-safe if you aren't super-duper sure the map is fully filled for all keys you'd ever look up, which should be the case if initialize_lebedev was ever called anywhere. But it wasn't, and that was Bad.

Now it's called exactly once (thanks, c++11's call_once! I do see that this isn't used anywhere else in the code, but I do see mutex is imported in several files, so I don't think I'm adding any new deps here).

The hangs should be gone, though I'll have to churn through another several thousand runs to likely be sure (as, again, it is a very rare kind of hang). This will take me a few days to confirm, but given all debugging efforts point to this being the problem, I'm like 99% confident this will do the trick.

That said, as far as I can tell, besides one print function the resulting order_ that's assigned to is never used. Maybe a candidate to be axed in the future?

Description

Actually invokes initialize_lebedev before accessing lebedev_mapping_ to ensure the mapping has values, and prevents a deadlock when running in parallel

User API & Changelog headlines

Prevents a nasty, rare hang

Dev notes & details

See the main PR body

Questions

What does order_ actually do in SphericalGrid? It never appears to be used anywhere except one print function that also appears unused.

Checklist

No tests needed -- no observable changes and it's hard to write a test for a very rare race condition :)

Status

Ready for review
Ready for merge

Previously, initialize_lebedev was never called, and in fact was getting optimized out of the module completely upon compilation. When lebedev_mappping_[] is then accessed across multiple OpenMP threads, the std::map is empty, and a deadlock can happen where two threads try to access-write (since [key] fills if key is not found), and the slightly slower thread ends up in a Bad State where it thinks there is a value but ends up infinitely looping on the lookup (the program will hang on []). This only happens once every several thousand runs, and only when running with a high degree of parallelism in a system with many atoms. I cannot induce it in captivity, but I have observed it in the wild. Anyway, [] accesses on std::map aren't thread-safe if you aren't super-duper sure the map is fully filled for all keys you'd ever look up, which *should* be the case if initialize_lebedev was ever called anywhere. But it wasn't, and that was Bad. Now it's called exactly once (thanks, c++11's call_once!). The hangs should be gone, though I'll have to churn through another several thousand runs to likely be sure (as, again, it is a very rare kind of hang). That said, as far as I can tell, besides one print function the resulting order_ that's assigned to is never *used*. Maybe a candidate to be axed in the future?

susilehtola · 2022-10-05T08:58:53Z

lebedev_mapping_ is a member of SphericalGrid, so initialize_lebedev() should be called in the constructor of SphericalGrid. No need to add mutexes etc.

More worryingly, I also see another lebedev_mapping_ in cubature.cc

psi4/psi4/src/psi4/libfock/cubature.cc

Line 5065 in ac8f87a

std::map<int, int> SphericalGrid::lebedev_mapping_;

Or have I misunderstood something basic? 😆

loriab · 2022-10-05T13:36:25Z

Thanks for the thorough report and fix. The seemingly extra declaration @susilehtola pointed out confused me, too, though perhaps https://stackoverflow.com/a/17392441 is the answer. The map looks to be from the boost-to-avoid-c++11-standard era, so it could be improved (https://github.com/psi4/psi4archive/blame/1.0.x/src/lib/libfock/cubature.h#L302).

tallakahath · 2022-10-05T13:41:10Z

The "other" lebedev_mapping_: That's just the declaration for it. The problem is that `SphericalGrid::build` is a _static_ method so the constructor isn't called when that happens. Someone might yank out the `new` there and still statically access `lebedev_mapping_` and we are back where we started. There's no guarantee of construction, so I stapled it into the one place it's actually used. (And, again, I don't understand why it's used at all as it never seems to show up downstream) Even if it was in the constructor though (which would fire on the _new_) it would still need the mutex to ensure it's initialized once as lebedev_mapping_ is also (purposefully) static. No sense doing the rebuild of the map on every single object instantiation.

…

On Wed, Oct 5, 2022, 4:59 AM Susi Lehtola ***@***.***> wrote: lebedev_mapping_ is a member of SphericalGrid, so initialize_lebedev() should be called in the constructor of SphericalGrid. No need to add mutexes etc. More worryingly, I also see another lebedev_mapping_ in cubature.cc https://github.com/psi4/psi4/blob/ac8f87a1dd3fdda2aabc3318713d6e5ce00e2c70/psi4/src/psi4/libfock/cubature.cc#L5065 — Reply to this email directly, view it on GitHub <#2736 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABTN7JCEAMFQBY2H6OKMPYLWBU7NRANCNFSM6AAAAAAQ5DTWOY> . You are receiving this because you authored the thread.Message ID: ***@***.***>

susilehtola · 2022-10-05T15:37:34Z

The initialization code is pretty simple, however. One could just replace the initialize_lebedev function by an initializer list in

psi4/psi4/src/psi4/libfock/cubature.cc

Line 5065 in ac8f87a

std::map<int, int> SphericalGrid::lebedev_mapping_;

This would also avoid the need for the mutexes.

The code is hacky also elsewhere; lebedev_error duplicates the same information. For what it's worth, I will be developing a general C++ Becke quadrature library that would eliminate these parts of the code.

tallakahath · 2022-10-05T16:13:41Z

I will admit I considered that, then decided I didn't want to reformat the map into the initializer because it was long and I was lazy.

My bigger concern is still: what was the consequence of order_ being incorrectly 0 for so long, and who actually uses this property up-stream? The git commit where this was added says Hack DFT grids to retain indexing data for ISA but the only files touched were:

src/lib/libfock/cubature.cc.REMOVED.git-id
src/lib/libfock/cubature.h
src/lib/libfock/gridblocker.h

And I can't find anywhere that seems to use order_ or the order() method of SphericalGrid.

Unfortunately, trying to look further back in the history of cubature.cc results in several instances of psi4/src/psi4/libfock/cubature.cc.REMOVED.git-id which (I think?) means things were stripped from the repo at some point, and the resulting commit stuff is all assigned to the wrong person (e.g., nearly all of cubature.cc is assigned to one commit 7e4ecf9 from dgasmith).

print_details is the only place the spherical_grids_ are ever apparently used, and I can comment out the block where they print nicely. That function is only ever used in v.cc and only if print_ > 2. I don't think this is used much, and is creating unnecessary complexity, and doesn't feature in any tests I can find.

I'm running through tests now to see if just... removing this entirely breaks anything. Short of print_details being accessible through some print settings, it doesn't appear that any of this is even exposed on the python side anyway?

susilehtola · 2022-10-05T16:17:59Z

I'm running through tests now to see if just... removing this entirely breaks anything. Short of print_details being accessible through some print settings, it doesn't appear that any of this is even exposed on the python side anyway?

Haha yes spaghetti code!

loriab · 2022-10-05T16:22:32Z

Anyone wanting to look back in psi4 history beyond 2016 (when history was rewritten to shrink the codebase by 90% by getting rid of boost tarballs) would need to use the psi4/psi4archive repo, https://github.com/psi4/psi4archive/blame/1.0.x/src/lib/libfock/cubature.h#L302 .

Yes, acknowledged spaghetti :-(

JonathonMisiewicz · 2022-10-05T16:34:42Z

My bigger concern is still: what was the consequence of order_ being incorrectly 0 for so long, and who actually uses this property up-stream? The git commit where this was added says Hack DFT grids to retain indexing data for ISA but the only files touched were:
src/lib/libfock/cubature.cc.REMOVED.git-id
src/lib/libfock/cubature.h
src/lib/libfock/gridblocker.h
And I can't find anywhere that seems to use order_ or the order() method of SphericalGrid.

The original coder probably added it anticipating it would be useful, but it seems to have never been used. I'm happy to see unused code burn.

tallakahath · 2022-10-05T16:41:52Z

100% tests passed, 0 tests failed out of 520

The original coder probably added it anticipating it would be useful, but it seems to have never been used. I'm happy to see unused code burn.

Burn, baby, burn. Patch incoming...

tallakahath · 2022-10-05T16:57:22Z

RadialGrid also looks ripe for burning -- smoketests passed, running full tests now.

tallakahath · 2022-10-05T17:57:54Z

100% tests passed, 0 tests failed out of 520

Shipping it.

Both RadialGrid and SphericalGrid were constructed but never used by anything except a somewhat obscure print situation. All real logic handled by Mngr classes. Cleaning up to make things easier to hack on in the future. All tests pass, so whatever this code was added for was possibly never implemented, or unimportant. If something obscure in 3rd party ISA tools breaks, though, this is the commit to (probably) look at...

psi4/src/psi4/libfock/cubature.cc

tallakahath · 2022-10-06T20:58:17Z

I think I have a version that maintains the nice printing functionality (and shoves the necessary info into a struct, and actually uses the nicer LebedevGridMgr::findOrderByNPoints function and other LebedevGridMgr functions) but I can't actually test if this builds until I get home due to #reasons. I'll shove it upstream on my local branch and see if it breaks in the test builder.

tallakahath · 2022-10-06T20:59:17Z

(Note that I expect this will probably fail to build or run but I cannot debug with an incremental build until late tonight at best, probably tomorrow)

tallakahath · 2022-10-06T22:46:56Z

Clearly I can't get this working right now and I won't be able to touch this from a useful computer until tomorrow. Sorry.

I think if rs[A] is so big that it's length won't fit in a regular old `int`, we have a bigger problem...

tallakahath · 2022-10-07T18:33:39Z

I now have a version that compiles at home, though I'm not getting the narrowing warning-as-error that the auto-builder is. Still I've patched for it. This should work now...

hokru · 2022-10-08T14:05:04Z

This PR does more than the description tells me. Why are we removing so much?

hokru · 2022-10-08T14:09:27Z

psi4/src/psi4/libfock/cubature.cc

-
-    return std::shared_ptr<RadialGrid>(grid);
-}
-std::shared_ptr<RadialGrid> RadialGrid::build_treutler(int npoints, double alpha, int Z) {


This is the only Treutler version with correct Eta values. It does give slightly different results compared to the currently used "Treutler" where Eta is always 1. Though I had no time to try to confirm if this one is fully correct.

Even if it's correct, it was never actually used anywhere. The only use of the RadialGrid built was to provide data in print_details -- it's not accessed by anything or anyone outside of that function, at least as far as all of the psi4 tests can determine. In addition, the only calls to RadialGrid::build specifically were of the type RadialGrid::build("Unknown, ... which never invoked build_becke or build_treutler.

Yes, i tested it in a local build. So many parts of cubature.cc are messy unfortunately. Could you move the Eta table to line 2045 as a comment and add a note that the eta mapping is missing, please?
When I added the table long time ago I wasn't aware these routines were never called.

susilehtola · 2022-10-10T14:43:25Z

psi4/src/psi4/libfock/cubature.h

+    int npoints_;
+    /// Alpha scale (for user reference)
+    double alpha_;
+    /// Nodes (including alpha)
+    double* r_;
+    /// Weights (including alpha and r^2)
+    double* w_;


Please do not introduce raw pointers, when their use can easily be avoided. You can just duplicate the data into std::vector<double>.

Fixed (I think? My C++ is very rusty, if you haven't noticed...).

JonathonMisiewicz · 2022-10-10T15:12:02Z

psi4/src/psi4/libfock/cubature.cc

-            int Lsphere = spherical_grids_[A][R]->order();
-            printer->Printf("    Node: %4zu, R = %11.3E, WR = %11.3E, Nsphere = %6d, Lsphere = %6d\n", R, Rval, Wval,
-                            Nsphere, Lsphere);
+        printer->Printf("    Atom: %4zu, Nrad = %6d, Alpha = %11.3E:\n", A,


Just to be clear on this point: radial_grids_ exists now as a convenience struct for the sake of this printout, but the data being printed actually existed and is used elsewhere. Is that right, or am I missing something?

Yes, radial_grids_ is a convenience struct. I checked, and it looks like, e.g., r and wr do get used downstream, or at least, used in the construction of the MassPoints that get push_back'd into the atomic_grids_ of MolecularGrid. I suppose I haven't checked what happens, long-term, with those, but hopefully something...

Okay, in that case, does it make more sense to not print this at all but rely on the user to call MassPoint.print() or MassPoint.print_details() when necessary?

I looked into this and my conclusion was that grid pruning and cutoff thresholds also acts on the MassPoints, so recovering things like the initial lebedev grid or number of grid points isn't possible. The MassPoints don't really carry around sufficient information to reconstruct the state printed in print_details. But I'm not very familiar with the codebase OR what's happening under-the-hood here, so another set of eyes to check my assertion that this Cannot Be Done would be appreciated.

However, if ya'll are OK ditching print_details all together, then all these problems go away. The replacement structs were added b/c of @susilehtola 's requests to not break the existing print information.

Any information that isn't actually used except in printing, I'm okay removing.

tallakahath · 2022-10-10T19:28:13Z

I've made a new PR #2743 to JUST fix the hang bug, and then we can keep bikeshedding here on what we actually want to do about RadialGrid and SphericalGrid but have a fix in the codebase that actually stops the nasty hang that took me down this rabbit hole in the first place.

@susilehtola this initializes with an initializer list like you requested
@hokru this loses NO information over the prior state-of-the-world like you suggested was desirable (if I was reading correctly)
@JonathonMisiewicz I guess you won't like it because spaghetti didn't go away but it's short and sweet, at least.

Obviously this patch would undo all of that patch but I think the scope of this discussion exceeded the original PR statement.

JonathonMisiewicz · 2022-10-10T19:29:47Z

This is turning into a spaghetti PR, so breakup appreciated.

JonathonMisiewicz · 2022-10-11T18:04:05Z

Thanks for getting the most important part of the PR cleared away. As far as I can see, the next steps are:

Fix the merge conflicts (annoying rebase)
Determine what information in print_details is actually used. I want to retain print information that is actually used, but I can do away with the rest. Once we have that, we can assess if keeping a modified RadialGrid is the best way to keep the useful print information. (I think we've established that r and wr are used, but order is not?)

This part of the code is not pretty, and I commend you for dealing with it at all.

tallakahath · 2022-10-11T18:37:42Z

Fix the merge conflicts (annoying rebase)

I'm not sure how much rebasing you want to do, but I think things are back to the state-of-the-world that we were discussing prior to #2743

Determine what information in print_details is actually used. I want to retain print information that is actually used, but I can do away with the rest. Once we have that, we can assess if keeping a modified RadialGrid is the best way to keep the useful print information. (I think we've established that r and wr are used, but order is not?)

"Used" is a funny word here. No tests when this print is removed, so it's not so important that someone wrote a test to preserve it. Furthermore, the order_ was always wrong. The other elements weren't, though, but you have to set a pretty high print level to actually see anything (v.cc requires print_ > 2).

If by "used" you mean that the piece of information represented in the print is used downstream: r and wr are the primary candidates here. That said, they aren't used as-is -- they're dotted through other elements found by a lookup in LebedevGridMgr and RadialGridMgr to make MassPoints (see

psi4/psi4/src/psi4/libfock/cubature.cc

Lines 3814 to 3822 in 9db9100

    
           for (int j = 0; j < numAngPts; j++) { 
        
               MassPoint mp = {r[i] * anggrid[j].x, r[i] * anggrid[j].y, r[i] * anggrid[j].z, 
        
                               wr[i] * anggrid[j].w}; 
        
               mp = std_orientation.MoveIntoPosition(mp, A); 
        
               mp.w *= nuc.computeNuclearWeight(mp, A, stratmannCutoff); 
        
               if (std::abs(mp.w) > weightcut) { 
        
                   atomic_grids_[A].push_back(mp); 
        
               }

)

I'm not sure if we care about diagnostic print info on the intermediaries (the radii of shells around atoms we care about, and the number of angular points we end up picking for each shell). Maybe we do? If so, a struct seems like the most reasonable solution. I suppose, all-in-all, it's not so much data to carry around. Appending that data onto every MassPoint however feels like overkill (and that's the only other place to put it).

But we definitely should clean up/remove the rest of the old RadialGrid and SphericalGrid, as these seem to only cause confusion about where changes should be made. RadialGridMgr and LebedevGridMgr are the classes of Actual Importance.

JonathonMisiewicz · 2022-10-11T18:45:48Z

1 sounds good.

On 2., yes, agreed that high-level, having both Grid and GridMgr is confusing, and that's the issue to deal with. If theGrid classes are only used for print purposes, then while the print functionality belongs somewhere, it's not there. Keeping this old code around doesn't make that task easier doesn't make print setup in the correct place any easier. We should then just be rid of the classes entirely, printing included.

Are the Grid classes used for something other than printing?

Too many changes since original review

tallakahath · 2022-10-26T17:30:49Z

Are the Grid classes used for something other than printing?

Grid classes are only used for printing. Hence why I moved things to skinnier structs, so that it's clear to future devs that "nothing goes here but storage" (and maybe even with a big comment about "THIS CONTENT IS NOT REFERENCED EXCEPT BY A PRINT". There's no reasonable way I can see to migrate this into the MassPoints.

I've been quiet b/c I am hoping the other primary devs will chime in here -- I'm not sure how to move forwards given conflicting desires (and given that the immediate urgency of "this causes program hangs" has been fixed).

hokru · 2022-10-26T18:24:00Z

This print_details is more a debug printing for devs than for users. The quadrature points themselves can be printed with the BlockOPoints objects if needed. The r and wr are intermediates.

I would be in favour of removing the classes at the cost of removing the detailed printing.
It's better to re-think how to print this information instead of trying to fit it into these classes.

susilehtola reviewed Oct 6, 2022

View reviewed changes

psi4/src/psi4/libfock/cubature.cc Outdated Show resolved Hide resolved

Merge branch 'psi4:master' into initialize_lebedev

b8cad0f

JonathonMisiewicz previously approved these changes Oct 6, 2022

View reviewed changes

Revive printing of advanced grid information

af4d409

Elizabeth added 5 commits October 6, 2022 17:36

Whoops, forgot to uncomment in the print_details in the header

ed29cb5

And one more obvious typo

4d6000c

Last attempt at singular whack-a-mole

b30769f

OK one more

dceb5ea

I hope I can squash these commits together when it's Go Time

15e7ece

Codebase unlikes the implicit narrowing

ab799f2

I think if rs[A] is so big that it's length won't fit in a regular old `int`, we have a bigger problem...

hokru reviewed Oct 8, 2022

View reviewed changes

Adding record of correct eta mapping per @hokru's request

c5de09f

susilehtola reviewed Oct 10, 2022

View reviewed changes

JonathonMisiewicz reviewed Oct 10, 2022

View reviewed changes

Per @susilehtola 's request

7fd286b

tallakahath changed the title ~~Ensure that initialize_lebedev is called~~ Ensure that initialize_lebedev is called --> Clean up spaghetti in cubature Oct 10, 2022

tallakahath mentioned this pull request Oct 10, 2022

Initialize lebedev_mapping_ and make use thread-safe -- and do ONLY THAT #2743

Merged

Merge branch 'upstream_master' into initialize_lebedev

587558a

tallakahath mentioned this pull request Oct 26, 2022

quadrature grid pruning with functions broken #2754

Closed

JonathonMisiewicz mentioned this pull request Nov 3, 2022

SphericalGrid::order Contract Broken #2765

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Ensure that initialize_lebedev is called --> Clean up spaghetti in cubature #2736

Ensure that initialize_lebedev is called --> Clean up spaghetti in cubature #2736

tallakahath commented Oct 5, 2022 •

edited

susilehtola commented Oct 5, 2022 •

edited

loriab commented Oct 5, 2022

tallakahath commented Oct 5, 2022 via email

susilehtola commented Oct 5, 2022

tallakahath commented Oct 5, 2022

susilehtola commented Oct 5, 2022

loriab commented Oct 5, 2022

JonathonMisiewicz commented Oct 5, 2022

tallakahath commented Oct 5, 2022

tallakahath commented Oct 5, 2022

tallakahath commented Oct 5, 2022

tallakahath commented Oct 6, 2022

tallakahath commented Oct 6, 2022

tallakahath commented Oct 6, 2022

tallakahath commented Oct 7, 2022

hokru commented Oct 8, 2022

hokru Oct 8, 2022

tallakahath Oct 8, 2022

hokru Oct 8, 2022

tallakahath Oct 10, 2022

susilehtola Oct 10, 2022

tallakahath Oct 10, 2022

JonathonMisiewicz Oct 10, 2022

tallakahath Oct 10, 2022

JonathonMisiewicz Oct 10, 2022

tallakahath Oct 10, 2022 •

edited

JonathonMisiewicz Oct 10, 2022

tallakahath commented Oct 10, 2022

JonathonMisiewicz commented Oct 10, 2022

JonathonMisiewicz commented Oct 11, 2022

tallakahath commented Oct 11, 2022

JonathonMisiewicz commented Oct 11, 2022 •

edited

tallakahath commented Oct 26, 2022

hokru commented Oct 26, 2022

Ensure that initialize_lebedev is called --> Clean up spaghetti in cubature #2736

Are you sure you want to change the base?

Ensure that initialize_lebedev is called --> Clean up spaghetti in cubature #2736

Conversation

tallakahath commented Oct 5, 2022 • edited

Description

User API & Changelog headlines

Dev notes & details

Questions

Checklist

Status

susilehtola commented Oct 5, 2022 • edited

loriab commented Oct 5, 2022

tallakahath commented Oct 5, 2022 via email

susilehtola commented Oct 5, 2022

tallakahath commented Oct 5, 2022

susilehtola commented Oct 5, 2022

loriab commented Oct 5, 2022

JonathonMisiewicz commented Oct 5, 2022

tallakahath commented Oct 5, 2022

tallakahath commented Oct 5, 2022

tallakahath commented Oct 5, 2022

tallakahath commented Oct 6, 2022

tallakahath commented Oct 6, 2022

tallakahath commented Oct 6, 2022

tallakahath commented Oct 7, 2022

hokru commented Oct 8, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tallakahath Oct 10, 2022 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tallakahath commented Oct 10, 2022

JonathonMisiewicz commented Oct 10, 2022

JonathonMisiewicz commented Oct 11, 2022

tallakahath commented Oct 11, 2022

JonathonMisiewicz commented Oct 11, 2022 • edited

tallakahath commented Oct 26, 2022

hokru commented Oct 26, 2022

tallakahath commented Oct 5, 2022 •

edited

susilehtola commented Oct 5, 2022 •

edited

tallakahath Oct 10, 2022 •

edited

JonathonMisiewicz commented Oct 11, 2022 •

edited