Surface finding faster (vector->map)

[jpbrodsky]

This change improves performance when the number of G4LogicalBorderSurfaces is large.

The G4LogicalBorderSurfaceTable is searched frequently in optical photon simulations. When this table is large because many border surfaces have been defined, the performance of that search can dominate the performance of the simulation as a whole. The table can easily become large since border surfaces are defined between physical volumes and it's common for there to be many PhysVolume copies of the same logical volume.

This change improves performance of this search by changing the type of the G4LogicalBorderSurfaceTable from
std::vector<G4LogicalBorderSurface*>
to
std::map<std::pair<const G4VPhysicalVolume*,const G4VPhysicalVolume*>, G4LogicalBorderSurface*>
since that matches the most common way this table is searched.

This change affects one area outside the internal behavior of G4LogicalBorderSurface, which is in the GDML output of border surfaces. I have some concerns because the GDML search through the border surfaces searches only by the first PhysVolume in the pair that defines the border surface. Might this result in GDML output that misses the fact that multiple BorderSurfaces can be defined from the same PhysVolume? I have left the current behavior of the GDML output unchanged, but it should be examined.

[arekfu]

Hi, sorry to jump in on the discussion (I am not really qualified to review this PR), but have you tried using an std::unordered_map instead of an std::map? I would expect even larger speedups.

[jpbrodsky]

Hi, sorry to jump in on the discussion (I am not really qualified to review this PR), but have you tried using an std::unordered_map instead of an std::map? I would expect even larger speedups.

That's a good point. I haven't tested an unordered map but it should be superior.

Edit: this would require defining a hash function for std::pair<const G4VPhysicalVolume*,const G4VPhysicalVolume*>. This is probably not hard, but it's outside my area of expertise. If you're not careful, your hash function can be the new bottleneck, since it'll be called on every access of unordered_map elements. If it's a little slow, you can have good scaling with the unordered map but such a significant hit from the hash function that you lose any benefit.

[jrmadsen]

I am not 100% percent familiar with this section of code but from glancing at your changes, it is likely that

template <typename _Key, typename _Obj> using uomap = std::unordered_map<_Key, _Obj>;

uomap<const G4VPhysicalVolume*, uomap<const G4VPhysicalVolume*, const G4LogicalBorderSurface*>>

// ignoring checks this would be
auto surface = borderSurfaceTable[volume1][volume2]; 

is the best solution w.r.t. performance. I expect if you tried to implement unordered_map on that pair, you would have to create a custom hash operator or make it a pointer. You'd have to insert twice into the map (e.g. create entries for borderSurfaceTable[volume1][volume2]; and borderSurfaceTable[volume2][volume1]; but you are only storing pointers so that's isn't a huge memory overhead. If you run that and find a tangible speed-up, let me know and will talk to the geom working group coordinator.

[jpbrodsky]

You'd have to insert twice into the map (e.g. create entries for borderSurfaceTable[volume1][volume2]; and borderSurfaceTable[volume2][volume1];

In the current design, surfaces are defined one-way FROM one volume TO another. In many cases, only one direction is necessary, since one volume is an opaque reflector. So, this is not an issue with your implementation using nested unordered maps.

I'll give nested unordered maps a shot and see if it works well.

[gcosmo]

The proposed changes are included in the development branch and will be available in the next Geant4 release.
Thanks for the contribution!

[hanswenzel]

Sorry I missed this

using std::map instead of std::vector to define G4LogicalBorderSurfaceTable really is a change in the the API

using G4LogicalBorderSurfaceTable
= std::map<std::pair<const G4VPhysicalVolume*,
const G4VPhysicalVolume*>, G4LogicalBorderSurface*>;

instead of previously:

using G4LogicalBorderSurfaceTable = std::vector<G4LogicalBorderSurface*>;

so the method
static const G4LogicalBorderSurfaceTable* GetSurfaceTable();
in G4LogicalBorderSurface.hh

returns a completely different object. Instead of an ordered container (vector) where one can use the index (e.g. Opticks depends on this) as a pointer one gets an unordered container that behaves completely differently. I would suggest to have a new memeber function GetSurfaceMap that returns a G4LogicalBorderSurfaceMap. As is this really is a change of the API

[gcosmo]

Hi @hanswenzel, you're right, it is an API change although done for a good reason, for which I am sure that also Opticks would benefit for cases where many optical surfaces are defined. Is there anything impeding to update Opticks to adapt to the new signature compatible with the new Geant4 release?
I suggest to close this PR and if you have any modification to suggest, please, open a MR directly in the development repository, as you should have access to it. Thanks.